]> SALOME platform Git repositories - modules/hydro.git/blob - src/HYDROData/HYDROData_DTM.cxx
Salome HOME
Merge branch 'BR_2017_PORTING' into BR_2018_V8_5
[modules/hydro.git] / src / HYDROData / HYDROData_DTM.cxx
1
2 #include <HYDROData_DTM.h>
3 #include <HYDROData_Profile.h>
4
5 #include <Geom2d_BSplineCurve.hxx>
6 #include <Geom2dAPI_Interpolate.hxx>
7 #include <TColgp_HArray1OfPnt2d.hxx>
8 #include <TColgp_Array1OfVec2d.hxx>
9 #include <TColStd_HArray1OfBoolean.hxx>
10 #include <TopoDS.hxx>
11 #include <TopoDS_Edge.hxx>
12 #include <TopoDS_Wire.hxx>
13 #include <TopExp_Explorer.hxx>
14 #include <BRep_Tool.hxx>
15 #include <gp_Ax3.hxx>
16 #include <Geom_Line.hxx>
17 #include <Geom2d_Line.hxx>
18 #include <Geom2d_TrimmedCurve.hxx>
19 #include <Geom_BSplineCurve.hxx>
20 #include <Geom2d_BSplineCurve.hxx>
21 #include <GeomAPI_Interpolate.hxx>
22 #include <TColStd_Array1OfReal.hxx>
23 #include <TColStd_Array1OfInteger.hxx>
24 #include <TColgp_Array1OfPnt.hxx>
25 #include <TColgp_Array1OfVec.hxx>
26 #include <TColgp_HArray1OfPnt.hxx>
27 #include <Geom2dAPI_InterCurveCurve.hxx>
28 #include <Geom2dAPI_ProjectPointOnCurve.hxx>
29 #include <Geom2dAdaptor_Curve.hxx>
30 #include <GCPnts_AbscissaPoint.hxx>
31 #include <BRepBuilderAPI_MakeEdge.hxx>
32 #include <limits>
33
34 #include <BRepLib_MakeEdge.hxx>
35 #include <BRepLib_MakeWire.hxx>
36 #include <BRep_Builder.hxx>
37 #include <GeomProjLib.hxx>
38 #include <Geom_TrimmedCurve.hxx>
39 #include <BRepTools_WireExplorer.hxx>
40 #include <TopTools_IndexedMapOfShape.hxx>
41 #include <BRepBuilderAPI_MakeFace.hxx>
42 #include <TopExp.hxx>
43 #include <TopTools_IndexedMapOfOrientedShape.hxx>
44 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
45
46 #include <BRepLib_MakeEdge.hxx>
47 #include <BRepLib_MakeWire.hxx>
48 #include <BRep_Builder.hxx>
49 #include <ShapeAnalysis_Wire.hxx>
50 #include <BRepAlgo_NormalProjection.hxx>
51 #include <ShapeUpgrade_UnifySameDomain.hxx>
52 #include <BRepBuilderAPI_MakePolygon.hxx>
53 #include <BOPAlgo_Builder.hxx>
54 #include <BRepAdaptor_Curve.hxx>
55 #include <GeomProjLib.hxx>
56 #include <gp_Pln.hxx>
57 #include <TopTools_IndexedDataMapOfShapeShape.hxx>
58 #include <TopTools_ListIteratorOfListOfShape.hxx>
59 #include <TopTools_SequenceOfShape.hxx>
60 #include <assert.h>
61 #include <float.h>
62
63 IMPLEMENT_STANDARD_RTTIEXT( HYDROData_DTM, HYDROData_Bathymetry )
64
65 HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv, const gp_Vec2d& theProfileDir, double theDeltaZ, double theMaxZ )
66   : myXcurv( theXCurv ), myProfileDir( theProfileDir ), myDeltaZ( theDeltaZ ), myMaxZ (theMaxZ)
67 {
68 }
69
70 HYDROData_DTM::CurveUZ::~CurveUZ()
71 {
72 }
73
74 double HYDROData_DTM::CurveUZ::Xcurv() const
75 {
76   return myXcurv;
77 }
78
79 gp_Vec2d HYDROData_DTM::CurveUZ::ProfileDir() const
80 {
81   return myProfileDir;
82 }
83
84 double HYDROData_DTM::CurveUZ::DeltaZ() const
85 {
86   return myDeltaZ;
87 }
88
89 double HYDROData_DTM::CurveUZ::MaxZ() const
90 {
91   return myMaxZ;
92 }
93
94 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator + ( const CurveUZ& c ) const
95 {
96   HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv(), ProfileDir() + c.ProfileDir(), DeltaZ() + c.DeltaZ(), MaxZ() + c.MaxZ() );
97   size_t n = size(), n1 = c.size();
98   if( n!=n1 )
99   {
100     std::cout << "Warning: different number of points in curves: " << n << ", " << n1 << std::endl;
101   }
102   int q = n < n1 ? n : n1;
103   res.reserve( q );
104   for( int i=0; i<q; i++ )
105   {
106     PointUZ p;
107     p.U = operator[]( i ).U + c[i].U;
108     p.Z = operator[]( i ).Z;
109     res.push_back( p );
110   }
111   return res;
112 }
113
114 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator * ( double d ) const
115 {
116   HYDROData_DTM::CurveUZ res( Xcurv()*d, ProfileDir()*d, DeltaZ()*d, MaxZ()*d );
117   size_t n = size();
118   res.reserve( n );
119   for( int i=0; i<n; i++ )
120   {
121     PointUZ p;
122     p.U = operator[]( i ).U * d;
123     p.Z = operator[]( i ).Z;
124     res.push_back( p );
125   }
126   return res;
127 }
128
129
130
131
132 HYDROData_DTM::HYDROData_DTM()
133 {
134 }
135
136 HYDROData_DTM::~HYDROData_DTM()
137 {
138 }
139
140 HYDROData_SequenceOfObjects HYDROData_DTM::GetProfiles() const
141 {
142   return GetReferenceObjects( DataTag_Profiles );
143 }
144
145 void HYDROData_DTM::SetProfiles( const HYDROData_SequenceOfObjects& theProfiles )
146 {
147   SetReferenceObjects( theProfiles, DataTag_Profiles );
148   Changed( Geom_3d );
149 }
150
151 double HYDROData_DTM::GetDDZ() const
152 {
153   return GetDouble( DataTag_DDZ );
154 }
155
156 void HYDROData_DTM::SetDDZ( double theDDZ )
157 {
158   SetDouble( DataTag_DDZ, theDDZ );
159   Changed( Geom_3d );
160 }
161
162 double HYDROData_DTM::GetSpatialStep() const
163 {
164   return GetDouble( DataTag_SpatialStep );
165 }
166
167 void HYDROData_DTM::SetSpatialStep( double theSpatialStep )
168 {
169   SetDouble( DataTag_SpatialStep, theSpatialStep );
170   Changed( Geom_3d );
171 }
172
173 void HYDROData_DTM::PointsToWire(const AltitudePoints& pnts, TopoDS_Wire& W )
174 {
175  
176   BRepBuilderAPI_MakePolygon PM;
177   for (int i = 0; i < pnts.size(); i++)
178     PM.Add(gp_Pnt(pnts[i].X, pnts[i].Y, pnts[i].Z));
179   
180   W = PM.Wire();
181 }
182
183 void HYDROData_DTM::PointsToEdge(const AltitudePoints& pnts, TopoDS_Edge& E )
184
185   Handle(TColgp_HArray1OfPnt) gpPoints = new TColgp_HArray1OfPnt(1, (int)pnts.size());
186
187   for (int i = 0; i < pnts.size(); i++)
188     gpPoints->SetValue(i+1, gp_Pnt(pnts[i].X, pnts[i].Y, pnts[i].Z));
189
190   GeomAPI_Interpolate anInterpolator(gpPoints, Standard_False,1.0e-6);
191   anInterpolator.Perform() ;
192   if (anInterpolator.IsDone()) 
193   {
194     Handle(Geom_BSplineCurve) C = anInterpolator.Curve();
195     E = BRepBuilderAPI_MakeEdge(C).Edge();
196   }
197 }
198
199 TopTools_IndexedMapOfOrientedShape HYDROData_DTM::Create3DShape(const AltitudePoints& left,
200                                                                 const AltitudePoints& right,
201                                                                 const std::vector<AltitudePoints>& main_profiles)
202 {  
203   TopTools_IndexedMapOfOrientedShape ll;
204   //TopoDS_Wire LWire, RWire;
205   //PointsToWire(left, LWire);
206   //PointsToWire(right, RWire);
207   TopoDS_Edge LEdge, REdge;
208   PointsToEdge(left, LEdge);
209   PointsToEdge(right, REdge);
210   if (!LEdge.IsNull())
211     ll.Add(LEdge.Oriented(TopAbs_FORWARD));
212
213   for (int k = 0; k < main_profiles.size(); k++)
214   {
215     TopoDS_Wire W;
216     PointsToWire(main_profiles[k], W);
217     TopAbs_Orientation Ori = TopAbs_INTERNAL;
218     if (k == 0 || k == main_profiles.size() - 1)
219       Ori = TopAbs_FORWARD;
220     ll.Add(W.Oriented(Ori));
221   }
222
223   if (!REdge.IsNull())
224     ll.Add(REdge.Oriented(TopAbs_FORWARD)); 
225   //yes, add subshapes in this order (left + profiles + right)
226   //otherwise the projected wire will be non-manifold
227
228   return ll;
229 }
230
231
232 void HYDROData_DTM::Update()
233 {
234   AltitudePoints points;
235   TopoDS_Shape Out3dPres;
236   TopoDS_Shape Out2dPres;
237   TopoDS_Shape OutLeftB;
238   TopoDS_Shape OutRightB;
239   TopoDS_Shape OutInlet;
240   TopoDS_Shape OutOutlet;
241
242   HYDROData_SequenceOfObjects objs = GetProfiles();  
243   double ddz = GetDDZ();
244   double step = GetSpatialStep();
245   std::set<int> InvInd;
246   bool WireIntersections; //__TODO
247   CreateProfilesFromDTM( objs, ddz, step, points, Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, true, true, InvInd, -1, WireIntersections );
248   SetAltitudePoints( points );  
249   
250   SetShape( DataTag_LeftBankShape, OutLeftB);
251   SetShape( DataTag_RightBankShape, OutRightB);
252   SetShape( DataTag_InletShape, OutInlet);
253   SetShape( DataTag_OutletShape, OutOutlet );
254   SetShape( DataTag_3DShape, Out3dPres );
255   SetShape( DataTag_2DShape, Out2dPres );
256
257   HYDROData_Bathymetry::Update();
258 }
259
260 void HYDROData_DTM::GetPresentationShapes( TopoDS_Shape& Out3dPres,
261                                            TopoDS_Shape& Out2dPres,
262                                            TopoDS_Shape& OutLeftB,
263                                            TopoDS_Shape& OutRightB,
264                                            TopoDS_Shape& OutInlet,
265                                            TopoDS_Shape& OutOutlet )
266 {
267   //without update!
268   OutLeftB = GetShape( DataTag_LeftBankShape);
269   OutRightB = GetShape( DataTag_RightBankShape);
270   OutInlet = GetShape( DataTag_InletShape);
271   OutOutlet = GetShape( DataTag_OutletShape );
272   Out3dPres = GetShape( DataTag_3DShape );
273   Out2dPres = GetShape( DataTag_2DShape );
274 }
275 void HYDROData_DTM::CreateProfilesFromDTM (const HYDROData_SequenceOfObjects& InpProfiles,
276                                            double ddz,
277                                            double step, 
278                                            AltitudePoints& points,
279                                            TopoDS_Shape& Out3dPres,
280                                            TopoDS_Shape& Out2dPres,
281                                            TopoDS_Shape& OutLeftB,
282                                            TopoDS_Shape& OutRightB,
283                                            TopoDS_Shape& OutInlet,
284                                            TopoDS_Shape& OutOutlet,
285                                            bool Create3dPres,
286                                            bool Create2dPres,
287                                            std::set<int>& InvInd,
288                                            int thePntsLimit,
289                                            bool& WireIntersections)
290 {
291   int aLower = InpProfiles.Lower(), anUpper = InpProfiles.Upper();
292   size_t n = anUpper - aLower + 1;
293
294   std::vector<Handle(HYDROData_Profile)> profiles;
295   profiles.reserve( n ); 
296   for( int i=aLower; i<=anUpper; i++ )
297   {
298     Handle(HYDROData_Profile) aProfile = Handle(HYDROData_Profile)::DownCast( InpProfiles.Value( i ) );
299     if( !aProfile.IsNull() )
300       profiles.push_back( aProfile );
301   }
302   const double EPS = 1E-3;
303   AltitudePoints left;
304   AltitudePoints right;
305   std::vector<AltitudePoints> main_profiles;
306
307   if( thePntsLimit > 0 )
308   {
309     int aNbPoints = EstimateNbPoints( profiles, ddz, step );
310     if( aNbPoints < 0 || aNbPoints > thePntsLimit )
311       return;
312   }
313
314   if( ddz>EPS && step>EPS )
315     CreateProfiles(profiles, ddz, step, left, right, points, main_profiles, 
316     Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, Create3dPres, Create2dPres, InvInd, WireIntersections );
317 }
318
319 bool HYDROData_DTM::GetPlanarFaceFromBanks( const TopoDS_Edge& LB, const TopoDS_Edge& RB, TopoDS_Face& outF,
320   TopTools_SequenceOfShape* Boundr)
321 {
322   BRep_Builder BB;
323   TopoDS_Face F;
324   Handle_Geom_Plane refpl = new Geom_Plane(gp_Pnt(0,0,0), gp_Dir(0,0,1));
325
326   TopoDS_Vertex VFI, VLI, VFO, VLO;
327   TopoDS_Edge prLB;
328   TopoDS_Edge prRB;
329
330
331   BRepAdaptor_Curve LBAD(LB);
332   Handle_Geom_Curve LBPC = GeomProjLib::ProjectOnPlane(LBAD.Curve().Curve(), refpl, gp_Dir(0, 0, -1), 1 );
333   prLB = BRepLib_MakeEdge(LBPC).Edge();
334
335   BRepAdaptor_Curve RBAD(RB);
336   Handle_Geom_Curve RBPC = GeomProjLib::ProjectOnPlane(RBAD.Curve().Curve(), refpl, gp_Dir(0, 0, -1), 1 );
337   prRB = BRepLib_MakeEdge(RBPC).Edge();
338
339   TopExp::Vertices(prLB, VFI, VFO, 1);
340   TopExp::Vertices(prRB, VLI, VLO, 1);
341   TopoDS_Edge prIL = BRepLib_MakeEdge(VFI, VLI).Edge();
342   TopoDS_Edge prOL = BRepLib_MakeEdge(VFO, VLO).Edge();
343   TopoDS_Wire prW = BRepLib_MakeWire(prLB, prIL, prOL, prRB).Wire();
344   outF = BRepBuilderAPI_MakeFace(refpl->Pln(), prW, 1).Face();
345
346   if (Boundr)
347   {
348     Boundr->Append(prLB);
349     Boundr->Append(prIL);
350     Boundr->Append(prOL);
351     Boundr->Append(prRB);
352   }
353
354   ShapeAnalysis_Wire WA(prW, outF, Precision::Confusion());
355   bool res = WA.CheckSelfIntersection();
356   return !res;
357 }
358
359 void HYDROData_DTM::CreateProfiles(const std::vector<Handle(HYDROData_Profile)>& theProfiles,
360                                    double theDDZ,
361                                    double theSpatialStep,
362                                    AltitudePoints& theOutLeft,
363                                    AltitudePoints& theOutRight,
364                                    AltitudePoints& theOutPoints,
365                                    std::vector<AltitudePoints>& theOutMainProfiles,
366                                    TopoDS_Shape& Out3dPres,
367                                    TopoDS_Shape& Out2dPres,
368                                    TopoDS_Shape& OutLeftB,
369                                    TopoDS_Shape& OutRightB,
370                                    TopoDS_Shape& OutInlet,
371                                    TopoDS_Shape& OutOutlet,
372                                    bool Create3dPres,
373                                    bool Create2dPres,
374                                    std::set<int>& InvInd,
375                                    bool& ProjStat)
376 {
377   if (theProfiles.empty())
378     return;
379   theOutPoints = Interpolate( theProfiles, theDDZ, theSpatialStep, theOutLeft, theOutRight, theOutMainProfiles, InvInd );
380   //note that if Create3dPres is false => Create2dPres flag is meaningless!
381   if( theOutPoints.empty() )
382     return;
383
384   if (Create3dPres)
385   {
386     TopTools_IndexedMapOfOrientedShape ll = Create3DShape( theOutLeft, theOutRight, theOutMainProfiles);
387     
388     if (ll.IsEmpty())
389       return;
390     BRep_Builder BB;
391     TopoDS_Compound cmp;
392     BB.MakeCompound(cmp);
393     for (int i = 1; i <= ll.Extent(); i++)
394       BB.Add(cmp, ll(i));
395
396     Out3dPres = cmp;
397
398     //same order as in HYDROData_DTM::Update()
399     OutLeftB = ll(1);
400     OutRightB = ll(ll.Extent());
401     OutInlet = ll(2);
402     OutOutlet = ll(ll.Extent() - 1);
403
404     if (Create2dPres)
405     {
406       TopoDS_Face outF;
407       ProjStat = GetPlanarFaceFromBanks(TopoDS::Edge(OutLeftB), TopoDS::Edge(OutRightB), outF, NULL);
408       Out2dPres = outF;
409     };
410   }
411 }
412
413
414
415
416 void HYDROData_DTM::GetProperties( const Handle(HYDROData_Profile)& theProfile,
417                     gp_Pnt& theLowestPoint, gp_Vec2d& theDir,
418                     double& theZMin, double& theZMax )
419 {
420   theLowestPoint = theProfile->GetBottomPoint();
421   
422   gp_XY aLeft, aRight;
423   theProfile->GetLeftPoint( aLeft, true, true );
424   theProfile->GetRightPoint( aRight, true, true );
425   double x = aRight.X()-aLeft.X();
426   double y = aRight.Y()-aLeft.Y();
427   theDir = gp_Vec2d( x, y );
428
429   HYDROData_Profile::ProfilePoints points = theProfile->GetProfilePoints();
430   int lo = points.Lower();
431   int up = points.Upper();
432   theZMin = DBL_MAX;
433   theZMax = -theZMin;
434   for( int i=lo; i<=up; i++ )
435   {
436     double z = points.Value( i ).Z();
437     if( z>theZMax )
438       theZMax = z;
439     if( z<theZMin )
440       theZMin = z;
441   }
442 }
443
444 inline gp_Pnt2d To2D( const gp_Pnt& thePnt, const gp_Trsf& theTr,
445                       double& theUMin, double& theUMax )
446 {
447   gp_Pnt p = thePnt.Transformed( theTr );
448   double u = p.X();
449   double z = p.Z();
450   if( u<theUMin )
451     theUMin = u;
452   if( u>theUMax )
453     theUMax = u;
454   return gp_Pnt2d( u, z );
455 }
456
457 Handle(TColgp_HArray1OfPnt2d) To2D( const TColgp_Array1OfPnt& thePoints,
458                                     const gp_Trsf& theTr,
459                                     double& theUMin, double& theUMax )
460 {
461   int low = thePoints.Lower(), up = thePoints.Upper();
462   Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( low, up );
463   for( int i=low; i<=up; i++ )
464     points->SetValue( i, To2D( thePoints.Value( i ), theTr, theUMin, theUMax ) );
465   return points;
466 }
467
468 Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve, 
469                                 Standard_Real theFirst, Standard_Real theLast,
470                                 const gp_Trsf& theTr,
471                                 double& theUMin, double& theUMax )
472 {
473   if( theCurve->IsKind( STANDARD_TYPE( Geom_Line ) ) )
474   {
475     gp_Pnt aFirstPnt, aLastPnt;
476     theCurve->D0( theFirst, aFirstPnt );
477     theCurve->D0( theLast, aLastPnt );
478
479     gp_Pnt2d
480       aFirst2d = To2D( aFirstPnt, theTr, theUMin, theUMax ),
481       aLast2d = To2D( aLastPnt, theTr, theUMin, theUMax );
482
483     gp_Vec2d dir( aFirst2d, aLast2d );
484     Handle(Geom2d_Line) aLine2d = new Geom2d_Line( aFirst2d, gp_Dir2d( dir.X(), dir.Y() ) );
485     return new Geom2d_TrimmedCurve( aLine2d, 0, aLast2d.Distance( aFirst2d ) );
486   }
487
488   if( theCurve->IsKind( STANDARD_TYPE( Geom_BSplineCurve ) ) )
489   {
490     Handle(Geom_BSplineCurve) aSpline = Handle(Geom_BSplineCurve)::DownCast( theCurve );
491
492     Handle(TColgp_HArray1OfPnt2d) poles = To2D( aSpline->Poles(), theTr, theUMin, theUMax );
493     const TColStd_Array1OfReal& knots = aSpline->Knots();
494     const TColStd_Array1OfInteger& multiplicities = aSpline->Multiplicities();
495     int aDegree = aSpline->Degree();
496
497     return new Geom2d_BSplineCurve( poles->Array1(), knots, multiplicities, aDegree );
498   }
499
500   return Handle(Geom2d_Curve)();
501 }
502
503 #include <GCE2d_MakeSegment.hxx>
504 #include <Geom2dAPI_InterCurveCurve.hxx>
505 bool IsCooriented( const Handle(HYDROData_Profile)& theProfile1,
506                    const Handle(HYDROData_Profile)& theProfile2 ) 
507 {
508   if( theProfile1==theProfile2 )
509     return true;
510   
511   gp_XY lp1, rp1, lp2, rp2;
512   theProfile1->GetLeftPoint(lp1);
513   theProfile1->GetRightPoint(rp1);
514   theProfile2->GetLeftPoint(lp2);
515   theProfile2->GetRightPoint(rp2);
516
517   GCE2d_MakeSegment s1(lp1, lp2);
518   GCE2d_MakeSegment s2(rp1, rp2);
519
520   Geom2dAPI_InterCurveCurve inter;
521   inter.Init(s1.Value(), s2.Value());
522   if (inter.NbPoints() == 0)
523     return true;
524   else
525     return false;
526 }
527
528 Handle(Geom2d_BSplineCurve) HYDROData_DTM::CreateHydraulicAxis( 
529   const std::vector<Handle(HYDROData_Profile)>& theProfiles,
530   std::vector<double>& theDistances )
531 {
532   size_t n = theProfiles.size();
533   if( n==1 )
534     return Handle_Geom2d_BSplineCurve();
535
536   Handle_Geom2d_BSplineCurve aResult;
537
538   Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( 1, (int)n );
539   TColgp_Array1OfVec2d tangents( 1, (int)n );
540   Handle(TColStd_HArray1OfBoolean) flags = new TColStd_HArray1OfBoolean( 1, (int)n );
541
542   // Stage 1. Orient all profiles to be co-oriented with the first profile
543   theProfiles[0]->Update();
544   for( size_t i = 1; i < n; i++ )
545   {
546     Handle(HYDROData_Profile) aProfile = theProfiles[i];
547     Handle(HYDROData_Profile) aPrevProfile = theProfiles[i-1];
548
549     if( !IsCooriented( aProfile, aPrevProfile ) )
550     {
551       gp_XY lp, rp;
552       aProfile->GetLeftPoint( lp, true );
553       aProfile->GetRightPoint( rp, true );
554       aProfile->SetLeftPoint( rp, true );
555       aProfile->SetRightPoint( lp, true );
556     }
557     aProfile->Update();
558   }
559
560   // Stage 2. Calculate normals so that each normal "points" to the next profile
561   for( size_t i = 0; i < n; i++ )
562   {
563     Handle(HYDROData_Profile) aProfile = theProfiles[i];
564     Handle(HYDROData_Profile) aNextProfile = i==n-1 ? theProfiles[i-1] : theProfiles[i+1];
565
566     gp_Pnt aLowest;
567     gp_Vec2d aNormal;
568     double zmin, zmax;
569
570     gp_XYZ curP = aProfile->GetBottomPoint(true);
571     gp_XY curP2d = gp_XY(curP.X(), curP.Y());
572
573     gp_XYZ nextP = aNextProfile->GetBottomPoint(true);
574     gp_XY nextP2d = gp_XY(nextP.X(), nextP.Y());
575
576     gp_Vec2d aPrTangent;
577     GetProperties( aProfile, aLowest, aPrTangent, zmin, zmax );
578     aNormal.SetCoord( -aPrTangent.Y(), aPrTangent.X() );
579
580     gp_Vec2d aDirToNextProfile(nextP2d.X() - curP2d.X(), nextP2d.Y() - curP2d.Y() );
581     if( i==n-1 )
582       aDirToNextProfile.Reverse();
583
584     if (aNormal.Dot(aDirToNextProfile) < 0)
585       aNormal.Reverse();
586
587     aNormal.Normalize();
588
589     points->SetValue( (int)(i+1), gp_Pnt2d( aLowest.X(), aLowest.Y() ) );
590     tangents.SetValue( (int)(i+1), aNormal );
591     flags->SetValue( (int)(i+1), Standard_True );
592   }
593
594   Geom2dAPI_Interpolate anInterpolator( points, Standard_False, Standard_False );
595   anInterpolator.Load( tangents, flags );
596   anInterpolator.Perform();
597   if( anInterpolator.IsDone() )
598   {
599     aResult = anInterpolator.Curve();
600
601     //fill the distances vector
602     Geom2dAdaptor_Curve anAdaptor( aResult );
603
604     theDistances.clear();
605     theDistances.reserve( n );
606     Standard_Real aParamFirst = anAdaptor.FirstParameter(), aParamLast = anAdaptor.LastParameter();
607     for( size_t i = 1; i <= n; i++ )
608     {
609       gp_Pnt2d aPnt = points->Value( (Standard_Integer)i );
610       Geom2dAPI_ProjectPointOnCurve aProject( aPnt, aResult );
611       Standard_Real aParam = aProject.LowerDistanceParameter();
612       double aDistance = GCPnts_AbscissaPoint::Length( anAdaptor, aParamFirst, aParam );  
613       theDistances.push_back( aDistance );
614     }
615   }
616   return aResult;
617 }
618
619 std::vector<Handle(Geom2d_Curve)> HYDROData_DTM::ProfileToParametric( 
620   const Handle(HYDROData_Profile)& theProfile,
621   double& theUMin, double& theUMax, gp_Vec2d& theDir )
622 {
623   std::vector<Handle(Geom2d_Curve)> curves;
624   
625   // Transformation of the coordinate systems
626   gp_Pnt aLowest;
627   double zmin, zmax;
628   GetProperties( theProfile, aLowest, theDir, zmin, zmax );
629
630   gp_Ax3 aStd3d( gp_Pnt( 0, 0, 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( 1, 0, 0 ) );
631   gp_Ax3 aLocal( gp_Pnt( aLowest.X(), aLowest.Y(), 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( theDir.X(), theDir.Y(), 0 ) );
632
633   gp_Trsf aTransf;
634   aTransf.SetTransformation( aStd3d, aLocal );
635
636   // Iteration via edges
637   TopoDS_Wire aWire = TopoDS::Wire( theProfile->GetShape3D() );
638   TopExp_Explorer anExp( aWire, TopAbs_EDGE );
639   for( ; anExp.More(); anExp.Next() )
640   {
641     // Extract an edge from wire
642     TopoDS_Edge anEdge = TopoDS::Edge( anExp.Current() );
643
644     // Extract a curve corresponding to the edge
645     TopLoc_Location aLoc;
646     Standard_Real aFirst, aLast;
647     Handle(Geom_Curve) aCurve = BRep_Tool::Curve( anEdge, aLoc, aFirst, aLast );
648
649     // Convert the curve to 2d CS
650     Handle(Geom2d_Curve) aCurve2d = CurveTo2D( aCurve, aFirst, aLast, aTransf, theUMin, theUMax );
651     if( !aCurve2d.IsNull() )
652       curves.push_back( aCurve2d );
653   }
654   return curves;
655 }
656
657
658 bool CalcMidWidth( const std::set<double>& intersections, double& theMid, double& theWid )
659 {
660   double umin = DBL_MAX,
661          umax = -umin;
662
663   size_t n = intersections.size();
664   if( n <= 0 )
665     return false;
666
667   std::set<double>::const_iterator it = intersections.begin(), last = intersections.end();
668   for( ; it!=last; it++ )
669   {
670     double u = *it;
671     if( u<umin )
672       umin = u;
673     if( u>umax )
674       umax = u;
675   }
676   theMid = ( umin+umax )/2;
677   theWid = umax-umin;
678   return true;
679 }
680
681 void HYDROData_DTM::ProfileDiscretization( const Handle(HYDROData_Profile)& theProfile, 
682                                            double theXCurv, double theMinZ, double theMaxZ, double theTopZ, double theDDZ,
683                                            CurveUZ& theMidPointCurve,
684                                            CurveUZ& theWidthCurve,                                           
685                                            int& intersection_nb,
686                                            double theTolerance)
687 {
688   double aDblMax = DBL_MAX,
689          aUMin = aDblMax,
690          aUMax = -aUMin,
691          aVMax = 1000000;
692   
693   gp_Vec2d aProfileDir;
694   std::vector<Handle(Geom2d_Curve)> curves = ProfileToParametric( theProfile, aUMin, aUMax, aProfileDir );
695   size_t n = curves.size();
696
697   if( n==0 )
698     return;
699
700   // we add the "virtual" vertical lines to simulate the intersection with profile 
701   gp_Pnt2d aFirst, aLast;
702   curves[0]->D0( curves[0]->FirstParameter(), aFirst );
703   curves[n-1]->D0( curves[n-1]->LastParameter(), aLast );
704   Handle(Geom2d_Line) aV1 = new Geom2d_Line( aFirst, gp_Dir2d( 0, 1 ) );
705   Handle(Geom2d_TrimmedCurve) aT1 = new Geom2d_TrimmedCurve( aV1, 0.0, aVMax );
706   
707   Handle(Geom2d_Line) aV2 = new Geom2d_Line( aLast, gp_Dir2d( 0, 1 ) );
708   Handle(Geom2d_TrimmedCurve) aT2 = new Geom2d_TrimmedCurve( aV2, 0.0, aVMax );
709
710   curves.push_back( aT1 );
711   curves.push_back( aT2 );
712   
713   int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 );
714   theMidPointCurve = CurveUZ( theXCurv, aProfileDir, theMinZ, theTopZ);
715   theMidPointCurve.reserve( psize );
716   theWidthCurve = CurveUZ( theXCurv, aProfileDir, theMinZ, theTopZ );
717   theWidthCurve.reserve( psize );
718
719   n = curves.size();
720   // for each discrete value of z we search intersection with profile
721   for( double z1 = theMinZ; z1 <= theMaxZ; z1 += theDDZ )
722   {
723     Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z1 ), gp_Dir2d( 1, 0 ) );
724     std::set<double> intersections;
725     for( size_t i = 0; i < n; i++ )
726     {
727       Handle(Geom2d_Curve) aCurve = curves[i];
728       Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
729       for( int k=1, m=anIntersect.NbPoints(); k<=m; k++ )
730         intersections.insert( anIntersect.Point( k ).X() );
731     }
732
733     intersection_nb = intersections.size();
734     if( intersection_nb >= 1 )
735     {
736       double u_mid, u_wid;
737       if( !CalcMidWidth( intersections, u_mid, u_wid ) )
738         continue;
739
740       double z = z1 - theMinZ;
741       PointUZ p_mid;
742       p_mid.U = u_mid;
743       p_mid.Z = z;
744       theMidPointCurve.push_back( p_mid );
745
746       PointUZ p_wid;
747       p_wid.U = u_wid;
748       p_wid.Z = z;
749       theWidthCurve.push_back( p_wid );
750     }
751   }
752 }
753
754 void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCurveB, 
755                                  int theNbSteps, std::vector<CurveUZ>& theInterpolation,
756                                  bool isAddSecond )
757 {
758   theInterpolation.clear();
759   int d = isAddSecond ? 2 : 1;
760   theInterpolation.reserve( theNbSteps+d );
761   double dt = 1.0 / double( theNbSteps + 1 );
762   double t = dt;
763   theInterpolation.push_back( theCurveA );
764   for( int i=0; i<theNbSteps; i++, t+=dt )
765   {
766     CurveUZ anInterp = theCurveA*(1-t) + theCurveB*t;
767     theInterpolation.push_back( anInterp );
768   }
769   if( isAddSecond )
770     theInterpolation.push_back( theCurveB );
771 }
772 #include <BRepLib_MakeEdge2d.hxx>
773 void HYDROData_DTM::CurveTo3D( const Handle(Geom2d_BSplineCurve)& theHydraulicAxis,
774                                const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
775                                AltitudePoints& thePoints )
776 {
777   Geom2dAdaptor_Curve anAdaptor( theHydraulicAxis );
778   TopoDS_Edge E2d = BRepLib_MakeEdge2d(theHydraulicAxis).Edge();
779   GCPnts_AbscissaPoint ap( anAdaptor, theMidCurve.Xcurv(), anAdaptor.FirstParameter() );  
780   double aParam = ap.Parameter();
781
782   gp_Pnt2d point;
783   anAdaptor.D0( aParam, point );
784   gp_Vec2d profile_dir = theMidCurve.ProfileDir();
785   //gp_Dir tangent_n( -profile_dir.Y(), profile_dir.X(), dz );
786   profile_dir.Normalize();
787   
788   size_t n = theMidCurve.size();
789   std::map<double, AltitudePoint> sorted_points;
790   bool isOnTop = false;
791   for( size_t i=0; i<n; i++ ) // build the two banks of the interpolated profile, from bottom to top
792   {
793     double param1 = theMidCurve[i].U - theWidthCurve[i].U / 2;
794     double param2 = theMidCurve[i].U + theWidthCurve[i].U / 2;
795
796     gp_Pnt2d p1 = point.Translated( param1 * profile_dir);
797     gp_Pnt2d p2 = point.Translated( param2 * profile_dir);
798
799     bool arrivedOnTop = false;
800     double z = 0;
801     if (theMidCurve[i].Z <= theMidCurve.MaxZ())
802       z = theMidCurve[i].Z + theMidCurve.DeltaZ();
803     else
804       {
805         z = theMidCurve.MaxZ() + theMidCurve.DeltaZ(); // limit z to linear interpolation between maxima on extremity profiles
806         arrivedOnTop = true; // do not keep points after this one
807       }
808     if (!isOnTop)
809       {
810         AltitudePoint p3d_1( p1.X(), p1.Y(), z ), p3d_2( p2.X(), p2.Y(), z );
811
812         sorted_points[param1] = p3d_1;
813         sorted_points[param2] = p3d_2;
814       }
815     //if (arrivedOnTop)
816     //  isOnTop =true; // do not keep points after this one (commented: leads to strange limits of 2D shape)
817   }
818
819   thePoints.reserve( sorted_points.size() );
820   const double EPS = 1E-12;
821   std::map<double, AltitudePoint>::const_iterator it = sorted_points.begin(), last = sorted_points.end();
822   for( ; it!=last; it++ )
823     if( thePoints.empty() || thePoints.back().SquareDistance( it->second ) > EPS )
824       thePoints.push_back( it->second );
825 }
826
827 //inline double max( double a, double b )
828 //{
829 //  if( a>b )
830 //    return a;
831 //  else
832 //    return b;
833 //}
834 //
835 //inline double min( double a, double b )
836 //{
837 //  if( a<b )
838 //    return a;
839 //  else
840 //    return b;
841 //}
842
843 #include <BRepLib_MakeWire.hxx>
844
845 std::vector<HYDROData_Bathymetry::AltitudePoints> HYDROData_DTM::Interpolate
846   ( const Handle(Geom2d_BSplineCurve)& theHydraulicAxis,
847     const Handle(HYDROData_Profile)& theProfileA,
848     double theXCurvA,
849     const Handle(HYDROData_Profile)& theProfileB,
850     double theXCurvB,
851     double theDDZ, int theNbSteps, bool isAddSecond,
852     int& inter_nb_1, int& inter_nb_2)
853 {
854   double zminA, zmaxA, zminB, zmaxB;
855   gp_Pnt lowestA, lowestB;
856   gp_Vec2d dirA, dirB;
857
858   GetProperties( theProfileA, lowestA, dirA, zminA, zmaxA ); 
859   GetProperties( theProfileB, lowestB, dirB, zminB, zmaxB ); 
860
861   
862   double hmax = zmaxA-zminA > zmaxB-zminB ? zmaxA-zminA : zmaxB-zminB;
863
864   //double dz = zminB - zminA;
865   //double zmin = min( zminA, zminB );
866   //double zmax = max( zmaxA, zmaxB );
867
868   CurveUZ midA(0, gp_Vec2d(), 0, 0), midB(0, gp_Vec2d(), 0, 0);
869   CurveUZ widA(0, gp_Vec2d(), 0, 0), widB(0, gp_Vec2d(), 0, 0);
870
871   ProfileDiscretization( theProfileA, theXCurvA, zminA, zminA+hmax, zmaxA-zminA, theDDZ, midA, widA, inter_nb_1 );
872   ProfileDiscretization( theProfileB, theXCurvB, zminB, zminB+hmax, zmaxB-zminB, theDDZ, midB, widB, inter_nb_2 );
873
874   std::vector<CurveUZ> mid, wid;
875   Interpolate( midA, midB, theNbSteps, mid, isAddSecond );
876   Interpolate( widA, widB, theNbSteps, wid, isAddSecond );
877
878   size_t p = mid.size();
879   size_t q = p>0 ? 2*mid[0].size() : 1;
880   std::vector<AltitudePoints> points;
881   points.resize( p );
882
883   for( size_t i=0; i<p; i++ )
884   {
885     points[i].reserve( q );
886     CurveTo3D( theHydraulicAxis, mid[i], wid[i], points[i] );
887   }
888
889   return points;
890 }
891
892 HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
893   ( const std::vector<Handle(HYDROData_Profile)>& theProfiles,
894     double theDDZ, double theSpatialStep,
895     AltitudePoints& theLeft,
896     AltitudePoints& theRight,
897     std::vector<AltitudePoints>& theMainProfiles,
898     std::set<int>& invalInd)
899 {
900   AltitudePoints points;
901   size_t n = theProfiles.size();
902   if( n<=1 )
903     return points;
904
905   std::vector<double> distances;
906   Handle(Geom2d_BSplineCurve) aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
907   if( aHydraulicAxis.IsNull() )
908     return points;
909
910   theMainProfiles.reserve( n );
911
912   for( size_t i=0, n1=n-1; i<n1; i++ )
913   {
914     double aDistance = distances[i+1]-distances[i];
915     int aNbSteps = int(aDistance/theSpatialStep);
916     bool isAddSecond = i==n1-1;
917
918     // 1. Calculate interpolated profiles
919     int inter_nb_1, inter_nb_2;
920     std::vector<AltitudePoints> local_points = Interpolate( aHydraulicAxis, theProfiles[i], distances[i], 
921       theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond, inter_nb_1, inter_nb_2 );
922     int lps = local_points.size();
923
924     if (inter_nb_1 > 2)
925       invalInd.insert(i);
926
927     if (inter_nb_2 > 2)
928       invalInd.insert(i+1);
929
930     // 2. Put all points into the global container
931     for( size_t j=0; j<lps; j++ )
932     {
933       const AltitudePoints& lp = local_points[j];
934       if( i==0 && j==0 )
935         points.reserve( lp.size() * n );
936       for( size_t k=0, ks=lp.size(); k<ks; k++ )
937         points.push_back( lp[k] );
938     }
939
940     // 3. Get left/right banks' points
941     if( i==0 )
942     {
943       theLeft.reserve( lps * n );
944       theRight.reserve( lps * n );
945     }
946     for( size_t j=0; j<lps; j++ )
947     {
948       const AltitudePoints& lp = local_points[j];
949       theLeft.push_back( lp[0] );
950       theRight.push_back( lp[lp.size()-1] );
951     }
952
953     // 4. Get main profiles points
954     theMainProfiles.push_back( local_points[0] );
955     if( isAddSecond )
956       theMainProfiles.push_back( local_points[lps-1] );
957   }
958   return points;
959 }
960
961 int HYDROData_DTM::EstimateNbPoints( const std::vector<Handle(HYDROData_Profile)>& theProfiles,
962                                      double theDDZ, double theSpatialStep )
963 {
964   size_t n = theProfiles.size();
965   if( n<=1 )
966     return 0;
967   if( theDDZ<1E-6 || theSpatialStep<1E-6 )
968     return 1 << 20;
969
970   std::vector<double> distances;
971   Handle(Geom2d_BSplineCurve) aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
972   if( aHydraulicAxis.IsNull() )
973     return 0;
974
975   double aCompleteDistance = distances[n-1];
976   int aNbSteps = int( aCompleteDistance / theSpatialStep ) + 1;
977   gp_Pnt aLowest;
978   gp_Vec2d aDir;
979   double aZMin, aZMax;
980   GetProperties( theProfiles[0], aLowest, aDir, aZMin, aZMax );
981   int aNbZSteps = (aZMax-aZMin)/theDDZ;
982
983   if( aNbSteps > ( 1<<16 ) || aNbZSteps > ( 1<<16 ) )
984     return 1 << 20;
985
986   return aNbSteps * aNbZSteps;
987 }