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 bool CalcMidWidth( const std::vector<gp_Pnt2d>& intersections, double& theMid, double& theWid )
 393 {
 394   double umin = std::numeric_limits<double>::max(),
 395          umax = -umin;
 396
 397   size_t n = intersections.size();
 398   if( n <= 1 )
 399     return false;
 400
 401   for( size_t i = 0; i < n; i++ )
 402   {
 403     double u = intersections[i].X();
 404     if( u<umin )
 405       umin = u;
 406     if( u>umax )
 407       umax = u;
 408   }
 409   theMid = ( umin+umax )/2;
 410   theWid = umax-umin;
 411   return true;
 412 }
 413
 414 void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& theProfile,
 415                                            double theXCurv, double theMinZ, double theMaxZ, double theDDZ,
 416                                            CurveUZ& theMidPointCurve,
 417                                            CurveUZ& theWidthCurve,
 418                                            double theTolerance )
 419 {
 420   double aDblMax = std::numeric_limits<double>::max(),
 421          aUMin = aDblMax,
 422          aUMax = -aUMin,
 423          aVMax = 1000000;
 424
 425   std::vector<Handle_Geom2d_Curve> curves = ProfileToParametric( theProfile, aUMin, aUMax );
 426   size_t n = curves.size();
 427
 428   if( n==0 )
 429     return;
 430
 431   // we add the "virtual" vertical lines to simulate the intersection with profile
 432   gp_Pnt2d aFirst, aLast;
 433   curves[0]->D0( curves[0]->FirstParameter(), aFirst );
 434   curves[n-1]->D0( curves[n-1]->LastParameter(), aLast );
 435   Handle(Geom2d_Line) aV1 = new Geom2d_Line( aFirst, gp_Dir2d( 0, 1 ) );
 436   Handle(Geom2d_TrimmedCurve) aT1 = new Geom2d_TrimmedCurve( aV1, 0.0, aVMax );
 437
 438   Handle(Geom2d_Line) aV2 = new Geom2d_Line( aLast, gp_Dir2d( 0, 1 ) );
 439   Handle(Geom2d_TrimmedCurve) aT2 = new Geom2d_TrimmedCurve( aV2, 0.0, aVMax );
 440
 441   curves.push_back( aT1 );
 442   curves.push_back( aT2 );
 443
 444   int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 );
 445   theMidPointCurve = CurveUZ( theXCurv );
 446   theMidPointCurve.reserve( psize );
 447   theWidthCurve = CurveUZ( theXCurv );
 448   theWidthCurve.reserve( psize );
 449
 450   n = curves.size();
 451   // for each discrete value of z we search intersection with profile
 452   for( double z = theMinZ; z <= theMaxZ; z += theDDZ )
 453   {
 454     Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z ), gp_Dir2d( 1, 0 ) );
 455     std::vector<gp_Pnt2d> intersections;
 456     for( size_t i = 0; i < n; i++ )
 457     {
 458       Handle_Geom2d_Curve aCurve = curves[i];
 459       Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
 460       for( int k=1, m=anIntersect.NbPoints(); k<=m; k++ )
 461         intersections.push_back( anIntersect.Point( k ) );
 462     }
 463
 464     if( intersections.size() >= 2 )
 465     {
 466       double u_mid, u_wid;
 467       if( !CalcMidWidth( intersections, u_mid, u_wid ) )
 468         continue;
 469
 470       PointUZ p_mid;
 471       p_mid.U = u_mid;
 472       p_mid.Z = z;
 473       theMidPointCurve.push_back( p_mid );
 474
 475       PointUZ p_wid;
 476       p_wid.U = u_wid;
 477       p_wid.Z = z;
 478       theWidthCurve.push_back( p_wid );
 479     }
 480   }
 481 }
 482
 483 void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCurveB,
 484                                  int theNbSteps, std::vector<CurveUZ>& theInterpolation,
 485                                  bool isAddSecond )
 486 {
 487   theInterpolation.clear();
 488   int d = isAddSecond ? 2 : 1;
 489   theInterpolation.reserve( theNbSteps+d );
 490   double dt = 1.0 / double( theNbSteps + 1 );
 491   double t = dt;
 492   theInterpolation.push_back( theCurveA );
 493   for( int i=0; i<theNbSteps; i++, t+=dt )
 494   {
 495     CurveUZ anInterp = theCurveA*(1-t) + theCurveB*t;
 496     theInterpolation.push_back( anInterp );
 497   }
 498   if( isAddSecond )
 499     theInterpolation.push_back( theCurveB );
 500 }
 501 #include <BRepLib_MakeEdge2d.hxx>
 502 void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
 503                                const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
 504                                AltitudePoints& thePoints,
 505                                Bank* theLeft, Bank* theRight, double dz )
 506 {
 507   Geom2dAdaptor_Curve anAdaptor( theHydraulicAxis );
 508   TopoDS_Edge E2d = BRepLib_MakeEdge2d(theHydraulicAxis).Edge();
 509   GCPnts_AbscissaPoint ap( anAdaptor, theMidCurve.Xcurv(), anAdaptor.FirstParameter() );
 510   double aParam = ap.Parameter();
 511
 512   gp_Pnt2d point;
 513   gp_Vec2d tangent, profile_dir;
 514   anAdaptor.D1( aParam, point, tangent );
 515   profile_dir.SetCoord( tangent.Y(), -tangent.X() );
 516   profile_dir.Normalize();
 517   gp_Dir tangent_n( tangent.X(), tangent.Y(), -dz );
 518   //gp_Dir tangent_n( 0, 0, 1 );
 519   size_t n = theMidCurve.size();
 520   double min_param = 1E+15;
 521   double max_param = -1E+15;
 522   double z1, z2;
 523   for( size_t i=0; i<n; i++ )
 524   {
 525     double param1 = theMidCurve[i].U - theWidthCurve[i].U / 2;
 526     double param2 = theMidCurve[i].U + theWidthCurve[i].U / 2;
 527
 528     gp_Pnt2d p1 = point.Translated( param1 * profile_dir);
 529     gp_Pnt2d p2 = point.Translated( param2 * profile_dir);
 530
 531     double z = theMidCurve[i].Z;
 532
 533     if( param1 < min_param )
 534     {
 535       min_param = param1;
 536       z1 = z;
 537     }
 538     if( param2 < min_param )
 539     {
 540       min_param = param2;
 541       z1 = z;
 542     }
 543     if( param1 > max_param )
 544     {
 545       max_param = param1;
 546       z2 = z;
 547     }
 548     if( param2 > max_param )
 549     {
 550       max_param = param2;
 551       z2 = z;
 552     }
 553
 554     AltitudePoint p3d_1( p1.X(), p1.Y(), z ), p3d_2( p2.X(), p2.Y(), z );
 555     thePoints.push_back( p3d_1 );
 556     thePoints.push_back( p3d_2 );
 557   }
 558
 559   if( theLeft )
 560   {
 561     gp_Pnt2d left2d = point/*.Translated( min_param * profile_dir )*/;
 562     gp_Pnt left( left2d.X(), left2d.Y(), z1 );
 563     theLeft->push_back( left, tangent_n );
 564   }
 565   if( theRight )
 566   {
 567     gp_Pnt2d right2d = point.Translated( max_param * profile_dir );
 568     gp_Pnt right( right2d.X(), right2d.Y(), z2 );
 569     theRight->push_back( right, tangent_n );
 570   }
 571 }
 572
 573 inline double max( double a, double b )
 574 {
 575   if( a>b )
 576     return a;
 577   else
 578     return b;
 579 }
 580
 581 #include <BRepLib_MakeWire.hxx>
 582
 583 HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
 584   ( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
 585     const Handle_HYDROData_Profile& theProfileA,
 586     double theXCurvA,
 587     const Handle_HYDROData_Profile& theProfileB,
 588     double theXCurvB,
 589     double theDDZ, int theNbSteps, bool isAddSecond,
 590     Bank* theLeft, Bank* theRight )
 591 {
 592   double zminA, zmaxA, zminB, zmaxB;
 593   gp_Pnt lowestA, lowestB;
 594   gp_Vec2d dirA, dirB;
 595
 596   GetProperties( theProfileA, lowestA, dirA, false, zminA, zmaxA );
 597   GetProperties( theProfileB, lowestB, dirB, false, zminB, zmaxB );
 598
 599   double dz = zminB - zminA;
 600
 601   double zmin = max( zminA, zminB );
 602   double zmax = max( zmaxA, zmaxB );
 603
 604   CurveUZ midA(0), midB(0);
 605   CurveUZ widA(0), widB(0);
 606
 607   ProfileDiscretization( theProfileA, theXCurvA, zmin, zmax, theDDZ, midA, widA );
 608   ProfileDiscretization( theProfileB, theXCurvB, zmin, zmax, theDDZ, midB, widB );
 609
 610   std::vector<CurveUZ> mid, wid;
 611   Interpolate( midA, midB, theNbSteps, mid, isAddSecond );
 612   Interpolate( widA, widB, theNbSteps, wid, isAddSecond );
 613
 614   size_t p = mid.size();
 615   size_t q = p>0 ? 2*mid[0].size() : 1;
 616   AltitudePoints points;
 617   points.reserve( p*q );
 618   for( size_t i=0; i<p; i++ )
 619     CurveTo3D( theHydraulicAxis, mid[i], wid[i], points, theLeft, theRight, dz );
 620
 621   BRepLib_MakeWire WM;
 622   for (int i =0; i < theLeft->myPoints.size() - 1; i++)
 623     WM.Add(BRepLib_MakeEdge(theLeft->myPoints[i], theLeft->myPoints[i+1]).Edge());
 624   TopoDS_Wire W = WM.Wire();
 625   return points;
 626 }
 627
 628 HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
 629   ( const std::vector<Handle_HYDROData_Profile>& theProfiles,
 630     double theDDZ, double theSpatialStep,
 631     Bank* theLeft,
 632     Bank* theRight )
 633 {
 634   AltitudePoints points;
 635   size_t n = theProfiles.size();
 636   if( n<=1 )
 637     return points;
 638
 639   std::vector<double> distances;
 640   Handle_Geom2d_BSplineCurve aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
 641   if( aHydraulicAxis.IsNull() )
 642     return points;
 643
 644   int aNbStepsComplete = 0;
 645   for( size_t i=0, n1=n-1; i<n1; i++ )
 646   {
 647     double aDistance = distances[i+1]-distances[i];
 648     aNbStepsComplete += ( int(aDistance/theSpatialStep) /*+ 1*/ );
 649   }
 650   if( theLeft )
 651     theLeft->reserve( aNbStepsComplete );
 652   if( theRight )
 653     theRight->reserve( aNbStepsComplete );
 654
 655   for( size_t i=0, n1=n-1; i<n1; i++ )
 656   {
 657     double aDistance = distances[i+1]-distances[i];
 658     int aNbSteps = int(aDistance/theSpatialStep) /*+ 1*/;
 659     bool isAddSecond = i==n1-1;
 660
 661     AltitudePoints local_points = Interpolate( aHydraulicAxis, theProfiles[i], distances[i],
 662       theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond, theLeft, theRight );
 663
 664     if( i==0 )
 665       points.reserve( local_points.size() * ( n-1 ) );
 666
 667     for( size_t j=0, m=local_points.size(); j<m; j++ )
 668       points.push_back( local_points[j] );
 669   }
 670   return points;
 671 }
 672
 673 void HYDROData_DTM::CreateBankShapes( TopoDS_Edge& theLeft, TopoDS_Edge& theRight ) const
 674 {
 675   theLeft = myLeft.createEdge3d();
 676   theRight = myRight.createEdge3d();
 677 }