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   TCollection_AsciiString aFilePath = GetFilePath();
  64   if ( !aFilePath.IsEmpty() )
  65   {
  66     aResList << QString( "%1.ImportFromFile( \"%2\" );" )
  67                 .arg( aBathymetryName ).arg( aFilePath.ToCString() );
  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.Length() * 3 - 1 );
  87
  88   AltitudePoints::Iterator anIter( thePoints );
  89   for ( int i = 0 ; anIter.More(); ++i, anIter.Next() )
  90   {
  91     const AltitudePoint& aPoint = anIter.Value();
  92
  93     aCoordsArray->SetValue( i * 3, aPoint.X() );
  94     aCoordsArray->SetValue( i * 3 + 1, aPoint.Y() );
  95     aCoordsArray->SetValue( i * 3 + 2, aPoint.Z() );
  96   }
  97
  98   SetToUpdate( true );
  99 }
 100
 101 HYDROData_Bathymetry::AltitudePoints HYDROData_Bathymetry::GetAltitudePoints() const
 102 {
 103   AltitudePoints aPoints;
 104
 105   TDF_Label aLabel = myLab.FindChild( DataTag_AltitudePoints, false );
 106   if ( aLabel.IsNull() )
 107     return aPoints;
 108
 109   Handle(TDataStd_RealArray) aCoordsArray;
 110   if ( !aLabel.FindAttribute( TDataStd_RealArray::GetID(), aCoordsArray ) )
 111     return aPoints;
 112
 113   for ( int i = aCoordsArray->Lower(), n = aCoordsArray->Upper(); i <= n; )
 114   {
 115     if ( i + 3 > n + 1 )
 116       break;
 117
 118     AltitudePoint aPoint;
 119     aPoint.SetX( aCoordsArray->Value( i++ ) );
 120     aPoint.SetY( aCoordsArray->Value( i++ ) );
 121     aPoint.SetZ( aCoordsArray->Value( i++ ) );
 122
 123     aPoints.Append( aPoint );
 124   }
 125
 126   return aPoints;
 127 }
 128
 129 void HYDROData_Bathymetry::RemoveAltitudePoints()
 130 {
 131   TDF_Label aLabel = myLab.FindChild( DataTag_AltitudePoints, false );
 132   if ( !aLabel.IsNull() )
 133   {
 134     aLabel.ForgetAllAttributes();
 135     SetToUpdate( true );
 136   }
 137 }
 138
 139 void interpolateAltitudeForPoints( const gp_XY&                               thePoint,
 140                                    const HYDROData_Bathymetry::AltitudePoint& theFirstPoint,
 141                                    const HYDROData_Bathymetry::AltitudePoint& theSecPoint,
 142                                    HYDROData_Bathymetry::AltitudePoint&       theResPoint,
 143                                    const bool&                                theIsVertical )
 144 {
 145   double aCoordX = thePoint.X();
 146   double aCoordY = thePoint.Y();
 147
 148   if ( theIsVertical )
 149   {
 150     aCoordX = theFirstPoint.X();
 151
 152     if ( !ValuesEquals( theFirstPoint.X(), theSecPoint.X() ) )
 153     {
 154       // Recalculate X coordinate by equation of line from two points
 155       aCoordX = ( ( ( thePoint.Y() - theFirstPoint.Y() ) * ( theSecPoint.X() - theFirstPoint.X() ) ) /
 156                   ( theSecPoint.Y() - theFirstPoint.Y() ) ) + theFirstPoint.X();
 157     }
 158   }
 159   else
 160   {
 161     aCoordY = theFirstPoint.Y();
 162
 163     if ( !ValuesEquals( theFirstPoint.Y(), theSecPoint.Y() ) )
 164     {
 165       // Recalculate y by equation of line from two points
 166       aCoordY = ( ( ( thePoint.X() - theFirstPoint.X() ) * ( theSecPoint.Y() - theFirstPoint.Y() ) ) /
 167                   ( theSecPoint.X() - theFirstPoint.X() ) ) + theFirstPoint.Y();
 168     }
 169   }
 170
 171   theResPoint.SetX( aCoordX );
 172   theResPoint.SetY( aCoordY );
 173
 174   // Calculate coefficient for interpolation
 175   double aLength = Sqrt( Pow( theSecPoint.Y() - theFirstPoint.Y(), 2 ) +
 176                          Pow( theSecPoint.X() - theFirstPoint.X(), 2 ) );
 177
 178   double aInterCoeff = 0;
 179   if ( aLength != 0 )
 180    aInterCoeff = ( theSecPoint.Z() - theFirstPoint.Z() ) / aLength;
 181
 182
 183   double aNewLength = Sqrt( Pow( theResPoint.Y() - theFirstPoint.Y(), 2 ) +
 184                             Pow( theResPoint.X() - theFirstPoint.X(), 2 ) );
 185
 186   // Calculate interpolated value
 187   double aResVal = theFirstPoint.Z() + aInterCoeff * aNewLength;
 188
 189   theResPoint.SetZ( aResVal );
 190 }
 191
 192 double HYDROData_Bathymetry::GetAltitudeForPoint( const gp_XY& thePoint ) const
 193 {
 194   double anInvalidAltitude = GetInvalidAltitude();
 195   double aResAltitude = anInvalidAltitude;
 196
 197   AltitudePoints anAltitudePoints = GetAltitudePoints();
 198   if ( anAltitudePoints.IsEmpty() )
 199     return aResAltitude;
 200
 201   QPolygonF aBoundingRect;
 202
 203   // Boundary plane
 204   // [ 0 (top-left) ]          [ 1 (top-right) ]
 205   //                  thePoint
 206   // [ 2 (bot-left) ]          [ 3 (bot-right) ]
 207   AltitudePoint aBounds[ 4 ] = { AltitudePoint( -DBL_MAX, -DBL_MAX, anInvalidAltitude ),
 208                                  AltitudePoint(  DBL_MAX, -DBL_MAX, anInvalidAltitude ),
 209                                  AltitudePoint( -DBL_MAX,  DBL_MAX, anInvalidAltitude ),
 210                                  AltitudePoint(  DBL_MAX,  DBL_MAX, anInvalidAltitude ) };
 211
 212   AltitudePoints::Iterator anIter( anAltitudePoints );
 213   for ( ; anIter.More(); anIter.Next() )
 214   {
 215     const AltitudePoint& aPoint = anIter.Value();
 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 TCollection_AsciiString& theFilePath )
 329 {
 330   TDataStd_AsciiString::Set( myLab.FindChild( DataTag_FilePath ), theFilePath );
 331 }
 332
 333 TCollection_AsciiString HYDROData_Bathymetry::GetFilePath() const
 334 {
 335   TCollection_AsciiString aRes;
 336
 337   TDF_Label aLabel = myLab.FindChild( DataTag_FilePath, false );
 338   if ( !aLabel.IsNull() )
 339   {
 340     Handle(TDataStd_AsciiString) anAsciiStr;
 341     if ( aLabel.FindAttribute( TDataStd_AsciiString::GetID(), anAsciiStr ) )
 342       aRes = anAsciiStr->Get();
 343   }
 344
 345   return aRes;
 346 }
 347
 348 void HYDROData_Bathymetry::SetAltitudesInverted( const bool theIsInverted,
 349                                                  const bool theIsUpdate )
 350 {
 351   bool anIsAltitudesInverted = IsAltitudesInverted();
 352   if ( anIsAltitudesInverted == theIsInverted )
 353     return;
 354
 355   TDataStd_Integer::Set( myLab.FindChild( DataTag_AltitudesInverted ), (Standard_Integer)theIsInverted );
 356
 357   SetToUpdate( true );
 358
 359   if ( !theIsUpdate )
 360     return;
 361
 362   // Update altitude points
 363   AltitudePoints anAltitudePoints = GetAltitudePoints();
 364   if ( anAltitudePoints.IsEmpty() )
 365     return;
 366
 367   AltitudePoints::Iterator anIter( anAltitudePoints );
 368   for ( ; anIter.More(); anIter.Next() )
 369   {
 370     AltitudePoint& aPoint = anIter.ChangeValue();
 371     aPoint.SetZ( aPoint.Z() * -1 );
 372   }
 373
 374   SetAltitudePoints( anAltitudePoints );
 375 }
 376
 377 bool HYDROData_Bathymetry::IsAltitudesInverted() const
 378 {
 379   bool aRes = false;
 380
 381   TDF_Label aLabel = myLab.FindChild( DataTag_AltitudesInverted, false );
 382   if ( !aLabel.IsNull() )
 383   {
 384     Handle(TDataStd_Integer) anIntVal;
 385     if ( aLabel.FindAttribute( TDataStd_Integer::GetID(), anIntVal ) )
 386       aRes = (bool)anIntVal->Get();
 387   }
 388
 389   return aRes;
 390 }
 391
 392 bool HYDROData_Bathymetry::ImportFromFile( const TCollection_AsciiString& theFileName )
 393 {
 394   // Try to open the file
 395   QFile aFile( theFileName.ToCString() );
 396   if ( !aFile.exists() || !aFile.open( QIODevice::ReadOnly ) )
 397     return false;
 398
 399   bool aRes = false;
 400
 401   QString aFileSuf = QFileInfo( aFile ).suffix().toLower();
 402
 403   AltitudePoints aPoints;
 404
 405   // Try to import the file
 406   if ( aFileSuf == "xyz" )
 407     aRes = importFromXYZFile( aFile, aPoints );
 408
 409   // Close the file
 410   aFile.close();
 411
 412   if ( aRes )
 413   {
 414     // Update file path and altitude points of this Bathymetry
 415     SetFilePath( theFileName );
 416     SetAltitudePoints( aPoints );
 417   }
 418
 419   return aRes && !aPoints.IsEmpty();
 420 }
 421
 422 bool HYDROData_Bathymetry::importFromXYZFile( QFile&          theFile,
 423                                               AltitudePoints& thePoints ) const
 424 {
 425   if ( !theFile.isOpen() )
 426     return false;
 427
 428   // Strings in file is written as:
 429   //  1. X(float) Y(float) Z(float)
 430   //  2. X(float) Y(float) Z(float)
 431   //  ...
 432
 433 #ifdef _TIMER
 434   OSD_Timer aTimer;
 435   aTimer.Start();
 436 #endif
 437
 438   bool anIsAltitudesInverted = IsAltitudesInverted();
 439   while ( !theFile.atEnd() )
 440   {
 441     QString aLine = theFile.readLine().simplified();
 442     if ( aLine.isEmpty() )
 443       continue;
 444
 445     QStringList aValues = aLine.split( ' ', QString::SkipEmptyParts );
 446     if ( aValues.length() < 3 )
 447       return false;
 448
 449     AltitudePoint aPoint;
 450
 451     QString anX = aValues.value( 0 );
 452     QString anY = aValues.value( 1 );
 453     QString aZ  = aValues.value( 2 );
 454
 455     bool isXOk = false, isYOk = false, isZOk = false;
 456
 457     aPoint.SetX( anX.toDouble( &isXOk ) );
 458     aPoint.SetY( anY.toDouble( &isYOk ) );
 459     aPoint.SetZ( aZ.toDouble( &isZOk ) );
 460
 461     if ( !isXOk || !isYOk || !isZOk )
 462       return false;
 463
 464     if ( boost::math::isnan( aPoint.X() ) || boost::math::isinf( aPoint.X() ) ||
 465          boost::math::isnan( aPoint.Y() ) || boost::math::isinf( aPoint.Y() ) ||
 466          boost::math::isnan( aPoint.Z() ) || boost::math::isinf( aPoint.Z() ) )
 467       return false;
 468
 469     // Invert the z value if requested
 470     if ( anIsAltitudesInverted )
 471       aPoint.SetZ( -aPoint.Z() );
 472
 473     thePoints.Append( aPoint );
 474   }
 475
 476 #ifdef _TIMER
 477   aTimer.Stop();
 478   std::ofstream stream( "W:/HYDRO/WORK/log.txt", std::ofstream::out );
 479   aTimer.Show( stream );
 480 #endif
 481
 482   return true;
 483 }
 484
 485
 486 bool HYDROData_Bathymetry::CreateBoundaryPolyline() const
 487 {
 488   Handle(HYDROData_Document) aDocument = HYDROData_Document::Document( myLab );
 489   Handle_HYDROData_PolylineXY aResult =
 490     Handle_HYDROData_PolylineXY::DownCast( aDocument->CreateObject( KIND_POLYLINEXY ) );
 491
 492   if( aResult.IsNull() )
 493     return false;
 494
 495   //search free name
 496   QString aPolylinePref = GetName() + "_Boundary";
 497   QString aPolylineName = HYDROData_Tool::GenerateObjectName( aDocument, aPolylinePref );
 498   aResult->SetName( aPolylineName );
 499
 500   double Xmin = 0.0, Xmax = 0.0, Ymin = 0.0, Ymax = 0.0;
 501   bool isFirst = true;
 502   AltitudePoints aPoints = GetAltitudePoints();
 503
 504   AltitudePoints::Iterator anIter( aPoints );
 505   for ( ; anIter.More(); anIter.Next() )
 506   {
 507     const AltitudePoint& aPoint = anIter.Value();
 508
 509     double x = aPoint.X(), y = aPoint.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 }