X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FHYDROData%2FHYDROData_DTM.cxx;h=207530f8df59e3a0c37bc80081e3315a3197bbcc;hb=bd81def133c235351f4c18834c6660450f16b9ee;hp=4db3a948b63f60fa52fc018b1ea6fa81c9ff54f3;hpb=c4738cc6d2fb23c60a1922847991007b3ac4e6d7;p=modules%2Fhydro.git diff --git a/src/HYDROData/HYDROData_DTM.cxx b/src/HYDROData/HYDROData_DTM.cxx index 4db3a948..207530f8 100644 --- a/src/HYDROData/HYDROData_DTM.cxx +++ b/src/HYDROData/HYDROData_DTM.cxx @@ -18,21 +18,52 @@ #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 ) +#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 +#include +IMPLEMENT_STANDARD_RTTIEXT( HYDROData_DTM, HYDROData_Bathymetry ) -HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv ) - : myXcurv( theXCurv ) +HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv, const gp_Vec2d& theProfileDir, double theDeltaZ, double theMaxZ ) + : myXcurv( theXCurv ), myProfileDir( theProfileDir ), myDeltaZ( theDeltaZ ), myMaxZ (theMaxZ) { } @@ -45,12 +76,32 @@ double HYDROData_DTM::CurveUZ::Xcurv() const return myXcurv; } +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() ); - 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()) + { + Handle(Geom_BSplineCurve) C = anInterpolator.Curve(); + E = BRepBuilderAPI_MakeEdge(C).Edge(); + } } +TopTools_IndexedMapOfOrientedShape HYDROData_DTM::Create3DShape(const AltitudePoints& left, + const AltitudePoints& right, + const std::vector& main_profiles) +{ + TopTools_IndexedMapOfOrientedShape ll; + //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; + 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 (!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 + + return ll; +} + + void HYDROData_DTM::Update() { - HYDROData_SequenceOfObjects objs = GetProfiles(); - int aLower = objs.Lower(), anUpper = objs.Upper(); - size_t n = anUpper-aLower+1; + AltitudePoints points; + TopoDS_Shape Out3dPres; + TopoDS_Shape Out2dPres; + TopoDS_Shape OutLeftB; + TopoDS_Shape OutRightB; + TopoDS_Shape OutInlet; + TopoDS_Shape OutOutlet; + + HYDROData_SequenceOfObjects objs = GetProfiles(); + double ddz = GetDDZ(); + double step = GetSpatialStep(); + 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); + SetShape( DataTag_OutletShape, OutOutlet ); + SetShape( DataTag_3DShape, Out3dPres ); + SetShape( DataTag_2DShape, Out2dPres ); + + HYDROData_Bathymetry::Update(); +} - std::vector profiles; - profiles.reserve( n ); +void HYDROData_DTM::GetPresentationShapes( TopoDS_Shape& Out3dPres, + TopoDS_Shape& Out2dPres, + TopoDS_Shape& OutLeftB, + TopoDS_Shape& OutRightB, + TopoDS_Shape& OutInlet, + TopoDS_Shape& OutOutlet ) +{ + //without update! + OutLeftB = GetShape( DataTag_LeftBankShape); + OutRightB = GetShape( DataTag_RightBankShape); + OutInlet = GetShape( DataTag_InletShape); + OutOutlet = GetShape( DataTag_OutletShape ); + Out3dPres = GetShape( DataTag_3DShape ); + Out2dPres = GetShape( DataTag_2DShape ); +} +void HYDROData_DTM::CreateProfilesFromDTM (const HYDROData_SequenceOfObjects& InpProfiles, + double ddz, + double step, + AltitudePoints& points, + TopoDS_Shape& Out3dPres, + TopoDS_Shape& Out2dPres, + TopoDS_Shape& OutLeftB, + TopoDS_Shape& OutRightB, + TopoDS_Shape& OutInlet, + TopoDS_Shape& OutOutlet, + bool Create3dPres, + bool Create2dPres, + std::set& InvInd, + int thePntsLimit, + bool& WireIntersections) +{ + int aLower = InpProfiles.Lower(), anUpper = InpProfiles.Upper(); + size_t n = anUpper - aLower + 1; + + std::vector profiles; + profiles.reserve( n ); for( int i=aLower; i<=anUpper; i++ ) { - Handle(HYDROData_Profile) aProfile = Handle(HYDROData_Profile)::DownCast( objs.Value( i ) ); + Handle(HYDROData_Profile) aProfile = Handle(HYDROData_Profile)::DownCast( InpProfiles.Value( i ) ); if( !aProfile.IsNull() ) profiles.push_back( aProfile ); } - - double ddz = GetDDZ(); - double step = GetSpatialStep(); const double EPS = 1E-3; - AltitudePoints points; + AltitudePoints left; + 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 ) - points = Interpolate( profiles, ddz, step ); - SetAltitudePoints( points ); + CreateProfiles(profiles, ddz, step, left, right, points, main_profiles, + Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, Create3dPres, Create2dPres, InvInd, WireIntersections ); } +bool HYDROData_DTM::GetPlanarFaceFromBanks( const TopoDS_Edge& LB, const TopoDS_Edge& RB, TopoDS_Face& outF, + TopTools_SequenceOfShape* Boundr) +{ + 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; + 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(); + 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) + { + 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, + double theDDZ, + double theSpatialStep, + AltitudePoints& theOutLeft, + AltitudePoints& theOutRight, + AltitudePoints& theOutPoints, + std::vector& theOutMainProfiles, + TopoDS_Shape& Out3dPres, + TopoDS_Shape& Out2dPres, + TopoDS_Shape& OutLeftB, + TopoDS_Shape& OutRightB, + TopoDS_Shape& OutInlet, + TopoDS_Shape& OutOutlet, + bool Create3dPres, + bool Create2dPres, + 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); + + if (ll.IsEmpty()) + return; + BRep_Builder BB; + TopoDS_Compound cmp; + BB.MakeCompound(cmp); + for (int i = 1; i <= ll.Extent(); i++) + BB.Add(cmp, ll(i)); + + Out3dPres = cmp; + + //same order as in HYDROData_DTM::Update() + OutLeftB = ll(1); + OutRightB = ll(ll.Extent()); + OutInlet = ll(2); + OutOutlet = ll(ll.Extent() - 1); + + if (Create2dPres) + { + TopoDS_Face outF; + ProjStat = GetPlanarFaceFromBanks(TopoDS::Edge(OutLeftB), TopoDS::Edge(OutRightB), outF, NULL); + Out2dPres = outF; + }; + } +} -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(); @@ -160,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++ ) { @@ -220,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 ) ); } @@ -239,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 ); @@ -281,7 +606,7 @@ Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis( Standard_Real aParamFirst = anAdaptor.FirstParameter(), aParamLast = anAdaptor.LastParameter(); for( size_t i = 1; i <= n; i++ ) { - gp_Pnt2d aPnt = points->Value( i ); + gp_Pnt2d aPnt = points->Value( (Standard_Integer)i ); Geom2dAPI_ProjectPointOnCurve aProject( aPnt, aResult ); Standard_Real aParam = aProject.LowerDistanceParameter(); double aDistance = GCPnts_AbscissaPoint::Length( anAdaptor, aParamFirst, aParam ); @@ -291,20 +616,19 @@ Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis( return aResult; } -std::vector HYDROData_DTM::ProfileToParametric( - const Handle_HYDROData_Profile& theProfile, - double& theUMin, double& theUMax ) +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; - gp_Vec2d aDir; double zmin, zmax; - GetProperties( theProfile, aLowest, aDir, 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( aDir.X(), aDir.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 ); @@ -331,50 +655,43 @@ std::vector HYDROData_DTM::ProfileToParametric( } -double CalcGC( const std::vector& intersections ) +bool CalcMidWidth( const std::set& intersections, double& theMid, double& theWid ) { - double u = 0; - size_t n = intersections.size(); - for( size_t i = 0; i < n; i++ ) - { - u += intersections[i].X(); - } - u /= n; - return u; -} - -double CalcWidth( const std::vector& intersections ) -{ - double umin = std::numeric_limits::max(), + double umin = DBL_MAX, umax = -umin; size_t n = intersections.size(); - if( n <= 1 ) - return 0; + 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 ) umax = u; } - return umax-umin; + theMid = ( umin+umax )/2; + theWid = umax-umin; + 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, - double theTolerance ) + CurveUZ& theWidthCurve, + int& intersection_nb, + double theTolerance) { - double aDblMax = std::numeric_limits::max(), + double aDblMax = DBL_MAX, aUMin = aDblMax, aUMax = -aUMin, aVMax = 1000000; - std::vector curves = ProfileToParametric( theProfile, aUMin, aUMax ); + gp_Vec2d aProfileDir; + std::vector curves = ProfileToParametric( theProfile, aUMin, aUMax, aProfileDir ); size_t n = curves.size(); if( n==0 ) @@ -394,34 +711,38 @@ void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& thePr curves.push_back( aT2 ); int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 ); - theMidPointCurve = CurveUZ( theXCurv ); + theMidPointCurve = CurveUZ( theXCurv, aProfileDir, theMinZ, theTopZ); theMidPointCurve.reserve( psize ); - theWidthCurve = CurveUZ( theXCurv ); + 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() ); } - if( intersections.size() >= 2 ) + intersection_nb = intersections.size(); + if( intersection_nb >= 1 ) { - double u_mid = CalcGC( intersections ); + 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; theMidPointCurve.push_back( p_mid ); - double u_wid = CalcWidth( intersections ); PointUZ p_wid; p_wid.U = u_wid; p_wid.Z = z; @@ -448,69 +769,107 @@ void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCur if( isAddSecond ) theInterpolation.push_back( theCurveB ); } - -void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxis, +#include +void HYDROData_DTM::CurveTo3D( const Handle(Geom2d_BSplineCurve)& theHydraulicAxis, const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve, AltitudePoints& thePoints ) { Geom2dAdaptor_Curve anAdaptor( theHydraulicAxis ); + TopoDS_Edge E2d = BRepLib_MakeEdge2d(theHydraulicAxis).Edge(); GCPnts_AbscissaPoint ap( anAdaptor, theMidCurve.Xcurv(), anAdaptor.FirstParameter() ); double aParam = ap.Parameter(); gp_Pnt2d point; - gp_Vec2d tangent, profile_dir; - anAdaptor.D1( aParam, point, tangent ); - profile_dir.SetCoord( tangent.Y(), -tangent.X() ); + anAdaptor.D0( aParam, point ); + gp_Vec2d profile_dir = theMidCurve.ProfileDir(); + //gp_Dir tangent_n( -profile_dir.Y(), profile_dir.X(), dz ); profile_dir.Normalize(); - + size_t n = theMidCurve.size(); - for( size_t i=0; i sorted_points; + bool isOnTop = false; + for( size_t i=0; ib ) - return a; - else - return b; + thePoints.reserve( sorted_points.size() ); + const double EPS = 1E-12; + std::map::const_iterator it = sorted_points.begin(), last = sorted_points.end(); + for( ; it!=last; it++ ) + if( thePoints.empty() || thePoints.back().SquareDistance( it->second ) > EPS ) + thePoints.push_back( it->second ); } -HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate - ( const Handle_Geom2d_BSplineCurve& theHydraulicAxis, - const Handle_HYDROData_Profile& theProfileA, +//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, double theXCurvA, - const Handle_HYDROData_Profile& theProfileB, + const Handle(HYDROData_Profile)& theProfileB, double theXCurvB, - double theDDZ, int theNbSteps, bool isAddSecond ) + double theDDZ, int theNbSteps, bool isAddSecond, + int& inter_nb_1, int& inter_nb_2) { double zminA, zmaxA, zminB, zmaxB; 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 zmin = max( zminA, zminB ); - double zmax = max( zmaxA, zmaxB ); + + double hmax = zmaxA-zminA > zmaxB-zminB ? zmaxA-zminA : zmaxB-zminB; + + //double dz = zminB - zminA; + //double zmin = min( zminA, zminB ); + //double zmax = max( zmaxA, zmaxB ); - CurveUZ midA(0), midB(0); - CurveUZ widA(0), widB(0); + 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 ); - ProfileDiscretization( theProfileB, theXCurvB, zmin, zmax, theDDZ, midB, widB ); + 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 ); @@ -518,16 +877,25 @@ HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate size_t p = mid.size(); size_t q = p>0 ? 2*mid[0].size() : 1; - AltitudePoints points; - points.reserve( p*q ); + std::vector points; + points.resize( p ); + for( size_t i=0; i& theProfiles, - double theDDZ, double theSpatialStep ) + ( const std::vector& theProfiles, + double theDDZ, double theSpatialStep, + AltitudePoints& theLeft, + AltitudePoints& theRight, + std::vector& theMainProfiles, + std::set& invalInd) { AltitudePoints points; size_t n = theProfiles.size(); @@ -535,21 +903,85 @@ 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; + theMainProfiles.reserve( n ); + for( size_t i=0, n1=n-1; i local_points = Interpolate( aHydraulicAxis, theProfiles[i], distances[i], + 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) + invalInd.insert(i); + + if (inter_nb_2 > 2) + invalInd.insert(i+1); + + // 2. Put all points into the global container + for( size_t j=0; j& 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; +}