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