src/HYDROData/HYDROData_Bathymetry.cxx

   1
   2 #include "HYDROData_Bathymetry.h"
   3 #include "HYDROData_Document.h"
   4 #include "HYDROData_Tool.h"
   5 #include "HYDROData_PolylineXY.h"
   6
   7 #include <boost/math/special_functions/fpclassify.hpp>
   8
   9 #include <gp_XY.hxx>
  10 #include <gp_XYZ.hxx>
  11
  12 #include <TDataStd_RealArray.hxx>
  13 #include <TDataStd_AsciiString.hxx>
  14 #include <TDataStd_Integer.hxx>
  15
  16 #include <QFile>
  17 #include <QFileInfo>
  18 #include <QPointF>
  19 #include <QPolygonF>
  20 #include <QStringList>
  21
  22 #include <math.h>
  23
  24 #define _TIMER
  25 #ifdef _TIMER
  26 #include <OSD_Timer.hxx>
  27 #endif
  28
  29 #define PYTHON_BATHYMETRY_ID "KIND_BATHYMETRY"
  30
  31
  32 IMPLEMENT_STANDARD_HANDLE(HYDROData_Bathymetry, HYDROData_IAltitudeObject)
  33 IMPLEMENT_STANDARD_RTTIEXT(HYDROData_Bathymetry, HYDROData_IAltitudeObject)
  34
  35 HYDROData_Bathymetry::HYDROData_Bathymetry()
  36 : HYDROData_IAltitudeObject()
  37 {
  38 }
  39
  40 HYDROData_Bathymetry::~HYDROData_Bathymetry()
  41 {
  42 }
  43
  44 QStringList HYDROData_Bathymetry::DumpToPython( MapOfTreatedObjects& theTreatedObjects ) const
  45 {
  46   QStringList aResList;
  47
  48   Handle(HYDROData_Document) aDocument = HYDROData_Document::Document( myLab );
  49   if ( aDocument.IsNull() )
  50     return aResList;
  51
  52   QString aDocName = aDocument->GetDocPyName();
  53   QString aBathymetryName = GetName();
  54
  55   aResList << QString( "%1 = %2.CreateObject( %3 );" )
  56               .arg( aBathymetryName ).arg( aDocName ).arg( PYTHON_BATHYMETRY_ID );
  57   aResList << QString( "%1.SetName( \"%2\" );" )
  58               .arg( aBathymetryName ).arg( aBathymetryName );
  59
  60   aResList << QString( "%1.SetAltitudesInverted( %2 );" )
  61               .arg( aBathymetryName ).arg( IsAltitudesInverted() );
  62
  63   QString aFilePath = GetFilePath();
  64   if ( !aFilePath.isEmpty() )
  65   {
  66     aResList << QString( "%1.ImportFromFile( \"%2\" );" )
  67                 .arg( aBathymetryName ).arg( aFilePath );
  68   }
  69   else
  70   {
  71     // TODO : bathymetry is composed from other bathymetry(ies)
  72   }
  73
  74   return aResList;
  75 }
  76
  77 void HYDROData_Bathymetry::SetAltitudePoints( const AltitudePoints& thePoints )
  78 {
  79   RemoveAltitudePoints();
  80
  81   if ( thePoints.isEmpty() )
  82     return;
  83
  84   // Save coordinates
  85   Handle(TDataStd_RealArray) aCoordsArray =
  86     TDataStd_RealArray::Set( myLab.FindChild( DataTag_AltitudePoints ), 0, thePoints.size() * 3 - 1 );
  87
  88   AltitudePoints::const_iterator aListItBeg =  thePoints.constBegin();
  89   AltitudePoints::const_iterator aListItEnd =  thePoints.constEnd();
  90   for ( int i = 0 ; aListItBeg != aListItEnd; ++i, ++aListItBeg )
  91   {
  92     const AltitudePoint& aPoint = *aListItBeg;
  93
  94     aCoordsArray->SetValue( i * 3, aPoint.X() );
  95     aCoordsArray->SetValue( i * 3 + 1, aPoint.Y() );
  96     aCoordsArray->SetValue( i * 3 + 2, aPoint.Z() );
  97   }
  98
  99   SetToUpdate( true );
 100 }
 101
 102 HYDROData_Bathymetry::AltitudePoints HYDROData_Bathymetry::GetAltitudePoints() const
 103 {
 104   AltitudePoints aPoints;
 105
 106   TDF_Label aLabel = myLab.FindChild( DataTag_AltitudePoints, false );
 107   if ( aLabel.IsNull() )
 108     return aPoints;
 109
 110   Handle(TDataStd_RealArray) aCoordsArray;
 111   if ( !aLabel.FindAttribute( TDataStd_RealArray::GetID(), aCoordsArray ) )
 112     return aPoints;
 113
 114   for ( int i = aCoordsArray->Lower(), n = aCoordsArray->Upper(); i <= n; )
 115   {
 116     if ( i + 3 > n + 1 )
 117       break;
 118
 119     AltitudePoint aPoint;
 120     aPoint.SetX( aCoordsArray->Value( i++ ) );
 121     aPoint.SetY( aCoordsArray->Value( i++ ) );
 122     aPoint.SetZ( aCoordsArray->Value( i++ ) );
 123
 124     aPoints << aPoint;
 125   }
 126
 127   return aPoints;
 128 }
 129
 130 void HYDROData_Bathymetry::RemoveAltitudePoints()
 131 {
 132   TDF_Label aLabel = myLab.FindChild( DataTag_AltitudePoints, false );
 133   if ( !aLabel.IsNull() )
 134   {
 135     aLabel.ForgetAllAttributes();
 136     SetToUpdate( true );
 137   }
 138 }
 139
 140 void interpolateAltitudeForPoints( const gp_XY&                               thePoint,
 141                                    const HYDROData_Bathymetry::AltitudePoint& theFirstPoint,
 142                                    const HYDROData_Bathymetry::AltitudePoint& theSecPoint,
 143                                    HYDROData_Bathymetry::AltitudePoint&       theResPoint,
 144                                    const bool&                                theIsVertical )
 145 {
 146   double aCoordX = thePoint.X();
 147   double aCoordY = thePoint.Y();
 148
 149   if ( theIsVertical )
 150   {
 151     aCoordX = theFirstPoint.X();
 152
 153     if ( !ValuesEquals( theFirstPoint.X(), theSecPoint.X() ) )
 154     {
 155       // Recalculate X coordinate by equation of line from two points
 156       aCoordX = ( ( ( thePoint.Y() - theFirstPoint.Y() ) * ( theSecPoint.X() - theFirstPoint.X() ) ) /
 157                   ( theSecPoint.Y() - theFirstPoint.Y() ) ) + theFirstPoint.X();
 158     }
 159   }
 160   else
 161   {
 162     aCoordY = theFirstPoint.Y();
 163
 164     if ( !ValuesEquals( theFirstPoint.Y(), theSecPoint.Y() ) )
 165     {
 166       // Recalculate y by equation of line from two points
 167       aCoordY = ( ( ( thePoint.X() - theFirstPoint.X() ) * ( theSecPoint.Y() - theFirstPoint.Y() ) ) /
 168                   ( theSecPoint.X() - theFirstPoint.X() ) ) + theFirstPoint.Y();
 169     }
 170   }
 171
 172   theResPoint.SetX( aCoordX );
 173   theResPoint.SetY( aCoordY );
 174
 175   // Calculate coefficient for interpolation
 176   double aLength = Sqrt( Pow( theSecPoint.Y() - theFirstPoint.Y(), 2 ) +
 177                          Pow( theSecPoint.X() - theFirstPoint.X(), 2 ) );
 178
 179   double aInterCoeff = 0;
 180   if ( aLength != 0 )
 181    aInterCoeff = ( theSecPoint.Z() - theFirstPoint.Z() ) / aLength;
 182
 183
 184   double aNewLength = Sqrt( Pow( theResPoint.Y() - theFirstPoint.Y(), 2 ) +
 185                             Pow( theResPoint.X() - theFirstPoint.X(), 2 ) );
 186
 187   // Calculate interpolated value
 188   double aResVal = theFirstPoint.Z() + aInterCoeff * aNewLength;
 189
 190   theResPoint.SetZ( aResVal );
 191 }
 192
 193 double HYDROData_Bathymetry::GetAltitudeForPoint( const gp_XY& thePoint ) const
 194 {
 195   double anInvalidAltitude = GetInvalidAltitude();
 196   double aResAltitude = anInvalidAltitude;
 197
 198   AltitudePoints anAltitudePoints = GetAltitudePoints();
 199   if ( anAltitudePoints.isEmpty() )
 200     return aResAltitude;
 201
 202   QPolygonF aBoundingRect;
 203
 204   // Boundary plane
 205   // [ 0 (top-left) ]          [ 1 (top-right) ]
 206   //                  thePoint
 207   // [ 2 (bot-left) ]          [ 3 (bot-right) ]
 208   AltitudePoint aBounds[ 4 ] = { AltitudePoint( -DBL_MAX, -DBL_MAX, anInvalidAltitude ),
 209                                  AltitudePoint(  DBL_MAX, -DBL_MAX, anInvalidAltitude ),
 210                                  AltitudePoint( -DBL_MAX,  DBL_MAX, anInvalidAltitude ),
 211                                  AltitudePoint(  DBL_MAX,  DBL_MAX, anInvalidAltitude ) };
 212
 213   AltitudePoints::const_iterator aListItBeg = anAltitudePoints.constBegin();
 214   AltitudePoints::const_iterator aListItEnd = anAltitudePoints.constEnd();
 215   for ( ; aListItBeg != aListItEnd; ++aListItBeg )
 216   {
 217     const AltitudePoint& aPoint = *aListItBeg;
 218
 219     double aDeltaX = Abs( aPoint.X() ) - Abs( thePoint.X() );
 220     double aDeltaY = Abs( aPoint.Y() ) - Abs( thePoint.Y() );
 221
 222     if ( ValuesEquals( aDeltaX, 0.0 ) ) // Both left and right sides
 223     {
 224       if ( ValuesEquals( aDeltaY, 0.0 ) ) // Both top and bottom sides
 225       {
 226         aResAltitude = aPoint.Z();
 227         return aResAltitude;
 228       }
 229       else if ( aDeltaY < 0 ) // top side
 230       {
 231         // top border
 232         if ( ValuesMoreEquals( aPoint.X(), aBounds[ 0 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 0 ].Y() ) )
 233           aBounds[ 0 ] = aPoint;
 234         if ( ValuesLessEquals( aPoint.X(), aBounds[ 1 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 1 ].Y() ) )
 235           aBounds[ 1 ] = aPoint;
 236       }
 237       else
 238       {
 239         // bottom border
 240         if ( ValuesMoreEquals( aPoint.X(), aBounds[ 2 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 2 ].Y() ) )
 241           aBounds[ 2 ] = aPoint;
 242         if ( ValuesLessEquals( aPoint.X(), aBounds[ 3 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 3 ].Y() ) )
 243           aBounds[ 3 ] = aPoint;
 244       }
 245     }
 246     else if ( aDeltaX < 0 ) // left side
 247     {
 248       if ( ValuesEquals( aDeltaY, 0.0 ) )
 249       {
 250         // Left border
 251         if ( ValuesMoreEquals( aPoint.X(), aBounds[ 0 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 0 ].Y() ) )
 252           aBounds[ 0 ] = aPoint;
 253         if ( ValuesMoreEquals( aPoint.X(), aBounds[ 2 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 2 ].Y() ) )
 254           aBounds[ 2 ] = aPoint;
 255       }
 256       else if ( aDeltaY < 0 )
 257       {
 258         // top left corner
 259         if ( ValuesMoreEquals( aPoint.X(), aBounds[ 0 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 0 ].Y() ) )
 260           aBounds[ 0 ] = aPoint;
 261       }
 262       else
 263       {
 264         // bottom left corner
 265         if ( ValuesMoreEquals( aPoint.X(), aBounds[ 2 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 2 ].Y() ) )
 266           aBounds[ 2 ] = aPoint;
 267       }
 268     }
 269     else // right side
 270     {
 271       if ( ValuesEquals( aDeltaY, 0.0 ) )
 272       {
 273         // Right border
 274         if ( ValuesLessEquals( aPoint.X(), aBounds[ 1 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 1 ].Y() ) )
 275           aBounds[ 1 ] = aPoint;
 276         if ( ValuesLessEquals( aPoint.X(), aBounds[ 3 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 3 ].Y() ) )
 277           aBounds[ 3 ] = aPoint;
 278       }
 279       else if ( aDeltaY < 0 )
 280       {
 281         // top right corner
 282         if ( ValuesLessEquals( aPoint.X(), aBounds[ 1 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 1 ].Y() ) )
 283           aBounds[ 1 ] = aPoint;
 284       }
 285       else
 286       {
 287         // bottom right corner
 288         if ( ValuesLessEquals( aPoint.X(), aBounds[ 3 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 3 ].Y() ) )
 289           aBounds[ 3 ] = aPoint;
 290       }
 291     }
 292
 293     // Update bounding rectangle of our global grid
 294     aBoundingRect << QPointF( aPoint.X(), aPoint.Y() );
 295   }
 296
 297   const double LIMIT = 1E300;
 298   if( fabs( aBounds[ 0 ].X() ) > LIMIT || fabs( aBounds[ 0 ].Y() ) > LIMIT ||
 299       fabs( aBounds[ 1 ].X() ) > LIMIT || fabs( aBounds[ 1 ].Y() ) > LIMIT ||
 300       fabs( aBounds[ 2 ].X() ) > LIMIT || fabs( aBounds[ 2 ].Y() ) > LIMIT ||
 301       fabs( aBounds[ 3 ].X() ) > LIMIT || fabs( aBounds[ 3 ].Y() ) > LIMIT )
 302     return anInvalidAltitude;
 303
 304
 305   // Check if requested point is inside of our bounding rectangle
 306   if ( !aBoundingRect.boundingRect().contains( thePoint.X(), thePoint.Y() ) )
 307     return aResAltitude;
 308
 309   // Calculate result altitude for point
 310   AltitudePoint aFirstPoint( aBounds[ 0 ] ), aSecPoint( aBounds[ 1 ] );
 311
 312   // At first we merge top and bottom borders
 313   if ( aBounds[ 0 ].Y() != aBounds[ 2 ].Y() || aBounds[ 0 ].X() != aBounds[ 2 ].X() )
 314     interpolateAltitudeForPoints( thePoint, aBounds[ 0 ], aBounds[ 2 ], aFirstPoint, true );
 315
 316   if ( aBounds[ 1 ].Y() != aBounds[ 3 ].Y() || aBounds[ 1 ].X() != aBounds[ 3 ].X() )
 317     interpolateAltitudeForPoints( thePoint, aBounds[ 1 ], aBounds[ 3 ], aSecPoint, true );
 318
 319   AltitudePoint aResPoint( aFirstPoint );
 320
 321   // At last we merge left and right borders
 322   if ( aFirstPoint.Y() != aSecPoint.Y() || aFirstPoint.X() != aSecPoint.X() )
 323     interpolateAltitudeForPoints( thePoint, aFirstPoint, aSecPoint, aResPoint, false );
 324
 325   aResAltitude = aResPoint.Z();
 326
 327   return aResAltitude;
 328 }
 329
 330 void HYDROData_Bathymetry::SetFilePath(const QString& theFilePath)
 331 {
 332   TCollection_AsciiString anAsciiStr( theFilePath.toStdString().c_str() );
 333   TDataStd_AsciiString::Set( myLab.FindChild( DataTag_FilePath ), anAsciiStr );
 334 }
 335
 336 QString HYDROData_Bathymetry::GetFilePath() const
 337 {
 338   QString aRes;
 339
 340   TDF_Label aLabel = myLab.FindChild( DataTag_FilePath, false );
 341   if ( !aLabel.IsNull() )
 342   {
 343     Handle(TDataStd_AsciiString) anAsciiStr;
 344     if ( aLabel.FindAttribute( TDataStd_AsciiString::GetID(), anAsciiStr ) )
 345       aRes = QString( anAsciiStr->Get().ToCString() );
 346   }
 347
 348   return aRes;
 349 }
 350
 351 void HYDROData_Bathymetry::SetAltitudesInverted( const bool theIsInverted,
 352                                                  const bool theIsUpdate )
 353 {
 354   bool anIsAltitudesInverted = IsAltitudesInverted();
 355   if ( anIsAltitudesInverted == theIsInverted )
 356     return;
 357
 358   TDataStd_Integer::Set( myLab.FindChild( DataTag_AltitudesInverted ), (Standard_Integer)theIsInverted );
 359
 360   SetToUpdate( true );
 361
 362   if ( !theIsUpdate )
 363     return;
 364
 365   // Update altitude points
 366   AltitudePoints anAltitudePoints = GetAltitudePoints();
 367   if ( anAltitudePoints.isEmpty() )
 368     return;
 369
 370   AltitudePoints::iterator aListItBeg = anAltitudePoints.begin();
 371   AltitudePoints::iterator aListItEnd = anAltitudePoints.end();
 372   for ( ; aListItBeg != aListItEnd; ++aListItBeg )
 373   {
 374     AltitudePoint& aPoint = *aListItBeg;
 375     aPoint.SetZ( aPoint.Z() * -1 );
 376   }
 377
 378   SetAltitudePoints( anAltitudePoints );
 379 }
 380
 381 bool HYDROData_Bathymetry::IsAltitudesInverted() const
 382 {
 383   bool aRes = false;
 384
 385   TDF_Label aLabel = myLab.FindChild( DataTag_AltitudesInverted, false );
 386   if ( !aLabel.IsNull() )
 387   {
 388     Handle(TDataStd_Integer) anIntVal;
 389     if ( aLabel.FindAttribute( TDataStd_Integer::GetID(), anIntVal ) )
 390       aRes = (bool)anIntVal->Get();
 391   }
 392
 393   return aRes;
 394 }
 395
 396 bool HYDROData_Bathymetry::ImportFromFile( const QString& theFileName )
 397 {
 398   // Try to open the file
 399   QFile aFile( theFileName );
 400   if ( !aFile.exists() || !aFile.open( QIODevice::ReadOnly ) )
 401     return false;
 402
 403   bool aRes = false;
 404
 405   QString aFileSuf = QFileInfo( aFile ).suffix().toLower();
 406
 407   AltitudePoints aPoints;
 408
 409   // Try to import the file
 410   if ( aFileSuf == "xyz" )
 411     aRes = importFromXYZFile( aFile, aPoints );
 412
 413   // Close the file
 414   aFile.close();
 415
 416   if ( aRes )
 417   {
 418     // Update file path and altitude points of this Bathymetry
 419     SetFilePath( theFileName );
 420     SetAltitudePoints( aPoints );
 421   }
 422
 423   return aRes && !aPoints.isEmpty();
 424 }
 425
 426 bool HYDROData_Bathymetry::importFromXYZFile( QFile&          theFile,
 427                                               AltitudePoints& thePoints ) const
 428 {
 429   if ( !theFile.isOpen() )
 430     return false;
 431
 432   // Strings in file is written as:
 433   //  1. X(float) Y(float) Z(float)
 434   //  2. X(float) Y(float) Z(float)
 435   //  ...
 436
 437 #ifdef _TIMER
 438   OSD_Timer aTimer;
 439   aTimer.Start();
 440 #endif
 441
 442   bool anIsAltitudesInverted = IsAltitudesInverted();
 443   while ( !theFile.atEnd() )
 444   {
 445     QString aLine = theFile.readLine().simplified();
 446     if ( aLine.isEmpty() )
 447       continue;
 448
 449     QStringList aValues = aLine.split( ' ', QString::SkipEmptyParts );
 450     if ( aValues.length() < 3 )
 451       return false;
 452
 453     AltitudePoint aPoint;
 454
 455     QString anX = aValues.value( 0 );
 456     QString anY = aValues.value( 1 );
 457     QString aZ  = aValues.value( 2 );
 458
 459     bool isXOk = false, isYOk = false, isZOk = false;
 460
 461     aPoint.SetX( anX.toDouble( &isXOk ) );
 462     aPoint.SetY( anY.toDouble( &isYOk ) );
 463     aPoint.SetZ( aZ.toDouble( &isZOk ) );
 464
 465     if ( !isXOk || !isYOk || !isZOk )
 466       return false;
 467
 468     if ( boost::math::isnan( aPoint.X() ) || boost::math::isinf( aPoint.X() ) ||
 469          boost::math::isnan( aPoint.Y() ) || boost::math::isinf( aPoint.Y() ) ||
 470          boost::math::isnan( aPoint.Z() ) || boost::math::isinf( aPoint.Z() ) )
 471       return false;
 472
 473     // Invert the z value if requested
 474     if ( anIsAltitudesInverted )
 475       aPoint.SetZ( -aPoint.Z() );
 476
 477     thePoints << aPoint;
 478   }
 479
 480 #ifdef _TIMER
 481   aTimer.Stop();
 482   std::ofstream stream( "W:/HYDRO/WORK/log.txt", std::ofstream::out );
 483   aTimer.Show( stream );
 484 #endif
 485
 486   return true;
 487 }
 488
 489
 490 bool HYDROData_Bathymetry::CreateBoundaryPolyline() const
 491 {
 492   Handle(HYDROData_Document) aDocument = HYDROData_Document::Document( myLab );
 493   Handle_HYDROData_PolylineXY aResult =
 494     Handle_HYDROData_PolylineXY::DownCast( aDocument->CreateObject( KIND_POLYLINEXY ) );
 495
 496   if( aResult.IsNull() )
 497     return false;
 498
 499   //search free name
 500   QString aPolylinePref = GetName() + "_Boundary";
 501   QString aPolylineName = HYDROData_Tool::GenerateObjectName( aDocument, aPolylinePref );
 502   aResult->SetName( aPolylineName );
 503
 504   double Xmin = 0.0, Xmax = 0.0, Ymin = 0.0, Ymax = 0.0;
 505   bool isFirst = true;
 506   AltitudePoints aPoints = GetAltitudePoints();
 507   foreach( AltitudePoint aPnt, aPoints )
 508   {
 509     double x = aPnt.X(), y = aPnt.Y();
 510     if( isFirst || x<Xmin )
 511       Xmin = x;
 512     if( isFirst || x>Xmax )
 513       Xmax = x;
 514     if( isFirst || y<Ymin )
 515       Ymin = y;
 516     if( isFirst || y>Ymax )
 517       Ymax = y;
 518     isFirst = false;
 519   }
 520
 521   aResult->AddSection( "bound", HYDROData_IPolyline::SECTION_POLYLINE, true );
 522   aResult->AddPoint( 0, HYDROData_IPolyline::Point( Xmin, Ymin ) );
 523   aResult->AddPoint( 0, HYDROData_IPolyline::Point( Xmin, Ymax ) );
 524   aResult->AddPoint( 0, HYDROData_IPolyline::Point( Xmax, Ymax ) );
 525   aResult->AddPoint( 0, HYDROData_IPolyline::Point( Xmax, Ymin ) );
 526   aResult->Update();
 527
 528   return true;
 529 }