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 QStringList HYDROData_StreamAltitude::DumpToPython( MapOfTreatedObjects& theTreatedObjects ) const
64 QStringList aResList = dumpObjectCreation( theTreatedObjects );
65 QString aName = GetObjPyName();
69 aResList << QString( "" );
70 aResList << QString( "%1.Update();" ).arg( aName );
71 aResList << QString( "" );
76 Standard_Real getAltitudeFromProfile( const Handle(HYDROData_Profile)& theProfile,
77 const Standard_Real& theLeftDist,
78 const Standard_Real& theRightDist )
80 Standard_Real aResAlt = 0.0;
82 gp_XY aFirstPoint, aLastPoint;
83 if ( !theProfile->GetLeftPoint( aFirstPoint ) ||
84 !theProfile->GetRightPoint( aLastPoint ) )
87 gp_Pnt aPnt1( aFirstPoint.X(), aFirstPoint.Y(), 0 );
88 gp_Pnt aPnt2( aLastPoint.X(), aLastPoint.Y(), 0 );
90 Standard_Real aProfileDist = aPnt1.Distance( aPnt2 );
92 Standard_Real aCoeff = aProfileDist / ( theLeftDist + theRightDist );
94 gp_Pnt anIntPoint( aPnt1.XYZ() + ( aCoeff * theLeftDist ) * gp_Dir( gp_Vec( aPnt1, aPnt2 ) ).XYZ() );
96 gp_Lin aPointLine( anIntPoint, gp::DZ() );
100 HYDROData_Profile::ProfilePoints aProfilePoints = theProfile->GetProfilePoints();
101 for ( int i = 1, n = aProfilePoints.Length(); i <= n; ++i )
103 gp_Pnt aProfPoint( aProfilePoints.Value( i ) );
105 Standard_Real aDist = aPointLine.Distance( aProfPoint );
106 if ( aDist <= gp::Resolution() )
108 // We found the intersected point
109 aResAlt = aProfPoint.Z();
113 gp_Lin aNormal = aPointLine.Normal( aProfPoint );
116 aPrevNormal = aNormal;
117 aPrevPoint = aProfPoint;
121 if ( aPrevNormal.Direction().Dot( aNormal.Direction() ) < 0 )
123 // We found the intersected edge
124 gp_Lin anEdgeLine( aPrevPoint, gp_Dir( gp_Vec( aPrevPoint, aProfPoint ) ) );
126 Extrema_ExtElC anExtrema( aPointLine, anEdgeLine, Precision::Angular() );
127 if ( !anExtrema.IsParallel() )
129 Extrema_POnCurv aFirstPnt, aSecPnt;
130 anExtrema.Points( 1, aFirstPnt, aSecPnt );
132 const gp_Pnt& anIntPnt = aSecPnt.Value();
133 aResAlt = anIntPnt.Z();
139 aPrevNormal = aNormal;
140 aPrevPoint = aProfPoint;
146 // =================================== workaround ==============================
147 // till fix of Classifier2d bug
149 TopAbs_State ComputePointState(const gp_XY& thePnt2d, const TopoDS_Face& theFace)
151 TopAbs_State aStatus(TopAbs_UNKNOWN);
154 TopoDS_Iterator it(theFace);
155 for(;it.More();it.Next()) {
156 aWire = TopoDS::Wire(it.Value());
159 if(nb > 1 || aWire.IsNull()) return aStatus;
162 BRepBuilderAPI_FindPlane fndPlane (theFace, Precision::Confusion());
164 aPlane = fndPlane.Plane()->Pln();
167 aNormal = gp_Vec(aPlane.Axis().Direction());
168 if(theFace.Orientation() == TopAbs_REVERSED)
170 gp_Pnt aPoint = gp_Pnt (thePnt2d.X(), thePnt2d.Y(), 0);
171 TColgp_SequenceOfVec aSeq;
172 TopTools_MapOfShape aMap;
173 it.Initialize(aWire);
174 for (;it.More(); it.Next()) {
175 const TopoDS_Vertex& aVx = TopExp::FirstVertex(TopoDS::Edge(it.Value()), Standard_True);
178 const gp_Pnt& aCurPnt = BRep_Tool::Pnt(TopoDS::Vertex(aVx));
179 if(aPoint.IsEqual(aCurPnt,Precision::Confusion())) {
183 gp_Vec aVec (aPoint, aCurPnt);
186 Standard_Real anAng(0.0);
187 for(int i = 1;i < aSeq.Length();i++) {
188 const gp_Vec& aV1 = aSeq.Value(i);
189 const gp_Vec& aV2 = aSeq.Value(i+1);
190 anAng += aV1.AngleWithRef(aV2, aNormal);
192 anAng += aSeq.Value(aSeq.Length()).AngleWithRef(aSeq.Value(1), aNormal);
193 if(abs(anAng) > Precision::Angular())
196 aStatus = TopAbs_OUT;
200 bool HYDROData_StreamAltitude::getBoundaryProfilesForPoint(
201 const gp_XY& thePoint,
202 Handle(HYDROData_Profile)& theLeftProfile,
203 Handle(HYDROData_Profile)& theRightProfile ) const
205 Handle(HYDROData_Stream) aStream =
206 Handle(HYDROData_Stream)::DownCast( GetFatherObject() );
207 if ( aStream.IsNull() )
210 HYDROData_SequenceOfObjects aStreamProfiles = aStream->GetProfiles();
211 if ( aStreamProfiles.Length() < 2 )
214 Handle(HYDROData_Profile) aPrevProfile;
215 gp_Pnt aPrevPnt1, aPrevPnt2;
216 for ( int i = 1, n = aStreamProfiles.Length(); i <= n; ++i )
218 Handle(HYDROData_Profile) aProfile =
219 Handle(HYDROData_Profile)::DownCast( aStreamProfiles.Value( i ) );
220 if ( aProfile.IsNull() )
223 gp_XY aFirstPoint, aLastPoint;
224 if ( !aProfile->GetLeftPoint( aFirstPoint ) || !aProfile->GetRightPoint( aLastPoint ) )
227 gp_Pnt aPnt1( aFirstPoint.X(), aFirstPoint.Y(), 0 );
228 gp_Pnt aPnt2( aLastPoint.X(), aLastPoint.Y(), 0 );
230 if ( !aPrevProfile.IsNull() )
232 BRepBuilderAPI_MakeEdge aLeftMakeEdge( aPrevPnt1, aPrevPnt2 );
233 BRepBuilderAPI_MakeEdge aBotMakeEdge( aPrevPnt2, aPnt2 );
234 BRepBuilderAPI_MakeEdge aRightMakeEdge( aPnt2, aPnt1 );
235 BRepBuilderAPI_MakeEdge aTopMakeEdge( aPnt1, aPrevPnt1 );
237 BRepBuilderAPI_MakeWire aMakeWire( aLeftMakeEdge.Edge(), aBotMakeEdge.Edge(),
238 aRightMakeEdge.Edge(), aTopMakeEdge.Edge() );
240 BRepBuilderAPI_MakeFace aMakeFace( aMakeWire.Wire() );
242 TopoDS_Face aProfilesFace = aMakeFace.Face();
244 TopoDS_Compound aCmp;
246 aBB.MakeCompound(aCmp);
247 aBB.Add(aCmp, aProfilesFace);
248 gp_Pnt aPnt (thePoint.X(), thePoint.Y(), 0.);
249 BRepBuilderAPI_MakeVertex aMk(aPnt);
250 aBB.Add(aCmp, aMk.Vertex());
251 BRepTools::Write(aCmp, "ProfileFace.brep");
253 BRepTopAdaptor_FClass2d aClassifier( aProfilesFace, Precision::Confusion() );
254 TopAbs_State aPointState =
256 ComputePointState(thePoint, aProfilesFace);
258 aClassifier.Perform( gp_Pnt2d( thePoint ), Standard_False );
261 cout << "Point status is = " << aPointState <<endl;
263 if ( aPointState != TopAbs_OUT )
265 theLeftProfile = aPrevProfile;
266 theRightProfile = aProfile;
271 aPrevProfile = aProfile;
276 return !theLeftProfile.IsNull() && !theRightProfile.IsNull();
279 double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) const
281 double aResAltitude = GetInvalidAltitude();
283 Handle(HYDROData_Stream) aStream =
284 Handle(HYDROData_Stream)::DownCast( GetFatherObject() );
285 if ( aStream.IsNull() )
288 TopoDS_Shape aStreamShape = aStream->GetTopShape();
289 if ( aStreamShape.IsNull() )
292 TopExp_Explorer aStreamFaceExp( aStreamShape, TopAbs_FACE );
293 if ( !aStreamFaceExp.More() )
296 // Get only face because of 2d profile wires is in compound
297 TopoDS_Face aStreamFace = TopoDS::Face( aStreamFaceExp.Current() );
299 // Check if point is inside of stream presentation
300 BRepTopAdaptor_FClass2d aClassifier( aStreamFace, Precision::Confusion() );
301 TopAbs_State aPointState =
303 ComputePointState(thePoint, aStreamFace);
305 aClassifier.Perform( gp_Pnt2d( thePoint ), Standard_False );
309 cout << "Point status is = " << aPointState <<endl;
310 TopoDS_Compound aCmp;
312 aBB.MakeCompound(aCmp);
313 aBB.Add(aCmp, aStreamFace);
314 gp_Pnt aPnt (thePoint.X(), thePoint.Y(), 0.);
315 BRepBuilderAPI_MakeVertex aMk(aPnt);
316 aBB.Add(aCmp, aMk.Vertex());
317 BRepTools::Write(aCmp, "FCL2d.brep");
319 if ( aPointState == TopAbs_OUT )
322 // Find the two profiles between which the point is lies
323 Handle(HYDROData_Profile) aLeftProfile, aRightProfile;
324 if ( !getBoundaryProfilesForPoint( thePoint, aLeftProfile, aRightProfile ) )
327 // Find the projections of point to borders of stream
328 gp_XYZ aPointToTest( thePoint.X(), thePoint.Y(), 0.0 );
330 TopoDS_Edge aStreamLeftEdge = TopoDS::Edge( aStream->GetLeftShape() );
331 TopoDS_Edge aStreamRightEdge = TopoDS::Edge( aStream->GetRightShape() );
333 Standard_Real aFirst = 0.0, aLast = 0.0;
334 Handle(Geom_Curve) anEdgeLeftCurve = BRep_Tool::Curve( aStreamLeftEdge, aFirst, aLast );
335 Handle(Geom_Curve) anEdgeRightCurve = BRep_Tool::Curve( aStreamRightEdge, aFirst, aLast );
337 GeomAPI_ProjectPointOnCurve aLeftProject( aPointToTest, anEdgeLeftCurve );
338 GeomAPI_ProjectPointOnCurve aRightProject( aPointToTest, anEdgeRightCurve );
340 Standard_Real aLeftDist = aLeftProject.LowerDistance();
341 Standard_Real aRightDist = aRightProject.LowerDistance();
343 // Find the altitude in profiles
344 Standard_Real aLeftAlt = getAltitudeFromProfile( aLeftProfile, aLeftDist, aRightDist );
345 Standard_Real aRightAlt = getAltitudeFromProfile( aRightProfile, aLeftDist, aRightDist );
347 // Interpolate altitudes
348 Standard_Real aFirstCoeff = aLeftDist / ( aLeftDist + aRightDist );
349 Standard_Real aSecCoeff = aRightDist / ( aLeftDist + aRightDist );
351 aResAltitude = aLeftAlt * aFirstCoeff + aRightAlt * aSecCoeff;