]> SALOME platform Git repositories - modules/hydro.git/blob - src/HYDROData/HYDROData_StreamAltitude.cxx
Salome HOME
17.01.2014. Redesign of workaround for calcualtion point state.
[modules/hydro.git] / src / HYDROData / HYDROData_StreamAltitude.cxx
1
2 #include "HYDROData_StreamAltitude.h"
3
4 #include "HYDROData_Document.h"
5 #include "HYDROData_Profile.h"
6 #include "HYDROData_Stream.h"
7 #include "HYDROData_ShapesTool.h"
8 #include <HYDROData_Tool.h>
9 #include <BRep_Tool.hxx>
10
11 #include <BRepBuilderAPI_MakeEdge.hxx>
12 #include <BRepBuilderAPI_MakeFace.hxx>
13 #include <BRepBuilderAPI_MakeWire.hxx>
14
15 #include <Extrema_ExtElC.hxx>
16
17 #include <GeomAPI_ProjectPointOnCurve.hxx>
18
19 #include <gp_Lin.hxx>
20
21 #include <Precision.hxx>
22
23 #include <TopExp_Explorer.hxx>
24
25 #include <TopoDS.hxx>
26 #include <TopoDS_Wire.hxx>
27 #include <TopoDS_Edge.hxx>
28
29 #include <TopTools_SequenceOfShape.hxx>
30
31 #include <QStringList>
32
33 #ifdef DEB_CLASS2D
34 #include <BRepTools.hxx>
35 #include <BRep_Builder.hxx>
36 #include <BRepBuilderAPI_MakeVertex.hxx>
37 #endif
38 IMPLEMENT_STANDARD_HANDLE(HYDROData_StreamAltitude, HYDROData_IAltitudeObject)
39 IMPLEMENT_STANDARD_RTTIEXT(HYDROData_StreamAltitude, HYDROData_IAltitudeObject)
40
41 HYDROData_StreamAltitude::HYDROData_StreamAltitude()
42 : HYDROData_IAltitudeObject()
43 {
44 }
45
46 HYDROData_StreamAltitude::~HYDROData_StreamAltitude()
47 {
48 }
49
50 Standard_Real getAltitudeFromProfile( const Handle(HYDROData_Profile)& theProfile,
51                                       const Standard_Real&             theLeftDist,
52                                       const Standard_Real&             theRightDist )
53 {
54   Standard_Real aResAlt = 0.0;
55
56   gp_XY aFirstPoint, aLastPoint;
57   if ( !theProfile->GetLeftPoint( aFirstPoint ) ||
58        !theProfile->GetRightPoint( aLastPoint ) )
59     return aResAlt;
60
61   gp_Pnt aPnt1( aFirstPoint.X(), aFirstPoint.Y(), 0 );
62   gp_Pnt aPnt2( aLastPoint.X(),  aLastPoint.Y(),  0 );
63
64   Standard_Real aProfileDist = aPnt1.Distance( aPnt2 );
65
66   Standard_Real aCoeff = aProfileDist / ( theLeftDist + theRightDist );
67
68   gp_Pnt anIntPoint( aPnt1.XYZ() + ( aCoeff * theLeftDist ) * gp_Dir( gp_Vec( aPnt1, aPnt2 ) ).XYZ() );
69
70   gp_Lin aPointLine( anIntPoint, gp::DZ() );
71
72   gp_Pnt aPrevPoint;
73   gp_Lin aPrevNormal;
74   HYDROData_Profile::ProfilePoints aProfilePoints = theProfile->GetProfilePoints();
75   for ( int i = 1, n = aProfilePoints.Length(); i <= n; ++i )
76   {
77     gp_Pnt aProfPoint( aProfilePoints.Value( i ) );
78
79     Standard_Real aDist = aPointLine.Distance( aProfPoint );
80     if ( aDist <= gp::Resolution() )
81     {
82       // We found the intersected point
83       aResAlt = aProfPoint.Z();
84       break;
85     }
86    
87     gp_Lin aNormal = aPointLine.Normal( aProfPoint );
88     if ( i == 1 )
89     {
90       aPrevNormal = aNormal;
91       aPrevPoint = aProfPoint;
92       continue;
93     }
94
95     if ( aPrevNormal.Direction().Dot( aNormal.Direction() ) < 0 )
96     {
97       // We found the intersected edge
98       gp_Lin anEdgeLine( aPrevPoint, gp_Dir( gp_Vec( aPrevPoint, aProfPoint ) ) );
99      
100       Extrema_ExtElC anExtrema( aPointLine, anEdgeLine, Precision::Angular() );
101       if ( !anExtrema.IsParallel() )
102       {
103         Extrema_POnCurv aFirstPnt, aSecPnt;
104         anExtrema.Points( 1, aFirstPnt, aSecPnt );
105
106         const gp_Pnt& anIntPnt = aSecPnt.Value();
107         aResAlt = anIntPnt.Z();
108
109         break;
110       }
111     }
112
113     aPrevNormal = aNormal;
114     aPrevPoint = aProfPoint;
115   }
116
117   return aResAlt;
118 }
119
120 bool HYDROData_StreamAltitude::getBoundaryProfilesForPoint(
121   const gp_XY&               thePoint,
122   Handle(HYDROData_Profile)& theLeftProfile,
123   Handle(HYDROData_Profile)& theRightProfile ) const
124 {
125   Handle(HYDROData_Stream) aStream =
126     Handle(HYDROData_Stream)::DownCast( GetFatherObject() );
127   if ( aStream.IsNull() )
128     return false;
129
130   HYDROData_SequenceOfObjects aStreamProfiles = aStream->GetProfiles();
131   if ( aStreamProfiles.Length() < 2 )
132     return false;
133
134   Handle(HYDROData_Profile) aPrevProfile;
135   gp_Pnt aPrevPnt1, aPrevPnt2;
136   for ( int i = 1, n = aStreamProfiles.Length(); i <= n; ++i )
137   {
138     Handle(HYDROData_Profile) aProfile =
139       Handle(HYDROData_Profile)::DownCast( aStreamProfiles.Value( i ) );
140     if ( aProfile.IsNull() )
141       continue;
142
143     gp_XY aFirstPoint, aLastPoint;
144     if ( !aProfile->GetLeftPoint( aFirstPoint ) || !aProfile->GetRightPoint( aLastPoint ) )
145       continue;
146
147     gp_Pnt aPnt1( aFirstPoint.X(), aFirstPoint.Y(), 0 );
148     gp_Pnt aPnt2( aLastPoint.X(),  aLastPoint.Y(),  0 );
149
150     if ( !aPrevProfile.IsNull() )
151     {
152       BRepBuilderAPI_MakeEdge aLeftMakeEdge( aPrevPnt1, aPrevPnt2 );
153       BRepBuilderAPI_MakeEdge aBotMakeEdge( aPrevPnt2, aPnt2 );
154       BRepBuilderAPI_MakeEdge aRightMakeEdge( aPnt2, aPnt1 );
155       BRepBuilderAPI_MakeEdge aTopMakeEdge( aPnt1, aPrevPnt1 );
156
157       BRepBuilderAPI_MakeWire aMakeWire( aLeftMakeEdge.Edge(), aBotMakeEdge.Edge(), 
158                                          aRightMakeEdge.Edge(), aTopMakeEdge.Edge() );
159
160       BRepBuilderAPI_MakeFace aMakeFace( aMakeWire.Wire() );
161         
162       TopoDS_Face aProfilesFace = aMakeFace.Face();
163 #ifdef DEB_CLASS2D
164           TopoDS_Compound aCmp;
165       BRep_Builder aBB;
166       aBB.MakeCompound(aCmp);
167           aBB.Add(aCmp, aProfilesFace);
168           gp_Pnt aPnt (thePoint.X(), thePoint.Y(), 0.);
169           BRepBuilderAPI_MakeVertex aMk(aPnt);
170           aBB.Add(aCmp, aMk.Vertex());
171           BRepTools::Write(aCmp, "ProfileFace.brep");
172 #endif     
173
174           TopAbs_State aPointState =  HYDROData_Tool::ComputePointState(thePoint, aProfilesFace);
175
176 #ifdef DEB_CLASS2D
177           cout << "Point status is = " << aPointState <<endl;
178 #endif
179       if ( aPointState != TopAbs_OUT )
180       {
181         theLeftProfile = aPrevProfile;
182         theRightProfile = aProfile;
183         break;
184       }
185     }
186
187     aPrevProfile = aProfile;
188     aPrevPnt1 = aPnt1;
189     aPrevPnt2 = aPnt2;
190   }
191
192   return !theLeftProfile.IsNull() && !theRightProfile.IsNull();
193 }
194
195 double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) const
196 {
197   double aResAltitude = GetInvalidAltitude();
198
199   Handle(HYDROData_Stream) aStream =
200     Handle(HYDROData_Stream)::DownCast( GetFatherObject() );
201   if ( aStream.IsNull() )
202     return aResAltitude;
203
204   TopoDS_Shape aStreamShape = aStream->GetTopShape();
205   if ( aStreamShape.IsNull() )
206     return aResAltitude;
207
208   TopExp_Explorer aStreamFaceExp( aStreamShape, TopAbs_FACE );
209   if ( !aStreamFaceExp.More() )
210     return aResAltitude;
211
212   // Get only face because of 2d profile wires is in compound
213   TopoDS_Face aStreamFace = TopoDS::Face( aStreamFaceExp.Current() );
214
215   // Check if point is inside of stream presentation
216   TopAbs_State aPointState = HYDROData_Tool::ComputePointState(thePoint, aStreamFace);
217
218 #ifdef DEB_CLASS2D
219       cout << "Point status is = " << aPointState <<endl;
220           TopoDS_Compound aCmp;
221       BRep_Builder aBB;
222       aBB.MakeCompound(aCmp);
223           aBB.Add(aCmp, aStreamFace);
224           gp_Pnt aPnt (thePoint.X(), thePoint.Y(), 0.);
225           BRepBuilderAPI_MakeVertex aMk(aPnt);
226           aBB.Add(aCmp, aMk.Vertex());
227           BRepTools::Write(aCmp, "FCL2d.brep");
228 #endif
229   if ( aPointState == TopAbs_OUT )
230     return aResAltitude;
231
232   // Find the two profiles between which the point is lies
233   Handle(HYDROData_Profile) aLeftProfile, aRightProfile;
234   if ( !getBoundaryProfilesForPoint( thePoint, aLeftProfile, aRightProfile ) )
235     return aResAltitude;
236
237   // Find the projections of point to borders of stream
238   gp_XYZ aPointToTest( thePoint.X(), thePoint.Y(), 0.0 );
239
240   TopoDS_Edge aStreamLeftEdge = TopoDS::Edge( aStream->GetLeftShape() );
241   TopoDS_Edge aStreamRightEdge = TopoDS::Edge( aStream->GetRightShape() );
242
243   Standard_Real aFirst = 0.0, aLast = 0.0;
244   Handle(Geom_Curve) anEdgeLeftCurve = BRep_Tool::Curve( aStreamLeftEdge, aFirst, aLast );
245   Handle(Geom_Curve) anEdgeRightCurve = BRep_Tool::Curve( aStreamRightEdge, aFirst, aLast );
246
247   GeomAPI_ProjectPointOnCurve aLeftProject( aPointToTest, anEdgeLeftCurve );
248   GeomAPI_ProjectPointOnCurve aRightProject( aPointToTest, anEdgeRightCurve );
249
250   Standard_Real aLeftDist = aLeftProject.LowerDistance();
251   Standard_Real aRightDist = aRightProject.LowerDistance();
252
253   // Find the altitude in profiles
254   Standard_Real aLeftAlt = getAltitudeFromProfile( aLeftProfile, aLeftDist, aRightDist );
255   Standard_Real aRightAlt = getAltitudeFromProfile( aRightProfile, aLeftDist, aRightDist );
256
257   // Interpolate altitudes
258   Standard_Real aFirstCoeff = aLeftDist / ( aLeftDist + aRightDist );
259   Standard_Real aSecCoeff = aRightDist / ( aLeftDist + aRightDist );
260
261   aResAltitude = aLeftAlt * aFirstCoeff + aRightAlt * aSecCoeff;
262
263   return aResAltitude;
264 }
265
266
267
268