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