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>
45 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
47 #include <BRepLib_MakeEdge.hxx>
48 #include <BRepLib_MakeWire.hxx>
49 #include <BRep_Builder.hxx>
50 #include <ShapeAnalysis_Wire.hxx>
51 #include <BRepAlgo_NormalProjection.hxx>
52 #include <ShapeUpgrade_UnifySameDomain.hxx>
53 #include <BRepBuilderAPI_MakePolygon.hxx>
54 #include <BOPAlgo_Builder.hxx>
56 #include <TopTools_IndexedDataMapOfShapeShape.hxx>
57 #include <TopTools_ListIteratorOfListOfShape.hxx>
58 #include <TopTools_SequenceOfShape.hxx>
61 IMPLEMENT_STANDARD_HANDLE( HYDROData_DTM, HYDROData_Bathymetry )
62 IMPLEMENT_STANDARD_RTTIEXT( HYDROData_DTM, HYDROData_Bathymetry )
64 HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv, const gp_Vec2d& theProfileDir, double theDeltaZ )
65 : myXcurv( theXCurv ), myProfileDir( theProfileDir ), myDeltaZ( theDeltaZ )
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 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator + ( const CurveUZ& c ) const
90 HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv(), ProfileDir() + c.ProfileDir(), DeltaZ() + c.DeltaZ() );
91 size_t n = size(), n1 = c.size();
94 std::cout << "Warning: different number of points in curves: " << n << ", " << n1 << std::endl;
97 for( int i=0; i<n; i++ )
100 p.U = operator[]( i ).U + c[i].U;
101 p.Z = operator[]( i ).Z;
107 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator * ( double d ) const
109 HYDROData_DTM::CurveUZ res( Xcurv()*d, ProfileDir()*d, DeltaZ()*d );
112 for( int i=0; i<n; i++ )
115 p.U = operator[]( i ).U * d;
116 p.Z = operator[]( i ).Z;
125 HYDROData_DTM::HYDROData_DTM()
129 HYDROData_DTM::~HYDROData_DTM()
133 HYDROData_SequenceOfObjects HYDROData_DTM::GetProfiles() const
135 return GetReferenceObjects( DataTag_Profiles );
138 void HYDROData_DTM::SetProfiles( const HYDROData_SequenceOfObjects& theProfiles )
140 SetReferenceObjects( theProfiles, DataTag_Profiles );
144 double HYDROData_DTM::GetDDZ() const
146 return GetDouble( DataTag_DDZ );
149 void HYDROData_DTM::SetDDZ( double theDDZ )
151 SetDouble( DataTag_DDZ, theDDZ );
155 double HYDROData_DTM::GetSpatialStep() const
157 return GetDouble( DataTag_SpatialStep );
160 void HYDROData_DTM::SetSpatialStep( double theSpatialStep )
162 SetDouble( DataTag_SpatialStep, theSpatialStep );
166 void HYDROData_DTM::PointsToWire(const AltitudePoints& pnts, TopoDS_Wire& W )
169 BRepBuilderAPI_MakePolygon PM;
170 for (int i = 0; i < pnts.size(); i++)
171 PM.Add(gp_Pnt(pnts[i].X, pnts[i].Y, pnts[i].Z));
176 void HYDROData_DTM::PointsToEdge(const AltitudePoints& pnts, TopoDS_Edge& E )
178 Handle(TColgp_HArray1OfPnt) gpPoints = new TColgp_HArray1OfPnt(1, (int)pnts.size());
180 for (int i = 0; i < pnts.size(); i++)
181 gpPoints->SetValue(i+1, gp_Pnt(pnts[i].X, pnts[i].Y, pnts[i].Z));
183 GeomAPI_Interpolate anInterpolator(gpPoints, Standard_False,1.0e-6);
184 anInterpolator.Perform() ;
185 if (anInterpolator.IsDone())
187 Handle(Geom_Curve) C = anInterpolator.Curve();
188 E = BRepBuilderAPI_MakeEdge(C).Edge();
192 TopTools_IndexedMapOfOrientedShape HYDROData_DTM::Create3DShape(const AltitudePoints& left,
193 const AltitudePoints& right,
194 const std::vector<AltitudePoints>& main_profiles)
196 TopTools_IndexedMapOfOrientedShape ll;
197 //TopoDS_Wire LWire, RWire;
198 //PointsToWire(left, LWire);
199 //PointsToWire(right, RWire);
200 TopoDS_Edge LEdge, REdge;
201 PointsToEdge(left, LEdge);
202 PointsToEdge(right, REdge);
204 ll.Add(LEdge.Oriented(TopAbs_FORWARD));
206 for (int k = 0; k < main_profiles.size(); k++)
209 PointsToWire(main_profiles[k], W);
210 TopAbs_Orientation Ori = TopAbs_INTERNAL;
211 if (k == 0 || k == main_profiles.size() - 1)
212 Ori = TopAbs_FORWARD;
213 ll.Add(W.Oriented(Ori));
217 ll.Add(REdge.Oriented(TopAbs_FORWARD));
218 //yes, add subshapes in this order (left + profiles + right)
219 //otherwise the projected wire will be non-manifold
225 void HYDROData_DTM::Update()
227 AltitudePoints points;
228 TopoDS_Shape Out3dPres;
229 TopoDS_Shape Out2dPres;
230 TopoDS_Shape OutLeftB;
231 TopoDS_Shape OutRightB;
232 TopoDS_Shape OutInlet;
233 TopoDS_Shape OutOutlet;
235 HYDROData_SequenceOfObjects objs = GetProfiles();
236 double ddz = GetDDZ();
237 double step = GetSpatialStep();
238 std::set<int> InvInd;
239 bool WireIntersections; //__TODO
240 CreateProfilesFromDTM( objs, ddz, step, points, Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, true, true, InvInd, -1, WireIntersections );
241 SetAltitudePoints( points );
243 SetShape( DataTag_LeftBankShape, OutLeftB);
244 SetShape( DataTag_RightBankShape, OutRightB);
245 SetShape( DataTag_InletShape, OutInlet);
246 SetShape( DataTag_OutletShape, OutOutlet );
247 SetShape( DataTag_3DShape, Out3dPres );
248 SetShape( DataTag_2DShape, Out2dPres );
250 HYDROData_Bathymetry::Update();
253 void HYDROData_DTM::GetPresentationShapes( TopoDS_Shape& Out3dPres,
254 TopoDS_Shape& Out2dPres,
255 TopoDS_Shape& OutLeftB,
256 TopoDS_Shape& OutRightB,
257 TopoDS_Shape& OutInlet,
258 TopoDS_Shape& OutOutlet )
261 OutLeftB = GetShape( DataTag_LeftBankShape);
262 OutRightB = GetShape( DataTag_RightBankShape);
263 OutInlet = GetShape( DataTag_InletShape);
264 OutOutlet = GetShape( DataTag_OutletShape );
265 Out3dPres = GetShape( DataTag_3DShape );
266 Out2dPres = GetShape( DataTag_2DShape );
268 void HYDROData_DTM::CreateProfilesFromDTM (const HYDROData_SequenceOfObjects& InpProfiles,
271 AltitudePoints& points,
272 TopoDS_Shape& Out3dPres,
273 TopoDS_Shape& Out2dPres,
274 TopoDS_Shape& OutLeftB,
275 TopoDS_Shape& OutRightB,
276 TopoDS_Shape& OutInlet,
277 TopoDS_Shape& OutOutlet,
280 std::set<int>& InvInd,
282 bool& WireIntersections)
284 int aLower = InpProfiles.Lower(), anUpper = InpProfiles.Upper();
285 size_t n = anUpper - aLower + 1;
287 std::vector<Handle_HYDROData_Profile> profiles;
288 profiles.reserve( n );
289 for( int i=aLower; i<=anUpper; i++ )
291 Handle(HYDROData_Profile) aProfile = Handle(HYDROData_Profile)::DownCast( InpProfiles.Value( i ) );
292 if( !aProfile.IsNull() )
293 profiles.push_back( aProfile );
295 const double EPS = 1E-3;
297 AltitudePoints right;
298 std::vector<AltitudePoints> main_profiles;
300 if( thePntsLimit > 0 )
302 int aNbPoints = EstimateNbPoints( profiles, ddz, step );
303 if( aNbPoints < 0 || aNbPoints > thePntsLimit )
307 if( ddz>EPS && step>EPS )
308 CreateProfiles(profiles, ddz, step, left, right, points, main_profiles,
309 Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, Create3dPres, Create2dPres, InvInd, WireIntersections );
312 bool HYDROData_DTM::GetPlanarFaceFromBanks( const TopoDS_Edge& LB, const TopoDS_Edge& RB, TopoDS_Face& outF,
313 TopTools_SequenceOfShape* Boundr)
317 Handle_Geom_Plane refpl = new Geom_Plane(gp_Pnt(0,0,0), gp_Dir(0,0,1));
318 BB.MakeFace(F, refpl, Precision::Confusion());
319 BRepAlgo_NormalProjection nproj(F);
322 nproj.SetDefaultParams();
327 //TopoDS_Shape projBanks = nproj.Projection();
328 TopoDS_Vertex VFI, VLI, VFO, VLO;
329 TopoDS_Edge prLB = TopoDS::Edge(nproj.Generated(LB).First());
330 TopoDS_Edge prRB = TopoDS::Edge(nproj.Generated(RB).First());
331 TopExp::Vertices(prLB, VFI, VFO, 1);
332 TopExp::Vertices(prRB, VLI, VLO, 1);
333 TopoDS_Edge prIL = BRepLib_MakeEdge(VFI, VLI).Edge();
334 TopoDS_Edge prOL = BRepLib_MakeEdge(VFO, VLO).Edge();
335 TopoDS_Wire prW = BRepLib_MakeWire(prLB, prIL, prOL, prRB).Wire();
336 outF = BRepBuilderAPI_MakeFace(refpl->Pln(), prW, 1).Face();
340 Boundr->Append(prLB);
341 Boundr->Append(prIL);
342 Boundr->Append(prOL);
343 Boundr->Append(prRB);
346 ShapeAnalysis_Wire WA(prW, outF, Precision::Confusion());
347 bool res = WA.CheckSelfIntersection();
351 void HYDROData_DTM::CreateProfiles(const std::vector<Handle_HYDROData_Profile>& theProfiles,
353 double theSpatialStep,
354 AltitudePoints& theOutLeft,
355 AltitudePoints& theOutRight,
356 AltitudePoints& theOutPoints,
357 std::vector<AltitudePoints>& theOutMainProfiles,
358 TopoDS_Shape& Out3dPres,
359 TopoDS_Shape& Out2dPres,
360 TopoDS_Shape& OutLeftB,
361 TopoDS_Shape& OutRightB,
362 TopoDS_Shape& OutInlet,
363 TopoDS_Shape& OutOutlet,
366 std::set<int>& InvInd,
369 if (theProfiles.empty())
371 theOutPoints = Interpolate( theProfiles, theDDZ, theSpatialStep, theOutLeft, theOutRight, theOutMainProfiles, InvInd );
372 //note that if Create3dPres is false => Create2dPres flag is meaningless!
373 if( theOutPoints.empty() )
378 TopTools_IndexedMapOfOrientedShape ll = Create3DShape( theOutLeft, theOutRight, theOutMainProfiles);
384 BB.MakeCompound(cmp);
385 for (int i = 1; i <= ll.Extent(); i++)
390 //same order as in HYDROData_DTM::Update()
392 OutRightB = ll(ll.Extent());
394 OutOutlet = ll(ll.Extent() - 1);
399 ProjStat = GetPlanarFaceFromBanks(TopoDS::Edge(OutLeftB), TopoDS::Edge(OutRightB), outF, NULL);
408 void HYDROData_DTM::GetProperties( const Handle_HYDROData_Profile& theProfile,
409 gp_Pnt& theLowestPoint, gp_Vec2d& theDir,
410 double& theZMin, double& theZMax )
412 theLowestPoint = theProfile->GetBottomPoint();
415 theProfile->GetLeftPoint( aLeft, true, true );
416 theProfile->GetRightPoint( aRight, true, true );
417 double x = aRight.X()-aLeft.X();
418 double y = aRight.Y()-aLeft.Y();
419 theDir = gp_Vec2d( x, y );
421 HYDROData_Profile::ProfilePoints points = theProfile->GetProfilePoints();
422 int lo = points.Lower();
423 int up = points.Upper();
424 theZMin = std::numeric_limits<double>::max();
426 for( int i=lo; i<=up; i++ )
428 double z = points.Value( i ).Z();
436 inline gp_Pnt2d To2D( const gp_Pnt& thePnt, const gp_Trsf& theTr,
437 double& theUMin, double& theUMax )
439 gp_Pnt p = thePnt.Transformed( theTr );
446 return gp_Pnt2d( u, z );
449 Handle(TColgp_HArray1OfPnt2d) To2D( const TColgp_Array1OfPnt& thePoints,
450 const gp_Trsf& theTr,
451 double& theUMin, double& theUMax )
453 int low = thePoints.Lower(), up = thePoints.Upper();
454 Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( low, up );
455 for( int i=low; i<=up; i++ )
456 points->SetValue( i, To2D( thePoints.Value( i ), theTr, theUMin, theUMax ) );
460 Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve,
461 Standard_Real theFirst, Standard_Real theLast,
462 const gp_Trsf& theTr,
463 double& theUMin, double& theUMax )
465 if( theCurve->IsKind( STANDARD_TYPE( Geom_Line ) ) )
467 gp_Pnt aFirstPnt, aLastPnt;
468 theCurve->D0( theFirst, aFirstPnt );
469 theCurve->D0( theLast, aLastPnt );
472 aFirst2d = To2D( aFirstPnt, theTr, theUMin, theUMax ),
473 aLast2d = To2D( aLastPnt, theTr, theUMin, theUMax );
475 gp_Vec2d dir( aFirst2d, aLast2d );
476 Handle_Geom2d_Line aLine2d = new Geom2d_Line( aFirst2d, gp_Dir2d( dir.X(), dir.Y() ) );
477 return new Geom2d_TrimmedCurve( aLine2d, 0, aLast2d.Distance( aFirst2d ) );
480 if( theCurve->IsKind( STANDARD_TYPE( Geom_BSplineCurve ) ) )
482 Handle(Geom_BSplineCurve) aSpline = Handle(Geom_BSplineCurve)::DownCast( theCurve );
484 Handle(TColgp_HArray1OfPnt2d) poles = To2D( aSpline->Poles(), theTr, theUMin, theUMax );
485 const TColStd_Array1OfReal& knots = aSpline->Knots();
486 const TColStd_Array1OfInteger& multiplicities = aSpline->Multiplicities();
487 int aDegree = aSpline->Degree();
489 return new Geom2d_BSplineCurve( poles->Array1(), knots, multiplicities, aDegree );
492 return Handle(Geom2d_Curve)();
495 #include <GCE2d_MakeSegment.hxx>
496 #include <Geom2dAPI_InterCurveCurve.hxx>
497 bool IsCooriented( const Handle_HYDROData_Profile& theProfile1,
498 const Handle_HYDROData_Profile& theProfile2 )
500 if( theProfile1==theProfile2 )
503 gp_XY lp1, rp1, lp2, rp2;
504 theProfile1->GetLeftPoint(lp1);
505 theProfile1->GetRightPoint(rp1);
506 theProfile2->GetLeftPoint(lp2);
507 theProfile2->GetRightPoint(rp2);
509 GCE2d_MakeSegment s1(lp1, lp2);
510 GCE2d_MakeSegment s2(rp1, rp2);
512 Geom2dAPI_InterCurveCurve inter;
514 if (inter.NbPoints() == 0)
520 Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis(
521 const std::vector<Handle_HYDROData_Profile>& theProfiles,
522 std::vector<double>& theDistances )
524 size_t n = theProfiles.size();
526 return Handle_Geom2d_BSplineCurve();
528 Handle_Geom2d_BSplineCurve aResult;
530 Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( 1, (int)n );
531 TColgp_Array1OfVec2d tangents( 1, (int)n );
532 Handle(TColStd_HArray1OfBoolean) flags = new TColStd_HArray1OfBoolean( 1, (int)n );
534 for( size_t i = 1; i <= n; i++ )
536 Handle_HYDROData_Profile aProfile = theProfiles[i-1];
537 Handle_HYDROData_Profile aPrevProfile = i==1 ? theProfiles[i-1] : theProfiles[i-2];
538 Handle_HYDROData_Profile aNextProfile = i==n ? theProfiles[i-1] : theProfiles[i];
540 if( !IsCooriented( aProfile, aNextProfile ) )
543 aProfile->GetLeftPoint( lp, true );
544 aProfile->GetRightPoint( rp, true );
545 aProfile->SetLeftPoint( rp, true );
546 aProfile->SetRightPoint( lp, true );
554 gp_XYZ curP = aProfile->GetBottomPoint();
555 gp_XY curP2d = gp_XY(curP.X(), curP.Y());
559 nextP = aPrevProfile->GetBottomPoint(true);
561 nextP = aNextProfile->GetBottomPoint(true);
563 gp_XY nextP2d = gp_XY(nextP.X(), nextP.Y());
566 GetProperties( aProfile, aLowest, aPrTangent, zmin, zmax );
567 aNormal.SetCoord( -aPrTangent.Y(), aPrTangent.X() );
569 gp_Vec2d aDirToNextProfile(nextP2d.X() - curP2d.X(), nextP2d.Y() - curP2d.Y() );
571 aDirToNextProfile.Reverse();
572 if (aNormal.Dot(aDirToNextProfile) < 0)
577 points->SetValue( (int)i, gp_Pnt2d( aLowest.X(), aLowest.Y() ) );
578 tangents.SetValue( (int)i, aNormal );
579 flags->SetValue( (int)i, Standard_True );
582 Geom2dAPI_Interpolate anInterpolator( points, Standard_False, Standard_False );
583 anInterpolator.Load( tangents, flags );
584 anInterpolator.Perform();
585 if( anInterpolator.IsDone() )
587 aResult = anInterpolator.Curve();
589 //fill the distances vector
590 Geom2dAdaptor_Curve anAdaptor( aResult );
592 theDistances.clear();
593 theDistances.reserve( n );
594 Standard_Real aParamFirst = anAdaptor.FirstParameter(), aParamLast = anAdaptor.LastParameter();
595 for( size_t i = 1; i <= n; i++ )
597 gp_Pnt2d aPnt = points->Value( (Standard_Integer)i );
598 Geom2dAPI_ProjectPointOnCurve aProject( aPnt, aResult );
599 Standard_Real aParam = aProject.LowerDistanceParameter();
600 double aDistance = GCPnts_AbscissaPoint::Length( anAdaptor, aParamFirst, aParam );
601 theDistances.push_back( aDistance );
607 std::vector<Handle_Geom2d_Curve> HYDROData_DTM::ProfileToParametric(
608 const Handle_HYDROData_Profile& theProfile,
609 double& theUMin, double& theUMax, gp_Vec2d& theDir )
611 std::vector<Handle_Geom2d_Curve> curves;
613 // Transformation of the coordinate systems
616 GetProperties( theProfile, aLowest, theDir, zmin, zmax );
618 gp_Ax3 aStd3d( gp_Pnt( 0, 0, 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( 1, 0, 0 ) );
619 gp_Ax3 aLocal( gp_Pnt( aLowest.X(), aLowest.Y(), 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( theDir.X(), theDir.Y(), 0 ) );
622 aTransf.SetTransformation( aStd3d, aLocal );
624 // Iteration via edges
625 TopoDS_Wire aWire = TopoDS::Wire( theProfile->GetShape3D() );
626 TopExp_Explorer anExp( aWire, TopAbs_EDGE );
627 for( ; anExp.More(); anExp.Next() )
629 // Extract an edge from wire
630 TopoDS_Edge anEdge = TopoDS::Edge( anExp.Current() );
632 // Extract a curve corresponding to the edge
633 TopLoc_Location aLoc;
634 Standard_Real aFirst, aLast;
635 Handle(Geom_Curve) aCurve = BRep_Tool::Curve( anEdge, aLoc, aFirst, aLast );
637 // Convert the curve to 2d CS
638 Handle(Geom2d_Curve) aCurve2d = CurveTo2D( aCurve, aFirst, aLast, aTransf, theUMin, theUMax );
639 if( !aCurve2d.IsNull() )
640 curves.push_back( aCurve2d );
646 bool CalcMidWidth( const std::set<double>& intersections, double& theMid, double& theWid )
648 double umin = std::numeric_limits<double>::max(),
651 size_t n = intersections.size();
655 std::set<double>::const_iterator it = intersections.begin(), last = intersections.end();
656 for( ; it!=last; it++ )
664 theMid = ( umin+umax )/2;
669 void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& theProfile,
670 double theXCurv, double theMinZ, double theMaxZ, double theDDZ,
671 CurveUZ& theMidPointCurve,
672 CurveUZ& theWidthCurve,
673 int& intersection_nb,
676 double aDblMax = std::numeric_limits<double>::max(),
681 gp_Vec2d aProfileDir;
682 std::vector<Handle_Geom2d_Curve> curves = ProfileToParametric( theProfile, aUMin, aUMax, aProfileDir );
683 size_t n = curves.size();
688 // we add the "virtual" vertical lines to simulate the intersection with profile
689 gp_Pnt2d aFirst, aLast;
690 curves[0]->D0( curves[0]->FirstParameter(), aFirst );
691 curves[n-1]->D0( curves[n-1]->LastParameter(), aLast );
692 Handle(Geom2d_Line) aV1 = new Geom2d_Line( aFirst, gp_Dir2d( 0, 1 ) );
693 Handle(Geom2d_TrimmedCurve) aT1 = new Geom2d_TrimmedCurve( aV1, 0.0, aVMax );
695 Handle(Geom2d_Line) aV2 = new Geom2d_Line( aLast, gp_Dir2d( 0, 1 ) );
696 Handle(Geom2d_TrimmedCurve) aT2 = new Geom2d_TrimmedCurve( aV2, 0.0, aVMax );
698 curves.push_back( aT1 );
699 curves.push_back( aT2 );
701 int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 );
702 theMidPointCurve = CurveUZ( theXCurv, aProfileDir, theMinZ );
703 theMidPointCurve.reserve( psize );
704 theWidthCurve = CurveUZ( theXCurv, aProfileDir, theMinZ );
705 theWidthCurve.reserve( psize );
708 // for each discrete value of z we search intersection with profile
709 for( double z1 = theMinZ; z1 <= theMaxZ; z1 += theDDZ )
711 Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z1 ), gp_Dir2d( 1, 0 ) );
712 std::set<double> intersections;
713 for( size_t i = 0; i < n; i++ )
715 Handle_Geom2d_Curve aCurve = curves[i];
716 Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
717 for( int k=1, m=anIntersect.NbPoints(); k<=m; k++ )
718 intersections.insert( anIntersect.Point( k ).X() );
721 intersection_nb = intersections.size();
722 if( intersection_nb >= 1 )
725 if( !CalcMidWidth( intersections, u_mid, u_wid ) )
728 double z = z1 - theMinZ;
732 theMidPointCurve.push_back( p_mid );
737 theWidthCurve.push_back( p_wid );
742 void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCurveB,
743 int theNbSteps, std::vector<CurveUZ>& theInterpolation,
746 theInterpolation.clear();
747 int d = isAddSecond ? 2 : 1;
748 theInterpolation.reserve( theNbSteps+d );
749 double dt = 1.0 / double( theNbSteps + 1 );
751 theInterpolation.push_back( theCurveA );
752 for( int i=0; i<theNbSteps; i++, t+=dt )
754 CurveUZ anInterp = theCurveA*(1-t) + theCurveB*t;
755 theInterpolation.push_back( anInterp );
758 theInterpolation.push_back( theCurveB );
760 #include <BRepLib_MakeEdge2d.hxx>
761 void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
762 const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
763 AltitudePoints& thePoints )
765 Geom2dAdaptor_Curve anAdaptor( theHydraulicAxis );
766 TopoDS_Edge E2d = BRepLib_MakeEdge2d(theHydraulicAxis).Edge();
767 GCPnts_AbscissaPoint ap( anAdaptor, theMidCurve.Xcurv(), anAdaptor.FirstParameter() );
768 double aParam = ap.Parameter();
771 anAdaptor.D0( aParam, point );
772 gp_Vec2d profile_dir = theMidCurve.ProfileDir();
773 //gp_Dir tangent_n( -profile_dir.Y(), profile_dir.X(), dz );
774 profile_dir.Normalize();
776 size_t n = theMidCurve.size();
777 std::map<double, AltitudePoint> sorted_points;
778 for( size_t i=0; i<n; i++ )
780 double param1 = theMidCurve[i].U - theWidthCurve[i].U / 2;
781 double param2 = theMidCurve[i].U + theWidthCurve[i].U / 2;
783 gp_Pnt2d p1 = point.Translated( param1 * profile_dir);
784 gp_Pnt2d p2 = point.Translated( param2 * profile_dir);
786 double z = theMidCurve[i].Z + theMidCurve.DeltaZ();
788 AltitudePoint p3d_1( p1.X(), p1.Y(), z ), p3d_2( p2.X(), p2.Y(), z );
790 sorted_points[param1] = p3d_1;
791 sorted_points[param2] = p3d_2;
794 thePoints.reserve( sorted_points.size() );
795 const double EPS = 1E-12;
796 std::map<double, AltitudePoint>::const_iterator it = sorted_points.begin(), last = sorted_points.end();
797 for( ; it!=last; it++ )
798 if( thePoints.empty() || thePoints.back().SquareDistance( it->second ) > EPS )
799 thePoints.push_back( it->second );
802 inline double max( double a, double b )
810 inline double min( double a, double b )
818 #include <BRepLib_MakeWire.hxx>
820 std::vector<HYDROData_Bathymetry::AltitudePoints> HYDROData_DTM::Interpolate
821 ( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
822 const Handle_HYDROData_Profile& theProfileA,
824 const Handle_HYDROData_Profile& theProfileB,
826 double theDDZ, int theNbSteps, bool isAddSecond,
827 int& inter_nb_1, int& inter_nb_2)
829 double zminA, zmaxA, zminB, zmaxB;
830 gp_Pnt lowestA, lowestB;
833 GetProperties( theProfileA, lowestA, dirA, zminA, zmaxA );
834 GetProperties( theProfileB, lowestB, dirB, zminB, zmaxB );
837 double hmax = max( zmaxA-zminA, zmaxB-zminB );
839 //double dz = zminB - zminA;
840 //double zmin = min( zminA, zminB );
841 //double zmax = max( zmaxA, zmaxB );
843 CurveUZ midA(0, gp_Vec2d(), 0), midB(0, gp_Vec2d(), 0);
844 CurveUZ widA(0, gp_Vec2d(), 0), widB(0, gp_Vec2d(), 0);
846 ProfileDiscretization( theProfileA, theXCurvA, zminA, zminA+hmax, theDDZ, midA, widA, inter_nb_1 );
847 ProfileDiscretization( theProfileB, theXCurvB, zminB, zminB+hmax, theDDZ, midB, widB, inter_nb_2 );
849 std::vector<CurveUZ> mid, wid;
850 Interpolate( midA, midB, theNbSteps, mid, isAddSecond );
851 Interpolate( widA, widB, theNbSteps, wid, isAddSecond );
853 size_t p = mid.size();
854 size_t q = p>0 ? 2*mid[0].size() : 1;
855 std::vector<AltitudePoints> points;
858 for( size_t i=0; i<p; i++ )
860 points[i].reserve( q );
861 CurveTo3D( theHydraulicAxis, mid[i], wid[i], points[i] );
867 HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
868 ( const std::vector<Handle_HYDROData_Profile>& theProfiles,
869 double theDDZ, double theSpatialStep,
870 AltitudePoints& theLeft,
871 AltitudePoints& theRight,
872 std::vector<AltitudePoints>& theMainProfiles,
873 std::set<int>& invalInd)
875 AltitudePoints points;
876 size_t n = theProfiles.size();
880 std::vector<double> distances;
881 Handle_Geom2d_BSplineCurve aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
882 if( aHydraulicAxis.IsNull() )
885 theMainProfiles.reserve( n );
887 for( size_t i=0, n1=n-1; i<n1; i++ )
889 double aDistance = distances[i+1]-distances[i];
890 int aNbSteps = int(aDistance/theSpatialStep);
891 bool isAddSecond = i==n1-1;
893 // 1. Calculate interpolated profiles
894 int inter_nb_1, inter_nb_2;
895 std::vector<AltitudePoints> local_points = Interpolate( aHydraulicAxis, theProfiles[i], distances[i],
896 theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond, inter_nb_1, inter_nb_2 );
897 int lps = local_points.size();
903 invalInd.insert(i+1);
905 // 2. Put all points into the global container
906 for( size_t j=0; j<lps; j++ )
908 const AltitudePoints& lp = local_points[j];
910 points.reserve( lp.size() * n );
911 for( size_t k=0, ks=lp.size(); k<ks; k++ )
912 points.push_back( lp[k] );
915 // 3. Get left/right banks' points
918 theLeft.reserve( lps * n );
919 theRight.reserve( lps * n );
921 for( size_t j=0; j<lps; j++ )
923 const AltitudePoints& lp = local_points[j];
924 theLeft.push_back( lp[0] );
925 theRight.push_back( lp[lp.size()-1] );
928 // 4. Get main profiles points
929 theMainProfiles.push_back( local_points[0] );
931 theMainProfiles.push_back( local_points[lps-1] );
936 int HYDROData_DTM::EstimateNbPoints( const std::vector<Handle_HYDROData_Profile>& theProfiles,
937 double theDDZ, double theSpatialStep )
939 size_t n = theProfiles.size();
942 if( theDDZ<1E-6 || theSpatialStep<1E-6 )
945 std::vector<double> distances;
946 Handle_Geom2d_BSplineCurve aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
947 if( aHydraulicAxis.IsNull() )
950 double aCompleteDistance = distances[n-1];
951 int aNbSteps = int( aCompleteDistance / theSpatialStep ) + 1;
955 GetProperties( theProfiles[0], aLowest, aDir, aZMin, aZMax );
956 int aNbZSteps = (aZMax-aZMin)/theDDZ;
958 if( aNbSteps > ( 1<<16 ) || aNbZSteps > ( 1<<16 ) )
961 return aNbSteps * aNbZSteps;