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