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