debug of DTM creation on different cases

diff --git a/src/HYDROData/HYDROData_DTM.cxx b/src/HYDROData/HYDROData_DTM.cxx
index 16f82e4c78d2bbcc2b49e3a902cf0f4cde674cd2..e2bbbdeaed916028e1f55157a1d2422810ed6887 100644 (file)
--- a/src/HYDROData/HYDROData_DTM.cxx
+++ b/src/HYDROData/HYDROData_DTM.cxx
@@ -53,8 +53,8 @@ IMPLEMENT_STANDARD_HANDLE( HYDROData_DTM, HYDROData_Bathymetry )
 IMPLEMENT_STANDARD_RTTIEXT( HYDROData_DTM, HYDROData_Bathymetry )
 
 
-HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv, const gp_Vec2d& theProfileDir )
-  : myXcurv( theXCurv ), myProfileDir( theProfileDir )
+HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv, const gp_Vec2d& theProfileDir, double theDeltaZ )
+  : myXcurv( theXCurv ), myProfileDir( theProfileDir ), myDeltaZ( theDeltaZ )
 {
 }
 
@@ -72,9 +72,14 @@ gp_Vec2d HYDROData_DTM::CurveUZ::ProfileDir() const
   return myProfileDir;
 }
 
+double HYDROData_DTM::CurveUZ::DeltaZ() const
+{
+  return myDeltaZ;
+}
+
 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator + ( const CurveUZ& c ) const
 {
-  HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv(), ProfileDir() + c.ProfileDir() );
+  HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv(), ProfileDir() + c.ProfileDir(), DeltaZ() + c.DeltaZ() );
   size_t n = size();
   res.reserve( n );
   for( int i=0; i<n; i++ )
@@ -89,7 +94,7 @@ HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator + ( const CurveUZ& c ) c
 
 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator * ( double d ) const
 {
-  HYDROData_DTM::CurveUZ res( Xcurv()*d, ProfileDir()*d );
+  HYDROData_DTM::CurveUZ res( Xcurv()*d, ProfileDir()*d, DeltaZ()*d );
   size_t n = size();
   res.reserve( n );
   for( int i=0; i<n; i++ )
@@ -541,9 +546,8 @@ std::vector<Handle_Geom2d_Curve> HYDROData_DTM::ProfileToParametric(
   double zmin, zmax;
   GetProperties( theProfile, aLowest, theDir, false, zmin, zmax );
 
-  aLowest.SetZ( 0 );
   gp_Ax3 aStd3d( gp_Pnt( 0, 0, 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( 1, 0, 0 ) );
-  gp_Ax3 aLocal( aLowest, gp_Dir( 0, 0, 1 ), gp_Dir( theDir.X(), theDir.Y(), 0 ) );
+  gp_Ax3 aLocal( gp_Pnt( aLowest.X(), aLowest.Y(), 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( theDir.X(), theDir.Y(), 0 ) );
 
   gp_Trsf aTransf;
   aTransf.SetTransformation( aStd3d, aLocal );
@@ -570,18 +574,19 @@ std::vector<Handle_Geom2d_Curve> HYDROData_DTM::ProfileToParametric(
 }
 
 
-bool CalcMidWidth( const std::vector<gp_Pnt2d>& intersections, double& theMid, double& theWid )
+bool CalcMidWidth( const std::set<double>& intersections, double& theMid, double& theWid )
 {
   double umin = std::numeric_limits<double>::max(),
          umax = -umin;
 
   size_t n = intersections.size();
-  if( n <= 1 )
+  if( n <= 0 )
     return false;
 
-  for( size_t i = 0; i < n; i++ )
+  std::set<double>::const_iterator it = intersections.begin(), last = intersections.end();
+  for( ; it!=last; it++ )
   {
-    double u = intersections[i].X();
+    double u = *it;
     if( u<umin )
       umin = u;
     if( u>umax )
@@ -625,32 +630,33 @@ void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& thePr
   curves.push_back( aT2 );
   
   int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 );
-  theMidPointCurve = CurveUZ( theXCurv, aProfileDir );
+  theMidPointCurve = CurveUZ( theXCurv, aProfileDir, theMinZ );
   theMidPointCurve.reserve( psize );
-  theWidthCurve = CurveUZ( theXCurv, aProfileDir );
+  theWidthCurve = CurveUZ( theXCurv, aProfileDir, theMinZ );
   theWidthCurve.reserve( psize );
 
   n = curves.size();
   // for each discrete value of z we search intersection with profile
-  for( double z = theMinZ; z <= theMaxZ; z += theDDZ )
+  for( double z1 = theMinZ; z1 <= theMaxZ; z1 += theDDZ )
   {
-    Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z ), gp_Dir2d( 1, 0 ) );
-    std::vector<gp_Pnt2d> intersections;
+    Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z1 ), gp_Dir2d( 1, 0 ) );
+    std::set<double> intersections;
     for( size_t i = 0; i < n; i++ )
     {
       Handle_Geom2d_Curve aCurve = curves[i];
       Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
       for( int k=1, m=anIntersect.NbPoints(); k<=m; k++ )
-        intersections.push_back( anIntersect.Point( k ) );
+        intersections.insert( anIntersect.Point( k ).X() );
     }
 
     intersection_nb = intersections.size();
-    if( intersection_nb >= 2 )
+    if( intersection_nb >= 1 )
     {
       double u_mid, u_wid;
       if( !CalcMidWidth( intersections, u_mid, u_wid ) )
         continue;
 
+      double z = z1 - theMinZ;
       PointUZ p_mid;
       p_mid.U = u_mid;
       p_mid.Z = z;
@@ -685,7 +691,7 @@ void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCur
 #include <BRepLib_MakeEdge2d.hxx>
 void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
                                const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
-                               AltitudePoints& thePoints, double dz )
+                               AltitudePoints& thePoints )
 {
   Geom2dAdaptor_Curve anAdaptor( theHydraulicAxis );
   TopoDS_Edge E2d = BRepLib_MakeEdge2d(theHydraulicAxis).Edge();
@@ -695,7 +701,7 @@ void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxi
   gp_Pnt2d point;
   anAdaptor.D0( aParam, point );
   gp_Vec2d profile_dir = theMidCurve.ProfileDir();
-  gp_Dir tangent_n( -profile_dir.Y(), profile_dir.X(), dz );
+  //gp_Dir tangent_n( -profile_dir.Y(), profile_dir.X(), dz );
   profile_dir.Normalize();
   
   size_t n = theMidCurve.size();
@@ -708,7 +714,7 @@ void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxi
     gp_Pnt2d p1 = point.Translated( param1 * profile_dir);
     gp_Pnt2d p2 = point.Translated( param2 * profile_dir);
 
-    double z = theMidCurve[i].Z;
+    double z = theMidCurve[i].Z + theMidCurve.DeltaZ();
 
     AltitudePoint p3d_1( p1.X(), p1.Y(), z ), p3d_2( p2.X(), p2.Y(), z );
 
@@ -730,6 +736,14 @@ inline double max( double a, double b )
     return b;
 }
 
+inline double min( double a, double b )
+{
+  if( a<b )
+    return a;
+  else
+    return b;
+}
+
 #include <BRepLib_MakeWire.hxx>
 
 std::vector<HYDROData_Bathymetry::AltitudePoints> HYDROData_DTM::Interpolate
@@ -748,16 +762,18 @@ std::vector<HYDROData_Bathymetry::AltitudePoints> HYDROData_DTM::Interpolate
   GetProperties( theProfileA, lowestA, dirA, false, zminA, zmaxA ); 
   GetProperties( theProfileB, lowestB, dirB, false, zminB, zmaxB ); 
 
-  double dz = zminB - zminA;
+  
+  double hmax = max( zmaxA-zminA, zmaxB-zminB );
 
-  double zmin = max( zminA, zminB );
-  double zmax = max( zmaxA, zmaxB );
+  //double dz = zminB - zminA;
+  //double zmin = min( zminA, zminB );
+  //double zmax = max( zmaxA, zmaxB );
 
-  CurveUZ midA(0, gp_Vec2d()), midB(0, gp_Vec2d());
-  CurveUZ widA(0, gp_Vec2d()), widB(0, gp_Vec2d());
+  CurveUZ midA(0, gp_Vec2d(), 0), midB(0, gp_Vec2d(), 0);
+  CurveUZ widA(0, gp_Vec2d(), 0), widB(0, gp_Vec2d(), 0);
 
-  ProfileDiscretization( theProfileA, theXCurvA, zmin, zmax, theDDZ, midA, widA, inter_nb_1 ); 
-  ProfileDiscretization( theProfileB, theXCurvB, zmin, zmax, theDDZ, midB, widB, inter_nb_2 );
+  ProfileDiscretization( theProfileA, theXCurvA, zminA, zminA+hmax, theDDZ, midA, widA, inter_nb_1 ); 
+  ProfileDiscretization( theProfileB, theXCurvB, zminB, zminB+hmax, theDDZ, midB, widB, inter_nb_2 );
 
   std::vector<CurveUZ> mid, wid;
   Interpolate( midA, midB, theNbSteps, mid, isAddSecond );
@@ -771,7 +787,7 @@ std::vector<HYDROData_Bathymetry::AltitudePoints> HYDROData_DTM::Interpolate
   for( size_t i=0; i<p; i++ )
   {
     points[i].reserve( q );
-    CurveTo3D( theHydraulicAxis, mid[i], wid[i], points[i], dz );
+    CurveTo3D( theHydraulicAxis, mid[i], wid[i], points[i] );
   }
 
   return points;

diff --git a/src/HYDROData/HYDROData_DTM.h b/src/HYDROData/HYDROData_DTM.h
index 9d1cd477ccd397e4c0d98cb74bce33274449c8dc..e10b286b61d5b7c8f4980281ddc820e9e28ab04a 100644 (file)
--- a/src/HYDROData/HYDROData_DTM.h
+++ b/src/HYDROData/HYDROData_DTM.h
@@ -85,11 +85,12 @@ public:
   class CurveUZ : public std::vector<PointUZ>
   {
   public:
-    CurveUZ( double theXcurv, const gp_Vec2d& theProfileDir = gp_Vec2d() );
+    CurveUZ( double theXcurv, const gp_Vec2d& theProfileDir, double theDeltaZ );
     ~CurveUZ();
 
     double Xcurv() const;
-    gp_Vec2d ProfileDir() const; 
+    gp_Vec2d ProfileDir() const;
+    double DeltaZ() const;
 
     CurveUZ operator + ( const CurveUZ& ) const;
     CurveUZ operator * ( double ) const;
@@ -97,6 +98,7 @@ public:
   private:
     double myXcurv;
     gp_Vec2d myProfileDir;
+    double myDeltaZ;
   };
 
 protected:
@@ -129,7 +131,7 @@ protected:
 
   static void CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
                          const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
-                         AltitudePoints& thePoints, double dz );
+                         AltitudePoints& thePoints );
   
   static void Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCurveB, 
                            int theNbSteps, std::vector<CurveUZ>& theInterpolation,

diff --git a/src/HYDRO_tests/test_HYDROData_DTM.cxx b/src/HYDRO_tests/test_HYDROData_DTM.cxx
index 153f3c44292eda3e1e7cfec5875706b422ef2577..2050a73bc93a2a5e4e0a1f20b3e09825fd450c17 100644 (file)
--- a/src/HYDRO_tests/test_HYDROData_DTM.cxx
+++ b/src/HYDRO_tests/test_HYDROData_DTM.cxx
@@ -338,7 +338,7 @@ void test_HYDROData_DTM::test_profile_discretization_polyline()
   aProfile->SetLeftPoint( gp_XY( 10, 10 ) );
   aProfile->SetRightPoint( gp_XY( 20, 20 ) );
 
-  HYDROData_DTM::CurveUZ aMid( 0.0 ), aWid( 0.0 );
+  HYDROData_DTM::CurveUZ aMid( 0.0, gp_Vec2d(), 0 ), aWid( 0.0, gp_Vec2d(), 0 );
   int dummy = 0;
   HYDROData_DTM::ProfileDiscretization( aProfile, 0.0, 0.0, 5.0, 0.5, aMid, aWid, dummy );
   CPPUNIT_ASSERT_EQUAL( 11, (int)aMid.size() );
@@ -376,7 +376,7 @@ Handle(HYDROData_Document) aDoc = HYDROData_Document::Document(1);
   aProfile->SetLeftPoint( gp_XY( 10, 10 ) );
   aProfile->SetRightPoint( gp_XY( 20, 20 ) );
 
-  HYDROData_DTM::CurveUZ aMid( 0.0 ), aWid( 0.0 );
+  HYDROData_DTM::CurveUZ aMid( 0.0, gp_Vec2d(), 0 ), aWid( 0.0, gp_Vec2d(), 0 );
   int dummy = 0 ;
   HYDROData_DTM::ProfileDiscretization( aProfile, 0.0, 0.0, 5.0, 0.5, aMid, aWid, dummy );
   CPPUNIT_ASSERT_EQUAL( 11, (int)aMid.size() );
@@ -431,7 +431,7 @@ void operator << ( std::ostream& s, const HYDROData_DTM::CurveUZ& c )
 
 void test_HYDROData_DTM::test_curves_interpolation()
 {
-  HYDROData_DTM::CurveUZ A(1.0), B(2.0);
+  HYDROData_DTM::CurveUZ A(1.0, gp_Vec2d(), 0), B(2.0, gp_Vec2d(), 0);
   A.push_back( HYDROData_DTM::PointUZ( 0, 0 ) );
   A.push_back( HYDROData_DTM::PointUZ( 1, 1 ) );
   A.push_back( HYDROData_DTM::PointUZ( 2, 2 ) );
@@ -498,13 +498,13 @@ void test_HYDROData_DTM::test_curve_to_3d()
 
   Handle_Geom2d_BSplineCurve HA = HYDROData_DTM::CreateHydraulicAxis( profiles, distances );
   HYDROData_DTM::AltitudePoints points;
-  HYDROData_DTM::CurveUZ mid( 5.0, gp_Vec2d(-10,10) );
+  HYDROData_DTM::CurveUZ mid( 5.0, gp_Vec2d(-10,10), 0 );
   mid.push_back( HYDROData_DTM::PointUZ( 0, 5 ) );
   mid.push_back( HYDROData_DTM::PointUZ( 1, 6 ) );
-  HYDROData_DTM::CurveUZ wid( 5.0, gp_Vec2d(-10,10) );
+  HYDROData_DTM::CurveUZ wid( 5.0, gp_Vec2d(-10,10), 0 );
   wid.push_back( HYDROData_DTM::PointUZ( 2, 5 ) );
   wid.push_back( HYDROData_DTM::PointUZ( 6, 6 ) );
-  HYDROData_DTM::CurveTo3D( HA, mid, wid, points, 0.0 );
+  HYDROData_DTM::CurveTo3D( HA, mid, wid, points );
 
   CPPUNIT_ASSERT_EQUAL( 4, (int)points.size() );
   CPPUNIT_ASSERT_EQUAL( HYDROData_DTM::AltitudePoint( 16.380, -2.186, 6.0 ), points[0] );