X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FHYDROData%2FHYDROData_StreamAltitude.cxx;h=2360eddd0fce93aa735dd9b99d12bfe530995f6e;hb=e35b4caf3341d2cece2b7921f349acd17fb89351;hp=656bbdc1670580363f499ea5031078759bf1ffac;hpb=15327a2433ee39a181928d06a2acb5104bfbc979;p=modules%2Fhydro.git diff --git a/src/HYDROData/HYDROData_StreamAltitude.cxx b/src/HYDROData/HYDROData_StreamAltitude.cxx index 656bbdc1..2360eddd 100644 --- a/src/HYDROData/HYDROData_StreamAltitude.cxx +++ b/src/HYDROData/HYDROData_StreamAltitude.cxx @@ -1,15 +1,31 @@ +// 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 #include -#include - #include #include #include @@ -22,6 +38,7 @@ #include +#include #include #include @@ -29,25 +46,22 @@ #include #include +#include + +#include #include -#define CLASS2D 1 -#ifdef CLASS2D -#include -#include -#include -#include -#include -#include -#include -#undef _NCollection_MapHasher -#endif +#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() @@ -59,283 +73,348 @@ HYDROData_StreamAltitude::~HYDROData_StreamAltitude() { } -Standard_Real getAltitudeFromProfile( const Handle(HYDROData_Profile)& theProfile, - const Standard_Real& theLeftDist, - const Standard_Real& theRightDist ) -{ - Standard_Real aResAlt = 0.0; +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_XY aFirstPoint, aLastPoint; - if ( !theProfile->GetLeftPoint( aFirstPoint ) || - !theProfile->GetRightPoint( aLastPoint ) ) - return aResAlt; + BRepBuilderAPI_MakeWire aMakeWire( aLeftMakeEdge.Edge(), aBotMakeEdge.Edge(), + aRightMakeEdge.Edge(), aTopMakeEdge.Edge() ); - gp_Pnt aPnt1( aFirstPoint.X(), aFirstPoint.Y(), 0 ); - gp_Pnt aPnt2( aLastPoint.X(), aLastPoint.Y(), 0 ); - - 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 <=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 " <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 <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 <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; }