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