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, const gp_Vec2d& theProfileDir )
39 : myXcurv( theXCurv ), myProfileDir( theProfileDir )
43 HYDROData_DTM::CurveUZ::~CurveUZ()
47 double HYDROData_DTM::CurveUZ::Xcurv() const
52 gp_Vec2d HYDROData_DTM::CurveUZ::ProfileDir() const
57 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator + ( const CurveUZ& c ) const
59 HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv(), ProfileDir() + c.ProfileDir() );
62 for( int i=0; i<n; i++ )
65 p.U = operator[]( i ).U + c[i].U;
66 p.Z = operator[]( i ).Z;
72 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator * ( double d ) const
74 HYDROData_DTM::CurveUZ res( Xcurv()*d, ProfileDir()*d );
77 for( int i=0; i<n; i++ )
80 p.U = operator[]( i ).U * d;
81 p.Z = operator[]( i ).Z;
90 HYDROData_DTM::HYDROData_DTM()
94 HYDROData_DTM::~HYDROData_DTM()
98 HYDROData_SequenceOfObjects HYDROData_DTM::GetProfiles() const
100 return GetReferenceObjects( DataTag_Profiles );
103 void HYDROData_DTM::SetProfiles( const HYDROData_SequenceOfObjects& theProfiles )
105 SetReferenceObjects( theProfiles, DataTag_Profiles );
108 double HYDROData_DTM::GetDDZ() const
110 return GetDouble( DataTag_DDZ );
113 void HYDROData_DTM::SetDDZ( double theDDZ )
115 SetDouble( DataTag_DDZ, theDDZ );
118 double HYDROData_DTM::GetSpatialStep() const
120 return GetDouble( DataTag_SpatialStep );
123 void HYDROData_DTM::SetSpatialStep( double theSpatialStep )
125 SetDouble( DataTag_SpatialStep, theSpatialStep );
128 void HYDROData_DTM::Update()
130 HYDROData_SequenceOfObjects objs = GetProfiles();
131 int aLower = objs.Lower(), anUpper = objs.Upper();
132 size_t n = anUpper-aLower+1;
134 std::vector<Handle_HYDROData_Profile> profiles;
135 profiles.reserve( n );
136 for( int i=aLower; i<=anUpper; i++ )
138 Handle(HYDROData_Profile) aProfile = Handle(HYDROData_Profile)::DownCast( objs.Value( i ) );
139 if( !aProfile.IsNull() )
140 profiles.push_back( aProfile );
143 double ddz = GetDDZ();
144 double step = GetSpatialStep();
145 const double EPS = 1E-3;
146 AltitudePoints points;
148 AltitudePoints right;
149 std::vector<AltitudePoints> main_profiles;
151 if( ddz>EPS && step>EPS )
152 points = Interpolate( profiles, ddz, step, left, right, main_profiles );
154 SetAltitudePoints( points );
155 //SetTopShape(); //2d
168 void HYDROData_DTM::GetProperties( const Handle_HYDROData_Profile& theProfile,
169 gp_Pnt& theLowestPoint, gp_Vec2d& theDir,
171 double& theZMin, double& theZMax )
173 theLowestPoint = theProfile->GetBottomPoint();
176 theProfile->GetLeftPoint( aLeft, true, true );
177 theProfile->GetRightPoint( aRight, true, true );
178 double x = aRight.X()-aLeft.X();
179 double y = aRight.Y()-aLeft.Y();
181 theDir = gp_Vec2d( -y, x );
183 theDir = gp_Vec2d( x, y );
185 HYDROData_Profile::ProfilePoints points = theProfile->GetProfilePoints();
186 int lo = points.Lower();
187 int up = points.Upper();
188 theZMin = std::numeric_limits<double>::max();
190 for( int i=lo; i<=up; i++ )
192 double z = points.Value( i ).Z();
200 inline gp_Pnt2d To2D( const gp_Pnt& thePnt, const gp_Trsf& theTr,
201 double& theUMin, double& theUMax )
203 gp_Pnt p = thePnt.Transformed( theTr );
210 return gp_Pnt2d( u, z );
213 Handle(TColgp_HArray1OfPnt2d) To2D( const TColgp_Array1OfPnt& thePoints,
214 const gp_Trsf& theTr,
215 double& theUMin, double& theUMax )
217 int low = thePoints.Lower(), up = thePoints.Upper();
218 Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( low, up );
219 for( int i=low; i<=up; i++ )
220 points->SetValue( i, To2D( thePoints.Value( i ), theTr, theUMin, theUMax ) );
224 Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve,
225 Standard_Real theFirst, Standard_Real theLast,
226 const gp_Trsf& theTr,
227 double& theUMin, double& theUMax )
229 if( theCurve->IsKind( STANDARD_TYPE( Geom_Line ) ) )
231 gp_Pnt aFirstPnt, aLastPnt;
232 theCurve->D0( theFirst, aFirstPnt );
233 theCurve->D0( theLast, aLastPnt );
236 aFirst2d = To2D( aFirstPnt, theTr, theUMin, theUMax ),
237 aLast2d = To2D( aLastPnt, theTr, theUMin, theUMax );
239 gp_Vec2d dir( aFirst2d, aLast2d );
240 Handle_Geom2d_Line aLine2d = new Geom2d_Line( aFirst2d, gp_Dir2d( dir.X(), dir.Y() ) );
241 return new Geom2d_TrimmedCurve( aLine2d, 0, aLast2d.Distance( aFirst2d ) );
244 if( theCurve->IsKind( STANDARD_TYPE( Geom_BSplineCurve ) ) )
246 Handle(Geom_BSplineCurve) aSpline = Handle(Geom_BSplineCurve)::DownCast( theCurve );
248 Handle(TColgp_HArray1OfPnt2d) poles = To2D( aSpline->Poles(), theTr, theUMin, theUMax );
249 const TColStd_Array1OfReal& knots = aSpline->Knots();
250 const TColStd_Array1OfInteger& multiplicities = aSpline->Multiplicities();
251 int aDegree = aSpline->Degree();
253 return new Geom2d_BSplineCurve( poles->Array1(), knots, multiplicities, aDegree );
256 return Handle(Geom2d_Curve)();
259 Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis(
260 const std::vector<Handle_HYDROData_Profile>& theProfiles,
261 std::vector<double>& theDistances )
263 size_t n = theProfiles.size();
264 Handle_Geom2d_BSplineCurve aResult;
266 Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( 1, (int)n );
267 TColgp_Array1OfVec2d tangents( 1, (int)n );
268 Handle(TColStd_HArray1OfBoolean) flags = new TColStd_HArray1OfBoolean( 1, (int)n );
270 for( size_t i = 1; i <= n; i++ )
272 Handle_HYDROData_Profile aProfile = theProfiles[i-1];
278 GetProperties( aProfile, aLowest, aTangent, true, zmin, zmax );
279 aTangent.Normalize();
281 points->SetValue( (int)i, gp_Pnt2d( aLowest.X(), aLowest.Y() ) );
282 tangents.SetValue( (int)i, aTangent );
283 flags->SetValue( (int)i, Standard_True );
286 Geom2dAPI_Interpolate anInterpolator( points, Standard_False, Standard_False );
287 anInterpolator.Load( tangents, flags );
288 anInterpolator.Perform();
289 if( anInterpolator.IsDone() )
291 aResult = anInterpolator.Curve();
293 //fill the distances vector
294 Geom2dAdaptor_Curve anAdaptor( aResult );
296 theDistances.clear();
297 theDistances.reserve( n );
298 Standard_Real aParamFirst = anAdaptor.FirstParameter(), aParamLast = anAdaptor.LastParameter();
299 for( size_t i = 1; i <= n; i++ )
301 gp_Pnt2d aPnt = points->Value( (Standard_Integer)i );
302 Geom2dAPI_ProjectPointOnCurve aProject( aPnt, aResult );
303 Standard_Real aParam = aProject.LowerDistanceParameter();
304 double aDistance = GCPnts_AbscissaPoint::Length( anAdaptor, aParamFirst, aParam );
305 theDistances.push_back( aDistance );
311 std::vector<Handle_Geom2d_Curve> HYDROData_DTM::ProfileToParametric(
312 const Handle_HYDROData_Profile& theProfile,
313 double& theUMin, double& theUMax, gp_Vec2d& theDir )
315 std::vector<Handle_Geom2d_Curve> curves;
317 // Transformation of the coordinate systems
320 GetProperties( theProfile, aLowest, theDir, false, zmin, zmax );
322 gp_Ax3 aStd3d( gp_Pnt( 0, 0, 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( 1, 0, 0 ) );
323 gp_Ax3 aLocal( aLowest, gp_Dir( 0, 0, 1 ), gp_Dir( theDir.X(), theDir.Y(), 0 ) );
326 aTransf.SetTransformation( aStd3d, aLocal );
328 // Iteration via edges
329 TopoDS_Wire aWire = TopoDS::Wire( theProfile->GetShape3D() );
330 TopExp_Explorer anExp( aWire, TopAbs_EDGE );
331 for( ; anExp.More(); anExp.Next() )
333 // Extract an edge from wire
334 TopoDS_Edge anEdge = TopoDS::Edge( anExp.Current() );
336 // Extract a curve corresponding to the edge
337 TopLoc_Location aLoc;
338 Standard_Real aFirst, aLast;
339 Handle(Geom_Curve) aCurve = BRep_Tool::Curve( anEdge, aLoc, aFirst, aLast );
341 // Convert the curve to 2d CS
342 Handle(Geom2d_Curve) aCurve2d = CurveTo2D( aCurve, aFirst, aLast, aTransf, theUMin, theUMax );
343 if( !aCurve2d.IsNull() )
344 curves.push_back( aCurve2d );
350 bool CalcMidWidth( const std::vector<gp_Pnt2d>& intersections, double& theMid, double& theWid )
352 double umin = std::numeric_limits<double>::max(),
355 size_t n = intersections.size();
359 for( size_t i = 0; i < n; i++ )
361 double u = intersections[i].X();
367 theMid = ( umin+umax )/2;
372 void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& theProfile,
373 double theXCurv, double theMinZ, double theMaxZ, double theDDZ,
374 CurveUZ& theMidPointCurve,
375 CurveUZ& theWidthCurve,
376 double theTolerance )
378 double aDblMax = std::numeric_limits<double>::max(),
383 gp_Vec2d aProfileDir;
384 std::vector<Handle_Geom2d_Curve> curves = ProfileToParametric( theProfile, aUMin, aUMax, aProfileDir );
385 size_t n = curves.size();
390 // we add the "virtual" vertical lines to simulate the intersection with profile
391 gp_Pnt2d aFirst, aLast;
392 curves[0]->D0( curves[0]->FirstParameter(), aFirst );
393 curves[n-1]->D0( curves[n-1]->LastParameter(), aLast );
394 Handle(Geom2d_Line) aV1 = new Geom2d_Line( aFirst, gp_Dir2d( 0, 1 ) );
395 Handle(Geom2d_TrimmedCurve) aT1 = new Geom2d_TrimmedCurve( aV1, 0.0, aVMax );
397 Handle(Geom2d_Line) aV2 = new Geom2d_Line( aLast, gp_Dir2d( 0, 1 ) );
398 Handle(Geom2d_TrimmedCurve) aT2 = new Geom2d_TrimmedCurve( aV2, 0.0, aVMax );
400 curves.push_back( aT1 );
401 curves.push_back( aT2 );
403 int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 );
404 theMidPointCurve = CurveUZ( theXCurv, aProfileDir );
405 theMidPointCurve.reserve( psize );
406 theWidthCurve = CurveUZ( theXCurv, aProfileDir );
407 theWidthCurve.reserve( psize );
410 // for each discrete value of z we search intersection with profile
411 for( double z = theMinZ; z <= theMaxZ; z += theDDZ )
413 Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z ), gp_Dir2d( 1, 0 ) );
414 std::vector<gp_Pnt2d> intersections;
415 for( size_t i = 0; i < n; i++ )
417 Handle_Geom2d_Curve aCurve = curves[i];
418 Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
419 for( int k=1, m=anIntersect.NbPoints(); k<=m; k++ )
420 intersections.push_back( anIntersect.Point( k ) );
423 if( intersections.size() >= 2 )
426 if( !CalcMidWidth( intersections, u_mid, u_wid ) )
432 theMidPointCurve.push_back( p_mid );
437 theWidthCurve.push_back( p_wid );
442 void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCurveB,
443 int theNbSteps, std::vector<CurveUZ>& theInterpolation,
446 theInterpolation.clear();
447 int d = isAddSecond ? 2 : 1;
448 theInterpolation.reserve( theNbSteps+d );
449 double dt = 1.0 / double( theNbSteps + 1 );
451 theInterpolation.push_back( theCurveA );
452 for( int i=0; i<theNbSteps; i++, t+=dt )
454 CurveUZ anInterp = theCurveA*(1-t) + theCurveB*t;
455 theInterpolation.push_back( anInterp );
458 theInterpolation.push_back( theCurveB );
460 #include <BRepLib_MakeEdge2d.hxx>
461 void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
462 const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
463 AltitudePoints& thePoints, double dz )
465 Geom2dAdaptor_Curve anAdaptor( theHydraulicAxis );
466 TopoDS_Edge E2d = BRepLib_MakeEdge2d(theHydraulicAxis).Edge();
467 GCPnts_AbscissaPoint ap( anAdaptor, theMidCurve.Xcurv(), anAdaptor.FirstParameter() );
468 double aParam = ap.Parameter();
471 anAdaptor.D0( aParam, point );
472 gp_Vec2d profile_dir = theMidCurve.ProfileDir();
473 gp_Dir tangent_n( -profile_dir.Y(), profile_dir.X(), dz );
474 profile_dir.Normalize();
476 size_t n = theMidCurve.size();
477 std::map<double, AltitudePoint> sorted_points;
478 for( size_t i=0; i<n; i++ )
480 double param1 = theMidCurve[i].U - theWidthCurve[i].U / 2;
481 double param2 = theMidCurve[i].U + theWidthCurve[i].U / 2;
483 gp_Pnt2d p1 = point.Translated( param1 * profile_dir);
484 gp_Pnt2d p2 = point.Translated( param2 * profile_dir);
486 double z = theMidCurve[i].Z;
488 AltitudePoint p3d_1( p1.X(), p1.Y(), z ), p3d_2( p2.X(), p2.Y(), z );
490 sorted_points[param1] = p3d_1;
491 sorted_points[param2] = p3d_2;
494 thePoints.reserve( sorted_points.size() );
495 std::map<double, AltitudePoint>::const_iterator it = sorted_points.begin(), last = sorted_points.end();
496 for( ; it!=last; it++ )
497 thePoints.push_back( it->second );
500 inline double max( double a, double b )
508 #include <BRepLib_MakeWire.hxx>
510 std::vector<HYDROData_Bathymetry::AltitudePoints> HYDROData_DTM::Interpolate
511 ( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
512 const Handle_HYDROData_Profile& theProfileA,
514 const Handle_HYDROData_Profile& theProfileB,
516 double theDDZ, int theNbSteps, bool isAddSecond )
518 double zminA, zmaxA, zminB, zmaxB;
519 gp_Pnt lowestA, lowestB;
522 GetProperties( theProfileA, lowestA, dirA, false, zminA, zmaxA );
523 GetProperties( theProfileB, lowestB, dirB, false, zminB, zmaxB );
525 double dz = zminB - zminA;
527 double zmin = max( zminA, zminB );
528 double zmax = max( zmaxA, zmaxB );
530 CurveUZ midA(0, gp_Vec2d()), midB(0, gp_Vec2d());
531 CurveUZ widA(0, gp_Vec2d()), widB(0, gp_Vec2d());
533 ProfileDiscretization( theProfileA, theXCurvA, zmin, zmax, theDDZ, midA, widA );
534 ProfileDiscretization( theProfileB, theXCurvB, zmin, zmax, theDDZ, midB, widB );
536 std::vector<CurveUZ> mid, wid;
537 Interpolate( midA, midB, theNbSteps, mid, isAddSecond );
538 Interpolate( widA, widB, theNbSteps, wid, isAddSecond );
540 size_t p = mid.size();
541 size_t q = p>0 ? 2*mid[0].size() : 1;
542 std::vector<AltitudePoints> points;
545 for( size_t i=0; i<p; i++ )
547 points[i].reserve( q );
548 CurveTo3D( theHydraulicAxis, mid[i], wid[i], points[i], dz );
554 HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
555 ( const std::vector<Handle_HYDROData_Profile>& theProfiles,
556 double theDDZ, double theSpatialStep,
557 AltitudePoints& theLeft,
558 AltitudePoints& theRight,
559 std::vector<AltitudePoints>& theMainProfiles )
561 AltitudePoints points;
562 size_t n = theProfiles.size();
566 std::vector<double> distances;
567 Handle_Geom2d_BSplineCurve aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
568 if( aHydraulicAxis.IsNull() )
571 theMainProfiles.reserve( n );
573 for( size_t i=0, n1=n-1; i<n1; i++ )
575 double aDistance = distances[i+1]-distances[i];
576 int aNbSteps = int(aDistance/theSpatialStep);
577 bool isAddSecond = i==n1-1;
579 // 1. Calculate interpolated profiles
580 std::vector<AltitudePoints> local_points = Interpolate( aHydraulicAxis, theProfiles[i], distances[i],
581 theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond );
582 int lps = local_points.size();
584 // 2. Put all points into the global container
585 for( size_t j=0; j<lps; j++ )
587 const AltitudePoints& lp = local_points[j];
589 points.reserve( lp.size() * n );
590 for( size_t k=0, ks=lp.size(); k<ks; k++ )
591 points.push_back( lp[k] );
594 // 3. Get left/right banks' points
597 theLeft.reserve( lps * n );
598 theRight.reserve( lps * n );
600 for( size_t j=0; j<lps; j++ )
602 const AltitudePoints& lp = local_points[j];
603 theLeft.push_back( lp[0] );
604 theRight.push_back( lp[lp.size()-1] );
607 // 4. Get main profiles points
608 theMainProfiles.push_back( local_points[0] );
610 theMainProfiles.push_back( local_points[lps-1] );