]> SALOME platform Git repositories - modules/hydro.git/blob - src/HYDROData/HYDROData_StreamAltitude.cxx
Salome HOME
merge BR_PORTING_OCCT_7
[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 #include <cmath>
55
56 #define _DEVDEBUG_
57 #include "HYDRO_trace.hxx"
58
59 #ifdef DEB_CLASS2D
60 #include <BRepTools.hxx>
61 #include <BRep_Builder.hxx>
62 #include <BRepBuilderAPI_MakeVertex.hxx>
63 #endif
64
65 IMPLEMENT_STANDARD_RTTIEXT(HYDROData_StreamAltitude, HYDROData_IAltitudeObject)
66
67 HYDROData_StreamAltitude::HYDROData_StreamAltitude()
68 : HYDROData_IAltitudeObject()
69 {
70 }
71
72 HYDROData_StreamAltitude::~HYDROData_StreamAltitude()
73 {
74 }
75
76 bool IsPointBetweenEdges( const gp_Pnt& aFirstPnt1, const gp_Pnt& aLastPnt1,
77                           const gp_Pnt& aFirstPnt2, const gp_Pnt& aLastPnt2,
78                           const gp_Pnt& thePoint) {
79   BRepBuilderAPI_MakeEdge aLeftMakeEdge( aFirstPnt1, aLastPnt1 );
80   BRepBuilderAPI_MakeEdge aBotMakeEdge( aLastPnt1, aLastPnt2 );
81   BRepBuilderAPI_MakeEdge aRightMakeEdge( aLastPnt2, aFirstPnt2 );
82   BRepBuilderAPI_MakeEdge aTopMakeEdge( aFirstPnt2, aFirstPnt1 );
83
84   BRepBuilderAPI_MakeWire aMakeWire( aLeftMakeEdge.Edge(), aBotMakeEdge.Edge(), 
85                                      aRightMakeEdge.Edge(), aTopMakeEdge.Edge() );
86
87   BRepBuilderAPI_MakeFace aMakeFace( aMakeWire.Wire() );
88         
89   TopoDS_Face aFace = aMakeFace.Face();
90 #ifdef DEB_CLASS2D
91   TopoDS_Compound aCmp;
92   BRep_Builder aBB;
93   aBB.MakeCompound(aCmp);
94   aBB.Add(aCmp, aFace);
95   BRepBuilderAPI_MakeVertex aMk(thePoint);
96   aBB.Add(aCmp, aMk.Vertex());
97   BRepTools::Write(aCmp, "ProfileFace.brep");
98 #endif     
99
100   gp_XY anXY( thePoint.X(), thePoint.Y() );
101   TopAbs_State aPointState =  HYDROData_Tool::ComputePointState(anXY, aFace);
102
103 #ifdef DEB_CLASS2D
104           cout << "Point status is = " << aPointState <<endl;
105 #endif
106   return aPointState != TopAbs_OUT;
107 }
108
109 /** compute edge1^(firstPnt1,point) and edge2^(firstPnt2,point):
110  *  if the z values are opposite, the point is between the edges.
111  *  We must also check if the point is near the edges
112  *  (inside the circle defined by the edges: approximation)
113  *  to discriminate false positives in sinuous cases
114  */
115 bool IsPointBetweenEdges2( const gp_Pnt& aFirstPnt1, const gp_Pnt& aLastPnt1,
116                            const gp_Pnt& aFirstPnt2, const gp_Pnt& aLastPnt2,
117                            const gp_Pnt& thePoint) {
118   double x1 = aLastPnt1.X() - aFirstPnt1.X(); // v1
119   double y1 = aLastPnt1.Y() - aFirstPnt1.Y();
120   double x2 = aLastPnt2.X() - aFirstPnt2.X(); // v2
121   double y2 = aLastPnt2.Y() - aFirstPnt2.Y();
122   double xa = thePoint.X() - aFirstPnt1.X();  // va
123   double ya = thePoint.Y() - aFirstPnt1.Y();
124   double xb = thePoint.X() - aFirstPnt2.X();  // vb
125   double yb = thePoint.Y() - aFirstPnt2.Y();
126   double z1 = x1*ya -xa*y1; // v1^va: z component
127   double z2 = x2*yb -xb*y2; // v2^vb: z component
128   bool isBetween = true;
129   if (((z1<0) && (z2<0)) || ((z1>=0) && (z2>=0)))
130     {
131       isBetween = false;
132     }
133   if (isBetween)
134     {
135       double xg = (aFirstPnt1.X() + aLastPnt1.X() + aFirstPnt2.X() + aLastPnt2.X())/4;
136       double yg = (aFirstPnt1.Y() + aLastPnt1.Y() + aFirstPnt2.Y() + aLastPnt2.Y())/4;
137       double df1 = (aFirstPnt1.X()-xg)*(aFirstPnt1.X()-xg) + (aFirstPnt1.Y()-yg)*(aFirstPnt1.Y()-yg);
138       double dl1 = (aLastPnt1.X()-xg)*(aLastPnt1.X()-xg) + (aLastPnt1.Y()-yg)*(aLastPnt1.Y()-yg);
139       double df2 = (aFirstPnt2.X()-xg)*(aFirstPnt2.X()-xg) + (aFirstPnt2.Y()-yg)*(aFirstPnt2.Y()-yg);
140       double dl2 = (aLastPnt2.X()-xg)*(aLastPnt2.X()-xg) + (aLastPnt2.Y()-yg)*(aLastPnt2.Y()-yg);
141       double r2 = std::max(df1,dl1);
142       r2 = std::max(r2,df2);
143       r2 = std::max(r2,dl2);
144       double d2 = (thePoint.X()-xg)*(thePoint.X()-xg) + (thePoint.Y()-yg)*(thePoint.Y()-yg);
145       if (d2 > r2)
146         isBetween = false;
147     }
148   return isBetween;
149 }
150
151 Standard_Real getAltitudeFromWire( const TopoDS_Wire& theWire,
152                                    const Standard_Real& theLeftDist,
153                                    const Standard_Real& theRightDist )
154 {
155   Standard_Real aResAlt = 0.0;
156
157   gp_XY aFirstPoint, aLastPoint;
158   TopoDS_Vertex aFirstVertex, aLastVertex;
159   TopExp::Vertices( theWire, aFirstVertex, aLastVertex );
160
161   gp_Pnt aPnt1( BRep_Tool::Pnt( aFirstVertex ) );
162   aPnt1.SetZ( 0 );
163   gp_Pnt aPnt2( BRep_Tool::Pnt( aLastVertex ) );
164   aPnt2.SetZ( 0 );
165   
166   Standard_Real aProfileDist = aPnt1.Distance( aPnt2 );
167
168   Standard_Real aCoeff = aProfileDist / ( theLeftDist + theRightDist );
169
170   gp_Pnt anIntPoint( aPnt1.XYZ() + ( aCoeff * theLeftDist ) * gp_Dir( gp_Vec( aPnt1, aPnt2 ) ).XYZ() );
171
172   return HYDROData_Tool::GetAltitudeForWire( theWire,
173                                              gp_XY(anIntPoint.X(), anIntPoint.Y()),
174                                              1E-2,
175                                              1E-2,
176                                              HYDROData_IAltitudeObject::GetInvalidAltitude() );
177 }
178
179 bool HYDROData_StreamAltitude::getBoundaryWiresForPoint(
180   const gp_XY& thePoint,
181   TopoDS_Wire& theLeftWire,
182   TopoDS_Wire& theRightWire ) const
183 {
184   gp_Pnt aTestPnt( thePoint.X(),  thePoint.Y(), 0 );
185
186   Handle(HYDROData_Object) anObject =
187     Handle(HYDROData_Object)::DownCast( GetFatherObject() );
188   if ( anObject.IsNull() ) {
189     return false;
190   }
191
192   if ( anObject->GetKind() == KIND_STREAM ) {
193     Handle(HYDROData_Stream) aStream = Handle(HYDROData_Stream)::DownCast( anObject );
194     if ( aStream.IsNull() )
195       return false;
196
197     HYDROData_SequenceOfObjects aStreamProfiles = aStream->GetProfiles();
198     if ( aStreamProfiles.Length() < 2 )
199       return false;
200
201     Handle(HYDROData_Profile) aPrevProfile;
202     gp_Pnt aPrevPnt1, aPrevPnt2;
203     for ( int i = 1, n = aStreamProfiles.Length(); i <= n; ++i )
204     {
205       Handle(HYDROData_Profile) aProfile =
206         Handle(HYDROData_Profile)::DownCast( aStreamProfiles.Value( i ) );
207       if ( aProfile.IsNull() )
208         continue;
209
210       gp_XY aFirstPoint, aLastPoint;
211       if ( !aProfile->GetLeftPoint( aFirstPoint, false ) ||
212            !aProfile->GetRightPoint( aLastPoint, false ) )
213         continue;
214       
215       gp_Pnt aPnt1( aFirstPoint.X(), aFirstPoint.Y(), 0 );
216       gp_Pnt aPnt2( aLastPoint.X(),  aLastPoint.Y(),  0 );
217
218       if ( !aPrevProfile.IsNull() )
219       {        
220         if ( IsPointBetweenEdges2( aPrevPnt1, aPrevPnt2, aPnt1, aPnt2, aTestPnt ) )
221         {
222           DEBTRACE("point is between profiles "<< i-1 << " and " <<i);
223           theLeftWire = TopoDS::Wire( aPrevProfile->GetShape3D() );
224           theRightWire = TopoDS::Wire( aProfile->GetShape3D() );
225           break;
226         }
227       }
228
229       aPrevProfile = aProfile;
230       aPrevPnt1 = aPnt1;
231       aPrevPnt2 = aPnt2;
232     }
233   } else if ( anObject->GetKind() == KIND_CHANNEL ) {
234     Handle(HYDROData_Channel) aChannel = Handle(HYDROData_Channel)::DownCast( anObject );
235     if ( aChannel.IsNull() )
236       return false;
237
238     TopTools_ListOfShape aWiresList;
239     TopExp_Explorer anExp( aChannel->GetShape3D(), TopAbs_WIRE );
240     for ( ; anExp.More(); anExp.Next() ) {
241       if(!anExp.Current().IsNull()) {
242         const TopoDS_Wire& aWire = TopoDS::Wire( anExp.Current() );
243         aWiresList.Append( aWire );
244       }
245     }
246
247     if ( aWiresList.Extent() < 2 )
248       return false;
249
250     TopoDS_Wire aPrevWire;
251     gp_Pnt aPrevPnt1, aPrevPnt2;
252
253     TopTools_ListIteratorOfListOfShape anIt( aWiresList );
254     for ( ; anIt.More(); anIt.Next() ) {
255       TopoDS_Wire& aWire = TopoDS::Wire( anIt.Value() );
256       if ( aWire.IsNull() )
257         continue;
258
259       TopoDS_Vertex aFirstVertex, aLastVertex;
260       TopExp::Vertices( aWire, aFirstVertex, aLastVertex );
261     
262       gp_Pnt aPnt1( BRep_Tool::Pnt( aFirstVertex ) );
263       aPnt1.SetZ( 0 );
264       gp_Pnt aPnt2( BRep_Tool::Pnt( aLastVertex ) );
265       aPnt2.SetZ( 0 );
266
267       if ( !aPrevWire.IsNull() ) {
268         if ( IsPointBetweenEdges2( aPrevPnt1, aPrevPnt2, aPnt1, aPnt2, aTestPnt ) ) {
269           theLeftWire = aPrevWire;
270           theRightWire = aWire;
271           break;
272         }
273       }
274
275       aPrevWire = aWire;
276       aPrevPnt1 = aPnt1;
277       aPrevPnt2 = aPnt2;
278     }
279   }
280
281   return !theLeftWire.IsNull() && !theRightWire.IsNull();
282 }
283
284 double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint,
285                                                       int theMethod) const
286 {
287   DEBTRACE("GetAltitudeForPoint p(" << thePoint.X() << ", " << thePoint.Y() << ")");
288   double aResAltitude = GetInvalidAltitude();
289
290   Handle(HYDROData_Object) anObject =
291     Handle(HYDROData_Object)::DownCast( GetFatherObject() );
292   if ( anObject.IsNull() )
293   {
294     DEBTRACE(" --- aStream.IsNull()");
295     return aResAltitude;
296   }
297   TopoDS_Shape aTopShape = anObject->GetTopShape();
298   if ( aTopShape.IsNull() )
299   {
300     DEBTRACE(" --- aTopShape.IsNull()");
301     return aResAltitude;
302   }
303
304   TopExp_Explorer aFaceExp( aTopShape, TopAbs_FACE );
305   if ( !aFaceExp.More() )
306   {
307     DEBTRACE(" --- !aFaceExp.More()");
308     return aResAltitude;
309   }
310
311   // Get only face because of 2d profile wires is in compound
312   TopoDS_Face aFace = TopoDS::Face( aFaceExp.Current() );
313
314   // Check if point is inside of stream/channel presentation
315   TopAbs_State aPointState = HYDROData_Tool::ComputePointState(thePoint, aFace);
316
317 #ifdef DEB_CLASS2D
318       cout << "Point status is = " << aPointState <<endl;
319           TopoDS_Compound aCmp;
320       BRep_Builder aBB;
321       aBB.MakeCompound(aCmp);
322           aBB.Add(aCmp, aFace);
323           gp_Pnt aPnt (thePoint.X(), thePoint.Y(), 0.);
324           BRepBuilderAPI_MakeVertex aMk(aPnt);
325           aBB.Add(aCmp, aMk.Vertex());
326           BRepTools::Write(aCmp, "FCL2d.brep");
327 #endif
328   if ( aPointState == TopAbs_OUT )
329   {
330         DEBTRACE(" --- aPointState == TopAbs_OUT");
331     return aResAltitude;
332   }
333
334   TopoDS_Edge aLeftEdge, aRightEdge;
335
336   if ( anObject->GetKind() == KIND_STREAM ) {
337     Handle(HYDROData_Stream) aStream =
338       Handle(HYDROData_Stream)::DownCast( GetFatherObject() );
339     if ( !aStream.IsNull() ) {
340       aLeftEdge = TopoDS::Edge( aStream->GetLeftShape() );
341       aRightEdge = TopoDS::Edge( aStream->GetRightShape() );
342     }
343   } else if ( anObject->GetKind() == KIND_CHANNEL ) {
344     Handle(HYDROData_Channel) aChannel =
345       Handle(HYDROData_Channel)::DownCast( GetFatherObject() );
346     if ( !aChannel.IsNull() ) {
347       aLeftEdge = TopoDS::Edge( aChannel->GetLeftShape() );
348       aRightEdge = TopoDS::Edge( aChannel->GetRightShape() );
349     }
350   }
351   
352   // Find the two profiles between which the point is lies
353   TopoDS_Wire aLeftWire, aRightWire;
354   if ( !getBoundaryWiresForPoint( thePoint, aLeftWire, aRightWire ) )
355   {
356       DEBTRACE(" --- !getBoundaryProfilesForPoint( thePoint, aLeftProfile, aRightProfile )");
357       return aResAltitude;
358   }
359
360   // Find the projections of point to borders of stream
361   gp_XYZ aPointToTest( thePoint.X(), thePoint.Y(), 0.0 );
362   
363   Standard_Real aFirst = 0.0, aLast = 0.0;
364   Handle(Geom_Curve) anEdgeLeftCurve = BRep_Tool::Curve( aLeftEdge, aFirst, aLast );
365   Handle(Geom_Curve) anEdgeRightCurve = BRep_Tool::Curve( aRightEdge, aFirst, aLast );
366
367   GeomAPI_ProjectPointOnCurve aLeftProject( aPointToTest, anEdgeLeftCurve );
368   GeomAPI_ProjectPointOnCurve aRightProject( aPointToTest, anEdgeRightCurve );
369
370   int aNbPoints1 = aRightProject.NbPoints();
371   int aNbPoints2 = aLeftProject.NbPoints();
372   if ( aNbPoints1 < 1 || aNbPoints2 < 1)
373     {
374       DEBTRACE(" --- projections " << aNbPoints1 << " " << aNbPoints2);
375       return aResAltitude;
376     }
377   Standard_Real aLeftDist = aLeftProject.LowerDistance();
378   Standard_Real aRightDist = aRightProject.LowerDistance();
379
380   // Find the altitude in profiles
381   Standard_Real aLeftAlt, aRightAlt;
382   gp_Pnt aLeftProfileP1, aLeftProfileP2, aRightProfileP1, aRightProfileP2;
383   aLeftAlt = getAltitudeFromWire( aLeftWire, aLeftDist, aRightDist );
384   aRightAlt = getAltitudeFromWire( aRightWire, aLeftDist, aRightDist );
385   
386   TopoDS_Vertex aFirstVertex, aLastVertex;
387   TopExp::Vertices( aLeftWire, aFirstVertex, aLastVertex );
388   aLeftProfileP1 = BRep_Tool::Pnt( aFirstVertex );
389   aLeftProfileP2 = BRep_Tool::Pnt( aLastVertex );
390
391   TopExp::Vertices( aRightWire, aFirstVertex, aLastVertex );
392   aRightProfileP1 = BRep_Tool::Pnt( aFirstVertex );
393   aRightProfileP2 = BRep_Tool::Pnt( aLastVertex );
394
395   // Interpolate altitudes
396   // Left profile line ( the segment between the first and the last profile point )
397   aLeftProfileP1.SetZ( 0 );
398   aLeftProfileP2.SetZ( 0 );
399   gp_Vec aLeftProfileVec( aLeftProfileP1, aLeftProfileP2 );
400   Handle(Geom_Line) aLeftProfileLine = new Geom_Line( gp_Ax1( aLeftProfileP1, aLeftProfileVec ) );
401   // Right profile line
402   aRightProfileP1.SetZ( 0 );
403   aRightProfileP2.SetZ( 0 );
404   gp_Vec aRightProfileVec( aRightProfileP1, aRightProfileP2 );
405   Handle(Geom_Line) aRightProfileLine = new Geom_Line( gp_Ax1( aRightProfileP1, aRightProfileVec ) );
406   // The point projections on the left and right profiles
407   GeomAPI_ProjectPointOnCurve aLeftProfileProject( aPointToTest, aLeftProfileLine );
408   GeomAPI_ProjectPointOnCurve aRightProfileProject( aPointToTest, aRightProfileLine );
409   // The point distance to the left and right profiles
410   Standard_Real aLeftProfileDist = aLeftProfileProject.LowerDistance();
411   Standard_Real aRightProfileDist = aRightProfileProject.LowerDistance();
412   // The coefficients
413   Standard_Real aFirstCoeff = aLeftProfileDist / ( aLeftProfileDist + aRightProfileDist );
414   Standard_Real aSecCoeff = aRightProfileDist / ( aLeftProfileDist + aRightProfileDist );
415
416   aResAltitude = aLeftAlt * ( 1 - aFirstCoeff ) + aRightAlt * ( 1 - aSecCoeff );
417   DEBTRACE("aResAltitude=" << aResAltitude);
418   return aResAltitude;
419 }
420
421
422
423