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