]> SALOME platform Git repositories - modules/hydro.git/blob - src/HYDROData/HYDROData_StreamAltitude.cxx
Salome HOME
Merge branch 'BR_LAND_COVER_MAP' into BR_quadtree
[modules/hydro.git] / src / HYDROData / HYDROData_StreamAltitude.cxx
1 // Copyright (C) 2014-2015  EDF-R&D
2 // This library is free software; you can redistribute it and/or
3 // modify it under the terms of the GNU Lesser General Public
4 // License as published by the Free Software Foundation; either
5 // version 2.1 of the License, or (at your option) any later version.
6 //
7 // This library is distributed in the hope that it will be useful,
8 // but WITHOUT ANY WARRANTY; without even the implied warranty of
9 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
10 // Lesser General Public License for more details.
11 //
12 // You should have received a copy of the GNU Lesser General Public
13 // License along with this library; if not, write to the Free Software
14 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
15 //
16 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
17 //
18
19 #include "HYDROData_StreamAltitude.h"
20
21 #include "HYDROData_Channel.h"
22 #include "HYDROData_Document.h"
23 #include "HYDROData_Profile.h"
24 #include "HYDROData_Stream.h"
25 #include "HYDROData_ShapesTool.h"
26 #include <HYDROData_Tool.h>
27 #include <BRep_Tool.hxx>
28
29 #include <BRepBuilderAPI_MakeEdge.hxx>
30 #include <BRepBuilderAPI_MakeFace.hxx>
31 #include <BRepBuilderAPI_MakeWire.hxx>
32
33 #include <Extrema_ExtElC.hxx>
34
35 #include <GeomAPI_ProjectPointOnCurve.hxx>
36
37 #include <gp_Lin.hxx>
38
39 #include <Precision.hxx>
40
41 #include <TopExp.hxx>
42 #include <TopExp_Explorer.hxx>
43
44 #include <TopoDS.hxx>
45 #include <TopoDS_Wire.hxx>
46 #include <TopoDS_Edge.hxx>
47
48 #include <TopTools_SequenceOfShape.hxx>
49 #include <TopTools_ListIteratorOfListOfShape.hxx>
50
51 #include <Geom_Line.hxx>
52
53 #include <QStringList>
54
55 #define _DEVDEBUG_
56 #include "HYDRO_trace.hxx"
57
58 #ifdef DEB_CLASS2D
59 #include <BRepTools.hxx>
60 #include <BRep_Builder.hxx>
61 #include <BRepBuilderAPI_MakeVertex.hxx>
62 #endif
63 IMPLEMENT_STANDARD_HANDLE(HYDROData_StreamAltitude, HYDROData_IAltitudeObject)
64 IMPLEMENT_STANDARD_RTTIEXT(HYDROData_StreamAltitude, HYDROData_IAltitudeObject)
65
66 HYDROData_StreamAltitude::HYDROData_StreamAltitude()
67 : HYDROData_IAltitudeObject()
68 {
69 }
70
71 HYDROData_StreamAltitude::~HYDROData_StreamAltitude()
72 {
73 }
74
75 bool IsPointBetweenEdges( const gp_Pnt& aFirstPnt1, const gp_Pnt& aLastPnt1,
76                           const gp_Pnt& aFirstPnt2, const gp_Pnt& aLastPnt2,
77                           const gp_Pnt& thePoint) {
78   BRepBuilderAPI_MakeEdge aLeftMakeEdge( aFirstPnt1, aLastPnt1 );
79   BRepBuilderAPI_MakeEdge aBotMakeEdge( aLastPnt1, aLastPnt2 );
80   BRepBuilderAPI_MakeEdge aRightMakeEdge( aLastPnt2, aFirstPnt2 );
81   BRepBuilderAPI_MakeEdge aTopMakeEdge( aFirstPnt2, aFirstPnt1 );
82
83   BRepBuilderAPI_MakeWire aMakeWire( aLeftMakeEdge.Edge(), aBotMakeEdge.Edge(), 
84                                      aRightMakeEdge.Edge(), aTopMakeEdge.Edge() );
85
86   BRepBuilderAPI_MakeFace aMakeFace( aMakeWire.Wire() );
87         
88   TopoDS_Face aFace = aMakeFace.Face();
89 #ifdef DEB_CLASS2D
90   TopoDS_Compound aCmp;
91   BRep_Builder aBB;
92   aBB.MakeCompound(aCmp);
93   aBB.Add(aCmp, aFace);
94   BRepBuilderAPI_MakeVertex aMk(thePoint);
95   aBB.Add(aCmp, aMk.Vertex());
96   BRepTools::Write(aCmp, "ProfileFace.brep");
97 #endif     
98
99   gp_XY anXY( thePoint.X(), thePoint.Y() );
100   TopAbs_State aPointState =  HYDROData_Tool::ComputePointState(anXY, aFace);
101
102 #ifdef DEB_CLASS2D
103           cout << "Point status is = " << aPointState <<endl;
104 #endif
105   return aPointState != TopAbs_OUT;
106 }
107
108 Standard_Real getAltitudeFromWire( const TopoDS_Wire& theWire,
109                                    const Standard_Real& theLeftDist,
110                                    const Standard_Real& theRightDist )
111 {
112   Standard_Real aResAlt = 0.0;
113
114   gp_XY aFirstPoint, aLastPoint;
115   TopoDS_Vertex aFirstVertex, aLastVertex;
116   TopExp::Vertices( theWire, aFirstVertex, aLastVertex );
117
118   gp_Pnt aPnt1( BRep_Tool::Pnt( aFirstVertex ) );
119   aPnt1.SetZ( 0 );
120   gp_Pnt aPnt2( BRep_Tool::Pnt( aLastVertex ) );
121   aPnt2.SetZ( 0 );
122   
123   Standard_Real aProfileDist = aPnt1.Distance( aPnt2 );
124
125   Standard_Real aCoeff = aProfileDist / ( theLeftDist + theRightDist );
126
127   gp_Pnt anIntPoint( aPnt1.XYZ() + ( aCoeff * theLeftDist ) * gp_Dir( gp_Vec( aPnt1, aPnt2 ) ).XYZ() );
128
129   return HYDROData_Tool::GetAltitudeForWire( theWire,
130                                              gp_XY(anIntPoint.X(), anIntPoint.Y()),
131                                              1E-2,
132                                              1E-2,
133                                              HYDROData_IAltitudeObject::GetInvalidAltitude() );
134 }
135
136 bool HYDROData_StreamAltitude::getBoundaryWiresForPoint(
137   const gp_XY& thePoint,
138   TopoDS_Wire& theLeftWire,
139   TopoDS_Wire& theRightWire ) const
140 {
141   gp_Pnt aTestPnt( thePoint.X(),  thePoint.Y(), 0 );
142
143   Handle(HYDROData_Object) anObject =
144     Handle(HYDROData_Object)::DownCast( GetFatherObject() );
145   if ( anObject.IsNull() ) {
146     return false;
147   }
148
149   if ( anObject->GetKind() == KIND_STREAM ) {
150     Handle(HYDROData_Stream) aStream = Handle(HYDROData_Stream)::DownCast( anObject );
151     if ( aStream.IsNull() )
152       return false;
153
154     HYDROData_SequenceOfObjects aStreamProfiles = aStream->GetProfiles();
155     if ( aStreamProfiles.Length() < 2 )
156       return false;
157
158     Handle(HYDROData_Profile) aPrevProfile;
159     gp_Pnt aPrevPnt1, aPrevPnt2;
160     for ( int i = 1, n = aStreamProfiles.Length(); i <= n; ++i )
161     {
162       Handle(HYDROData_Profile) aProfile =
163         Handle(HYDROData_Profile)::DownCast( aStreamProfiles.Value( i ) );
164       if ( aProfile.IsNull() )
165         continue;
166
167       gp_XY aFirstPoint, aLastPoint;
168       if ( !aProfile->GetLeftPoint( aFirstPoint, false ) ||
169            !aProfile->GetRightPoint( aLastPoint, false ) )
170         continue;
171       
172       gp_Pnt aPnt1( aFirstPoint.X(), aFirstPoint.Y(), 0 );
173       gp_Pnt aPnt2( aLastPoint.X(),  aLastPoint.Y(),  0 );
174
175       if ( !aPrevProfile.IsNull() )
176       {        
177         if ( IsPointBetweenEdges( aPrevPnt1, aPrevPnt2, aPnt1, aPnt2, aTestPnt ) )
178         {
179           theLeftWire = TopoDS::Wire( aPrevProfile->GetShape3D() );
180           theRightWire = TopoDS::Wire( aProfile->GetShape3D() );
181           break;
182         }
183       }
184
185       aPrevProfile = aProfile;
186       aPrevPnt1 = aPnt1;
187       aPrevPnt2 = aPnt2;
188     }
189   } else if ( anObject->GetKind() == KIND_CHANNEL ) {
190     Handle(HYDROData_Channel) aChannel = Handle(HYDROData_Channel)::DownCast( anObject );
191     if ( aChannel.IsNull() )
192       return false;
193
194     TopTools_ListOfShape aWiresList;
195     TopExp_Explorer anExp( aChannel->GetShape3D(), TopAbs_WIRE );
196     for ( ; anExp.More(); anExp.Next() ) {
197       if(!anExp.Current().IsNull()) {
198         const TopoDS_Wire& aWire = TopoDS::Wire( anExp.Current() );
199         aWiresList.Append( aWire );
200       }
201     }
202
203     if ( aWiresList.Extent() < 2 )
204       return false;
205
206     TopoDS_Wire aPrevWire;
207     gp_Pnt aPrevPnt1, aPrevPnt2;
208
209     TopTools_ListIteratorOfListOfShape anIt( aWiresList );
210     for ( ; anIt.More(); anIt.Next() ) {
211       TopoDS_Wire& aWire = TopoDS::Wire( anIt.Value() );
212       if ( aWire.IsNull() )
213         continue;
214
215       TopoDS_Vertex aFirstVertex, aLastVertex;
216       TopExp::Vertices( aWire, aFirstVertex, aLastVertex );
217     
218       gp_Pnt aPnt1( BRep_Tool::Pnt( aFirstVertex ) );
219       aPnt1.SetZ( 0 );
220       gp_Pnt aPnt2( BRep_Tool::Pnt( aLastVertex ) );
221       aPnt2.SetZ( 0 );
222
223       if ( !aPrevWire.IsNull() ) {
224         if ( IsPointBetweenEdges( aPrevPnt1, aPrevPnt2, aPnt1, aPnt2, aTestPnt ) ) {
225           theLeftWire = aPrevWire;
226           theRightWire = aWire;
227           break;
228         }
229       }
230
231       aPrevWire = aWire;
232       aPrevPnt1 = aPnt1;
233       aPrevPnt2 = aPnt2;
234     }
235   }
236
237   return !theLeftWire.IsNull() && !theRightWire.IsNull();
238 }
239
240 double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) const
241 {
242   DEBTRACE("GetAltitudeForPoint p(" << thePoint.X() << ", " << thePoint.Y() << ")");
243   double aResAltitude = GetInvalidAltitude();
244
245   Handle(HYDROData_Object) anObject =
246     Handle(HYDROData_Object)::DownCast( GetFatherObject() );
247   if ( anObject.IsNull() )
248   {
249     DEBTRACE("aStream.IsNull()");
250     return aResAltitude;
251   }
252   TopoDS_Shape aTopShape = anObject->GetTopShape();
253   if ( aTopShape.IsNull() )
254   {
255     DEBTRACE("aTopShape.IsNull()");
256     return aResAltitude;
257   }
258
259   TopExp_Explorer aFaceExp( aTopShape, TopAbs_FACE );
260   if ( !aFaceExp.More() )
261   {
262     DEBTRACE("!aFaceExp.More()");
263     return aResAltitude;
264   }
265
266   // Get only face because of 2d profile wires is in compound
267   TopoDS_Face aFace = TopoDS::Face( aFaceExp.Current() );
268
269   // Check if point is inside of stream/channel presentation
270   TopAbs_State aPointState = HYDROData_Tool::ComputePointState(thePoint, aFace);
271
272 #ifdef DEB_CLASS2D
273       cout << "Point status is = " << aPointState <<endl;
274           TopoDS_Compound aCmp;
275       BRep_Builder aBB;
276       aBB.MakeCompound(aCmp);
277           aBB.Add(aCmp, aFace);
278           gp_Pnt aPnt (thePoint.X(), thePoint.Y(), 0.);
279           BRepBuilderAPI_MakeVertex aMk(aPnt);
280           aBB.Add(aCmp, aMk.Vertex());
281           BRepTools::Write(aCmp, "FCL2d.brep");
282 #endif
283   if ( aPointState == TopAbs_OUT )
284   {
285         DEBTRACE("aPointState == TopAbs_OUT");
286     return aResAltitude;
287   }
288
289   TopoDS_Edge aLeftEdge, aRightEdge;
290
291   if ( anObject->GetKind() == KIND_STREAM ) {
292     Handle(HYDROData_Stream) aStream =
293       Handle(HYDROData_Stream)::DownCast( GetFatherObject() );
294     if ( !aStream.IsNull() ) {
295       aLeftEdge = TopoDS::Edge( aStream->GetLeftShape() );
296       aRightEdge = TopoDS::Edge( aStream->GetRightShape() );
297     }
298   } else if ( anObject->GetKind() == KIND_CHANNEL ) {
299     Handle(HYDROData_Channel) aChannel =
300       Handle(HYDROData_Channel)::DownCast( GetFatherObject() );
301     if ( !aChannel.IsNull() ) {
302       aLeftEdge = TopoDS::Edge( aChannel->GetLeftShape() );
303       aRightEdge = TopoDS::Edge( aChannel->GetRightShape() );
304     }
305   }
306   
307   // Find the two profiles between which the point is lies
308   TopoDS_Wire aLeftWire, aRightWire;
309   if ( !getBoundaryWiresForPoint( thePoint, aLeftWire, aRightWire ) )
310   {
311       DEBTRACE("!getBoundaryProfilesForPoint( thePoint, aLeftProfile, aRightProfile )");
312       return aResAltitude;
313   }
314
315   // Find the projections of point to borders of stream
316   gp_XYZ aPointToTest( thePoint.X(), thePoint.Y(), 0.0 );
317   
318   Standard_Real aFirst = 0.0, aLast = 0.0;
319   Handle(Geom_Curve) anEdgeLeftCurve = BRep_Tool::Curve( aLeftEdge, aFirst, aLast );
320   Handle(Geom_Curve) anEdgeRightCurve = BRep_Tool::Curve( aRightEdge, aFirst, aLast );
321
322   GeomAPI_ProjectPointOnCurve aLeftProject( aPointToTest, anEdgeLeftCurve );
323   GeomAPI_ProjectPointOnCurve aRightProject( aPointToTest, anEdgeRightCurve );
324
325   int aNbPoints1 = aRightProject.NbPoints();
326   int aNbPoints2 = aLeftProject.NbPoints();
327   if ( aNbPoints1 < 1 || aNbPoints2 < 1)
328      return aResAltitude;
329   Standard_Real aLeftDist = aLeftProject.LowerDistance();
330   Standard_Real aRightDist = aRightProject.LowerDistance();
331
332   // Find the altitude in profiles
333   Standard_Real aLeftAlt, aRightAlt;
334   gp_Pnt aLeftProfileP1, aLeftProfileP2, aRightProfileP1, aRightProfileP2;
335   aLeftAlt = getAltitudeFromWire( aLeftWire, aLeftDist, aRightDist );
336   aRightAlt = getAltitudeFromWire( aRightWire, aLeftDist, aRightDist );
337   
338   TopoDS_Vertex aFirstVertex, aLastVertex;
339   TopExp::Vertices( aLeftWire, aFirstVertex, aLastVertex );
340   aLeftProfileP1 = BRep_Tool::Pnt( aFirstVertex );
341   aLeftProfileP2 = BRep_Tool::Pnt( aLastVertex );
342
343   TopExp::Vertices( aRightWire, aFirstVertex, aLastVertex );
344   aRightProfileP1 = BRep_Tool::Pnt( aFirstVertex );
345   aRightProfileP2 = BRep_Tool::Pnt( aLastVertex );
346
347   // Interpolate altitudes
348   // Left profile line ( the segment between the firts and the last profile point )
349   aLeftProfileP1.SetZ( 0 );
350   aLeftProfileP2.SetZ( 0 );
351   gp_Vec aLeftProfileVec( aLeftProfileP1, aLeftProfileP2 );
352   Handle(Geom_Line) aLeftProfileLine = new Geom_Line( gp_Ax1( aLeftProfileP1, aLeftProfileVec ) );
353   // Right profile line
354   aRightProfileP1.SetZ( 0 );
355   aRightProfileP2.SetZ( 0 );
356   gp_Vec aRightProfileVec( aRightProfileP1, aRightProfileP2 );
357   Handle(Geom_Line) aRightProfileLine = new Geom_Line( gp_Ax1( aRightProfileP1, aRightProfileVec ) );
358   // The point projections on the left and right profiles
359   GeomAPI_ProjectPointOnCurve aLeftProfileProject( aPointToTest, aLeftProfileLine );
360   GeomAPI_ProjectPointOnCurve aRightProfileProject( aPointToTest, aRightProfileLine );
361   // The point distance to the left and right profiles
362   Standard_Real aLeftProfileDist = aLeftProfileProject.LowerDistance();
363   Standard_Real aRightProfileDist = aRightProfileProject.LowerDistance();
364   // The coefficients
365   Standard_Real aFirstCoeff = aLeftProfileDist / ( aLeftProfileDist + aRightProfileDist );
366   Standard_Real aSecCoeff = aRightProfileDist / ( aLeftProfileDist + aRightProfileDist );
367
368   aResAltitude = aLeftAlt * ( 1 - aFirstCoeff ) + aRightAlt * ( 1 - aSecCoeff );
369   DEBTRACE("aResAltitude=" << aResAltitude);
370   return aResAltitude;
371 }
372
373
374
375