src/HYDROData/HYDROData_Bathymetry.cxx

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