From: asl Date: Wed, 7 Sep 2016 07:35:30 +0000 (+0300) Subject: implementation of algorithms for DTM X-Git-Tag: v1.6~80 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=5124dac22315ff6a5e6cce6509a5f86ea15e1b78;p=modules%2Fhydro.git implementation of algorithms for DTM --- diff --git a/src/HYDROData/HYDROData_DTM.cxx b/src/HYDROData/HYDROData_DTM.cxx index a80edf22..384f672a 100644 --- a/src/HYDROData/HYDROData_DTM.cxx +++ b/src/HYDROData/HYDROData_DTM.cxx @@ -21,10 +21,59 @@ #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 HYDROData_DTM::CurveUZ::Xcurv() const +{ + return myXcurv; +} + +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; iGetBottomPoint(); @@ -52,24 +101,34 @@ void GetProperties( const Handle_HYDROData_Profile& theProfile, theDir = gp_Vec2d( x, y ); } -inline gp_Pnt2d To2D( const gp_Pnt& thePnt, const gp_Trsf& theTr ) +inline gp_Pnt2d To2D( const gp_Pnt& thePnt, const gp_Trsf& theTr, + double& theUMin, double& theUMax ) { gp_Pnt p = thePnt.Transformed( theTr ); - return gp_Pnt2d( p.X(), p.Z() ); + double u = p.X(); + double z = p.Z(); + if( utheUMax ) + theUMax = u; + return gp_Pnt2d( u, z ); } -Handle(TColgp_HArray1OfPnt2d) To2D( const TColgp_Array1OfPnt& thePoints, const gp_Trsf& theTr ) +Handle(TColgp_HArray1OfPnt2d) To2D( const TColgp_Array1OfPnt& thePoints, + const gp_Trsf& theTr, + double& theUMin, double& theUMax ) { int low = thePoints.Lower(), up = thePoints.Upper(); Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( low, up ); for( int i=low; i<=up; i++ ) - points->SetValue( i, To2D( thePoints.Value( i ), theTr ) ); + points->SetValue( i, To2D( thePoints.Value( i ), theTr, theUMin, theUMax ) ); return points; } Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve, Standard_Real theFirst, Standard_Real theLast, - const gp_Trsf& theTr ) + const gp_Trsf& theTr, + double& theUMin, double& theUMax ) { if( theCurve->IsKind( STANDARD_TYPE( Geom_Line ) ) ) { @@ -78,8 +137,8 @@ Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve, theCurve->D0( theLast, aLastPnt ); gp_Pnt2d - aFirst2d = To2D( aFirstPnt, theTr ), - aLast2d = To2D( aLastPnt, theTr ); + aFirst2d = To2D( aFirstPnt, theTr, theUMin, theUMax ), + 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() ) ); @@ -90,7 +149,7 @@ Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve, { Handle(Geom_BSplineCurve) aSpline = Handle(Geom_BSplineCurve)::DownCast( theCurve ); - Handle(TColgp_HArray1OfPnt2d) poles = To2D( aSpline->Poles(), theTr ); + Handle(TColgp_HArray1OfPnt2d) poles = To2D( aSpline->Poles(), theTr, theUMin, theUMax ); const TColStd_Array1OfReal& knots = aSpline->Knots(); const TColStd_Array1OfInteger& multiplicities = aSpline->Multiplicities(); int aDegree = aSpline->Degree(); @@ -102,7 +161,8 @@ Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve, } Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis( - const std::vector& theProfiles ) + const std::vector& theProfiles, + std::vector& theDistances ) { size_t n = theProfiles.size(); Handle_Geom2d_BSplineCurve aResult; @@ -114,9 +174,12 @@ Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis( for( size_t i = 1; i <= n; i++ ) { Handle_HYDROData_Profile aProfile = theProfiles[i-1]; + aProfile->Update(); + gp_Pnt aLowest; gp_Vec2d aTangent; - GetProperties( aProfile, aLowest, aTangent, true ); + double zmin, zmax; + GetProperties( aProfile, aLowest, aTangent, true, zmin, zmax ); aTangent.Normalize(); points->SetValue( (int)i, gp_Pnt2d( aLowest.X(), aLowest.Y() ) ); @@ -128,19 +191,24 @@ Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis( anInterpolator.Load( tangents, flags ); anInterpolator.Perform(); if( anInterpolator.IsDone() ) + { aResult = anInterpolator.Curve(); + //TODO: fill the distances vector + } 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 ) { std::vector curves; - theProfile->Update(); // Transformation of the coordinate systems gp_Pnt aLowest; gp_Vec2d aDir; - GetProperties( theProfile, aLowest, aDir, false ); + double zmin, zmax; + GetProperties( theProfile, aLowest, aDir, 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 ) ); @@ -162,30 +230,195 @@ std::vector HYDROData_DTM::ProfileToParametric( const Handl Handle(Geom_Curve) aCurve = BRep_Tool::Curve( anEdge, aLoc, aFirst, aLast ); // Convert the curve to 2d CS - Handle(Geom2d_Curve) aCurve2d = CurveTo2D( aCurve, aFirst, aLast, aTransf ); + Handle(Geom2d_Curve) aCurve2d = CurveTo2D( aCurve, aFirst, aLast, aTransf, theUMin, theUMax ); if( !aCurve2d.IsNull() ) curves.push_back( aCurve2d ); } return curves; } + +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 ) +{ + double umin = std::numeric_limits::max(), + umax = -umin; + + size_t n = intersections.size(); + if( n <= 1 ) + return 0; + + for( size_t i = 0; i < n; i++ ) + { + double u = intersections[i].X(); + if( uumax ) + umax = u; + } + return umax-umin; +} + void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& theProfile, - double theMinZ, double theMaxZ, double theDDZ, + double theXCurv, double theMinZ, double theMaxZ, double theDDZ, CurveUZ& theMidPointCurve, - CurveUZ& theWidthCurve ) + CurveUZ& theWidthCurve, + double theTolerance ) { - /*for( double z = theMinZ; z<=theMaxZ; z += theDDZ ) + double aDblMax = std::numeric_limits::max(), + aUMin = aDblMax, + aUMax = -aUMin; + + std::vector curves = ProfileToParametric( theProfile, aUMin, aUMax ); + size_t n = curves.size(); + + if( n==0 ) + return; + + // we add the "virtual" vertical lines to simulate the intersection with profile + gp_Pnt2d aFirst, aLast; + 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_Line) aV2 = new Geom2d_Line( aLast, gp_Dir2d( 0, 1 ) ); + Handle(Geom2d_TrimmedCurve) aT2 = new Geom2d_TrimmedCurve( aV2, 0.0, aDblMax ); + + curves.push_back( aT1 ); + curves.push_back( aT2 ); + + int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 ); + theMidPointCurve = CurveUZ( theXCurv ); + theMidPointCurve.reserve( psize ); + theWidthCurve = CurveUZ( theXCurv ); + 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 ) { + Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z ), gp_Dir2d( 1, 0 ) ); + std::vector 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 ) + { + double u_mid = CalcGC( intersections ); + 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; + theWidthCurve.push_back( p_wid ); + } + } } void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCurveB, - int theNbSteps, std::vector& theInterpolation ) + int theNbSteps, std::vector& theInterpolation, + bool isAddSecond ) { + theInterpolation.clear(); + int d = isAddSecond ? 2 : 1; + theInterpolation.reserve( theNbSteps+d ); + double dt = 1.0 / double( theNbSteps + 1 ); + double t = dt; + theInterpolation.push_back( theCurveA ); + for( int i=0; i& theInterpolation ) +void HYDROData_DTM::CurveTo3d( const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve, + AltitudePoints& thePoints ) { + //TODO +} + +inline double max( double a, double b ) +{ + if( a>b ) + return a; + else + 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 ) +{ + 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 ); + + double zmin = max( zminA, zminB ); + double zmax = max( zmaxA, zmaxB ); + + CurveUZ midA(0), midB(0); + CurveUZ widA(0), widB(0); + + ProfileDiscretization( theProfileA, theXCurvA, zmin, zmax, theDDZ, midA, widA ); + ProfileDiscretization( theProfileB, theXCurvB, zmin, zmax, theDDZ, midB, widB ); + + 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; i& theProfiles, + double theDDZ, double theSpatialStep, AltitudePoints& thePoints ) +{ + size_t n = theProfiles.size(); + if( n<=1 ) + return; + + std::vector distances; + Handle_Geom2d_BSplineCurve aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances ); + if( aHydraulicAxis.IsNull() ) + return; + + for( size_t i=0, n1=n-1; i { public: - CurveUZ( double theXCurv ); + CurveUZ( double theXcurv ); ~CurveUZ(); + double Xcurv() const; + + CurveUZ operator + ( const CurveUZ& ) const; + CurveUZ operator * ( double ) const; + private: - double Xcurv; + double myXcurv; }; - static Handle_Geom2d_BSplineCurve CreateHydraulicAxis( const std::vector& theProfiles ); + static Handle_Geom2d_BSplineCurve CreateHydraulicAxis( + const std::vector& theProfiles, + std::vector& theDistances ); - static std::vector ProfileToParametric( const Handle_HYDROData_Profile& theProfile ); + static std::vector ProfileToParametric( const Handle_HYDROData_Profile& theProfile, + double& theUMin, double& theUMax ); - static void ProfileDiscretization( const Handle_HYDROData_Profile& theProfile, - double theMinZ, double theMaxZ, double theDDZ, + static void ProfileDiscretization( const Handle_HYDROData_Profile& theProfile, + double theXCurv, double theMinZ, double theMaxZ, double theDDZ, CurveUZ& theMidPointCurve, - CurveUZ& theWidthCurve ); + CurveUZ& theWidthCurve, + double theTolerance = 1E-6 ); + static void CurveTo3d( const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve, AltitudePoints& thePoints ); + static void Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCurveB, - int theNbSteps, std::vector& theInterpolation ); - - static void CurveTo3d( const CurveUZ& theCurve, const CurveUZ& theCurveB, - int theNbSteps, std::vector& theInterpolation ); + int theNbSteps, std::vector& theInterpolation, + bool isAddSecond ); + + static void Interpolate( const Handle_HYDROData_Profile& theProfileA, + double theXCurvA, + const Handle_HYDROData_Profile& theProfileB, + double theXCurvB, + double theDDZ, int theNbSteps, bool isAddSecond, + AltitudePoints& thePoints ); + + static void Interpolate( const std::vector& theProfiles, + double theDDZ, double theSpatialStep, AltitudePoints& thePoints ); }; #endif diff --git a/src/HYDRO_tests/test_HYDROData_DTM.cxx b/src/HYDRO_tests/test_HYDROData_DTM.cxx index d1189c74..5786bef3 100644 --- a/src/HYDRO_tests/test_HYDROData_DTM.cxx +++ b/src/HYDRO_tests/test_HYDROData_DTM.cxx @@ -66,11 +66,17 @@ void test_HYDROData_DTM::test_profile_conversion_to_2d() aProfile2->SetLeftPoint( gp_XY( 10, 10 ) ); aProfile2->SetRightPoint( gp_XY( 20, 20 ) ); - std::vector curves1 = HYDROData_DTM::ProfileToParametric( aProfile1 ); - std::vector curves2 = HYDROData_DTM::ProfileToParametric( aProfile2 ); + double aUMin1 = std::numeric_limits::max(), + aUMax1 = -aUMin1, + aUMin2 = aUMin1, + aUMax2 = aUMax1; + std::vector curves1 = HYDROData_DTM::ProfileToParametric( aProfile1, aUMin1, aUMax1 ); + std::vector curves2 = HYDROData_DTM::ProfileToParametric( aProfile2, aUMin2, aUMax2 ); gp_Pnt2d aFirst, aLast; CPPUNIT_ASSERT_EQUAL( 3, (int)curves1.size() ); + CPPUNIT_ASSERT_DOUBLES_EQUAL( -5.303, aUMin1, EPS ); + CPPUNIT_ASSERT_DOUBLES_EQUAL( 8.839, aUMax1, EPS ); curves1[0]->D0( curves1[0]->FirstParameter(), aFirst ); curves1[0]->D0( curves1[0]->LastParameter(), aLast ); CPPUNIT_ASSERT_DOUBLES_EQUAL( -5.303, aFirst.X(), EPS ); @@ -91,6 +97,8 @@ void test_HYDROData_DTM::test_profile_conversion_to_2d() CPPUNIT_ASSERT_DOUBLES_EQUAL( 4.0, aLast.Y(), EPS ); CPPUNIT_ASSERT_EQUAL( 1, (int)curves2.size() ); + CPPUNIT_ASSERT_DOUBLES_EQUAL( -5.303, aUMin2, EPS ); + CPPUNIT_ASSERT_DOUBLES_EQUAL( 8.839, aUMax2, EPS ); Handle(Geom2d_BSplineCurve) aBSpline = Handle(Geom2d_BSplineCurve)::DownCast( curves2[0] ); CPPUNIT_ASSERT_EQUAL( false, (bool)aBSpline.IsNull() ); const TColgp_Array1OfPnt2d& poles = aBSpline->Poles();