Calculation of the point's altitude from bathymetry (feature 9).

diff --git a/src/HYDROData/HYDROData_Bathymetry.cxx b/src/HYDROData/HYDROData_Bathymetry.cxx
index 31bd8cc7435de7c709beda555bbf18e51073eec5..afbb8d0460dee0cd24be5ed88fd06f0e98fda951 100644 (file)
--- a/src/HYDROData/HYDROData_Bathymetry.cxx
+++ b/src/HYDROData/HYDROData_Bathymetry.cxx
@@ -1,5 +1,6 @@
 
 #include "HYDROData_Bathymetry.h"
+#include "HYDROData_Tool.h"
 
 #include <gp_XY.hxx>
 #include <gp_XYZ.hxx>
@@ -9,8 +10,12 @@
 #include <QFile>
 #include <QFileInfo>
 #include <QPointF>
+#include <QPolygonF>
 #include <QStringList>
 
+#define INVALID_ALTITUDE_VALUE -9999.0
+
+
 IMPLEMENT_STANDARD_HANDLE(HYDROData_Bathymetry, HYDROData_Object)
 IMPLEMENT_STANDARD_RTTIEXT(HYDROData_Bathymetry, HYDROData_Object)
 
@@ -77,6 +82,59 @@ void HYDROData_Bathymetry::RemoveAltitudePoints()
   aLab.ForgetAllAttributes();
 }
 
+void interpolateAltitudeForPoints( const gp_XY&                               thePoint,
+                                   const HYDROData_Bathymetry::AltitudePoint& theFirstPoint,
+                                   const HYDROData_Bathymetry::AltitudePoint& theSecPoint,
+                                   HYDROData_Bathymetry::AltitudePoint&       theResPoint,
+                                   const bool&                                theIsVertical )
+{
+  double aCoordX = thePoint.X();
+  double aCoordY = thePoint.Y();
+
+  if ( theIsVertical )
+  {
+    aCoordX = theFirstPoint.X();
+
+    if ( !ValuesEquals( theFirstPoint.X(), theSecPoint.X() ) )
+    {
+      // Recalculate X coordinate by equation of line from two points
+      aCoordX = ( ( ( thePoint.Y() - theFirstPoint.Y() ) * ( theSecPoint.X() - theFirstPoint.X() ) ) /
+                  ( theSecPoint.Y() - theFirstPoint.Y() ) ) + theFirstPoint.X();
+    }
+  }
+  else
+  {
+    aCoordY = theFirstPoint.Y();
+
+    if ( !ValuesEquals( theFirstPoint.Y(), theSecPoint.Y() ) )
+    {
+      // Recalculate y by equation of line from two points
+      aCoordY = ( ( ( thePoint.X() - theFirstPoint.X() ) * ( theSecPoint.Y() - theFirstPoint.Y() ) ) /
+                  ( theSecPoint.X() - theFirstPoint.X() ) ) + theFirstPoint.Y();
+    }
+  }
+
+  theResPoint.SetX( aCoordX );
+  theResPoint.SetY( aCoordY );
+
+  // Calculate coefficient for interpolation
+  double aLength = Sqrt( Pow( theSecPoint.Y() - theFirstPoint.Y(), 2 ) +
+                         Pow( theSecPoint.X() - theFirstPoint.X(), 2 ) );
+
+  double aInterCoeff = 0;
+  if ( aLength != 0 )
+   aInterCoeff = ( theSecPoint.Z() - theFirstPoint.Z() ) / aLength;
+
+
+  double aNewLength = Sqrt( Pow( theResPoint.Y() - theFirstPoint.Y(), 2 ) +
+                            Pow( theResPoint.X() - theFirstPoint.X(), 2 ) );
+
+  // Calculate interpolated value
+  double aResVal = theFirstPoint.Z() + aInterCoeff * aNewLength;
+
+  theResPoint.SetZ( aResVal );
+}
+
 double HYDROData_Bathymetry::GetAltitudeForPoint( const QPointF& thePoint ) const
 {
   gp_XY aGpPoint( thePoint.x(), thePoint.y() );
@@ -88,8 +146,125 @@ double HYDROData_Bathymetry::GetAltitudeForPoint( const gp_XY& thePoint ) const
   double aResAltitude = -9999.90;
   
   AltitudePoints anAltitudePoints = GetAltitudePoints();
+  if ( anAltitudePoints.isEmpty() )
+    return aResAltitude;
+
+  QPolygonF aBoundingRect;
+
+  // Boundary plane
+  // [ 0 (top-left) ]          [ 1 (top-right) ]
+  //                  thePoint
+  // [ 2 (bot-left) ]          [ 3 (bot-right) ] 
+  AltitudePoint aBounds[ 4 ] = { AltitudePoint( -DBL_MAX, -DBL_MAX, INVALID_ALTITUDE_VALUE ),
+                                 AltitudePoint(  DBL_MAX, -DBL_MAX, INVALID_ALTITUDE_VALUE ),
+                                 AltitudePoint( -DBL_MAX,  DBL_MAX, INVALID_ALTITUDE_VALUE ),
+                                 AltitudePoint(  DBL_MAX,  DBL_MAX, INVALID_ALTITUDE_VALUE ) }; 
+
+  AltitudePoints::const_iterator aListItBeg = anAltitudePoints.constBegin();
+  AltitudePoints::const_iterator aListItEnd = anAltitudePoints.constEnd();
+  for ( ; aListItBeg != aListItEnd; ++aListItBeg )
+  {
+    const AltitudePoint& aPoint = *aListItBeg;
+
+    double aDeltaX = Abs( aPoint.X() ) - Abs( thePoint.X() );
+    double aDeltaY = Abs( aPoint.Y() ) - Abs( thePoint.Y() );
+
+    if ( ValuesEquals( aDeltaX, 0.0 ) ) // Both left and right sides
+    {
+      if ( ValuesEquals( aDeltaY, 0.0 ) ) // Both top and bottom sides
+      {
+        aResAltitude = aPoint.Z();
+        return aResAltitude;
+      }
+      else if ( aDeltaY < 0 ) // top side
+      {
+        // top border
+        if ( ValuesMoreEquals( aPoint.X(), aBounds[ 0 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 0 ].Y() ) )
+          aBounds[ 0 ] = aPoint;
+        if ( ValuesLessEquals( aPoint.X(), aBounds[ 1 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 1 ].Y() ) )
+          aBounds[ 1 ] = aPoint;
+      }
+      else
+      {
+        // bottom border
+        if ( ValuesMoreEquals( aPoint.X(), aBounds[ 2 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 2 ].Y() ) )
+          aBounds[ 2 ] = aPoint;
+        if ( ValuesLessEquals( aPoint.X(), aBounds[ 3 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 3 ].Y() ) )
+          aBounds[ 3 ] = aPoint;
+      }
+    }
+    else if ( aDeltaX < 0 ) // left side
+    {
+      if ( ValuesEquals( aDeltaY, 0.0 ) )
+      {
+        // Left border
+        if ( ValuesMoreEquals( aPoint.X(), aBounds[ 0 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 0 ].Y() ) )
+          aBounds[ 0 ] = aPoint;
+        if ( ValuesMoreEquals( aPoint.X(), aBounds[ 2 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 2 ].Y() ) )
+          aBounds[ 2 ] = aPoint;
+      }
+      else if ( aDeltaY < 0 )
+      {
+        // top left corner
+        if ( ValuesMoreEquals( aPoint.X(), aBounds[ 0 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 0 ].Y() ) )
+          aBounds[ 0 ] = aPoint;
+      }
+      else
+      {
+        // bottom left corner
+        if ( ValuesMoreEquals( aPoint.X(), aBounds[ 2 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 2 ].Y() ) )
+          aBounds[ 2 ] = aPoint;
+      }
+    }
+    else // right side
+    {
+      if ( ValuesEquals( aDeltaY, 0.0 ) )
+      {
+        // Right border
+        if ( ValuesLessEquals( aPoint.X(), aBounds[ 1 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 1 ].Y() ) )
+          aBounds[ 1 ] = aPoint;
+        if ( ValuesLessEquals( aPoint.X(), aBounds[ 3 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 3 ].Y() ) )
+          aBounds[ 3 ] = aPoint;
+      }
+      else if ( aDeltaY < 0 )
+      {
+        // top right corner
+        if ( ValuesLessEquals( aPoint.X(), aBounds[ 1 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 1 ].Y() ) )
+          aBounds[ 1 ] = aPoint;
+      }
+      else
+      {
+        // bottom right corner
+        if ( ValuesLessEquals( aPoint.X(), aBounds[ 3 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 3 ].Y() ) )
+          aBounds[ 3 ] = aPoint;
+      }
+    }
+
+    // Update bounding rectangle of our global grid
+    aBoundingRect << QPointF( aPoint.X(), aPoint.Y() );
+  }
 
-  // TODO : implement
+  // Check if requested point is inside of our bounding rectangle
+  if ( !aBoundingRect.boundingRect().contains( thePoint.X(), thePoint.Y() ) )
+    return INVALID_ALTITUDE_VALUE;
+
+  // Calculate result altitude for point
+  AltitudePoint aFirstPoint( aBounds[ 0 ] ), aSecPoint( aBounds[ 1 ] );
+
+  // At first we merge top and bottom borders
+  if ( aBounds[ 0 ].Y() != aBounds[ 2 ].Y() || aBounds[ 0 ].X() != aBounds[ 2 ].X() )
+    interpolateAltitudeForPoints( thePoint, aBounds[ 0 ], aBounds[ 2 ], aFirstPoint, true );
+
+  if ( aBounds[ 1 ].Y() != aBounds[ 3 ].Y() || aBounds[ 1 ].X() != aBounds[ 3 ].X() )
+    interpolateAltitudeForPoints( thePoint, aBounds[ 1 ], aBounds[ 3 ], aSecPoint, true );
+
+  AltitudePoint aResPoint( aFirstPoint );
+
+  // At last we merge left and right borders
+  if ( aFirstPoint.Y() != aSecPoint.Y() || aFirstPoint.X() != aSecPoint.X() )
+    interpolateAltitudeForPoints( thePoint, aFirstPoint, aSecPoint, aResPoint, false );
+    
+  aResAltitude = aResPoint.Z();
 
   return aResAltitude;
 }
@@ -136,7 +311,7 @@ bool HYDROData_Bathymetry::importFromXYZFile( QFile&          theFile,
 
   while ( !theFile.atEnd() )
   {
-    QString aLine = theFile.readLine();
+    QString aLine = theFile.readLine().simplified();
     if ( aLine.isEmpty() )
       continue;