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