Salome HOME
ImportFromFile()
[modules/hydro.git] / src / HYDROData / HYDROData_DTM.cxx
index 384f672a1348247ef936ac23b5f9523756fb6397..4db3a948b63f60fa52fc018b1ea6fa81c9ff54f3 100644 (file)
@@ -22,6 +22,9 @@
 #include <TColStd_Array1OfInteger.hxx>
 #include <TColgp_Array1OfPnt.hxx>
 #include <Geom2dAPI_InterCurveCurve.hxx>
+#include <Geom2dAPI_ProjectPointOnCurve.hxx>
+#include <Geom2dAdaptor_Curve.hxx>
+#include <GCPnts_AbscissaPoint.hxx>
 #include <limits>
 
 IMPLEMENT_STANDARD_HANDLE( HYDROData_DTM, HYDROData_Bathymetry )
@@ -82,8 +85,70 @@ HYDROData_DTM::~HYDROData_DTM()
 {
 }
 
+HYDROData_SequenceOfObjects HYDROData_DTM::GetProfiles() const
+{
+  return GetReferenceObjects( DataTag_Profiles );
+}
+
+void HYDROData_DTM::SetProfiles( const HYDROData_SequenceOfObjects& theProfiles )
+{
+  SetReferenceObjects( theProfiles, DataTag_Profiles );
+}
+
+double HYDROData_DTM::GetDDZ() const
+{
+  return GetDouble( DataTag_DDZ );
+}
+
+void HYDROData_DTM::SetDDZ( double theDDZ )
+{
+  SetDouble( DataTag_DDZ, theDDZ );
+}
+
+double HYDROData_DTM::GetSpatialStep() const
+{
+  return GetDouble( DataTag_SpatialStep );
+}
+
+void HYDROData_DTM::SetSpatialStep( double theSpatialStep )
+{
+  SetDouble( DataTag_SpatialStep, theSpatialStep );
+}
+
+void HYDROData_DTM::Update()
+{
+  HYDROData_SequenceOfObjects objs = GetProfiles();
+  int aLower = objs.Lower(), anUpper = objs.Upper();
+  size_t n = anUpper-aLower+1;
+
+  std::vector<Handle_HYDROData_Profile> profiles;
+  profiles.reserve( n );
+  for( int i=aLower; i<=anUpper; i++ )
+  {
+    Handle(HYDROData_Profile) aProfile = Handle(HYDROData_Profile)::DownCast( objs.Value( i ) );
+    if( !aProfile.IsNull() )
+      profiles.push_back( aProfile );
+  }
+
+  double ddz = GetDDZ();
+  double step = GetSpatialStep();
+  const double EPS = 1E-3;
+  AltitudePoints points;
+  if( ddz>EPS && step>EPS )
+    points = Interpolate( profiles, ddz, step );
+  SetAltitudePoints( points );
+}
+
 
-void GetProperties( const Handle_HYDROData_Profile& theProfile,
+
+
+
+
+
+
+
+
+void HYDROData_DTM::GetProperties( const Handle_HYDROData_Profile& theProfile,
                     gp_Pnt& theLowestPoint, gp_Vec2d& theDir,
                     bool isNormalDir,
                     double& theZMin, double& theZMax )
@@ -99,6 +164,20 @@ void GetProperties( const Handle_HYDROData_Profile& theProfile,
     theDir = gp_Vec2d( -y, x );
   else
     theDir = gp_Vec2d( x, y );
+
+  HYDROData_Profile::ProfilePoints points = theProfile->GetProfilePoints();
+  int lo = points.Lower();
+  int up = points.Upper();
+  theZMin = std::numeric_limits<double>::max();
+  theZMax = -theZMin;
+  for( int i=lo; i<=up; i++ )
+  {
+    double z = points.Value( i ).Z();
+    if( z>theZMax )
+      theZMax = z;
+    if( z<theZMin )
+      theZMin = z;
+  }
 }
 
 inline gp_Pnt2d To2D( const gp_Pnt& thePnt, const gp_Trsf& theTr,
@@ -193,7 +272,21 @@ Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis(
   if( anInterpolator.IsDone() )
   {
     aResult = anInterpolator.Curve();
-    //TODO: fill the distances vector
+
+    //fill the distances vector
+    Geom2dAdaptor_Curve anAdaptor( aResult );
+
+    theDistances.clear();
+    theDistances.reserve( n );
+    Standard_Real aParamFirst = anAdaptor.FirstParameter(), aParamLast = anAdaptor.LastParameter();
+    for( size_t i = 1; i <= n; i++ )
+    {
+      gp_Pnt2d aPnt = points->Value( i );
+      Geom2dAPI_ProjectPointOnCurve aProject( aPnt, aResult );
+      Standard_Real aParam = aProject.LowerDistanceParameter();
+      double aDistance = GCPnts_AbscissaPoint::Length( anAdaptor, aParamFirst, aParam );  
+      theDistances.push_back( aDistance );
+    }
   }
   return aResult;
 }
@@ -278,7 +371,8 @@ void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& thePr
 {
   double aDblMax = std::numeric_limits<double>::max(),
          aUMin = aDblMax,
-         aUMax = -aUMin;
+         aUMax = -aUMin,
+         aVMax = 1000000;
   
   std::vector<Handle_Geom2d_Curve> curves = ProfileToParametric( theProfile, aUMin, aUMax );
   size_t n = curves.size();
@@ -291,10 +385,10 @@ void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& thePr
   curves[0]->D0( curves[0]->FirstParameter(), aFirst );
   curves[n-1]->D0( curves[n-1]->LastParameter(), aLast );
   Handle(Geom2d_Line) aV1 = new Geom2d_Line( aFirst, gp_Dir2d( 0, 1 ) );
-  Handle(Geom2d_TrimmedCurve) aT1 = new Geom2d_TrimmedCurve( aV1, 0.0, aDblMax );
+  Handle(Geom2d_TrimmedCurve) aT1 = new Geom2d_TrimmedCurve( aV1, 0.0, aVMax );
   
   Handle(Geom2d_Line) aV2 = new Geom2d_Line( aLast, gp_Dir2d( 0, 1 ) );
-  Handle(Geom2d_TrimmedCurve) aT2 = new Geom2d_TrimmedCurve( aV2, 0.0, aDblMax );
+  Handle(Geom2d_TrimmedCurve) aT2 = new Geom2d_TrimmedCurve( aV2, 0.0, aVMax );
 
   curves.push_back( aT1 );
   curves.push_back( aT2 );
@@ -315,7 +409,7 @@ void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& thePr
     {
       Handle_Geom2d_Curve aCurve = curves[i];
       Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
-      for( int k=0, m=anIntersect.NbPoints(); k<m; k++ )
+      for( int k=1, m=anIntersect.NbPoints(); k<=m; k++ )
         intersections.push_back( anIntersect.Point( k ) );
     }
 
@@ -355,10 +449,35 @@ void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCur
     theInterpolation.push_back( theCurveB );
 }
 
-void HYDROData_DTM::CurveTo3d( const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
+void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
+                               const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
                                AltitudePoints& thePoints )
 {
-  //TODO
+  Geom2dAdaptor_Curve anAdaptor( theHydraulicAxis );
+  GCPnts_AbscissaPoint ap( anAdaptor, theMidCurve.Xcurv(), anAdaptor.FirstParameter() );  
+  double aParam = ap.Parameter();
+
+  gp_Pnt2d point;
+  gp_Vec2d tangent, profile_dir;
+  anAdaptor.D1( aParam, point, tangent );
+  profile_dir.SetCoord( tangent.Y(), -tangent.X() );
+  profile_dir.Normalize();
+
+  size_t n = theMidCurve.size();
+  for( size_t i=0; i<n; i++ )
+  {
+    double param1 = theMidCurve[i].U - theWidthCurve[i].U / 2;
+    double param2 = theMidCurve[i].U + theWidthCurve[i].U / 2;
+
+    gp_Pnt2d p1 = point.Translated( param1 * profile_dir / 2 );
+    gp_Pnt2d p2 = point.Translated( param2 * profile_dir / 2 );
+
+    double z = theMidCurve[i].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 );
+  }
 }
 
 inline double max( double a, double b )
@@ -369,12 +488,13 @@ inline double max( double a, double b )
     return b;
 }
 
-void HYDROData_DTM::Interpolate( const Handle_HYDROData_Profile& theProfileA,
-                                 double theXCurvA,
-                                 const Handle_HYDROData_Profile& theProfileB,
-                                 double theXCurvB,
-                                 double theDDZ, int theNbSteps, bool isAddSecond,
-                                 AltitudePoints& thePoints )
+HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
+  ( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
+    const Handle_HYDROData_Profile& theProfileA,
+    double theXCurvA,
+    const Handle_HYDROData_Profile& theProfileB,
+    double theXCurvB,
+    double theDDZ, int theNbSteps, bool isAddSecond )
 {
   double zminA, zmaxA, zminB, zmaxB;
   gp_Pnt lowestA, lowestB;
@@ -396,29 +516,40 @@ void HYDROData_DTM::Interpolate( const Handle_HYDROData_Profile& theProfileA,
   Interpolate( midA, midB, theNbSteps, mid, isAddSecond );
   Interpolate( widA, widB, theNbSteps, wid, isAddSecond );
 
-  size_t q = mid.size();
-  for( size_t i=0; i<q; i++ )
-    CurveTo3d( mid[i], wid[i], thePoints );
+  size_t p = mid.size();
+  size_t q = p>0 ? 2*mid[0].size() : 1;
+  AltitudePoints points;
+  points.reserve( p*q );
+  for( size_t i=0; i<p; i++ )
+    CurveTo3D( theHydraulicAxis, mid[i], wid[i], points );
+  return points;
 }
 
-void HYDROData_DTM::Interpolate( const std::vector<Handle_HYDROData_Profile>& theProfiles,
-                                 double theDDZ, double theSpatialStep, AltitudePoints& thePoints )
+HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
+  ( const std::vector<Handle_HYDROData_Profile>& theProfiles,
+    double theDDZ, double theSpatialStep )
 {
+  AltitudePoints points;
   size_t n = theProfiles.size();
   if( n<=1 )
-    return;
+    return points;
 
   std::vector<double> distances;
   Handle_Geom2d_BSplineCurve aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
   if( aHydraulicAxis.IsNull() )
-    return;
+    return points;
 
   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;
-    Interpolate( theProfiles[i], distances[i], 
-      theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond, thePoints );
+    AltitudePoints local_points = Interpolate( aHydraulicAxis, theProfiles[i], distances[i], 
+      theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond );
+    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] );
   }
+  return points;
 }