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 //@MZN zzz IMPLEMENT_STANDARD_HANDLE( HYDROData_DTM, HYDROData_Bathymetry )
57 //@MZN zzz IMPLEMENT_STANDARD_RTTIEXT( HYDROData_DTM, HYDROData_Bathymetry )
59 HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv, const gp_Vec2d& theProfileDir, double theDeltaZ )
60 : myXcurv( theXCurv ), myProfileDir( theProfileDir ), myDeltaZ( theDeltaZ )
64 HYDROData_DTM::CurveUZ::~CurveUZ()
68 double HYDROData_DTM::CurveUZ::Xcurv() const
73 gp_Vec2d HYDROData_DTM::CurveUZ::ProfileDir() const
78 double HYDROData_DTM::CurveUZ::DeltaZ() const
83 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator + ( const CurveUZ& c ) const
85 HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv(), ProfileDir() + c.ProfileDir(), DeltaZ() + c.DeltaZ() );
86 size_t n = size(), n1 = c.size();
89 std::cout << "Warning: different number of points in curves: " << n << ", " << n1 << std::endl;
92 for( int i=0; i<n; i++ )
95 p.U = operator[]( i ).U + c[i].U;
96 p.Z = operator[]( i ).Z;
102 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator * ( double d ) const
104 HYDROData_DTM::CurveUZ res( Xcurv()*d, ProfileDir()*d, DeltaZ()*d );
107 for( int i=0; i<n; i++ )
110 p.U = operator[]( i ).U * d;
111 p.Z = operator[]( i ).Z;
120 HYDROData_DTM::HYDROData_DTM()
124 HYDROData_DTM::~HYDROData_DTM()
128 HYDROData_SequenceOfObjects HYDROData_DTM::GetProfiles() const
130 return GetReferenceObjects( DataTag_Profiles );
133 void HYDROData_DTM::SetProfiles( const HYDROData_SequenceOfObjects& theProfiles )
135 SetReferenceObjects( theProfiles, DataTag_Profiles );
139 double HYDROData_DTM::GetDDZ() const
141 return GetDouble( DataTag_DDZ );
144 void HYDROData_DTM::SetDDZ( double theDDZ )
146 SetDouble( DataTag_DDZ, theDDZ );
150 double HYDROData_DTM::GetSpatialStep() const
152 return GetDouble( DataTag_SpatialStep );
155 void HYDROData_DTM::SetSpatialStep( double theSpatialStep )
157 SetDouble( DataTag_SpatialStep, theSpatialStep );
161 void HYDROData_DTM::PointToWire(const AltitudePoints& pnts, TopoDS_Wire& W )
163 /*BRepLib_MakeWire WM;
166 for (int i = 0; i < pnts.size() - 1; i++)
168 gp_Pnt p1(pnts[i].X, pnts[i].Y, pnts[i].Z);
169 gp_Pnt p2(pnts[i+1].X, pnts[i+1].Y, pnts[i+1].Z);
170 WM.Add(BRepLib_MakeEdge(p1, p2).Edge());
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 TopTools_IndexedMapOfOrientedShape HYDROData_DTM::Create3DShape(const AltitudePoints& left,
182 const AltitudePoints& right,
183 const std::vector<AltitudePoints>& main_profiles)
185 TopTools_IndexedMapOfOrientedShape ll;
186 TopoDS_Wire LWire, RWire;
187 PointToWire(left, LWire);
188 PointToWire(right, RWire);
190 ll.Add(LWire.Oriented(TopAbs_FORWARD));
192 for (int k = 0; k < main_profiles.size(); k++)
195 PointToWire(main_profiles[k], W);
196 TopAbs_Orientation Ori = TopAbs_INTERNAL;
197 if (k == 0 || k == main_profiles.size() - 1)
198 Ori = TopAbs_FORWARD;
199 ll.Add(W.Oriented(Ori));
203 ll.Add(RWire.Oriented(TopAbs_FORWARD));
204 //yes, add subshapes in this order (left + profiles + right)
205 //otherwise the projected wire will be non-manifold
211 void HYDROData_DTM::Update()
213 AltitudePoints points;
214 TopoDS_Shape Out3dPres;
215 TopoDS_Shape Out2dPres;
216 TopoDS_Shape OutLeftB;
217 TopoDS_Shape OutRightB;
218 TopoDS_Shape OutInlet;
219 TopoDS_Shape OutOutlet;
221 HYDROData_SequenceOfObjects objs = GetProfiles();
222 double ddz = GetDDZ();
223 double step = GetSpatialStep();
224 std::set<int> InvInd;
225 bool WireIntersections; //__TODO
226 CreateProfilesFromDTM( objs, ddz, step, points, Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, true, true, InvInd, -1, WireIntersections );
227 SetAltitudePoints( points );
229 SetShape( DataTag_LeftBankShape, OutLeftB);
230 SetShape( DataTag_RightBankShape, OutRightB);
231 SetShape( DataTag_InletShape, OutInlet);
232 SetShape( DataTag_OutletShape, OutOutlet );
233 SetShape( DataTag_3DShape, Out3dPres );
234 SetShape( DataTag_2DShape, Out2dPres );
236 HYDROData_Bathymetry::Update();
239 void HYDROData_DTM::GetPresentationShapes( TopoDS_Shape& Out3dPres,
240 TopoDS_Shape& Out2dPres,
241 TopoDS_Shape& OutLeftB,
242 TopoDS_Shape& OutRightB,
243 TopoDS_Shape& OutInlet,
244 TopoDS_Shape& OutOutlet )
247 OutLeftB = GetShape( DataTag_LeftBankShape);
248 OutRightB = GetShape( DataTag_RightBankShape);
249 OutInlet = GetShape( DataTag_InletShape);
250 OutOutlet = GetShape( DataTag_OutletShape );
251 Out3dPres = GetShape( DataTag_3DShape );
252 Out2dPres = GetShape( DataTag_2DShape );
254 void HYDROData_DTM::CreateProfilesFromDTM (const HYDROData_SequenceOfObjects& InpProfiles,
257 AltitudePoints& points,
258 TopoDS_Shape& Out3dPres,
259 TopoDS_Shape& Out2dPres,
260 TopoDS_Shape& OutLeftB,
261 TopoDS_Shape& OutRightB,
262 TopoDS_Shape& OutInlet,
263 TopoDS_Shape& OutOutlet,
266 std::set<int>& InvInd,
268 bool& WireIntersections)
270 int aLower = InpProfiles.Lower(), anUpper = InpProfiles.Upper();
271 size_t n = anUpper - aLower + 1;
273 std::vector<Handle(HYDROData_Profile)> profiles;
274 profiles.reserve( n );
275 for( int i=aLower; i<=anUpper; i++ )
277 Handle(HYDROData_Profile) aProfile = Handle(HYDROData_Profile)::DownCast( InpProfiles.Value( i ) );
278 if( !aProfile.IsNull() )
279 profiles.push_back( aProfile );
281 const double EPS = 1E-3;
283 AltitudePoints right;
284 std::vector<AltitudePoints> main_profiles;
286 if( thePntsLimit > 0 )
288 int aNbPoints = EstimateNbPoints( profiles, ddz, step );
289 if( aNbPoints < 0 || aNbPoints > thePntsLimit )
293 if( ddz>EPS && step>EPS )
294 CreateProfiles(profiles, ddz, step, left, right, points, main_profiles,
295 Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, Create3dPres, Create2dPres, InvInd, WireIntersections );
298 void HYDROData_DTM::ProjWireOnPlane(const TopoDS_Shape& inpWire, const Handle(Geom_Plane)& RefPlane,
299 TopTools_DataMapOfShapeListOfShape* E2PE)
303 //project shape (edges) on planar face
305 BB.MakeFace(F, RefPlane, Precision::Confusion());
306 BRepAlgo_NormalProjection nproj(F);
308 nproj.SetDefaultParams();
313 TopoDS_Shape projRes = nproj.Projection();
315 // unite all vertexes/edges from projected result
316 BOPAlgo_Builder anAlgo;
317 TopExp_Explorer exp(projRes, TopAbs_EDGE);
318 for (;exp.More(); exp.Next())
320 const TopoDS_Edge& E = TopoDS::Edge(exp.Current());
321 if (E.Orientation() != TopAbs_INTERNAL)
322 anAlgo.AddArgument(E);
325 int stat = anAlgo.ErrorStatus();
326 TopoDS_Shape projResConn = anAlgo.Shape();
328 // make wire => vertexes and edges should be untouched after this operation!
329 exp.Init(projResConn, TopAbs_EDGE);
330 //TopTools_ListOfShape llE;
332 //BRepLib_MakeWire WM;
334 //for (;exp.More();exp.Next())
335 // llE.Append(exp.Current());
338 //outWire = WM.Wire();
340 //outWire.Orientation(inpWire.Orientation()); //take from the original wire
342 //history mode: edge to projected edges
345 TopExp_Explorer ex(inpWire, TopAbs_EDGE);
346 for (;ex.More();ex.Next())
348 const TopoDS_Edge& CE = TopoDS::Edge(ex.Current());
349 TopTools_ListOfShape NEL;
350 const TopTools_ListOfShape& LS = nproj.Generated(CE);
351 TopTools_ListIteratorOfListOfShape it(LS);
352 for (;it.More();it.Next())
354 const TopoDS_Shape& PCE = it.Value();
355 TopTools_ListOfShape PLS = anAlgo.Modified(PCE);
358 TopTools_ListIteratorOfListOfShape itp(PLS);
359 for (;itp.More();itp.Next())
360 NEL.Append(itp.Value());
367 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
368 bool HYDROData_DTM::Get2dFaceFrom3dPres(const TopoDS_Compound& cmp, TopoDS_Face& outF,
369 TopTools_SequenceOfShape* Boundr, std::set<int> ind )
371 //ind : set of indices (starts with 0). index == number of boundary (inlet, outlet, etc..)
373 //if Boundr is not null => this method will return sequence of boundary wires (inlet, outlet...)
375 Handle(Geom_Plane) refpl = new Geom_Plane(gp_Pnt(0,0,0), gp_Dir(0,0,1));
376 TopTools_DataMapOfShapeListOfShape E2PE;
377 ProjWireOnPlane(cmp, refpl, &E2PE);
378 TopTools_ListOfShape ELL;
380 TopoDS_Iterator it(cmp);
382 for (;it.More(); it.Next())
384 const TopoDS_Wire& W = TopoDS::Wire(it.Value());
385 if (W.Orientation() != TopAbs_INTERNAL)
388 TopExp_Explorer ex(W, TopAbs_EDGE);
389 TopTools_ListOfShape LEpW;
390 TopTools_ListOfShape LEpWDummy;
391 for (;ex.More();ex.Next())
393 const TopoDS_Edge& CE = TopoDS::Edge(ex.Current());
394 TopTools_ListOfShape LS = E2PE.Find(CE);
398 if (ind.count(i) != 0)
402 const TopoDS_Wire& WMW = WM.Wire();
403 //assume that wire is a straight line,
404 //take first and last vertex and make simple edge (RE)
405 TopoDS_Vertex VF, VL;
406 TopExp::Vertices(WMW, VF, VL);
407 TopoDS_Edge RE = BRepLib_MakeEdge(VF, VL).Edge();
410 LEpWDummy = LEpW; //LEpW will be nullified after appending to ELL
415 LEpWDummy.Append(RE);
418 //TODO: in the new version of OCCT, USD can process separate wires
419 //ShapeUpgrade_UnifySameDomain USD(WMW, 1U, 0U, 1U); //concat bsplines
421 //const TopoDS_Shape& RSU = USD.Shape();
422 //TopExp_Explorer exp(RSU, TopAbs_EDGE);
423 //TopTools_ListOfShape DummyL;
424 //for (;it.More();it.Next())
425 // DummyL.Append(it.Value());
426 //if (DummyL.Extent() == 1)
427 // ELL.Append(DummyL.First()); //if one edge => accept this result
429 // ELL.Append(LEpW); //else put 'as is'
439 //make inlet, outlet, left/tight banks [wires]
440 //shouldn't change topology of the edges
441 BRepLib_MakeWire IWM;
443 Boundr->Append(IWM.Wire());
450 BRepLib_MakeWire WME;
453 const TopoDS_Wire& resW = WME.Wire();
454 BRepBuilderAPI_MakeFace mf(refpl, resW, true); //check inside is true by def
457 ShapeAnalysis_Wire WA(resW, outF, Precision::Confusion());
458 bool res = WA.CheckSelfIntersection(); //TODO check that this return false if OK
461 ///!!! the internal wires cant be added with 'internal' ori.
462 // it's possible to do with brep builder yet the result will not be correct!
463 // more proper way is to use BOP operation here.
464 /*for (int i = 1; i <= IntW.Extent(); i++)
467 const TopoDS_Wire& W = TopoDS::Wire(IntW(i));
468 ProjWireOnPlane(W, refpl, outIW);
473 void HYDROData_DTM::CreateProfiles(const std::vector<Handle(HYDROData_Profile)>& theProfiles,
475 double theSpatialStep,
476 AltitudePoints& theOutLeft,
477 AltitudePoints& theOutRight,
478 AltitudePoints& theOutPoints,
479 std::vector<AltitudePoints>& theOutMainProfiles,
480 TopoDS_Shape& Out3dPres,
481 TopoDS_Shape& Out2dPres,
482 TopoDS_Shape& OutLeftB,
483 TopoDS_Shape& OutRightB,
484 TopoDS_Shape& OutInlet,
485 TopoDS_Shape& OutOutlet,
488 std::set<int>& InvInd,
489 bool& WireIntersections)
491 if (theProfiles.empty())
493 theOutPoints = Interpolate( theProfiles, theDDZ, theSpatialStep, theOutLeft, theOutRight, theOutMainProfiles, InvInd );
494 //note that if Create3dPres is false => Create2dPres flag is meaningless!
497 TopTools_IndexedMapOfOrientedShape ll = Create3DShape( theOutLeft, theOutRight, theOutMainProfiles);
503 BB.MakeCompound(cmp);
504 for (int i = 1; i <= ll.Extent(); i++)
509 //same order as in HYDROData_DTM::Update()
511 OutRightB = ll(ll.Extent());
513 OutOutlet = ll(ll.Extent() - 1);
518 WireIntersections = Get2dFaceFrom3dPres(cmp, outF); //__TODO
527 void HYDROData_DTM::GetProperties( const Handle(HYDROData_Profile)& theProfile,
528 gp_Pnt& theLowestPoint, gp_Vec2d& theDir,
530 double& theZMin, double& theZMax )
532 theLowestPoint = theProfile->GetBottomPoint();
535 theProfile->GetLeftPoint( aLeft, true, true );
536 theProfile->GetRightPoint( aRight, true, true );
537 double x = aRight.X()-aLeft.X();
538 double y = aRight.Y()-aLeft.Y();
540 theDir = gp_Vec2d( -y, x );
542 theDir = gp_Vec2d( x, y );
544 HYDROData_Profile::ProfilePoints points = theProfile->GetProfilePoints();
545 int lo = points.Lower();
546 int up = points.Upper();
547 theZMin = std::numeric_limits<double>::max();
549 for( int i=lo; i<=up; i++ )
551 double z = points.Value( i ).Z();
559 inline gp_Pnt2d To2D( const gp_Pnt& thePnt, const gp_Trsf& theTr,
560 double& theUMin, double& theUMax )
562 gp_Pnt p = thePnt.Transformed( theTr );
569 return gp_Pnt2d( u, z );
572 Handle(TColgp_HArray1OfPnt2d) To2D( const TColgp_Array1OfPnt& thePoints,
573 const gp_Trsf& theTr,
574 double& theUMin, double& theUMax )
576 int low = thePoints.Lower(), up = thePoints.Upper();
577 Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( low, up );
578 for( int i=low; i<=up; i++ )
579 points->SetValue( i, To2D( thePoints.Value( i ), theTr, theUMin, theUMax ) );
583 Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve,
584 Standard_Real theFirst, Standard_Real theLast,
585 const gp_Trsf& theTr,
586 double& theUMin, double& theUMax )
588 if( theCurve->IsKind( STANDARD_TYPE( Geom_Line ) ) )
590 gp_Pnt aFirstPnt, aLastPnt;
591 theCurve->D0( theFirst, aFirstPnt );
592 theCurve->D0( theLast, aLastPnt );
595 aFirst2d = To2D( aFirstPnt, theTr, theUMin, theUMax ),
596 aLast2d = To2D( aLastPnt, theTr, theUMin, theUMax );
598 gp_Vec2d dir( aFirst2d, aLast2d );
599 Handle(Geom2d_Line) aLine2d = new Geom2d_Line( aFirst2d, gp_Dir2d( dir.X(), dir.Y() ) );
600 return new Geom2d_TrimmedCurve( aLine2d, 0, aLast2d.Distance( aFirst2d ) );
603 if( theCurve->IsKind( STANDARD_TYPE( Geom_BSplineCurve ) ) )
605 Handle(Geom_BSplineCurve) aSpline = Handle(Geom_BSplineCurve)::DownCast( theCurve );
607 Handle(TColgp_HArray1OfPnt2d) poles = To2D( aSpline->Poles(), theTr, theUMin, theUMax );
608 const TColStd_Array1OfReal& knots = aSpline->Knots();
609 const TColStd_Array1OfInteger& multiplicities = aSpline->Multiplicities();
610 int aDegree = aSpline->Degree();
612 return new Geom2d_BSplineCurve( poles->Array1(), knots, multiplicities, aDegree );
615 return Handle(Geom2d_Curve)();
618 Handle(Geom2d_BSplineCurve) HYDROData_DTM::CreateHydraulicAxis(
619 const std::vector<Handle(HYDROData_Profile)>& theProfiles,
620 std::vector<double>& theDistances )
622 size_t n = theProfiles.size();
623 Handle(Geom2d_BSplineCurve) aResult;
625 Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( 1, (int)n );
626 TColgp_Array1OfVec2d tangents( 1, (int)n );
627 Handle(TColStd_HArray1OfBoolean) flags = new TColStd_HArray1OfBoolean( 1, (int)n );
629 for( size_t i = 1; i <= n; i++ )
631 Handle(HYDROData_Profile) aProfile = theProfiles[i-1];
637 GetProperties( aProfile, aLowest, aTangent, true, zmin, zmax );
638 aTangent.Normalize();
640 points->SetValue( (int)i, gp_Pnt2d( aLowest.X(), aLowest.Y() ) );
641 tangents.SetValue( (int)i, aTangent );
642 flags->SetValue( (int)i, Standard_True );
645 Geom2dAPI_Interpolate anInterpolator( points, Standard_False, Standard_False );
646 anInterpolator.Load( tangents, flags );
647 anInterpolator.Perform();
648 if( anInterpolator.IsDone() )
650 aResult = anInterpolator.Curve();
652 //fill the distances vector
653 Geom2dAdaptor_Curve anAdaptor( aResult );
655 theDistances.clear();
656 theDistances.reserve( n );
657 Standard_Real aParamFirst = anAdaptor.FirstParameter(), aParamLast = anAdaptor.LastParameter();
658 for( size_t i = 1; i <= n; i++ )
660 gp_Pnt2d aPnt = points->Value( (Standard_Integer)i );
661 Geom2dAPI_ProjectPointOnCurve aProject( aPnt, aResult );
662 Standard_Real aParam = aProject.LowerDistanceParameter();
663 double aDistance = GCPnts_AbscissaPoint::Length( anAdaptor, aParamFirst, aParam );
664 theDistances.push_back( aDistance );
670 std::vector<Handle(Geom2d_Curve)> HYDROData_DTM::ProfileToParametric(
671 const Handle(HYDROData_Profile)& theProfile,
672 double& theUMin, double& theUMax, gp_Vec2d& theDir )
674 std::vector<Handle(Geom2d_Curve)> curves;
676 // Transformation of the coordinate systems
679 GetProperties( theProfile, aLowest, theDir, false, zmin, zmax );
681 gp_Ax3 aStd3d( gp_Pnt( 0, 0, 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( 1, 0, 0 ) );
682 gp_Ax3 aLocal( gp_Pnt( aLowest.X(), aLowest.Y(), 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( theDir.X(), theDir.Y(), 0 ) );
685 aTransf.SetTransformation( aStd3d, aLocal );
687 // Iteration via edges
688 TopoDS_Wire aWire = TopoDS::Wire( theProfile->GetShape3D() );
689 TopExp_Explorer anExp( aWire, TopAbs_EDGE );
690 for( ; anExp.More(); anExp.Next() )
692 // Extract an edge from wire
693 TopoDS_Edge anEdge = TopoDS::Edge( anExp.Current() );
695 // Extract a curve corresponding to the edge
696 TopLoc_Location aLoc;
697 Standard_Real aFirst, aLast;
698 Handle(Geom_Curve) aCurve = BRep_Tool::Curve( anEdge, aLoc, aFirst, aLast );
700 // Convert the curve to 2d CS
701 Handle(Geom2d_Curve) aCurve2d = CurveTo2D( aCurve, aFirst, aLast, aTransf, theUMin, theUMax );
702 if( !aCurve2d.IsNull() )
703 curves.push_back( aCurve2d );
709 bool CalcMidWidth( const std::set<double>& intersections, double& theMid, double& theWid )
711 double umin = std::numeric_limits<double>::max(),
714 size_t n = intersections.size();
718 std::set<double>::const_iterator it = intersections.begin(), last = intersections.end();
719 for( ; it!=last; it++ )
727 theMid = ( umin+umax )/2;
732 void HYDROData_DTM::ProfileDiscretization( const Handle(HYDROData_Profile)& theProfile,
733 double theXCurv, double theMinZ, double theMaxZ, double theDDZ,
734 CurveUZ& theMidPointCurve,
735 CurveUZ& theWidthCurve,
736 int& intersection_nb,
739 double aDblMax = std::numeric_limits<double>::max(),
744 gp_Vec2d aProfileDir;
745 std::vector<Handle(Geom2d_Curve)> curves = ProfileToParametric( theProfile, aUMin, aUMax, aProfileDir );
746 size_t n = curves.size();
751 // we add the "virtual" vertical lines to simulate the intersection with profile
752 gp_Pnt2d aFirst, aLast;
753 curves[0]->D0( curves[0]->FirstParameter(), aFirst );
754 curves[n-1]->D0( curves[n-1]->LastParameter(), aLast );
755 Handle(Geom2d_Line) aV1 = new Geom2d_Line( aFirst, gp_Dir2d( 0, 1 ) );
756 Handle(Geom2d_TrimmedCurve) aT1 = new Geom2d_TrimmedCurve( aV1, 0.0, aVMax );
758 Handle(Geom2d_Line) aV2 = new Geom2d_Line( aLast, gp_Dir2d( 0, 1 ) );
759 Handle(Geom2d_TrimmedCurve) aT2 = new Geom2d_TrimmedCurve( aV2, 0.0, aVMax );
761 curves.push_back( aT1 );
762 curves.push_back( aT2 );
764 int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 );
765 theMidPointCurve = CurveUZ( theXCurv, aProfileDir, theMinZ );
766 theMidPointCurve.reserve( psize );
767 theWidthCurve = CurveUZ( theXCurv, aProfileDir, theMinZ );
768 theWidthCurve.reserve( psize );
771 // for each discrete value of z we search intersection with profile
772 for( double z1 = theMinZ; z1 <= theMaxZ; z1 += theDDZ )
774 Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z1 ), gp_Dir2d( 1, 0 ) );
775 std::set<double> intersections;
776 for( size_t i = 0; i < n; i++ )
778 Handle(Geom2d_Curve) aCurve = curves[i];
779 Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
780 for( int k=1, m=anIntersect.NbPoints(); k<=m; k++ )
781 intersections.insert( anIntersect.Point( k ).X() );
784 intersection_nb = intersections.size();
785 if( intersection_nb >= 1 )
788 if( !CalcMidWidth( intersections, u_mid, u_wid ) )
791 double z = z1 - theMinZ;
795 theMidPointCurve.push_back( p_mid );
800 theWidthCurve.push_back( p_wid );
805 void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCurveB,
806 int theNbSteps, std::vector<CurveUZ>& theInterpolation,
809 theInterpolation.clear();
810 int d = isAddSecond ? 2 : 1;
811 theInterpolation.reserve( theNbSteps+d );
812 double dt = 1.0 / double( theNbSteps + 1 );
814 theInterpolation.push_back( theCurveA );
815 for( int i=0; i<theNbSteps; i++, t+=dt )
817 CurveUZ anInterp = theCurveA*(1-t) + theCurveB*t;
818 theInterpolation.push_back( anInterp );
821 theInterpolation.push_back( theCurveB );
823 #include <BRepLib_MakeEdge2d.hxx>
824 void HYDROData_DTM::CurveTo3D( const Handle(Geom2d_BSplineCurve)& theHydraulicAxis,
825 const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
826 AltitudePoints& thePoints )
828 Geom2dAdaptor_Curve anAdaptor( theHydraulicAxis );
829 TopoDS_Edge E2d = BRepLib_MakeEdge2d(theHydraulicAxis).Edge();
830 GCPnts_AbscissaPoint ap( anAdaptor, theMidCurve.Xcurv(), anAdaptor.FirstParameter() );
831 double aParam = ap.Parameter();
834 anAdaptor.D0( aParam, point );
835 gp_Vec2d profile_dir = theMidCurve.ProfileDir();
836 //gp_Dir tangent_n( -profile_dir.Y(), profile_dir.X(), dz );
837 profile_dir.Normalize();
839 size_t n = theMidCurve.size();
840 std::map<double, AltitudePoint> sorted_points;
841 for( size_t i=0; i<n; i++ )
843 double param1 = theMidCurve[i].U - theWidthCurve[i].U / 2;
844 double param2 = theMidCurve[i].U + theWidthCurve[i].U / 2;
846 gp_Pnt2d p1 = point.Translated( param1 * profile_dir);
847 gp_Pnt2d p2 = point.Translated( param2 * profile_dir);
849 double z = theMidCurve[i].Z + theMidCurve.DeltaZ();
851 AltitudePoint p3d_1( p1.X(), p1.Y(), z ), p3d_2( p2.X(), p2.Y(), z );
853 sorted_points[param1] = p3d_1;
854 sorted_points[param2] = p3d_2;
857 thePoints.reserve( sorted_points.size() );
858 const double EPS = 1E-12;
859 std::map<double, AltitudePoint>::const_iterator it = sorted_points.begin(), last = sorted_points.end();
860 for( ; it!=last; it++ )
861 if( thePoints.empty() || thePoints.back().SquareDistance( it->second ) > EPS )
862 thePoints.push_back( it->second );
865 inline double max( double a, double b )
873 inline double min( double a, double b )
881 #include <BRepLib_MakeWire.hxx>
883 std::vector<HYDROData_Bathymetry::AltitudePoints> HYDROData_DTM::Interpolate
884 ( const Handle(Geom2d_BSplineCurve)& theHydraulicAxis,
885 const Handle(HYDROData_Profile)& theProfileA,
887 const Handle(HYDROData_Profile)& theProfileB,
889 double theDDZ, int theNbSteps, bool isAddSecond,
890 int& inter_nb_1, int& inter_nb_2)
892 double zminA, zmaxA, zminB, zmaxB;
893 gp_Pnt lowestA, lowestB;
896 GetProperties( theProfileA, lowestA, dirA, false, zminA, zmaxA );
897 GetProperties( theProfileB, lowestB, dirB, false, zminB, zmaxB );
900 double hmax = max( zmaxA-zminA, zmaxB-zminB );
902 //double dz = zminB - zminA;
903 //double zmin = min( zminA, zminB );
904 //double zmax = max( zmaxA, zmaxB );
906 CurveUZ midA(0, gp_Vec2d(), 0), midB(0, gp_Vec2d(), 0);
907 CurveUZ widA(0, gp_Vec2d(), 0), widB(0, gp_Vec2d(), 0);
909 ProfileDiscretization( theProfileA, theXCurvA, zminA, zminA+hmax, theDDZ, midA, widA, inter_nb_1 );
910 ProfileDiscretization( theProfileB, theXCurvB, zminB, zminB+hmax, theDDZ, midB, widB, inter_nb_2 );
912 std::vector<CurveUZ> mid, wid;
913 Interpolate( midA, midB, theNbSteps, mid, isAddSecond );
914 Interpolate( widA, widB, theNbSteps, wid, isAddSecond );
916 size_t p = mid.size();
917 size_t q = p>0 ? 2*mid[0].size() : 1;
918 std::vector<AltitudePoints> points;
921 for( size_t i=0; i<p; i++ )
923 points[i].reserve( q );
924 CurveTo3D( theHydraulicAxis, mid[i], wid[i], points[i] );
930 HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
931 ( const std::vector<Handle(HYDROData_Profile)>& theProfiles,
932 double theDDZ, double theSpatialStep,
933 AltitudePoints& theLeft,
934 AltitudePoints& theRight,
935 std::vector<AltitudePoints>& theMainProfiles,
936 std::set<int>& invalInd)
938 AltitudePoints points;
939 size_t n = theProfiles.size();
943 std::vector<double> distances;
944 Handle(Geom2d_BSplineCurve) aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
945 if( aHydraulicAxis.IsNull() )
948 theMainProfiles.reserve( n );
950 for( size_t i=0, n1=n-1; i<n1; i++ )
952 double aDistance = distances[i+1]-distances[i];
953 int aNbSteps = int(aDistance/theSpatialStep);
954 bool isAddSecond = i==n1-1;
956 // 1. Calculate interpolated profiles
957 int inter_nb_1, inter_nb_2;
958 std::vector<AltitudePoints> local_points = Interpolate( aHydraulicAxis, theProfiles[i], distances[i],
959 theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond, inter_nb_1, inter_nb_2 );
960 int lps = local_points.size();
966 invalInd.insert(i+1);
968 // 2. Put all points into the global container
969 for( size_t j=0; j<lps; j++ )
971 const AltitudePoints& lp = local_points[j];
973 points.reserve( lp.size() * n );
974 for( size_t k=0, ks=lp.size(); k<ks; k++ )
975 points.push_back( lp[k] );
978 // 3. Get left/right banks' points
981 theLeft.reserve( lps * n );
982 theRight.reserve( lps * n );
984 for( size_t j=0; j<lps; j++ )
986 const AltitudePoints& lp = local_points[j];
987 theLeft.push_back( lp[0] );
988 theRight.push_back( lp[lp.size()-1] );
991 // 4. Get main profiles points
992 theMainProfiles.push_back( local_points[0] );
994 theMainProfiles.push_back( local_points[lps-1] );
999 int HYDROData_DTM::EstimateNbPoints( const std::vector<Handle(HYDROData_Profile)>& theProfiles,
1000 double theDDZ, double theSpatialStep )
1002 size_t n = theProfiles.size();
1005 if( theDDZ<1E-6 || theSpatialStep<1E-6 )
1008 std::vector<double> distances;
1009 Handle(Geom2d_BSplineCurve) aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
1010 if( aHydraulicAxis.IsNull() )
1013 double aCompleteDistance = distances[n-1];
1014 int aNbSteps = int( aCompleteDistance / theSpatialStep ) + 1;
1017 double aZMin, aZMax;
1018 GetProperties( theProfiles[0], aLowest, aDir, true, aZMin, aZMax );
1019 int aNbZSteps = (aZMax-aZMin)/theDDZ;
1021 if( aNbSteps > ( 1<<16 ) || aNbZSteps > ( 1<<16 ) )
1024 return aNbSteps * aNbZSteps;