2 #include "HYDROData_StreamAltitude.h"
4 #include "HYDROData_Document.h"
5 #include "HYDROData_Profile.h"
6 #include "HYDROData_Stream.h"
7 #include "HYDROData_ShapesTool.h"
9 #include <BRep_Tool.hxx>
11 #include <BRepTopAdaptor_FClass2d.hxx>
13 #include <BRepBuilderAPI_MakeEdge.hxx>
14 #include <BRepBuilderAPI_MakeFace.hxx>
15 #include <BRepBuilderAPI_MakeWire.hxx>
17 #include <Extrema_ExtElC.hxx>
19 #include <GeomAPI_ProjectPointOnCurve.hxx>
23 #include <Precision.hxx>
25 #include <TopExp_Explorer.hxx>
28 #include <TopoDS_Wire.hxx>
29 #include <TopoDS_Edge.hxx>
31 #include <TopTools_SequenceOfShape.hxx>
33 #include <QStringList>
38 #include <Geom_Plane.hxx>
39 #include <BRepBuilderAPI_FindPlane.hxx>
40 #include <TopoDS_Vertex.hxx>
41 #include <TColgp_SequenceOfVec.hxx>
42 #include <TopTools_ShapeMapHasher.hxx>
43 #undef _NCollection_MapHasher
46 #include <BRepTools.hxx>
47 #include <BRep_Builder.hxx>
48 #include <BRepBuilderAPI_MakeVertex.hxx>
50 IMPLEMENT_STANDARD_HANDLE(HYDROData_StreamAltitude, HYDROData_IAltitudeObject)
51 IMPLEMENT_STANDARD_RTTIEXT(HYDROData_StreamAltitude, HYDROData_IAltitudeObject)
53 HYDROData_StreamAltitude::HYDROData_StreamAltitude()
54 : HYDROData_IAltitudeObject()
58 HYDROData_StreamAltitude::~HYDROData_StreamAltitude()
62 Standard_Real getAltitudeFromProfile( const Handle(HYDROData_Profile)& theProfile,
63 const Standard_Real& theLeftDist,
64 const Standard_Real& theRightDist )
66 Standard_Real aResAlt = 0.0;
68 gp_XY aFirstPoint, aLastPoint;
69 if ( !theProfile->GetLeftPoint( aFirstPoint ) ||
70 !theProfile->GetRightPoint( aLastPoint ) )
73 gp_Pnt aPnt1( aFirstPoint.X(), aFirstPoint.Y(), 0 );
74 gp_Pnt aPnt2( aLastPoint.X(), aLastPoint.Y(), 0 );
76 Standard_Real aProfileDist = aPnt1.Distance( aPnt2 );
78 Standard_Real aCoeff = aProfileDist / ( theLeftDist + theRightDist );
80 gp_Pnt anIntPoint( aPnt1.XYZ() + ( aCoeff * theLeftDist ) * gp_Dir( gp_Vec( aPnt1, aPnt2 ) ).XYZ() );
82 gp_Lin aPointLine( anIntPoint, gp::DZ() );
86 HYDROData_Profile::ProfilePoints aProfilePoints = theProfile->GetProfilePoints();
87 for ( int i = 1, n = aProfilePoints.Length(); i <= n; ++i )
89 gp_Pnt aProfPoint( aProfilePoints.Value( i ) );
91 Standard_Real aDist = aPointLine.Distance( aProfPoint );
92 if ( aDist <= gp::Resolution() )
94 // We found the intersected point
95 aResAlt = aProfPoint.Z();
99 gp_Lin aNormal = aPointLine.Normal( aProfPoint );
102 aPrevNormal = aNormal;
103 aPrevPoint = aProfPoint;
107 if ( aPrevNormal.Direction().Dot( aNormal.Direction() ) < 0 )
109 // We found the intersected edge
110 gp_Lin anEdgeLine( aPrevPoint, gp_Dir( gp_Vec( aPrevPoint, aProfPoint ) ) );
112 Extrema_ExtElC anExtrema( aPointLine, anEdgeLine, Precision::Angular() );
113 if ( !anExtrema.IsParallel() )
115 Extrema_POnCurv aFirstPnt, aSecPnt;
116 anExtrema.Points( 1, aFirstPnt, aSecPnt );
118 const gp_Pnt& anIntPnt = aSecPnt.Value();
119 aResAlt = anIntPnt.Z();
125 aPrevNormal = aNormal;
126 aPrevPoint = aProfPoint;
132 // =================================== workaround ==============================
133 // till fix of Classifier2d bug
135 TopAbs_State ComputePointState(const gp_XY& thePnt2d, const TopoDS_Face& theFace)
137 TopAbs_State aStatus(TopAbs_UNKNOWN);
140 TopoDS_Iterator it(theFace);
141 for(;it.More();it.Next()) {
142 aWire = TopoDS::Wire(it.Value());
145 if(nb > 1 || aWire.IsNull()) return aStatus;
148 BRepBuilderAPI_FindPlane fndPlane (theFace, Precision::Confusion());
150 aPlane = fndPlane.Plane()->Pln();
153 aNormal = gp_Vec(aPlane.Axis().Direction());
154 if(theFace.Orientation() == TopAbs_REVERSED)
156 gp_Pnt aPoint = gp_Pnt (thePnt2d.X(), thePnt2d.Y(), 0);
157 TColgp_SequenceOfVec aSeq;
158 TopTools_MapOfShape aMap;
159 it.Initialize(aWire);
160 for (;it.More(); it.Next()) {
161 const TopoDS_Vertex& aVx = TopExp::FirstVertex(TopoDS::Edge(it.Value()), Standard_True);
164 const gp_Pnt& aCurPnt = BRep_Tool::Pnt(TopoDS::Vertex(aVx));
165 if(aPoint.IsEqual(aCurPnt,Precision::Confusion())) {
169 gp_Vec aVec (aPoint, aCurPnt);
172 Standard_Real anAng(0.0);
173 for(int i = 1;i < aSeq.Length();i++) {
174 const gp_Vec& aV1 = aSeq.Value(i);
175 const gp_Vec& aV2 = aSeq.Value(i+1);
176 anAng += aV1.AngleWithRef(aV2, aNormal);
178 anAng += aSeq.Value(aSeq.Length()).AngleWithRef(aSeq.Value(1), aNormal);
179 if(abs(anAng) > Precision::Angular())
182 aStatus = TopAbs_OUT;
186 bool HYDROData_StreamAltitude::getBoundaryProfilesForPoint(
187 const gp_XY& thePoint,
188 Handle(HYDROData_Profile)& theLeftProfile,
189 Handle(HYDROData_Profile)& theRightProfile ) const
191 Handle(HYDROData_Stream) aStream =
192 Handle(HYDROData_Stream)::DownCast( GetFatherObject() );
193 if ( aStream.IsNull() )
196 HYDROData_SequenceOfObjects aStreamProfiles = aStream->GetProfiles();
197 if ( aStreamProfiles.Length() < 2 )
200 Handle(HYDROData_Profile) aPrevProfile;
201 gp_Pnt aPrevPnt1, aPrevPnt2;
202 for ( int i = 1, n = aStreamProfiles.Length(); i <= n; ++i )
204 Handle(HYDROData_Profile) aProfile =
205 Handle(HYDROData_Profile)::DownCast( aStreamProfiles.Value( i ) );
206 if ( aProfile.IsNull() )
209 gp_XY aFirstPoint, aLastPoint;
210 if ( !aProfile->GetLeftPoint( aFirstPoint ) || !aProfile->GetRightPoint( aLastPoint ) )
213 gp_Pnt aPnt1( aFirstPoint.X(), aFirstPoint.Y(), 0 );
214 gp_Pnt aPnt2( aLastPoint.X(), aLastPoint.Y(), 0 );
216 if ( !aPrevProfile.IsNull() )
218 BRepBuilderAPI_MakeEdge aLeftMakeEdge( aPrevPnt1, aPrevPnt2 );
219 BRepBuilderAPI_MakeEdge aBotMakeEdge( aPrevPnt2, aPnt2 );
220 BRepBuilderAPI_MakeEdge aRightMakeEdge( aPnt2, aPnt1 );
221 BRepBuilderAPI_MakeEdge aTopMakeEdge( aPnt1, aPrevPnt1 );
223 BRepBuilderAPI_MakeWire aMakeWire( aLeftMakeEdge.Edge(), aBotMakeEdge.Edge(),
224 aRightMakeEdge.Edge(), aTopMakeEdge.Edge() );
226 BRepBuilderAPI_MakeFace aMakeFace( aMakeWire.Wire() );
228 TopoDS_Face aProfilesFace = aMakeFace.Face();
230 TopoDS_Compound aCmp;
232 aBB.MakeCompound(aCmp);
233 aBB.Add(aCmp, aProfilesFace);
234 gp_Pnt aPnt (thePoint.X(), thePoint.Y(), 0.);
235 BRepBuilderAPI_MakeVertex aMk(aPnt);
236 aBB.Add(aCmp, aMk.Vertex());
237 BRepTools::Write(aCmp, "ProfileFace.brep");
239 BRepTopAdaptor_FClass2d aClassifier( aProfilesFace, Precision::Confusion() );
240 TopAbs_State aPointState =
242 ComputePointState(thePoint, aProfilesFace);
244 aClassifier.Perform( gp_Pnt2d( thePoint ), Standard_False );
247 cout << "Point status is = " << aPointState <<endl;
249 if ( aPointState != TopAbs_OUT )
251 theLeftProfile = aPrevProfile;
252 theRightProfile = aProfile;
257 aPrevProfile = aProfile;
262 return !theLeftProfile.IsNull() && !theRightProfile.IsNull();
265 double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) const
267 double aResAltitude = GetInvalidAltitude();
269 Handle(HYDROData_Stream) aStream =
270 Handle(HYDROData_Stream)::DownCast( GetFatherObject() );
271 if ( aStream.IsNull() )
274 TopoDS_Shape aStreamShape = aStream->GetTopShape();
275 if ( aStreamShape.IsNull() )
278 TopExp_Explorer aStreamFaceExp( aStreamShape, TopAbs_FACE );
279 if ( !aStreamFaceExp.More() )
282 // Get only face because of 2d profile wires is in compound
283 TopoDS_Face aStreamFace = TopoDS::Face( aStreamFaceExp.Current() );
285 // Check if point is inside of stream presentation
286 BRepTopAdaptor_FClass2d aClassifier( aStreamFace, Precision::Confusion() );
287 TopAbs_State aPointState =
289 ComputePointState(thePoint, aStreamFace);
291 aClassifier.Perform( gp_Pnt2d( thePoint ), Standard_False );
295 cout << "Point status is = " << aPointState <<endl;
296 TopoDS_Compound aCmp;
298 aBB.MakeCompound(aCmp);
299 aBB.Add(aCmp, aStreamFace);
300 gp_Pnt aPnt (thePoint.X(), thePoint.Y(), 0.);
301 BRepBuilderAPI_MakeVertex aMk(aPnt);
302 aBB.Add(aCmp, aMk.Vertex());
303 BRepTools::Write(aCmp, "FCL2d.brep");
305 if ( aPointState == TopAbs_OUT )
308 // Find the two profiles between which the point is lies
309 Handle(HYDROData_Profile) aLeftProfile, aRightProfile;
310 if ( !getBoundaryProfilesForPoint( thePoint, aLeftProfile, aRightProfile ) )
313 // Find the projections of point to borders of stream
314 gp_XYZ aPointToTest( thePoint.X(), thePoint.Y(), 0.0 );
316 TopoDS_Edge aStreamLeftEdge = TopoDS::Edge( aStream->GetLeftShape() );
317 TopoDS_Edge aStreamRightEdge = TopoDS::Edge( aStream->GetRightShape() );
319 Standard_Real aFirst = 0.0, aLast = 0.0;
320 Handle(Geom_Curve) anEdgeLeftCurve = BRep_Tool::Curve( aStreamLeftEdge, aFirst, aLast );
321 Handle(Geom_Curve) anEdgeRightCurve = BRep_Tool::Curve( aStreamRightEdge, aFirst, aLast );
323 GeomAPI_ProjectPointOnCurve aLeftProject( aPointToTest, anEdgeLeftCurve );
324 GeomAPI_ProjectPointOnCurve aRightProject( aPointToTest, anEdgeRightCurve );
326 Standard_Real aLeftDist = aLeftProject.LowerDistance();
327 Standard_Real aRightDist = aRightProject.LowerDistance();
329 // Find the altitude in profiles
330 Standard_Real aLeftAlt = getAltitudeFromProfile( aLeftProfile, aLeftDist, aRightDist );
331 Standard_Real aRightAlt = getAltitudeFromProfile( aRightProfile, aLeftDist, aRightDist );
333 // Interpolate altitudes
334 Standard_Real aFirstCoeff = aLeftDist / ( aLeftDist + aRightDist );
335 Standard_Real aSecCoeff = aRightDist / ( aLeftDist + aRightDist );
337 aResAltitude = aLeftAlt * aFirstCoeff + aRightAlt * aSecCoeff;