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