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