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