]> SALOME platform Git repositories - modules/hydro.git/blob - src/HYDROData/HYDROData_TopoCurve.cxx
Salome HOME
refs #623: Fatal error as result of splitting polyline
[modules/hydro.git] / src / HYDROData / HYDROData_TopoCurve.cxx
1 // Copyright (C) 2007-2015  CEA/DEN, EDF R&D, OPEN CASCADE
2 //
3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
5 //
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
10 //
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14 // Lesser General Public License for more details.
15 //
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
19 //
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 //
22
23 #include <HYDROData_TopoCurve.h>
24
25 #include <Approx_Curve3d.hxx>
26 #include <BRep_Builder.hxx>
27 #include <BRepAdaptor_Curve.hxx>
28 #include <BRepAdaptor_HCurve.hxx>
29 #include <BRepBuilderAPI_MakeEdge.hxx>
30 #include <Extrema_ExtCC.hxx>
31 #include <Extrema_ExtPC.hxx>
32 #include <GeomAPI_Interpolate.hxx>
33 #include <Geom_BSplineCurve.hxx>
34 #include <Precision.hxx>
35 #include <ShapeAnalysis_TransferParametersProj.hxx>
36 #include <ShapeBuild_Edge.hxx>
37 #include <TColgp_Array1OfVec.hxx>
38 #include <TColgp_HArray1OfPnt.hxx>
39 #include <TColStd_HArray1OfBoolean.hxx>
40 #include <TopExp_Explorer.hxx>
41 #include <TopoDS.hxx>
42 #include <TopoDS_Wire.hxx>
43 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
44 #include <TopTools_ListOfShape.hxx>
45
46 //! The type is intended to traverse the container
47 //! either from the begin to the end or vice versa.
48 template<typename ContainerType, typename IteratorType>
49 class Iterator
50 {
51 private:
52   IteratorType myIterator; //!< The iterator.
53   IteratorType myIteratorLimit; //!< The iterator limit.
54
55   //! The pointer to the method to traverse the next item.
56   IteratorType& (IteratorType::*myNext)();
57
58 public:
59   //! The constructor.
60   Iterator(
61     const ContainerType& theContainer,
62     const bool theIsForward)
63   {
64     if (theIsForward)
65     {
66       myIterator = theContainer.begin();
67       myIteratorLimit = theContainer.end();
68       myNext = &IteratorType::operator++;
69     }
70     else
71     {
72       myIterator = --theContainer.end();
73       myIteratorLimit = --theContainer.begin();
74       myNext = &IteratorType::operator--;
75     }
76   }
77
78   //! Returna 'true' if the container contains not yet traversed item.
79   bool More() const
80   {
81     return myIterator != myIteratorLimit;
82   }
83
84   //! Traverses to the next item.
85   IteratorType& operator ++()
86   {
87     return (myIterator.*myNext)();
88   }
89
90   //! Returns the iterator.
91   IteratorType& operator *() {return myIterator;}
92 };
93
94 // Inserts the value after the position.
95 template<typename ItemType> static void InsertAfter(
96   const typename std::list<ItemType>::iterator& thePosition,
97   const ItemType& theValue,
98   std::list<ItemType>& theList)
99 {
100   typename std::list<ItemType>::iterator aEIt2 = thePosition;
101   if (++aEIt2 != theList.end())
102   {
103     theList.insert(aEIt2, theValue);
104   }
105   else
106   {
107     theList.push_back(theValue);
108   }
109 }
110
111 // Converts the curve to a smooth cubic B-spline using the deflection.
112 static Handle(Geom_BSplineCurve) BSpline(
113   const BRepAdaptor_Curve& theCurve, const double theDeflection)
114 {
115   Handle(BRepAdaptor_HCurve) aCurve = new BRepAdaptor_HCurve(theCurve);
116   Approx_Curve3d aConverter(aCurve, theDeflection, GeomAbs_C1, 4, 3);
117   Handle(Geom_BSplineCurve) aBSpline;
118   return aConverter.HasResult() ? aConverter.Curve() : aBSpline;
119 }
120
121 // Replaces the vertex of the edge considering the edge orientation.
122 static TopoDS_Edge ReplaceVertex(
123   const TopoDS_Edge& theEdge, const bool theIsEndVertex)
124 {
125   TopoDS_Vertex aVertices[] = {
126     TopExp::FirstVertex(theEdge, Standard_True),
127     TopExp::LastVertex(theEdge, Standard_True)};
128   aVertices[theIsEndVertex ? 1 : 0].EmptyCopy();
129   TopoDS_Edge aNewEdge = TopoDS::Edge(theEdge.Oriented(TopAbs_FORWARD));
130   aNewEdge =
131     ShapeBuild_Edge().CopyReplaceVertices(aNewEdge, aVertices[0], aVertices[1]);
132   aNewEdge.Orientation(theEdge.Orientation());
133   return aNewEdge;
134 }
135
136 // Projects the point to the curve.
137 static double ProjectPointToCurve(
138   const gp_XYZ& thePoint,
139   const Adaptor3d_Curve& theCurve,
140   double& theParameter)
141 {
142   // Calculate the nearest curve internal extremum.
143   Extrema_ExtPC aAlgo(thePoint, theCurve);
144   int aMinEN = -2;
145   double aMinSqDist = DBL_MAX;
146   if (aAlgo.IsDone())
147   {
148     const int aECount = aAlgo.NbExt();
149     for (int aEN = 1; aEN <= aECount; ++aEN)
150     {
151       const gp_XYZ& aP = aAlgo.Point(aEN).Value().XYZ();
152       const double aSqDist = (thePoint - aP).SquareModulus();
153       if (aMinSqDist > aSqDist)
154       {
155         aMinSqDist = aSqDist;
156         aMinEN = aEN;
157       }
158     }
159   }
160
161   // Calculate the nearest curve end extremum.
162   const double aParams[] =
163     {theCurve.FirstParameter(), theCurve.LastParameter()};
164   const gp_XYZ aEnds[] =
165     {theCurve.Value(aParams[0]).XYZ(), theCurve.Value(aParams[1]).XYZ()};
166   const double aSqDists[] = {
167     (thePoint - aEnds[0]).SquareModulus(),
168     (thePoint - aEnds[1]).SquareModulus()};
169   for (int aEI = 0; aEI < 2; ++aEI)
170   {
171     if (aMinSqDist > aSqDists[aEI])
172     {
173       aMinSqDist = aSqDists[aEI];
174       aMinEN = -aEI;
175     }
176   }
177
178   if (aMinEN <= 0)
179   {
180     theParameter = aParams[-aMinEN];
181     return aMinSqDist;
182   }
183
184   const Extrema_POnCurv& aPOnC = aAlgo.Point(aMinEN);
185   const gp_XYZ& aP = aPOnC.Value().XYZ();
186   theParameter = aPOnC.Parameter();
187   for (int aEI = 0; aEI < 2; ++aEI)
188   {
189     if (Abs(theParameter - aParams[aEI]) < Precision::PConfusion() ||
190       (aP - aEnds[aEI]).SquareModulus() < Precision::SquareConfusion())
191     {
192       theParameter = aParams[aEI];
193     }
194   }
195   return aMinSqDist;
196 }
197
198 // Projects the point to the edge.
199 static double ProjectPointToEdge(
200   const gp_XYZ& thePoint, const TopoDS_Edge& theEdge, double& theParameter)
201 {
202   return ProjectPointToCurve(thePoint, BRepAdaptor_Curve(theEdge), theParameter);
203 }
204
205 // Adds the parameter to the curve parameter list.
206 static int AddParameter(
207   const Adaptor3d_Curve& theCurve,
208   const double theParameter,
209   std::list<double>& theParameters)
210 {
211   // Check the coincidence.
212   std::list<double> aEndParams;
213   aEndParams.push_back(theCurve.FirstParameter());
214   aEndParams.push_back(theCurve.LastParameter());
215   std::list<double>* aParams[] = {&theParameters, &aEndParams};
216   const gp_XYZ aPoint = theCurve.Value(theParameter).XYZ();
217   for (int aLI = 0; aLI < 2; ++aLI)
218   {
219     std::list<double>::iterator aPIt = aParams[aLI]->begin();
220     std::list<double>::iterator aLastPIt = aParams[aLI]->end();
221     for (int aPI = 0; aPIt != aLastPIt; ++aPI, ++aPIt)
222     {
223       const double aParam = *aPIt;
224       if (Abs(theParameter - aParam) < Precision::PConfusion() ||
225         (theCurve.Value(aParam).XYZ() - aPoint).SquareModulus() <=
226           Precision::SquareConfusion())
227       {
228         int aIntCount = 0;
229         if (aLI != 0)
230         {
231           if (aPI == 0)
232           {
233             theParameters.push_front(aEndParams.front());
234           }
235           else
236           {
237             theParameters.push_back(aEndParams.back());
238           }
239           ++aIntCount;
240         }
241         return aIntCount;
242       }
243     }
244   }
245
246   // Calculate the position to insert.
247   std::list<double>::iterator aPIt = theParameters.begin();
248   std::list<double>::iterator aLastPIt = theParameters.end();
249   if (aPIt != aLastPIt && *aPIt < theParameter)
250   {
251     for (++aPIt; aPIt != aLastPIt && *aPIt < theParameter; ++aPIt);
252     if (aPIt != aLastPIt)
253     {
254       theParameters.insert(aPIt, theParameter);
255     }
256     else
257     {
258       theParameters.push_back(theParameter);
259     }
260   }
261   else
262   {
263     theParameters.push_front(theParameter);
264   }
265   return 1;
266 }
267
268 // Intersects the first curve by the second one and
269 // adds the intersection parameters to the ordered list.
270 static int IntersectCurve(
271   const Adaptor3d_Curve& theCurve1,
272   const Adaptor3d_Curve& theCurve2,
273   std::list<double>& theParameters)
274 {
275   // Process the ends.
276   int aIntCount = 0;
277   const gp_XYZ aEndPs[] = {
278     theCurve2.Value(theCurve2.FirstParameter()).XYZ(),
279     theCurve2.Value(theCurve2.LastParameter()).XYZ()};
280   for (int aPI = 0; aPI < 2; ++aPI)
281   {
282     double aParameter;
283     if (ProjectPointToCurve(aEndPs[aPI], theCurve1, aParameter) <=
284       Precision::SquareConfusion())
285     {
286       aIntCount += AddParameter(theCurve1, aParameter, theParameters);
287     }
288   }
289
290   // Process the internal extremums.
291   Extrema_ExtCC aAlgo(theCurve1, theCurve2);
292   if (aAlgo.IsDone())
293   {
294     const int aECount = aAlgo.NbExt();
295     for (int aEN = 1; aEN <= aECount; ++aEN)
296     {
297       Extrema_POnCurv aP1, aP2;
298       aAlgo.Points(aEN, aP1, aP2);
299       if (aP1.Value().SquareDistance(aP2.Value()) <=
300         Precision::SquareConfusion())
301       {
302         aIntCount += AddParameter(theCurve1, aP1.Parameter(), theParameters);
303       }
304     }
305   }
306   return aIntCount;
307 }
308
309 // Intersects the first edge by the second one and
310 // adds the intersection parameters to the ordered list.
311 static int IntersectEdge(
312   const TopoDS_Edge& theEdge1,
313   const TopoDS_Edge& theEdge2,
314   std::list<double>& theParameters)
315 {
316   BRepAdaptor_Curve aCurve1 = BRepAdaptor_Curve(theEdge1);
317   BRepAdaptor_Curve aCurve2 = BRepAdaptor_Curve(theEdge2);
318   return IntersectCurve(aCurve1, aCurve2, theParameters);
319 }
320
321 // Returns the curve tangent in the position: 0 - start, 1 - end.
322 static gp_XYZ Tangent(const Adaptor3d_Curve& theCurve, const int thePosition)
323 {
324   const Standard_Real aParam = (thePosition == 0) ?
325     theCurve.FirstParameter() : theCurve.LastParameter();
326   gp_Pnt aP;
327   gp_Vec aV;
328   theCurve.D1(aParam, aP, aV);
329   Standard_Real aNorm = aV.Magnitude();
330   aNorm = (aNorm >= Precision::PConfusion()) ? aNorm : 0;
331   return ((1 / aNorm) * aV).XYZ();
332 }
333
334 // Returns the edge tangent in the position: 0 - start, 1 - end.
335 static gp_XYZ Tangent(const TopoDS_Edge& theEdge, const int thePosition)
336 {
337   BRepAdaptor_Curve aCurve(theEdge);
338   return Tangent(BRepAdaptor_Curve(theEdge), thePosition);
339 }
340
341 static bool Interpolate(
342   const gp_XYZ thePoint1,
343   const gp_XYZ thePoint2,
344   const gp_XYZ theTangent1,
345   const gp_XYZ theTangent2,
346   Handle(Geom_BSplineCurve)& theBSpline)
347 {
348   Handle(TColgp_HArray1OfPnt) aPs = new TColgp_HArray1OfPnt(1, 2);
349   TColgp_Array1OfVec aTs(1, 2);
350   Handle(TColStd_HArray1OfBoolean) aTFs = new TColStd_HArray1OfBoolean(1, 2);
351   aPs->SetValue(1, thePoint1);
352   aPs->SetValue(2, thePoint2);
353   aTs.SetValue(1, theTangent1);
354   aTs.SetValue(2, theTangent2);
355   aTFs->SetValue(1, Standard_True);
356   aTFs->SetValue(2, Standard_True);
357   GeomAPI_Interpolate aInterpolator(aPs, Standard_False, 0);
358   aInterpolator.Load(aTs, aTFs, Standard_False);
359   aInterpolator.Perform();
360   const bool aResult = (aInterpolator.IsDone() == Standard_True);
361   if (aResult)
362   {
363     theBSpline = aInterpolator.Curve();
364   }
365   return aResult;
366 }
367
368 bool HYDROData_TopoCurve::Initialize(const TopoDS_Wire& theWire)
369 {
370   // Check for nonemptiness.
371   myEdges.clear();
372   TopTools_IndexedDataMapOfShapeListOfShape aVertexToEdges;
373   TopExp::MapShapesAndAncestors(theWire,
374     TopAbs_VERTEX, TopAbs_EDGE, aVertexToEdges);
375   const int aVCount = aVertexToEdges.Extent();
376   if (aVCount == 0)
377   {
378     return false;
379   }
380
381   // Check for 1 manifoldness.
382   bool isClosed = false;
383   {
384     int aEndCount = 0;
385     for (int aVN = 1; aVN <= aVCount; ++aVN)
386     {
387       const int aEdgeCount = aVertexToEdges.FindFromIndex(aVN).Extent();
388       if (aEdgeCount == 1)
389       {
390         ++aEndCount;
391       }
392       if (aEdgeCount > 2)
393       {
394         return false;
395       }
396     }
397     isClosed = (aEndCount == 0);
398     if (!isClosed && aEndCount != 2)
399     {
400       return false;
401     }
402   }
403
404   // Find the start.
405   int aVN = 1;
406   if (!isClosed)
407   {
408     for (; aVertexToEdges.FindFromIndex(aVN).Extent() == 2; ++aVN);
409     if (!aVertexToEdges.FindKey(aVN).IsEqual(TopExp::FirstVertex(
410       TopoDS::Edge(aVertexToEdges.FindFromIndex(aVN).First()), Standard_True)))
411     {
412       for (; aVertexToEdges.FindFromIndex(aVN).Extent() == 2; ++aVN);
413     }
414   }
415   else
416   {
417     TopTools_ListOfShape& aEdges = aVertexToEdges.ChangeFromIndex(1);
418     if (!aVertexToEdges.FindKey(aVN).IsEqual(
419       TopExp::FirstVertex(TopoDS::Edge(aEdges.First()), Standard_True)))
420     {
421       const TopoDS_Shape aEdge = aEdges.First();
422       aEdges.First() = aEdges.Last();
423       aEdges.Last() = aEdge;
424     }
425   }
426
427   // Calculate the edge order.
428   TopTools_ListOfShape* aEdges = &aVertexToEdges.ChangeFromIndex(aVN);
429   while (!aEdges->IsEmpty())
430   {
431     const TopoDS_Edge aEdge = TopoDS::Edge(aEdges->First());
432     aEdges->RemoveFirst();
433     myEdges.push_back(aEdge);
434     int aVN2 = aVertexToEdges.FindIndex(TopExp::FirstVertex(aEdge));
435     if (aVN2 == aVN)
436     {
437       aVN2 = aVertexToEdges.FindIndex(TopExp::LastVertex(aEdge));
438     }
439     aVN = aVN2;
440
441     aEdges = &aVertexToEdges.ChangeFromIndex(aVN2);
442     const TopoDS_Edge aEdge2 = TopoDS::Edge(aEdges->First());
443     if (aEdge2.IsEqual(aEdge))
444     {
445       aEdges->RemoveFirst();
446     }
447     else
448     {
449       aEdges->Clear();
450       aEdges->Append(aEdge2);
451     }
452   }
453
454   // Check for connectedness and free vertex.
455   return aVCount - myEdges.size() == (isClosed ? 0 : 1);
456 }
457
458 TopoDS_Wire HYDROData_TopoCurve::Wire() const
459 {
460   TopoDS_Wire aWire;
461   BRep_Builder aBulder;
462   aBulder.MakeWire(aWire);
463   std::list<TopoDS_Edge>::const_iterator aEItLast = myEdges.end();
464   std::list<TopoDS_Edge>::const_iterator aEIt = myEdges.begin();
465   for (; aEIt != aEItLast; ++aEIt)
466   {
467     aBulder.Add(aWire, *aEIt);
468   }
469   return aWire;
470 }
471
472 void HYDROData_TopoCurve::Cut(
473   const std::list<TopoDS_Edge>::iterator& theEdgePosition,
474   const double theParameter,
475   HYDROData_TopoCurve& theCurve)
476 {
477   // Locate the edge.
478   std::list<TopoDS_Edge>::iterator aFirstEIt = myEdges.begin();
479   std::list<TopoDS_Edge>::iterator aEIt = aFirstEIt;
480   for (; aEIt != theEdgePosition; ++aEIt);
481
482   // Cut the edge.
483   TopoDS_Edge aEdge = *aEIt;
484   BRepAdaptor_Curve aCurve(aEdge);
485   int aParamI = -1;
486   const double aEdgeEndParams[] =
487     {aCurve.FirstParameter(), aCurve.LastParameter()};
488   if (Abs(theParameter - aEdgeEndParams[0]) < Precision::PConfusion())
489   {
490     aParamI = 0;
491   }
492   else if (Abs(theParameter - aEdgeEndParams[1]) < Precision::PConfusion())
493   {
494     aParamI = 1;
495   }
496   const TopAbs_Orientation aOrient = aEdge.Orientation();
497   if (aOrient == TopAbs_REVERSED)
498   {
499     aParamI ^= 1;
500   }
501   if (aParamI < 0)
502   {
503     aEdge.Orientation(TopAbs_FORWARD);
504     TopoDS_Vertex aSplitV1, aSplitV2;
505     BRep_Builder().MakeVertex(
506       aSplitV1, aCurve.Value(theParameter), Precision::Confusion());
507     BRep_Builder().MakeVertex(
508       aSplitV2, aCurve.Value(theParameter), Precision::Confusion());
509     TopoDS_Edge aEParts[] = {
510       ShapeBuild_Edge().CopyReplaceVertices(aEdge, TopoDS_Vertex(),
511         TopoDS::Vertex(aSplitV1.Oriented(TopAbs_REVERSED))),
512       ShapeBuild_Edge().CopyReplaceVertices(aEdge, aSplitV2, TopoDS_Vertex())};
513     ShapeBuild_Edge().CopyPCurves(aEParts[0], aEdge);
514     ShapeBuild_Edge().CopyPCurves(aEParts[1], aEdge);
515     BRep_Builder().SameRange(aEParts[0], Standard_False);
516     BRep_Builder().SameRange(aEParts[1], Standard_False);
517     BRep_Builder().SameParameter(aEParts[0], Standard_False);
518     BRep_Builder().SameParameter(aEParts[1], Standard_False);
519     ShapeAnalysis_TransferParametersProj aSATPP(aEdge, TopoDS_Face());
520     aSATPP.SetMaxTolerance(Precision::Confusion());
521     aSATPP.TransferRange(aEParts[0],
522       aEdgeEndParams[0], theParameter, Standard_False);
523     aSATPP.TransferRange(aEParts[1],
524       theParameter, aEdgeEndParams[1], Standard_False);
525     aEParts[0].Orientation(aOrient);
526     aEParts[1].Orientation(aOrient);
527
528     const int aFirstPI = (aOrient != TopAbs_REVERSED) ? 0 : 1;
529     *aEIt = aEParts[aFirstPI];
530     InsertAfter(aEIt, aEParts[1 - aFirstPI], myEdges);
531     ++aEIt;
532   }
533   else if (aParamI > 0)
534   {
535     ++aEIt;
536   }
537
538   // Calculate the curve parts.
539   std::list<TopoDS_Edge>::iterator aLastEIt = myEdges.end();
540   if (aEIt != aFirstEIt && aEIt != aLastEIt)
541   {
542     std::list<TopoDS_Edge>* aEdges = !IsClosed() ? &theCurve.myEdges : &myEdges;
543     aEdges->splice(aEdges->begin(), myEdges, aEIt, aLastEIt);
544   }
545 }
546
547 void HYDROData_TopoCurve::Cut(
548   const std::list<TopoDS_Edge>::const_iterator& theEdgePosition,
549   const double theParameter,
550   HYDROData_TopoCurve& theCurve1,
551   HYDROData_TopoCurve& theCurve2) const
552 {
553   theCurve1 = *this;
554   std::list<TopoDS_Edge>::const_iterator aEPos = myEdges.begin();
555   std::list<TopoDS_Edge>::iterator aEPos1 = theCurve1.myEdges.begin();
556   for (; aEPos != theEdgePosition; ++aEPos1, ++aEPos);
557   theCurve1.Cut(aEPos1, theParameter, theCurve2);
558 }
559
560 void HYDROData_TopoCurve::Cut(
561   const std::deque<std::list<double> >& theParameters,
562   std::deque<HYDROData_TopoCurve>& theCurves) const
563 {
564   HYDROData_TopoCurve aCurves[2];
565   aCurves[0] = *this;
566   int aCI = 0;
567   std::list<TopoDS_Edge>::iterator aEIt = aCurves[0].myEdges.begin();
568   std::deque<std::list<double> >::const_iterator aPLIt = theParameters.begin();
569   for (std::deque<std::list<double> >::const_iterator aLastPLIt =
570     theParameters.end(); aPLIt != aLastPLIt; ++aPLIt)
571   {
572     const int aSize = aCurves[aCI].myEdges.size();
573     for (Iterator<std::list<double>, std::list<double>::const_iterator> aPIt(
574       *aPLIt, (aEIt->Orientation() != TopAbs_REVERSED)); aPIt.More(); ++aPIt)
575     {
576       const int aCI1 = 1 - aCI;
577       aCurves[aCI].Cut(aEIt, **aPIt, aCurves[aCI1]);
578       if (!aCurves[aCI1].IsEmpty())
579       {
580         theCurves.push_back(HYDROData_TopoCurve());
581         theCurves.back().append(aCurves[aCI]);
582         aEIt = aCurves[aCI1].myEdges.begin();
583         aCI = aCI1;
584       }
585       else
586       {
587         aEIt = aCurves[aCI].myEdges.begin();
588       }
589     }
590     if (aCurves[aCI].myEdges.size() == aSize)
591     {
592       ++aEIt;
593     }
594   }
595   theCurves.push_back(aCurves[aCI]);
596 }
597
598 double HYDROData_TopoCurve::Project(
599   const gp_XYZ& thePoint,
600   std::list<TopoDS_Edge>::const_iterator& theEdgeIterator,
601   double& theParameter) const
602 {
603   double aMinSqDist = DBL_MAX;
604   std::list<TopoDS_Edge>::const_iterator aEIt = myEdges.begin();
605   std::list<TopoDS_Edge>::const_iterator aLastEIt = myEdges.end();
606   for (; aEIt != aLastEIt; ++aEIt)
607   {
608     double aParam;
609     const double aSqDist = ProjectPointToEdge(thePoint, *aEIt, aParam);
610     if (aMinSqDist > aSqDist)
611     {
612       aMinSqDist = aSqDist;
613       theEdgeIterator = aEIt;
614       theParameter = aParam;
615     }
616   }
617   return aMinSqDist;
618 }
619
620 int HYDROData_TopoCurve::Intersect(
621   const TopoDS_Wire& theWire,
622   std::deque<std::list<double> >& theParameters) const
623 {
624   int aIntCount = 0;
625   theParameters.resize(myEdges.size());
626   std::list<TopoDS_Edge>::const_iterator aEIt = myEdges.begin();
627   std::list<TopoDS_Edge>::const_iterator aLastEIt = myEdges.end();
628   std::deque<std::list<double> >::iterator aPIt = theParameters.begin();
629   for (; aEIt != aLastEIt; ++aPIt, ++aEIt)
630   {
631     const TopoDS_Edge& aEdge = *aEIt;
632     std::list<double>& aParams = *aPIt;
633     TopExp_Explorer aEIt2(theWire, TopAbs_EDGE);
634     for (; aEIt2.More(); aEIt2.Next())
635     {
636       aIntCount += IntersectEdge(aEdge,TopoDS::Edge(aEIt2.Current()), aParams);
637     }
638   }
639   return aIntCount;
640 }
641
642 void HYDROData_TopoCurve::CloseCurve()
643 {
644   const TopoDS_Vertex aVertex = TopoDS::Vertex(TopExp::LastVertex(
645     myEdges.back(), Standard_True).Oriented(TopAbs_FORWARD));
646   TopoDS_Edge& aEdge = myEdges.front();
647   const TopoDS_Edge aForwardEdge = TopoDS::Edge(aEdge.Oriented(TopAbs_FORWARD));
648   aEdge = TopoDS::Edge(ShapeBuild_Edge().CopyReplaceVertices(
649     aForwardEdge, aVertex, TopoDS_Vertex()).Oriented(aEdge.Orientation()));
650 }
651
652 void HYDROData_TopoCurve::Merge(
653   const int thePosition, HYDROData_TopoCurve& theCurve)
654 {
655   if (thePosition == 0)
656   {
657     const TopoDS_Vertex aVertex = TopoDS::Vertex(TopExp::LastVertex(
658       theCurve.myEdges.back(), Standard_True).Oriented(TopAbs_FORWARD));
659     TopoDS_Edge& aEdge = myEdges.front();
660     aEdge = TopoDS::Edge(ShapeBuild_Edge().CopyReplaceVertices(
661         TopoDS::Edge(aEdge.Oriented(TopAbs_FORWARD)), aVertex, TopoDS_Vertex()).
662       Oriented(aEdge.Orientation()));
663     prepend(theCurve);
664   }
665   else
666   {
667     const TopoDS_Vertex aVertex = TopoDS::Vertex(TopExp::FirstVertex(
668       theCurve.myEdges.front(), Standard_True).Oriented(TopAbs_REVERSED));
669     TopoDS_Edge& aEdge = myEdges.back();
670     aEdge = TopoDS::Edge(ShapeBuild_Edge().CopyReplaceVertices(
671         TopoDS::Edge(aEdge.Oriented(TopAbs_FORWARD)), TopoDS_Vertex(), aVertex).
672       Oriented(aEdge.Orientation()));
673     append(theCurve);
674   }
675 }
676
677 void HYDROData_TopoCurve::Merge(
678   const double theTolerance, std::deque<HYDROData_TopoCurve>& theCurves)
679 {
680   // Process the curve closeness.
681   const double aSqTol = theTolerance * theTolerance;
682   const gp_XYZ aPs[] = {
683     BRep_Tool::Pnt(TopExp::FirstVertex(myEdges.front(), Standard_True)).XYZ(),
684     BRep_Tool::Pnt(TopExp::LastVertex(myEdges.back(), Standard_True)).XYZ()};
685   bool isClosed = IsClosed();
686   if (!isClosed && (aPs[0] - aPs[1]).SquareModulus() <= aSqTol)
687   {
688     CloseCurve();
689     isClosed = true;
690   }
691
692   // Find the merge places.
693   HYDROData_TopoCurve* aCurves[] = {NULL, NULL};
694   int aOrder = 0;
695   if (!isClosed)
696   {
697     std::deque<HYDROData_TopoCurve>::iterator aCIt = theCurves.begin();
698     std::deque<HYDROData_TopoCurve>::iterator aLastCIt = theCurves.end();
699     for (; aCIt != aLastCIt; ++aCIt)
700     {
701       HYDROData_TopoCurve& aCurve = *aCIt;
702       if (aCurve.IsEmpty() || aCurve.IsClosed())
703       {
704         continue;
705       }
706
707       const gp_XYZ aP1 = BRep_Tool::Pnt(
708         TopExp::FirstVertex(aCurve.myEdges.front(), Standard_True)).XYZ();
709       if (aCurves[0] == NULL && (aPs[1] - aP1).SquareModulus() <= aSqTol)
710       {
711         aCurves[0] = &aCurve;
712       }
713
714       const gp_XYZ aP2 = BRep_Tool::Pnt(
715         TopExp::LastVertex(aCurve.myEdges.back(), Standard_True)).XYZ();
716       if (aCurves[1] == NULL && (aPs[0] - aP2).SquareModulus() <= aSqTol)
717       {
718         aCurves[1] = &aCurve;
719         aOrder = (aCurves[0] == NULL) ? 1 : 0;
720       }
721     }
722   }
723
724   if (aCurves[0] == NULL && aCurves[1] == NULL)
725   {
726     theCurves.push_back(HYDROData_TopoCurve());
727     theCurves.back().append(*this);
728   }
729   else if (aCurves[1] == NULL)
730   {
731     aCurves[0]->Merge(0, *this);
732   }
733   else if (aCurves[0] == NULL)
734   {
735     aCurves[1]->Merge(1, *this);
736   }
737   else
738   {
739     aCurves[aOrder]->Merge(aOrder, *this);
740     if (aCurves[0] != aCurves[1])
741     {
742       aCurves[aOrder]->Merge(aOrder, *aCurves[1 - aOrder]);
743     }
744     else
745     {
746       aCurves[aOrder]->CloseCurve();
747     }
748   }
749 }
750
751 bool HYDROData_TopoCurve::Connect(
752   const Standard_Real theTolerance, std::deque<HYDROData_TopoCurve>& theCurves)
753 {
754   const double aSqTol = theTolerance * theTolerance;
755   const gp_XYZ aPs[] = {
756     BRep_Tool::Pnt(TopExp::FirstVertex(myEdges.front(), Standard_True)).XYZ(),
757     BRep_Tool::Pnt(TopExp::LastVertex(myEdges.back(), Standard_True)).XYZ()};
758   bool isClosed = IsClosed();
759   if (!isClosed && (aPs[0] - aPs[1]).SquareModulus() <= aSqTol)
760   {
761     CloseCurve();
762     isClosed = true;
763   }
764
765   if (!isClosed)
766   {
767     std::deque<HYDROData_TopoCurve>::iterator aCIt = theCurves.begin();
768     std::deque<HYDROData_TopoCurve>::iterator aLastCIt = theCurves.end();
769     for (; aCIt != aLastCIt; ++aCIt)
770     {
771       HYDROData_TopoCurve& aCurve2 = *aCIt;
772       if (aCurve2.IsEmpty() || aCurve2.IsClosed())
773       {
774         continue;
775       }
776
777       const TopoDS_Edge* aEdges2[] =
778         {&aCurve2.myEdges.front(), &aCurve2.myEdges.back()};
779       const gp_XYZ aPs2[] = {
780         BRep_Tool::Pnt(TopExp::FirstVertex(*aEdges2[0], Standard_True)).XYZ(),
781         BRep_Tool::Pnt(TopExp::LastVertex(*aEdges2[1], Standard_True)).XYZ()};
782       const double aSqDists[] =
783         {(aPs[1] - aPs2[0]).SquareModulus(), (aPs[0] - aPs2[1]).SquareModulus()};
784       const int aOrder = (aSqDists[0] <= aSqDists[1]) ? 0 : 1;
785       if (aSqDists[aOrder] > aSqTol)
786       {
787         const TopoDS_Edge& aEdge =
788           (aOrder == 0) ? myEdges.back() : myEdges.front();
789         const gp_XYZ aPs3[] = {aPs[1 - aOrder], aPs2[aOrder]};
790         const gp_XYZ aTs[] =
791           {Tangent(aEdge, 1 - aOrder), Tangent(*aEdges2[aOrder], aOrder)};
792         Handle(Geom_BSplineCurve) aBSpline;
793         if (!Interpolate(aPs3[aOrder], aPs3[1 - aOrder],
794           aTs[aOrder], aTs[1 - aOrder], aBSpline))
795         {
796           return false;
797         }
798
799         HYDROData_TopoCurve aECurve = BRepBuilderAPI_MakeEdge(aBSpline).Edge();
800         aCurve2.Merge(aOrder, aECurve);
801       }
802       aCurve2.Merge(aOrder, *this);
803       if (aSqDists[1 - aOrder] <= aSqTol)
804       {
805         aCurve2.CloseCurve();
806       }
807       return true;
808     }
809   }
810
811   theCurves.push_back(HYDROData_TopoCurve());
812   theCurves.back().append(*this);
813   return true;
814 }
815
816 int HYDROData_TopoCurve::BSplinePiecewiseCurve(
817   const double theDeflection, HYDROData_TopoCurve& theCurve) const
818 {
819   int aPieceCount = 0;
820   std::list<TopoDS_Edge>::const_iterator aLastEIt = myEdges.end();
821   std::list<TopoDS_Edge>::const_iterator aEIt = myEdges.begin();
822   TopoDS_Vertex aEndVertex;
823   TopoDS_Edge aPrevEdge;
824   for (; aEIt != aLastEIt; ++aEIt)
825   {
826     Handle(Geom_BSplineCurve) aBSpline =
827       ::BSpline(BRepAdaptor_Curve(*aEIt), theDeflection);
828     if (aBSpline.IsNull())
829     {
830       return 0;
831     }
832
833     if (aEIt->Orientation() == TopAbs_REVERSED)
834     {
835       aBSpline->Reverse();
836     }
837
838     TopoDS_Edge aEdge;
839     BRep_Builder().MakeEdge(aEdge, aBSpline, Precision::Confusion());
840     TopoDS_Vertex aVertex;
841     BRep_Builder().MakeVertex(aVertex,
842       aBSpline->Value(aBSpline->FirstParameter()), Precision::Confusion());
843     if (!aPrevEdge.IsNull())
844     {
845       BRep_Builder().Add(aPrevEdge, aVertex.Oriented(TopAbs_REVERSED));
846     }
847     else
848     {
849       aEndVertex = aVertex;
850     }
851     BRep_Builder().Add(aEdge, aVertex);
852     theCurve.myEdges.push_back(aEdge);
853     aPieceCount += aBSpline->NbKnots() - 1;
854     aPrevEdge = aEdge;
855   }
856
857   if (!IsClosed())
858   {
859     BRepAdaptor_Curve aCurve(aPrevEdge);
860     BRep_Builder().MakeVertex(aEndVertex,
861       aCurve.Value(aCurve.LastParameter()), Precision::Confusion());
862   }
863   BRep_Builder().Add(aPrevEdge, aEndVertex.Oriented(TopAbs_REVERSED));
864   return aPieceCount;
865 }
866
867 bool HYDROData_TopoCurve::ValuesInKnots(std::list<gp_XYZ>& theValues) const
868 {
869   std::list<TopoDS_Edge>::const_iterator aLastEIt = myEdges.end();
870   std::list<TopoDS_Edge>::const_iterator aEIt = myEdges.begin();
871   for (; aEIt != aLastEIt; ++aEIt)
872   {
873     Handle(Geom_BSplineCurve) aCurve;
874     {
875       TopLoc_Location aLoc;
876       double aParams[2];
877       aCurve = Handle(Geom_BSplineCurve)::
878         DownCast(BRep_Tool::Curve(*aEIt, aLoc, aParams[0], aParams[1]));
879       if (!aLoc.IsIdentity() || aEIt->Orientation() != TopAbs_FORWARD ||
880         aCurve.IsNull())
881       {
882         return false;
883       }
884     }
885
886     for (int aNbKnots = aCurve->NbKnots(), aKN = 1; aKN < aNbKnots; ++aKN)
887     {
888       theValues.push_back(aCurve->Value(aCurve->Knot(aKN)).XYZ());
889     }
890   }
891
892   if (!IsClosed())
893   {
894     TopLoc_Location aLoc;
895     double aParams[2];
896     Handle(Geom_Curve) aCurve =
897       BRep_Tool::Curve(myEdges.back(), aLoc, aParams[0], aParams[1]);
898     theValues.push_back(aCurve->Value(aParams[1]).XYZ());
899   }
900   return true;
901 }