#include <Geom_Line.hxx>
#include <QStringList>
+#include <cmath>
+
+#define _DEVDEBUG_
+#include "HYDRO_trace.hxx"
#ifdef DEB_CLASS2D
#include <BRepTools.hxx>
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 )
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 " <<i);
theLeftWire = TopoDS::Wire( aPrevProfile->GetShape3D() );
theRightWire = TopoDS::Wire( aProfile->GetShape3D() );
break;
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;
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() );
BRepTools::Write(aCmp, "FCL2d.brep");
#endif
if ( aPointState == TopAbs_OUT )
+ {
+ DEBTRACE(" --- aPointState == TopAbs_OUT");
return aResAltitude;
+ }
TopoDS_Edge aLeftEdge, aRightEdge;
// 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 );
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();
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 );
Standard_Real aSecCoeff = aRightProfileDist / ( aLeftProfileDist + aRightProfileDist );
aResAltitude = aLeftAlt * ( 1 - aFirstCoeff ) + aRightAlt * ( 1 - aSecCoeff );
-
+ DEBTRACE("aResAltitude=" << aResAltitude);
return aResAltitude;
}