]> SALOME platform Git repositories - modules/hydro.git/blob - src/HYDROData/HYDROData_TopoCurve.cxx
Salome HOME
Merge branch 'BR_v14_rc' into BR_quadtree
[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   const bool isClosed = IsClosed();
502   if (aParamI < 0)
503   {
504     aEdge.Orientation(TopAbs_FORWARD);
505     TopoDS_Vertex aSplitV1, aSplitV2;
506     BRep_Builder().MakeVertex(
507       aSplitV1, aCurve.Value(theParameter), Precision::Confusion());
508     BRep_Builder().MakeVertex(
509       aSplitV2, aCurve.Value(theParameter), Precision::Confusion());
510     TopoDS_Edge aEParts[] = {
511       ShapeBuild_Edge().CopyReplaceVertices(aEdge, TopoDS_Vertex(),
512         TopoDS::Vertex(aSplitV1.Oriented(TopAbs_REVERSED))),
513       ShapeBuild_Edge().CopyReplaceVertices(aEdge, aSplitV2, TopoDS_Vertex())};
514     ShapeBuild_Edge().CopyPCurves(aEParts[0], aEdge);
515     ShapeBuild_Edge().CopyPCurves(aEParts[1], aEdge);
516     BRep_Builder().SameRange(aEParts[0], Standard_False);
517     BRep_Builder().SameRange(aEParts[1], Standard_False);
518     BRep_Builder().SameParameter(aEParts[0], Standard_False);
519     BRep_Builder().SameParameter(aEParts[1], Standard_False);
520     ShapeAnalysis_TransferParametersProj aSATPP(aEdge, TopoDS_Face());
521     aSATPP.SetMaxTolerance(Precision::Confusion());
522     aSATPP.TransferRange(aEParts[0],
523       aEdgeEndParams[0], theParameter, Standard_False);
524     aSATPP.TransferRange(aEParts[1],
525       theParameter, aEdgeEndParams[1], Standard_False);
526     aEParts[0].Orientation(aOrient);
527     aEParts[1].Orientation(aOrient);
528
529     const int aFirstPI = (aOrient != TopAbs_REVERSED) ? 0 : 1;
530     *aEIt = aEParts[aFirstPI];
531     InsertAfter(aEIt, aEParts[1 - aFirstPI], myEdges);
532     ++aEIt;
533   }
534   else
535   {
536     TopoDS_Edge aNewEdge = ReplaceVertex(aEdge, (aParamI == 0) ? false : true);
537     *aEIt = aNewEdge;
538     if (aParamI > 0)
539     {
540       ++aEIt;
541     }
542   }
543
544   // Calculate the curve parts.
545   std::list<TopoDS_Edge>::iterator aLastEIt = myEdges.end();
546   if (aEIt != aFirstEIt && aEIt != aLastEIt)
547   {
548     std::list<TopoDS_Edge>* aEdges = !isClosed ? &theCurve.myEdges : &myEdges;
549     aEdges->splice(aEdges->begin(), myEdges, aEIt, aLastEIt);
550   }
551 }
552
553 void HYDROData_TopoCurve::Cut(
554   const std::list<TopoDS_Edge>::const_iterator& theEdgePosition,
555   const double theParameter,
556   HYDROData_TopoCurve& theCurve1,
557   HYDROData_TopoCurve& theCurve2) const
558 {
559   theCurve1 = *this;
560   std::list<TopoDS_Edge>::const_iterator aEPos = myEdges.begin();
561   std::list<TopoDS_Edge>::iterator aEPos1 = theCurve1.myEdges.begin();
562   for (; aEPos != theEdgePosition; ++aEPos1, ++aEPos);
563   theCurve1.Cut(aEPos1, theParameter, theCurve2);
564 }
565
566 void HYDROData_TopoCurve::Cut(
567   const std::deque<std::list<double> >& theParameters,
568   std::deque<HYDROData_TopoCurve>& theCurves) const
569 {
570   HYDROData_TopoCurve aCurves[2];
571   aCurves[0] = *this;
572   int aCI = 0;
573   std::list<TopoDS_Edge>::iterator aEIt = aCurves[0].myEdges.begin();
574   std::deque<std::list<double> >::const_iterator aPLIt = theParameters.begin();
575   for (std::deque<std::list<double> >::const_iterator aLastPLIt =
576     theParameters.end(); aPLIt != aLastPLIt; ++aPLIt)
577   {
578     TopoDS_Edge aNextEdge;
579     {
580       std::list<TopoDS_Edge>::iterator aNextEIt = aEIt;
581       ++aNextEIt;
582       if (aNextEIt != aCurves[aCI].myEdges.end())
583       {
584         aNextEdge = *aNextEIt;
585       }
586     }
587
588     for (Iterator<std::list<double>, std::list<double>::const_iterator> aPIt(
589       *aPLIt, (aEIt->Orientation() != TopAbs_REVERSED)); aPIt.More(); ++aPIt)
590     {
591       const int aCI1 = 1 - aCI;
592       aCurves[aCI].Cut(aEIt, **aPIt, aCurves[aCI1]);
593       if (!aCurves[aCI1].IsEmpty())
594       {
595         theCurves.push_back(HYDROData_TopoCurve());
596         theCurves.back().append(aCurves[aCI]);
597         aEIt = aCurves[aCI1].myEdges.begin();
598         aCI = aCI1;
599       }
600       else
601       {
602         aEIt = aCurves[aCI].myEdges.begin();
603       }
604     }
605
606     if (!aNextEdge.IsNull() && !aEIt->IsEqual(aNextEdge))
607     {
608       ++aEIt;
609     }
610   }
611   theCurves.push_back(aCurves[aCI]);
612 }
613
614 double HYDROData_TopoCurve::Project(
615   const gp_XYZ& thePoint,
616   std::list<TopoDS_Edge>::const_iterator& theEdgeIterator,
617   double& theParameter) const
618 {
619   double aMinSqDist = DBL_MAX;
620   std::list<TopoDS_Edge>::const_iterator aEIt = myEdges.begin();
621   std::list<TopoDS_Edge>::const_iterator aLastEIt = myEdges.end();
622   for (; aEIt != aLastEIt; ++aEIt)
623   {
624     double aParam;
625     const double aSqDist = ProjectPointToEdge(thePoint, *aEIt, aParam);
626     if (aMinSqDist > aSqDist)
627     {
628       aMinSqDist = aSqDist;
629       theEdgeIterator = aEIt;
630       theParameter = aParam;
631     }
632   }
633   return aMinSqDist;
634 }
635
636 int HYDROData_TopoCurve::Intersect(
637   const TopoDS_Wire& theWire,
638   std::deque<std::list<double> >& theParameters) const
639 {
640   int aIntCount = 0;
641   theParameters.resize(myEdges.size());
642   std::list<TopoDS_Edge>::const_iterator aEIt = myEdges.begin();
643   std::list<TopoDS_Edge>::const_iterator aLastEIt = myEdges.end();
644   std::deque<std::list<double> >::iterator aPIt = theParameters.begin();
645   for (; aEIt != aLastEIt; ++aPIt, ++aEIt)
646   {
647     const TopoDS_Edge& aEdge = *aEIt;
648     std::list<double>& aParams = *aPIt;
649     TopExp_Explorer aEIt2(theWire, TopAbs_EDGE);
650     for (; aEIt2.More(); aEIt2.Next())
651     {
652       aIntCount += IntersectEdge(aEdge,TopoDS::Edge(aEIt2.Current()), aParams);
653     }
654   }
655   return aIntCount;
656 }
657
658 void HYDROData_TopoCurve::CloseCurve()
659 {
660   const TopoDS_Vertex aVertex = TopoDS::Vertex(TopExp::LastVertex(
661     myEdges.back(), Standard_True).Oriented(TopAbs_FORWARD));
662   TopoDS_Edge& aEdge = myEdges.front();
663   const TopoDS_Edge aForwardEdge = TopoDS::Edge(aEdge.Oriented(TopAbs_FORWARD));
664   aEdge = TopoDS::Edge(ShapeBuild_Edge().CopyReplaceVertices(
665     aForwardEdge, aVertex, TopoDS_Vertex()).Oriented(aEdge.Orientation()));
666 }
667
668 void HYDROData_TopoCurve::Merge(
669   const int thePosition, HYDROData_TopoCurve& theCurve)
670 {
671   if (thePosition == 0)
672   {
673     const TopoDS_Vertex aVertex = TopoDS::Vertex(TopExp::LastVertex(
674       theCurve.myEdges.back(), Standard_True).Oriented(TopAbs_FORWARD));
675     TopoDS_Edge& aEdge = myEdges.front();
676     aEdge = TopoDS::Edge(ShapeBuild_Edge().CopyReplaceVertices(
677         TopoDS::Edge(aEdge.Oriented(TopAbs_FORWARD)), aVertex, TopoDS_Vertex()).
678       Oriented(aEdge.Orientation()));
679     prepend(theCurve);
680   }
681   else
682   {
683     const TopoDS_Vertex aVertex = TopoDS::Vertex(TopExp::FirstVertex(
684       theCurve.myEdges.front(), Standard_True).Oriented(TopAbs_REVERSED));
685     TopoDS_Edge& aEdge = myEdges.back();
686     aEdge = TopoDS::Edge(ShapeBuild_Edge().CopyReplaceVertices(
687         TopoDS::Edge(aEdge.Oriented(TopAbs_FORWARD)), TopoDS_Vertex(), aVertex).
688       Oriented(aEdge.Orientation()));
689     append(theCurve);
690   }
691 }
692
693 void HYDROData_TopoCurve::Merge(
694   const double theTolerance, std::deque<HYDROData_TopoCurve>& theCurves)
695 {
696   // Process the curve closeness.
697   const double aSqTol = theTolerance * theTolerance;
698   const gp_XYZ aPs[] = {
699     BRep_Tool::Pnt(TopExp::FirstVertex(myEdges.front(), Standard_True)).XYZ(),
700     BRep_Tool::Pnt(TopExp::LastVertex(myEdges.back(), Standard_True)).XYZ()};
701   bool isClosed = IsClosed();
702   if (!isClosed && (aPs[0] - aPs[1]).SquareModulus() <= aSqTol)
703   {
704     CloseCurve();
705     isClosed = true;
706   }
707
708   // Find the merge places.
709   HYDROData_TopoCurve* aCurves[] = {NULL, NULL};
710   int aOrder = 0;
711   if (!isClosed)
712   {
713     std::deque<HYDROData_TopoCurve>::iterator aCIt = theCurves.begin();
714     std::deque<HYDROData_TopoCurve>::iterator aLastCIt = theCurves.end();
715     for (; aCIt != aLastCIt; ++aCIt)
716     {
717       HYDROData_TopoCurve& aCurve = *aCIt;
718       if (aCurve.IsEmpty() || aCurve.IsClosed())
719       {
720         continue;
721       }
722
723       const gp_XYZ aP1 = BRep_Tool::Pnt(
724         TopExp::FirstVertex(aCurve.myEdges.front(), Standard_True)).XYZ();
725       if (aCurves[0] == NULL && (aPs[1] - aP1).SquareModulus() <= aSqTol)
726       {
727         aCurves[0] = &aCurve;
728       }
729
730       const gp_XYZ aP2 = BRep_Tool::Pnt(
731         TopExp::LastVertex(aCurve.myEdges.back(), Standard_True)).XYZ();
732       if (aCurves[1] == NULL && (aPs[0] - aP2).SquareModulus() <= aSqTol)
733       {
734         aCurves[1] = &aCurve;
735         aOrder = (aCurves[0] == NULL) ? 1 : 0;
736       }
737     }
738   }
739
740   if (aCurves[0] == NULL && aCurves[1] == NULL)
741   {
742     theCurves.push_back(HYDROData_TopoCurve());
743     theCurves.back().append(*this);
744   }
745   else if (aCurves[1] == NULL)
746   {
747     aCurves[0]->Merge(0, *this);
748   }
749   else if (aCurves[0] == NULL)
750   {
751     aCurves[1]->Merge(1, *this);
752   }
753   else
754   {
755     aCurves[aOrder]->Merge(aOrder, *this);
756     if (aCurves[0] != aCurves[1])
757     {
758       aCurves[aOrder]->Merge(aOrder, *aCurves[1 - aOrder]);
759     }
760     else
761     {
762       aCurves[aOrder]->CloseCurve();
763     }
764   }
765 }
766
767 bool HYDROData_TopoCurve::Connect(
768   const Standard_Real theTolerance, std::deque<HYDROData_TopoCurve>& theCurves)
769 {
770   const double aSqTol = theTolerance * theTolerance;
771   const gp_XYZ aPs[] = {
772     BRep_Tool::Pnt(TopExp::FirstVertex(myEdges.front(), Standard_True)).XYZ(),
773     BRep_Tool::Pnt(TopExp::LastVertex(myEdges.back(), Standard_True)).XYZ()};
774   bool isClosed = IsClosed();
775   if (!isClosed && (aPs[0] - aPs[1]).SquareModulus() <= aSqTol)
776   {
777     CloseCurve();
778     isClosed = true;
779   }
780
781   if (!isClosed)
782   {
783     std::deque<HYDROData_TopoCurve>::iterator aCIt = theCurves.begin();
784     std::deque<HYDROData_TopoCurve>::iterator aLastCIt = theCurves.end();
785     for (; aCIt != aLastCIt; ++aCIt)
786     {
787       HYDROData_TopoCurve& aCurve2 = *aCIt;
788       if (aCurve2.IsEmpty() || aCurve2.IsClosed())
789       {
790         continue;
791       }
792
793       const TopoDS_Edge* aEdges2[] =
794         {&aCurve2.myEdges.front(), &aCurve2.myEdges.back()};
795       const gp_XYZ aPs2[] = {
796         BRep_Tool::Pnt(TopExp::FirstVertex(*aEdges2[0], Standard_True)).XYZ(),
797         BRep_Tool::Pnt(TopExp::LastVertex(*aEdges2[1], Standard_True)).XYZ()};
798       const double aSqDists[] =
799         {(aPs[1] - aPs2[0]).SquareModulus(), (aPs[0] - aPs2[1]).SquareModulus()};
800       const int aOrder = (aSqDists[0] <= aSqDists[1]) ? 0 : 1;
801       if (aSqDists[aOrder] > aSqTol)
802       {
803         const TopoDS_Edge& aEdge =
804           (aOrder == 0) ? myEdges.back() : myEdges.front();
805         const gp_XYZ aPs3[] = {aPs[1 - aOrder], aPs2[aOrder]};
806         const gp_XYZ aTs[] =
807           {Tangent(aEdge, 1 - aOrder), Tangent(*aEdges2[aOrder], aOrder)};
808         Handle(Geom_BSplineCurve) aBSpline;
809         if (!Interpolate(aPs3[aOrder], aPs3[1 - aOrder],
810           aTs[aOrder], aTs[1 - aOrder], aBSpline))
811         {
812           return false;
813         }
814
815         HYDROData_TopoCurve aECurve = BRepBuilderAPI_MakeEdge(aBSpline).Edge();
816         aCurve2.Merge(aOrder, aECurve);
817       }
818       aCurve2.Merge(aOrder, *this);
819       if (aSqDists[1 - aOrder] <= aSqTol)
820       {
821         aCurve2.CloseCurve();
822       }
823       return true;
824     }
825   }
826
827   theCurves.push_back(HYDROData_TopoCurve());
828   theCurves.back().append(*this);
829   return true;
830 }
831
832 int HYDROData_TopoCurve::BSplinePiecewiseCurve(
833   const double theDeflection, HYDROData_TopoCurve& theCurve) const
834 {
835   int aPieceCount = 0;
836   std::list<TopoDS_Edge>::const_iterator aLastEIt = myEdges.end();
837   std::list<TopoDS_Edge>::const_iterator aEIt = myEdges.begin();
838   TopoDS_Vertex aEndVertex;
839   TopoDS_Edge aPrevEdge;
840   for (; aEIt != aLastEIt; ++aEIt)
841   {
842     Handle(Geom_BSplineCurve) aBSpline =
843       ::BSpline(BRepAdaptor_Curve(*aEIt), theDeflection);
844     if (aBSpline.IsNull())
845     {
846       return 0;
847     }
848
849     if (aEIt->Orientation() == TopAbs_REVERSED)
850     {
851       aBSpline->Reverse();
852     }
853
854     TopoDS_Edge aEdge;
855     BRep_Builder().MakeEdge(aEdge, aBSpline, Precision::Confusion());
856     TopoDS_Vertex aVertex;
857     BRep_Builder().MakeVertex(aVertex,
858       aBSpline->Value(aBSpline->FirstParameter()), Precision::Confusion());
859     if (!aPrevEdge.IsNull())
860     {
861       BRep_Builder().Add(aPrevEdge, aVertex.Oriented(TopAbs_REVERSED));
862     }
863     else
864     {
865       aEndVertex = aVertex;
866     }
867     BRep_Builder().Add(aEdge, aVertex);
868     theCurve.myEdges.push_back(aEdge);
869     aPieceCount += aBSpline->NbKnots() - 1;
870     aPrevEdge = aEdge;
871   }
872
873   if (!IsClosed())
874   {
875     BRepAdaptor_Curve aCurve(aPrevEdge);
876     BRep_Builder().MakeVertex(aEndVertex,
877       aCurve.Value(aCurve.LastParameter()), Precision::Confusion());
878   }
879   BRep_Builder().Add(aPrevEdge, aEndVertex.Oriented(TopAbs_REVERSED));
880   return aPieceCount;
881 }
882
883 bool HYDROData_TopoCurve::ValuesInKnots(std::list<gp_XYZ>& theValues) const
884 {
885   std::list<TopoDS_Edge>::const_iterator aLastEIt = myEdges.end();
886   std::list<TopoDS_Edge>::const_iterator aEIt = myEdges.begin();
887   for (; aEIt != aLastEIt; ++aEIt)
888   {
889     Handle(Geom_BSplineCurve) aCurve;
890     {
891       TopLoc_Location aLoc;
892       double aParams[2];
893       aCurve = Handle(Geom_BSplineCurve)::
894         DownCast(BRep_Tool::Curve(*aEIt, aLoc, aParams[0], aParams[1]));
895       if (!aLoc.IsIdentity() || aEIt->Orientation() != TopAbs_FORWARD ||
896         aCurve.IsNull())
897       {
898         return false;
899       }
900     }
901
902     for (int aNbKnots = aCurve->NbKnots(), aKN = 1; aKN < aNbKnots; ++aKN)
903     {
904       theValues.push_back(aCurve->Value(aCurve->Knot(aKN)).XYZ());
905     }
906   }
907
908   if (!IsClosed())
909   {
910     TopLoc_Location aLoc;
911     double aParams[2];
912     Handle(Geom_Curve) aCurve =
913       BRep_Tool::Curve(myEdges.back(), aLoc, aParams[0], aParams[1]);
914     theValues.push_back(aCurve->Value(aParams[1]).XYZ());
915   }
916   return true;
917 }