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.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 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 anIter( anAltitudePoints );
 369   for ( ; anIter.More(); anIter.Next() )
 370   {
 371     AltitudePoint& aPoint = anIter.ChangeValue();
 372     aPoint.SetZ( aPoint.Z() * -1 );
 373   }
 374
 375   SetAltitudePoints( anAltitudePoints );
 376 }
 377
 378 bool HYDROData_Bathymetry::IsAltitudesInverted() const
 379 {
 380   bool aRes = false;
 381
 382   TDF_Label aLabel = myLab.FindChild( DataTag_AltitudesInverted, false );
 383   if ( !aLabel.IsNull() )
 384   {
 385     Handle(TDataStd_Integer) anIntVal;
 386     if ( aLabel.FindAttribute( TDataStd_Integer::GetID(), anIntVal ) )
 387       aRes = (bool)anIntVal->Get();
 388   }
 389
 390   return aRes;
 391 }
 392
 393 bool HYDROData_Bathymetry::ImportFromFile( const QString& theFileName )
 394 {
 395   // Try to open the file
 396   QFile aFile( theFileName );
 397   if ( !aFile.exists() || !aFile.open( QIODevice::ReadOnly ) )
 398     return false;
 399
 400   bool aRes = false;
 401
 402   QString aFileSuf = QFileInfo( aFile ).suffix().toLower();
 403
 404   AltitudePoints aPoints;
 405
 406   // Try to import the file
 407   if ( aFileSuf == "xyz" )
 408     aRes = importFromXYZFile( aFile, aPoints );
 409
 410   // Close the file
 411   aFile.close();
 412
 413   if ( aRes )
 414   {
 415     // Update file path and altitude points of this Bathymetry
 416     SetFilePath( theFileName );
 417     SetAltitudePoints( aPoints );
 418   }
 419
 420   return aRes && !aPoints.IsEmpty();
 421 }
 422
 423 bool HYDROData_Bathymetry::importFromXYZFile( QFile&          theFile,
 424                                               AltitudePoints& thePoints ) const
 425 {
 426   if ( !theFile.isOpen() )
 427     return false;
 428
 429   // Strings in file is written as:
 430   //  1. X(float) Y(float) Z(float)
 431   //  2. X(float) Y(float) Z(float)
 432   //  ...
 433
 434 #ifdef _TIMER
 435   OSD_Timer aTimer;
 436   aTimer.Start();
 437 #endif
 438
 439   bool anIsAltitudesInverted = IsAltitudesInverted();
 440   while ( !theFile.atEnd() )
 441   {
 442     QString aLine = theFile.readLine().simplified();
 443     if ( aLine.isEmpty() )
 444       continue;
 445
 446     QStringList aValues = aLine.split( ' ', QString::SkipEmptyParts );
 447     if ( aValues.length() < 3 )
 448       return false;
 449
 450     AltitudePoint aPoint;
 451
 452     QString anX = aValues.value( 0 );
 453     QString anY = aValues.value( 1 );
 454     QString aZ  = aValues.value( 2 );
 455
 456     bool isXOk = false, isYOk = false, isZOk = false;
 457
 458     aPoint.SetX( anX.toDouble( &isXOk ) );
 459     aPoint.SetY( anY.toDouble( &isYOk ) );
 460     aPoint.SetZ( aZ.toDouble( &isZOk ) );
 461
 462     if ( !isXOk || !isYOk || !isZOk )
 463       return false;
 464
 465     if ( boost::math::isnan( aPoint.X() ) || boost::math::isinf( aPoint.X() ) ||
 466          boost::math::isnan( aPoint.Y() ) || boost::math::isinf( aPoint.Y() ) ||
 467          boost::math::isnan( aPoint.Z() ) || boost::math::isinf( aPoint.Z() ) )
 468       return false;
 469
 470     // Invert the z value if requested
 471     if ( anIsAltitudesInverted )
 472       aPoint.SetZ( -aPoint.Z() );
 473
 474     thePoints.Append( aPoint );
 475   }
 476
 477 #ifdef _TIMER
 478   aTimer.Stop();
 479   std::ofstream stream( "W:/HYDRO/WORK/log.txt", std::ofstream::out );
 480   aTimer.Show( stream );
 481 #endif
 482
 483   return true;
 484 }
 485
 486
 487 bool HYDROData_Bathymetry::CreateBoundaryPolyline() const
 488 {
 489   Handle(HYDROData_Document) aDocument = HYDROData_Document::Document( myLab );
 490   Handle_HYDROData_PolylineXY aResult =
 491     Handle_HYDROData_PolylineXY::DownCast( aDocument->CreateObject( KIND_POLYLINEXY ) );
 492
 493   if( aResult.IsNull() )
 494     return false;
 495
 496   //search free name
 497   QString aPolylinePref = GetName() + "_Boundary";
 498   QString aPolylineName = HYDROData_Tool::GenerateObjectName( aDocument, aPolylinePref );
 499   aResult->SetName( aPolylineName );
 500
 501   double Xmin = 0.0, Xmax = 0.0, Ymin = 0.0, Ymax = 0.0;
 502   bool isFirst = true;
 503   AltitudePoints aPoints = GetAltitudePoints();
 504
 505   AltitudePoints::Iterator anIter( aPoints );
 506   for ( ; anIter.More(); anIter.Next() )
 507   {
 508     const AltitudePoint& aPoint = anIter.Value();
 509
 510     double x = aPoint.X(), y = aPoint.Y();
 511     if( isFirst || x<Xmin )
 512       Xmin = x;
 513     if( isFirst || x>Xmax )
 514       Xmax = x;
 515     if( isFirst || y<Ymin )
 516       Ymin = y;
 517     if( isFirst || y>Ymax )
 518       Ymax = y;
 519     isFirst = false;
 520   }
 521
 522   aResult->AddSection( "bound", HYDROData_IPolyline::SECTION_POLYLINE, true );
 523   aResult->AddPoint( 0, HYDROData_IPolyline::Point( Xmin, Ymin ) );
 524   aResult->AddPoint( 0, HYDROData_IPolyline::Point( Xmin, Ymax ) );
 525   aResult->AddPoint( 0, HYDROData_IPolyline::Point( Xmax, Ymax ) );
 526   aResult->AddPoint( 0, HYDROData_IPolyline::Point( Xmax, Ymin ) );
 527   aResult->Update();
 528
 529   return true;
 530 }