X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FHYDROData%2FHYDROData_DTM.cxx;h=207530f8df59e3a0c37bc80081e3315a3197bbcc;hb=bd81def133c235351f4c18834c6660450f16b9ee;hp=6678a8b90035bef1d3d4e70e2c7bec180cd0bc87;hpb=c34437a61318ce4c9aa56d4cbfeaea0e639dd2d3;p=modules%2Fhydro.git

diff --git a/src/HYDROData/HYDROData_DTM.cxx b/src/HYDROData/HYDROData_DTM.cxx
index 6678a8b9..207530f8 100644
--- a/src/HYDROData/HYDROData_DTM.cxx
+++ b/src/HYDROData/HYDROData_DTM.cxx
@@ -36,7 +36,6 @@
 #include <BRep_Builder.hxx>
 #include <GeomProjLib.hxx>
 #include <Geom_TrimmedCurve.hxx>
-#include <Geom_Plane.hxx>
 #include <BRepTools_WireExplorer.hxx>
 #include <TopTools_IndexedMapOfShape.hxx>
 #include <BRepBuilderAPI_MakeFace.hxx>
@@ -59,12 +58,12 @@
 #include <TopTools_ListIteratorOfListOfShape.hxx>
 #include <TopTools_SequenceOfShape.hxx>
 #include <assert.h>
+#include <float.h>
 
-IMPLEMENT_STANDARD_HANDLE( HYDROData_DTM, HYDROData_Bathymetry )
 IMPLEMENT_STANDARD_RTTIEXT( HYDROData_DTM, HYDROData_Bathymetry )
 
-HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv, const gp_Vec2d& theProfileDir, double theDeltaZ )
-  : myXcurv( theXCurv ), myProfileDir( theProfileDir ), myDeltaZ( theDeltaZ )
+HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv, const gp_Vec2d& theProfileDir, double theDeltaZ, double theMaxZ )
+  : myXcurv( theXCurv ), myProfileDir( theProfileDir ), myDeltaZ( theDeltaZ ), myMaxZ (theMaxZ)
 {
 }
 
@@ -87,16 +86,22 @@ double HYDROData_DTM::CurveUZ::DeltaZ() const
   return myDeltaZ;
 }
 
+double HYDROData_DTM::CurveUZ::MaxZ() const
+{
+  return myMaxZ;
+}
+
 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator + ( const CurveUZ& c ) const
 {
-  HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv(), ProfileDir() + c.ProfileDir(), DeltaZ() + c.DeltaZ() );
+  HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv(), ProfileDir() + c.ProfileDir(), DeltaZ() + c.DeltaZ(), MaxZ() + c.MaxZ() );
   size_t n = size(), n1 = c.size();
   if( n!=n1 )
   {
     std::cout << "Warning: different number of points in curves: " << n << ", " << n1 << std::endl;
   }
-  res.reserve( n );
-  for( int i=0; i<n; i++ )
+  int q = n < n1 ? n : n1;
+  res.reserve( q );
+  for( int i=0; i<q; i++ )
   {
     PointUZ p;
     p.U = operator[]( i ).U + c[i].U;
@@ -108,7 +113,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, DeltaZ()*d );
+  HYDROData_DTM::CurveUZ res( Xcurv()*d, ProfileDir()*d, DeltaZ()*d, MaxZ()*d );
   size_t n = size();
   res.reserve( n );
   for( int i=0; i<n; i++ )
@@ -186,7 +191,7 @@ void HYDROData_DTM::PointsToEdge(const AltitudePoints& pnts, TopoDS_Edge& E )
   anInterpolator.Perform() ;
   if (anInterpolator.IsDone()) 
   {
-    Handle(Geom_Curve) C = anInterpolator.Curve();
+    Handle(Geom_BSplineCurve) C = anInterpolator.Curve();
     E = BRepBuilderAPI_MakeEdge(C).Edge();
   }
 }
@@ -286,7 +291,7 @@ void HYDROData_DTM::CreateProfilesFromDTM (const HYDROData_SequenceOfObjects& In
   int aLower = InpProfiles.Lower(), anUpper = InpProfiles.Upper();
   size_t n = anUpper - aLower + 1;
 
-  std::vector<Handle_HYDROData_Profile> profiles;
+  std::vector<Handle(HYDROData_Profile)> profiles;
   profiles.reserve( n ); 
   for( int i=aLower; i<=anUpper; i++ )
   {
@@ -351,7 +356,7 @@ bool HYDROData_DTM::GetPlanarFaceFromBanks( const TopoDS_Edge& LB, const TopoDS_
   return !res;
 }
 
-void HYDROData_DTM::CreateProfiles(const std::vector<Handle_HYDROData_Profile>& theProfiles,
+void HYDROData_DTM::CreateProfiles(const std::vector<Handle(HYDROData_Profile)>& theProfiles,
                                    double theDDZ,
                                    double theSpatialStep,
                                    AltitudePoints& theOutLeft,
@@ -408,7 +413,7 @@ void HYDROData_DTM::CreateProfiles(const std::vector<Handle_HYDROData_Profile>&
 
 
 
-void HYDROData_DTM::GetProperties( const Handle_HYDROData_Profile& theProfile,
+void HYDROData_DTM::GetProperties( const Handle(HYDROData_Profile)& theProfile,
                     gp_Pnt& theLowestPoint, gp_Vec2d& theDir,
                     double& theZMin, double& theZMax )
 {
@@ -424,7 +429,7 @@ void HYDROData_DTM::GetProperties( const Handle_HYDROData_Profile& theProfile,
   HYDROData_Profile::ProfilePoints points = theProfile->GetProfilePoints();
   int lo = points.Lower();
   int up = points.Upper();
-  theZMin = std::numeric_limits<double>::max();
+  theZMin = DBL_MAX;
   theZMax = -theZMin;
   for( int i=lo; i<=up; i++ )
   {
@@ -476,7 +481,7 @@ Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve,
       aLast2d = To2D( aLastPnt, theTr, theUMin, theUMax );
 
     gp_Vec2d dir( aFirst2d, aLast2d );
-    Handle_Geom2d_Line aLine2d = new Geom2d_Line( aFirst2d, gp_Dir2d( dir.X(), dir.Y() ) );
+    Handle(Geom2d_Line) aLine2d = new Geom2d_Line( aFirst2d, gp_Dir2d( dir.X(), dir.Y() ) );
     return new Geom2d_TrimmedCurve( aLine2d, 0, aLast2d.Distance( aFirst2d ) );
   }
 
@@ -497,8 +502,8 @@ Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve,
 
 #include <GCE2d_MakeSegment.hxx>
 #include <Geom2dAPI_InterCurveCurve.hxx>
-bool IsCooriented( const Handle_HYDROData_Profile& theProfile1,
-                   const Handle_HYDROData_Profile& theProfile2 ) 
+bool IsCooriented( const Handle(HYDROData_Profile)& theProfile1,
+                   const Handle(HYDROData_Profile)& theProfile2 ) 
 {
   if( theProfile1==theProfile2 )
     return true;
@@ -513,15 +518,15 @@ bool IsCooriented( const Handle_HYDROData_Profile& theProfile1,
   GCE2d_MakeSegment s2(rp1, rp2);
 
   Geom2dAPI_InterCurveCurve inter;
-  inter.Init(s1, s2);
+  inter.Init(s1.Value(), s2.Value());
   if (inter.NbPoints() == 0)
     return true;
   else
     return false;
 }
 
-Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis( 
-  const std::vector<Handle_HYDROData_Profile>& theProfiles,
+Handle(Geom2d_BSplineCurve) HYDROData_DTM::CreateHydraulicAxis( 
+  const std::vector<Handle(HYDROData_Profile)>& theProfiles,
   std::vector<double>& theDistances )
 {
   size_t n = theProfiles.size();
@@ -538,8 +543,8 @@ Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis(
   theProfiles[0]->Update();
   for( size_t i = 1; i < n; i++ )
   {
-    Handle_HYDROData_Profile aProfile = theProfiles[i];
-    Handle_HYDROData_Profile aPrevProfile = theProfiles[i-1];
+    Handle(HYDROData_Profile) aProfile = theProfiles[i];
+    Handle(HYDROData_Profile) aPrevProfile = theProfiles[i-1];
 
     if( !IsCooriented( aProfile, aPrevProfile ) )
     {
@@ -555,8 +560,8 @@ Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis(
   // Stage 2. Calculate normals so that each normal "points" to the next profile
   for( size_t i = 0; i < n; i++ )
   {
-    Handle_HYDROData_Profile aProfile = theProfiles[i];
-    Handle_HYDROData_Profile aNextProfile = i==n-1 ? theProfiles[i-1] : theProfiles[i+1];
+    Handle(HYDROData_Profile) aProfile = theProfiles[i];
+    Handle(HYDROData_Profile) aNextProfile = i==n-1 ? theProfiles[i-1] : theProfiles[i+1];
 
     gp_Pnt aLowest;
     gp_Vec2d aNormal;
@@ -611,11 +616,11 @@ Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis(
   return aResult;
 }
 
-std::vector<Handle_Geom2d_Curve> HYDROData_DTM::ProfileToParametric( 
-  const Handle_HYDROData_Profile& theProfile,
+std::vector<Handle(Geom2d_Curve)> HYDROData_DTM::ProfileToParametric( 
+  const Handle(HYDROData_Profile)& theProfile,
   double& theUMin, double& theUMax, gp_Vec2d& theDir )
 {
-  std::vector<Handle_Geom2d_Curve> curves;
+  std::vector<Handle(Geom2d_Curve)> curves;
   
   // Transformation of the coordinate systems
   gp_Pnt aLowest;
@@ -652,7 +657,7 @@ std::vector<Handle_Geom2d_Curve> HYDROData_DTM::ProfileToParametric(
 
 bool CalcMidWidth( const std::set<double>& intersections, double& theMid, double& theWid )
 {
-  double umin = std::numeric_limits<double>::max(),
+  double umin = DBL_MAX,
          umax = -umin;
 
   size_t n = intersections.size();
@@ -673,20 +678,20 @@ bool CalcMidWidth( const std::set<double>& intersections, double& theMid, double
   return true;
 }
 
-void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& theProfile, 
-                                           double theXCurv, double theMinZ, double theMaxZ, double theDDZ,
+void HYDROData_DTM::ProfileDiscretization( const Handle(HYDROData_Profile)& theProfile, 
+                                           double theXCurv, double theMinZ, double theMaxZ, double theTopZ, double theDDZ,
                                            CurveUZ& theMidPointCurve,
                                            CurveUZ& theWidthCurve,                                           
                                            int& intersection_nb,
                                            double theTolerance)
 {
-  double aDblMax = std::numeric_limits<double>::max(),
+  double aDblMax = DBL_MAX,
          aUMin = aDblMax,
          aUMax = -aUMin,
          aVMax = 1000000;
   
   gp_Vec2d aProfileDir;
-  std::vector<Handle_Geom2d_Curve> curves = ProfileToParametric( theProfile, aUMin, aUMax, aProfileDir );
+  std::vector<Handle(Geom2d_Curve)> curves = ProfileToParametric( theProfile, aUMin, aUMax, aProfileDir );
   size_t n = curves.size();
 
   if( n==0 )
@@ -706,9 +711,9 @@ void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& thePr
   curves.push_back( aT2 );
   
   int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 );
-  theMidPointCurve = CurveUZ( theXCurv, aProfileDir, theMinZ );
+  theMidPointCurve = CurveUZ( theXCurv, aProfileDir, theMinZ, theTopZ);
   theMidPointCurve.reserve( psize );
-  theWidthCurve = CurveUZ( theXCurv, aProfileDir, theMinZ );
+  theWidthCurve = CurveUZ( theXCurv, aProfileDir, theMinZ, theTopZ );
   theWidthCurve.reserve( psize );
 
   n = curves.size();
@@ -719,7 +724,7 @@ void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& thePr
     std::set<double> intersections;
     for( size_t i = 0; i < n; i++ )
     {
-      Handle_Geom2d_Curve aCurve = curves[i];
+      Handle(Geom2d_Curve) aCurve = curves[i];
       Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
       for( int k=1, m=anIntersect.NbPoints(); k<=m; k++ )
         intersections.insert( anIntersect.Point( k ).X() );
@@ -765,7 +770,7 @@ void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCur
     theInterpolation.push_back( theCurveB );
 }
 #include <BRepLib_MakeEdge2d.hxx>
-void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
+void HYDROData_DTM::CurveTo3D( const Handle(Geom2d_BSplineCurve)& theHydraulicAxis,
                                const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
                                AltitudePoints& thePoints )
 {
@@ -782,7 +787,8 @@ void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxi
   
   size_t n = theMidCurve.size();
   std::map<double, AltitudePoint> sorted_points;
-  for( size_t i=0; i<n; i++ )
+  bool isOnTop = false;
+  for( size_t i=0; i<n; i++ ) // build the two banks of the interpolated profile, from bottom to top
   {
     double param1 = theMidCurve[i].U - theWidthCurve[i].U / 2;
     double param2 = theMidCurve[i].U + theWidthCurve[i].U / 2;
@@ -790,12 +796,24 @@ 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 + theMidCurve.DeltaZ();
-
-    AltitudePoint p3d_1( p1.X(), p1.Y(), z ), p3d_2( p2.X(), p2.Y(), z );
-
-    sorted_points[param1] = p3d_1;
-    sorted_points[param2] = p3d_2;
+    bool arrivedOnTop = false;
+    double z = 0;
+    if (theMidCurve[i].Z <= theMidCurve.MaxZ())
+      z = theMidCurve[i].Z + theMidCurve.DeltaZ();
+    else
+      {
+        z = theMidCurve.MaxZ() + theMidCurve.DeltaZ(); // limit z to linear interpolation between maxima on extremity profiles
+        arrivedOnTop = true; // do not keep points after this one
+      }
+    if (!isOnTop)
+      {
+        AltitudePoint p3d_1( p1.X(), p1.Y(), z ), p3d_2( p2.X(), p2.Y(), z );
+
+        sorted_points[param1] = p3d_1;
+        sorted_points[param2] = p3d_2;
+      }
+    //if (arrivedOnTop)
+    //  isOnTop =true; // do not keep points after this one (commented: leads to strange limits of 2D shape)
   }
 
   thePoints.reserve( sorted_points.size() );
@@ -806,29 +824,29 @@ void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxi
       thePoints.push_back( it->second );
 }
 
-inline double max( double a, double b )
-{
-  if( a>b )
-    return a;
-  else
-    return b;
-}
-
-inline double min( double a, double b )
-{
-  if( a<b )
-    return a;
-  else
-    return b;
-}
+//inline double max( double a, double b )
+//{
+//  if( a>b )
+//    return a;
+//  else
+//    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
-  ( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
-    const Handle_HYDROData_Profile& theProfileA,
+  ( const Handle(Geom2d_BSplineCurve)& theHydraulicAxis,
+    const Handle(HYDROData_Profile)& theProfileA,
     double theXCurvA,
-    const Handle_HYDROData_Profile& theProfileB,
+    const Handle(HYDROData_Profile)& theProfileB,
     double theXCurvB,
     double theDDZ, int theNbSteps, bool isAddSecond,
     int& inter_nb_1, int& inter_nb_2)
@@ -841,17 +859,17 @@ std::vector<HYDROData_Bathymetry::AltitudePoints> HYDROData_DTM::Interpolate
   GetProperties( theProfileB, lowestB, dirB, zminB, zmaxB ); 
 
   
-  double hmax = max( zmaxA-zminA, zmaxB-zminB );
+  double hmax = zmaxA-zminA > zmaxB-zminB ? zmaxA-zminA : zmaxB-zminB;
 
   //double dz = zminB - zminA;
   //double zmin = min( zminA, zminB );
   //double zmax = max( zmaxA, zmaxB );
 
-  CurveUZ midA(0, gp_Vec2d(), 0), midB(0, gp_Vec2d(), 0);
-  CurveUZ widA(0, gp_Vec2d(), 0), widB(0, gp_Vec2d(), 0);
+  CurveUZ midA(0, gp_Vec2d(), 0, 0), midB(0, gp_Vec2d(), 0, 0);
+  CurveUZ widA(0, gp_Vec2d(), 0, 0), widB(0, gp_Vec2d(), 0, 0);
 
-  ProfileDiscretization( theProfileA, theXCurvA, zminA, zminA+hmax, theDDZ, midA, widA, inter_nb_1 ); 
-  ProfileDiscretization( theProfileB, theXCurvB, zminB, zminB+hmax, theDDZ, midB, widB, inter_nb_2 );
+  ProfileDiscretization( theProfileA, theXCurvA, zminA, zminA+hmax, zmaxA-zminA, theDDZ, midA, widA, inter_nb_1 );
+  ProfileDiscretization( theProfileB, theXCurvB, zminB, zminB+hmax, zmaxB-zminB, theDDZ, midB, widB, inter_nb_2 );
 
   std::vector<CurveUZ> mid, wid;
   Interpolate( midA, midB, theNbSteps, mid, isAddSecond );
@@ -872,7 +890,7 @@ std::vector<HYDROData_Bathymetry::AltitudePoints> HYDROData_DTM::Interpolate
 }
 
 HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
-  ( const std::vector<Handle_HYDROData_Profile>& theProfiles,
+  ( const std::vector<Handle(HYDROData_Profile)>& theProfiles,
     double theDDZ, double theSpatialStep,
     AltitudePoints& theLeft,
     AltitudePoints& theRight,
@@ -885,7 +903,7 @@ HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
     return points;
 
   std::vector<double> distances;
-  Handle_Geom2d_BSplineCurve aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
+  Handle(Geom2d_BSplineCurve) aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
   if( aHydraulicAxis.IsNull() )
     return points;
 
@@ -940,7 +958,7 @@ HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
   return points;
 }
 
-int HYDROData_DTM::EstimateNbPoints( const std::vector<Handle_HYDROData_Profile>& theProfiles,
+int HYDROData_DTM::EstimateNbPoints( const std::vector<Handle(HYDROData_Profile)>& theProfiles,
                                      double theDDZ, double theSpatialStep )
 {
   size_t n = theProfiles.size();
@@ -950,7 +968,7 @@ int HYDROData_DTM::EstimateNbPoints( const std::vector<Handle_HYDROData_Profile>
     return 1 << 20;
 
   std::vector<double> distances;
-  Handle_Geom2d_BSplineCurve aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
+  Handle(Geom2d_BSplineCurve) aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
   if( aHydraulicAxis.IsNull() )
     return 0;