Salome HOME
refs #651: wrapping land cover map by Python
[modules/hydro.git] / src / HYDROData / HYDROData_Projection.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_Projection.h>
20 #include <BRepAdaptor_Curve.hxx>
21 #include <BRep_Builder.hxx>
22 #include <BRepBuilderAPI_Sewing.hxx>
23 #include <BRepLib.hxx>
24 #include <BRepLib_MakeEdge.hxx>
25 #include <BRepLib_MakeFace.hxx>
26 #include <BRepExtrema_DistShapeShape.hxx>
27 #include <Bnd_Box.hxx>
28 #include <BRepBndLib.hxx>
29 #include <gp_Pnt.hxx>
30 #include <gp_Lin.hxx>
31 #include <HLRAlgo_Projector.hxx>
32 #include <HLRBRep_Algo.hxx>
33 #include <HLRBRep_HLRToShape.hxx>
34 #include <IntCurvesFace_ShapeIntersector.hxx>
35 #include <Precision.hxx>
36 #include <TopoDS_Edge.hxx>
37 #include <TopoDS_Shape.hxx>
38 #include <TopoDS_Vertex.hxx>
39 #include <TopoDS_Wire.hxx>
40 #include <TopoDS.hxx>
41 #include <TopoDS_Iterator.hxx>
42 #include <TopExp.hxx>
43 #include <TopTools_ListIteratorOfListOfShape.hxx>
44 #include <TopTools_MapOfShape.hxx>
45 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
46 #include <TopTools_DataMapOfShapeReal.hxx>
47
48 HYDROData_Make3dMesh::HYDROData_Make3dMesh( const TopoDS_Shape& aShape,
49                                             const Standard_Real Tolerance )
50 {
51   myIntersector.Load(aShape, Tolerance);
52 }
53
54 Standard_Boolean HYDROData_Make3dMesh::GetHighestOriginal(const Standard_Real aX,
55                                                           const Standard_Real aY,
56                                                           gp_Pnt& HighestPoint)
57 {
58   gp_Pnt ProjPnt(aX, aY, 0.);
59   gp_Lin Ray(ProjPnt, -gp::DZ());
60   myIntersector.PerformNearest(Ray, -Precision::Infinite(), 0.);
61   if (!myIntersector.IsDone() || myIntersector.NbPnt() == 0)
62     return Standard_False;
63
64   HighestPoint = myIntersector.Pnt(1);
65   return Standard_True;
66 }
67
68 //=======================================================================
69 //function : FindClosestDirection
70 //purpose  : auxiliary function
71 //=======================================================================
72 TopoDS_Edge HYDROData_Projection::FindClosestDirection( const TopoDS_Vertex& StartVertex,
73                                                         const TopTools_ListOfShape& Candidates,
74                                                         const gp_Vec& StartDir )
75 {
76   Standard_Real MinAngle = M_PI;
77   TopoDS_Edge theEdge;
78
79   TopTools_ListIteratorOfListOfShape itl(Candidates);
80   for (; itl.More(); itl.Next())
81   {
82     TopoDS_Edge anEdge = TopoDS::Edge(itl.Value());
83     BRepAdaptor_Curve BAcurve(anEdge);
84     gp_Pnt aPnt, aPnt2;
85     gp_Vec aDir, aDir2;
86     
87     Standard_Real Par, Par2 = 0.; //to avoid warning
88     Standard_Boolean ToReverse = Standard_False;
89     
90     TopoDS_Vertex V1, V2;
91     TopExp::Vertices(anEdge, V1, V2);
92     if (StartVertex.IsSame(V1))
93       Par = 0.99*BAcurve.FirstParameter() + 0.01*BAcurve.LastParameter();
94     else
95     {
96       Par = 0.01*BAcurve.FirstParameter() + 0.99*BAcurve.LastParameter();
97       ToReverse = Standard_True;
98     }
99     if (V1.IsSame(V2))
100     {
101       Par2 = 0.01*BAcurve.FirstParameter() + 0.99*BAcurve.LastParameter();
102     }
103     
104     BAcurve.D1(Par, aPnt, aDir);
105     if (ToReverse)
106       aDir.Reverse();
107     
108     if (V1.IsSame(V2))
109     {
110       BAcurve.D1(Par2, aPnt2, aDir2);
111       aDir2.Reverse();
112     }
113     
114     Standard_Real anAngle = StartDir.Angle(aDir);
115     if (anAngle < MinAngle)
116     {
117       MinAngle = anAngle;
118       theEdge = anEdge;
119     }
120     
121     if (V1.IsSame(V2))
122     {
123       anAngle = StartDir.Angle(aDir2);
124       if (anAngle < MinAngle)
125       {
126         MinAngle = anAngle;
127         theEdge = anEdge;
128       }
129     }     
130   }
131
132   return theEdge;
133 }
134
135 //=======================================================================
136 //function : BuildOutWire
137 //purpose  : auxiliary function
138 //=======================================================================
139 TopoDS_Wire HYDROData_Projection::BuildOutWire( const TopTools_IndexedDataMapOfShapeListOfShape& VEmap,
140                                                 const TopoDS_Vertex& StartVertex,
141                                                 TopoDS_Edge&   StartEdge)
142 {
143   Standard_Real AngularTol = 0.1; //0.01; //1.e-4; //Precision::Angular();
144   TopTools_MapOfShape BoundEdges;
145
146   TopoDS_Wire theWire;
147   BRep_Builder BB;
148   BB.MakeWire(theWire);
149
150   TopoDS_Vertex V1, V2;
151   TopExp::Vertices(StartEdge, V1, V2);
152   if (StartVertex.IsSame(V1))
153     StartEdge.Orientation(TopAbs_FORWARD);
154   else
155     StartEdge.Orientation(TopAbs_REVERSED);
156   BB.Add(theWire, StartEdge);
157   BoundEdges.Add(StartEdge);
158
159   TopoDS_Vertex CurVertex = (StartVertex.IsSame(V1))? V2 : V1;
160   TopoDS_Edge CurEdge = StartEdge;
161   gp_Vec Normal(0., 0., 1.);
162   for (;;)
163   {
164     if (CurVertex.IsSame(StartVertex))
165     {
166       theWire.Closed(Standard_True);
167       break;
168     }
169     
170     const TopTools_ListOfShape& Candidates = VEmap.FindFromKey(CurVertex);
171     Standard_Real MinAngle = M_PI;
172     
173     BRepAdaptor_Curve CurCurve(CurEdge);
174     TopExp::Vertices(CurEdge, V1, V2);
175     Standard_Real CurPar;
176     if (CurVertex.IsSame(V2))
177       CurPar = CurCurve.LastParameter();
178     else
179       CurPar = CurCurve.FirstParameter();
180     gp_Pnt aPnt;
181     gp_Vec CurDir;
182     CurCurve.D1(CurPar, aPnt, CurDir);
183     if (CurDir.SquareMagnitude() < Precision::PConfusion())
184     {
185       CurPar = (CurVertex.IsSame(V2))?
186         0.01*CurCurve.FirstParameter() + 0.99*CurCurve.LastParameter() :
187         0.99*CurCurve.FirstParameter() + 0.01*CurCurve.LastParameter();
188       CurCurve.D1(CurPar, aPnt, CurDir);
189     }
190     if (CurVertex.IsSame(V1))
191       CurDir.Reverse();
192     
193     TopoDS_Edge theEdge;
194     TopTools_ListIteratorOfListOfShape itl(Candidates);
195     for (; itl.More(); itl.Next())
196     {
197       TopoDS_Edge anEdge = TopoDS::Edge(itl.Value());
198       if (BoundEdges.Contains(anEdge)) //if (anEdge.IsSame(CurEdge))
199         continue;
200       
201       TopExp::Vertices(anEdge, V1, V2);
202       if (V1.IsSame(V2))
203         continue;
204       
205       BRepAdaptor_Curve BAcurve(anEdge);
206       gp_Vec aDir;
207       Standard_Real aPar = (CurVertex.IsSame(V1))?
208         0.99*BAcurve.FirstParameter() + 0.01*BAcurve.LastParameter() :
209         0.01*BAcurve.FirstParameter() + 0.99*BAcurve.LastParameter();
210         //BAcurve.FirstParameter() : BAcurve.LastParameter();
211       
212       BAcurve.D1(aPar, aPnt, aDir);
213       if (CurVertex.IsSame(V2))
214         aDir.Reverse();
215       
216       Standard_Real anAngle = CurDir.AngleWithRef(aDir, Normal);
217       if (Abs(anAngle + M_PI) < AngularTol) //anAngle = -PI
218         continue;
219       if (anAngle < MinAngle)
220       {
221         MinAngle = anAngle;
222         theEdge = anEdge;
223       }
224     }
225     if (theEdge.IsNull())
226       break;
227     
228     TopExp::Vertices(theEdge, V1, V2);
229     if (CurVertex.IsSame(V1))
230       theEdge.Orientation(TopAbs_FORWARD);
231     else
232       theEdge.Orientation(TopAbs_REVERSED);
233     BB.Add(theWire, theEdge);
234     BoundEdges.Add(theEdge);
235     
236     CurVertex = (CurVertex.IsSame(V1))? V2 : V1;
237     CurEdge = theEdge;
238   }
239   
240   return theWire;
241 }
242
243 //=======================================================================
244 //Function : MakeProjection
245 //purpose  : this function makes the "shadow" of a shape on the plane XOY
246 //=======================================================================
247 TopoDS_Face HYDROData_Projection::MakeProjection( const TopoDS_Shape& aShape )
248 {
249   if ( aShape.IsNull() )
250     return TopoDS_Face();
251
252   HLRAlgo_Projector theProjector(gp::XOY());
253
254   Handle(HLRBRep_Algo) aHLRAlgo = new HLRBRep_Algo();
255   aHLRAlgo->Add(aShape);
256   aHLRAlgo->Projector(theProjector);
257   aHLRAlgo->Update();
258   aHLRAlgo->Hide();
259   HLRBRep_HLRToShape aHLRToShape(aHLRAlgo);
260
261   TopoDS_Shape SharpEdges, OutLines;
262   TopoDS_Compound Total;
263   BRep_Builder BB;
264   BB.MakeCompound(Total);
265
266   SharpEdges = aHLRToShape.VCompound();
267   OutLines = aHLRToShape.OutLineVCompound();
268
269   TopoDS_Iterator itc;
270   if (!SharpEdges.IsNull())
271   {
272     for (itc.Initialize(SharpEdges); itc.More(); itc.Next())
273       BB.Add(Total, itc.Value());
274   }
275   if (!OutLines.IsNull())
276   {
277     for (itc.Initialize(OutLines); itc.More(); itc.Next())
278       BB.Add(Total, itc.Value());
279   }
280
281   BRepBuilderAPI_Sewing aSewing;
282   Standard_Real tol = 1.0e-06;
283   Standard_Boolean NonManifoldMode = Standard_False;
284   aSewing.Init(tol, Standard_True,Standard_True,Standard_True,NonManifoldMode);
285   aSewing.SetFloatingEdgesMode(Standard_True);
286
287   aSewing.Add( Total );
288   aSewing.Perform();
289
290   TopoDS_Shape SewedEdges = aSewing.SewedShape();
291
292   BRepLib::BuildCurves3d( SewedEdges );
293
294   Bnd_Box theBox;
295   BRepBndLib::Add(SewedEdges, theBox);
296   Standard_Real Xmin, Ymin, Zmin, Xmax, Ymax, Zmax;
297   theBox.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
298
299   gp_Pnt P1(Xmin, Ymin, 0.), P2(Xmax, Ymin, 0.);
300   TopoDS_Edge LowerBound = BRepLib_MakeEdge(P1, P2);
301   
302   Standard_Real Deviation = 1.e-7; //1.e-5;
303   BRepExtrema_DistShapeShape DSS(SewedEdges, LowerBound, Deviation);
304
305   TopTools_IndexedDataMapOfShapeListOfShape StartShapes;
306   TopTools_DataMapOfShapeReal EdgesWithPars;
307   Standard_Integer i;
308   for (i = 1; i <= DSS.NbSolution(); i++)
309   {
310     if (!StartShapes.Contains( DSS.SupportOnShape1(i) ))
311     {
312       TopTools_ListOfShape aList;
313       aList.Append( DSS.SupportOnShape2(i) );
314       StartShapes.Add( DSS.SupportOnShape1(i), aList );
315     }
316     else
317     {
318       TopTools_ListOfShape& aList = StartShapes.ChangeFromKey( DSS.SupportOnShape1(i) );
319       aList.Append( DSS.SupportOnShape2(i) );
320     }
321     if (DSS.SupportTypeShape1(i) == BRepExtrema_IsOnEdge)
322     {
323       Standard_Real aPar;
324       DSS.ParOnEdgeS1(i, aPar);
325       EdgesWithPars.Bind( DSS.SupportOnShape1(i), aPar );
326     }
327   }
328   
329   TopoDS_Shape aStartShape;
330   for (i = 1; i <= StartShapes.Extent(); i++)
331   {
332     const TopoDS_Shape& aShape = StartShapes.FindKey(i);
333     if (aShape.ShapeType() == TopAbs_EDGE)
334     {
335       aStartShape = aShape;
336       break;
337     }
338   }
339   if (aStartShape.IsNull())
340     aStartShape = StartShapes.FindKey(1); //it is a vertex
341   
342   TopoDS_Vertex StartVertex;
343   TopoDS_Edge   StartEdge;
344   gp_Vec StartDir(1.,0.,0.);
345   TopTools_IndexedDataMapOfShapeListOfShape VEmap;
346   TopExp::MapShapesAndAncestors(SewedEdges, TopAbs_VERTEX, TopAbs_EDGE, VEmap);
347   if (aStartShape.ShapeType() == TopAbs_EDGE)
348   {
349     StartEdge = TopoDS::Edge(aStartShape);
350     Standard_Real StartPar = EdgesWithPars(StartEdge);
351     BRepAdaptor_Curve BAcurve(StartEdge);
352     gp_Pnt aPnt;
353     gp_Vec aDir;
354     BAcurve.D1(StartPar, aPnt, aDir);
355     StartVertex = (StartDir * aDir > 0.)?
356       TopExp::FirstVertex(StartEdge) : TopExp::LastVertex(StartEdge);
357   }
358   else // it is a vertex
359   {
360     StartVertex = TopoDS::Vertex(aStartShape);
361     const TopTools_ListOfShape& Candidates = VEmap.FindFromKey(StartVertex);
362     StartEdge = FindClosestDirection(StartVertex, Candidates, StartDir);
363   }
364
365   TopoDS_Wire OutWire = BuildOutWire(VEmap, StartVertex, StartEdge);
366   TopoDS_Face ProjFace = BRepLib_MakeFace(OutWire, Standard_True); //only plane
367   
368   return ProjFace;
369 }
370