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