X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FHYDROData%2FHYDROData_DTM.cxx;h=df831b619a110a4ba5dd9b61549e7c3be98774fd;hb=2fcf12a721b77f222a9e597c65dc2cafe35db97c;hp=384f672a1348247ef936ac23b5f9523756fb6397;hpb=5124dac22315ff6a5e6cce6509a5f86ea15e1b78;p=modules%2Fhydro.git diff --git a/src/HYDROData/HYDROData_DTM.cxx b/src/HYDROData/HYDROData_DTM.cxx index 384f672a..df831b61 100644 --- a/src/HYDROData/HYDROData_DTM.cxx +++ b/src/HYDROData/HYDROData_DTM.cxx @@ -18,18 +18,49 @@ #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 +#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 ) - : myXcurv( theXCurv ) +HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv, const gp_Vec2d& theProfileDir, double theDeltaZ ) + : myXcurv( theXCurv ), myProfileDir( theProfileDir ), myDeltaZ( theDeltaZ ) { } @@ -42,10 +73,24 @@ 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; +} + HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator + ( const CurveUZ& c ) const { - HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv() ); - size_t n = size(); + HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv(), ProfileDir() + c.ProfileDir(), DeltaZ() + c.DeltaZ() ); + size_t n = size(), n1 = c.size(); + if( n!=n1 ) + { + std::cout << "Warning: different number of points in curves: " << n << ", " << n1 << std::endl; + } res.reserve( n ); for( int i=0; i& main_profiles) +{ + TopTools_IndexedMapOfOrientedShape ll; + TopoDS_Wire LWire, RWire; + PointToWire(left, LWire); + PointToWire(right, RWire); + if (!LWire.IsNull()) + ll.Add(LWire.Oriented(TopAbs_FORWARD)); + + for (int k = 0; k < main_profiles.size(); k++) + { + TopoDS_Wire W; + PointToWire(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)); + //yes, add subshapes in this order (left + profiles + right) + //otherwise the projected wire will be non-manifold + + return ll; +} + + +void HYDROData_DTM::Update() +{ + 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(); +} + +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( InpProfiles.Value( i ) ); + if( !aProfile.IsNull() ) + profiles.push_back( aProfile ); + } + const double EPS = 1E-3; + 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 ) + CreateProfiles(profiles, ddz, step, left, right, points, main_profiles, + Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, Create3dPres, Create2dPres, InvInd, WireIntersections ); +} + +void HYDROData_DTM::ProjWireOnPlane(const TopoDS_Shape& inpWire, const Handle_Geom_Plane& RefPlane, + TopTools_DataMapOfShapeListOfShape* E2PE) +{ + BRep_Builder BB; + + //project shape (edges) on planar face + TopoDS_Face F; + BB.MakeFace(F, RefPlane, Precision::Confusion()); + BRepAlgo_NormalProjection nproj(F); + nproj.Add(inpWire); + nproj.SetDefaultParams(); + nproj.Build(); + if(!nproj.IsDone()) + return; + + TopoDS_Shape projRes = nproj.Projection(); + + // unite all vertexes/edges from projected result + BOPAlgo_Builder anAlgo; + TopExp_Explorer exp(projRes, TopAbs_EDGE); + for (;exp.More(); exp.Next()) + { + const TopoDS_Edge& E = TopoDS::Edge(exp.Current()); + if (E.Orientation() != TopAbs_INTERNAL) + anAlgo.AddArgument(E); + } + anAlgo.Perform(); + int stat = anAlgo.ErrorStatus(); + TopoDS_Shape projResConn = anAlgo.Shape(); + + // make wire => vertexes and edges should be untouched after this operation! + exp.Init(projResConn, TopAbs_EDGE); + //TopTools_ListOfShape llE; + //TopoDS_Wire RW; + //BRepLib_MakeWire WM; + + //for (;exp.More();exp.Next()) + // llE.Append(exp.Current()); + // + //WM.Add(llE); + //outWire = WM.Wire(); + + //outWire.Orientation(inpWire.Orientation()); //take from the original wire + + //history mode: edge to projected edges + if (E2PE) + { + TopExp_Explorer ex(inpWire, TopAbs_EDGE); + for (;ex.More();ex.Next()) + { + const TopoDS_Edge& CE = TopoDS::Edge(ex.Current()); + TopTools_ListOfShape NEL; + const TopTools_ListOfShape& LS = nproj.Generated(CE); + TopTools_ListIteratorOfListOfShape it(LS); + for (;it.More();it.Next()) + { + const TopoDS_Shape& PCE = it.Value(); + TopTools_ListOfShape PLS = anAlgo.Modified(PCE); + if (PLS.IsEmpty()) + PLS.Append(PCE); + TopTools_ListIteratorOfListOfShape itp(PLS); + for (;itp.More();itp.Next()) + NEL.Append(itp.Value()); + } + + E2PE->Bind(CE, NEL); + } + } +} +#include +bool HYDROData_DTM::Get2dFaceFrom3dPres(const TopoDS_Compound& cmp, TopoDS_Face& outF, + TopTools_SequenceOfShape* Boundr, std::set ind ) +{ + //ind : set of indices (starts with 0). index == number of boundary (inlet, outlet, etc..) + //in compound cmp. + //if Boundr is not null => this method will return sequence of boundary wires (inlet, outlet...) + + Handle_Geom_Plane refpl = new Geom_Plane(gp_Pnt(0,0,0), gp_Dir(0,0,1)); + TopTools_DataMapOfShapeListOfShape E2PE; + ProjWireOnPlane(cmp, refpl, &E2PE); + TopTools_ListOfShape ELL; + + TopoDS_Iterator it(cmp); + int i = 0; + for (;it.More(); it.Next()) + { + const TopoDS_Wire& W = TopoDS::Wire(it.Value()); + if (W.Orientation() != TopAbs_INTERNAL) + { + TopoDS_Wire PW; + TopExp_Explorer ex(W, TopAbs_EDGE); + TopTools_ListOfShape LEpW; + TopTools_ListOfShape LEpWDummy; + for (;ex.More();ex.Next()) + { + const TopoDS_Edge& CE = TopoDS::Edge(ex.Current()); + TopTools_ListOfShape LS = E2PE.Find(CE); + LEpW.Append(LS); + } + + if (ind.count(i) != 0) + { + BRepLib_MakeWire WM; + WM.Add(LEpW); + const TopoDS_Wire& WMW = WM.Wire(); + //assume that wire is a straight line, + //take first and last vertex and make simple edge (RE) + TopoDS_Vertex VF, VL; + TopExp::Vertices(WMW, VF, VL); + TopoDS_Edge RE = BRepLib_MakeEdge(VF, VL).Edge(); + if (RE.IsNull()) + { + LEpWDummy = LEpW; //LEpW will be nullified after appending to ELL + ELL.Append(LEpW); + } + else + { + LEpWDummy.Append(RE); + ELL.Append(RE); + } + //TODO: in the new version of OCCT, USD can process separate wires + //ShapeUpgrade_UnifySameDomain USD(WMW, 1U, 0U, 1U); //concat bsplines + //USD.Build(); + //const TopoDS_Shape& RSU = USD.Shape(); + //TopExp_Explorer exp(RSU, TopAbs_EDGE); + //TopTools_ListOfShape DummyL; + //for (;it.More();it.Next()) + // DummyL.Append(it.Value()); + //if (DummyL.Extent() == 1) + // ELL.Append(DummyL.First()); //if one edge => accept this result + //else + // ELL.Append(LEpW); //else put 'as is' + } + else + { + LEpWDummy = LEpW; + ELL.Append(LEpW); + } + + if (Boundr) + { + //make inlet, outlet, left/tight banks [wires] + //shouldn't change topology of the edges + BRepLib_MakeWire IWM; + IWM.Add(LEpWDummy); + Boundr->Append(IWM.Wire()); + } + } + i++; + } + + //make primary wire + BRepLib_MakeWire WME; + WME.Add(ELL); + + const TopoDS_Wire& resW = WME.Wire(); + BRepBuilderAPI_MakeFace mf(refpl, resW, true); //check inside is true by def + outF = mf.Face(); + + ShapeAnalysis_Wire WA(resW, outF, Precision::Confusion()); + bool res = WA.CheckSelfIntersection(); //TODO check that this return false if OK + return res; + + ///!!! 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++) + { + TopoDS_Wire outIW; + const TopoDS_Wire& W = TopoDS::Wire(IntW(i)); + ProjWireOnPlane(W, refpl, outIW); + BB.Add(outF, outIW); + }*/ +} + +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& WireIntersections) +{ + if (theProfiles.empty()) + return; + theOutPoints = Interpolate( theProfiles, theDDZ, theSpatialStep, theOutLeft, theOutRight, theOutMainProfiles, InvInd ); + //note that if Create3dPres is false => Create2dPres flag is meaningless! + 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; + WireIntersections = Get2dFaceFrom3dPres(cmp, outF); //__TODO + Out2dPres = outF; + }; + } +} + -void 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 ) @@ -99,6 +543,20 @@ void GetProperties( const Handle_HYDROData_Profile& theProfile, theDir = gp_Vec2d( -y, x ); else theDir = gp_Vec2d( x, y ); + + HYDROData_Profile::ProfilePoints points = theProfile->GetProfilePoints(); + int lo = points.Lower(); + int up = points.Upper(); + theZMin = std::numeric_limits::max(); + theZMax = -theZMin; + for( int i=lo; i<=up; i++ ) + { + double z = points.Value( i ).Z(); + if( z>theZMax ) + theZMax = z; + if( zValue( (Standard_Integer)i ); + Geom2dAPI_ProjectPointOnCurve aProject( aPnt, aResult ); + Standard_Real aParam = aProject.LowerDistanceParameter(); + double aDistance = GCPnts_AbscissaPoint::Length( anAdaptor, aParamFirst, aParam ); + theDistances.push_back( aDistance ); + } } return aResult; } std::vector HYDROData_DTM::ProfileToParametric( const Handle_HYDROData_Profile& theProfile, - double& theUMin, double& theUMax ) + double& theUMin, double& theUMax, gp_Vec2d& theDir ) { 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, false, 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 ); @@ -238,49 +709,43 @@ std::vector HYDROData_DTM::ProfileToParametric( } -double CalcGC( const std::vector& intersections ) -{ - 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 ) +bool CalcMidWidth( const std::set& intersections, double& theMid, double& theWid ) { double umin = std::numeric_limits::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, CurveUZ& theMidPointCurve, - CurveUZ& theWidthCurve, - double theTolerance ) + CurveUZ& theWidthCurve, + int& intersection_nb, + double theTolerance) { double aDblMax = std::numeric_limits::max(), aUMin = aDblMax, - aUMax = -aUMin; + 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 ) @@ -291,43 +756,47 @@ void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& thePr curves[0]->D0( curves[0]->FirstParameter(), aFirst ); curves[n-1]->D0( curves[n-1]->LastParameter(), aLast ); Handle(Geom2d_Line) aV1 = new Geom2d_Line( aFirst, gp_Dir2d( 0, 1 ) ); - Handle(Geom2d_TrimmedCurve) aT1 = new Geom2d_TrimmedCurve( aV1, 0.0, aDblMax ); + Handle(Geom2d_TrimmedCurve) aT1 = new Geom2d_TrimmedCurve( aV1, 0.0, aVMax ); Handle(Geom2d_Line) aV2 = new Geom2d_Line( aLast, gp_Dir2d( 0, 1 ) ); - Handle(Geom2d_TrimmedCurve) aT2 = new Geom2d_TrimmedCurve( aV2, 0.0, aDblMax ); + Handle(Geom2d_TrimmedCurve) aT2 = new Geom2d_TrimmedCurve( aV2, 0.0, aVMax ); curves.push_back( aT1 ); curves.push_back( aT2 ); int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 ); - theMidPointCurve = CurveUZ( theXCurv ); + theMidPointCurve = CurveUZ( theXCurv, aProfileDir, theMinZ ); theMidPointCurve.reserve( psize ); - theWidthCurve = CurveUZ( theXCurv ); + theWidthCurve = CurveUZ( theXCurv, aProfileDir, theMinZ ); 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]; Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance ); - for( int k=0, m=anIntersect.NbPoints(); k= 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; @@ -354,11 +823,46 @@ void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCur if( isAddSecond ) theInterpolation.push_back( theCurveB ); } - -void HYDROData_DTM::CurveTo3d( const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve, +#include +void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxis, + const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve, AltitudePoints& thePoints ) { - //TODO + 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; + 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(); + std::map sorted_points; + for( size_t i=0; i::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 ); } inline double max( double a, double b ) @@ -369,12 +873,24 @@ inline double max( double a, double b ) return b; } -void HYDROData_DTM::Interpolate( const Handle_HYDROData_Profile& theProfileA, - double theXCurvA, - const Handle_HYDROData_Profile& theProfileB, - double theXCurvB, - double theDDZ, int theNbSteps, bool isAddSecond, - AltitudePoints& thePoints ) +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, + double theXCurvB, + double theDDZ, int theNbSteps, bool isAddSecond, + int& inter_nb_1, int& inter_nb_2) { double zminA, zmaxA, zminB, zmaxB; gp_Pnt lowestA, lowestB; @@ -383,42 +899,130 @@ void HYDROData_DTM::Interpolate( const Handle_HYDROData_Profile& theProfileA, GetProperties( theProfileA, lowestA, dirA, false, zminA, zmaxA ); GetProperties( theProfileB, lowestB, dirB, false, zminB, zmaxB ); - double zmin = max( zminA, zminB ); - double zmax = max( zmaxA, zmaxB ); + + double hmax = max( 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), midB(0, gp_Vec2d(), 0); + CurveUZ widA(0, gp_Vec2d(), 0), widB(0, gp_Vec2d(), 0); - ProfileDiscretization( theProfileA, theXCurvA, zmin, zmax, theDDZ, midA, widA ); - ProfileDiscretization( theProfileB, theXCurvB, zmin, zmax, theDDZ, midB, widB ); + ProfileDiscretization( theProfileA, theXCurvA, zminA, zminA+hmax, theDDZ, midA, widA, inter_nb_1 ); + ProfileDiscretization( theProfileB, theXCurvB, zminB, zminB+hmax, theDDZ, midB, widB, inter_nb_2 ); std::vector mid, wid; Interpolate( midA, midB, theNbSteps, mid, isAddSecond ); Interpolate( widA, widB, theNbSteps, wid, isAddSecond ); - size_t q = mid.size(); - for( size_t i=0; i0 ? 2*mid[0].size() : 1; + std::vector points; + points.resize( p ); + + for( size_t i=0; i& theProfiles, - double theDDZ, double theSpatialStep, AltitudePoints& thePoints ) +HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate + ( 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(); if( n<=1 ) - return; + return points; std::vector distances; Handle_Geom2d_BSplineCurve aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances ); if( aHydraulicAxis.IsNull() ) - return; + 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, true, aZMin, aZMax ); + int aNbZSteps = (aZMax-aZMin)/theDDZ; + + if( aNbSteps > ( 1<<16 ) || aNbZSteps > ( 1<<16 ) ) + return 1 << 20; + + return aNbSteps * aNbZSteps; }