]> SALOME platform Git repositories - modules/hydro.git/blob - src/HYDROData/HYDROData_DTM.cxx
Salome HOME
Merge commit '0d956a814' into pre/IMPS_2016
[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   for( size_t i = 1; i <= n; i++ )
538   {
539     Handle_HYDROData_Profile aProfile = theProfiles[i-1];
540     Handle_HYDROData_Profile aPrevProfile = i==1 ? theProfiles[i-1] : theProfiles[i-2];
541     Handle_HYDROData_Profile aNextProfile = i==n ? theProfiles[i-1] : theProfiles[i];
542
543     if( !IsCooriented( aProfile, aNextProfile ) )
544     {
545       gp_XY lp, rp;
546       aProfile->GetLeftPoint( lp, true );
547       aProfile->GetRightPoint( rp, true );
548       aProfile->SetLeftPoint( rp, true );
549       aProfile->SetRightPoint( lp, true );
550     }
551     aProfile->Update();
552
553     gp_Pnt aLowest;
554     gp_Vec2d aNormal;
555     double zmin, zmax;
556
557     gp_XYZ curP = aProfile->GetBottomPoint();
558     gp_XY curP2d = gp_XY(curP.X(), curP.Y());
559
560     gp_XYZ nextP;
561     if( i==n )
562       nextP = aPrevProfile->GetBottomPoint(true);
563     else
564       nextP = aNextProfile->GetBottomPoint(true);
565
566     gp_XY nextP2d = gp_XY(nextP.X(), nextP.Y());
567
568     gp_Vec2d aPrTangent;
569     GetProperties( aProfile, aLowest, aPrTangent, zmin, zmax );
570     aNormal.SetCoord( -aPrTangent.Y(), aPrTangent.X() );
571
572     gp_Vec2d aDirToNextProfile(nextP2d.X() - curP2d.X(), nextP2d.Y() - curP2d.Y() );
573     if( i==n )
574       aDirToNextProfile.Reverse();
575     if (aNormal.Dot(aDirToNextProfile) < 0)
576       aNormal.Reverse();
577
578     aNormal.Normalize();
579
580     points->SetValue( (int)i, gp_Pnt2d( aLowest.X(), aLowest.Y() ) );
581     tangents.SetValue( (int)i, aNormal );
582     flags->SetValue( (int)i, Standard_True );
583   }
584
585   Geom2dAPI_Interpolate anInterpolator( points, Standard_False, Standard_False );
586   anInterpolator.Load( tangents, flags );
587   anInterpolator.Perform();
588   if( anInterpolator.IsDone() )
589   {
590     aResult = anInterpolator.Curve();
591
592     //fill the distances vector
593     Geom2dAdaptor_Curve anAdaptor( aResult );
594
595     theDistances.clear();
596     theDistances.reserve( n );
597     Standard_Real aParamFirst = anAdaptor.FirstParameter(), aParamLast = anAdaptor.LastParameter();
598     for( size_t i = 1; i <= n; i++ )
599     {
600       gp_Pnt2d aPnt = points->Value( (Standard_Integer)i );
601       Geom2dAPI_ProjectPointOnCurve aProject( aPnt, aResult );
602       Standard_Real aParam = aProject.LowerDistanceParameter();
603       double aDistance = GCPnts_AbscissaPoint::Length( anAdaptor, aParamFirst, aParam );  
604       theDistances.push_back( aDistance );
605     }
606   }
607   return aResult;
608 }
609
610 std::vector<Handle_Geom2d_Curve> HYDROData_DTM::ProfileToParametric( 
611   const Handle_HYDROData_Profile& theProfile,
612   double& theUMin, double& theUMax, gp_Vec2d& theDir )
613 {
614   std::vector<Handle_Geom2d_Curve> curves;
615   
616   // Transformation of the coordinate systems
617   gp_Pnt aLowest;
618   double zmin, zmax;
619   GetProperties( theProfile, aLowest, theDir, zmin, zmax );
620
621   gp_Ax3 aStd3d( gp_Pnt( 0, 0, 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( 1, 0, 0 ) );
622   gp_Ax3 aLocal( gp_Pnt( aLowest.X(), aLowest.Y(), 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( theDir.X(), theDir.Y(), 0 ) );
623
624   gp_Trsf aTransf;
625   aTransf.SetTransformation( aStd3d, aLocal );
626
627   // Iteration via edges
628   TopoDS_Wire aWire = TopoDS::Wire( theProfile->GetShape3D() );
629   TopExp_Explorer anExp( aWire, TopAbs_EDGE );
630   for( ; anExp.More(); anExp.Next() )
631   {
632     // Extract an edge from wire
633     TopoDS_Edge anEdge = TopoDS::Edge( anExp.Current() );
634
635     // Extract a curve corresponding to the edge
636     TopLoc_Location aLoc;
637     Standard_Real aFirst, aLast;
638     Handle(Geom_Curve) aCurve = BRep_Tool::Curve( anEdge, aLoc, aFirst, aLast );
639
640     // Convert the curve to 2d CS
641     Handle(Geom2d_Curve) aCurve2d = CurveTo2D( aCurve, aFirst, aLast, aTransf, theUMin, theUMax );
642     if( !aCurve2d.IsNull() )
643       curves.push_back( aCurve2d );
644   }
645   return curves;
646 }
647
648
649 bool CalcMidWidth( const std::set<double>& intersections, double& theMid, double& theWid )
650 {
651   double umin = std::numeric_limits<double>::max(),
652          umax = -umin;
653
654   size_t n = intersections.size();
655   if( n <= 0 )
656     return false;
657
658   std::set<double>::const_iterator it = intersections.begin(), last = intersections.end();
659   for( ; it!=last; it++ )
660   {
661     double u = *it;
662     if( u<umin )
663       umin = u;
664     if( u>umax )
665       umax = u;
666   }
667   theMid = ( umin+umax )/2;
668   theWid = umax-umin;
669   return true;
670 }
671
672 void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& theProfile, 
673                                            double theXCurv, double theMinZ, double theMaxZ, double theDDZ,
674                                            CurveUZ& theMidPointCurve,
675                                            CurveUZ& theWidthCurve,                                           
676                                            int& intersection_nb,
677                                            double theTolerance)
678 {
679   double aDblMax = std::numeric_limits<double>::max(),
680          aUMin = aDblMax,
681          aUMax = -aUMin,
682          aVMax = 1000000;
683   
684   gp_Vec2d aProfileDir;
685   std::vector<Handle_Geom2d_Curve> curves = ProfileToParametric( theProfile, aUMin, aUMax, aProfileDir );
686   size_t n = curves.size();
687
688   if( n==0 )
689     return;
690
691   // we add the "virtual" vertical lines to simulate the intersection with profile 
692   gp_Pnt2d aFirst, aLast;
693   curves[0]->D0( curves[0]->FirstParameter(), aFirst );
694   curves[n-1]->D0( curves[n-1]->LastParameter(), aLast );
695   Handle(Geom2d_Line) aV1 = new Geom2d_Line( aFirst, gp_Dir2d( 0, 1 ) );
696   Handle(Geom2d_TrimmedCurve) aT1 = new Geom2d_TrimmedCurve( aV1, 0.0, aVMax );
697   
698   Handle(Geom2d_Line) aV2 = new Geom2d_Line( aLast, gp_Dir2d( 0, 1 ) );
699   Handle(Geom2d_TrimmedCurve) aT2 = new Geom2d_TrimmedCurve( aV2, 0.0, aVMax );
700
701   curves.push_back( aT1 );
702   curves.push_back( aT2 );
703   
704   int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 );
705   theMidPointCurve = CurveUZ( theXCurv, aProfileDir, theMinZ );
706   theMidPointCurve.reserve( psize );
707   theWidthCurve = CurveUZ( theXCurv, aProfileDir, theMinZ );
708   theWidthCurve.reserve( psize );
709
710   n = curves.size();
711   // for each discrete value of z we search intersection with profile
712   for( double z1 = theMinZ; z1 <= theMaxZ; z1 += theDDZ )
713   {
714     Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z1 ), gp_Dir2d( 1, 0 ) );
715     std::set<double> intersections;
716     for( size_t i = 0; i < n; i++ )
717     {
718       Handle_Geom2d_Curve aCurve = curves[i];
719       Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
720       for( int k=1, m=anIntersect.NbPoints(); k<=m; k++ )
721         intersections.insert( anIntersect.Point( k ).X() );
722     }
723
724     intersection_nb = intersections.size();
725     if( intersection_nb >= 1 )
726     {
727       double u_mid, u_wid;
728       if( !CalcMidWidth( intersections, u_mid, u_wid ) )
729         continue;
730
731       double z = z1 - theMinZ;
732       PointUZ p_mid;
733       p_mid.U = u_mid;
734       p_mid.Z = z;
735       theMidPointCurve.push_back( p_mid );
736
737       PointUZ p_wid;
738       p_wid.U = u_wid;
739       p_wid.Z = z;
740       theWidthCurve.push_back( p_wid );
741     }
742   }
743 }
744
745 void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCurveB, 
746                                  int theNbSteps, std::vector<CurveUZ>& theInterpolation,
747                                  bool isAddSecond )
748 {
749   theInterpolation.clear();
750   int d = isAddSecond ? 2 : 1;
751   theInterpolation.reserve( theNbSteps+d );
752   double dt = 1.0 / double( theNbSteps + 1 );
753   double t = dt;
754   theInterpolation.push_back( theCurveA );
755   for( int i=0; i<theNbSteps; i++, t+=dt )
756   {
757     CurveUZ anInterp = theCurveA*(1-t) + theCurveB*t;
758     theInterpolation.push_back( anInterp );
759   }
760   if( isAddSecond )
761     theInterpolation.push_back( theCurveB );
762 }
763 #include <BRepLib_MakeEdge2d.hxx>
764 void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
765                                const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
766                                AltitudePoints& thePoints )
767 {
768   Geom2dAdaptor_Curve anAdaptor( theHydraulicAxis );
769   TopoDS_Edge E2d = BRepLib_MakeEdge2d(theHydraulicAxis).Edge();
770   GCPnts_AbscissaPoint ap( anAdaptor, theMidCurve.Xcurv(), anAdaptor.FirstParameter() );  
771   double aParam = ap.Parameter();
772
773   gp_Pnt2d point;
774   anAdaptor.D0( aParam, point );
775   gp_Vec2d profile_dir = theMidCurve.ProfileDir();
776   //gp_Dir tangent_n( -profile_dir.Y(), profile_dir.X(), dz );
777   profile_dir.Normalize();
778   
779   size_t n = theMidCurve.size();
780   std::map<double, AltitudePoint> sorted_points;
781   for( size_t i=0; i<n; i++ )
782   {
783     double param1 = theMidCurve[i].U - theWidthCurve[i].U / 2;
784     double param2 = theMidCurve[i].U + theWidthCurve[i].U / 2;
785
786     gp_Pnt2d p1 = point.Translated( param1 * profile_dir);
787     gp_Pnt2d p2 = point.Translated( param2 * profile_dir);
788
789     double z = theMidCurve[i].Z + theMidCurve.DeltaZ();
790
791     AltitudePoint p3d_1( p1.X(), p1.Y(), z ), p3d_2( p2.X(), p2.Y(), z );
792
793     sorted_points[param1] = p3d_1;
794     sorted_points[param2] = p3d_2;
795   }
796
797   thePoints.reserve( sorted_points.size() );
798   const double EPS = 1E-12;
799   std::map<double, AltitudePoint>::const_iterator it = sorted_points.begin(), last = sorted_points.end();
800   for( ; it!=last; it++ )
801     if( thePoints.empty() || thePoints.back().SquareDistance( it->second ) > EPS )
802       thePoints.push_back( it->second );
803 }
804
805 inline double max( double a, double b )
806 {
807   if( a>b )
808     return a;
809   else
810     return b;
811 }
812
813 inline double min( double a, double b )
814 {
815   if( a<b )
816     return a;
817   else
818     return b;
819 }
820
821 #include <BRepLib_MakeWire.hxx>
822
823 std::vector<HYDROData_Bathymetry::AltitudePoints> HYDROData_DTM::Interpolate
824   ( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
825     const Handle_HYDROData_Profile& theProfileA,
826     double theXCurvA,
827     const Handle_HYDROData_Profile& theProfileB,
828     double theXCurvB,
829     double theDDZ, int theNbSteps, bool isAddSecond,
830     int& inter_nb_1, int& inter_nb_2)
831 {
832   double zminA, zmaxA, zminB, zmaxB;
833   gp_Pnt lowestA, lowestB;
834   gp_Vec2d dirA, dirB;
835
836   GetProperties( theProfileA, lowestA, dirA, zminA, zmaxA ); 
837   GetProperties( theProfileB, lowestB, dirB, zminB, zmaxB ); 
838
839   
840   double hmax = max( zmaxA-zminA, zmaxB-zminB );
841
842   //double dz = zminB - zminA;
843   //double zmin = min( zminA, zminB );
844   //double zmax = max( zmaxA, zmaxB );
845
846   CurveUZ midA(0, gp_Vec2d(), 0), midB(0, gp_Vec2d(), 0);
847   CurveUZ widA(0, gp_Vec2d(), 0), widB(0, gp_Vec2d(), 0);
848
849   ProfileDiscretization( theProfileA, theXCurvA, zminA, zminA+hmax, theDDZ, midA, widA, inter_nb_1 ); 
850   ProfileDiscretization( theProfileB, theXCurvB, zminB, zminB+hmax, theDDZ, midB, widB, inter_nb_2 );
851
852   std::vector<CurveUZ> mid, wid;
853   Interpolate( midA, midB, theNbSteps, mid, isAddSecond );
854   Interpolate( widA, widB, theNbSteps, wid, isAddSecond );
855
856   size_t p = mid.size();
857   size_t q = p>0 ? 2*mid[0].size() : 1;
858   std::vector<AltitudePoints> points;
859   points.resize( p );
860
861   for( size_t i=0; i<p; i++ )
862   {
863     points[i].reserve( q );
864     CurveTo3D( theHydraulicAxis, mid[i], wid[i], points[i] );
865   }
866
867   return points;
868 }
869
870 HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
871   ( const std::vector<Handle_HYDROData_Profile>& theProfiles,
872     double theDDZ, double theSpatialStep,
873     AltitudePoints& theLeft,
874     AltitudePoints& theRight,
875     std::vector<AltitudePoints>& theMainProfiles,
876     std::set<int>& invalInd)
877 {
878   AltitudePoints points;
879   size_t n = theProfiles.size();
880   if( n<=1 )
881     return points;
882
883   std::vector<double> distances;
884   Handle_Geom2d_BSplineCurve aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
885   if( aHydraulicAxis.IsNull() )
886     return points;
887
888   theMainProfiles.reserve( n );
889
890   for( size_t i=0, n1=n-1; i<n1; i++ )
891   {
892     double aDistance = distances[i+1]-distances[i];
893     int aNbSteps = int(aDistance/theSpatialStep);
894     bool isAddSecond = i==n1-1;
895
896     // 1. Calculate interpolated profiles
897     int inter_nb_1, inter_nb_2;
898     std::vector<AltitudePoints> local_points = Interpolate( aHydraulicAxis, theProfiles[i], distances[i], 
899       theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond, inter_nb_1, inter_nb_2 );
900     int lps = local_points.size();
901
902     if (inter_nb_1 > 2)
903       invalInd.insert(i);
904
905     if (inter_nb_2 > 2)
906       invalInd.insert(i+1);
907
908     // 2. Put all points into the global container
909     for( size_t j=0; j<lps; j++ )
910     {
911       const AltitudePoints& lp = local_points[j];
912       if( i==0 && j==0 )
913         points.reserve( lp.size() * n );
914       for( size_t k=0, ks=lp.size(); k<ks; k++ )
915         points.push_back( lp[k] );
916     }
917
918     // 3. Get left/right banks' points
919     if( i==0 )
920     {
921       theLeft.reserve( lps * n );
922       theRight.reserve( lps * n );
923     }
924     for( size_t j=0; j<lps; j++ )
925     {
926       const AltitudePoints& lp = local_points[j];
927       theLeft.push_back( lp[0] );
928       theRight.push_back( lp[lp.size()-1] );
929     }
930
931     // 4. Get main profiles points
932     theMainProfiles.push_back( local_points[0] );
933     if( isAddSecond )
934       theMainProfiles.push_back( local_points[lps-1] );
935   }
936   return points;
937 }
938
939 int HYDROData_DTM::EstimateNbPoints( const std::vector<Handle_HYDROData_Profile>& theProfiles,
940                                      double theDDZ, double theSpatialStep )
941 {
942   size_t n = theProfiles.size();
943   if( n<=1 )
944     return 0;
945   if( theDDZ<1E-6 || theSpatialStep<1E-6 )
946     return 1 << 20;
947
948   std::vector<double> distances;
949   Handle_Geom2d_BSplineCurve aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
950   if( aHydraulicAxis.IsNull() )
951     return 0;
952
953   double aCompleteDistance = distances[n-1];
954   int aNbSteps = int( aCompleteDistance / theSpatialStep ) + 1;
955   gp_Pnt aLowest;
956   gp_Vec2d aDir;
957   double aZMin, aZMax;
958   GetProperties( theProfiles[0], aLowest, aDir, aZMin, aZMax );
959   int aNbZSteps = (aZMax-aZMin)/theDDZ;
960
961   if( aNbSteps > ( 1<<16 ) || aNbZSteps > ( 1<<16 ) )
962     return 1 << 20;
963
964   return aNbSteps * aNbZSteps;
965 }