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