X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FHYDROData%2FHYDROData_DTM.cxx;h=207530f8df59e3a0c37bc80081e3315a3197bbcc;hb=bd81def133c235351f4c18834c6660450f16b9ee;hp=9ac801dcfb81fb2c8eb5b5aea323e07d5808f646;hpb=4bb18a8f971e38d25612f346158d7c0ec99fb47a;p=modules%2Fhydro.git diff --git a/src/HYDROData/HYDROData_DTM.cxx b/src/HYDROData/HYDROData_DTM.cxx index 9ac801dc..207530f8 100644 --- a/src/HYDROData/HYDROData_DTM.cxx +++ b/src/HYDROData/HYDROData_DTM.cxx @@ -36,25 +36,34 @@ #include #include #include -#include #include #include #include #include #include +#include #include #include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include - - -IMPLEMENT_STANDARD_HANDLE( HYDROData_DTM, HYDROData_Bathymetry ) IMPLEMENT_STANDARD_RTTIEXT( HYDROData_DTM, HYDROData_Bathymetry ) - -HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv, const gp_Vec2d& theProfileDir ) - : myXcurv( theXCurv ), myProfileDir( theProfileDir ) +HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv, const gp_Vec2d& theProfileDir, double theDeltaZ, double theMaxZ ) + : myXcurv( theXCurv ), myProfileDir( theProfileDir ), myDeltaZ( theDeltaZ ), myMaxZ (theMaxZ) { } @@ -72,12 +81,27 @@ gp_Vec2d HYDROData_DTM::CurveUZ::ProfileDir() const return myProfileDir; } +double HYDROData_DTM::CurveUZ::DeltaZ() const +{ + return myDeltaZ; +} + +double HYDROData_DTM::CurveUZ::MaxZ() const +{ + return myMaxZ; +} + HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator + ( const CurveUZ& c ) const { - HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv(), ProfileDir() + c.ProfileDir() ); - size_t n = size(); - res.reserve( n ); - for( int i=0; iSetValue(i+1, gp_Pnt(pnts[i].X, pnts[i].Y, pnts[i].Z)); + + GeomAPI_Interpolate anInterpolator(gpPoints, Standard_False,1.0e-6); + anInterpolator.Perform() ; + if (anInterpolator.IsDone()) { - gp_Pnt p1(pnts[i].X, pnts[i].Y, pnts[i].Z); - gp_Pnt p2(pnts[i+1].X, pnts[i+1].Y, pnts[i+1].Z); - WM.Add(BRepLib_MakeEdge(p1, p2).Edge()); + Handle(Geom_BSplineCurve) C = anInterpolator.Curve(); + E = BRepBuilderAPI_MakeEdge(C).Edge(); } - if (WM.IsDone()) - W = WM.Wire(); } TopTools_IndexedMapOfOrientedShape HYDROData_DTM::Create3DShape(const AltitudePoints& left, @@ -163,24 +201,27 @@ TopTools_IndexedMapOfOrientedShape HYDROData_DTM::Create3DShape(const AltitudePo const std::vector& main_profiles) { TopTools_IndexedMapOfOrientedShape ll; - TopoDS_Wire LWire, RWire; - PointToWire(left, LWire); - PointToWire(right, RWire); - if (!LWire.IsNull()) - ll.Add(LWire.Oriented(TopAbs_FORWARD)); + //TopoDS_Wire LWire, RWire; + //PointsToWire(left, LWire); + //PointsToWire(right, RWire); + TopoDS_Edge LEdge, REdge; + PointsToEdge(left, LEdge); + PointsToEdge(right, REdge); + if (!LEdge.IsNull()) + ll.Add(LEdge.Oriented(TopAbs_FORWARD)); for (int k = 0; k < main_profiles.size(); k++) { TopoDS_Wire W; - PointToWire(main_profiles[k], W); + PointsToWire(main_profiles[k], W); TopAbs_Orientation Ori = TopAbs_INTERNAL; if (k == 0 || k == main_profiles.size() - 1) Ori = TopAbs_FORWARD; ll.Add(W.Oriented(Ori)); } - if (!RWire.IsNull()) - ll.Add(RWire.Oriented(TopAbs_FORWARD)); + if (!REdge.IsNull()) + ll.Add(REdge.Oriented(TopAbs_FORWARD)); //yes, add subshapes in this order (left + profiles + right) //otherwise the projected wire will be non-manifold @@ -201,9 +242,11 @@ void HYDROData_DTM::Update() HYDROData_SequenceOfObjects objs = GetProfiles(); double ddz = GetDDZ(); double step = GetSpatialStep(); - CreateProfilesFromDTM( objs, ddz, step, points, Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, true, true, std::set() ); + std::set InvInd; + bool WireIntersections; //__TODO + CreateProfilesFromDTM( objs, ddz, step, points, Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, true, true, InvInd, -1, WireIntersections ); SetAltitudePoints( points ); - + SetShape( DataTag_LeftBankShape, OutLeftB); SetShape( DataTag_RightBankShape, OutRightB); SetShape( DataTag_InletShape, OutInlet); @@ -229,7 +272,6 @@ void HYDROData_DTM::GetPresentationShapes( TopoDS_Shape& Out3dPres, Out3dPres = GetShape( DataTag_3DShape ); Out2dPres = GetShape( DataTag_2DShape ); } - void HYDROData_DTM::CreateProfilesFromDTM (const HYDROData_SequenceOfObjects& InpProfiles, double ddz, double step, @@ -242,13 +284,14 @@ void HYDROData_DTM::CreateProfilesFromDTM (const HYDROData_SequenceOfObjects& In TopoDS_Shape& OutOutlet, bool Create3dPres, bool Create2dPres, - std::set& InvInd) + std::set& InvInd, + int thePntsLimit, + bool& WireIntersections) { - int aLower = InpProfiles.Lower(), anUpper = InpProfiles.Upper(); size_t n = anUpper - aLower + 1; - std::vector profiles; + std::vector profiles; profiles.reserve( n ); for( int i=aLower; i<=anUpper; i++ ) { @@ -261,72 +304,59 @@ void HYDROData_DTM::CreateProfilesFromDTM (const HYDROData_SequenceOfObjects& In AltitudePoints right; std::vector main_profiles; + if( thePntsLimit > 0 ) + { + int aNbPoints = EstimateNbPoints( profiles, ddz, step ); + if( aNbPoints < 0 || aNbPoints > thePntsLimit ) + return; + } + if( ddz>EPS && step>EPS ) CreateProfiles(profiles, ddz, step, left, right, points, main_profiles, - Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, Create3dPres, Create2dPres, InvInd ); + Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, Create3dPres, Create2dPres, InvInd, WireIntersections ); } -void HYDROData_DTM::ProjWireOnPlane(const TopoDS_Wire& inpWire, const Handle_Geom_Plane& RefPlane, TopoDS_Wire& outWire) +bool HYDROData_DTM::GetPlanarFaceFromBanks( const TopoDS_Edge& LB, const TopoDS_Edge& RB, TopoDS_Face& outF, + TopTools_SequenceOfShape* Boundr) { - //its also possible to use BrepAlgo_NormalProjection here! - BRepTools_WireExplorer ex(TopoDS::Wire(inpWire.Oriented(TopAbs_FORWARD))); - BRepLib_MakeWire WM; - for (;ex.More();ex.Next()) - { - const TopoDS_Edge& CE = ex.Current(); - double f, l; - Handle(Geom_Curve) C3d = BRep_Tool::Curve(CE, f, l); - Handle(Geom_Curve) ProjectedCurve = GeomProjLib::ProjectOnPlane(new Geom_TrimmedCurve(C3d, f, l), RefPlane, RefPlane->Position().Direction(), Standard_True); - TopoDS_Edge ProjEdge = BRepLib_MakeEdge(ProjectedCurve); - WM.Add(ProjEdge); //auto sharing between edges if vertex is coincident - } - outWire = WM.Wire(); - outWire.Orientation(inpWire.Orientation()); //take from the original wire -} + BRep_Builder BB; + TopoDS_Face F; + Handle_Geom_Plane refpl = new Geom_Plane(gp_Pnt(0,0,0), gp_Dir(0,0,1)); + TopoDS_Vertex VFI, VLI, VFO, VLO; + TopoDS_Edge prLB; + TopoDS_Edge prRB; -void HYDROData_DTM::Get2dFaceFrom3dPres(const TopoDS_Compound& cmp, TopoDS_Face& outF ) -{ - Handle_Geom_Plane refpl = new Geom_Plane(gp_Pnt(0,0,0), gp_Dir(0,0,1)); - BRepLib_MakeWire WM; - TopoDS_Iterator it(cmp); - //TopTools_IndexedMapOfShape IntW; - for (;it.More(); it.Next()) - { - const TopoDS_Wire& W = TopoDS::Wire(it.Value()); - if (W.Orientation() != TopAbs_INTERNAL) - { - //use list of edges to protect againts non-manifold cases. - //auto sharing between edges will be added automatically - TopTools_IndexedMapOfShape ME; - TopTools_ListOfShape LE; - TopExp::MapShapes(W, TopAbs_EDGE, ME); - for (int i = 1; i <= ME.Extent(); i++) - LE.Append(ME(i)); - WM.Add(LE); - } - //else - // IntW.Add(W); - } - TopoDS_Wire outW; - ProjWireOnPlane(WM.Wire(), refpl, outW); - BRepBuilderAPI_MakeFace mf(refpl, outW); //check inside is true by def - outF = mf.Face(); + BRepAdaptor_Curve LBAD(LB); + Handle_Geom_Curve LBPC = GeomProjLib::ProjectOnPlane(LBAD.Curve().Curve(), refpl, gp_Dir(0, 0, -1), 1 ); + prLB = BRepLib_MakeEdge(LBPC).Edge(); + + BRepAdaptor_Curve RBAD(RB); + Handle_Geom_Curve RBPC = GeomProjLib::ProjectOnPlane(RBAD.Curve().Curve(), refpl, gp_Dir(0, 0, -1), 1 ); + prRB = BRepLib_MakeEdge(RBPC).Edge(); - ///!!! the internal wires cant be added with 'internal' ori. - // it's possible to do with brep builder yet the result will not be correct! - // more proper way is to use BOP operation here. - /*for (int i = 1; i <= IntW.Extent(); i++) + TopExp::Vertices(prLB, VFI, VFO, 1); + TopExp::Vertices(prRB, VLI, VLO, 1); + TopoDS_Edge prIL = BRepLib_MakeEdge(VFI, VLI).Edge(); + TopoDS_Edge prOL = BRepLib_MakeEdge(VFO, VLO).Edge(); + TopoDS_Wire prW = BRepLib_MakeWire(prLB, prIL, prOL, prRB).Wire(); + outF = BRepBuilderAPI_MakeFace(refpl->Pln(), prW, 1).Face(); + + if (Boundr) { - TopoDS_Wire outIW; - const TopoDS_Wire& W = TopoDS::Wire(IntW(i)); - ProjWireOnPlane(W, refpl, outIW); - BB.Add(outF, outIW); - }*/ + Boundr->Append(prLB); + Boundr->Append(prIL); + Boundr->Append(prOL); + Boundr->Append(prRB); + } + + ShapeAnalysis_Wire WA(prW, outF, Precision::Confusion()); + bool res = WA.CheckSelfIntersection(); + return !res; } -void HYDROData_DTM::CreateProfiles(const std::vector& theProfiles, +void HYDROData_DTM::CreateProfiles(const std::vector& theProfiles, double theDDZ, double theSpatialStep, AltitudePoints& theOutLeft, @@ -341,12 +371,16 @@ void HYDROData_DTM::CreateProfiles(const std::vector& TopoDS_Shape& OutOutlet, bool Create3dPres, bool Create2dPres, - std::set& InvInd) + std::set& InvInd, + bool& ProjStat) { if (theProfiles.empty()) return; theOutPoints = Interpolate( theProfiles, theDDZ, theSpatialStep, theOutLeft, theOutRight, theOutMainProfiles, InvInd ); //note that if Create3dPres is false => Create2dPres flag is meaningless! + if( theOutPoints.empty() ) + return; + if (Create3dPres) { TopTools_IndexedMapOfOrientedShape ll = Create3DShape( theOutLeft, theOutRight, theOutMainProfiles); @@ -370,7 +404,7 @@ void HYDROData_DTM::CreateProfiles(const std::vector& if (Create2dPres) { TopoDS_Face outF; - Get2dFaceFrom3dPres(cmp, outF); + ProjStat = GetPlanarFaceFromBanks(TopoDS::Edge(OutLeftB), TopoDS::Edge(OutRightB), outF, NULL); Out2dPres = outF; }; } @@ -379,9 +413,8 @@ void HYDROData_DTM::CreateProfiles(const std::vector& -void HYDROData_DTM::GetProperties( const Handle_HYDROData_Profile& theProfile, +void HYDROData_DTM::GetProperties( const Handle(HYDROData_Profile)& theProfile, gp_Pnt& theLowestPoint, gp_Vec2d& theDir, - bool isNormalDir, double& theZMin, double& theZMax ) { theLowestPoint = theProfile->GetBottomPoint(); @@ -391,15 +424,12 @@ void HYDROData_DTM::GetProperties( const Handle_HYDROData_Profile& theProfile, theProfile->GetRightPoint( aRight, true, true ); double x = aRight.X()-aLeft.X(); double y = aRight.Y()-aLeft.Y(); - if( isNormalDir ) - theDir = gp_Vec2d( -y, x ); - else - theDir = gp_Vec2d( x, y ); + theDir = gp_Vec2d( x, y ); HYDROData_Profile::ProfilePoints points = theProfile->GetProfilePoints(); int lo = points.Lower(); int up = points.Upper(); - theZMin = std::numeric_limits::max(); + theZMin = DBL_MAX; theZMax = -theZMin; for( int i=lo; i<=up; i++ ) { @@ -451,7 +481,7 @@ Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve, aLast2d = To2D( aLastPnt, theTr, theUMin, theUMax ); gp_Vec2d dir( aFirst2d, aLast2d ); - Handle_Geom2d_Line aLine2d = new Geom2d_Line( aFirst2d, gp_Dir2d( dir.X(), dir.Y() ) ); + Handle(Geom2d_Line) aLine2d = new Geom2d_Line( aFirst2d, gp_Dir2d( dir.X(), dir.Y() ) ); return new Geom2d_TrimmedCurve( aLine2d, 0, aLast2d.Distance( aFirst2d ) ); } @@ -470,31 +500,95 @@ Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve, return Handle(Geom2d_Curve)(); } -Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis( - const std::vector& theProfiles, +#include +#include +bool IsCooriented( const Handle(HYDROData_Profile)& theProfile1, + const Handle(HYDROData_Profile)& theProfile2 ) +{ + if( theProfile1==theProfile2 ) + return true; + + gp_XY lp1, rp1, lp2, rp2; + theProfile1->GetLeftPoint(lp1); + theProfile1->GetRightPoint(rp1); + theProfile2->GetLeftPoint(lp2); + theProfile2->GetRightPoint(rp2); + + GCE2d_MakeSegment s1(lp1, lp2); + GCE2d_MakeSegment s2(rp1, rp2); + + Geom2dAPI_InterCurveCurve inter; + inter.Init(s1.Value(), s2.Value()); + if (inter.NbPoints() == 0) + return true; + else + return false; +} + +Handle(Geom2d_BSplineCurve) HYDROData_DTM::CreateHydraulicAxis( + const std::vector& theProfiles, std::vector& theDistances ) { size_t n = theProfiles.size(); + if( n==1 ) + return Handle_Geom2d_BSplineCurve(); + Handle_Geom2d_BSplineCurve aResult; Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( 1, (int)n ); TColgp_Array1OfVec2d tangents( 1, (int)n ); Handle(TColStd_HArray1OfBoolean) flags = new TColStd_HArray1OfBoolean( 1, (int)n ); - for( size_t i = 1; i <= n; i++ ) + // Stage 1. Orient all profiles to be co-oriented with the first profile + theProfiles[0]->Update(); + for( size_t i = 1; i < n; i++ ) { - Handle_HYDROData_Profile aProfile = theProfiles[i-1]; + Handle(HYDROData_Profile) aProfile = theProfiles[i]; + Handle(HYDROData_Profile) aPrevProfile = theProfiles[i-1]; + + if( !IsCooriented( aProfile, aPrevProfile ) ) + { + gp_XY lp, rp; + aProfile->GetLeftPoint( lp, true ); + aProfile->GetRightPoint( rp, true ); + aProfile->SetLeftPoint( rp, true ); + aProfile->SetRightPoint( lp, true ); + } aProfile->Update(); + } + + // Stage 2. Calculate normals so that each normal "points" to the next profile + for( size_t i = 0; i < n; i++ ) + { + Handle(HYDROData_Profile) aProfile = theProfiles[i]; + Handle(HYDROData_Profile) aNextProfile = i==n-1 ? theProfiles[i-1] : theProfiles[i+1]; gp_Pnt aLowest; - gp_Vec2d aTangent; + gp_Vec2d aNormal; double zmin, zmax; - GetProperties( aProfile, aLowest, aTangent, true, zmin, zmax ); - aTangent.Normalize(); - points->SetValue( (int)i, gp_Pnt2d( aLowest.X(), aLowest.Y() ) ); - tangents.SetValue( (int)i, aTangent ); - flags->SetValue( (int)i, Standard_True ); + gp_XYZ curP = aProfile->GetBottomPoint(true); + gp_XY curP2d = gp_XY(curP.X(), curP.Y()); + + gp_XYZ nextP = aNextProfile->GetBottomPoint(true); + gp_XY nextP2d = gp_XY(nextP.X(), nextP.Y()); + + gp_Vec2d aPrTangent; + GetProperties( aProfile, aLowest, aPrTangent, zmin, zmax ); + aNormal.SetCoord( -aPrTangent.Y(), aPrTangent.X() ); + + gp_Vec2d aDirToNextProfile(nextP2d.X() - curP2d.X(), nextP2d.Y() - curP2d.Y() ); + if( i==n-1 ) + aDirToNextProfile.Reverse(); + + if (aNormal.Dot(aDirToNextProfile) < 0) + aNormal.Reverse(); + + aNormal.Normalize(); + + points->SetValue( (int)(i+1), gp_Pnt2d( aLowest.X(), aLowest.Y() ) ); + tangents.SetValue( (int)(i+1), aNormal ); + flags->SetValue( (int)(i+1), Standard_True ); } Geom2dAPI_Interpolate anInterpolator( points, Standard_False, Standard_False ); @@ -522,19 +616,19 @@ Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis( return aResult; } -std::vector HYDROData_DTM::ProfileToParametric( - const Handle_HYDROData_Profile& theProfile, +std::vector HYDROData_DTM::ProfileToParametric( + const Handle(HYDROData_Profile)& theProfile, double& theUMin, double& theUMax, gp_Vec2d& theDir ) { - std::vector curves; + std::vector curves; // Transformation of the coordinate systems gp_Pnt aLowest; double zmin, zmax; - GetProperties( theProfile, aLowest, theDir, false, zmin, zmax ); + GetProperties( theProfile, aLowest, theDir, zmin, zmax ); gp_Ax3 aStd3d( gp_Pnt( 0, 0, 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( 1, 0, 0 ) ); - gp_Ax3 aLocal( aLowest, gp_Dir( 0, 0, 1 ), gp_Dir( theDir.X(), theDir.Y(), 0 ) ); + gp_Ax3 aLocal( gp_Pnt( aLowest.X(), aLowest.Y(), 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( theDir.X(), theDir.Y(), 0 ) ); gp_Trsf aTransf; aTransf.SetTransformation( aStd3d, aLocal ); @@ -561,18 +655,19 @@ std::vector HYDROData_DTM::ProfileToParametric( } -bool CalcMidWidth( const std::vector& intersections, double& theMid, double& theWid ) +bool CalcMidWidth( const std::set& intersections, double& theMid, double& theWid ) { - double umin = std::numeric_limits::max(), + double umin = DBL_MAX, umax = -umin; size_t n = intersections.size(); - if( n <= 1 ) + if( n <= 0 ) return false; - for( size_t i = 0; i < n; i++ ) + std::set::const_iterator it = intersections.begin(), last = intersections.end(); + for( ; it!=last; it++ ) { - double u = intersections[i].X(); + double u = *it; if( uumax ) @@ -583,20 +678,20 @@ bool CalcMidWidth( const std::vector& intersections, double& theMid, d return true; } -void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& theProfile, - double theXCurv, double theMinZ, double theMaxZ, double theDDZ, +void HYDROData_DTM::ProfileDiscretization( const Handle(HYDROData_Profile)& theProfile, + double theXCurv, double theMinZ, double theMaxZ, double theTopZ, double theDDZ, CurveUZ& theMidPointCurve, CurveUZ& theWidthCurve, int& intersection_nb, double theTolerance) { - double aDblMax = std::numeric_limits::max(), + double aDblMax = DBL_MAX, aUMin = aDblMax, aUMax = -aUMin, aVMax = 1000000; gp_Vec2d aProfileDir; - std::vector curves = ProfileToParametric( theProfile, aUMin, aUMax, aProfileDir ); + std::vector curves = ProfileToParametric( theProfile, aUMin, aUMax, aProfileDir ); size_t n = curves.size(); if( n==0 ) @@ -616,32 +711,33 @@ void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& thePr curves.push_back( aT2 ); int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 ); - theMidPointCurve = CurveUZ( theXCurv, aProfileDir ); + theMidPointCurve = CurveUZ( theXCurv, aProfileDir, theMinZ, theTopZ); theMidPointCurve.reserve( psize ); - theWidthCurve = CurveUZ( theXCurv, aProfileDir ); + theWidthCurve = CurveUZ( theXCurv, aProfileDir, theMinZ, theTopZ ); theWidthCurve.reserve( psize ); n = curves.size(); // for each discrete value of z we search intersection with profile - for( double z = theMinZ; z <= theMaxZ; z += theDDZ ) + for( double z1 = theMinZ; z1 <= theMaxZ; z1 += theDDZ ) { - Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z ), gp_Dir2d( 1, 0 ) ); - std::vector intersections; + Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z1 ), gp_Dir2d( 1, 0 ) ); + std::set intersections; for( size_t i = 0; i < n; i++ ) { - Handle_Geom2d_Curve aCurve = curves[i]; + Handle(Geom2d_Curve) aCurve = curves[i]; Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance ); for( int k=1, m=anIntersect.NbPoints(); k<=m; k++ ) - intersections.push_back( anIntersect.Point( k ) ); + intersections.insert( anIntersect.Point( k ).X() ); } intersection_nb = intersections.size(); - if( intersection_nb >= 2 ) + if( intersection_nb >= 1 ) { double u_mid, u_wid; if( !CalcMidWidth( intersections, u_mid, u_wid ) ) continue; + double z = z1 - theMinZ; PointUZ p_mid; p_mid.U = u_mid; p_mid.Z = z; @@ -674,9 +770,9 @@ void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCur theInterpolation.push_back( theCurveB ); } #include -void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxis, +void HYDROData_DTM::CurveTo3D( const Handle(Geom2d_BSplineCurve)& theHydraulicAxis, const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve, - AltitudePoints& thePoints, double dz ) + AltitudePoints& thePoints ) { Geom2dAdaptor_Curve anAdaptor( theHydraulicAxis ); TopoDS_Edge E2d = BRepLib_MakeEdge2d(theHydraulicAxis).Edge(); @@ -686,12 +782,13 @@ void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxi gp_Pnt2d point; anAdaptor.D0( aParam, point ); gp_Vec2d profile_dir = theMidCurve.ProfileDir(); - gp_Dir tangent_n( -profile_dir.Y(), profile_dir.X(), dz ); + //gp_Dir tangent_n( -profile_dir.Y(), profile_dir.X(), dz ); profile_dir.Normalize(); size_t n = theMidCurve.size(); std::map sorted_points; - for( size_t i=0; i::const_iterator it = sorted_points.begin(), last = sorted_points.end(); for( ; it!=last; it++ ) - thePoints.push_back( it->second ); -} - -inline double max( double a, double b ) -{ - if( a>b ) - return a; - else - return b; -} + if( thePoints.empty() || thePoints.back().SquareDistance( it->second ) > EPS ) + thePoints.push_back( it->second ); +} + +//inline double max( double a, double b ) +//{ +// if( a>b ) +// return a; +// else +// return b; +//} +// +//inline double min( double a, double b ) +//{ +// if( a std::vector HYDROData_DTM::Interpolate - ( const Handle_Geom2d_BSplineCurve& theHydraulicAxis, - const Handle_HYDROData_Profile& theProfileA, + ( const Handle(Geom2d_BSplineCurve)& theHydraulicAxis, + const Handle(HYDROData_Profile)& theProfileA, double theXCurvA, - const Handle_HYDROData_Profile& theProfileB, + const Handle(HYDROData_Profile)& theProfileB, double theXCurvB, double theDDZ, int theNbSteps, bool isAddSecond, int& inter_nb_1, int& inter_nb_2) @@ -736,19 +855,21 @@ std::vector HYDROData_DTM::Interpolate gp_Pnt lowestA, lowestB; gp_Vec2d dirA, dirB; - GetProperties( theProfileA, lowestA, dirA, false, zminA, zmaxA ); - GetProperties( theProfileB, lowestB, dirB, false, zminB, zmaxB ); + GetProperties( theProfileA, lowestA, dirA, zminA, zmaxA ); + GetProperties( theProfileB, lowestB, dirB, zminB, zmaxB ); - double dz = zminB - zminA; + + double hmax = zmaxA-zminA > zmaxB-zminB ? zmaxA-zminA : zmaxB-zminB; - double zmin = max( zminA, zminB ); - double zmax = max( zmaxA, zmaxB ); + //double dz = zminB - zminA; + //double zmin = min( zminA, zminB ); + //double zmax = max( zmaxA, zmaxB ); - CurveUZ midA(0, gp_Vec2d()), midB(0, gp_Vec2d()); - CurveUZ widA(0, gp_Vec2d()), widB(0, gp_Vec2d()); + CurveUZ midA(0, gp_Vec2d(), 0, 0), midB(0, gp_Vec2d(), 0, 0); + CurveUZ widA(0, gp_Vec2d(), 0, 0), widB(0, gp_Vec2d(), 0, 0); - ProfileDiscretization( theProfileA, theXCurvA, zmin, zmax, theDDZ, midA, widA, inter_nb_1 ); - ProfileDiscretization( theProfileB, theXCurvB, zmin, zmax, theDDZ, midB, widB, inter_nb_2 ); + ProfileDiscretization( theProfileA, theXCurvA, zminA, zminA+hmax, zmaxA-zminA, theDDZ, midA, widA, inter_nb_1 ); + ProfileDiscretization( theProfileB, theXCurvB, zminB, zminB+hmax, zmaxB-zminB, theDDZ, midB, widB, inter_nb_2 ); std::vector mid, wid; Interpolate( midA, midB, theNbSteps, mid, isAddSecond ); @@ -762,14 +883,14 @@ std::vector HYDROData_DTM::Interpolate for( size_t i=0; i& theProfiles, + ( const std::vector& theProfiles, double theDDZ, double theSpatialStep, AltitudePoints& theLeft, AltitudePoints& theRight, @@ -782,7 +903,7 @@ HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate return points; std::vector distances; - Handle_Geom2d_BSplineCurve aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances ); + Handle(Geom2d_BSplineCurve) aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances ); if( aHydraulicAxis.IsNull() ) return points; @@ -800,10 +921,10 @@ HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond, inter_nb_1, inter_nb_2 ); int lps = local_points.size(); - if (inter_nb_1 >= 2) + if (inter_nb_1 > 2) invalInd.insert(i); - if (inter_nb_2 >= 2) + if (inter_nb_2 > 2) invalInd.insert(i+1); // 2. Put all points into the global container @@ -836,3 +957,31 @@ HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate } return points; } + +int HYDROData_DTM::EstimateNbPoints( const std::vector& theProfiles, + double theDDZ, double theSpatialStep ) +{ + size_t n = theProfiles.size(); + if( n<=1 ) + return 0; + if( theDDZ<1E-6 || theSpatialStep<1E-6 ) + return 1 << 20; + + std::vector distances; + Handle(Geom2d_BSplineCurve) aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances ); + if( aHydraulicAxis.IsNull() ) + return 0; + + double aCompleteDistance = distances[n-1]; + int aNbSteps = int( aCompleteDistance / theSpatialStep ) + 1; + gp_Pnt aLowest; + gp_Vec2d aDir; + double aZMin, aZMax; + GetProperties( theProfiles[0], aLowest, aDir, aZMin, aZMax ); + int aNbZSteps = (aZMax-aZMin)/theDDZ; + + if( aNbSteps > ( 1<<16 ) || aNbZSteps > ( 1<<16 ) ) + return 1 << 20; + + return aNbSteps * aNbZSteps; +}