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 PYTHON_BATHYMETRY_ID "KIND_BATHYMETRY"
  19
  20
  21 IMPLEMENT_STANDARD_HANDLE(HYDROData_Bathymetry, HYDROData_IAltitudeObject)
  22 IMPLEMENT_STANDARD_RTTIEXT(HYDROData_Bathymetry, HYDROData_IAltitudeObject)
  23
  24 HYDROData_Bathymetry::HYDROData_Bathymetry()
  25 : HYDROData_IAltitudeObject()
  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   aResList << QString( "%1.SetName( \"%2\" );" )
  47               .arg( aBathymetryName ).arg( aBathymetryName );
  48
  49   QString aFilePath = GetFilePath();
  50   if ( !aFilePath.isEmpty() )
  51   {
  52     aResList << QString( "%1.ImportFromFile( \"%2\" );" )
  53                 .arg( aBathymetryName ).arg( aFilePath );
  54   }
  55   else
  56   {
  57     // TODO : bathymetry is composed from other bathymetry(ies)
  58   }
  59
  60   return aResList;
  61 }
  62
  63 void HYDROData_Bathymetry::SetAltitudePoints( const AltitudePoints& thePoints )
  64 {
  65   RemoveAltitudePoints();
  66
  67   if ( thePoints.isEmpty() )
  68     return;
  69
  70   // Save coordinates
  71   Handle(TDataStd_RealArray) aCoordsArray =
  72     TDataStd_RealArray::Set( myLab.FindChild( DataTag_AltitudePoints ), 0, thePoints.size() * 3 - 1 );
  73
  74   AltitudePoints::const_iterator aListItBeg =  thePoints.constBegin();
  75   AltitudePoints::const_iterator aListItEnd =  thePoints.constEnd();
  76   for ( int i = 0 ; aListItBeg != aListItEnd; ++i, ++aListItBeg )
  77   {
  78     const AltitudePoint& aPoint = *aListItBeg;
  79
  80     aCoordsArray->SetValue( i * 3, aPoint.X() );
  81     aCoordsArray->SetValue( i * 3 + 1, aPoint.Y() );
  82     aCoordsArray->SetValue( i * 3 + 2, aPoint.Z() );
  83   }
  84 }
  85
  86 HYDROData_Bathymetry::AltitudePoints HYDROData_Bathymetry::GetAltitudePoints() const
  87 {
  88   AltitudePoints aPoints;
  89
  90   Handle(TDataStd_RealArray) aCoordsArray;
  91   if ( !myLab.FindChild( DataTag_AltitudePoints ).FindAttribute( TDataStd_RealArray::GetID(), aCoordsArray ) )
  92     return aPoints;
  93
  94   int aLowerIdx = aCoordsArray->Lower();
  95   int anUpperIdx = aCoordsArray->Upper();
  96   for ( int i = aCoordsArray->Lower(), n = aCoordsArray->Upper(); i <= n; )
  97   {
  98     if ( i + 3 > n )
  99       break;
 100
 101     AltitudePoint aPoint;
 102     aPoint.SetX( aCoordsArray->Value( i++ ) );
 103     aPoint.SetY( aCoordsArray->Value( i++ ) );
 104     aPoint.SetZ( aCoordsArray->Value( i++ ) );
 105
 106     aPoints << aPoint;
 107   }
 108
 109   return aPoints;
 110 }
 111
 112 void HYDROData_Bathymetry::RemoveAltitudePoints()
 113 {
 114   TDF_Label aLab = myLab.FindChild( DataTag_AltitudePoints );
 115   aLab.ForgetAllAttributes();
 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 aResAltitude = GetInvalidAltitude();
 174
 175   AltitudePoints anAltitudePoints = GetAltitudePoints();
 176   if ( anAltitudePoints.isEmpty() )
 177     return aResAltitude;
 178
 179   QPolygonF aBoundingRect;
 180
 181   // Boundary plane
 182   // [ 0 (top-left) ]          [ 1 (top-right) ]
 183   //                  thePoint
 184   // [ 2 (bot-left) ]          [ 3 (bot-right) ]
 185   AltitudePoint aBounds[ 4 ] = { AltitudePoint( -DBL_MAX, -DBL_MAX, GetInvalidAltitude() ),
 186                                  AltitudePoint(  DBL_MAX, -DBL_MAX, GetInvalidAltitude() ),
 187                                  AltitudePoint( -DBL_MAX,  DBL_MAX, GetInvalidAltitude() ),
 188                                  AltitudePoint(  DBL_MAX,  DBL_MAX, GetInvalidAltitude() ) };
 189
 190   AltitudePoints::const_iterator aListItBeg = anAltitudePoints.constBegin();
 191   AltitudePoints::const_iterator aListItEnd = anAltitudePoints.constEnd();
 192   for ( ; aListItBeg != aListItEnd; ++aListItBeg )
 193   {
 194     const AltitudePoint& aPoint = *aListItBeg;
 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   // Check if requested point is inside of our bounding rectangle
 275   if ( !aBoundingRect.boundingRect().contains( thePoint.X(), thePoint.Y() ) )
 276     return aResAltitude;
 277
 278   // Calculate result altitude for point
 279   AltitudePoint aFirstPoint( aBounds[ 0 ] ), aSecPoint( aBounds[ 1 ] );
 280
 281   // At first we merge top and bottom borders
 282   if ( aBounds[ 0 ].Y() != aBounds[ 2 ].Y() || aBounds[ 0 ].X() != aBounds[ 2 ].X() )
 283     interpolateAltitudeForPoints( thePoint, aBounds[ 0 ], aBounds[ 2 ], aFirstPoint, true );
 284
 285   if ( aBounds[ 1 ].Y() != aBounds[ 3 ].Y() || aBounds[ 1 ].X() != aBounds[ 3 ].X() )
 286     interpolateAltitudeForPoints( thePoint, aBounds[ 1 ], aBounds[ 3 ], aSecPoint, true );
 287
 288   AltitudePoint aResPoint( aFirstPoint );
 289
 290   // At last we merge left and right borders
 291   if ( aFirstPoint.Y() != aSecPoint.Y() || aFirstPoint.X() != aSecPoint.X() )
 292     interpolateAltitudeForPoints( thePoint, aFirstPoint, aSecPoint, aResPoint, false );
 293
 294   aResAltitude = aResPoint.Z();
 295
 296   return aResAltitude;
 297 }
 298
 299 void HYDROData_Bathymetry::SetFilePath(const QString& theFilePath)
 300 {
 301   TCollection_AsciiString anAsciiStr( theFilePath.toStdString().c_str() );
 302   TDataStd_AsciiString::Set( myLab.FindChild( DataTag_FilePath ), anAsciiStr );
 303 }
 304
 305 QString HYDROData_Bathymetry::GetFilePath() const
 306 {
 307   QString aRes;
 308
 309   Handle(TDataStd_AsciiString) anAsciiStr;
 310   if ( myLab.FindChild( DataTag_FilePath ).FindAttribute( TDataStd_AsciiString::GetID(), anAsciiStr ) )
 311     aRes = QString( anAsciiStr->Get().ToCString() );
 312
 313   return aRes;
 314 }
 315
 316 bool HYDROData_Bathymetry::ImportFromFile( const QString& theFileName )
 317 {
 318   // Try to open the file
 319   QFile aFile( theFileName );
 320   if ( !aFile.exists() || !aFile.open( QIODevice::ReadOnly ) )
 321     return false;
 322
 323   bool aRes = false;
 324
 325   QString aFileSuf = QFileInfo( aFile ).suffix().toLower();
 326
 327   AltitudePoints aPoints;
 328
 329   // Try to import the file
 330   if ( aFileSuf == "xyz" )
 331     aRes = importFromXYZFile( aFile, aPoints );
 332
 333   // Close the file
 334   aFile.close();
 335
 336   if ( aRes )
 337   {
 338     // Update file path and altitude points of this Bathymetry
 339     SetFilePath( theFileName );
 340     SetAltitudePoints( aPoints );
 341   }
 342
 343   return aRes && !aPoints.isEmpty();
 344 }
 345
 346 bool HYDROData_Bathymetry::importFromXYZFile( QFile&          theFile,
 347                                               AltitudePoints& thePoints )
 348 {
 349   if ( !theFile.isOpen() )
 350     return false;
 351
 352   // Strings in file is written as:
 353   //  1. X(float) Y(float) Z(float)
 354   //  2. X(float) Y(float) Z(float)
 355   //  ...
 356
 357   while ( !theFile.atEnd() )
 358   {
 359     QString aLine = theFile.readLine().simplified();
 360     if ( aLine.isEmpty() )
 361       continue;
 362
 363     QStringList aValues = aLine.split( QRegExp( "\\s+" ), QString::SkipEmptyParts );
 364     if ( aValues.length() < 3 )
 365       return false;
 366
 367     AltitudePoint aPoint;
 368
 369     QString anX = aValues.value( 0 );
 370     QString anY = aValues.value( 1 );
 371     QString aZ  = aValues.value( 2 );
 372
 373     bool isXOk = false, isYOk = false, isZOk = false;
 374
 375     aPoint.SetX( anX.toDouble( &isXOk ) );
 376     aPoint.SetY( anY.toDouble( &isYOk ) );
 377     aPoint.SetZ(  aZ.toDouble( &isZOk ) );
 378
 379     if ( !isXOk || !isYOk || !isZOk )
 380       return false;
 381
 382     thePoints << aPoint;
 383   }
 384
 385   return true;
 386 }
 387
 388
 389
 390