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 }