src/HYDROData/HYDROData_DTM.cxx

   1
   2 #include <HYDROData_DTM.h>
   3 #include <HYDROData_Profile.h>
   4
   5 #include <Geom2d_BSplineCurve.hxx>
   6 #include <Geom2dAPI_Interpolate.hxx>
   7 #include <TColgp_HArray1OfPnt2d.hxx>
   8 #include <TColgp_Array1OfVec2d.hxx>
   9 #include <TColStd_HArray1OfBoolean.hxx>
  10 #include <TopoDS.hxx>
  11 #include <TopoDS_Edge.hxx>
  12 #include <TopoDS_Wire.hxx>
  13 #include <TopExp_Explorer.hxx>
  14 #include <BRep_Tool.hxx>
  15 #include <gp_Ax3.hxx>
  16 #include <Geom_Line.hxx>
  17 #include <Geom2d_Line.hxx>
  18 #include <Geom2d_TrimmedCurve.hxx>
  19 #include <Geom_BSplineCurve.hxx>
  20 #include <Geom2d_BSplineCurve.hxx>
  21 #include <GeomAPI_Interpolate.hxx>
  22 #include <TColStd_Array1OfReal.hxx>
  23 #include <TColStd_Array1OfInteger.hxx>
  24 #include <TColgp_Array1OfPnt.hxx>
  25 #include <TColgp_Array1OfVec.hxx>
  26 #include <TColgp_HArray1OfPnt.hxx>
  27 #include <Geom2dAPI_InterCurveCurve.hxx>
  28 #include <Geom2dAPI_ProjectPointOnCurve.hxx>
  29 #include <Geom2dAdaptor_Curve.hxx>
  30 #include <GCPnts_AbscissaPoint.hxx>
  31 #include <BRepBuilderAPI_MakeEdge.hxx>
  32 #include <limits>
  33
  34 IMPLEMENT_STANDARD_HANDLE( HYDROData_DTM, HYDROData_Bathymetry )
  35 IMPLEMENT_STANDARD_RTTIEXT( HYDROData_DTM, HYDROData_Bathymetry )
  36
  37
  38 HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv, const gp_Vec2d& theProfileDir )
  39   : myXcurv( theXCurv ), myProfileDir( theProfileDir )
  40 {
  41 }
  42
  43 HYDROData_DTM::CurveUZ::~CurveUZ()
  44 {
  45 }
  46
  47 double HYDROData_DTM::CurveUZ::Xcurv() const
  48 {
  49   return myXcurv;
  50 }
  51
  52 gp_Vec2d HYDROData_DTM::CurveUZ::ProfileDir() const
  53 {
  54   return myProfileDir;
  55 }
  56
  57 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator + ( const CurveUZ& c ) const
  58 {
  59   HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv(), ProfileDir() + c.ProfileDir() );
  60   size_t n = size();
  61   res.reserve( n );
  62   for( int i=0; i<n; i++ )
  63   {
  64     PointUZ p;
  65     p.U = operator[]( i ).U + c[i].U;
  66     p.Z = operator[]( i ).Z;
  67     res.push_back( p );
  68   }
  69   return res;
  70 }
  71
  72 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator * ( double d ) const
  73 {
  74   HYDROData_DTM::CurveUZ res( Xcurv()*d, ProfileDir()*d );
  75   size_t n = size();
  76   res.reserve( n );
  77   for( int i=0; i<n; i++ )
  78   {
  79     PointUZ p;
  80     p.U = operator[]( i ).U * d;
  81     p.Z = operator[]( i ).Z;
  82     res.push_back( p );
  83   }
  84   return res;
  85 }
  86
  87
  88
  89
  90 HYDROData_DTM::HYDROData_DTM()
  91 {
  92 }
  93
  94 HYDROData_DTM::~HYDROData_DTM()
  95 {
  96 }
  97
  98 HYDROData_SequenceOfObjects HYDROData_DTM::GetProfiles() const
  99 {
 100   return GetReferenceObjects( DataTag_Profiles );
 101 }
 102
 103 void HYDROData_DTM::SetProfiles( const HYDROData_SequenceOfObjects& theProfiles )
 104 {
 105   SetReferenceObjects( theProfiles, DataTag_Profiles );
 106 }
 107
 108 double HYDROData_DTM::GetDDZ() const
 109 {
 110   return GetDouble( DataTag_DDZ );
 111 }
 112
 113 void HYDROData_DTM::SetDDZ( double theDDZ )
 114 {
 115   SetDouble( DataTag_DDZ, theDDZ );
 116 }
 117
 118 double HYDROData_DTM::GetSpatialStep() const
 119 {
 120   return GetDouble( DataTag_SpatialStep );
 121 }
 122
 123 void HYDROData_DTM::SetSpatialStep( double theSpatialStep )
 124 {
 125   SetDouble( DataTag_SpatialStep, theSpatialStep );
 126 }
 127
 128 void HYDROData_DTM::Update()
 129 {
 130   HYDROData_SequenceOfObjects objs = GetProfiles();
 131   int aLower = objs.Lower(), anUpper = objs.Upper();
 132   size_t n = anUpper-aLower+1;
 133
 134   std::vector<Handle_HYDROData_Profile> profiles;
 135   profiles.reserve( n );
 136   for( int i=aLower; i<=anUpper; i++ )
 137   {
 138     Handle(HYDROData_Profile) aProfile = Handle(HYDROData_Profile)::DownCast( objs.Value( i ) );
 139     if( !aProfile.IsNull() )
 140       profiles.push_back( aProfile );
 141   }
 142
 143   double ddz = GetDDZ();
 144   double step = GetSpatialStep();
 145   const double EPS = 1E-3;
 146   AltitudePoints points;
 147   AltitudePoints left;
 148   AltitudePoints right;
 149   std::vector<AltitudePoints> main_profiles;
 150
 151   if( ddz>EPS && step>EPS )
 152     points = Interpolate( profiles, ddz, step, left, right, main_profiles );
 153
 154   SetAltitudePoints( points );
 155   //SetTopShape(); //2d
 156   //SetShape3D(); //3d
 157 }
 158
 159
 160
 161
 162
 163
 164
 165
 166
 167
 168 void HYDROData_DTM::GetProperties( const Handle_HYDROData_Profile& theProfile,
 169                     gp_Pnt& theLowestPoint, gp_Vec2d& theDir,
 170                     bool isNormalDir,
 171                     double& theZMin, double& theZMax )
 172 {
 173   theLowestPoint = theProfile->GetBottomPoint();
 174
 175   gp_XY aLeft, aRight;
 176   theProfile->GetLeftPoint( aLeft, true, true );
 177   theProfile->GetRightPoint( aRight, true, true );
 178   double x = aRight.X()-aLeft.X();
 179   double y = aRight.Y()-aLeft.Y();
 180   if( isNormalDir )
 181     theDir = gp_Vec2d( -y, x );
 182   else
 183     theDir = gp_Vec2d( x, y );
 184
 185   HYDROData_Profile::ProfilePoints points = theProfile->GetProfilePoints();
 186   int lo = points.Lower();
 187   int up = points.Upper();
 188   theZMin = std::numeric_limits<double>::max();
 189   theZMax = -theZMin;
 190   for( int i=lo; i<=up; i++ )
 191   {
 192     double z = points.Value( i ).Z();
 193     if( z>theZMax )
 194       theZMax = z;
 195     if( z<theZMin )
 196       theZMin = z;
 197   }
 198 }
 199
 200 inline gp_Pnt2d To2D( const gp_Pnt& thePnt, const gp_Trsf& theTr,
 201                       double& theUMin, double& theUMax )
 202 {
 203   gp_Pnt p = thePnt.Transformed( theTr );
 204   double u = p.X();
 205   double z = p.Z();
 206   if( u<theUMin )
 207     theUMin = u;
 208   if( u>theUMax )
 209     theUMax = u;
 210   return gp_Pnt2d( u, z );
 211 }
 212
 213 Handle(TColgp_HArray1OfPnt2d) To2D( const TColgp_Array1OfPnt& thePoints,
 214                                     const gp_Trsf& theTr,
 215                                     double& theUMin, double& theUMax )
 216 {
 217   int low = thePoints.Lower(), up = thePoints.Upper();
 218   Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( low, up );
 219   for( int i=low; i<=up; i++ )
 220     points->SetValue( i, To2D( thePoints.Value( i ), theTr, theUMin, theUMax ) );
 221   return points;
 222 }
 223
 224 Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve,
 225                                 Standard_Real theFirst, Standard_Real theLast,
 226                                 const gp_Trsf& theTr,
 227                                 double& theUMin, double& theUMax )
 228 {
 229   if( theCurve->IsKind( STANDARD_TYPE( Geom_Line ) ) )
 230   {
 231     gp_Pnt aFirstPnt, aLastPnt;
 232     theCurve->D0( theFirst, aFirstPnt );
 233     theCurve->D0( theLast, aLastPnt );
 234
 235     gp_Pnt2d
 236       aFirst2d = To2D( aFirstPnt, theTr, theUMin, theUMax ),
 237       aLast2d = To2D( aLastPnt, theTr, theUMin, theUMax );
 238
 239     gp_Vec2d dir( aFirst2d, aLast2d );
 240     Handle_Geom2d_Line aLine2d = new Geom2d_Line( aFirst2d, gp_Dir2d( dir.X(), dir.Y() ) );
 241     return new Geom2d_TrimmedCurve( aLine2d, 0, aLast2d.Distance( aFirst2d ) );
 242   }
 243
 244   if( theCurve->IsKind( STANDARD_TYPE( Geom_BSplineCurve ) ) )
 245   {
 246     Handle(Geom_BSplineCurve) aSpline = Handle(Geom_BSplineCurve)::DownCast( theCurve );
 247
 248     Handle(TColgp_HArray1OfPnt2d) poles = To2D( aSpline->Poles(), theTr, theUMin, theUMax );
 249     const TColStd_Array1OfReal& knots = aSpline->Knots();
 250     const TColStd_Array1OfInteger& multiplicities = aSpline->Multiplicities();
 251     int aDegree = aSpline->Degree();
 252
 253     return new Geom2d_BSplineCurve( poles->Array1(), knots, multiplicities, aDegree );
 254   }
 255
 256   return Handle(Geom2d_Curve)();
 257 }
 258
 259 Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis(
 260   const std::vector<Handle_HYDROData_Profile>& theProfiles,
 261   std::vector<double>& theDistances )
 262 {
 263   size_t n = theProfiles.size();
 264   Handle_Geom2d_BSplineCurve aResult;
 265
 266   Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( 1, (int)n );
 267   TColgp_Array1OfVec2d tangents( 1, (int)n );
 268   Handle(TColStd_HArray1OfBoolean) flags = new TColStd_HArray1OfBoolean( 1, (int)n );
 269
 270   for( size_t i = 1; i <= n; i++ )
 271   {
 272     Handle_HYDROData_Profile aProfile = theProfiles[i-1];
 273     aProfile->Update();
 274
 275     gp_Pnt aLowest;
 276     gp_Vec2d aTangent;
 277     double zmin, zmax;
 278     GetProperties( aProfile, aLowest, aTangent, true, zmin, zmax );
 279     aTangent.Normalize();
 280
 281     points->SetValue( (int)i, gp_Pnt2d( aLowest.X(), aLowest.Y() ) );
 282     tangents.SetValue( (int)i, aTangent );
 283     flags->SetValue( (int)i, Standard_True );
 284   }
 285
 286   Geom2dAPI_Interpolate anInterpolator( points, Standard_False, Standard_False );
 287   anInterpolator.Load( tangents, flags );
 288   anInterpolator.Perform();
 289   if( anInterpolator.IsDone() )
 290   {
 291     aResult = anInterpolator.Curve();
 292
 293     //fill the distances vector
 294     Geom2dAdaptor_Curve anAdaptor( aResult );
 295
 296     theDistances.clear();
 297     theDistances.reserve( n );
 298     Standard_Real aParamFirst = anAdaptor.FirstParameter(), aParamLast = anAdaptor.LastParameter();
 299     for( size_t i = 1; i <= n; i++ )
 300     {
 301       gp_Pnt2d aPnt = points->Value( (Standard_Integer)i );
 302       Geom2dAPI_ProjectPointOnCurve aProject( aPnt, aResult );
 303       Standard_Real aParam = aProject.LowerDistanceParameter();
 304       double aDistance = GCPnts_AbscissaPoint::Length( anAdaptor, aParamFirst, aParam );
 305       theDistances.push_back( aDistance );
 306     }
 307   }
 308   return aResult;
 309 }
 310
 311 std::vector<Handle_Geom2d_Curve> HYDROData_DTM::ProfileToParametric(
 312   const Handle_HYDROData_Profile& theProfile,
 313   double& theUMin, double& theUMax, gp_Vec2d& theDir )
 314 {
 315   std::vector<Handle_Geom2d_Curve> curves;
 316
 317   // Transformation of the coordinate systems
 318   gp_Pnt aLowest;
 319   double zmin, zmax;
 320   GetProperties( theProfile, aLowest, theDir, false, zmin, zmax );
 321
 322   gp_Ax3 aStd3d( gp_Pnt( 0, 0, 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( 1, 0, 0 ) );
 323   gp_Ax3 aLocal( aLowest, gp_Dir( 0, 0, 1 ), gp_Dir( theDir.X(), theDir.Y(), 0 ) );
 324
 325   gp_Trsf aTransf;
 326   aTransf.SetTransformation( aStd3d, aLocal );
 327
 328   // Iteration via edges
 329   TopoDS_Wire aWire = TopoDS::Wire( theProfile->GetShape3D() );
 330   TopExp_Explorer anExp( aWire, TopAbs_EDGE );
 331   for( ; anExp.More(); anExp.Next() )
 332   {
 333     // Extract an edge from wire
 334     TopoDS_Edge anEdge = TopoDS::Edge( anExp.Current() );
 335
 336     // Extract a curve corresponding to the edge
 337     TopLoc_Location aLoc;
 338     Standard_Real aFirst, aLast;
 339     Handle(Geom_Curve) aCurve = BRep_Tool::Curve( anEdge, aLoc, aFirst, aLast );
 340
 341     // Convert the curve to 2d CS
 342     Handle(Geom2d_Curve) aCurve2d = CurveTo2D( aCurve, aFirst, aLast, aTransf, theUMin, theUMax );
 343     if( !aCurve2d.IsNull() )
 344       curves.push_back( aCurve2d );
 345   }
 346   return curves;
 347 }
 348
 349
 350 bool CalcMidWidth( const std::vector<gp_Pnt2d>& intersections, double& theMid, double& theWid )
 351 {
 352   double umin = std::numeric_limits<double>::max(),
 353          umax = -umin;
 354
 355   size_t n = intersections.size();
 356   if( n <= 1 )
 357     return false;
 358
 359   for( size_t i = 0; i < n; i++ )
 360   {
 361     double u = intersections[i].X();
 362     if( u<umin )
 363       umin = u;
 364     if( u>umax )
 365       umax = u;
 366   }
 367   theMid = ( umin+umax )/2;
 368   theWid = umax-umin;
 369   return true;
 370 }
 371
 372 void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& theProfile,
 373                                            double theXCurv, double theMinZ, double theMaxZ, double theDDZ,
 374                                            CurveUZ& theMidPointCurve,
 375                                            CurveUZ& theWidthCurve,
 376                                            double theTolerance )
 377 {
 378   double aDblMax = std::numeric_limits<double>::max(),
 379          aUMin = aDblMax,
 380          aUMax = -aUMin,
 381          aVMax = 1000000;
 382
 383   gp_Vec2d aProfileDir;
 384   std::vector<Handle_Geom2d_Curve> curves = ProfileToParametric( theProfile, aUMin, aUMax, aProfileDir );
 385   size_t n = curves.size();
 386
 387   if( n==0 )
 388     return;
 389
 390   // we add the "virtual" vertical lines to simulate the intersection with profile
 391   gp_Pnt2d aFirst, aLast;
 392   curves[0]->D0( curves[0]->FirstParameter(), aFirst );
 393   curves[n-1]->D0( curves[n-1]->LastParameter(), aLast );
 394   Handle(Geom2d_Line) aV1 = new Geom2d_Line( aFirst, gp_Dir2d( 0, 1 ) );
 395   Handle(Geom2d_TrimmedCurve) aT1 = new Geom2d_TrimmedCurve( aV1, 0.0, aVMax );
 396
 397   Handle(Geom2d_Line) aV2 = new Geom2d_Line( aLast, gp_Dir2d( 0, 1 ) );
 398   Handle(Geom2d_TrimmedCurve) aT2 = new Geom2d_TrimmedCurve( aV2, 0.0, aVMax );
 399
 400   curves.push_back( aT1 );
 401   curves.push_back( aT2 );
 402
 403   int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 );
 404   theMidPointCurve = CurveUZ( theXCurv, aProfileDir );
 405   theMidPointCurve.reserve( psize );
 406   theWidthCurve = CurveUZ( theXCurv, aProfileDir );
 407   theWidthCurve.reserve( psize );
 408
 409   n = curves.size();
 410   // for each discrete value of z we search intersection with profile
 411   for( double z = theMinZ; z <= theMaxZ; z += theDDZ )
 412   {
 413     Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z ), gp_Dir2d( 1, 0 ) );
 414     std::vector<gp_Pnt2d> intersections;
 415     for( size_t i = 0; i < n; i++ )
 416     {
 417       Handle_Geom2d_Curve aCurve = curves[i];
 418       Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
 419       for( int k=1, m=anIntersect.NbPoints(); k<=m; k++ )
 420         intersections.push_back( anIntersect.Point( k ) );
 421     }
 422
 423     if( intersections.size() >= 2 )
 424     {
 425       double u_mid, u_wid;
 426       if( !CalcMidWidth( intersections, u_mid, u_wid ) )
 427         continue;
 428
 429       PointUZ p_mid;
 430       p_mid.U = u_mid;
 431       p_mid.Z = z;
 432       theMidPointCurve.push_back( p_mid );
 433
 434       PointUZ p_wid;
 435       p_wid.U = u_wid;
 436       p_wid.Z = z;
 437       theWidthCurve.push_back( p_wid );
 438     }
 439   }
 440 }
 441
 442 void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCurveB,
 443                                  int theNbSteps, std::vector<CurveUZ>& theInterpolation,
 444                                  bool isAddSecond )
 445 {
 446   theInterpolation.clear();
 447   int d = isAddSecond ? 2 : 1;
 448   theInterpolation.reserve( theNbSteps+d );
 449   double dt = 1.0 / double( theNbSteps + 1 );
 450   double t = dt;
 451   theInterpolation.push_back( theCurveA );
 452   for( int i=0; i<theNbSteps; i++, t+=dt )
 453   {
 454     CurveUZ anInterp = theCurveA*(1-t) + theCurveB*t;
 455     theInterpolation.push_back( anInterp );
 456   }
 457   if( isAddSecond )
 458     theInterpolation.push_back( theCurveB );
 459 }
 460 #include <BRepLib_MakeEdge2d.hxx>
 461 void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
 462                                const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
 463                                AltitudePoints& thePoints, double dz )
 464 {
 465   Geom2dAdaptor_Curve anAdaptor( theHydraulicAxis );
 466   TopoDS_Edge E2d = BRepLib_MakeEdge2d(theHydraulicAxis).Edge();
 467   GCPnts_AbscissaPoint ap( anAdaptor, theMidCurve.Xcurv(), anAdaptor.FirstParameter() );
 468   double aParam = ap.Parameter();
 469
 470   gp_Pnt2d point;
 471   anAdaptor.D0( aParam, point );
 472   gp_Vec2d profile_dir = theMidCurve.ProfileDir();
 473   gp_Dir tangent_n( -profile_dir.Y(), profile_dir.X(), dz );
 474   profile_dir.Normalize();
 475
 476   size_t n = theMidCurve.size();
 477   std::map<double, AltitudePoint> sorted_points;
 478   for( size_t i=0; i<n; i++ )
 479   {
 480     double param1 = theMidCurve[i].U - theWidthCurve[i].U / 2;
 481     double param2 = theMidCurve[i].U + theWidthCurve[i].U / 2;
 482
 483     gp_Pnt2d p1 = point.Translated( param1 * profile_dir);
 484     gp_Pnt2d p2 = point.Translated( param2 * profile_dir);
 485
 486     double z = theMidCurve[i].Z;
 487
 488     AltitudePoint p3d_1( p1.X(), p1.Y(), z ), p3d_2( p2.X(), p2.Y(), z );
 489
 490     sorted_points[param1] = p3d_1;
 491     sorted_points[param2] = p3d_2;
 492   }
 493
 494   thePoints.reserve( sorted_points.size() );
 495   std::map<double, AltitudePoint>::const_iterator it = sorted_points.begin(), last = sorted_points.end();
 496   for( ; it!=last; it++ )
 497     thePoints.push_back( it->second );
 498 }
 499
 500 inline double max( double a, double b )
 501 {
 502   if( a>b )
 503     return a;
 504   else
 505     return b;
 506 }
 507
 508 #include <BRepLib_MakeWire.hxx>
 509
 510 std::vector<HYDROData_Bathymetry::AltitudePoints> HYDROData_DTM::Interpolate
 511   ( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
 512     const Handle_HYDROData_Profile& theProfileA,
 513     double theXCurvA,
 514     const Handle_HYDROData_Profile& theProfileB,
 515     double theXCurvB,
 516     double theDDZ, int theNbSteps, bool isAddSecond )
 517 {
 518   double zminA, zmaxA, zminB, zmaxB;
 519   gp_Pnt lowestA, lowestB;
 520   gp_Vec2d dirA, dirB;
 521
 522   GetProperties( theProfileA, lowestA, dirA, false, zminA, zmaxA );
 523   GetProperties( theProfileB, lowestB, dirB, false, zminB, zmaxB );
 524
 525   double dz = zminB - zminA;
 526
 527   double zmin = max( zminA, zminB );
 528   double zmax = max( zmaxA, zmaxB );
 529
 530   CurveUZ midA(0, gp_Vec2d()), midB(0, gp_Vec2d());
 531   CurveUZ widA(0, gp_Vec2d()), widB(0, gp_Vec2d());
 532
 533   ProfileDiscretization( theProfileA, theXCurvA, zmin, zmax, theDDZ, midA, widA );
 534   ProfileDiscretization( theProfileB, theXCurvB, zmin, zmax, theDDZ, midB, widB );
 535
 536   std::vector<CurveUZ> mid, wid;
 537   Interpolate( midA, midB, theNbSteps, mid, isAddSecond );
 538   Interpolate( widA, widB, theNbSteps, wid, isAddSecond );
 539
 540   size_t p = mid.size();
 541   size_t q = p>0 ? 2*mid[0].size() : 1;
 542   std::vector<AltitudePoints> points;
 543   points.resize( p );
 544
 545   for( size_t i=0; i<p; i++ )
 546   {
 547     points[i].reserve( q );
 548     CurveTo3D( theHydraulicAxis, mid[i], wid[i], points[i], dz );
 549   }
 550
 551   return points;
 552 }
 553
 554 HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
 555   ( const std::vector<Handle_HYDROData_Profile>& theProfiles,
 556     double theDDZ, double theSpatialStep,
 557     AltitudePoints& theLeft,
 558     AltitudePoints& theRight,
 559     std::vector<AltitudePoints>& theMainProfiles )
 560 {
 561   AltitudePoints points;
 562   size_t n = theProfiles.size();
 563   if( n<=1 )
 564     return points;
 565
 566   std::vector<double> distances;
 567   Handle_Geom2d_BSplineCurve aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
 568   if( aHydraulicAxis.IsNull() )
 569     return points;
 570
 571   theMainProfiles.reserve( n );
 572
 573   for( size_t i=0, n1=n-1; i<n1; i++ )
 574   {
 575     double aDistance = distances[i+1]-distances[i];
 576     int aNbSteps = int(aDistance/theSpatialStep);
 577     bool isAddSecond = i==n1-1;
 578
 579     // 1. Calculate interpolated profiles
 580     std::vector<AltitudePoints> local_points = Interpolate( aHydraulicAxis, theProfiles[i], distances[i],
 581       theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond );
 582     int lps = local_points.size();
 583
 584     // 2. Put all points into the global container
 585     for( size_t j=0; j<lps; j++ )
 586     {
 587       const AltitudePoints& lp = local_points[j];
 588       if( i==0 && j==0 )
 589         points.reserve( lp.size() * n );
 590       for( size_t k=0, ks=lp.size(); k<ks; k++ )
 591         points.push_back( lp[k] );
 592     }
 593
 594     // 3. Get left/right banks' points
 595     if( i==0 )
 596     {
 597       theLeft.reserve( lps * n );
 598       theRight.reserve( lps * n );
 599     }
 600     for( size_t j=0; j<lps; j++ )
 601     {
 602       const AltitudePoints& lp = local_points[j];
 603       theLeft.push_back( lp[0] );
 604       theRight.push_back( lp[lp.size()-1] );
 605     }
 606
 607     // 4. Get main profiles points
 608     theMainProfiles.push_back( local_points[0] );
 609     if( isAddSecond )
 610       theMainProfiles.push_back( local_points[lps-1] );
 611   }
 612   return points;
 613 }