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