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