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