X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FHYDROData%2FHYDROData_StreamAltitude.cxx;h=2360eddd0fce93aa735dd9b99d12bfe530995f6e;hb=e35b4caf3341d2cece2b7921f349acd17fb89351;hp=60452989426cf7a635598b66210b67c6184f2ee6;hpb=5395b3a6f169302c35e0e9683915f1136ba99a0b;p=modules%2Fhydro.git diff --git a/src/HYDROData/HYDROData_StreamAltitude.cxx b/src/HYDROData/HYDROData_StreamAltitude.cxx index 60452989..2360eddd 100644 --- a/src/HYDROData/HYDROData_StreamAltitude.cxx +++ b/src/HYDROData/HYDROData_StreamAltitude.cxx @@ -51,13 +51,17 @@ #include #include +#include + +#define _DEVDEBUG_ +#include "HYDRO_trace.hxx" #ifdef DEB_CLASS2D #include #include #include #endif -IMPLEMENT_STANDARD_HANDLE(HYDROData_StreamAltitude, HYDROData_IAltitudeObject) + IMPLEMENT_STANDARD_RTTIEXT(HYDROData_StreamAltitude, HYDROData_IAltitudeObject) HYDROData_StreamAltitude::HYDROData_StreamAltitude() @@ -102,6 +106,48 @@ bool IsPointBetweenEdges( const gp_Pnt& aFirstPnt1, const gp_Pnt& aLastPnt1, return aPointState != TopAbs_OUT; } +/** compute edge1^(firstPnt1,point) and edge2^(firstPnt2,point): + * if the z values are opposite, the point is between the edges. + * We must also check if the point is near the edges + * (inside the circle defined by the edges: approximation) + * to discriminate false positives in sinuous cases + */ +bool IsPointBetweenEdges2( const gp_Pnt& aFirstPnt1, const gp_Pnt& aLastPnt1, + const gp_Pnt& aFirstPnt2, const gp_Pnt& aLastPnt2, + const gp_Pnt& thePoint) { + double x1 = aLastPnt1.X() - aFirstPnt1.X(); // v1 + double y1 = aLastPnt1.Y() - aFirstPnt1.Y(); + double x2 = aLastPnt2.X() - aFirstPnt2.X(); // v2 + double y2 = aLastPnt2.Y() - aFirstPnt2.Y(); + double xa = thePoint.X() - aFirstPnt1.X(); // va + double ya = thePoint.Y() - aFirstPnt1.Y(); + double xb = thePoint.X() - aFirstPnt2.X(); // vb + double yb = thePoint.Y() - aFirstPnt2.Y(); + double z1 = x1*ya -xa*y1; // v1^va: z component + double z2 = x2*yb -xb*y2; // v2^vb: z component + bool isBetween = true; + if (((z1<0) && (z2<0)) || ((z1>=0) && (z2>=0))) + { + isBetween = false; + } + if (isBetween) + { + double xg = (aFirstPnt1.X() + aLastPnt1.X() + aFirstPnt2.X() + aLastPnt2.X())/4; + double yg = (aFirstPnt1.Y() + aLastPnt1.Y() + aFirstPnt2.Y() + aLastPnt2.Y())/4; + double df1 = (aFirstPnt1.X()-xg)*(aFirstPnt1.X()-xg) + (aFirstPnt1.Y()-yg)*(aFirstPnt1.Y()-yg); + double dl1 = (aLastPnt1.X()-xg)*(aLastPnt1.X()-xg) + (aLastPnt1.Y()-yg)*(aLastPnt1.Y()-yg); + double df2 = (aFirstPnt2.X()-xg)*(aFirstPnt2.X()-xg) + (aFirstPnt2.Y()-yg)*(aFirstPnt2.Y()-yg); + double dl2 = (aLastPnt2.X()-xg)*(aLastPnt2.X()-xg) + (aLastPnt2.Y()-yg)*(aLastPnt2.Y()-yg); + double r2 = std::max(df1,dl1); + r2 = std::max(r2,df2); + r2 = std::max(r2,dl2); + double d2 = (thePoint.X()-xg)*(thePoint.X()-xg) + (thePoint.Y()-yg)*(thePoint.Y()-yg); + if (d2 > r2) + isBetween = false; + } + return isBetween; +} + Standard_Real getAltitudeFromWire( const TopoDS_Wire& theWire, const Standard_Real& theLeftDist, const Standard_Real& theRightDist ) @@ -171,8 +217,9 @@ bool HYDROData_StreamAltitude::getBoundaryWiresForPoint( if ( !aPrevProfile.IsNull() ) { - if ( IsPointBetweenEdges( aPrevPnt1, aPrevPnt2, aPnt1, aPnt2, aTestPnt ) ) + if ( IsPointBetweenEdges2( aPrevPnt1, aPrevPnt2, aPnt1, aPnt2, aTestPnt ) ) { + DEBTRACE("point is between profiles "<< i-1 << " and " <GetShape3D() ); theRightWire = TopoDS::Wire( aProfile->GetShape3D() ); break; @@ -218,7 +265,7 @@ bool HYDROData_StreamAltitude::getBoundaryWiresForPoint( aPnt2.SetZ( 0 ); if ( !aPrevWire.IsNull() ) { - if ( IsPointBetweenEdges( aPrevPnt1, aPrevPnt2, aPnt1, aPnt2, aTestPnt ) ) { + if ( IsPointBetweenEdges2( aPrevPnt1, aPrevPnt2, aPnt1, aPnt2, aTestPnt ) ) { theLeftWire = aPrevWire; theRightWire = aWire; break; @@ -234,22 +281,32 @@ bool HYDROData_StreamAltitude::getBoundaryWiresForPoint( return !theLeftWire.IsNull() && !theRightWire.IsNull(); } -double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) const +double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint, + int theMethod) const { + DEBTRACE("GetAltitudeForPoint p(" << thePoint.X() << ", " << thePoint.Y() << ")"); double aResAltitude = GetInvalidAltitude(); Handle(HYDROData_Object) anObject = Handle(HYDROData_Object)::DownCast( GetFatherObject() ); if ( anObject.IsNull() ) + { + DEBTRACE(" --- aStream.IsNull()"); return aResAltitude; - + } TopoDS_Shape aTopShape = anObject->GetTopShape(); if ( aTopShape.IsNull() ) + { + DEBTRACE(" --- aTopShape.IsNull()"); return aResAltitude; + } TopExp_Explorer aFaceExp( aTopShape, TopAbs_FACE ); if ( !aFaceExp.More() ) + { + DEBTRACE(" --- !aFaceExp.More()"); return aResAltitude; + } // Get only face because of 2d profile wires is in compound TopoDS_Face aFace = TopoDS::Face( aFaceExp.Current() ); @@ -269,7 +326,10 @@ double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) co BRepTools::Write(aCmp, "FCL2d.brep"); #endif if ( aPointState == TopAbs_OUT ) + { + DEBTRACE(" --- aPointState == TopAbs_OUT"); return aResAltitude; + } TopoDS_Edge aLeftEdge, aRightEdge; @@ -292,7 +352,10 @@ double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) co // Find the two profiles between which the point is lies TopoDS_Wire aLeftWire, aRightWire; if ( !getBoundaryWiresForPoint( thePoint, aLeftWire, aRightWire ) ) + { + DEBTRACE(" --- !getBoundaryProfilesForPoint( thePoint, aLeftProfile, aRightProfile )"); return aResAltitude; + } // Find the projections of point to borders of stream gp_XYZ aPointToTest( thePoint.X(), thePoint.Y(), 0.0 ); @@ -307,7 +370,10 @@ double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) co int aNbPoints1 = aRightProject.NbPoints(); int aNbPoints2 = aLeftProject.NbPoints(); if ( aNbPoints1 < 1 || aNbPoints2 < 1) - return aResAltitude; + { + DEBTRACE(" --- projections " << aNbPoints1 << " " << aNbPoints2); + return aResAltitude; + } Standard_Real aLeftDist = aLeftProject.LowerDistance(); Standard_Real aRightDist = aRightProject.LowerDistance(); @@ -327,7 +393,7 @@ double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) co aRightProfileP2 = BRep_Tool::Pnt( aLastVertex ); // Interpolate altitudes - // Left profile line ( the segment between the firts and the last profile point ) + // Left profile line ( the segment between the first and the last profile point ) aLeftProfileP1.SetZ( 0 ); aLeftProfileP2.SetZ( 0 ); gp_Vec aLeftProfileVec( aLeftProfileP1, aLeftProfileP2 ); @@ -348,7 +414,7 @@ double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) co Standard_Real aSecCoeff = aRightProfileDist / ( aLeftProfileDist + aRightProfileDist ); aResAltitude = aLeftAlt * ( 1 - aFirstCoeff ) + aRightAlt * ( 1 - aSecCoeff ); - + DEBTRACE("aResAltitude=" << aResAltitude); return aResAltitude; }