]> SALOME platform Git repositories - modules/hydro.git/blob - src/HYDROData/HYDROData_TopoCurve.cxx
Salome HOME
Merge branch 'BR_H2018_2' of https://codev-tuleap.cea.fr/plugins/git/salome/hydro...
[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 //#define _DEVDEBUG_
43 #include "HYDRO_trace.hxx"
44 #include <BRepTools.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 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       DEBTRACE("aParameter " << aParameter);
287       aIntCount += AddParameter(theCurve1, aParameter, theParameters);
288     }
289   }
290
291   // Process the internal extremums.
292   Extrema_ExtCC aAlgo(theCurve1, theCurve2); //, 1.e-6, 1.e-6);
293   if (aAlgo.IsDone())
294   {
295     const int aECount = aAlgo.NbExt();
296     for (int aEN = 1; aEN <= aECount; ++aEN)
297     {
298       Extrema_POnCurv aP1, aP2;
299       aAlgo.Points(aEN, aP1, aP2);
300       DEBTRACE("SquareDistance " << aP1.Value().SquareDistance(aP2.Value()));
301       if (aP1.Value().SquareDistance(aP2.Value()) <=
302         Precision::SquareConfusion())
303       {
304         DEBTRACE("aP1.Parameter() " << aP1.Parameter());
305         aIntCount += AddParameter(theCurve1, aP1.Parameter(), theParameters);
306       }
307     }
308   }
309   return aIntCount;
310 }
311
312 // Intersects the first edge by the second one and
313 // adds the intersection parameters to the ordered list.
314 static int IntersectEdge(
315   const TopoDS_Edge& theEdge1,
316   const TopoDS_Edge& theEdge2,
317   std::list<double>& theParameters)
318 {
319   BRepAdaptor_Curve aCurve1 = BRepAdaptor_Curve(theEdge1);
320   BRepAdaptor_Curve aCurve2 = BRepAdaptor_Curve(theEdge2);
321   return IntersectCurve(aCurve1, aCurve2, theParameters);
322 }
323
324 // Returns the curve tangent in the position: 0 - start, 1 - end.
325 static gp_XYZ Tangent(const Adaptor3d_Curve& theCurve, const int thePosition)
326 {
327   const Standard_Real aParam = (thePosition == 0) ?
328     theCurve.FirstParameter() : theCurve.LastParameter();
329   gp_Pnt aP;
330   gp_Vec aV;
331   theCurve.D1(aParam, aP, aV);
332   Standard_Real aNorm = aV.Magnitude();
333   aNorm = (aNorm >= Precision::PConfusion()) ? aNorm : 0;
334   return ((1 / aNorm) * aV).XYZ();
335 }
336
337 // Returns the edge tangent in the position: 0 - start, 1 - end.
338 static gp_XYZ Tangent(const TopoDS_Edge& theEdge, const int thePosition)
339 {
340   BRepAdaptor_Curve aCurve(theEdge);
341   return Tangent(BRepAdaptor_Curve(theEdge), thePosition);
342 }
343
344 static bool Interpolate(
345   const gp_XYZ thePoint1,
346   const gp_XYZ thePoint2,
347   const gp_XYZ theTangent1,
348   const gp_XYZ theTangent2,
349   Handle(Geom_BSplineCurve)& theBSpline)
350 {
351   Handle(TColgp_HArray1OfPnt) aPs = new TColgp_HArray1OfPnt(1, 2);
352   TColgp_Array1OfVec aTs(1, 2);
353   Handle(TColStd_HArray1OfBoolean) aTFs = new TColStd_HArray1OfBoolean(1, 2);
354   aPs->SetValue(1, thePoint1);
355   aPs->SetValue(2, thePoint2);
356   aTs.SetValue(1, theTangent1);
357   aTs.SetValue(2, theTangent2);
358   aTFs->SetValue(1, Standard_True);
359   aTFs->SetValue(2, Standard_True);
360   GeomAPI_Interpolate aInterpolator(aPs, Standard_False, 0);
361   aInterpolator.Load(aTs, aTFs, Standard_False);
362   aInterpolator.Perform();
363   const bool aResult = (aInterpolator.IsDone() == Standard_True);
364   if (aResult)
365   {
366     theBSpline = aInterpolator.Curve();
367   }
368   return aResult;
369 }
370
371 bool HYDROData_TopoCurve::Initialize(const TopoDS_Wire& theWire)
372 {
373   // Check for non-emptiness.
374   myEdges.clear();
375   TopTools_IndexedDataMapOfShapeListOfShape aVertexToEdges;
376   TopExp::MapShapesAndAncestors(theWire,
377     TopAbs_VERTEX, TopAbs_EDGE, aVertexToEdges);
378   const int aVCount = aVertexToEdges.Extent();
379   DEBTRACE("initialize VCount= "<< aVCount);
380   if (aVCount == 0)
381   {
382     return false;
383   }
384
385   // Check for 1 manifoldness.
386   bool isClosed = false;
387   {
388     int aEndCount = 0;
389     for (int aVN = 1; aVN <= aVCount; ++aVN)
390     {
391       const int aEdgeCount = aVertexToEdges.FindFromIndex(aVN).Extent();
392       if (aEdgeCount == 1)
393       {
394         ++aEndCount;
395       }
396       if (aEdgeCount > 2)
397       {
398         return false;
399       }
400     }
401     isClosed = (aEndCount == 0);
402     if (!isClosed && aEndCount != 2)
403     {
404       return false;
405     }
406   }
407
408   // Find the start.
409   int aVN = 1;
410   if (!isClosed)
411   {
412     for (; aVertexToEdges.FindFromIndex(aVN).Extent() == 2; ++aVN);
413     if (!aVertexToEdges.FindKey(aVN).IsEqual(TopExp::FirstVertex(
414       TopoDS::Edge(aVertexToEdges.FindFromIndex(aVN).First()), Standard_True)))
415     {
416       for (; aVertexToEdges.FindFromIndex(aVN).Extent() == 2; ++aVN);
417     }
418   }
419   else
420   {
421     TopTools_ListOfShape& aEdges = aVertexToEdges.ChangeFromIndex(1);
422     if (!aVertexToEdges.FindKey(aVN).IsEqual(
423       TopExp::FirstVertex(TopoDS::Edge(aEdges.First()), Standard_True)))
424     {
425       const TopoDS_Shape aEdge = aEdges.First();
426       aEdges.First() = aEdges.Last();
427       aEdges.Last() = aEdge;
428     }
429   }
430
431   // Calculate the edge order.
432   TopTools_ListOfShape* aEdges = &aVertexToEdges.ChangeFromIndex(aVN);
433   while (!aEdges->IsEmpty())
434   {
435     const TopoDS_Edge aEdge = TopoDS::Edge(aEdges->First());
436     aEdges->RemoveFirst();
437     myEdges.push_back(aEdge);
438     int aVN2 = aVertexToEdges.FindIndex(TopExp::FirstVertex(aEdge));
439     if (aVN2 == aVN)
440     {
441       aVN2 = aVertexToEdges.FindIndex(TopExp::LastVertex(aEdge));
442     }
443     aVN = aVN2;
444
445     aEdges = &aVertexToEdges.ChangeFromIndex(aVN2);
446     const TopoDS_Edge aEdge2 = TopoDS::Edge(aEdges->First());
447     if (aEdge2.IsEqual(aEdge))
448     {
449       aEdges->RemoveFirst();
450     }
451     else
452     {
453       aEdges->Clear();
454       aEdges->Append(aEdge2);
455     }
456   }
457
458   // Check for connectedness and free vertex.
459   return aVCount - myEdges.size() == (isClosed ? 0 : 1);
460 }
461
462 TopoDS_Wire HYDROData_TopoCurve::Wire() const
463 {
464   TopoDS_Wire aWire;
465   BRep_Builder aBuilder;
466   aBuilder.MakeWire(aWire);
467   std::list<TopoDS_Edge>::const_iterator aEItLast = myEdges.end();
468   std::list<TopoDS_Edge>::const_iterator aEIt = myEdges.begin();
469   for (; aEIt != aEItLast; ++aEIt)
470   {
471     aBuilder.Add(aWire, *aEIt);
472   }
473   return aWire;
474 }
475
476 bool HYDROData_TopoCurve::Cut(
477   const std::list<TopoDS_Edge>::iterator& theEdgePosition,
478   const double theParameter,
479   HYDROData_TopoCurve& theCurve)
480 {
481   bool aResult = false;
482
483   // Locate the edge.
484   std::list<TopoDS_Edge>::iterator aFirstEIt = myEdges.begin();
485   std::list<TopoDS_Edge>::iterator aEIt = aFirstEIt;
486   for (; aEIt != theEdgePosition; ++aEIt);
487
488   // Cut the edge.
489   TopoDS_Edge aEdge = *aEIt;
490   BRepAdaptor_Curve aCurve(aEdge);
491   int aParamI = -1;
492   const double aEdgeEndParams[] =
493     {aCurve.FirstParameter(), aCurve.LastParameter()};
494   if (Abs(theParameter - aEdgeEndParams[0]) < Precision::PConfusion())
495   {
496     aParamI = 0;
497   }
498   else if (Abs(theParameter - aEdgeEndParams[1]) < Precision::PConfusion())
499   {
500     aParamI = 1;
501   }
502   const TopAbs_Orientation aOrient = aEdge.Orientation();
503   if (aOrient == TopAbs_REVERSED)
504   {
505     aParamI ^= 1;
506   }
507   const bool isClosed = IsClosed();
508   DEBTRACE("aParamI: " << aParamI << " isClosed: "<< isClosed);
509   if (aParamI < 0)
510   {
511     aEdge.Orientation(TopAbs_FORWARD);
512     TopoDS_Vertex aSplitV1, aSplitV2;
513     BRep_Builder().MakeVertex(
514       aSplitV1, aCurve.Value(theParameter), Precision::Confusion());
515     BRep_Builder().MakeVertex(
516       aSplitV2, aCurve.Value(theParameter), Precision::Confusion());
517     TopoDS_Edge aEParts[] = {
518       ShapeBuild_Edge().CopyReplaceVertices(aEdge, TopoDS_Vertex(),
519         TopoDS::Vertex(aSplitV1.Oriented(TopAbs_REVERSED))),
520       ShapeBuild_Edge().CopyReplaceVertices(aEdge, aSplitV2, TopoDS_Vertex())};
521     ShapeBuild_Edge().CopyPCurves(aEParts[0], aEdge);
522     ShapeBuild_Edge().CopyPCurves(aEParts[1], aEdge);
523     BRep_Builder().SameRange(aEParts[0], Standard_False);
524     BRep_Builder().SameRange(aEParts[1], Standard_False);
525     BRep_Builder().SameParameter(aEParts[0], Standard_False);
526     BRep_Builder().SameParameter(aEParts[1], Standard_False);
527     ShapeAnalysis_TransferParametersProj aSATPP(aEdge, TopoDS_Face());
528     aSATPP.SetMaxTolerance(Precision::Confusion());
529     aSATPP.TransferRange(aEParts[0],
530       aEdgeEndParams[0], theParameter, Standard_False);
531     aSATPP.TransferRange(aEParts[1],
532       theParameter, aEdgeEndParams[1], Standard_False);
533     aEParts[0].Orientation(aOrient);
534     aEParts[1].Orientation(aOrient);
535
536     const int aFirstPI = (aOrient != TopAbs_REVERSED) ? 0 : 1;
537     *aEIt = aEParts[aFirstPI];
538     InsertAfter(aEIt, aEParts[1 - aFirstPI], myEdges);
539     ++aEIt;
540
541     aResult = true;
542   }
543   else
544   {
545     TopoDS_Edge aNewEdge = ReplaceVertex(aEdge, (aParamI == 0) ? false : true);
546     *aEIt = aNewEdge;
547     if (aParamI > 0)
548     {
549       ++aEIt;
550
551       std::list<TopoDS_Edge>::iterator aEdgePosition = theEdgePosition;
552       if (isClosed || ++aEdgePosition != myEdges.end())
553       {
554         aResult = true;
555       }
556     }
557     else
558     {
559       if (isClosed || theEdgePosition != aFirstEIt)
560       {
561         aResult = true;
562       }
563     }
564   }
565
566   // Calculate the curve parts.
567   std::list<TopoDS_Edge>::iterator aLastEIt = myEdges.end();
568   if (aEIt != aFirstEIt && aEIt != aLastEIt)
569   {
570     std::list<TopoDS_Edge>* aEdges = !isClosed ? &theCurve.myEdges : &myEdges;
571     aEdges->splice(aEdges->begin(), myEdges, aEIt, aLastEIt);
572   }
573
574   return aResult;
575 }
576
577 void HYDROData_TopoCurve::Cut(
578   const std::list<TopoDS_Edge>::const_iterator& theEdgePosition,
579   const double theParameter,
580   HYDROData_TopoCurve& theCurve1,
581   HYDROData_TopoCurve& theCurve2) const
582 {
583   theCurve1 = *this;
584   std::list<TopoDS_Edge>::const_iterator aEPos = myEdges.begin();
585   std::list<TopoDS_Edge>::iterator aEPos1 = theCurve1.myEdges.begin();
586   for (; aEPos != theEdgePosition; ++aEPos1, ++aEPos);
587   theCurve1.Cut(aEPos1, theParameter, theCurve2);
588 }
589
590 bool HYDROData_TopoCurve::Cut(
591   const std::deque<std::list<double> >& theParameters,
592   std::deque<HYDROData_TopoCurve>& theCurves) const
593 {
594   bool aResult = false;
595   HYDROData_TopoCurve aCurves[2];
596   aCurves[0] = *this;
597   int aCI = 0;
598   std::list<TopoDS_Edge>::iterator aEIt = aCurves[0].myEdges.begin();
599   std::deque<std::list<double> >::const_iterator aPLIt = theParameters.begin();
600   for (std::deque<std::list<double> >::const_iterator aLastPLIt =
601     theParameters.end(); aPLIt != aLastPLIt; ++aPLIt)
602   {
603     TopoDS_Edge aNextEdge;
604     {
605       std::list<TopoDS_Edge>::iterator aNextEIt = aEIt;
606       ++aNextEIt;
607       if (aNextEIt != aCurves[aCI].myEdges.end())
608       {
609         aNextEdge = *aNextEIt;
610       }
611     }
612
613     for (Iterator<std::list<double>, std::list<double>::const_iterator> aPIt(
614       *aPLIt, (aEIt->Orientation() != TopAbs_REVERSED)); aPIt.More(); ++aPIt)
615     {
616       const int aCI1 = 1 - aCI;
617       aResult |= aCurves[aCI].Cut(aEIt, **aPIt, aCurves[aCI1]);
618       if (!aCurves[aCI1].IsEmpty())
619       {
620         theCurves.push_back(HYDROData_TopoCurve());
621         theCurves.back().append(aCurves[aCI]);
622         aEIt = aCurves[aCI1].myEdges.begin();
623         aCI = aCI1;
624       }
625       else
626       {
627         aEIt = aCurves[aCI].myEdges.begin();
628       }
629     }
630
631     if (!aNextEdge.IsNull() && !aEIt->IsEqual(aNextEdge))
632     {
633       ++aEIt;
634     }
635   }
636   theCurves.push_back(aCurves[aCI]);
637   return aResult;
638 }
639
640 double HYDROData_TopoCurve::Project(
641   const gp_XYZ& thePoint,
642   std::list<TopoDS_Edge>::const_iterator& theEdgeIterator,
643   double& theParameter) const
644 {
645   double aMinSqDist = DBL_MAX;
646   std::list<TopoDS_Edge>::const_iterator aEIt = myEdges.begin();
647   std::list<TopoDS_Edge>::const_iterator aLastEIt = myEdges.end();
648   for (; aEIt != aLastEIt; ++aEIt)
649   {
650     double aParam;
651     const double aSqDist = ProjectPointToEdge(thePoint, *aEIt, aParam);
652     if (aMinSqDist > aSqDist)
653     {
654       aMinSqDist = aSqDist;
655       theEdgeIterator = aEIt;
656       theParameter = aParam;
657     }
658   }
659   return aMinSqDist;
660 }
661
662 int HYDROData_TopoCurve::Intersect(
663   const TopoDS_Wire& theWire,
664   std::deque<std::list<double> >& theParameters) const
665 {
666   std::string brepName = "theWireToIntersect";
667   brepName += ".brep";
668   BRepTools::Write( theWire, brepName.c_str() );
669
670   int aIntCount = 0;
671   theParameters.resize(myEdges.size());
672   DEBTRACE("myEdges.size() " << myEdges.size());
673   std::list<TopoDS_Edge>::const_iterator aEIt = myEdges.begin();
674   std::list<TopoDS_Edge>::const_iterator aLastEIt = myEdges.end();
675   std::deque<std::list<double> >::iterator aPIt = theParameters.begin();
676   for (; aEIt != aLastEIt; ++aPIt, ++aEIt)
677   {
678     DEBTRACE("---");
679     const TopoDS_Edge& aEdge = *aEIt;
680     std::list<double>& aParams = *aPIt;
681     TopExp_Explorer aEIt2(theWire, TopAbs_EDGE);
682     for (; aEIt2.More(); aEIt2.Next())
683     {
684       aIntCount += IntersectEdge(aEdge,TopoDS::Edge(aEIt2.Current()), aParams);
685       DEBTRACE("aParams.size() " << aParams.size());
686     }
687   }
688   DEBTRACE("aIntCount " << aIntCount);
689   return aIntCount;
690 }
691
692 void HYDROData_TopoCurve::CloseCurve()
693 {
694   const TopoDS_Vertex aVertex = TopoDS::Vertex(TopExp::LastVertex(
695     myEdges.back(), Standard_True).Oriented(TopAbs_FORWARD));
696   TopoDS_Edge& aEdge = myEdges.front();
697   const TopoDS_Edge aForwardEdge = TopoDS::Edge(aEdge.Oriented(TopAbs_FORWARD));
698   aEdge = TopoDS::Edge(ShapeBuild_Edge().CopyReplaceVertices(
699     aForwardEdge, aVertex, TopoDS_Vertex()).Oriented(aEdge.Orientation()));
700 }
701
702 void HYDROData_TopoCurve::Merge(
703   const int thePosition, HYDROData_TopoCurve& theCurve)
704 {
705   if (thePosition == 0)
706   {
707     const TopoDS_Vertex aVertex = TopoDS::Vertex(TopExp::LastVertex(
708       theCurve.myEdges.back(), Standard_True).Oriented(TopAbs_FORWARD));
709     TopoDS_Edge& aEdge = myEdges.front();
710     aEdge = TopoDS::Edge(ShapeBuild_Edge().CopyReplaceVertices(
711         TopoDS::Edge(aEdge.Oriented(TopAbs_FORWARD)), aVertex, TopoDS_Vertex()).
712       Oriented(aEdge.Orientation()));
713     prepend(theCurve);
714   }
715   else
716   {
717     const TopoDS_Vertex aVertex = TopoDS::Vertex(TopExp::FirstVertex(
718       theCurve.myEdges.front(), Standard_True).Oriented(TopAbs_REVERSED));
719     TopoDS_Edge& aEdge = myEdges.back();
720     aEdge = TopoDS::Edge(ShapeBuild_Edge().CopyReplaceVertices(
721         TopoDS::Edge(aEdge.Oriented(TopAbs_FORWARD)), TopoDS_Vertex(), aVertex).
722       Oriented(aEdge.Orientation()));
723     append(theCurve);
724   }
725 }
726
727 void HYDROData_TopoCurve::Merge(
728   const double theTolerance, std::deque<HYDROData_TopoCurve>& theCurves)
729 {
730   // Process the curve closeness.
731   const double aSqTol = theTolerance * theTolerance;
732   const gp_XYZ aPs[] = {
733     BRep_Tool::Pnt(TopExp::FirstVertex(myEdges.front(), Standard_True)).XYZ(),
734     BRep_Tool::Pnt(TopExp::LastVertex(myEdges.back(), Standard_True)).XYZ()};
735   bool isClosed = IsClosed();
736   if (!isClosed && (aPs[0] - aPs[1]).SquareModulus() <= aSqTol)
737   {
738     CloseCurve();
739     isClosed = true;
740   }
741
742   // Find the merge places.
743   HYDROData_TopoCurve* aCurves[] = {NULL, NULL};
744   int aOrder = 0;
745   if (!isClosed)
746   {
747     std::deque<HYDROData_TopoCurve>::iterator aCIt = theCurves.begin();
748     std::deque<HYDROData_TopoCurve>::iterator aLastCIt = theCurves.end();
749     for (; aCIt != aLastCIt; ++aCIt)
750     {
751       HYDROData_TopoCurve& aCurve = *aCIt;
752       if (aCurve.IsEmpty() || aCurve.IsClosed())
753       {
754         continue;
755       }
756
757       const gp_XYZ aP1 = BRep_Tool::Pnt(
758         TopExp::FirstVertex(aCurve.myEdges.front(), Standard_True)).XYZ();
759       if (aCurves[0] == NULL && (aPs[1] - aP1).SquareModulus() <= aSqTol)
760       {
761         aCurves[0] = &aCurve;
762       }
763
764       const gp_XYZ aP2 = BRep_Tool::Pnt(
765         TopExp::LastVertex(aCurve.myEdges.back(), Standard_True)).XYZ();
766       if (aCurves[1] == NULL && (aPs[0] - aP2).SquareModulus() <= aSqTol)
767       {
768         aCurves[1] = &aCurve;
769         aOrder = (aCurves[0] == NULL) ? 1 : 0;
770       }
771     }
772   }
773
774   if (aCurves[0] == NULL && aCurves[1] == NULL)
775   {
776     theCurves.push_back(HYDROData_TopoCurve());
777     theCurves.back().append(*this);
778   }
779   else if (aCurves[1] == NULL)
780   {
781     aCurves[0]->Merge(0, *this);
782   }
783   else if (aCurves[0] == NULL)
784   {
785     aCurves[1]->Merge(1, *this);
786   }
787   else
788   {
789     aCurves[aOrder]->Merge(aOrder, *this);
790     if (aCurves[0] != aCurves[1])
791     {
792       aCurves[aOrder]->Merge(aOrder, *aCurves[1 - aOrder]);
793     }
794     else
795     {
796       aCurves[aOrder]->CloseCurve();
797     }
798   }
799 }
800
801 bool HYDROData_TopoCurve::Connect(
802   const Standard_Real theTolerance, std::deque<HYDROData_TopoCurve>& theCurves)
803 {
804   const double aSqTol = theTolerance * theTolerance;
805   const gp_XYZ aPs[] = {
806     BRep_Tool::Pnt(TopExp::FirstVertex(myEdges.front(), Standard_True)).XYZ(),
807     BRep_Tool::Pnt(TopExp::LastVertex(myEdges.back(), Standard_True)).XYZ()};
808   bool isClosed = IsClosed();
809   if (!isClosed && (aPs[0] - aPs[1]).SquareModulus() <= aSqTol)
810   {
811     CloseCurve();
812     isClosed = true;
813   }
814
815   if (!isClosed)
816   {
817     std::deque<HYDROData_TopoCurve>::iterator aCIt = theCurves.begin();
818     std::deque<HYDROData_TopoCurve>::iterator aLastCIt = theCurves.end();
819     for (; aCIt != aLastCIt; ++aCIt)
820     {
821       HYDROData_TopoCurve& aCurve2 = *aCIt;
822       if (aCurve2.IsEmpty() || aCurve2.IsClosed())
823       {
824         continue;
825       }
826
827       const TopoDS_Edge* aEdges2[] =
828         {&aCurve2.myEdges.front(), &aCurve2.myEdges.back()};
829       const gp_XYZ aPs2[] = {
830         BRep_Tool::Pnt(TopExp::FirstVertex(*aEdges2[0], Standard_True)).XYZ(),
831         BRep_Tool::Pnt(TopExp::LastVertex(*aEdges2[1], Standard_True)).XYZ()};
832       const double aSqDists[] =
833         {(aPs[1] - aPs2[0]).SquareModulus(), (aPs[0] - aPs2[1]).SquareModulus()};
834       const int aOrder = (aSqDists[0] <= aSqDists[1]) ? 0 : 1;
835       if (aSqDists[aOrder] > aSqTol)
836       {
837         const TopoDS_Edge& aEdge =
838           (aOrder == 0) ? myEdges.back() : myEdges.front();
839         const gp_XYZ aPs3[] = {aPs[1 - aOrder], aPs2[aOrder]};
840         const gp_XYZ aTs[] =
841           {Tangent(aEdge, 1 - aOrder), Tangent(*aEdges2[aOrder], aOrder)};
842         Handle(Geom_BSplineCurve) aBSpline;
843         if (!Interpolate(aPs3[aOrder], aPs3[1 - aOrder],
844           aTs[aOrder], aTs[1 - aOrder], aBSpline))
845         {
846           return false;
847         }
848
849         HYDROData_TopoCurve aECurve = BRepBuilderAPI_MakeEdge(aBSpline).Edge();
850         aCurve2.Merge(aOrder, aECurve);
851       }
852       aCurve2.Merge(aOrder, *this);
853       if (aSqDists[1 - aOrder] <= aSqTol)
854       {
855         aCurve2.CloseCurve();
856       }
857       return true;
858     }
859   }
860
861   theCurves.push_back(HYDROData_TopoCurve());
862   theCurves.back().append(*this);
863   return true;
864 }
865
866 int HYDROData_TopoCurve::BSplinePiecewiseCurve(
867   const double theDeflection, HYDROData_TopoCurve& theCurve) const
868 {
869   int aPieceCount = 0;
870   std::list<TopoDS_Edge>::const_iterator aLastEIt = myEdges.end();
871   std::list<TopoDS_Edge>::const_iterator aEIt = myEdges.begin();
872   TopoDS_Vertex aEndVertex;
873   TopoDS_Edge aPrevEdge;
874   for (; aEIt != aLastEIt; ++aEIt)
875   {
876     Handle(Geom_BSplineCurve) aBSpline =
877       ::BSpline(BRepAdaptor_Curve(*aEIt), theDeflection);
878     if (aBSpline.IsNull())
879     {
880       return 0;
881     }
882
883     if (aEIt->Orientation() == TopAbs_REVERSED)
884     {
885       aBSpline->Reverse();
886     }
887
888     TopoDS_Edge aEdge;
889     BRep_Builder().MakeEdge(aEdge, aBSpline, Precision::Confusion());
890     TopoDS_Vertex aVertex;
891     BRep_Builder().MakeVertex(aVertex,
892       aBSpline->Value(aBSpline->FirstParameter()), Precision::Confusion());
893     if (!aPrevEdge.IsNull())
894     {
895       BRep_Builder().Add(aPrevEdge, aVertex.Oriented(TopAbs_REVERSED));
896     }
897     else
898     {
899       aEndVertex = aVertex;
900     }
901     BRep_Builder().Add(aEdge, aVertex);
902     theCurve.myEdges.push_back(aEdge);
903     aPieceCount += aBSpline->NbKnots() - 1;
904     aPrevEdge = aEdge;
905   }
906
907   if (!IsClosed())
908   {
909     BRepAdaptor_Curve aCurve(aPrevEdge);
910     BRep_Builder().MakeVertex(aEndVertex,
911       aCurve.Value(aCurve.LastParameter()), Precision::Confusion());
912   }
913   BRep_Builder().Add(aPrevEdge, aEndVertex.Oriented(TopAbs_REVERSED));
914   return aPieceCount;
915 }
916
917 bool HYDROData_TopoCurve::ValuesInKnots(std::list<gp_XYZ>& theValues) const
918 {
919   std::list<TopoDS_Edge>::const_iterator aLastEIt = myEdges.end();
920   std::list<TopoDS_Edge>::const_iterator aEIt = myEdges.begin();
921   for (; aEIt != aLastEIt; ++aEIt)
922   {
923     Handle(Geom_BSplineCurve) aCurve;
924     {
925       TopLoc_Location aLoc;
926       double aParams[2];
927       aCurve = Handle(Geom_BSplineCurve)::
928         DownCast(BRep_Tool::Curve(*aEIt, aLoc, aParams[0], aParams[1]));
929       if (!aLoc.IsIdentity() || aEIt->Orientation() != TopAbs_FORWARD ||
930         aCurve.IsNull())
931       {
932         return false;
933       }
934     }
935
936     for (int aNbKnots = aCurve->NbKnots(), aKN = 1; aKN < aNbKnots; ++aKN)
937     {
938       theValues.push_back(aCurve->Value(aCurve->Knot(aKN)).XYZ());
939     }
940   }
941
942   if (!IsClosed())
943   {
944     TopLoc_Location aLoc;
945     double aParams[2];
946     Handle(Geom_Curve) aCurve =
947       BRep_Tool::Curve(myEdges.back(), aLoc, aParams[0], aParams[1]);
948     theValues.push_back(aCurve->Value(aParams[1]).XYZ());
949   }
950   return true;
951 }