2 #include <HYDROData_DTM.h>
3 #include <HYDROData_Profile.h>
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>
11 #include <TopoDS_Edge.hxx>
12 #include <TopoDS_Wire.hxx>
13 #include <TopExp_Explorer.hxx>
14 #include <BRep_Tool.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>
34 #include <BRepLib_MakeEdge.hxx>
35 #include <BRepLib_MakeWire.hxx>
36 #include <BRep_Builder.hxx>
37 #include <GeomProjLib.hxx>
38 #include <Geom_TrimmedCurve.hxx>
39 #include <BRepTools_WireExplorer.hxx>
40 #include <TopTools_IndexedMapOfShape.hxx>
41 #include <BRepBuilderAPI_MakeFace.hxx>
43 #include <TopTools_IndexedMapOfOrientedShape.hxx>
44 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
46 #include <BRepLib_MakeEdge.hxx>
47 #include <BRepLib_MakeWire.hxx>
48 #include <BRep_Builder.hxx>
49 #include <ShapeAnalysis_Wire.hxx>
50 #include <BRepAlgo_NormalProjection.hxx>
51 #include <ShapeUpgrade_UnifySameDomain.hxx>
52 #include <BRepBuilderAPI_MakePolygon.hxx>
53 #include <BOPAlgo_Builder.hxx>
54 #include <BRepAdaptor_Curve.hxx>
55 #include <GeomProjLib.hxx>
57 #include <TopTools_IndexedDataMapOfShapeShape.hxx>
58 #include <TopTools_ListIteratorOfListOfShape.hxx>
59 #include <TopTools_SequenceOfShape.hxx>
62 IMPLEMENT_STANDARD_RTTIEXT( HYDROData_DTM, HYDROData_Bathymetry )
64 HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv, const gp_Vec2d& theProfileDir, double theDeltaZ, double theMaxZ )
65 : myXcurv( theXCurv ), myProfileDir( theProfileDir ), myDeltaZ( theDeltaZ ), myMaxZ (theMaxZ)
69 HYDROData_DTM::CurveUZ::~CurveUZ()
73 double HYDROData_DTM::CurveUZ::Xcurv() const
78 gp_Vec2d HYDROData_DTM::CurveUZ::ProfileDir() const
83 double HYDROData_DTM::CurveUZ::DeltaZ() const
88 double HYDROData_DTM::CurveUZ::MaxZ() const
93 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator + ( const CurveUZ& c ) const
95 HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv(), ProfileDir() + c.ProfileDir(), DeltaZ() + c.DeltaZ(), MaxZ() + c.MaxZ() );
96 size_t n = size(), n1 = c.size();
99 std::cout << "Warning: different number of points in curves: " << n << ", " << n1 << std::endl;
102 for( int i=0; i<n; i++ )
105 p.U = operator[]( i ).U + c[i].U;
106 p.Z = operator[]( i ).Z;
112 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator * ( double d ) const
114 HYDROData_DTM::CurveUZ res( Xcurv()*d, ProfileDir()*d, DeltaZ()*d, MaxZ()*d );
117 for( int i=0; i<n; i++ )
120 p.U = operator[]( i ).U * d;
121 p.Z = operator[]( i ).Z;
130 HYDROData_DTM::HYDROData_DTM()
134 HYDROData_DTM::~HYDROData_DTM()
138 HYDROData_SequenceOfObjects HYDROData_DTM::GetProfiles() const
140 return GetReferenceObjects( DataTag_Profiles );
143 void HYDROData_DTM::SetProfiles( const HYDROData_SequenceOfObjects& theProfiles )
145 SetReferenceObjects( theProfiles, DataTag_Profiles );
149 double HYDROData_DTM::GetDDZ() const
151 return GetDouble( DataTag_DDZ );
154 void HYDROData_DTM::SetDDZ( double theDDZ )
156 SetDouble( DataTag_DDZ, theDDZ );
160 double HYDROData_DTM::GetSpatialStep() const
162 return GetDouble( DataTag_SpatialStep );
165 void HYDROData_DTM::SetSpatialStep( double theSpatialStep )
167 SetDouble( DataTag_SpatialStep, theSpatialStep );
171 void HYDROData_DTM::PointsToWire(const AltitudePoints& pnts, TopoDS_Wire& W )
174 BRepBuilderAPI_MakePolygon PM;
175 for (int i = 0; i < pnts.size(); i++)
176 PM.Add(gp_Pnt(pnts[i].X, pnts[i].Y, pnts[i].Z));
181 void HYDROData_DTM::PointsToEdge(const AltitudePoints& pnts, TopoDS_Edge& E )
183 Handle(TColgp_HArray1OfPnt) gpPoints = new TColgp_HArray1OfPnt(1, (int)pnts.size());
185 for (int i = 0; i < pnts.size(); i++)
186 gpPoints->SetValue(i+1, gp_Pnt(pnts[i].X, pnts[i].Y, pnts[i].Z));
188 GeomAPI_Interpolate anInterpolator(gpPoints, Standard_False,1.0e-6);
189 anInterpolator.Perform() ;
190 if (anInterpolator.IsDone())
192 Handle(Geom_BSplineCurve) C = anInterpolator.Curve();
193 E = BRepBuilderAPI_MakeEdge(C).Edge();
197 TopTools_IndexedMapOfOrientedShape HYDROData_DTM::Create3DShape(const AltitudePoints& left,
198 const AltitudePoints& right,
199 const std::vector<AltitudePoints>& main_profiles)
201 TopTools_IndexedMapOfOrientedShape ll;
202 //TopoDS_Wire LWire, RWire;
203 //PointsToWire(left, LWire);
204 //PointsToWire(right, RWire);
205 TopoDS_Edge LEdge, REdge;
206 PointsToEdge(left, LEdge);
207 PointsToEdge(right, REdge);
209 ll.Add(LEdge.Oriented(TopAbs_FORWARD));
211 for (int k = 0; k < main_profiles.size(); k++)
214 PointsToWire(main_profiles[k], W);
215 TopAbs_Orientation Ori = TopAbs_INTERNAL;
216 if (k == 0 || k == main_profiles.size() - 1)
217 Ori = TopAbs_FORWARD;
218 ll.Add(W.Oriented(Ori));
222 ll.Add(REdge.Oriented(TopAbs_FORWARD));
223 //yes, add subshapes in this order (left + profiles + right)
224 //otherwise the projected wire will be non-manifold
230 void HYDROData_DTM::Update()
232 AltitudePoints points;
233 TopoDS_Shape Out3dPres;
234 TopoDS_Shape Out2dPres;
235 TopoDS_Shape OutLeftB;
236 TopoDS_Shape OutRightB;
237 TopoDS_Shape OutInlet;
238 TopoDS_Shape OutOutlet;
240 HYDROData_SequenceOfObjects objs = GetProfiles();
241 double ddz = GetDDZ();
242 double step = GetSpatialStep();
243 std::set<int> InvInd;
244 bool WireIntersections; //__TODO
245 CreateProfilesFromDTM( objs, ddz, step, points, Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, true, true, InvInd, -1, WireIntersections );
246 SetAltitudePoints( points );
248 SetShape( DataTag_LeftBankShape, OutLeftB);
249 SetShape( DataTag_RightBankShape, OutRightB);
250 SetShape( DataTag_InletShape, OutInlet);
251 SetShape( DataTag_OutletShape, OutOutlet );
252 SetShape( DataTag_3DShape, Out3dPres );
253 SetShape( DataTag_2DShape, Out2dPres );
255 HYDROData_Bathymetry::Update();
258 void HYDROData_DTM::GetPresentationShapes( TopoDS_Shape& Out3dPres,
259 TopoDS_Shape& Out2dPres,
260 TopoDS_Shape& OutLeftB,
261 TopoDS_Shape& OutRightB,
262 TopoDS_Shape& OutInlet,
263 TopoDS_Shape& OutOutlet )
266 OutLeftB = GetShape( DataTag_LeftBankShape);
267 OutRightB = GetShape( DataTag_RightBankShape);
268 OutInlet = GetShape( DataTag_InletShape);
269 OutOutlet = GetShape( DataTag_OutletShape );
270 Out3dPres = GetShape( DataTag_3DShape );
271 Out2dPres = GetShape( DataTag_2DShape );
273 void HYDROData_DTM::CreateProfilesFromDTM (const HYDROData_SequenceOfObjects& InpProfiles,
276 AltitudePoints& points,
277 TopoDS_Shape& Out3dPres,
278 TopoDS_Shape& Out2dPres,
279 TopoDS_Shape& OutLeftB,
280 TopoDS_Shape& OutRightB,
281 TopoDS_Shape& OutInlet,
282 TopoDS_Shape& OutOutlet,
285 std::set<int>& InvInd,
287 bool& WireIntersections)
289 int aLower = InpProfiles.Lower(), anUpper = InpProfiles.Upper();
290 size_t n = anUpper - aLower + 1;
292 std::vector<Handle(HYDROData_Profile)> profiles;
293 profiles.reserve( n );
294 for( int i=aLower; i<=anUpper; i++ )
296 Handle(HYDROData_Profile) aProfile = Handle(HYDROData_Profile)::DownCast( InpProfiles.Value( i ) );
297 if( !aProfile.IsNull() )
298 profiles.push_back( aProfile );
300 const double EPS = 1E-3;
302 AltitudePoints right;
303 std::vector<AltitudePoints> main_profiles;
305 if( thePntsLimit > 0 )
307 int aNbPoints = EstimateNbPoints( profiles, ddz, step );
308 if( aNbPoints < 0 || aNbPoints > thePntsLimit )
312 if( ddz>EPS && step>EPS )
313 CreateProfiles(profiles, ddz, step, left, right, points, main_profiles,
314 Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, Create3dPres, Create2dPres, InvInd, WireIntersections );
317 bool HYDROData_DTM::GetPlanarFaceFromBanks( const TopoDS_Edge& LB, const TopoDS_Edge& RB, TopoDS_Face& outF,
318 TopTools_SequenceOfShape* Boundr)
322 Handle_Geom_Plane refpl = new Geom_Plane(gp_Pnt(0,0,0), gp_Dir(0,0,1));
324 TopoDS_Vertex VFI, VLI, VFO, VLO;
329 BRepAdaptor_Curve LBAD(LB);
330 Handle_Geom_Curve LBPC = GeomProjLib::ProjectOnPlane(LBAD.Curve().Curve(), refpl, gp_Dir(0, 0, -1), 1 );
331 prLB = BRepLib_MakeEdge(LBPC).Edge();
333 BRepAdaptor_Curve RBAD(RB);
334 Handle_Geom_Curve RBPC = GeomProjLib::ProjectOnPlane(RBAD.Curve().Curve(), refpl, gp_Dir(0, 0, -1), 1 );
335 prRB = BRepLib_MakeEdge(RBPC).Edge();
337 TopExp::Vertices(prLB, VFI, VFO, 1);
338 TopExp::Vertices(prRB, VLI, VLO, 1);
339 TopoDS_Edge prIL = BRepLib_MakeEdge(VFI, VLI).Edge();
340 TopoDS_Edge prOL = BRepLib_MakeEdge(VFO, VLO).Edge();
341 TopoDS_Wire prW = BRepLib_MakeWire(prLB, prIL, prOL, prRB).Wire();
342 outF = BRepBuilderAPI_MakeFace(refpl->Pln(), prW, 1).Face();
346 Boundr->Append(prLB);
347 Boundr->Append(prIL);
348 Boundr->Append(prOL);
349 Boundr->Append(prRB);
352 ShapeAnalysis_Wire WA(prW, outF, Precision::Confusion());
353 bool res = WA.CheckSelfIntersection();
357 void HYDROData_DTM::CreateProfiles(const std::vector<Handle(HYDROData_Profile)>& theProfiles,
359 double theSpatialStep,
360 AltitudePoints& theOutLeft,
361 AltitudePoints& theOutRight,
362 AltitudePoints& theOutPoints,
363 std::vector<AltitudePoints>& theOutMainProfiles,
364 TopoDS_Shape& Out3dPres,
365 TopoDS_Shape& Out2dPres,
366 TopoDS_Shape& OutLeftB,
367 TopoDS_Shape& OutRightB,
368 TopoDS_Shape& OutInlet,
369 TopoDS_Shape& OutOutlet,
372 std::set<int>& InvInd,
375 if (theProfiles.empty())
377 theOutPoints = Interpolate( theProfiles, theDDZ, theSpatialStep, theOutLeft, theOutRight, theOutMainProfiles, InvInd );
378 //note that if Create3dPres is false => Create2dPres flag is meaningless!
379 if( theOutPoints.empty() )
384 TopTools_IndexedMapOfOrientedShape ll = Create3DShape( theOutLeft, theOutRight, theOutMainProfiles);
390 BB.MakeCompound(cmp);
391 for (int i = 1; i <= ll.Extent(); i++)
396 //same order as in HYDROData_DTM::Update()
398 OutRightB = ll(ll.Extent());
400 OutOutlet = ll(ll.Extent() - 1);
405 ProjStat = GetPlanarFaceFromBanks(TopoDS::Edge(OutLeftB), TopoDS::Edge(OutRightB), outF, NULL);
414 void HYDROData_DTM::GetProperties( const Handle(HYDROData_Profile)& theProfile,
415 gp_Pnt& theLowestPoint, gp_Vec2d& theDir,
416 double& theZMin, double& theZMax )
418 theLowestPoint = theProfile->GetBottomPoint();
421 theProfile->GetLeftPoint( aLeft, true, true );
422 theProfile->GetRightPoint( aRight, true, true );
423 double x = aRight.X()-aLeft.X();
424 double y = aRight.Y()-aLeft.Y();
425 theDir = gp_Vec2d( x, y );
427 HYDROData_Profile::ProfilePoints points = theProfile->GetProfilePoints();
428 int lo = points.Lower();
429 int up = points.Upper();
430 theZMin = std::numeric_limits<double>::max();
432 for( int i=lo; i<=up; i++ )
434 double z = points.Value( i ).Z();
442 inline gp_Pnt2d To2D( const gp_Pnt& thePnt, const gp_Trsf& theTr,
443 double& theUMin, double& theUMax )
445 gp_Pnt p = thePnt.Transformed( theTr );
452 return gp_Pnt2d( u, z );
455 Handle(TColgp_HArray1OfPnt2d) To2D( const TColgp_Array1OfPnt& thePoints,
456 const gp_Trsf& theTr,
457 double& theUMin, double& theUMax )
459 int low = thePoints.Lower(), up = thePoints.Upper();
460 Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( low, up );
461 for( int i=low; i<=up; i++ )
462 points->SetValue( i, To2D( thePoints.Value( i ), theTr, theUMin, theUMax ) );
466 Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve,
467 Standard_Real theFirst, Standard_Real theLast,
468 const gp_Trsf& theTr,
469 double& theUMin, double& theUMax )
471 if( theCurve->IsKind( STANDARD_TYPE( Geom_Line ) ) )
473 gp_Pnt aFirstPnt, aLastPnt;
474 theCurve->D0( theFirst, aFirstPnt );
475 theCurve->D0( theLast, aLastPnt );
478 aFirst2d = To2D( aFirstPnt, theTr, theUMin, theUMax ),
479 aLast2d = To2D( aLastPnt, theTr, theUMin, theUMax );
481 gp_Vec2d dir( aFirst2d, aLast2d );
482 Handle(Geom2d_Line) aLine2d = new Geom2d_Line( aFirst2d, gp_Dir2d( dir.X(), dir.Y() ) );
483 return new Geom2d_TrimmedCurve( aLine2d, 0, aLast2d.Distance( aFirst2d ) );
486 if( theCurve->IsKind( STANDARD_TYPE( Geom_BSplineCurve ) ) )
488 Handle(Geom_BSplineCurve) aSpline = Handle(Geom_BSplineCurve)::DownCast( theCurve );
490 Handle(TColgp_HArray1OfPnt2d) poles = To2D( aSpline->Poles(), theTr, theUMin, theUMax );
491 const TColStd_Array1OfReal& knots = aSpline->Knots();
492 const TColStd_Array1OfInteger& multiplicities = aSpline->Multiplicities();
493 int aDegree = aSpline->Degree();
495 return new Geom2d_BSplineCurve( poles->Array1(), knots, multiplicities, aDegree );
498 return Handle(Geom2d_Curve)();
501 #include <GCE2d_MakeSegment.hxx>
502 #include <Geom2dAPI_InterCurveCurve.hxx>
503 bool IsCooriented( const Handle(HYDROData_Profile)& theProfile1,
504 const Handle(HYDROData_Profile)& theProfile2 )
506 if( theProfile1==theProfile2 )
509 gp_XY lp1, rp1, lp2, rp2;
510 theProfile1->GetLeftPoint(lp1);
511 theProfile1->GetRightPoint(rp1);
512 theProfile2->GetLeftPoint(lp2);
513 theProfile2->GetRightPoint(rp2);
515 GCE2d_MakeSegment s1(lp1, lp2);
516 GCE2d_MakeSegment s2(rp1, rp2);
518 Geom2dAPI_InterCurveCurve inter;
519 inter.Init(s1.Value(), s2.Value());
520 if (inter.NbPoints() == 0)
526 Handle(Geom2d_BSplineCurve) HYDROData_DTM::CreateHydraulicAxis(
527 const std::vector<Handle(HYDROData_Profile)>& theProfiles,
528 std::vector<double>& theDistances )
530 size_t n = theProfiles.size();
532 return Handle_Geom2d_BSplineCurve();
534 Handle_Geom2d_BSplineCurve aResult;
536 Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( 1, (int)n );
537 TColgp_Array1OfVec2d tangents( 1, (int)n );
538 Handle(TColStd_HArray1OfBoolean) flags = new TColStd_HArray1OfBoolean( 1, (int)n );
540 // Stage 1. Orient all profiles to be co-oriented with the first profile
541 theProfiles[0]->Update();
542 for( size_t i = 1; i < n; i++ )
544 Handle(HYDROData_Profile) aProfile = theProfiles[i];
545 Handle(HYDROData_Profile) aPrevProfile = theProfiles[i-1];
547 if( !IsCooriented( aProfile, aPrevProfile ) )
550 aProfile->GetLeftPoint( lp, true );
551 aProfile->GetRightPoint( rp, true );
552 aProfile->SetLeftPoint( rp, true );
553 aProfile->SetRightPoint( lp, true );
558 // Stage 2. Calculate normals so that each normal "points" to the next profile
559 for( size_t i = 0; i < n; i++ )
561 Handle(HYDROData_Profile) aProfile = theProfiles[i];
562 Handle(HYDROData_Profile) aNextProfile = i==n-1 ? theProfiles[i-1] : theProfiles[i+1];
568 gp_XYZ curP = aProfile->GetBottomPoint(true);
569 gp_XY curP2d = gp_XY(curP.X(), curP.Y());
571 gp_XYZ nextP = aNextProfile->GetBottomPoint(true);
572 gp_XY nextP2d = gp_XY(nextP.X(), nextP.Y());
575 GetProperties( aProfile, aLowest, aPrTangent, zmin, zmax );
576 aNormal.SetCoord( -aPrTangent.Y(), aPrTangent.X() );
578 gp_Vec2d aDirToNextProfile(nextP2d.X() - curP2d.X(), nextP2d.Y() - curP2d.Y() );
580 aDirToNextProfile.Reverse();
582 if (aNormal.Dot(aDirToNextProfile) < 0)
587 points->SetValue( (int)(i+1), gp_Pnt2d( aLowest.X(), aLowest.Y() ) );
588 tangents.SetValue( (int)(i+1), aNormal );
589 flags->SetValue( (int)(i+1), Standard_True );
592 Geom2dAPI_Interpolate anInterpolator( points, Standard_False, Standard_False );
593 anInterpolator.Load( tangents, flags );
594 anInterpolator.Perform();
595 if( anInterpolator.IsDone() )
597 aResult = anInterpolator.Curve();
599 //fill the distances vector
600 Geom2dAdaptor_Curve anAdaptor( aResult );
602 theDistances.clear();
603 theDistances.reserve( n );
604 Standard_Real aParamFirst = anAdaptor.FirstParameter(), aParamLast = anAdaptor.LastParameter();
605 for( size_t i = 1; i <= n; i++ )
607 gp_Pnt2d aPnt = points->Value( (Standard_Integer)i );
608 Geom2dAPI_ProjectPointOnCurve aProject( aPnt, aResult );
609 Standard_Real aParam = aProject.LowerDistanceParameter();
610 double aDistance = GCPnts_AbscissaPoint::Length( anAdaptor, aParamFirst, aParam );
611 theDistances.push_back( aDistance );
617 std::vector<Handle(Geom2d_Curve)> HYDROData_DTM::ProfileToParametric(
618 const Handle(HYDROData_Profile)& theProfile,
619 double& theUMin, double& theUMax, gp_Vec2d& theDir )
621 std::vector<Handle(Geom2d_Curve)> curves;
623 // Transformation of the coordinate systems
626 GetProperties( theProfile, aLowest, theDir, zmin, zmax );
628 gp_Ax3 aStd3d( gp_Pnt( 0, 0, 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( 1, 0, 0 ) );
629 gp_Ax3 aLocal( gp_Pnt( aLowest.X(), aLowest.Y(), 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( theDir.X(), theDir.Y(), 0 ) );
632 aTransf.SetTransformation( aStd3d, aLocal );
634 // Iteration via edges
635 TopoDS_Wire aWire = TopoDS::Wire( theProfile->GetShape3D() );
636 TopExp_Explorer anExp( aWire, TopAbs_EDGE );
637 for( ; anExp.More(); anExp.Next() )
639 // Extract an edge from wire
640 TopoDS_Edge anEdge = TopoDS::Edge( anExp.Current() );
642 // Extract a curve corresponding to the edge
643 TopLoc_Location aLoc;
644 Standard_Real aFirst, aLast;
645 Handle(Geom_Curve) aCurve = BRep_Tool::Curve( anEdge, aLoc, aFirst, aLast );
647 // Convert the curve to 2d CS
648 Handle(Geom2d_Curve) aCurve2d = CurveTo2D( aCurve, aFirst, aLast, aTransf, theUMin, theUMax );
649 if( !aCurve2d.IsNull() )
650 curves.push_back( aCurve2d );
656 bool CalcMidWidth( const std::set<double>& intersections, double& theMid, double& theWid )
658 double umin = std::numeric_limits<double>::max(),
661 size_t n = intersections.size();
665 std::set<double>::const_iterator it = intersections.begin(), last = intersections.end();
666 for( ; it!=last; it++ )
674 theMid = ( umin+umax )/2;
679 void HYDROData_DTM::ProfileDiscretization( const Handle(HYDROData_Profile)& theProfile,
680 double theXCurv, double theMinZ, double theMaxZ, double theTopZ, double theDDZ,
681 CurveUZ& theMidPointCurve,
682 CurveUZ& theWidthCurve,
683 int& intersection_nb,
686 double aDblMax = std::numeric_limits<double>::max(),
691 gp_Vec2d aProfileDir;
692 std::vector<Handle(Geom2d_Curve)> curves = ProfileToParametric( theProfile, aUMin, aUMax, aProfileDir );
693 size_t n = curves.size();
698 // we add the "virtual" vertical lines to simulate the intersection with profile
699 gp_Pnt2d aFirst, aLast;
700 curves[0]->D0( curves[0]->FirstParameter(), aFirst );
701 curves[n-1]->D0( curves[n-1]->LastParameter(), aLast );
702 Handle(Geom2d_Line) aV1 = new Geom2d_Line( aFirst, gp_Dir2d( 0, 1 ) );
703 Handle(Geom2d_TrimmedCurve) aT1 = new Geom2d_TrimmedCurve( aV1, 0.0, aVMax );
705 Handle(Geom2d_Line) aV2 = new Geom2d_Line( aLast, gp_Dir2d( 0, 1 ) );
706 Handle(Geom2d_TrimmedCurve) aT2 = new Geom2d_TrimmedCurve( aV2, 0.0, aVMax );
708 curves.push_back( aT1 );
709 curves.push_back( aT2 );
711 int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 );
712 theMidPointCurve = CurveUZ( theXCurv, aProfileDir, theMinZ, theTopZ);
713 theMidPointCurve.reserve( psize );
714 theWidthCurve = CurveUZ( theXCurv, aProfileDir, theMinZ, theTopZ );
715 theWidthCurve.reserve( psize );
718 // for each discrete value of z we search intersection with profile
719 for( double z1 = theMinZ; z1 <= theMaxZ; z1 += theDDZ )
721 Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z1 ), gp_Dir2d( 1, 0 ) );
722 std::set<double> intersections;
723 for( size_t i = 0; i < n; i++ )
725 Handle(Geom2d_Curve) aCurve = curves[i];
726 Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
727 for( int k=1, m=anIntersect.NbPoints(); k<=m; k++ )
728 intersections.insert( anIntersect.Point( k ).X() );
731 intersection_nb = intersections.size();
732 if( intersection_nb >= 1 )
735 if( !CalcMidWidth( intersections, u_mid, u_wid ) )
738 double z = z1 - theMinZ;
742 theMidPointCurve.push_back( p_mid );
747 theWidthCurve.push_back( p_wid );
752 void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCurveB,
753 int theNbSteps, std::vector<CurveUZ>& theInterpolation,
756 theInterpolation.clear();
757 int d = isAddSecond ? 2 : 1;
758 theInterpolation.reserve( theNbSteps+d );
759 double dt = 1.0 / double( theNbSteps + 1 );
761 theInterpolation.push_back( theCurveA );
762 for( int i=0; i<theNbSteps; i++, t+=dt )
764 CurveUZ anInterp = theCurveA*(1-t) + theCurveB*t;
765 theInterpolation.push_back( anInterp );
768 theInterpolation.push_back( theCurveB );
770 #include <BRepLib_MakeEdge2d.hxx>
771 void HYDROData_DTM::CurveTo3D( const Handle(Geom2d_BSplineCurve)& theHydraulicAxis,
772 const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
773 AltitudePoints& thePoints )
775 Geom2dAdaptor_Curve anAdaptor( theHydraulicAxis );
776 TopoDS_Edge E2d = BRepLib_MakeEdge2d(theHydraulicAxis).Edge();
777 GCPnts_AbscissaPoint ap( anAdaptor, theMidCurve.Xcurv(), anAdaptor.FirstParameter() );
778 double aParam = ap.Parameter();
781 anAdaptor.D0( aParam, point );
782 gp_Vec2d profile_dir = theMidCurve.ProfileDir();
783 //gp_Dir tangent_n( -profile_dir.Y(), profile_dir.X(), dz );
784 profile_dir.Normalize();
786 size_t n = theMidCurve.size();
787 std::map<double, AltitudePoint> sorted_points;
788 bool isOnTop = false;
789 for( size_t i=0; i<n; i++ ) // build the two banks of the interpolated profile, from bottom to top
791 double param1 = theMidCurve[i].U - theWidthCurve[i].U / 2;
792 double param2 = theMidCurve[i].U + theWidthCurve[i].U / 2;
794 gp_Pnt2d p1 = point.Translated( param1 * profile_dir);
795 gp_Pnt2d p2 = point.Translated( param2 * profile_dir);
797 bool arrivedOnTop = false;
799 if (theMidCurve[i].Z <= theMidCurve.MaxZ())
800 z = theMidCurve[i].Z + theMidCurve.DeltaZ();
803 z = theMidCurve.MaxZ() + theMidCurve.DeltaZ(); // limit z to linear interpolation between maxima on extremity profiles
804 arrivedOnTop = true; // do not keep points after this one
808 AltitudePoint p3d_1( p1.X(), p1.Y(), z ), p3d_2( p2.X(), p2.Y(), z );
810 sorted_points[param1] = p3d_1;
811 sorted_points[param2] = p3d_2;
814 // isOnTop =true; // do not keep points after this one (commented: leads to strange limits of 2D shape)
817 thePoints.reserve( sorted_points.size() );
818 const double EPS = 1E-12;
819 std::map<double, AltitudePoint>::const_iterator it = sorted_points.begin(), last = sorted_points.end();
820 for( ; it!=last; it++ )
821 if( thePoints.empty() || thePoints.back().SquareDistance( it->second ) > EPS )
822 thePoints.push_back( it->second );
825 inline double max( double a, double b )
833 inline double min( double a, double b )
841 #include <BRepLib_MakeWire.hxx>
843 std::vector<HYDROData_Bathymetry::AltitudePoints> HYDROData_DTM::Interpolate
844 ( const Handle(Geom2d_BSplineCurve)& theHydraulicAxis,
845 const Handle(HYDROData_Profile)& theProfileA,
847 const Handle(HYDROData_Profile)& theProfileB,
849 double theDDZ, int theNbSteps, bool isAddSecond,
850 int& inter_nb_1, int& inter_nb_2)
852 double zminA, zmaxA, zminB, zmaxB;
853 gp_Pnt lowestA, lowestB;
856 GetProperties( theProfileA, lowestA, dirA, zminA, zmaxA );
857 GetProperties( theProfileB, lowestB, dirB, zminB, zmaxB );
860 double hmax = max( zmaxA-zminA, zmaxB-zminB );
862 //double dz = zminB - zminA;
863 //double zmin = min( zminA, zminB );
864 //double zmax = max( zmaxA, zmaxB );
866 CurveUZ midA(0, gp_Vec2d(), 0, 0), midB(0, gp_Vec2d(), 0, 0);
867 CurveUZ widA(0, gp_Vec2d(), 0, 0), widB(0, gp_Vec2d(), 0, 0);
869 ProfileDiscretization( theProfileA, theXCurvA, zminA, zminA+hmax, zmaxA-zminA, theDDZ, midA, widA, inter_nb_1 );
870 ProfileDiscretization( theProfileB, theXCurvB, zminB, zminB+hmax, zmaxB-zminB, theDDZ, midB, widB, inter_nb_2 );
872 std::vector<CurveUZ> mid, wid;
873 Interpolate( midA, midB, theNbSteps, mid, isAddSecond );
874 Interpolate( widA, widB, theNbSteps, wid, isAddSecond );
876 size_t p = mid.size();
877 size_t q = p>0 ? 2*mid[0].size() : 1;
878 std::vector<AltitudePoints> points;
881 for( size_t i=0; i<p; i++ )
883 points[i].reserve( q );
884 CurveTo3D( theHydraulicAxis, mid[i], wid[i], points[i] );
890 HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
891 ( const std::vector<Handle(HYDROData_Profile)>& theProfiles,
892 double theDDZ, double theSpatialStep,
893 AltitudePoints& theLeft,
894 AltitudePoints& theRight,
895 std::vector<AltitudePoints>& theMainProfiles,
896 std::set<int>& invalInd)
898 AltitudePoints points;
899 size_t n = theProfiles.size();
903 std::vector<double> distances;
904 Handle(Geom2d_BSplineCurve) aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
905 if( aHydraulicAxis.IsNull() )
908 theMainProfiles.reserve( n );
910 for( size_t i=0, n1=n-1; i<n1; i++ )
912 double aDistance = distances[i+1]-distances[i];
913 int aNbSteps = int(aDistance/theSpatialStep);
914 bool isAddSecond = i==n1-1;
916 // 1. Calculate interpolated profiles
917 int inter_nb_1, inter_nb_2;
918 std::vector<AltitudePoints> local_points = Interpolate( aHydraulicAxis, theProfiles[i], distances[i],
919 theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond, inter_nb_1, inter_nb_2 );
920 int lps = local_points.size();
926 invalInd.insert(i+1);
928 // 2. Put all points into the global container
929 for( size_t j=0; j<lps; j++ )
931 const AltitudePoints& lp = local_points[j];
933 points.reserve( lp.size() * n );
934 for( size_t k=0, ks=lp.size(); k<ks; k++ )
935 points.push_back( lp[k] );
938 // 3. Get left/right banks' points
941 theLeft.reserve( lps * n );
942 theRight.reserve( lps * n );
944 for( size_t j=0; j<lps; j++ )
946 const AltitudePoints& lp = local_points[j];
947 theLeft.push_back( lp[0] );
948 theRight.push_back( lp[lp.size()-1] );
951 // 4. Get main profiles points
952 theMainProfiles.push_back( local_points[0] );
954 theMainProfiles.push_back( local_points[lps-1] );
959 int HYDROData_DTM::EstimateNbPoints( const std::vector<Handle(HYDROData_Profile)>& theProfiles,
960 double theDDZ, double theSpatialStep )
962 size_t n = theProfiles.size();
965 if( theDDZ<1E-6 || theSpatialStep<1E-6 )
968 std::vector<double> distances;
969 Handle(Geom2d_BSplineCurve) aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
970 if( aHydraulicAxis.IsNull() )
973 double aCompleteDistance = distances[n-1];
974 int aNbSteps = int( aCompleteDistance / theSpatialStep ) + 1;
978 GetProperties( theProfiles[0], aLowest, aDir, aZMin, aZMax );
979 int aNbZSteps = (aZMax-aZMin)/theDDZ;
981 if( aNbSteps > ( 1<<16 ) || aNbZSteps > ( 1<<16 ) )
984 return aNbSteps * aNbZSteps;