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