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