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>
44 #include <BRepLib_MakeEdge.hxx>
45 #include <BRepLib_MakeWire.hxx>
46 #include <BRep_Builder.hxx>
47 #include <ShapeAnalysis_Wire.hxx>
48 #include <BRepAlgo_NormalProjection.hxx>
49 #include <ShapeUpgrade_UnifySameDomain.hxx>
50 #include <BRepBuilderAPI_MakePolygon.hxx>
51 #include <BOPAlgo_Builder.hxx>
53 #include <TopTools_IndexedDataMapOfShapeShape.hxx>
54 #include <TopTools_ListIteratorOfListOfShape.hxx>
56 IMPLEMENT_STANDARD_RTTIEXT( HYDROData_DTM, HYDROData_Bathymetry )
58 HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv, const gp_Vec2d& theProfileDir, double theDeltaZ )
59 : myXcurv( theXCurv ), myProfileDir( theProfileDir ), myDeltaZ( theDeltaZ )
63 HYDROData_DTM::CurveUZ::~CurveUZ()
67 double HYDROData_DTM::CurveUZ::Xcurv() const
72 gp_Vec2d HYDROData_DTM::CurveUZ::ProfileDir() const
77 double HYDROData_DTM::CurveUZ::DeltaZ() const
82 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator + ( const CurveUZ& c ) const
84 HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv(), ProfileDir() + c.ProfileDir(), DeltaZ() + c.DeltaZ() );
85 size_t n = size(), n1 = c.size();
88 std::cout << "Warning: different number of points in curves: " << n << ", " << n1 << std::endl;
91 for( int i=0; i<n; i++ )
94 p.U = operator[]( i ).U + c[i].U;
95 p.Z = operator[]( i ).Z;
101 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator * ( double d ) const
103 HYDROData_DTM::CurveUZ res( Xcurv()*d, ProfileDir()*d, DeltaZ()*d );
106 for( int i=0; i<n; i++ )
109 p.U = operator[]( i ).U * d;
110 p.Z = operator[]( i ).Z;
119 HYDROData_DTM::HYDROData_DTM()
123 HYDROData_DTM::~HYDROData_DTM()
127 HYDROData_SequenceOfObjects HYDROData_DTM::GetProfiles() const
129 return GetReferenceObjects( DataTag_Profiles );
132 void HYDROData_DTM::SetProfiles( const HYDROData_SequenceOfObjects& theProfiles )
134 SetReferenceObjects( theProfiles, DataTag_Profiles );
138 double HYDROData_DTM::GetDDZ() const
140 return GetDouble( DataTag_DDZ );
143 void HYDROData_DTM::SetDDZ( double theDDZ )
145 SetDouble( DataTag_DDZ, theDDZ );
149 double HYDROData_DTM::GetSpatialStep() const
151 return GetDouble( DataTag_SpatialStep );
154 void HYDROData_DTM::SetSpatialStep( double theSpatialStep )
156 SetDouble( DataTag_SpatialStep, theSpatialStep );
160 void HYDROData_DTM::PointToWire(const AltitudePoints& pnts, TopoDS_Wire& W )
162 /*BRepLib_MakeWire WM;
165 for (int i = 0; i < pnts.size() - 1; i++)
167 gp_Pnt p1(pnts[i].X, pnts[i].Y, pnts[i].Z);
168 gp_Pnt p2(pnts[i+1].X, pnts[i+1].Y, pnts[i+1].Z);
169 WM.Add(BRepLib_MakeEdge(p1, p2).Edge());
173 BRepBuilderAPI_MakePolygon PM;
174 for (int i = 0; i < pnts.size(); i++)
175 PM.Add(gp_Pnt(pnts[i].X, pnts[i].Y, pnts[i].Z));
180 TopTools_IndexedMapOfOrientedShape HYDROData_DTM::Create3DShape(const AltitudePoints& left,
181 const AltitudePoints& right,
182 const std::vector<AltitudePoints>& main_profiles)
184 TopTools_IndexedMapOfOrientedShape ll;
185 TopoDS_Wire LWire, RWire;
186 PointToWire(left, LWire);
187 PointToWire(right, RWire);
189 ll.Add(LWire.Oriented(TopAbs_FORWARD));
191 for (int k = 0; k < main_profiles.size(); k++)
194 PointToWire(main_profiles[k], W);
195 TopAbs_Orientation Ori = TopAbs_INTERNAL;
196 if (k == 0 || k == main_profiles.size() - 1)
197 Ori = TopAbs_FORWARD;
198 ll.Add(W.Oriented(Ori));
202 ll.Add(RWire.Oriented(TopAbs_FORWARD));
203 //yes, add subshapes in this order (left + profiles + right)
204 //otherwise the projected wire will be non-manifold
210 void HYDROData_DTM::Update()
212 AltitudePoints points;
213 TopoDS_Shape Out3dPres;
214 TopoDS_Shape Out2dPres;
215 TopoDS_Shape OutLeftB;
216 TopoDS_Shape OutRightB;
217 TopoDS_Shape OutInlet;
218 TopoDS_Shape OutOutlet;
220 HYDROData_SequenceOfObjects objs = GetProfiles();
221 double ddz = GetDDZ();
222 double step = GetSpatialStep();
223 std::set<int> InvInd;
224 bool WireIntersections; //__TODO
225 CreateProfilesFromDTM( objs, ddz, step, points, Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, true, true, InvInd, -1, WireIntersections );
226 SetAltitudePoints( points );
228 SetShape( DataTag_LeftBankShape, OutLeftB);
229 SetShape( DataTag_RightBankShape, OutRightB);
230 SetShape( DataTag_InletShape, OutInlet);
231 SetShape( DataTag_OutletShape, OutOutlet );
232 SetShape( DataTag_3DShape, Out3dPres );
233 SetShape( DataTag_2DShape, Out2dPres );
235 HYDROData_Bathymetry::Update();
238 void HYDROData_DTM::GetPresentationShapes( TopoDS_Shape& Out3dPres,
239 TopoDS_Shape& Out2dPres,
240 TopoDS_Shape& OutLeftB,
241 TopoDS_Shape& OutRightB,
242 TopoDS_Shape& OutInlet,
243 TopoDS_Shape& OutOutlet )
246 OutLeftB = GetShape( DataTag_LeftBankShape);
247 OutRightB = GetShape( DataTag_RightBankShape);
248 OutInlet = GetShape( DataTag_InletShape);
249 OutOutlet = GetShape( DataTag_OutletShape );
250 Out3dPres = GetShape( DataTag_3DShape );
251 Out2dPres = GetShape( DataTag_2DShape );
253 void HYDROData_DTM::CreateProfilesFromDTM (const HYDROData_SequenceOfObjects& InpProfiles,
256 AltitudePoints& points,
257 TopoDS_Shape& Out3dPres,
258 TopoDS_Shape& Out2dPres,
259 TopoDS_Shape& OutLeftB,
260 TopoDS_Shape& OutRightB,
261 TopoDS_Shape& OutInlet,
262 TopoDS_Shape& OutOutlet,
265 std::set<int>& InvInd,
267 bool& WireIntersections)
269 int aLower = InpProfiles.Lower(), anUpper = InpProfiles.Upper();
270 size_t n = anUpper - aLower + 1;
272 std::vector<Handle(HYDROData_Profile)> profiles;
273 profiles.reserve( n );
274 for( int i=aLower; i<=anUpper; i++ )
276 Handle(HYDROData_Profile) aProfile = Handle(HYDROData_Profile)::DownCast( InpProfiles.Value( i ) );
277 if( !aProfile.IsNull() )
278 profiles.push_back( aProfile );
280 const double EPS = 1E-3;
282 AltitudePoints right;
283 std::vector<AltitudePoints> main_profiles;
285 if( thePntsLimit > 0 )
287 int aNbPoints = EstimateNbPoints( profiles, ddz, step );
288 if( aNbPoints < 0 || aNbPoints > thePntsLimit )
292 if( ddz>EPS && step>EPS )
293 CreateProfiles(profiles, ddz, step, left, right, points, main_profiles,
294 Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, Create3dPres, Create2dPres, InvInd, WireIntersections );
297 void HYDROData_DTM::ProjWireOnPlane(const TopoDS_Shape& inpWire, const Handle(Geom_Plane)& RefPlane,
298 TopTools_DataMapOfShapeListOfShape* E2PE)
302 //project shape (edges) on planar face
304 BB.MakeFace(F, RefPlane, Precision::Confusion());
305 BRepAlgo_NormalProjection nproj(F);
307 nproj.SetDefaultParams();
312 TopoDS_Shape projRes = nproj.Projection();
314 // unite all vertexes/edges from projected result
315 BOPAlgo_Builder anAlgo;
316 TopExp_Explorer exp(projRes, TopAbs_EDGE);
317 for (;exp.More(); exp.Next())
319 const TopoDS_Edge& E = TopoDS::Edge(exp.Current());
320 if (E.Orientation() != TopAbs_INTERNAL)
321 anAlgo.AddArgument(E);
324 int stat = anAlgo.ErrorStatus();
325 TopoDS_Shape projResConn = anAlgo.Shape();
327 // make wire => vertexes and edges should be untouched after this operation!
328 exp.Init(projResConn, TopAbs_EDGE);
329 //TopTools_ListOfShape llE;
331 //BRepLib_MakeWire WM;
333 //for (;exp.More();exp.Next())
334 // llE.Append(exp.Current());
337 //outWire = WM.Wire();
339 //outWire.Orientation(inpWire.Orientation()); //take from the original wire
341 //history mode: edge to projected edges
344 TopExp_Explorer ex(inpWire, TopAbs_EDGE);
345 for (;ex.More();ex.Next())
347 const TopoDS_Edge& CE = TopoDS::Edge(ex.Current());
348 TopTools_ListOfShape NEL;
349 const TopTools_ListOfShape& LS = nproj.Generated(CE);
350 TopTools_ListIteratorOfListOfShape it(LS);
351 for (;it.More();it.Next())
353 const TopoDS_Shape& PCE = it.Value();
354 TopTools_ListOfShape PLS = anAlgo.Modified(PCE);
357 TopTools_ListIteratorOfListOfShape itp(PLS);
358 for (;itp.More();itp.Next())
359 NEL.Append(itp.Value());
366 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
367 bool HYDROData_DTM::Get2dFaceFrom3dPres(const TopoDS_Compound& cmp, TopoDS_Face& outF,
368 TopTools_SequenceOfShape* Boundr, std::set<int> ind )
370 //ind : set of indices (starts with 0). index == number of boundary (inlet, outlet, etc..)
372 //if Boundr is not null => this method will return sequence of boundary wires (inlet, outlet...)
374 Handle(Geom_Plane) refpl = new Geom_Plane(gp_Pnt(0,0,0), gp_Dir(0,0,1));
375 TopTools_DataMapOfShapeListOfShape E2PE;
376 ProjWireOnPlane(cmp, refpl, &E2PE);
377 TopTools_ListOfShape ELL;
379 TopoDS_Iterator it(cmp);
381 for (;it.More(); it.Next())
383 const TopoDS_Wire& W = TopoDS::Wire(it.Value());
384 if (W.Orientation() != TopAbs_INTERNAL)
387 TopExp_Explorer ex(W, TopAbs_EDGE);
388 TopTools_ListOfShape LEpW;
389 TopTools_ListOfShape LEpWDummy;
390 for (;ex.More();ex.Next())
392 const TopoDS_Edge& CE = TopoDS::Edge(ex.Current());
393 TopTools_ListOfShape LS = E2PE.Find(CE);
397 if (ind.count(i) != 0)
401 const TopoDS_Wire& WMW = WM.Wire();
402 //assume that wire is a straight line,
403 //take first and last vertex and make simple edge (RE)
404 TopoDS_Vertex VF, VL;
405 TopExp::Vertices(WMW, VF, VL);
406 TopoDS_Edge RE = BRepLib_MakeEdge(VF, VL).Edge();
409 LEpWDummy = LEpW; //LEpW will be nullified after appending to ELL
414 LEpWDummy.Append(RE);
417 //TODO: in the new version of OCCT, USD can process separate wires
418 //ShapeUpgrade_UnifySameDomain USD(WMW, 1U, 0U, 1U); //concat bsplines
420 //const TopoDS_Shape& RSU = USD.Shape();
421 //TopExp_Explorer exp(RSU, TopAbs_EDGE);
422 //TopTools_ListOfShape DummyL;
423 //for (;it.More();it.Next())
424 // DummyL.Append(it.Value());
425 //if (DummyL.Extent() == 1)
426 // ELL.Append(DummyL.First()); //if one edge => accept this result
428 // ELL.Append(LEpW); //else put 'as is'
438 //make inlet, outlet, left/tight banks [wires]
439 //shouldn't change topology of the edges
440 BRepLib_MakeWire IWM;
442 Boundr->Append(IWM.Wire());
449 BRepLib_MakeWire WME;
452 const TopoDS_Wire& resW = WME.Wire();
453 BRepBuilderAPI_MakeFace mf(refpl, resW, true); //check inside is true by def
456 ShapeAnalysis_Wire WA(resW, outF, Precision::Confusion());
457 bool res = WA.CheckSelfIntersection(); //TODO check that this return false if OK
460 ///!!! the internal wires cant be added with 'internal' ori.
461 // it's possible to do with brep builder yet the result will not be correct!
462 // more proper way is to use BOP operation here.
463 /*for (int i = 1; i <= IntW.Extent(); i++)
466 const TopoDS_Wire& W = TopoDS::Wire(IntW(i));
467 ProjWireOnPlane(W, refpl, outIW);
472 void HYDROData_DTM::CreateProfiles(const std::vector<Handle(HYDROData_Profile)>& theProfiles,
474 double theSpatialStep,
475 AltitudePoints& theOutLeft,
476 AltitudePoints& theOutRight,
477 AltitudePoints& theOutPoints,
478 std::vector<AltitudePoints>& theOutMainProfiles,
479 TopoDS_Shape& Out3dPres,
480 TopoDS_Shape& Out2dPres,
481 TopoDS_Shape& OutLeftB,
482 TopoDS_Shape& OutRightB,
483 TopoDS_Shape& OutInlet,
484 TopoDS_Shape& OutOutlet,
487 std::set<int>& InvInd,
488 bool& WireIntersections)
490 if (theProfiles.empty())
492 theOutPoints = Interpolate( theProfiles, theDDZ, theSpatialStep, theOutLeft, theOutRight, theOutMainProfiles, InvInd );
493 //note that if Create3dPres is false => Create2dPres flag is meaningless!
494 if( theOutPoints.empty() )
499 TopTools_IndexedMapOfOrientedShape ll = Create3DShape( theOutLeft, theOutRight, theOutMainProfiles);
505 BB.MakeCompound(cmp);
506 for (int i = 1; i <= ll.Extent(); i++)
511 //same order as in HYDROData_DTM::Update()
513 OutRightB = ll(ll.Extent());
515 OutOutlet = ll(ll.Extent() - 1);
520 WireIntersections = Get2dFaceFrom3dPres(cmp, outF); //__TODO
529 void HYDROData_DTM::GetProperties( const Handle(HYDROData_Profile)& theProfile,
530 gp_Pnt& theLowestPoint, gp_Vec2d& theDir,
531 double& theZMin, double& theZMax )
533 theLowestPoint = theProfile->GetBottomPoint();
536 theProfile->GetLeftPoint( aLeft, true, true );
537 theProfile->GetRightPoint( aRight, true, true );
538 double x = aRight.X()-aLeft.X();
539 double y = aRight.Y()-aLeft.Y();
540 theDir = gp_Vec2d( x, y );
542 HYDROData_Profile::ProfilePoints points = theProfile->GetProfilePoints();
543 int lo = points.Lower();
544 int up = points.Upper();
545 theZMin = std::numeric_limits<double>::max();
547 for( int i=lo; i<=up; i++ )
549 double z = points.Value( i ).Z();
557 inline gp_Pnt2d To2D( const gp_Pnt& thePnt, const gp_Trsf& theTr,
558 double& theUMin, double& theUMax )
560 gp_Pnt p = thePnt.Transformed( theTr );
567 return gp_Pnt2d( u, z );
570 Handle(TColgp_HArray1OfPnt2d) To2D( const TColgp_Array1OfPnt& thePoints,
571 const gp_Trsf& theTr,
572 double& theUMin, double& theUMax )
574 int low = thePoints.Lower(), up = thePoints.Upper();
575 Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( low, up );
576 for( int i=low; i<=up; i++ )
577 points->SetValue( i, To2D( thePoints.Value( i ), theTr, theUMin, theUMax ) );
581 Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve,
582 Standard_Real theFirst, Standard_Real theLast,
583 const gp_Trsf& theTr,
584 double& theUMin, double& theUMax )
586 if( theCurve->IsKind( STANDARD_TYPE( Geom_Line ) ) )
588 gp_Pnt aFirstPnt, aLastPnt;
589 theCurve->D0( theFirst, aFirstPnt );
590 theCurve->D0( theLast, aLastPnt );
593 aFirst2d = To2D( aFirstPnt, theTr, theUMin, theUMax ),
594 aLast2d = To2D( aLastPnt, theTr, theUMin, theUMax );
596 gp_Vec2d dir( aFirst2d, aLast2d );
597 Handle(Geom2d_Line) aLine2d = new Geom2d_Line( aFirst2d, gp_Dir2d( dir.X(), dir.Y() ) );
598 return new Geom2d_TrimmedCurve( aLine2d, 0, aLast2d.Distance( aFirst2d ) );
601 if( theCurve->IsKind( STANDARD_TYPE( Geom_BSplineCurve ) ) )
603 Handle(Geom_BSplineCurve) aSpline = Handle(Geom_BSplineCurve)::DownCast( theCurve );
605 Handle(TColgp_HArray1OfPnt2d) poles = To2D( aSpline->Poles(), theTr, theUMin, theUMax );
606 const TColStd_Array1OfReal& knots = aSpline->Knots();
607 const TColStd_Array1OfInteger& multiplicities = aSpline->Multiplicities();
608 int aDegree = aSpline->Degree();
610 return new Geom2d_BSplineCurve( poles->Array1(), knots, multiplicities, aDegree );
613 return Handle(Geom2d_Curve)();
616 #include <GCE2d_MakeSegment.hxx>
617 #include <Geom2dAPI_InterCurveCurve.hxx>
618 bool IsCooriented( const Handle(HYDROData_Profile)& theProfile1,
619 const Handle(HYDROData_Profile)& theProfile2 )
621 if( theProfile1==theProfile2 )
624 gp_XY lp1, rp1, lp2, rp2;
625 theProfile1->GetLeftPoint(lp1);
626 theProfile1->GetRightPoint(rp1);
627 theProfile2->GetLeftPoint(lp2);
628 theProfile2->GetRightPoint(rp2);
630 GCE2d_MakeSegment s1(lp1, lp2);
631 GCE2d_MakeSegment s2(rp1, rp2);
633 Geom2dAPI_InterCurveCurve inter;
634 inter.Init(s1.Value(), s2.Value());
635 if (inter.NbPoints() == 0)
641 Handle(Geom2d_BSplineCurve) HYDROData_DTM::CreateHydraulicAxis(
642 const std::vector<Handle(HYDROData_Profile)>& theProfiles,
643 std::vector<double>& theDistances )
645 size_t n = theProfiles.size();
647 return Handle_Geom2d_BSplineCurve();
649 Handle_Geom2d_BSplineCurve aResult;
651 Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( 1, (int)n );
652 TColgp_Array1OfVec2d tangents( 1, (int)n );
653 Handle(TColStd_HArray1OfBoolean) flags = new TColStd_HArray1OfBoolean( 1, (int)n );
655 for( size_t i = 1; i <= n; i++ )
657 Handle(HYDROData_Profile) aProfile = theProfiles[i-1];
658 Handle(HYDROData_Profile) aPrevProfile = i==1 ? theProfiles[i-1] : theProfiles[i-2];
659 Handle(HYDROData_Profile) aNextProfile = i==n ? theProfiles[i-1] : theProfiles[i];
661 if( !IsCooriented( aProfile, aNextProfile ) )
664 aProfile->GetLeftPoint( lp, true );
665 aProfile->GetRightPoint( rp, true );
666 aProfile->SetLeftPoint( rp, true );
667 aProfile->SetRightPoint( lp, true );
675 gp_XYZ curP = aProfile->GetBottomPoint();
676 gp_XY curP2d = gp_XY(curP.X(), curP.Y());
680 nextP = aPrevProfile->GetBottomPoint(true);
682 nextP = aNextProfile->GetBottomPoint(true);
684 gp_XY nextP2d = gp_XY(nextP.X(), nextP.Y());
687 GetProperties( aProfile, aLowest, aPrTangent, zmin, zmax );
688 aNormal.SetCoord( -aPrTangent.Y(), aPrTangent.X() );
690 gp_Vec2d aDirToNextProfile(nextP2d.X() - curP2d.X(), nextP2d.Y() - curP2d.Y() );
692 aDirToNextProfile.Reverse();
693 if (aNormal.Dot(aDirToNextProfile) < 0)
698 points->SetValue( (int)i, gp_Pnt2d( aLowest.X(), aLowest.Y() ) );
699 tangents.SetValue( (int)i, aNormal );
700 flags->SetValue( (int)i, Standard_True );
703 Geom2dAPI_Interpolate anInterpolator( points, Standard_False, Standard_False );
704 anInterpolator.Load( tangents, flags );
705 anInterpolator.Perform();
706 if( anInterpolator.IsDone() )
708 aResult = anInterpolator.Curve();
710 //fill the distances vector
711 Geom2dAdaptor_Curve anAdaptor( aResult );
713 theDistances.clear();
714 theDistances.reserve( n );
715 Standard_Real aParamFirst = anAdaptor.FirstParameter(), aParamLast = anAdaptor.LastParameter();
716 for( size_t i = 1; i <= n; i++ )
718 gp_Pnt2d aPnt = points->Value( (Standard_Integer)i );
719 Geom2dAPI_ProjectPointOnCurve aProject( aPnt, aResult );
720 Standard_Real aParam = aProject.LowerDistanceParameter();
721 double aDistance = GCPnts_AbscissaPoint::Length( anAdaptor, aParamFirst, aParam );
722 theDistances.push_back( aDistance );
728 std::vector<Handle(Geom2d_Curve)> HYDROData_DTM::ProfileToParametric(
729 const Handle(HYDROData_Profile)& theProfile,
730 double& theUMin, double& theUMax, gp_Vec2d& theDir )
732 std::vector<Handle(Geom2d_Curve)> curves;
734 // Transformation of the coordinate systems
737 GetProperties( theProfile, aLowest, theDir, zmin, zmax );
739 gp_Ax3 aStd3d( gp_Pnt( 0, 0, 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( 1, 0, 0 ) );
740 gp_Ax3 aLocal( gp_Pnt( aLowest.X(), aLowest.Y(), 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( theDir.X(), theDir.Y(), 0 ) );
743 aTransf.SetTransformation( aStd3d, aLocal );
745 // Iteration via edges
746 TopoDS_Wire aWire = TopoDS::Wire( theProfile->GetShape3D() );
747 TopExp_Explorer anExp( aWire, TopAbs_EDGE );
748 for( ; anExp.More(); anExp.Next() )
750 // Extract an edge from wire
751 TopoDS_Edge anEdge = TopoDS::Edge( anExp.Current() );
753 // Extract a curve corresponding to the edge
754 TopLoc_Location aLoc;
755 Standard_Real aFirst, aLast;
756 Handle(Geom_Curve) aCurve = BRep_Tool::Curve( anEdge, aLoc, aFirst, aLast );
758 // Convert the curve to 2d CS
759 Handle(Geom2d_Curve) aCurve2d = CurveTo2D( aCurve, aFirst, aLast, aTransf, theUMin, theUMax );
760 if( !aCurve2d.IsNull() )
761 curves.push_back( aCurve2d );
767 bool CalcMidWidth( const std::set<double>& intersections, double& theMid, double& theWid )
769 double umin = std::numeric_limits<double>::max(),
772 size_t n = intersections.size();
776 std::set<double>::const_iterator it = intersections.begin(), last = intersections.end();
777 for( ; it!=last; it++ )
785 theMid = ( umin+umax )/2;
790 void HYDROData_DTM::ProfileDiscretization( const Handle(HYDROData_Profile)& theProfile,
791 double theXCurv, double theMinZ, double theMaxZ, double theDDZ,
792 CurveUZ& theMidPointCurve,
793 CurveUZ& theWidthCurve,
794 int& intersection_nb,
797 double aDblMax = std::numeric_limits<double>::max(),
802 gp_Vec2d aProfileDir;
803 std::vector<Handle(Geom2d_Curve)> curves = ProfileToParametric( theProfile, aUMin, aUMax, aProfileDir );
804 size_t n = curves.size();
809 // we add the "virtual" vertical lines to simulate the intersection with profile
810 gp_Pnt2d aFirst, aLast;
811 curves[0]->D0( curves[0]->FirstParameter(), aFirst );
812 curves[n-1]->D0( curves[n-1]->LastParameter(), aLast );
813 Handle(Geom2d_Line) aV1 = new Geom2d_Line( aFirst, gp_Dir2d( 0, 1 ) );
814 Handle(Geom2d_TrimmedCurve) aT1 = new Geom2d_TrimmedCurve( aV1, 0.0, aVMax );
816 Handle(Geom2d_Line) aV2 = new Geom2d_Line( aLast, gp_Dir2d( 0, 1 ) );
817 Handle(Geom2d_TrimmedCurve) aT2 = new Geom2d_TrimmedCurve( aV2, 0.0, aVMax );
819 curves.push_back( aT1 );
820 curves.push_back( aT2 );
822 int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 );
823 theMidPointCurve = CurveUZ( theXCurv, aProfileDir, theMinZ );
824 theMidPointCurve.reserve( psize );
825 theWidthCurve = CurveUZ( theXCurv, aProfileDir, theMinZ );
826 theWidthCurve.reserve( psize );
829 // for each discrete value of z we search intersection with profile
830 for( double z1 = theMinZ; z1 <= theMaxZ; z1 += theDDZ )
832 Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z1 ), gp_Dir2d( 1, 0 ) );
833 std::set<double> intersections;
834 for( size_t i = 0; i < n; i++ )
836 Handle(Geom2d_Curve) aCurve = curves[i];
837 Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
838 for( int k=1, m=anIntersect.NbPoints(); k<=m; k++ )
839 intersections.insert( anIntersect.Point( k ).X() );
842 intersection_nb = intersections.size();
843 if( intersection_nb >= 1 )
846 if( !CalcMidWidth( intersections, u_mid, u_wid ) )
849 double z = z1 - theMinZ;
853 theMidPointCurve.push_back( p_mid );
858 theWidthCurve.push_back( p_wid );
863 void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCurveB,
864 int theNbSteps, std::vector<CurveUZ>& theInterpolation,
867 theInterpolation.clear();
868 int d = isAddSecond ? 2 : 1;
869 theInterpolation.reserve( theNbSteps+d );
870 double dt = 1.0 / double( theNbSteps + 1 );
872 theInterpolation.push_back( theCurveA );
873 for( int i=0; i<theNbSteps; i++, t+=dt )
875 CurveUZ anInterp = theCurveA*(1-t) + theCurveB*t;
876 theInterpolation.push_back( anInterp );
879 theInterpolation.push_back( theCurveB );
881 #include <BRepLib_MakeEdge2d.hxx>
882 void HYDROData_DTM::CurveTo3D( const Handle(Geom2d_BSplineCurve)& theHydraulicAxis,
883 const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
884 AltitudePoints& thePoints )
886 Geom2dAdaptor_Curve anAdaptor( theHydraulicAxis );
887 TopoDS_Edge E2d = BRepLib_MakeEdge2d(theHydraulicAxis).Edge();
888 GCPnts_AbscissaPoint ap( anAdaptor, theMidCurve.Xcurv(), anAdaptor.FirstParameter() );
889 double aParam = ap.Parameter();
892 anAdaptor.D0( aParam, point );
893 gp_Vec2d profile_dir = theMidCurve.ProfileDir();
894 //gp_Dir tangent_n( -profile_dir.Y(), profile_dir.X(), dz );
895 profile_dir.Normalize();
897 size_t n = theMidCurve.size();
898 std::map<double, AltitudePoint> sorted_points;
899 for( size_t i=0; i<n; i++ )
901 double param1 = theMidCurve[i].U - theWidthCurve[i].U / 2;
902 double param2 = theMidCurve[i].U + theWidthCurve[i].U / 2;
904 gp_Pnt2d p1 = point.Translated( param1 * profile_dir);
905 gp_Pnt2d p2 = point.Translated( param2 * profile_dir);
907 double z = theMidCurve[i].Z + theMidCurve.DeltaZ();
909 AltitudePoint p3d_1( p1.X(), p1.Y(), z ), p3d_2( p2.X(), p2.Y(), z );
911 sorted_points[param1] = p3d_1;
912 sorted_points[param2] = p3d_2;
915 thePoints.reserve( sorted_points.size() );
916 const double EPS = 1E-12;
917 std::map<double, AltitudePoint>::const_iterator it = sorted_points.begin(), last = sorted_points.end();
918 for( ; it!=last; it++ )
919 if( thePoints.empty() || thePoints.back().SquareDistance( it->second ) > EPS )
920 thePoints.push_back( it->second );
923 inline double max( double a, double b )
931 inline double min( double a, double b )
939 #include <BRepLib_MakeWire.hxx>
941 std::vector<HYDROData_Bathymetry::AltitudePoints> HYDROData_DTM::Interpolate
942 ( const Handle(Geom2d_BSplineCurve)& theHydraulicAxis,
943 const Handle(HYDROData_Profile)& theProfileA,
945 const Handle(HYDROData_Profile)& theProfileB,
947 double theDDZ, int theNbSteps, bool isAddSecond,
948 int& inter_nb_1, int& inter_nb_2)
950 double zminA, zmaxA, zminB, zmaxB;
951 gp_Pnt lowestA, lowestB;
954 GetProperties( theProfileA, lowestA, dirA, zminA, zmaxA );
955 GetProperties( theProfileB, lowestB, dirB, zminB, zmaxB );
958 double hmax = max( zmaxA-zminA, zmaxB-zminB );
960 //double dz = zminB - zminA;
961 //double zmin = min( zminA, zminB );
962 //double zmax = max( zmaxA, zmaxB );
964 CurveUZ midA(0, gp_Vec2d(), 0), midB(0, gp_Vec2d(), 0);
965 CurveUZ widA(0, gp_Vec2d(), 0), widB(0, gp_Vec2d(), 0);
967 ProfileDiscretization( theProfileA, theXCurvA, zminA, zminA+hmax, theDDZ, midA, widA, inter_nb_1 );
968 ProfileDiscretization( theProfileB, theXCurvB, zminB, zminB+hmax, theDDZ, midB, widB, inter_nb_2 );
970 std::vector<CurveUZ> mid, wid;
971 Interpolate( midA, midB, theNbSteps, mid, isAddSecond );
972 Interpolate( widA, widB, theNbSteps, wid, isAddSecond );
974 size_t p = mid.size();
975 size_t q = p>0 ? 2*mid[0].size() : 1;
976 std::vector<AltitudePoints> points;
979 for( size_t i=0; i<p; i++ )
981 points[i].reserve( q );
982 CurveTo3D( theHydraulicAxis, mid[i], wid[i], points[i] );
988 HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
989 ( const std::vector<Handle(HYDROData_Profile)>& theProfiles,
990 double theDDZ, double theSpatialStep,
991 AltitudePoints& theLeft,
992 AltitudePoints& theRight,
993 std::vector<AltitudePoints>& theMainProfiles,
994 std::set<int>& invalInd)
996 AltitudePoints points;
997 size_t n = theProfiles.size();
1001 std::vector<double> distances;
1002 Handle(Geom2d_BSplineCurve) aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
1003 if( aHydraulicAxis.IsNull() )
1006 theMainProfiles.reserve( n );
1008 for( size_t i=0, n1=n-1; i<n1; i++ )
1010 double aDistance = distances[i+1]-distances[i];
1011 int aNbSteps = int(aDistance/theSpatialStep);
1012 bool isAddSecond = i==n1-1;
1014 // 1. Calculate interpolated profiles
1015 int inter_nb_1, inter_nb_2;
1016 std::vector<AltitudePoints> local_points = Interpolate( aHydraulicAxis, theProfiles[i], distances[i],
1017 theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond, inter_nb_1, inter_nb_2 );
1018 int lps = local_points.size();
1024 invalInd.insert(i+1);
1026 // 2. Put all points into the global container
1027 for( size_t j=0; j<lps; j++ )
1029 const AltitudePoints& lp = local_points[j];
1031 points.reserve( lp.size() * n );
1032 for( size_t k=0, ks=lp.size(); k<ks; k++ )
1033 points.push_back( lp[k] );
1036 // 3. Get left/right banks' points
1039 theLeft.reserve( lps * n );
1040 theRight.reserve( lps * n );
1042 for( size_t j=0; j<lps; j++ )
1044 const AltitudePoints& lp = local_points[j];
1045 theLeft.push_back( lp[0] );
1046 theRight.push_back( lp[lp.size()-1] );
1049 // 4. Get main profiles points
1050 theMainProfiles.push_back( local_points[0] );
1052 theMainProfiles.push_back( local_points[lps-1] );
1057 int HYDROData_DTM::EstimateNbPoints( const std::vector<Handle(HYDROData_Profile)>& theProfiles,
1058 double theDDZ, double theSpatialStep )
1060 size_t n = theProfiles.size();
1063 if( theDDZ<1E-6 || theSpatialStep<1E-6 )
1066 std::vector<double> distances;
1067 Handle(Geom2d_BSplineCurve) aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
1068 if( aHydraulicAxis.IsNull() )
1071 double aCompleteDistance = distances[n-1];
1072 int aNbSteps = int( aCompleteDistance / theSpatialStep ) + 1;
1075 double aZMin, aZMax;
1076 GetProperties( theProfiles[0], aLowest, aDir, aZMin, aZMax );
1077 int aNbZSteps = (aZMax-aZMin)/theDDZ;
1079 if( aNbSteps > ( 1<<16 ) || aNbZSteps > ( 1<<16 ) )
1082 return aNbSteps * aNbZSteps;