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 IMPLEMENT_STANDARD_HANDLE( HYDROData_DTM, HYDROData_Bathymetry )
35 IMPLEMENT_STANDARD_RTTIEXT( HYDROData_DTM, HYDROData_Bathymetry )
38 HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv )
43 HYDROData_DTM::CurveUZ::~CurveUZ()
47 double HYDROData_DTM::CurveUZ::Xcurv() const
52 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator + ( const CurveUZ& c ) const
54 HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv() );
57 for( int i=0; i<n; i++ )
60 p.U = operator[]( i ).U + c[i].U;
61 p.Z = operator[]( i ).Z;
67 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator * ( double d ) const
69 HYDROData_DTM::CurveUZ res( Xcurv()*d );
72 for( int i=0; i<n; i++ )
75 p.U = operator[]( i ).U * d;
76 p.Z = operator[]( i ).Z;
82 void HYDROData_DTM::Bank::reserve( int theNbPoints )
84 myPoints.reserve( theNbPoints );
85 myDirs.reserve( theNbPoints );
88 void HYDROData_DTM::Bank::push_back( const gp_Pnt& thePnt, const gp_Dir& theTangent )
90 myPoints.push_back( thePnt );
91 myDirs.push_back( theTangent );
94 void HYDROData_DTM::Bank::clear()
100 TopoDS_Edge HYDROData_DTM::Bank::createEdge( bool is3d ) const
102 size_t n = myPoints.size();
103 Handle_Geom_BSplineCurve aCurve;
105 Handle(TColgp_HArray1OfPnt) points = new TColgp_HArray1OfPnt( 1, (int)n );
106 TColgp_Array1OfVec tangents( 1, (int)n );
107 Handle(TColStd_HArray1OfBoolean) flags = new TColStd_HArray1OfBoolean( 1, (int)n );
109 for( size_t i = 1; i <= n; i++ )
111 gp_Pnt aPnt = myPoints[i-1];
114 gp_Vec aVec = myDirs[i-1];
117 points->SetValue( (int)i, aPnt );
118 tangents.SetValue( (int)i, aVec );
119 flags->SetValue( (int)i, Standard_True );
122 GeomAPI_Interpolate anInterpolator( points, Standard_False, Standard_False );
123 anInterpolator.Load( tangents, flags );
124 anInterpolator.Perform();
125 if( anInterpolator.IsDone() )
127 aCurve = anInterpolator.Curve();
128 return BRepBuilderAPI_MakeEdge( aCurve ).Edge();
131 return TopoDS_Edge();
136 HYDROData_DTM::HYDROData_DTM()
140 HYDROData_DTM::~HYDROData_DTM()
144 HYDROData_SequenceOfObjects HYDROData_DTM::GetProfiles() const
146 return GetReferenceObjects( DataTag_Profiles );
149 void HYDROData_DTM::SetProfiles( const HYDROData_SequenceOfObjects& theProfiles )
151 SetReferenceObjects( theProfiles, DataTag_Profiles );
154 double HYDROData_DTM::GetDDZ() const
156 return GetDouble( DataTag_DDZ );
159 void HYDROData_DTM::SetDDZ( double theDDZ )
161 SetDouble( DataTag_DDZ, theDDZ );
164 double HYDROData_DTM::GetSpatialStep() const
166 return GetDouble( DataTag_SpatialStep );
169 void HYDROData_DTM::SetSpatialStep( double theSpatialStep )
171 SetDouble( DataTag_SpatialStep, theSpatialStep );
174 void HYDROData_DTM::Update()
176 HYDROData_SequenceOfObjects objs = GetProfiles();
177 int aLower = objs.Lower(), anUpper = objs.Upper();
178 size_t n = anUpper-aLower+1;
180 std::vector<Handle_HYDROData_Profile> profiles;
181 profiles.reserve( n );
182 for( int i=aLower; i<=anUpper; i++ )
184 Handle(HYDROData_Profile) aProfile = Handle(HYDROData_Profile)::DownCast( objs.Value( i ) );
185 if( !aProfile.IsNull() )
186 profiles.push_back( aProfile );
189 double ddz = GetDDZ();
190 double step = GetSpatialStep();
191 const double EPS = 1E-3;
192 AltitudePoints points;
196 if( ddz>EPS && step>EPS )
197 points = Interpolate( profiles, ddz, step, &myLeft, &myRight );
198 SetAltitudePoints( points );
210 void HYDROData_DTM::GetProperties( const Handle_HYDROData_Profile& theProfile,
211 gp_Pnt& theLowestPoint, gp_Vec2d& theDir,
213 double& theZMin, double& theZMax )
215 theLowestPoint = theProfile->GetBottomPoint();
218 theProfile->GetLeftPoint( aLeft, true, true );
219 theProfile->GetRightPoint( aRight, true, true );
220 double x = aRight.X()-aLeft.X();
221 double y = aRight.Y()-aLeft.Y();
223 theDir = gp_Vec2d( -y, x );
225 theDir = gp_Vec2d( x, y );
227 HYDROData_Profile::ProfilePoints points = theProfile->GetProfilePoints();
228 int lo = points.Lower();
229 int up = points.Upper();
230 theZMin = std::numeric_limits<double>::max();
232 for( int i=lo; i<=up; i++ )
234 double z = points.Value( i ).Z();
242 inline gp_Pnt2d To2D( const gp_Pnt& thePnt, const gp_Trsf& theTr,
243 double& theUMin, double& theUMax )
245 gp_Pnt p = thePnt.Transformed( theTr );
252 return gp_Pnt2d( u, z );
255 Handle(TColgp_HArray1OfPnt2d) To2D( const TColgp_Array1OfPnt& thePoints,
256 const gp_Trsf& theTr,
257 double& theUMin, double& theUMax )
259 int low = thePoints.Lower(), up = thePoints.Upper();
260 Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( low, up );
261 for( int i=low; i<=up; i++ )
262 points->SetValue( i, To2D( thePoints.Value( i ), theTr, theUMin, theUMax ) );
266 Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve,
267 Standard_Real theFirst, Standard_Real theLast,
268 const gp_Trsf& theTr,
269 double& theUMin, double& theUMax )
271 if( theCurve->IsKind( STANDARD_TYPE( Geom_Line ) ) )
273 gp_Pnt aFirstPnt, aLastPnt;
274 theCurve->D0( theFirst, aFirstPnt );
275 theCurve->D0( theLast, aLastPnt );
278 aFirst2d = To2D( aFirstPnt, theTr, theUMin, theUMax ),
279 aLast2d = To2D( aLastPnt, theTr, theUMin, theUMax );
281 gp_Vec2d dir( aFirst2d, aLast2d );
282 Handle_Geom2d_Line aLine2d = new Geom2d_Line( aFirst2d, gp_Dir2d( dir.X(), dir.Y() ) );
283 return new Geom2d_TrimmedCurve( aLine2d, 0, aLast2d.Distance( aFirst2d ) );
286 if( theCurve->IsKind( STANDARD_TYPE( Geom_BSplineCurve ) ) )
288 Handle(Geom_BSplineCurve) aSpline = Handle(Geom_BSplineCurve)::DownCast( theCurve );
290 Handle(TColgp_HArray1OfPnt2d) poles = To2D( aSpline->Poles(), theTr, theUMin, theUMax );
291 const TColStd_Array1OfReal& knots = aSpline->Knots();
292 const TColStd_Array1OfInteger& multiplicities = aSpline->Multiplicities();
293 int aDegree = aSpline->Degree();
295 return new Geom2d_BSplineCurve( poles->Array1(), knots, multiplicities, aDegree );
298 return Handle(Geom2d_Curve)();
301 Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis(
302 const std::vector<Handle_HYDROData_Profile>& theProfiles,
303 std::vector<double>& theDistances )
305 size_t n = theProfiles.size();
306 Handle_Geom2d_BSplineCurve aResult;
308 Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( 1, (int)n );
309 TColgp_Array1OfVec2d tangents( 1, (int)n );
310 Handle(TColStd_HArray1OfBoolean) flags = new TColStd_HArray1OfBoolean( 1, (int)n );
312 for( size_t i = 1; i <= n; i++ )
314 Handle_HYDROData_Profile aProfile = theProfiles[i-1];
320 GetProperties( aProfile, aLowest, aTangent, true, zmin, zmax );
321 aTangent.Normalize();
323 points->SetValue( (int)i, gp_Pnt2d( aLowest.X(), aLowest.Y() ) );
324 tangents.SetValue( (int)i, aTangent );
325 flags->SetValue( (int)i, Standard_True );
328 Geom2dAPI_Interpolate anInterpolator( points, Standard_False, Standard_False );
329 anInterpolator.Load( tangents, flags );
330 anInterpolator.Perform();
331 if( anInterpolator.IsDone() )
333 aResult = anInterpolator.Curve();
335 //fill the distances vector
336 Geom2dAdaptor_Curve anAdaptor( aResult );
338 theDistances.clear();
339 theDistances.reserve( n );
340 Standard_Real aParamFirst = anAdaptor.FirstParameter(), aParamLast = anAdaptor.LastParameter();
341 for( size_t i = 1; i <= n; i++ )
343 gp_Pnt2d aPnt = points->Value( (Standard_Integer)i );
344 Geom2dAPI_ProjectPointOnCurve aProject( aPnt, aResult );
345 Standard_Real aParam = aProject.LowerDistanceParameter();
346 double aDistance = GCPnts_AbscissaPoint::Length( anAdaptor, aParamFirst, aParam );
347 theDistances.push_back( aDistance );
353 std::vector<Handle_Geom2d_Curve> HYDROData_DTM::ProfileToParametric(
354 const Handle_HYDROData_Profile& theProfile,
355 double& theUMin, double& theUMax )
357 std::vector<Handle_Geom2d_Curve> curves;
359 // Transformation of the coordinate systems
363 GetProperties( theProfile, aLowest, aDir, false, zmin, zmax );
365 gp_Ax3 aStd3d( gp_Pnt( 0, 0, 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( 1, 0, 0 ) );
366 gp_Ax3 aLocal( aLowest, gp_Dir( 0, 0, 1 ), gp_Dir( aDir.X(), aDir.Y(), 0 ) );
369 aTransf.SetTransformation( aStd3d, aLocal );
371 // Iteration via edges
372 TopoDS_Wire aWire = TopoDS::Wire( theProfile->GetShape3D() );
373 TopExp_Explorer anExp( aWire, TopAbs_EDGE );
374 for( ; anExp.More(); anExp.Next() )
376 // Extract an edge from wire
377 TopoDS_Edge anEdge = TopoDS::Edge( anExp.Current() );
379 // Extract a curve corresponding to the edge
380 TopLoc_Location aLoc;
381 Standard_Real aFirst, aLast;
382 Handle(Geom_Curve) aCurve = BRep_Tool::Curve( anEdge, aLoc, aFirst, aLast );
384 // Convert the curve to 2d CS
385 Handle(Geom2d_Curve) aCurve2d = CurveTo2D( aCurve, aFirst, aLast, aTransf, theUMin, theUMax );
386 if( !aCurve2d.IsNull() )
387 curves.push_back( aCurve2d );
393 double CalcGC( const std::vector<gp_Pnt2d>& intersections )
396 size_t n = intersections.size();
397 for( size_t i = 0; i < n; i++ )
399 u += intersections[i].X();
405 double CalcWidth( const std::vector<gp_Pnt2d>& intersections )
407 double umin = std::numeric_limits<double>::max(),
410 size_t n = intersections.size();
414 for( size_t i = 0; i < n; i++ )
416 double u = intersections[i].X();
425 void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& theProfile,
426 double theXCurv, double theMinZ, double theMaxZ, double theDDZ,
427 CurveUZ& theMidPointCurve,
428 CurveUZ& theWidthCurve,
429 double theTolerance )
431 double aDblMax = std::numeric_limits<double>::max(),
436 std::vector<Handle_Geom2d_Curve> curves = ProfileToParametric( theProfile, aUMin, aUMax );
437 size_t n = curves.size();
442 // we add the "virtual" vertical lines to simulate the intersection with profile
443 gp_Pnt2d aFirst, aLast;
444 curves[0]->D0( curves[0]->FirstParameter(), aFirst );
445 curves[n-1]->D0( curves[n-1]->LastParameter(), aLast );
446 Handle(Geom2d_Line) aV1 = new Geom2d_Line( aFirst, gp_Dir2d( 0, 1 ) );
447 Handle(Geom2d_TrimmedCurve) aT1 = new Geom2d_TrimmedCurve( aV1, 0.0, aVMax );
449 Handle(Geom2d_Line) aV2 = new Geom2d_Line( aLast, gp_Dir2d( 0, 1 ) );
450 Handle(Geom2d_TrimmedCurve) aT2 = new Geom2d_TrimmedCurve( aV2, 0.0, aVMax );
452 curves.push_back( aT1 );
453 curves.push_back( aT2 );
455 int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 );
456 theMidPointCurve = CurveUZ( theXCurv );
457 theMidPointCurve.reserve( psize );
458 theWidthCurve = CurveUZ( theXCurv );
459 theWidthCurve.reserve( psize );
462 // for each discrete value of z we search intersection with profile
463 for( double z = theMinZ; z <= theMaxZ; z += theDDZ )
465 Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z ), gp_Dir2d( 1, 0 ) );
466 std::vector<gp_Pnt2d> intersections;
467 for( size_t i = 0; i < n; i++ )
469 Handle_Geom2d_Curve aCurve = curves[i];
470 Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
471 for( int k=1, m=anIntersect.NbPoints(); k<=m; k++ )
472 intersections.push_back( anIntersect.Point( k ) );
475 if( intersections.size() >= 2 )
477 double u_mid = CalcGC( intersections );
481 theMidPointCurve.push_back( p_mid );
483 double u_wid = CalcWidth( intersections );
487 theWidthCurve.push_back( p_wid );
492 void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCurveB,
493 int theNbSteps, std::vector<CurveUZ>& theInterpolation,
496 theInterpolation.clear();
497 int d = isAddSecond ? 2 : 1;
498 theInterpolation.reserve( theNbSteps+d );
499 double dt = 1.0 / double( theNbSteps + 1 );
501 theInterpolation.push_back( theCurveA );
502 for( int i=0; i<theNbSteps; i++, t+=dt )
504 CurveUZ anInterp = theCurveA*(1-t) + theCurveB*t;
505 theInterpolation.push_back( anInterp );
508 theInterpolation.push_back( theCurveB );
511 void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
512 const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
513 AltitudePoints& thePoints,
514 Bank* theLeft, Bank* theRight, double dz )
516 Geom2dAdaptor_Curve anAdaptor( theHydraulicAxis );
517 GCPnts_AbscissaPoint ap( anAdaptor, theMidCurve.Xcurv(), anAdaptor.FirstParameter() );
518 double aParam = ap.Parameter();
521 gp_Vec2d tangent, profile_dir;
522 anAdaptor.D1( aParam, point, tangent );
523 profile_dir.SetCoord( tangent.Y(), -tangent.X() );
524 profile_dir.Normalize();
525 gp_Dir tangent_n( tangent.X(), tangent.Y(), dz );
527 size_t n = theMidCurve.size();
528 double min_param = 1E+15;
529 double max_param = -1E+15;
531 for( size_t i=0; i<n; i++ )
533 double param1 = theMidCurve[i].U - theWidthCurve[i].U / 2;
534 double param2 = theMidCurve[i].U + theWidthCurve[i].U / 2;
536 gp_Pnt2d p1 = point.Translated( param1 * profile_dir / 2 );
537 gp_Pnt2d p2 = point.Translated( param2 * profile_dir / 2 );
539 double z = theMidCurve[i].Z;
541 if( param1 < min_param )
546 if( param2 < min_param )
551 if( param1 > max_param )
556 if( param2 > max_param )
562 AltitudePoint p3d_1( p1.X(), p1.Y(), z ), p3d_2( p2.X(), p2.Y(), z );
563 thePoints.push_back( p3d_1 );
564 thePoints.push_back( p3d_2 );
569 gp_Pnt2d left2d = point.Translated( min_param * profile_dir / 2 );
570 gp_Pnt left( left2d.X(), left2d.Y(), z1 );
571 theLeft->push_back( left, tangent_n );
575 gp_Pnt2d right2d = point.Translated( max_param * profile_dir / 2 );
576 gp_Pnt right( right2d.X(), right2d.Y(), z2 );
577 theRight->push_back( right, tangent_n );
581 inline double max( double a, double b )
589 HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
590 ( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
591 const Handle_HYDROData_Profile& theProfileA,
593 const Handle_HYDROData_Profile& theProfileB,
595 double theDDZ, int theNbSteps, bool isAddSecond,
596 Bank* theLeft, Bank* theRight )
598 double zminA, zmaxA, zminB, zmaxB;
599 gp_Pnt lowestA, lowestB;
602 GetProperties( theProfileA, lowestA, dirA, false, zminA, zmaxA );
603 GetProperties( theProfileB, lowestB, dirB, false, zminB, zmaxB );
605 double dz = zminB - zminA;
607 double zmin = max( zminA, zminB );
608 double zmax = max( zmaxA, zmaxB );
610 CurveUZ midA(0), midB(0);
611 CurveUZ widA(0), widB(0);
613 ProfileDiscretization( theProfileA, theXCurvA, zmin, zmax, theDDZ, midA, widA );
614 ProfileDiscretization( theProfileB, theXCurvB, zmin, zmax, theDDZ, midB, widB );
616 std::vector<CurveUZ> mid, wid;
617 Interpolate( midA, midB, theNbSteps, mid, isAddSecond );
618 Interpolate( widA, widB, theNbSteps, wid, isAddSecond );
620 size_t p = mid.size();
621 size_t q = p>0 ? 2*mid[0].size() : 1;
622 AltitudePoints points;
623 points.reserve( p*q );
624 for( size_t i=0; i<p; i++ )
625 CurveTo3D( theHydraulicAxis, mid[i], wid[i], points, theLeft, theRight, dz );
630 HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
631 ( const std::vector<Handle_HYDROData_Profile>& theProfiles,
632 double theDDZ, double theSpatialStep,
636 AltitudePoints points;
637 size_t n = theProfiles.size();
641 std::vector<double> distances;
642 Handle_Geom2d_BSplineCurve aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
643 if( aHydraulicAxis.IsNull() )
646 int aNbStepsComplete = 0;
647 for( size_t i=0, n1=n-1; i<n1; i++ )
649 double aDistance = distances[i+1]-distances[i];
650 aNbStepsComplete += ( int(aDistance/theSpatialStep) + 1 );
653 theLeft->reserve( aNbStepsComplete );
655 theRight->reserve( aNbStepsComplete );
657 for( size_t i=0, n1=n-1; i<n1; i++ )
659 double aDistance = distances[i+1]-distances[i];
660 int aNbSteps = int(aDistance/theSpatialStep) + 1;
661 bool isAddSecond = i==n1-1;
663 AltitudePoints local_points = Interpolate( aHydraulicAxis, theProfiles[i], distances[i],
664 theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond, theLeft, theRight );
667 points.reserve( local_points.size() * ( n-1 ) );
669 for( size_t j=0, m=local_points.size(); j<m; j++ )
670 points.push_back( local_points[j] );