presentation of DTM

diff --git a/src/HYDROData/HYDROData_DTM.cxx b/src/HYDROData/HYDROData_DTM.cxx
index 4db3a948b63f60fa52fc018b1ea6fa81c9ff54f3..579c09eb19386692d46980371ac2a602ff987db1 100644 (file)
--- a/src/HYDROData/HYDROData_DTM.cxx
+++ b/src/HYDROData/HYDROData_DTM.cxx
@@ -18,13 +18,17 @@
 #include <Geom2d_TrimmedCurve.hxx>
 #include <Geom_BSplineCurve.hxx>
 #include <Geom2d_BSplineCurve.hxx>
+#include <GeomAPI_Interpolate.hxx>
 #include <TColStd_Array1OfReal.hxx>
 #include <TColStd_Array1OfInteger.hxx>
 #include <TColgp_Array1OfPnt.hxx>
+#include <TColgp_Array1OfVec.hxx>
+#include <TColgp_HArray1OfPnt.hxx>
 #include <Geom2dAPI_InterCurveCurve.hxx>
 #include <Geom2dAPI_ProjectPointOnCurve.hxx>
 #include <Geom2dAdaptor_Curve.hxx>
 #include <GCPnts_AbscissaPoint.hxx>
+#include <BRepBuilderAPI_MakeEdge.hxx>
 #include <limits>
 
 IMPLEMENT_STANDARD_HANDLE( HYDROData_DTM, HYDROData_Bathymetry )
@@ -75,6 +79,58 @@ HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator * ( double d ) const
   return res;
 }
 
+void HYDROData_DTM::Bank::reserve( int theNbPoints )
+{
+  myPoints.reserve( theNbPoints );
+  myDirs.reserve( theNbPoints );
+}
+
+void HYDROData_DTM::Bank::push_back( const gp_Pnt& thePnt, const gp_Dir& theTangent )
+{
+  myPoints.push_back( thePnt );
+  myDirs.push_back( theTangent );
+}
+
+void HYDROData_DTM::Bank::clear()
+{
+  myPoints.clear();
+  myDirs.clear();
+}
+
+TopoDS_Edge HYDROData_DTM::Bank::createEdge( bool is3d ) const
+{
+  size_t n = myPoints.size();
+  Handle_Geom_BSplineCurve aCurve;
+
+  Handle(TColgp_HArray1OfPnt) points = new TColgp_HArray1OfPnt( 1, (int)n );
+  TColgp_Array1OfVec tangents( 1, (int)n );
+  Handle(TColStd_HArray1OfBoolean) flags = new TColStd_HArray1OfBoolean( 1, (int)n );
+
+  for( size_t i = 1; i <= n; i++ )
+  {
+    gp_Pnt aPnt = myPoints[i-1];
+    if( !is3d )
+      aPnt.SetZ( 0.0 );
+    gp_Vec aVec = myDirs[i-1];
+    if( !is3d )
+      aVec.SetZ( 0.0 );
+    points->SetValue( (int)i, aPnt );
+    tangents.SetValue( (int)i, aVec );
+    flags->SetValue( (int)i, Standard_True );
+  }
+
+  GeomAPI_Interpolate anInterpolator( points, Standard_False, Standard_False );
+  anInterpolator.Load( tangents, flags );
+  anInterpolator.Perform();
+  if( anInterpolator.IsDone() )
+  {
+    aCurve = anInterpolator.Curve();
+    return BRepBuilderAPI_MakeEdge( aCurve ).Edge();
+  }
+  else
+    return TopoDS_Edge();
+}
+
 
 
 HYDROData_DTM::HYDROData_DTM()
@@ -134,8 +190,11 @@ void HYDROData_DTM::Update()
   double step = GetSpatialStep();
   const double EPS = 1E-3;
   AltitudePoints points;
+  
+  myLeft.clear();
+  myRight.clear();
   if( ddz>EPS && step>EPS )
-    points = Interpolate( profiles, ddz, step );
+    points = Interpolate( profiles, ddz, step, &myLeft, &myRight );
   SetAltitudePoints( points );
 }
 
@@ -281,7 +340,7 @@ Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis(
     Standard_Real aParamFirst = anAdaptor.FirstParameter(), aParamLast = anAdaptor.LastParameter();
     for( size_t i = 1; i <= n; i++ )
     {
-      gp_Pnt2d aPnt = points->Value( i );
+      gp_Pnt2d aPnt = points->Value( (Standard_Integer)i );
       Geom2dAPI_ProjectPointOnCurve aProject( aPnt, aResult );
       Standard_Real aParam = aProject.LowerDistanceParameter();
       double aDistance = GCPnts_AbscissaPoint::Length( anAdaptor, aParamFirst, aParam );  
@@ -451,7 +510,8 @@ void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCur
 
 void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
                                const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
-                               AltitudePoints& thePoints )
+                               AltitudePoints& thePoints,
+                               Bank* theLeft, Bank* theRight, double dz )
 {
   Geom2dAdaptor_Curve anAdaptor( theHydraulicAxis );
   GCPnts_AbscissaPoint ap( anAdaptor, theMidCurve.Xcurv(), anAdaptor.FirstParameter() );  
@@ -462,8 +522,12 @@ void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxi
   anAdaptor.D1( aParam, point, tangent );
   profile_dir.SetCoord( tangent.Y(), -tangent.X() );
   profile_dir.Normalize();
+  gp_Dir tangent_n( tangent.X(), tangent.Y(), dz );
 
   size_t n = theMidCurve.size();
+  double min_param = 1E+15;
+  double max_param = -1E+15;
+  double z1, z2;
   for( size_t i=0; i<n; i++ )
   {
     double param1 = theMidCurve[i].U - theWidthCurve[i].U / 2;
@@ -474,10 +538,44 @@ void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxi
 
     double z = theMidCurve[i].Z;
 
+    if( param1 < min_param )
+    {
+      min_param = param1;
+      z1 = z;
+    }
+    if( param2 < min_param )
+    {
+      min_param = param2;
+      z1 = z;
+    }
+    if( param1 > max_param )
+    {
+      max_param = param1;
+      z2 = z;
+    }
+    if( param2 > max_param )
+    {
+      max_param = param2;
+      z2 = z;
+    }
+
     AltitudePoint p3d_1( p1.X(), p1.Y(), z ), p3d_2( p2.X(), p2.Y(), z );
     thePoints.push_back( p3d_1 );
     thePoints.push_back( p3d_2 );
   }
+
+  if( theLeft )
+  {
+    gp_Pnt2d left2d = point.Translated( min_param * profile_dir / 2 );
+    gp_Pnt left( left2d.X(), left2d.Y(), z1 );
+    theLeft->push_back( left, tangent_n );
+  }
+  if( theRight )
+  {
+    gp_Pnt2d right2d = point.Translated( max_param * profile_dir / 2 );
+    gp_Pnt right( right2d.X(), right2d.Y(), z2 );
+    theRight->push_back( right, tangent_n );
+  }
 }
 
 inline double max( double a, double b )
@@ -494,7 +592,8 @@ HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
     double theXCurvA,
     const Handle_HYDROData_Profile& theProfileB,
     double theXCurvB,
-    double theDDZ, int theNbSteps, bool isAddSecond )
+    double theDDZ, int theNbSteps, bool isAddSecond,
+    Bank* theLeft, Bank* theRight )
 {
   double zminA, zmaxA, zminB, zmaxB;
   gp_Pnt lowestA, lowestB;
@@ -503,6 +602,8 @@ 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 zmin = max( zminA, zminB );
   double zmax = max( zmaxA, zmaxB );
 
@@ -521,13 +622,16 @@ HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
   AltitudePoints points;
   points.reserve( p*q );
   for( size_t i=0; i<p; i++ )
-    CurveTo3D( theHydraulicAxis, mid[i], wid[i], points );
+    CurveTo3D( theHydraulicAxis, mid[i], wid[i], points, theLeft, theRight, dz );
+
   return points;
 }
 
 HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
   ( const std::vector<Handle_HYDROData_Profile>& theProfiles,
-    double theDDZ, double theSpatialStep )
+    double theDDZ, double theSpatialStep,
+    Bank* theLeft,
+    Bank* theRight )
 {
   AltitudePoints points;
   size_t n = theProfiles.size();
@@ -539,15 +643,29 @@ HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
   if( aHydraulicAxis.IsNull() )
     return points;
 
+  int aNbStepsComplete = 0;
+  for( size_t i=0, n1=n-1; i<n1; i++ )
+  {
+    double aDistance = distances[i+1]-distances[i];
+    aNbStepsComplete += ( int(aDistance/theSpatialStep) + 1 );
+  }
+  if( theLeft )
+    theLeft->reserve( aNbStepsComplete );
+  if( theRight )
+    theRight->reserve( aNbStepsComplete );
+  
   for( size_t i=0, n1=n-1; i<n1; i++ )
   {
     double aDistance = distances[i+1]-distances[i];
     int aNbSteps = int(aDistance/theSpatialStep) + 1;
     bool isAddSecond = i==n1-1;
+
     AltitudePoints local_points = Interpolate( aHydraulicAxis, theProfiles[i], distances[i], 
-      theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond );
+      theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond, theLeft, theRight );
+
     if( i==0 )
       points.reserve( local_points.size() * ( n-1 ) );
+
     for( size_t j=0, m=local_points.size(); j<m; j++ )
       points.push_back( local_points[j] );
   }

diff --git a/src/HYDROData/HYDROData_DTM.h b/src/HYDROData/HYDROData_DTM.h
index 87bb10affd899d042820ad12d51cb983d316a0e8..717a9ed588850a5c21719710444d86b466335110 100644 (file)
--- a/src/HYDROData/HYDROData_DTM.h
+++ b/src/HYDROData/HYDROData_DTM.h
@@ -21,12 +21,14 @@
 
 #include "HYDROData_Bathymetry.h"
 #include <vector>
+#include <gp_Pnt2d.hxx>
 
 class Handle_HYDROData_Profile;
 class Handle_Geom2d_BSplineCurve;
 class Handle_Geom2d_Curve;
 class gp_Pnt;
 class gp_Vec2d;
+class TopoDS_Edge;
 
 DEFINE_STANDARD_HANDLE( HYDROData_DTM, HYDROData_Bathymetry )
 
@@ -82,6 +84,18 @@ public:
   private:
     double myXcurv;
   };
+  class Bank
+  {
+  public:
+    void reserve( int theNbPoints );
+    void push_back( const gp_Pnt& thePnt, const gp_Dir& theTangent );
+    void clear();
+    TopoDS_Edge createEdge( bool is3d ) const;
+
+  private:
+    std::vector<gp_Pnt> myPoints;
+    std::vector<gp_Dir> myDirs;
+  };
 
 protected:
   friend class HYDROData_Iterator;
@@ -110,7 +124,10 @@ protected:
 
   static void CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
                          const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
-                         AltitudePoints& thePoints );
+                         AltitudePoints& thePoints,
+                         Bank* theLeftBank = 0,
+                         Bank* theRightBank = 0,
+                         double dz = 0 );
   
   static void Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCurveB, 
                            int theNbSteps, std::vector<CurveUZ>& theInterpolation,
@@ -121,10 +138,18 @@ protected:
                                      double theXCurvA,
                                      const Handle_HYDROData_Profile& theProfileB,
                                      double theXCurvB,
-                                     double theDDZ, int theNbSteps, bool isAddSecond );
+                                     double theDDZ, int theNbSteps, bool isAddSecond,
+                                     Bank* theLeftBank = 0,
+                                     Bank* theRightBank = 0 );
 
   static AltitudePoints Interpolate( const std::vector<Handle_HYDROData_Profile>& theProfiles,
-                                     double theDDZ, double theSpatialStep );
+                                     double theDDZ, double theSpatialStep,
+                                     Bank* theLeftBank = 0,
+                                     Bank* theRightBank = 0 );
+
+private:
+  Bank myLeft;
+  Bank myRight;
 };
 
 #endif

diff --git a/src/HYDRO_tests/test_HYDROData_Stream.cxx b/src/HYDRO_tests/test_HYDROData_Stream.cxx
index d354d17584fe1f9f16cdddad6eae3f6847f19b37..1466db52b9a919c46a285a6f591a942a67184c78 100644 (file)
--- a/src/HYDRO_tests/test_HYDROData_Stream.cxx
+++ b/src/HYDRO_tests/test_HYDROData_Stream.cxx
@@ -23,6 +23,7 @@
 #include <HYDROData_Stream.h>
 #include <HYDROData_IPolyline.h>
 #include <HYDROData_Profile.h>
+#include <HYDROData_PolylineXY.h>
 #include <TestViewer.h>
 #include <QApplication>
 
@@ -182,5 +183,39 @@ void test_HYDROData_Stream::test_dump()
   CPPUNIT_ASSERT_EQUAL( std::string( "stream_1.Update()" ),                                aScript1[8].toStdString() );
   CPPUNIT_ASSERT_EQUAL( std::string( "" ),                                                 aScript1[9].toStdString() );
 
+  // Case 2. With hydraulic axis
+
+  Handle(HYDROData_Stream) aStream2 = 
+    Handle(HYDROData_Stream)::DownCast( aDoc->CreateObject( KIND_STREAM ) );
+  aStream2->SetName( "stream_2" );
+
+  Handle(HYDROData_PolylineXY) anHAxis = 
+    Handle(HYDROData_PolylineXY)::DownCast( aDoc->CreateObject( KIND_POLYLINEXY ) );
+  anHAxis->SetName( "axis" );
+
+  aStream2->SetProfiles( profiles, false );
+  aStream2->SetDDZ( 0.2 );
+  aStream2->SetSpatialStep( 3.0 );
+  aStream2->SetReferenceObject( anHAxis, HYDROData_Stream::DataTag_HydraulicAxis );
+
+  objs.clear();
+  objs["p1"] = aProfile1;
+  objs["p2"] = aProfile2;
+  objs["axis"] = anHAxis;
+
+  QStringList aScript2 = aStream2->DumpToPython( "", objs );
+  CPPUNIT_ASSERT_EQUAL( 11, aScript2.size() );
+  CPPUNIT_ASSERT_EQUAL( std::string( "stream_2 = hydro_doc.CreateObject( KIND_STREAM )" ), aScript2[0].toStdString() );
+  CPPUNIT_ASSERT_EQUAL( std::string( "stream_2.SetName( \"stream_2\" )" ),                 aScript2[1].toStdString() );
+  CPPUNIT_ASSERT_EQUAL( std::string( "" ),                                                 aScript2[2].toStdString() );
+  CPPUNIT_ASSERT_EQUAL( std::string( "stream_2.SetHydraulicAxis( axis )" ),                aScript2[3].toStdString() );
+  CPPUNIT_ASSERT_EQUAL( std::string( "stream_2.AddProfile( p1 )" ),                        aScript2[4].toStdString() );
+  CPPUNIT_ASSERT_EQUAL( std::string( "stream_2.AddProfile( p2 )" ),                        aScript2[5].toStdString() );
+  CPPUNIT_ASSERT_EQUAL( std::string( "stream_2.SetDDZ( 0.200 )" ),                         aScript2[6].toStdString() );
+  CPPUNIT_ASSERT_EQUAL( std::string( "stream_2.SetSpatialStep( 3.000 )" ),                 aScript2[7].toStdString() );
+  CPPUNIT_ASSERT_EQUAL( std::string( "" ),                                                 aScript2[8].toStdString() );
+  CPPUNIT_ASSERT_EQUAL( std::string( "stream_2.Update()" ),                                aScript2[9].toStdString() );
+  CPPUNIT_ASSERT_EQUAL( std::string( "" ),                                                 aScript2[10].toStdString() );
+
   aDoc->Close();
 }