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