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 <Geom_Plane.hxx>
40 #include <BRepTools_WireExplorer.hxx>
41 #include <TopTools_IndexedMapOfShape.hxx>
42 #include <BRepBuilderAPI_MakeFace.hxx>
44 #include <TopTools_IndexedMapOfOrientedShape.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>
55 #include <TopTools_IndexedDataMapOfShapeShape.hxx>
56 #include <TopTools_ListIteratorOfListOfShape.hxx>
57 #include <TopTools_SequenceOfShape.hxx>
59 IMPLEMENT_STANDARD_HANDLE( HYDROData_DTM, HYDROData_Bathymetry )
60 IMPLEMENT_STANDARD_RTTIEXT( HYDROData_DTM, HYDROData_Bathymetry )
62 HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv, const gp_Vec2d& theProfileDir, double theDeltaZ )
63 : myXcurv( theXCurv ), myProfileDir( theProfileDir ), myDeltaZ( theDeltaZ )
67 HYDROData_DTM::CurveUZ::~CurveUZ()
71 double HYDROData_DTM::CurveUZ::Xcurv() const
76 gp_Vec2d HYDROData_DTM::CurveUZ::ProfileDir() const
81 double HYDROData_DTM::CurveUZ::DeltaZ() const
86 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator + ( const CurveUZ& c ) const
88 HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv(), ProfileDir() + c.ProfileDir(), DeltaZ() + c.DeltaZ() );
89 size_t n = size(), n1 = c.size();
92 std::cout << "Warning: different number of points in curves: " << n << ", " << n1 << std::endl;
95 for( int i=0; i<n; i++ )
98 p.U = operator[]( i ).U + c[i].U;
99 p.Z = operator[]( i ).Z;
105 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator * ( double d ) const
107 HYDROData_DTM::CurveUZ res( Xcurv()*d, ProfileDir()*d, DeltaZ()*d );
110 for( int i=0; i<n; i++ )
113 p.U = operator[]( i ).U * d;
114 p.Z = operator[]( i ).Z;
123 HYDROData_DTM::HYDROData_DTM()
127 HYDROData_DTM::~HYDROData_DTM()
131 HYDROData_SequenceOfObjects HYDROData_DTM::GetProfiles() const
133 return GetReferenceObjects( DataTag_Profiles );
136 void HYDROData_DTM::SetProfiles( const HYDROData_SequenceOfObjects& theProfiles )
138 SetReferenceObjects( theProfiles, DataTag_Profiles );
142 double HYDROData_DTM::GetDDZ() const
144 return GetDouble( DataTag_DDZ );
147 void HYDROData_DTM::SetDDZ( double theDDZ )
149 SetDouble( DataTag_DDZ, theDDZ );
153 double HYDROData_DTM::GetSpatialStep() const
155 return GetDouble( DataTag_SpatialStep );
158 void HYDROData_DTM::SetSpatialStep( double theSpatialStep )
160 SetDouble( DataTag_SpatialStep, theSpatialStep );
164 void HYDROData_DTM::PointToWire(const AltitudePoints& pnts, TopoDS_Wire& W )
166 /*BRepLib_MakeWire WM;
169 for (int i = 0; i < pnts.size() - 1; i++)
171 gp_Pnt p1(pnts[i].X, pnts[i].Y, pnts[i].Z);
172 gp_Pnt p2(pnts[i+1].X, pnts[i+1].Y, pnts[i+1].Z);
173 WM.Add(BRepLib_MakeEdge(p1, p2).Edge());
177 BRepBuilderAPI_MakePolygon PM;
178 for (int i = 0; i < pnts.size(); i++)
179 PM.Add(gp_Pnt(pnts[i].X, pnts[i].Y, pnts[i].Z));
184 TopTools_IndexedMapOfOrientedShape HYDROData_DTM::Create3DShape(const AltitudePoints& left,
185 const AltitudePoints& right,
186 const std::vector<AltitudePoints>& main_profiles)
188 TopTools_IndexedMapOfOrientedShape ll;
189 TopoDS_Wire LWire, RWire;
190 PointToWire(left, LWire);
191 PointToWire(right, RWire);
193 ll.Add(LWire.Oriented(TopAbs_FORWARD));
195 for (int k = 0; k < main_profiles.size(); k++)
198 PointToWire(main_profiles[k], W);
199 TopAbs_Orientation Ori = TopAbs_INTERNAL;
200 if (k == 0 || k == main_profiles.size() - 1)
201 Ori = TopAbs_FORWARD;
202 ll.Add(W.Oriented(Ori));
206 ll.Add(RWire.Oriented(TopAbs_FORWARD));
207 //yes, add subshapes in this order (left + profiles + right)
208 //otherwise the projected wire will be non-manifold
214 void HYDROData_DTM::Update()
216 AltitudePoints points;
217 TopoDS_Shape Out3dPres;
218 TopoDS_Shape Out2dPres;
219 TopoDS_Shape OutLeftB;
220 TopoDS_Shape OutRightB;
221 TopoDS_Shape OutInlet;
222 TopoDS_Shape OutOutlet;
224 HYDROData_SequenceOfObjects objs = GetProfiles();
225 double ddz = GetDDZ();
226 double step = GetSpatialStep();
227 std::set<int> InvInd;
228 bool WireIntersections; //__TODO
229 CreateProfilesFromDTM( objs, ddz, step, points, Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, true, true, InvInd, -1, WireIntersections );
230 SetAltitudePoints( points );
232 SetShape( DataTag_LeftBankShape, OutLeftB);
233 SetShape( DataTag_RightBankShape, OutRightB);
234 SetShape( DataTag_InletShape, OutInlet);
235 SetShape( DataTag_OutletShape, OutOutlet );
236 SetShape( DataTag_3DShape, Out3dPres );
237 SetShape( DataTag_2DShape, Out2dPres );
239 HYDROData_Bathymetry::Update();
242 void HYDROData_DTM::GetPresentationShapes( TopoDS_Shape& Out3dPres,
243 TopoDS_Shape& Out2dPres,
244 TopoDS_Shape& OutLeftB,
245 TopoDS_Shape& OutRightB,
246 TopoDS_Shape& OutInlet,
247 TopoDS_Shape& OutOutlet )
250 OutLeftB = GetShape( DataTag_LeftBankShape);
251 OutRightB = GetShape( DataTag_RightBankShape);
252 OutInlet = GetShape( DataTag_InletShape);
253 OutOutlet = GetShape( DataTag_OutletShape );
254 Out3dPres = GetShape( DataTag_3DShape );
255 Out2dPres = GetShape( DataTag_2DShape );
257 void HYDROData_DTM::CreateProfilesFromDTM (const HYDROData_SequenceOfObjects& InpProfiles,
260 AltitudePoints& points,
261 TopoDS_Shape& Out3dPres,
262 TopoDS_Shape& Out2dPres,
263 TopoDS_Shape& OutLeftB,
264 TopoDS_Shape& OutRightB,
265 TopoDS_Shape& OutInlet,
266 TopoDS_Shape& OutOutlet,
269 std::set<int>& InvInd,
271 bool& WireIntersections)
273 int aLower = InpProfiles.Lower(), anUpper = InpProfiles.Upper();
274 size_t n = anUpper - aLower + 1;
276 std::vector<Handle_HYDROData_Profile> profiles;
277 profiles.reserve( n );
278 for( int i=aLower; i<=anUpper; i++ )
280 Handle(HYDROData_Profile) aProfile = Handle(HYDROData_Profile)::DownCast( InpProfiles.Value( i ) );
281 if( !aProfile.IsNull() )
282 profiles.push_back( aProfile );
284 const double EPS = 1E-3;
286 AltitudePoints right;
287 std::vector<AltitudePoints> main_profiles;
289 if( thePntsLimit > 0 )
291 int aNbPoints = EstimateNbPoints( profiles, ddz, step );
292 if( aNbPoints < 0 || aNbPoints > thePntsLimit )
296 if( ddz>EPS && step>EPS )
297 CreateProfiles(profiles, ddz, step, left, right, points, main_profiles,
298 Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, Create3dPres, Create2dPres, InvInd, WireIntersections );
301 void HYDROData_DTM::ProjWireOnPlane(const TopoDS_Shape& inpWire, const Handle_Geom_Plane& RefPlane,
302 TopTools_DataMapOfShapeListOfShape* E2PE)
306 //project shape (edges) on planar face
308 BB.MakeFace(F, RefPlane, Precision::Confusion());
309 BRepAlgo_NormalProjection nproj(F);
311 nproj.SetDefaultParams();
316 TopoDS_Shape projRes = nproj.Projection();
318 // unite all vertexes/edges from projected result
319 BOPAlgo_Builder anAlgo;
320 TopExp_Explorer exp(projRes, TopAbs_EDGE);
321 for (;exp.More(); exp.Next())
323 const TopoDS_Edge& E = TopoDS::Edge(exp.Current());
324 if (E.Orientation() != TopAbs_INTERNAL)
325 anAlgo.AddArgument(E);
328 int stat = anAlgo.ErrorStatus();
329 TopoDS_Shape projResConn = anAlgo.Shape();
331 // make wire => vertexes and edges should be untouched after this operation!
332 exp.Init(projResConn, TopAbs_EDGE);
333 //TopTools_ListOfShape llE;
335 //BRepLib_MakeWire WM;
337 //for (;exp.More();exp.Next())
338 // llE.Append(exp.Current());
341 //outWire = WM.Wire();
343 //outWire.Orientation(inpWire.Orientation()); //take from the original wire
345 //history mode: edge to projected edges
348 TopExp_Explorer ex(inpWire, TopAbs_EDGE);
349 for (;ex.More();ex.Next())
351 const TopoDS_Edge& CE = TopoDS::Edge(ex.Current());
352 TopTools_ListOfShape NEL;
353 const TopTools_ListOfShape& LS = nproj.Generated(CE);
354 TopTools_ListIteratorOfListOfShape it(LS);
355 for (;it.More();it.Next())
357 const TopoDS_Shape& PCE = it.Value();
358 TopTools_ListOfShape PLS = anAlgo.Modified(PCE);
361 TopTools_ListIteratorOfListOfShape itp(PLS);
362 for (;itp.More();itp.Next())
363 NEL.Append(itp.Value());
370 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
371 bool HYDROData_DTM::Get2dFaceFrom3dPres(const TopoDS_Compound& cmp, TopoDS_Face& outF,
372 TopTools_SequenceOfShape* Boundr, std::set<int> ind )
374 //ind : set of indices (starts with 0). index == number of boundary (inlet, outlet, etc..)
376 //if Boundr is not null => this method will return sequence of boundary wires (inlet, outlet...)
378 Handle_Geom_Plane refpl = new Geom_Plane(gp_Pnt(0,0,0), gp_Dir(0,0,1));
379 TopTools_DataMapOfShapeListOfShape E2PE;
380 ProjWireOnPlane(cmp, refpl, &E2PE);
381 TopTools_ListOfShape ELL;
383 TopoDS_Iterator it(cmp);
385 for (;it.More(); it.Next())
387 const TopoDS_Wire& W = TopoDS::Wire(it.Value());
388 if (W.Orientation() != TopAbs_INTERNAL)
391 TopExp_Explorer ex(W, TopAbs_EDGE);
392 TopTools_ListOfShape LEpW;
393 TopTools_ListOfShape LEpWDummy;
394 for (;ex.More();ex.Next())
396 const TopoDS_Edge& CE = TopoDS::Edge(ex.Current());
397 TopTools_ListOfShape LS = E2PE.Find(CE);
401 if (ind.count(i) != 0)
405 const TopoDS_Wire& WMW = WM.Wire();
406 //assume that wire is a straight line,
407 //take first and last vertex and make simple edge (RE)
408 TopoDS_Vertex VF, VL;
409 TopExp::Vertices(WMW, VF, VL);
410 TopoDS_Edge RE = BRepLib_MakeEdge(VF, VL).Edge();
413 LEpWDummy = LEpW; //LEpW will be nullified after appending to ELL
418 LEpWDummy.Append(RE);
421 //TODO: in the new version of OCCT, USD can process separate wires
422 //ShapeUpgrade_UnifySameDomain USD(WMW, 1U, 0U, 1U); //concat bsplines
424 //const TopoDS_Shape& RSU = USD.Shape();
425 //TopExp_Explorer exp(RSU, TopAbs_EDGE);
426 //TopTools_ListOfShape DummyL;
427 //for (;it.More();it.Next())
428 // DummyL.Append(it.Value());
429 //if (DummyL.Extent() == 1)
430 // ELL.Append(DummyL.First()); //if one edge => accept this result
432 // ELL.Append(LEpW); //else put 'as is'
442 //make inlet, outlet, left/tight banks [wires]
443 //shouldn't change topology of the edges
444 BRepLib_MakeWire IWM;
446 Boundr->Append(IWM.Wire());
453 BRepLib_MakeWire WME;
456 const TopoDS_Wire& resW = WME.Wire();
457 BRepBuilderAPI_MakeFace mf(refpl, resW, true); //check inside is true by def
460 ShapeAnalysis_Wire WA(resW, outF, Precision::Confusion());
461 bool res = WA.CheckSelfIntersection(); //TODO check that this return false if OK
464 ///!!! the internal wires cant be added with 'internal' ori.
465 // it's possible to do with brep builder yet the result will not be correct!
466 // more proper way is to use BOP operation here.
467 /*for (int i = 1; i <= IntW.Extent(); i++)
470 const TopoDS_Wire& W = TopoDS::Wire(IntW(i));
471 ProjWireOnPlane(W, refpl, outIW);
476 void HYDROData_DTM::CreateProfiles(const std::vector<Handle_HYDROData_Profile>& theProfiles,
478 double theSpatialStep,
479 AltitudePoints& theOutLeft,
480 AltitudePoints& theOutRight,
481 AltitudePoints& theOutPoints,
482 std::vector<AltitudePoints>& theOutMainProfiles,
483 TopoDS_Shape& Out3dPres,
484 TopoDS_Shape& Out2dPres,
485 TopoDS_Shape& OutLeftB,
486 TopoDS_Shape& OutRightB,
487 TopoDS_Shape& OutInlet,
488 TopoDS_Shape& OutOutlet,
491 std::set<int>& InvInd,
492 bool& WireIntersections)
494 if (theProfiles.empty())
496 theOutPoints = Interpolate( theProfiles, theDDZ, theSpatialStep, theOutLeft, theOutRight, theOutMainProfiles, InvInd );
497 //note that if Create3dPres is false => Create2dPres flag is meaningless!
498 if( theOutPoints.empty() )
503 TopTools_IndexedMapOfOrientedShape ll = Create3DShape( theOutLeft, theOutRight, theOutMainProfiles);
509 BB.MakeCompound(cmp);
510 for (int i = 1; i <= ll.Extent(); i++)
515 //same order as in HYDROData_DTM::Update()
517 OutRightB = ll(ll.Extent());
519 OutOutlet = ll(ll.Extent() - 1);
524 WireIntersections = Get2dFaceFrom3dPres(cmp, outF); //__TODO
533 void HYDROData_DTM::GetProperties( const Handle_HYDROData_Profile& theProfile,
534 gp_Pnt& theLowestPoint, gp_Vec2d& theDir,
535 double& theZMin, double& theZMax )
537 theLowestPoint = theProfile->GetBottomPoint();
540 theProfile->GetLeftPoint( aLeft, true, true );
541 theProfile->GetRightPoint( aRight, true, true );
542 double x = aRight.X()-aLeft.X();
543 double y = aRight.Y()-aLeft.Y();
544 theDir = gp_Vec2d( x, y );
546 HYDROData_Profile::ProfilePoints points = theProfile->GetProfilePoints();
547 int lo = points.Lower();
548 int up = points.Upper();
549 theZMin = std::numeric_limits<double>::max();
551 for( int i=lo; i<=up; i++ )
553 double z = points.Value( i ).Z();
561 inline gp_Pnt2d To2D( const gp_Pnt& thePnt, const gp_Trsf& theTr,
562 double& theUMin, double& theUMax )
564 gp_Pnt p = thePnt.Transformed( theTr );
571 return gp_Pnt2d( u, z );
574 Handle(TColgp_HArray1OfPnt2d) To2D( const TColgp_Array1OfPnt& thePoints,
575 const gp_Trsf& theTr,
576 double& theUMin, double& theUMax )
578 int low = thePoints.Lower(), up = thePoints.Upper();
579 Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( low, up );
580 for( int i=low; i<=up; i++ )
581 points->SetValue( i, To2D( thePoints.Value( i ), theTr, theUMin, theUMax ) );
585 Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve,
586 Standard_Real theFirst, Standard_Real theLast,
587 const gp_Trsf& theTr,
588 double& theUMin, double& theUMax )
590 if( theCurve->IsKind( STANDARD_TYPE( Geom_Line ) ) )
592 gp_Pnt aFirstPnt, aLastPnt;
593 theCurve->D0( theFirst, aFirstPnt );
594 theCurve->D0( theLast, aLastPnt );
597 aFirst2d = To2D( aFirstPnt, theTr, theUMin, theUMax ),
598 aLast2d = To2D( aLastPnt, theTr, theUMin, theUMax );
600 gp_Vec2d dir( aFirst2d, aLast2d );
601 Handle_Geom2d_Line aLine2d = new Geom2d_Line( aFirst2d, gp_Dir2d( dir.X(), dir.Y() ) );
602 return new Geom2d_TrimmedCurve( aLine2d, 0, aLast2d.Distance( aFirst2d ) );
605 if( theCurve->IsKind( STANDARD_TYPE( Geom_BSplineCurve ) ) )
607 Handle(Geom_BSplineCurve) aSpline = Handle(Geom_BSplineCurve)::DownCast( theCurve );
609 Handle(TColgp_HArray1OfPnt2d) poles = To2D( aSpline->Poles(), theTr, theUMin, theUMax );
610 const TColStd_Array1OfReal& knots = aSpline->Knots();
611 const TColStd_Array1OfInteger& multiplicities = aSpline->Multiplicities();
612 int aDegree = aSpline->Degree();
614 return new Geom2d_BSplineCurve( poles->Array1(), knots, multiplicities, aDegree );
617 return Handle(Geom2d_Curve)();
620 #include <GCE2d_MakeSegment.hxx>
621 #include <Geom2dAPI_InterCurveCurve.hxx>
622 bool IsCooriented( const Handle_HYDROData_Profile& theProfile1,
623 const Handle_HYDROData_Profile& theProfile2 )
625 if( theProfile1==theProfile2 )
628 gp_XY lp1, rp1, lp2, rp2;
629 theProfile1->GetLeftPoint(lp1);
630 theProfile1->GetRightPoint(rp1);
631 theProfile2->GetLeftPoint(lp2);
632 theProfile2->GetRightPoint(rp2);
634 GCE2d_MakeSegment s1(lp1, lp2);
635 GCE2d_MakeSegment s2(rp1, rp2);
637 Geom2dAPI_InterCurveCurve inter;
639 if (inter.NbPoints() == 0)
645 Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis(
646 const std::vector<Handle_HYDROData_Profile>& theProfiles,
647 std::vector<double>& theDistances )
649 size_t n = theProfiles.size();
651 return Handle_Geom2d_BSplineCurve();
653 Handle_Geom2d_BSplineCurve aResult;
655 Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( 1, (int)n );
656 TColgp_Array1OfVec2d tangents( 1, (int)n );
657 Handle(TColStd_HArray1OfBoolean) flags = new TColStd_HArray1OfBoolean( 1, (int)n );
659 for( size_t i = 1; i <= n; i++ )
661 Handle_HYDROData_Profile aProfile = theProfiles[i-1];
662 Handle_HYDROData_Profile aPrevProfile = i==1 ? theProfiles[i-1] : theProfiles[i-2];
663 Handle_HYDROData_Profile aNextProfile = i==n ? theProfiles[i-1] : theProfiles[i];
665 if( !IsCooriented( aProfile, aNextProfile ) )
668 aProfile->GetLeftPoint( lp, true );
669 aProfile->GetRightPoint( rp, true );
670 aProfile->SetLeftPoint( rp, true );
671 aProfile->SetRightPoint( lp, true );
679 gp_XYZ curP = aProfile->GetBottomPoint();
680 gp_XY curP2d = gp_XY(curP.X(), curP.Y());
684 nextP = aPrevProfile->GetBottomPoint(true);
686 nextP = aNextProfile->GetBottomPoint(true);
688 gp_XY nextP2d = gp_XY(nextP.X(), nextP.Y());
691 GetProperties( aProfile, aLowest, aPrTangent, zmin, zmax );
692 aNormal.SetCoord( -aPrTangent.Y(), aPrTangent.X() );
694 gp_Vec2d aDirToNextProfile(nextP2d.X() - curP2d.X(), nextP2d.Y() - curP2d.Y() );
696 aDirToNextProfile.Reverse();
697 if (aNormal.Dot(aDirToNextProfile) < 0)
702 points->SetValue( (int)i, gp_Pnt2d( aLowest.X(), aLowest.Y() ) );
703 tangents.SetValue( (int)i, aNormal );
704 flags->SetValue( (int)i, Standard_True );
707 Geom2dAPI_Interpolate anInterpolator( points, Standard_False, Standard_False );
708 anInterpolator.Load( tangents, flags );
709 anInterpolator.Perform();
710 if( anInterpolator.IsDone() )
712 aResult = anInterpolator.Curve();
714 //fill the distances vector
715 Geom2dAdaptor_Curve anAdaptor( aResult );
717 theDistances.clear();
718 theDistances.reserve( n );
719 Standard_Real aParamFirst = anAdaptor.FirstParameter(), aParamLast = anAdaptor.LastParameter();
720 for( size_t i = 1; i <= n; i++ )
722 gp_Pnt2d aPnt = points->Value( (Standard_Integer)i );
723 Geom2dAPI_ProjectPointOnCurve aProject( aPnt, aResult );
724 Standard_Real aParam = aProject.LowerDistanceParameter();
725 double aDistance = GCPnts_AbscissaPoint::Length( anAdaptor, aParamFirst, aParam );
726 theDistances.push_back( aDistance );
732 std::vector<Handle_Geom2d_Curve> HYDROData_DTM::ProfileToParametric(
733 const Handle_HYDROData_Profile& theProfile,
734 double& theUMin, double& theUMax, gp_Vec2d& theDir )
736 std::vector<Handle_Geom2d_Curve> curves;
738 // Transformation of the coordinate systems
741 GetProperties( theProfile, aLowest, theDir, zmin, zmax );
743 gp_Ax3 aStd3d( gp_Pnt( 0, 0, 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( 1, 0, 0 ) );
744 gp_Ax3 aLocal( gp_Pnt( aLowest.X(), aLowest.Y(), 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( theDir.X(), theDir.Y(), 0 ) );
747 aTransf.SetTransformation( aStd3d, aLocal );
749 // Iteration via edges
750 TopoDS_Wire aWire = TopoDS::Wire( theProfile->GetShape3D() );
751 TopExp_Explorer anExp( aWire, TopAbs_EDGE );
752 for( ; anExp.More(); anExp.Next() )
754 // Extract an edge from wire
755 TopoDS_Edge anEdge = TopoDS::Edge( anExp.Current() );
757 // Extract a curve corresponding to the edge
758 TopLoc_Location aLoc;
759 Standard_Real aFirst, aLast;
760 Handle(Geom_Curve) aCurve = BRep_Tool::Curve( anEdge, aLoc, aFirst, aLast );
762 // Convert the curve to 2d CS
763 Handle(Geom2d_Curve) aCurve2d = CurveTo2D( aCurve, aFirst, aLast, aTransf, theUMin, theUMax );
764 if( !aCurve2d.IsNull() )
765 curves.push_back( aCurve2d );
771 bool CalcMidWidth( const std::set<double>& intersections, double& theMid, double& theWid )
773 double umin = std::numeric_limits<double>::max(),
776 size_t n = intersections.size();
780 std::set<double>::const_iterator it = intersections.begin(), last = intersections.end();
781 for( ; it!=last; it++ )
789 theMid = ( umin+umax )/2;
794 void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& theProfile,
795 double theXCurv, double theMinZ, double theMaxZ, double theDDZ,
796 CurveUZ& theMidPointCurve,
797 CurveUZ& theWidthCurve,
798 int& intersection_nb,
801 double aDblMax = std::numeric_limits<double>::max(),
806 gp_Vec2d aProfileDir;
807 std::vector<Handle_Geom2d_Curve> curves = ProfileToParametric( theProfile, aUMin, aUMax, aProfileDir );
808 size_t n = curves.size();
813 // we add the "virtual" vertical lines to simulate the intersection with profile
814 gp_Pnt2d aFirst, aLast;
815 curves[0]->D0( curves[0]->FirstParameter(), aFirst );
816 curves[n-1]->D0( curves[n-1]->LastParameter(), aLast );
817 Handle(Geom2d_Line) aV1 = new Geom2d_Line( aFirst, gp_Dir2d( 0, 1 ) );
818 Handle(Geom2d_TrimmedCurve) aT1 = new Geom2d_TrimmedCurve( aV1, 0.0, aVMax );
820 Handle(Geom2d_Line) aV2 = new Geom2d_Line( aLast, gp_Dir2d( 0, 1 ) );
821 Handle(Geom2d_TrimmedCurve) aT2 = new Geom2d_TrimmedCurve( aV2, 0.0, aVMax );
823 curves.push_back( aT1 );
824 curves.push_back( aT2 );
826 int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 );
827 theMidPointCurve = CurveUZ( theXCurv, aProfileDir, theMinZ );
828 theMidPointCurve.reserve( psize );
829 theWidthCurve = CurveUZ( theXCurv, aProfileDir, theMinZ );
830 theWidthCurve.reserve( psize );
833 // for each discrete value of z we search intersection with profile
834 for( double z1 = theMinZ; z1 <= theMaxZ; z1 += theDDZ )
836 Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z1 ), gp_Dir2d( 1, 0 ) );
837 std::set<double> intersections;
838 for( size_t i = 0; i < n; i++ )
840 Handle_Geom2d_Curve aCurve = curves[i];
841 Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
842 for( int k=1, m=anIntersect.NbPoints(); k<=m; k++ )
843 intersections.insert( anIntersect.Point( k ).X() );
846 intersection_nb = intersections.size();
847 if( intersection_nb >= 1 )
850 if( !CalcMidWidth( intersections, u_mid, u_wid ) )
853 double z = z1 - theMinZ;
857 theMidPointCurve.push_back( p_mid );
862 theWidthCurve.push_back( p_wid );
867 void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCurveB,
868 int theNbSteps, std::vector<CurveUZ>& theInterpolation,
871 theInterpolation.clear();
872 int d = isAddSecond ? 2 : 1;
873 theInterpolation.reserve( theNbSteps+d );
874 double dt = 1.0 / double( theNbSteps + 1 );
876 theInterpolation.push_back( theCurveA );
877 for( int i=0; i<theNbSteps; i++, t+=dt )
879 CurveUZ anInterp = theCurveA*(1-t) + theCurveB*t;
880 theInterpolation.push_back( anInterp );
883 theInterpolation.push_back( theCurveB );
885 #include <BRepLib_MakeEdge2d.hxx>
886 void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
887 const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
888 AltitudePoints& thePoints )
890 Geom2dAdaptor_Curve anAdaptor( theHydraulicAxis );
891 TopoDS_Edge E2d = BRepLib_MakeEdge2d(theHydraulicAxis).Edge();
892 GCPnts_AbscissaPoint ap( anAdaptor, theMidCurve.Xcurv(), anAdaptor.FirstParameter() );
893 double aParam = ap.Parameter();
896 anAdaptor.D0( aParam, point );
897 gp_Vec2d profile_dir = theMidCurve.ProfileDir();
898 //gp_Dir tangent_n( -profile_dir.Y(), profile_dir.X(), dz );
899 profile_dir.Normalize();
901 size_t n = theMidCurve.size();
902 std::map<double, AltitudePoint> sorted_points;
903 for( size_t i=0; i<n; i++ )
905 double param1 = theMidCurve[i].U - theWidthCurve[i].U / 2;
906 double param2 = theMidCurve[i].U + theWidthCurve[i].U / 2;
908 gp_Pnt2d p1 = point.Translated( param1 * profile_dir);
909 gp_Pnt2d p2 = point.Translated( param2 * profile_dir);
911 double z = theMidCurve[i].Z + theMidCurve.DeltaZ();
913 AltitudePoint p3d_1( p1.X(), p1.Y(), z ), p3d_2( p2.X(), p2.Y(), z );
915 sorted_points[param1] = p3d_1;
916 sorted_points[param2] = p3d_2;
919 thePoints.reserve( sorted_points.size() );
920 const double EPS = 1E-12;
921 std::map<double, AltitudePoint>::const_iterator it = sorted_points.begin(), last = sorted_points.end();
922 for( ; it!=last; it++ )
923 if( thePoints.empty() || thePoints.back().SquareDistance( it->second ) > EPS )
924 thePoints.push_back( it->second );
927 inline double max( double a, double b )
935 inline double min( double a, double b )
943 #include <BRepLib_MakeWire.hxx>
945 std::vector<HYDROData_Bathymetry::AltitudePoints> HYDROData_DTM::Interpolate
946 ( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
947 const Handle_HYDROData_Profile& theProfileA,
949 const Handle_HYDROData_Profile& theProfileB,
951 double theDDZ, int theNbSteps, bool isAddSecond,
952 int& inter_nb_1, int& inter_nb_2)
954 double zminA, zmaxA, zminB, zmaxB;
955 gp_Pnt lowestA, lowestB;
958 GetProperties( theProfileA, lowestA, dirA, zminA, zmaxA );
959 GetProperties( theProfileB, lowestB, dirB, zminB, zmaxB );
962 double hmax = max( zmaxA-zminA, zmaxB-zminB );
964 //double dz = zminB - zminA;
965 //double zmin = min( zminA, zminB );
966 //double zmax = max( zmaxA, zmaxB );
968 CurveUZ midA(0, gp_Vec2d(), 0), midB(0, gp_Vec2d(), 0);
969 CurveUZ widA(0, gp_Vec2d(), 0), widB(0, gp_Vec2d(), 0);
971 ProfileDiscretization( theProfileA, theXCurvA, zminA, zminA+hmax, theDDZ, midA, widA, inter_nb_1 );
972 ProfileDiscretization( theProfileB, theXCurvB, zminB, zminB+hmax, theDDZ, midB, widB, inter_nb_2 );
974 std::vector<CurveUZ> mid, wid;
975 Interpolate( midA, midB, theNbSteps, mid, isAddSecond );
976 Interpolate( widA, widB, theNbSteps, wid, isAddSecond );
978 size_t p = mid.size();
979 size_t q = p>0 ? 2*mid[0].size() : 1;
980 std::vector<AltitudePoints> points;
983 for( size_t i=0; i<p; i++ )
985 points[i].reserve( q );
986 CurveTo3D( theHydraulicAxis, mid[i], wid[i], points[i] );
992 HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
993 ( const std::vector<Handle_HYDROData_Profile>& theProfiles,
994 double theDDZ, double theSpatialStep,
995 AltitudePoints& theLeft,
996 AltitudePoints& theRight,
997 std::vector<AltitudePoints>& theMainProfiles,
998 std::set<int>& invalInd)
1000 AltitudePoints points;
1001 size_t n = theProfiles.size();
1005 std::vector<double> distances;
1006 Handle_Geom2d_BSplineCurve aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
1007 if( aHydraulicAxis.IsNull() )
1010 theMainProfiles.reserve( n );
1012 for( size_t i=0, n1=n-1; i<n1; i++ )
1014 double aDistance = distances[i+1]-distances[i];
1015 int aNbSteps = int(aDistance/theSpatialStep);
1016 bool isAddSecond = i==n1-1;
1018 // 1. Calculate interpolated profiles
1019 int inter_nb_1, inter_nb_2;
1020 std::vector<AltitudePoints> local_points = Interpolate( aHydraulicAxis, theProfiles[i], distances[i],
1021 theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond, inter_nb_1, inter_nb_2 );
1022 int lps = local_points.size();
1028 invalInd.insert(i+1);
1030 // 2. Put all points into the global container
1031 for( size_t j=0; j<lps; j++ )
1033 const AltitudePoints& lp = local_points[j];
1035 points.reserve( lp.size() * n );
1036 for( size_t k=0, ks=lp.size(); k<ks; k++ )
1037 points.push_back( lp[k] );
1040 // 3. Get left/right banks' points
1043 theLeft.reserve( lps * n );
1044 theRight.reserve( lps * n );
1046 for( size_t j=0; j<lps; j++ )
1048 const AltitudePoints& lp = local_points[j];
1049 theLeft.push_back( lp[0] );
1050 theRight.push_back( lp[lp.size()-1] );
1053 // 4. Get main profiles points
1054 theMainProfiles.push_back( local_points[0] );
1056 theMainProfiles.push_back( local_points[lps-1] );
1061 int HYDROData_DTM::EstimateNbPoints( const std::vector<Handle_HYDROData_Profile>& theProfiles,
1062 double theDDZ, double theSpatialStep )
1064 size_t n = theProfiles.size();
1067 if( theDDZ<1E-6 || theSpatialStep<1E-6 )
1070 std::vector<double> distances;
1071 Handle_Geom2d_BSplineCurve aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
1072 if( aHydraulicAxis.IsNull() )
1075 double aCompleteDistance = distances[n-1];
1076 int aNbSteps = int( aCompleteDistance / theSpatialStep ) + 1;
1079 double aZMin, aZMax;
1080 GetProperties( theProfiles[0], aLowest, aDir, aZMin, aZMax );
1081 int aNbZSteps = (aZMax-aZMin)/theDDZ;
1083 if( aNbSteps > ( 1<<16 ) || aNbZSteps > ( 1<<16 ) )
1086 return aNbSteps * aNbZSteps;