correction interpolation Z pour les streams: détection des deux profils encadrant...

diff --git a/src/HYDROData/HYDROData_StreamAltitude.cxx b/src/HYDROData/HYDROData_StreamAltitude.cxx
index 96b33ad8cb62e0ad1fb8b67114081823b35191f1..4811380bd3c7c82f15aecbb85ee4db4fa344f9c8 100644 (file)
--- a/src/HYDROData/HYDROData_StreamAltitude.cxx
+++ b/src/HYDROData/HYDROData_StreamAltitude.cxx
@@ -51,6 +51,7 @@
 #include <Geom_Line.hxx>
 
 #include <QStringList>
+#include <cmath>
 
 #define _DEVDEBUG_
 #include "HYDRO_trace.hxx"
@@ -105,6 +106,48 @@ bool IsPointBetweenEdges( const gp_Pnt& aFirstPnt1, const gp_Pnt& aLastPnt1,
   return aPointState != TopAbs_OUT;
 }
 
+/** compute edge1^(firstPnt1,point) and edge2^(firstPnt2,point):
+ *  if the z values are opposite, the point is between the edges.
+ *  We must also check if the point is near the edges
+ *  (inside the circle defined by the edges: approximation)
+ *  to discriminate false positives in sinuous cases
+ */
+bool IsPointBetweenEdges2( const gp_Pnt& aFirstPnt1, const gp_Pnt& aLastPnt1,
+                           const gp_Pnt& aFirstPnt2, const gp_Pnt& aLastPnt2,
+                           const gp_Pnt& thePoint) {
+  double x1 = aLastPnt1.X() - aFirstPnt1.X(); // v1
+  double y1 = aLastPnt1.Y() - aFirstPnt1.Y();
+  double x2 = aLastPnt2.X() - aFirstPnt2.X(); // v2
+  double y2 = aLastPnt2.Y() - aFirstPnt2.Y();
+  double xa = thePoint.X() - aFirstPnt1.X();  // va
+  double ya = thePoint.Y() - aFirstPnt1.Y();
+  double xb = thePoint.X() - aFirstPnt2.X();  // vb
+  double yb = thePoint.Y() - aFirstPnt2.Y();
+  double z1 = x1*ya -xa*y1; // v1^va: z component
+  double z2 = x2*yb -xb*y2; // v2^vb: z component
+  bool isBetween = true;
+  if (((z1<0) and (z2<0)) or ((z1>=0) and (z2>=0)))
+    {
+      isBetween = false;
+    }
+  if (isBetween)
+    {
+      double xg = (aFirstPnt1.X() + aLastPnt1.X() + aFirstPnt2.X() + aLastPnt2.X())/4;
+      double yg = (aFirstPnt1.Y() + aLastPnt1.Y() + aFirstPnt2.Y() + aLastPnt2.Y())/4;
+      double df1 = (aFirstPnt1.X()-xg)*(aFirstPnt1.X()-xg) + (aFirstPnt1.Y()-yg)*(aFirstPnt1.Y()-yg);
+      double dl1 = (aLastPnt1.X()-xg)*(aLastPnt1.X()-xg) + (aLastPnt1.Y()-yg)*(aLastPnt1.Y()-yg);
+      double df2 = (aFirstPnt2.X()-xg)*(aFirstPnt2.X()-xg) + (aFirstPnt2.Y()-yg)*(aFirstPnt2.Y()-yg);
+      double dl2 = (aLastPnt2.X()-xg)*(aLastPnt2.X()-xg) + (aLastPnt2.Y()-yg)*(aLastPnt2.Y()-yg);
+      double r2 = fmax(df1,dl1);
+      r2 = fmax(r2,df2);
+      r2 = fmax(r2,dl2);
+      double d2 = (thePoint.X()-xg)*(thePoint.X()-xg) + (thePoint.Y()-yg)*(thePoint.Y()-yg);
+      if (d2 > r2)
+        isBetween = false;
+    }
+  return isBetween;
+}
+
 Standard_Real getAltitudeFromWire( const TopoDS_Wire& theWire,
                                    const Standard_Real& theLeftDist,
                                    const Standard_Real& theRightDist )
@@ -174,8 +217,9 @@ bool HYDROData_StreamAltitude::getBoundaryWiresForPoint(
 
       if ( !aPrevProfile.IsNull() )
       {        
-        if ( IsPointBetweenEdges( aPrevPnt1, aPrevPnt2, aPnt1, aPnt2, aTestPnt ) )
+        if ( IsPointBetweenEdges2( aPrevPnt1, aPrevPnt2, aPnt1, aPnt2, aTestPnt ) )
         {
+          DEBTRACE("point is between profiles "<< i-1 << " and " <<i);
           theLeftWire = TopoDS::Wire( aPrevProfile->GetShape3D() );
           theRightWire = TopoDS::Wire( aProfile->GetShape3D() );
           break;
@@ -221,7 +265,7 @@ bool HYDROData_StreamAltitude::getBoundaryWiresForPoint(
       aPnt2.SetZ( 0 );
 
       if ( !aPrevWire.IsNull() ) {
-        if ( IsPointBetweenEdges( aPrevPnt1, aPrevPnt2, aPnt1, aPnt2, aTestPnt ) ) {
+        if ( IsPointBetweenEdges2( aPrevPnt1, aPrevPnt2, aPnt1, aPnt2, aTestPnt ) ) {
           theLeftWire = aPrevWire;
           theRightWire = aWire;
           break;
@@ -246,20 +290,20 @@ double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) co
     Handle(HYDROData_Object)::DownCast( GetFatherObject() );
   if ( anObject.IsNull() )
   {
-    DEBTRACE("aStream.IsNull()");
+    DEBTRACE(" --- aStream.IsNull()");
     return aResAltitude;
   }
   TopoDS_Shape aTopShape = anObject->GetTopShape();
   if ( aTopShape.IsNull() )
   {
-    DEBTRACE("aTopShape.IsNull()");
+    DEBTRACE(" --- aTopShape.IsNull()");
     return aResAltitude;
   }
 
   TopExp_Explorer aFaceExp( aTopShape, TopAbs_FACE );
   if ( !aFaceExp.More() )
   {
-    DEBTRACE("!aFaceExp.More()");
+    DEBTRACE(" --- !aFaceExp.More()");
     return aResAltitude;
   }
 
@@ -282,7 +326,7 @@ double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) co
 #endif
   if ( aPointState == TopAbs_OUT )
   {
-       DEBTRACE("aPointState == TopAbs_OUT");
+       DEBTRACE(" --- aPointState == TopAbs_OUT");
     return aResAltitude;
   }
 
@@ -308,7 +352,7 @@ double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) co
   TopoDS_Wire aLeftWire, aRightWire;
   if ( !getBoundaryWiresForPoint( thePoint, aLeftWire, aRightWire ) )
   {
-      DEBTRACE("!getBoundaryProfilesForPoint( thePoint, aLeftProfile, aRightProfile )");
+      DEBTRACE(" --- !getBoundaryProfilesForPoint( thePoint, aLeftProfile, aRightProfile )");
       return aResAltitude;
   }
 
@@ -325,7 +369,10 @@ double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) co
   int aNbPoints1 = aRightProject.NbPoints();
   int aNbPoints2 = aLeftProject.NbPoints();
   if ( aNbPoints1 < 1 || aNbPoints2 < 1)
-     return aResAltitude;
+    {
+      DEBTRACE(" --- projections " << aNbPoints1 << " " << aNbPoints2);
+      return aResAltitude;
+    }
   Standard_Real aLeftDist = aLeftProject.LowerDistance();
   Standard_Real aRightDist = aRightProject.LowerDistance();
 
@@ -345,7 +392,7 @@ double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) co
   aRightProfileP2 = BRep_Tool::Pnt( aLastVertex );
 
   // Interpolate altitudes
-  // Left profile line ( the segment between the firts and the last profile point )
+  // Left profile line ( the segment between the first and the last profile point )
   aLeftProfileP1.SetZ( 0 );
   aLeftProfileP2.SetZ( 0 );
   gp_Vec aLeftProfileVec( aLeftProfileP1, aLeftProfileP2 );