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::createEdge3d() 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];
112 gp_Vec aVec = myDirs[i-1];
113 points->SetValue( (int)i, aPnt );
114 tangents.SetValue( (int)i, aVec );
115 flags->SetValue( (int)i, Standard_True );
118 GeomAPI_Interpolate anInterpolator( points, Standard_False, Standard_False );
119 anInterpolator.Load( tangents, flags );
120 anInterpolator.Perform();
121 if( anInterpolator.IsDone() )
123 aCurve = anInterpolator.Curve();
124 return BRepBuilderAPI_MakeEdge( aCurve ).Edge();
127 return TopoDS_Edge();
132 HYDROData_DTM::HYDROData_DTM()
136 HYDROData_DTM::~HYDROData_DTM()
140 HYDROData_SequenceOfObjects HYDROData_DTM::GetProfiles() const
142 return GetReferenceObjects( DataTag_Profiles );
145 void HYDROData_DTM::SetProfiles( const HYDROData_SequenceOfObjects& theProfiles )
147 SetReferenceObjects( theProfiles, DataTag_Profiles );
150 double HYDROData_DTM::GetDDZ() const
152 return GetDouble( DataTag_DDZ );
155 void HYDROData_DTM::SetDDZ( double theDDZ )
157 SetDouble( DataTag_DDZ, theDDZ );
160 double HYDROData_DTM::GetSpatialStep() const
162 return GetDouble( DataTag_SpatialStep );
165 void HYDROData_DTM::SetSpatialStep( double theSpatialStep )
167 SetDouble( DataTag_SpatialStep, theSpatialStep );
170 void HYDROData_DTM::Update()
172 HYDROData_SequenceOfObjects objs = GetProfiles();
173 int aLower = objs.Lower(), anUpper = objs.Upper();
174 size_t n = anUpper-aLower+1;
176 std::vector<Handle_HYDROData_Profile> profiles;
177 profiles.reserve( n );
178 for( int i=aLower; i<=anUpper; i++ )
180 Handle(HYDROData_Profile) aProfile = Handle(HYDROData_Profile)::DownCast( objs.Value( i ) );
181 if( !aProfile.IsNull() )
182 profiles.push_back( aProfile );
185 double ddz = GetDDZ();
186 double step = GetSpatialStep();
187 const double EPS = 1E-3;
188 AltitudePoints points;
192 if( ddz>EPS && step>EPS )
193 points = Interpolate( profiles, ddz, step, &myLeft, &myRight );
194 SetAltitudePoints( points );
206 void HYDROData_DTM::GetProperties( const Handle_HYDROData_Profile& theProfile,
207 gp_Pnt& theLowestPoint, gp_Vec2d& theDir,
209 double& theZMin, double& theZMax )
211 theLowestPoint = theProfile->GetBottomPoint();
214 theProfile->GetLeftPoint( aLeft, true, true );
215 theProfile->GetRightPoint( aRight, true, true );
216 double x = aRight.X()-aLeft.X();
217 double y = aRight.Y()-aLeft.Y();
219 theDir = gp_Vec2d( -y, x );
221 theDir = gp_Vec2d( x, y );
223 HYDROData_Profile::ProfilePoints points = theProfile->GetProfilePoints();
224 int lo = points.Lower();
225 int up = points.Upper();
226 theZMin = std::numeric_limits<double>::max();
228 for( int i=lo; i<=up; i++ )
230 double z = points.Value( i ).Z();
238 inline gp_Pnt2d To2D( const gp_Pnt& thePnt, const gp_Trsf& theTr,
239 double& theUMin, double& theUMax )
241 gp_Pnt p = thePnt.Transformed( theTr );
248 return gp_Pnt2d( u, z );
251 Handle(TColgp_HArray1OfPnt2d) To2D( const TColgp_Array1OfPnt& thePoints,
252 const gp_Trsf& theTr,
253 double& theUMin, double& theUMax )
255 int low = thePoints.Lower(), up = thePoints.Upper();
256 Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( low, up );
257 for( int i=low; i<=up; i++ )
258 points->SetValue( i, To2D( thePoints.Value( i ), theTr, theUMin, theUMax ) );
262 Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve,
263 Standard_Real theFirst, Standard_Real theLast,
264 const gp_Trsf& theTr,
265 double& theUMin, double& theUMax )
267 if( theCurve->IsKind( STANDARD_TYPE( Geom_Line ) ) )
269 gp_Pnt aFirstPnt, aLastPnt;
270 theCurve->D0( theFirst, aFirstPnt );
271 theCurve->D0( theLast, aLastPnt );
274 aFirst2d = To2D( aFirstPnt, theTr, theUMin, theUMax ),
275 aLast2d = To2D( aLastPnt, theTr, theUMin, theUMax );
277 gp_Vec2d dir( aFirst2d, aLast2d );
278 Handle_Geom2d_Line aLine2d = new Geom2d_Line( aFirst2d, gp_Dir2d( dir.X(), dir.Y() ) );
279 return new Geom2d_TrimmedCurve( aLine2d, 0, aLast2d.Distance( aFirst2d ) );
282 if( theCurve->IsKind( STANDARD_TYPE( Geom_BSplineCurve ) ) )
284 Handle(Geom_BSplineCurve) aSpline = Handle(Geom_BSplineCurve)::DownCast( theCurve );
286 Handle(TColgp_HArray1OfPnt2d) poles = To2D( aSpline->Poles(), theTr, theUMin, theUMax );
287 const TColStd_Array1OfReal& knots = aSpline->Knots();
288 const TColStd_Array1OfInteger& multiplicities = aSpline->Multiplicities();
289 int aDegree = aSpline->Degree();
291 return new Geom2d_BSplineCurve( poles->Array1(), knots, multiplicities, aDegree );
294 return Handle(Geom2d_Curve)();
297 Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis(
298 const std::vector<Handle_HYDROData_Profile>& theProfiles,
299 std::vector<double>& theDistances )
301 size_t n = theProfiles.size();
302 Handle_Geom2d_BSplineCurve aResult;
304 Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( 1, (int)n );
305 TColgp_Array1OfVec2d tangents( 1, (int)n );
306 Handle(TColStd_HArray1OfBoolean) flags = new TColStd_HArray1OfBoolean( 1, (int)n );
308 for( size_t i = 1; i <= n; i++ )
310 Handle_HYDROData_Profile aProfile = theProfiles[i-1];
316 GetProperties( aProfile, aLowest, aTangent, true, zmin, zmax );
317 aTangent.Normalize();
319 points->SetValue( (int)i, gp_Pnt2d( aLowest.X(), aLowest.Y() ) );
320 tangents.SetValue( (int)i, aTangent );
321 flags->SetValue( (int)i, Standard_True );
324 Geom2dAPI_Interpolate anInterpolator( points, Standard_False, Standard_False );
325 anInterpolator.Load( tangents, flags );
326 anInterpolator.Perform();
327 if( anInterpolator.IsDone() )
329 aResult = anInterpolator.Curve();
331 //fill the distances vector
332 Geom2dAdaptor_Curve anAdaptor( aResult );
334 theDistances.clear();
335 theDistances.reserve( n );
336 Standard_Real aParamFirst = anAdaptor.FirstParameter(), aParamLast = anAdaptor.LastParameter();
337 for( size_t i = 1; i <= n; i++ )
339 gp_Pnt2d aPnt = points->Value( (Standard_Integer)i );
340 Geom2dAPI_ProjectPointOnCurve aProject( aPnt, aResult );
341 Standard_Real aParam = aProject.LowerDistanceParameter();
342 double aDistance = GCPnts_AbscissaPoint::Length( anAdaptor, aParamFirst, aParam );
343 theDistances.push_back( aDistance );
349 std::vector<Handle_Geom2d_Curve> HYDROData_DTM::ProfileToParametric(
350 const Handle_HYDROData_Profile& theProfile,
351 double& theUMin, double& theUMax )
353 std::vector<Handle_Geom2d_Curve> curves;
355 // Transformation of the coordinate systems
359 GetProperties( theProfile, aLowest, aDir, false, zmin, zmax );
361 gp_Ax3 aStd3d( gp_Pnt( 0, 0, 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( 1, 0, 0 ) );
362 gp_Ax3 aLocal( aLowest, gp_Dir( 0, 0, 1 ), gp_Dir( aDir.X(), aDir.Y(), 0 ) );
365 aTransf.SetTransformation( aStd3d, aLocal );
367 // Iteration via edges
368 TopoDS_Wire aWire = TopoDS::Wire( theProfile->GetShape3D() );
369 TopExp_Explorer anExp( aWire, TopAbs_EDGE );
370 for( ; anExp.More(); anExp.Next() )
372 // Extract an edge from wire
373 TopoDS_Edge anEdge = TopoDS::Edge( anExp.Current() );
375 // Extract a curve corresponding to the edge
376 TopLoc_Location aLoc;
377 Standard_Real aFirst, aLast;
378 Handle(Geom_Curve) aCurve = BRep_Tool::Curve( anEdge, aLoc, aFirst, aLast );
380 // Convert the curve to 2d CS
381 Handle(Geom2d_Curve) aCurve2d = CurveTo2D( aCurve, aFirst, aLast, aTransf, theUMin, theUMax );
382 if( !aCurve2d.IsNull() )
383 curves.push_back( aCurve2d );
389 double CalcGC( const std::vector<gp_Pnt2d>& intersections )
392 size_t n = intersections.size();
393 for( size_t i = 0; i < n; i++ )
395 u += intersections[i].X();
401 double CalcWidth( const std::vector<gp_Pnt2d>& intersections )
403 double umin = std::numeric_limits<double>::max(),
406 size_t n = intersections.size();
410 for( size_t i = 0; i < n; i++ )
412 double u = intersections[i].X();
421 void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& theProfile,
422 double theXCurv, double theMinZ, double theMaxZ, double theDDZ,
423 CurveUZ& theMidPointCurve,
424 CurveUZ& theWidthCurve,
425 double theTolerance )
427 double aDblMax = std::numeric_limits<double>::max(),
432 std::vector<Handle_Geom2d_Curve> curves = ProfileToParametric( theProfile, aUMin, aUMax );
433 size_t n = curves.size();
438 // we add the "virtual" vertical lines to simulate the intersection with profile
439 gp_Pnt2d aFirst, aLast;
440 curves[0]->D0( curves[0]->FirstParameter(), aFirst );
441 curves[n-1]->D0( curves[n-1]->LastParameter(), aLast );
442 Handle(Geom2d_Line) aV1 = new Geom2d_Line( aFirst, gp_Dir2d( 0, 1 ) );
443 Handle(Geom2d_TrimmedCurve) aT1 = new Geom2d_TrimmedCurve( aV1, 0.0, aVMax );
445 Handle(Geom2d_Line) aV2 = new Geom2d_Line( aLast, gp_Dir2d( 0, 1 ) );
446 Handle(Geom2d_TrimmedCurve) aT2 = new Geom2d_TrimmedCurve( aV2, 0.0, aVMax );
448 curves.push_back( aT1 );
449 curves.push_back( aT2 );
451 int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 );
452 theMidPointCurve = CurveUZ( theXCurv );
453 theMidPointCurve.reserve( psize );
454 theWidthCurve = CurveUZ( theXCurv );
455 theWidthCurve.reserve( psize );
458 // for each discrete value of z we search intersection with profile
459 for( double z = theMinZ; z <= theMaxZ; z += theDDZ )
461 Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z ), gp_Dir2d( 1, 0 ) );
462 std::vector<gp_Pnt2d> intersections;
463 for( size_t i = 0; i < n; i++ )
465 Handle_Geom2d_Curve aCurve = curves[i];
466 Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
467 for( int k=1, m=anIntersect.NbPoints(); k<=m; k++ )
468 intersections.push_back( anIntersect.Point( k ) );
471 if( intersections.size() >= 2 )
473 double u_mid = CalcGC( intersections );
477 theMidPointCurve.push_back( p_mid );
479 double u_wid = CalcWidth( intersections );
483 theWidthCurve.push_back( p_wid );
488 void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCurveB,
489 int theNbSteps, std::vector<CurveUZ>& theInterpolation,
492 theInterpolation.clear();
493 int d = isAddSecond ? 2 : 1;
494 theInterpolation.reserve( theNbSteps+d );
495 double dt = 1.0 / double( theNbSteps + 1 );
497 theInterpolation.push_back( theCurveA );
498 for( int i=0; i<theNbSteps; i++, t+=dt )
500 CurveUZ anInterp = theCurveA*(1-t) + theCurveB*t;
501 theInterpolation.push_back( anInterp );
504 theInterpolation.push_back( theCurveB );
507 void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
508 const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
509 AltitudePoints& thePoints,
510 Bank* theLeft, Bank* theRight, double dz )
512 Geom2dAdaptor_Curve anAdaptor( theHydraulicAxis );
513 GCPnts_AbscissaPoint ap( anAdaptor, theMidCurve.Xcurv(), anAdaptor.FirstParameter() );
514 double aParam = ap.Parameter();
517 gp_Vec2d tangent, profile_dir;
518 anAdaptor.D1( aParam, point, tangent );
519 profile_dir.SetCoord( tangent.Y(), -tangent.X() );
520 profile_dir.Normalize();
521 gp_Dir tangent_n( tangent.X(), tangent.Y(), dz );
523 size_t n = theMidCurve.size();
524 double min_param = 1E+15;
525 double max_param = -1E+15;
527 for( size_t i=0; i<n; i++ )
529 double param1 = theMidCurve[i].U - theWidthCurve[i].U / 2;
530 double param2 = theMidCurve[i].U + theWidthCurve[i].U / 2;
532 gp_Pnt2d p1 = point.Translated( param1 * profile_dir / 2 );
533 gp_Pnt2d p2 = point.Translated( param2 * profile_dir / 2 );
535 double z = theMidCurve[i].Z;
537 if( param1 < min_param )
542 if( param2 < min_param )
547 if( param1 > max_param )
552 if( param2 > max_param )
558 AltitudePoint p3d_1( p1.X(), p1.Y(), z ), p3d_2( p2.X(), p2.Y(), z );
559 thePoints.push_back( p3d_1 );
560 thePoints.push_back( p3d_2 );
565 gp_Pnt2d left2d = point.Translated( min_param * profile_dir / 2 );
566 gp_Pnt left( left2d.X(), left2d.Y(), z1 );
567 theLeft->push_back( left, tangent_n );
571 gp_Pnt2d right2d = point.Translated( max_param * profile_dir / 2 );
572 gp_Pnt right( right2d.X(), right2d.Y(), z2 );
573 theRight->push_back( right, tangent_n );
577 inline double max( double a, double b )
585 HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
586 ( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
587 const Handle_HYDROData_Profile& theProfileA,
589 const Handle_HYDROData_Profile& theProfileB,
591 double theDDZ, int theNbSteps, bool isAddSecond,
592 Bank* theLeft, Bank* theRight )
594 double zminA, zmaxA, zminB, zmaxB;
595 gp_Pnt lowestA, lowestB;
598 GetProperties( theProfileA, lowestA, dirA, false, zminA, zmaxA );
599 GetProperties( theProfileB, lowestB, dirB, false, zminB, zmaxB );
601 double dz = zminB - zminA;
603 double zmin = max( zminA, zminB );
604 double zmax = max( zmaxA, zmaxB );
606 CurveUZ midA(0), midB(0);
607 CurveUZ widA(0), widB(0);
609 ProfileDiscretization( theProfileA, theXCurvA, zmin, zmax, theDDZ, midA, widA );
610 ProfileDiscretization( theProfileB, theXCurvB, zmin, zmax, theDDZ, midB, widB );
612 std::vector<CurveUZ> mid, wid;
613 Interpolate( midA, midB, theNbSteps, mid, isAddSecond );
614 Interpolate( widA, widB, theNbSteps, wid, isAddSecond );
616 size_t p = mid.size();
617 size_t q = p>0 ? 2*mid[0].size() : 1;
618 AltitudePoints points;
619 points.reserve( p*q );
620 for( size_t i=0; i<p; i++ )
621 CurveTo3D( theHydraulicAxis, mid[i], wid[i], points, theLeft, theRight, dz );
626 HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
627 ( const std::vector<Handle_HYDROData_Profile>& theProfiles,
628 double theDDZ, double theSpatialStep,
632 AltitudePoints points;
633 size_t n = theProfiles.size();
637 std::vector<double> distances;
638 Handle_Geom2d_BSplineCurve aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
639 if( aHydraulicAxis.IsNull() )
642 int aNbStepsComplete = 0;
643 for( size_t i=0, n1=n-1; i<n1; i++ )
645 double aDistance = distances[i+1]-distances[i];
646 aNbStepsComplete += ( int(aDistance/theSpatialStep) + 1 );
649 theLeft->reserve( aNbStepsComplete );
651 theRight->reserve( aNbStepsComplete );
653 for( size_t i=0, n1=n-1; i<n1; i++ )
655 double aDistance = distances[i+1]-distances[i];
656 int aNbSteps = int(aDistance/theSpatialStep) + 1;
657 bool isAddSecond = i==n1-1;
659 AltitudePoints local_points = Interpolate( aHydraulicAxis, theProfiles[i], distances[i],
660 theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond, theLeft, theRight );
663 points.reserve( local_points.size() * ( n-1 ) );
665 for( size_t j=0, m=local_points.size(); j<m; j++ )
666 points.push_back( local_points[j] );
671 void HYDROData_DTM::CreateBankShapes( TopoDS_Edge& theLeft, TopoDS_Edge& theRight ) const
673 theLeft = myLeft.createEdge3d();
674 theRight = myRight.createEdge3d();