From: Paul RASCLE <paul.rascle@edf.fr>
Date: Wed, 5 Aug 2015 11:30:41 +0000 (+0200)
Subject: channel altitude interpolation
X-Git-Tag: BR_quadtree_20150925~9
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=ca1ea3d105fb2b5552e51746c87a1a09cf2a4caf;p=modules%2Fhydro.git

channel altitude interpolation
---

diff --git a/src/HYDROData/HYDROData_ChannelAltitude.cxx b/src/HYDROData/HYDROData_ChannelAltitude.cxx
index 1a7b515f..1b956fea 100644
--- a/src/HYDROData/HYDROData_ChannelAltitude.cxx
+++ b/src/HYDROData/HYDROData_ChannelAltitude.cxx
@@ -23,6 +23,8 @@
 #include "HYDROData_Channel.h"
 #include "HYDROData_Projection.h"
 #include "HYDROData_Polyline3D.h"
+#include "HYDROData_PolylineXY.h"
+#include "HYDROData_ProfileUZ.h"
 #include "HYDROData_Profile.h"
 
 #define _DEVDEBUG_
@@ -95,11 +97,11 @@ double HYDROData_ChannelAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) c
   TopoDS_Shape aShape = aGuideLine->GetShape3D();
   double middleZ = -9999;
   aGuideLine->GetMiddleZ(middleZ); // use the middle Z value of the 3d line to help the projection.
-  gp_Pnt P1(thePoint.X(), thePoint.Y(), middleZ);
   if (middleZ < -9000)
     {
       DEBTRACE("the middle Z value of the 3d line is incorrect");
     }
+  gp_Pnt P1(thePoint.X(), thePoint.Y(), middleZ);
   TopoDS_Shape aPoint = BRepBuilderAPI_MakeVertex(P1).Shape();
 
   if (aPoint.IsNull() || aShape.IsNull())
@@ -245,7 +247,45 @@ double HYDROData_ChannelAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) c
           const gp_Pnt &aPntProj = aDistShSh.PointOnShape2(i);
           TopoDS_Shape aProj = BRepBuilderAPI_MakeVertex(aPntProj).Shape();
           DEBTRACE("projection: (" << aPntProj.X() << ", " << aPntProj.Y() << ", " << aPntProj.Z() << ")");
-          return aPntProj.Z() + 2.;
+          gp_XY aProjXY = gp_XY(aPntProj.X(), aPntProj.Y());
+          aProjXY.Subtract(thePoint);
+          double distance = aProjXY.Modulus();
+          DEBTRACE("distance " << distance);
+          int i1 = 0;
+          gp_XY pt1 = gp_XY();
+          gp_XY pt2 = gp_XY();
+          HYDROData_ProfileUZ::PointsList aProfilePoints = aProfile->GetParametricPoints();
+          for ( int i = 1, aNbPoints = aProfilePoints.Size(); i <= aNbPoints; ++i )
+            {
+              const HYDROData_IPolyline::Point& aPolylinePoint = aProfilePoints.Value( i );
+              DEBTRACE("  profile point: " << aPolylinePoint.X() << " " << aPolylinePoint.Y());
+              if (aPolylinePoint.X() < distance)
+                {
+                  i1 = i;
+                  pt1 = aPolylinePoint;
+                }
+              if (aPolylinePoint.X() >= distance)
+                {
+                  pt2 = aPolylinePoint;
+                  break;
+                }
+            }
+          if ((i1 == 0) && (distance > 0))
+            {
+              DEBTRACE("Projection aborted : non centered profile");
+              return aResAltitude;
+            }
+          if (i1 == aProfilePoints.Size()) // distance >= profile width
+            {
+              aResAltitude = aPntProj.Z() + pt1.Y();
+            }
+          else
+            {
+              double z = pt1.Y() + (pt2.Y() - pt1.Y())*(distance -pt1.X())/(pt2.X()-pt1.X());
+              aResAltitude = aPntProj.Z() + z;
+            }
+          DEBTRACE("altitude: " << aResAltitude << " delta: " << aResAltitude - aPntProj.Z());
+          return aResAltitude;
         }
     }