+// Copyright (C) 2014-2015 EDF-R&D
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
#include "HYDROData_StreamAltitude.h"
+#include "HYDROData_Channel.h"
#include "HYDROData_Document.h"
#include "HYDROData_Profile.h"
#include "HYDROData_Stream.h"
#include "HYDROData_ShapesTool.h"
-
+#include <HYDROData_Tool.h>
#include <BRep_Tool.hxx>
-#include <BRepTopAdaptor_FClass2d.hxx>
-
#include <BRepBuilderAPI_MakeEdge.hxx>
#include <BRepBuilderAPI_MakeFace.hxx>
#include <BRepBuilderAPI_MakeWire.hxx>
#include <Precision.hxx>
+#include <TopExp.hxx>
#include <TopExp_Explorer.hxx>
#include <TopoDS.hxx>
#include <TopoDS_Edge.hxx>
#include <TopTools_SequenceOfShape.hxx>
+#include <TopTools_ListIteratorOfListOfShape.hxx>
+
+#include <Geom_Line.hxx>
#include <QStringList>
-#define CLASS2D 1
-#ifdef CLASS2D
-#include <gp_Pln.hxx>
-#include <TopExp.hxx>
-#include <Geom_Plane.hxx>
-#include <BRepBuilderAPI_FindPlane.hxx>
-#include <TopoDS_Vertex.hxx>
-#include <TColgp_SequenceOfVec.hxx>
-#include <TopTools_ShapeMapHasher.hxx>
-#undef _NCollection_MapHasher
-#endif
+#include <cmath>
+
+#define _DEVDEBUG_
+#include "HYDRO_trace.hxx"
+
#ifdef DEB_CLASS2D
#include <BRepTools.hxx>
#include <BRep_Builder.hxx>
{
}
-Standard_Real getAltitudeFromProfile( const Handle(HYDROData_Profile)& theProfile,
- const Standard_Real& theLeftDist,
- const Standard_Real& theRightDist )
-{
- Standard_Real aResAlt = 0.0;
-
- gp_XY aFirstPoint, aLastPoint;
- if ( !theProfile->GetLeftPoint( aFirstPoint ) ||
- !theProfile->GetRightPoint( aLastPoint ) )
- return aResAlt;
+bool IsPointBetweenEdges( const gp_Pnt& aFirstPnt1, const gp_Pnt& aLastPnt1,
+ const gp_Pnt& aFirstPnt2, const gp_Pnt& aLastPnt2,
+ const gp_Pnt& thePoint) {
+ BRepBuilderAPI_MakeEdge aLeftMakeEdge( aFirstPnt1, aLastPnt1 );
+ BRepBuilderAPI_MakeEdge aBotMakeEdge( aLastPnt1, aLastPnt2 );
+ BRepBuilderAPI_MakeEdge aRightMakeEdge( aLastPnt2, aFirstPnt2 );
+ BRepBuilderAPI_MakeEdge aTopMakeEdge( aFirstPnt2, aFirstPnt1 );
- gp_Pnt aPnt1( aFirstPoint.X(), aFirstPoint.Y(), 0 );
- gp_Pnt aPnt2( aLastPoint.X(), aLastPoint.Y(), 0 );
+ BRepBuilderAPI_MakeWire aMakeWire( aLeftMakeEdge.Edge(), aBotMakeEdge.Edge(),
+ aRightMakeEdge.Edge(), aTopMakeEdge.Edge() );
- Standard_Real aProfileDist = aPnt1.Distance( aPnt2 );
-
- Standard_Real aCoeff = aProfileDist / ( theLeftDist + theRightDist );
-
- gp_Pnt anIntPoint( aPnt1.XYZ() + ( aCoeff * theLeftDist ) * gp_Dir( gp_Vec( aPnt1, aPnt2 ) ).XYZ() );
-
- gp_Lin aPointLine( anIntPoint, gp::DZ() );
+ BRepBuilderAPI_MakeFace aMakeFace( aMakeWire.Wire() );
+
+ TopoDS_Face aFace = aMakeFace.Face();
+#ifdef DEB_CLASS2D
+ TopoDS_Compound aCmp;
+ BRep_Builder aBB;
+ aBB.MakeCompound(aCmp);
+ aBB.Add(aCmp, aFace);
+ BRepBuilderAPI_MakeVertex aMk(thePoint);
+ aBB.Add(aCmp, aMk.Vertex());
+ BRepTools::Write(aCmp, "ProfileFace.brep");
+#endif
+
+ gp_XY anXY( thePoint.X(), thePoint.Y() );
+ TopAbs_State aPointState = HYDROData_Tool::ComputePointState(anXY, aFace);
- gp_Pnt aPrevPoint;
- gp_Lin aPrevNormal;
- HYDROData_Profile::ProfilePoints aProfilePoints = theProfile->GetProfilePoints();
- for ( int i = 1, n = aProfilePoints.Length(); i <= n; ++i )
- {
- gp_Pnt aProfPoint( aProfilePoints.Value( i ) );
+#ifdef DEB_CLASS2D
+ cout << "Point status is = " << aPointState <<endl;
+#endif
+ return aPointState != TopAbs_OUT;
+}
- Standard_Real aDist = aPointLine.Distance( aProfPoint );
- if ( aDist <= gp::Resolution() )
+/** 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)))
{
- // We found the intersected point
- aResAlt = aProfPoint.Z();
- break;
+ isBetween = false;
}
-
- gp_Lin aNormal = aPointLine.Normal( aProfPoint );
- if ( i == 1 )
+ if (isBetween)
{
- aPrevNormal = aNormal;
- aPrevPoint = aProfPoint;
- continue;
+ 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;
+}
- if ( aPrevNormal.Direction().Dot( aNormal.Direction() ) < 0 )
- {
- // We found the intersected edge
- gp_Lin anEdgeLine( aPrevPoint, gp_Dir( gp_Vec( aPrevPoint, aProfPoint ) ) );
-
- Extrema_ExtElC anExtrema( aPointLine, anEdgeLine, Precision::Angular() );
- if ( !anExtrema.IsParallel() )
- {
- Extrema_POnCurv aFirstPnt, aSecPnt;
- anExtrema.Points( 1, aFirstPnt, aSecPnt );
-
- const gp_Pnt& anIntPnt = aSecPnt.Value();
- aResAlt = anIntPnt.Z();
-
- break;
- }
- }
+Standard_Real getAltitudeFromWire( const TopoDS_Wire& theWire,
+ const Standard_Real& theLeftDist,
+ const Standard_Real& theRightDist )
+{
+ Standard_Real aResAlt = 0.0;
- aPrevNormal = aNormal;
- aPrevPoint = aProfPoint;
- }
+ gp_XY aFirstPoint, aLastPoint;
+ TopoDS_Vertex aFirstVertex, aLastVertex;
+ TopExp::Vertices( theWire, aFirstVertex, aLastVertex );
+
+ gp_Pnt aPnt1( BRep_Tool::Pnt( aFirstVertex ) );
+ aPnt1.SetZ( 0 );
+ gp_Pnt aPnt2( BRep_Tool::Pnt( aLastVertex ) );
+ aPnt2.SetZ( 0 );
+
+ Standard_Real aProfileDist = aPnt1.Distance( aPnt2 );
- return aResAlt;
-}
+ Standard_Real aCoeff = aProfileDist / ( theLeftDist + theRightDist );
-// =================================== workaround ==============================
-// till fix of Classifier2d bug
-#ifdef CLASS2D
-TopAbs_State ComputePointState(const gp_XY& thePnt2d, const TopoDS_Face& theFace)
-{
- TopAbs_State aStatus(TopAbs_UNKNOWN);
- TopoDS_Wire aWire;
- int nb(0);
- TopoDS_Iterator it(theFace);
- for(;it.More();it.Next()) {
- aWire = TopoDS::Wire(it.Value());
- nb++;
- }
- if(nb > 1 || aWire.IsNull()) return aStatus;
- gp_Pln aPlane;
- gp_Vec aNormal;
- BRepBuilderAPI_FindPlane fndPlane (theFace, Precision::Confusion());
- if(fndPlane.Found())
- aPlane = fndPlane.Plane()->Pln();
- else
- return aStatus;
- aNormal = gp_Vec(aPlane.Axis().Direction());
- if(theFace.Orientation() == TopAbs_REVERSED)
- aNormal.Reverse();
- gp_Pnt aPoint = gp_Pnt (thePnt2d.X(), thePnt2d.Y(), 0);
- TColgp_SequenceOfVec aSeq;
- TopTools_MapOfShape aMap;
- it.Initialize(aWire);
- for (;it.More(); it.Next()) {
- const TopoDS_Vertex& aVx = TopExp::FirstVertex(TopoDS::Edge(it.Value()), Standard_True);
- if(!aMap.Add(aVx))
- continue;
- const gp_Pnt& aCurPnt = BRep_Tool::Pnt(TopoDS::Vertex(aVx));
- if(aPoint.IsEqual(aCurPnt,Precision::Confusion())) {
- aStatus = TopAbs_ON;
- return aStatus;
- }
- gp_Vec aVec (aPoint, aCurPnt);
- aSeq.Append(aVec);
- }
- Standard_Real anAng(0.0);
- for(int i = 1;i < aSeq.Length();i++) {
- const gp_Vec& aV1 = aSeq.Value(i);
- const gp_Vec& aV2 = aSeq.Value(i+1);
- anAng += aV1.AngleWithRef(aV2, aNormal);
- }
- anAng += aSeq.Value(aSeq.Length()).AngleWithRef(aSeq.Value(1), aNormal);
- if(abs(anAng) > Precision::Angular())
- aStatus = TopAbs_IN;
- else
- aStatus = TopAbs_OUT;
- return aStatus;
+ gp_Pnt anIntPoint( aPnt1.XYZ() + ( aCoeff * theLeftDist ) * gp_Dir( gp_Vec( aPnt1, aPnt2 ) ).XYZ() );
+
+ return HYDROData_Tool::GetAltitudeForWire( theWire,
+ gp_XY(anIntPoint.X(), anIntPoint.Y()),
+ 1E-2,
+ 1E-2,
+ HYDROData_IAltitudeObject::GetInvalidAltitude() );
}
-#endif
-bool HYDROData_StreamAltitude::getBoundaryProfilesForPoint(
- const gp_XY& thePoint,
- Handle(HYDROData_Profile)& theLeftProfile,
- Handle(HYDROData_Profile)& theRightProfile ) const
+
+bool HYDROData_StreamAltitude::getBoundaryWiresForPoint(
+ const gp_XY& thePoint,
+ TopoDS_Wire& theLeftWire,
+ TopoDS_Wire& theRightWire ) const
{
- Handle(HYDROData_Stream) aStream =
- Handle(HYDROData_Stream)::DownCast( GetFatherObject() );
- if ( aStream.IsNull() )
- return false;
+ gp_Pnt aTestPnt( thePoint.X(), thePoint.Y(), 0 );
- HYDROData_SequenceOfObjects aStreamProfiles = aStream->GetProfiles();
- if ( aStreamProfiles.Length() < 2 )
+ Handle(HYDROData_Object) anObject =
+ Handle(HYDROData_Object)::DownCast( GetFatherObject() );
+ if ( anObject.IsNull() ) {
return false;
+ }
- Handle(HYDROData_Profile) aPrevProfile;
- gp_Pnt aPrevPnt1, aPrevPnt2;
- for ( int i = 1, n = aStreamProfiles.Length(); i <= n; ++i )
- {
- Handle(HYDROData_Profile) aProfile =
- Handle(HYDROData_Profile)::DownCast( aStreamProfiles.Value( i ) );
- if ( aProfile.IsNull() )
- continue;
-
- gp_XY aFirstPoint, aLastPoint;
- if ( !aProfile->GetLeftPoint( aFirstPoint ) || !aProfile->GetRightPoint( aLastPoint ) )
- continue;
+ if ( anObject->GetKind() == KIND_STREAM ) {
+ Handle(HYDROData_Stream) aStream = Handle(HYDROData_Stream)::DownCast( anObject );
+ if ( aStream.IsNull() )
+ return false;
- gp_Pnt aPnt1( aFirstPoint.X(), aFirstPoint.Y(), 0 );
- gp_Pnt aPnt2( aLastPoint.X(), aLastPoint.Y(), 0 );
+ HYDROData_SequenceOfObjects aStreamProfiles = aStream->GetProfiles();
+ if ( aStreamProfiles.Length() < 2 )
+ return false;
- if ( !aPrevProfile.IsNull() )
+ Handle(HYDROData_Profile) aPrevProfile;
+ gp_Pnt aPrevPnt1, aPrevPnt2;
+ for ( int i = 1, n = aStreamProfiles.Length(); i <= n; ++i )
{
- BRepBuilderAPI_MakeEdge aLeftMakeEdge( aPrevPnt1, aPrevPnt2 );
- BRepBuilderAPI_MakeEdge aBotMakeEdge( aPrevPnt2, aPnt2 );
- BRepBuilderAPI_MakeEdge aRightMakeEdge( aPnt2, aPnt1 );
- BRepBuilderAPI_MakeEdge aTopMakeEdge( aPnt1, aPrevPnt1 );
-
- BRepBuilderAPI_MakeWire aMakeWire( aLeftMakeEdge.Edge(), aBotMakeEdge.Edge(),
- aRightMakeEdge.Edge(), aTopMakeEdge.Edge() );
+ Handle(HYDROData_Profile) aProfile =
+ Handle(HYDROData_Profile)::DownCast( aStreamProfiles.Value( i ) );
+ if ( aProfile.IsNull() )
+ continue;
+
+ gp_XY aFirstPoint, aLastPoint;
+ if ( !aProfile->GetLeftPoint( aFirstPoint, false ) ||
+ !aProfile->GetRightPoint( aLastPoint, false ) )
+ continue;
+
+ gp_Pnt aPnt1( aFirstPoint.X(), aFirstPoint.Y(), 0 );
+ gp_Pnt aPnt2( aLastPoint.X(), aLastPoint.Y(), 0 );
+
+ if ( !aPrevProfile.IsNull() )
+ {
+ 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;
+ }
+ }
- BRepBuilderAPI_MakeFace aMakeFace( aMakeWire.Wire() );
-
- TopoDS_Face aProfilesFace = aMakeFace.Face();
-#ifdef DEB_CLASS2D
- TopoDS_Compound aCmp;
- BRep_Builder aBB;
- aBB.MakeCompound(aCmp);
- aBB.Add(aCmp, aProfilesFace);
- gp_Pnt aPnt (thePoint.X(), thePoint.Y(), 0.);
- BRepBuilderAPI_MakeVertex aMk(aPnt);
- aBB.Add(aCmp, aMk.Vertex());
- BRepTools::Write(aCmp, "ProfileFace.brep");
-#endif
- BRepTopAdaptor_FClass2d aClassifier( aProfilesFace, Precision::Confusion() );
- TopAbs_State aPointState =
-#ifdef CLASS2D
- ComputePointState(thePoint, aProfilesFace);
-#else
- aClassifier.Perform( gp_Pnt2d( thePoint ), Standard_False );
-#endif
-#ifdef DEB_CLASS2D
- cout << "Point status is = " << aPointState <<endl;
-#endif
- if ( aPointState != TopAbs_OUT )
- {
- theLeftProfile = aPrevProfile;
- theRightProfile = aProfile;
- break;
+ aPrevProfile = aProfile;
+ aPrevPnt1 = aPnt1;
+ aPrevPnt2 = aPnt2;
+ }
+ } else if ( anObject->GetKind() == KIND_CHANNEL ) {
+ Handle(HYDROData_Channel) aChannel = Handle(HYDROData_Channel)::DownCast( anObject );
+ if ( aChannel.IsNull() )
+ return false;
+
+ TopTools_ListOfShape aWiresList;
+ TopExp_Explorer anExp( aChannel->GetShape3D(), TopAbs_WIRE );
+ for ( ; anExp.More(); anExp.Next() ) {
+ if(!anExp.Current().IsNull()) {
+ const TopoDS_Wire& aWire = TopoDS::Wire( anExp.Current() );
+ aWiresList.Append( aWire );
}
}
- aPrevProfile = aProfile;
- aPrevPnt1 = aPnt1;
- aPrevPnt2 = aPnt2;
+ if ( aWiresList.Extent() < 2 )
+ return false;
+
+ TopoDS_Wire aPrevWire;
+ gp_Pnt aPrevPnt1, aPrevPnt2;
+
+ TopTools_ListIteratorOfListOfShape anIt( aWiresList );
+ for ( ; anIt.More(); anIt.Next() ) {
+ TopoDS_Wire& aWire = TopoDS::Wire( anIt.Value() );
+ if ( aWire.IsNull() )
+ continue;
+
+ TopoDS_Vertex aFirstVertex, aLastVertex;
+ TopExp::Vertices( aWire, aFirstVertex, aLastVertex );
+
+ gp_Pnt aPnt1( BRep_Tool::Pnt( aFirstVertex ) );
+ aPnt1.SetZ( 0 );
+ gp_Pnt aPnt2( BRep_Tool::Pnt( aLastVertex ) );
+ aPnt2.SetZ( 0 );
+
+ if ( !aPrevWire.IsNull() ) {
+ if ( IsPointBetweenEdges2( aPrevPnt1, aPrevPnt2, aPnt1, aPnt2, aTestPnt ) ) {
+ theLeftWire = aPrevWire;
+ theRightWire = aWire;
+ break;
+ }
+ }
+
+ aPrevWire = aWire;
+ aPrevPnt1 = aPnt1;
+ aPrevPnt2 = aPnt2;
+ }
}
- return !theLeftProfile.IsNull() && !theRightProfile.IsNull();
+ 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_Stream) aStream =
- Handle(HYDROData_Stream)::DownCast( GetFatherObject() );
- if ( aStream.IsNull() )
+ Handle(HYDROData_Object) anObject =
+ Handle(HYDROData_Object)::DownCast( GetFatherObject() );
+ if ( anObject.IsNull() )
+ {
+ DEBTRACE(" --- aStream.IsNull()");
return aResAltitude;
-
- TopoDS_Shape aStreamShape = aStream->GetTopShape();
- if ( aStreamShape.IsNull() )
+ }
+ TopoDS_Shape aTopShape = anObject->GetTopShape();
+ if ( aTopShape.IsNull() )
+ {
+ DEBTRACE(" --- aTopShape.IsNull()");
return aResAltitude;
+ }
- TopExp_Explorer aStreamFaceExp( aStreamShape, TopAbs_FACE );
- if ( !aStreamFaceExp.More() )
+ 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 aStreamFace = TopoDS::Face( aStreamFaceExp.Current() );
-
- // Check if point is inside of stream presentation
- BRepTopAdaptor_FClass2d aClassifier( aStreamFace, Precision::Confusion() );
- TopAbs_State aPointState =
-#ifdef CLASS2D
- ComputePointState(thePoint, aStreamFace);
-#else
- aClassifier.Perform( gp_Pnt2d( thePoint ), Standard_False );
-#endif
+ TopoDS_Face aFace = TopoDS::Face( aFaceExp.Current() );
+
+ // Check if point is inside of stream/channel presentation
+ TopAbs_State aPointState = HYDROData_Tool::ComputePointState(thePoint, aFace);
#ifdef DEB_CLASS2D
cout << "Point status is = " << aPointState <<endl;
TopoDS_Compound aCmp;
BRep_Builder aBB;
aBB.MakeCompound(aCmp);
- aBB.Add(aCmp, aStreamFace);
+ aBB.Add(aCmp, aFace);
gp_Pnt aPnt (thePoint.X(), thePoint.Y(), 0.);
BRepBuilderAPI_MakeVertex aMk(aPnt);
aBB.Add(aCmp, aMk.Vertex());
BRepTools::Write(aCmp, "FCL2d.brep");
#endif
if ( aPointState == TopAbs_OUT )
+ {
+ DEBTRACE(" --- aPointState == TopAbs_OUT");
return aResAltitude;
+ }
+ TopoDS_Edge aLeftEdge, aRightEdge;
+
+ if ( anObject->GetKind() == KIND_STREAM ) {
+ Handle(HYDROData_Stream) aStream =
+ Handle(HYDROData_Stream)::DownCast( GetFatherObject() );
+ if ( !aStream.IsNull() ) {
+ aLeftEdge = TopoDS::Edge( aStream->GetLeftShape() );
+ aRightEdge = TopoDS::Edge( aStream->GetRightShape() );
+ }
+ } else if ( anObject->GetKind() == KIND_CHANNEL ) {
+ Handle(HYDROData_Channel) aChannel =
+ Handle(HYDROData_Channel)::DownCast( GetFatherObject() );
+ if ( !aChannel.IsNull() ) {
+ aLeftEdge = TopoDS::Edge( aChannel->GetLeftShape() );
+ aRightEdge = TopoDS::Edge( aChannel->GetRightShape() );
+ }
+ }
+
// Find the two profiles between which the point is lies
- Handle(HYDROData_Profile) aLeftProfile, aRightProfile;
- if ( !getBoundaryProfilesForPoint( thePoint, aLeftProfile, aRightProfile ) )
- return aResAltitude;
+ 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 );
-
- TopoDS_Edge aStreamLeftEdge = TopoDS::Edge( aStream->GetLeftShape() );
- TopoDS_Edge aStreamRightEdge = TopoDS::Edge( aStream->GetRightShape() );
-
+
Standard_Real aFirst = 0.0, aLast = 0.0;
- Handle(Geom_Curve) anEdgeLeftCurve = BRep_Tool::Curve( aStreamLeftEdge, aFirst, aLast );
- Handle(Geom_Curve) anEdgeRightCurve = BRep_Tool::Curve( aStreamRightEdge, aFirst, aLast );
+ Handle(Geom_Curve) anEdgeLeftCurve = BRep_Tool::Curve( aLeftEdge, aFirst, aLast );
+ Handle(Geom_Curve) anEdgeRightCurve = BRep_Tool::Curve( aRightEdge, aFirst, aLast );
GeomAPI_ProjectPointOnCurve aLeftProject( aPointToTest, anEdgeLeftCurve );
GeomAPI_ProjectPointOnCurve aRightProject( aPointToTest, anEdgeRightCurve );
+ int aNbPoints1 = aRightProject.NbPoints();
+ int aNbPoints2 = aLeftProject.NbPoints();
+ if ( aNbPoints1 < 1 || aNbPoints2 < 1)
+ {
+ DEBTRACE(" --- projections " << aNbPoints1 << " " << aNbPoints2);
+ return aResAltitude;
+ }
Standard_Real aLeftDist = aLeftProject.LowerDistance();
Standard_Real aRightDist = aRightProject.LowerDistance();
// Find the altitude in profiles
- Standard_Real aLeftAlt = getAltitudeFromProfile( aLeftProfile, aLeftDist, aRightDist );
- Standard_Real aRightAlt = getAltitudeFromProfile( aRightProfile, aLeftDist, aRightDist );
+ Standard_Real aLeftAlt, aRightAlt;
+ gp_Pnt aLeftProfileP1, aLeftProfileP2, aRightProfileP1, aRightProfileP2;
+ aLeftAlt = getAltitudeFromWire( aLeftWire, aLeftDist, aRightDist );
+ aRightAlt = getAltitudeFromWire( aRightWire, aLeftDist, aRightDist );
+
+ TopoDS_Vertex aFirstVertex, aLastVertex;
+ TopExp::Vertices( aLeftWire, aFirstVertex, aLastVertex );
+ aLeftProfileP1 = BRep_Tool::Pnt( aFirstVertex );
+ aLeftProfileP2 = BRep_Tool::Pnt( aLastVertex );
+
+ TopExp::Vertices( aRightWire, aFirstVertex, aLastVertex );
+ aRightProfileP1 = BRep_Tool::Pnt( aFirstVertex );
+ aRightProfileP2 = BRep_Tool::Pnt( aLastVertex );
// Interpolate altitudes
- Standard_Real aFirstCoeff = aLeftDist / ( aLeftDist + aRightDist );
- Standard_Real aSecCoeff = aRightDist / ( aLeftDist + aRightDist );
-
- aResAltitude = aLeftAlt * aFirstCoeff + aRightAlt * aSecCoeff;
-
+ // Left profile line ( the segment between the first and the last profile point )
+ aLeftProfileP1.SetZ( 0 );
+ aLeftProfileP2.SetZ( 0 );
+ gp_Vec aLeftProfileVec( aLeftProfileP1, aLeftProfileP2 );
+ Handle(Geom_Line) aLeftProfileLine = new Geom_Line( gp_Ax1( aLeftProfileP1, aLeftProfileVec ) );
+ // Right profile line
+ aRightProfileP1.SetZ( 0 );
+ aRightProfileP2.SetZ( 0 );
+ gp_Vec aRightProfileVec( aRightProfileP1, aRightProfileP2 );
+ Handle(Geom_Line) aRightProfileLine = new Geom_Line( gp_Ax1( aRightProfileP1, aRightProfileVec ) );
+ // The point projections on the left and right profiles
+ GeomAPI_ProjectPointOnCurve aLeftProfileProject( aPointToTest, aLeftProfileLine );
+ GeomAPI_ProjectPointOnCurve aRightProfileProject( aPointToTest, aRightProfileLine );
+ // The point distance to the left and right profiles
+ Standard_Real aLeftProfileDist = aLeftProfileProject.LowerDistance();
+ Standard_Real aRightProfileDist = aRightProfileProject.LowerDistance();
+ // The coefficients
+ Standard_Real aFirstCoeff = aLeftProfileDist / ( aLeftProfileDist + aRightProfileDist );
+ Standard_Real aSecCoeff = aRightProfileDist / ( aLeftProfileDist + aRightProfileDist );
+
+ aResAltitude = aLeftAlt * ( 1 - aFirstCoeff ) + aRightAlt * ( 1 - aSecCoeff );
+ DEBTRACE("aResAltitude=" << aResAltitude);
return aResAltitude;
}