src/HYDROData/HYDROData_Bathymetry.cxx

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