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();   // number of edges (an edge can be a spline with intermediate vertices)
 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   std::string brepName = "theWireToSplit";
 408   brepName += ".brep";
 409   BRepTools::Write( theWire, brepName.c_str() );
 410
 411   // Find the start, i.e. the end which is the first vertex (or just the first vertex when closed) ==> aVN
 412   int aVN = 1;
 413   if (!isClosed)
 414   {
 415     for (; aVertexToEdges.FindFromIndex(aVN).Extent() == 2; ++aVN);
 416     if (!aVertexToEdges.FindKey(aVN).IsEqual(TopExp::FirstVertex(
 417       TopoDS::Edge(aVertexToEdges.FindFromIndex(aVN).First()), Standard_True)))
 418     {
 419       for (; aVertexToEdges.FindFromIndex(aVN).Extent() == 2; ++aVN);
 420     }
 421   }
 422   else
 423   {
 424     TopTools_ListOfShape& aEdges = aVertexToEdges.ChangeFromIndex(1);
 425     if (!aVertexToEdges.FindKey(aVN).IsEqual(
 426       TopExp::FirstVertex(TopoDS::Edge(aEdges.First()), Standard_True)))
 427     {
 428       const TopoDS_Shape aEdge = aEdges.First();
 429       aEdges.First() = aEdges.Last();
 430       aEdges.Last() = aEdge;
 431     }
 432   }
 433
 434   // Calculate the edge order ==> list in myEdges
 435   TopTools_ListOfShape* aEdges = &aVertexToEdges.ChangeFromIndex(aVN); // 1 or 2 edges from first vertex
 436   while (!aEdges->IsEmpty())
 437   {
 438     const TopoDS_Edge aEdge = TopoDS::Edge(aEdges->First());
 439     aEdges->RemoveFirst();
 440     myEdges.push_back(aEdge);                                          // add an edge in the list
 441     int aVN2 = aVertexToEdges.FindIndex(TopExp::FirstVertex(aEdge));
 442     if (aVN2 == aVN)
 443     {
 444       aVN2 = aVertexToEdges.FindIndex(TopExp::LastVertex(aEdge));
 445     }
 446     aVN = aVN2;                                                        // next vertex
 447
 448     aEdges = &aVertexToEdges.ChangeFromIndex(aVN2);
 449     const TopoDS_Edge aEdge2 = TopoDS::Edge(aEdges->First());
 450     if (aEdge2.IsEqual(aEdge))
 451     {
 452       aEdges->RemoveFirst();
 453     }
 454     else
 455     {
 456       aEdges->Clear();
 457       aEdges->Append(aEdge2);                                          // next edge
 458     }
 459   }
 460
 461   // Check for connectedness and free vertex.
 462   return (aVCount - myEdges.size()) == (isClosed ? 0 : 1);
 463 }
 464
 465 TopoDS_Wire HYDROData_TopoCurve::Wire() const
 466 {
 467   TopoDS_Wire aWire;
 468   BRep_Builder aBuilder;
 469   aBuilder.MakeWire(aWire);
 470   std::list<TopoDS_Edge>::const_iterator aEItLast = myEdges.end();
 471   std::list<TopoDS_Edge>::const_iterator aEIt = myEdges.begin();
 472   for (; aEIt != aEItLast; ++aEIt)
 473   {
 474     aBuilder.Add(aWire, *aEIt);
 475   }
 476   return aWire;
 477 }
 478
 479 bool HYDROData_TopoCurve::Cut(
 480   const std::list<TopoDS_Edge>::iterator& theEdgePosition,
 481   const double theParameter,
 482   HYDROData_TopoCurve& theCurve)
 483 {
 484   bool aResult = false;
 485
 486   // Locate the edge.
 487   std::list<TopoDS_Edge>::iterator aFirstEIt = myEdges.begin();
 488   std::list<TopoDS_Edge>::iterator aEIt = aFirstEIt;
 489   for (; aEIt != theEdgePosition; ++aEIt);
 490
 491   // Cut the edge.
 492   TopoDS_Edge aEdge = *aEIt;
 493   BRepAdaptor_Curve aCurve(aEdge);
 494   int aParamI = -1;
 495   const double aEdgeEndParams[] =
 496     {aCurve.FirstParameter(), aCurve.LastParameter()};
 497   if (Abs(theParameter - aEdgeEndParams[0]) < Precision::PConfusion())
 498   {
 499     aParamI = 0;
 500   }
 501   else if (Abs(theParameter - aEdgeEndParams[1]) < Precision::PConfusion())
 502   {
 503     aParamI = 1;
 504   }
 505   const TopAbs_Orientation aOrient = aEdge.Orientation();
 506   if (aOrient == TopAbs_REVERSED)
 507   {
 508     aParamI ^= 1;
 509   }
 510   const bool isClosed = IsClosed();
 511   DEBTRACE("aParamI: " << aParamI << " isClosed: "<< isClosed);
 512   if (aParamI < 0)
 513   {
 514     aEdge.Orientation(TopAbs_FORWARD);
 515     TopoDS_Vertex aSplitV1, aSplitV2;
 516     BRep_Builder().MakeVertex(
 517       aSplitV1, aCurve.Value(theParameter), Precision::Confusion());
 518     BRep_Builder().MakeVertex(
 519       aSplitV2, aCurve.Value(theParameter), Precision::Confusion());
 520     TopoDS_Edge aEParts[] = {
 521       ShapeBuild_Edge().CopyReplaceVertices(aEdge, TopoDS_Vertex(),
 522         TopoDS::Vertex(aSplitV1.Oriented(TopAbs_REVERSED))),
 523       ShapeBuild_Edge().CopyReplaceVertices(aEdge, aSplitV2, TopoDS_Vertex())};
 524     ShapeBuild_Edge().CopyPCurves(aEParts[0], aEdge);
 525     ShapeBuild_Edge().CopyPCurves(aEParts[1], aEdge);
 526     BRep_Builder().SameRange(aEParts[0], Standard_False);
 527     BRep_Builder().SameRange(aEParts[1], Standard_False);
 528     BRep_Builder().SameParameter(aEParts[0], Standard_False);
 529     BRep_Builder().SameParameter(aEParts[1], Standard_False);
 530     ShapeAnalysis_TransferParametersProj aSATPP(aEdge, TopoDS_Face());
 531     aSATPP.SetMaxTolerance(Precision::Confusion());
 532     aSATPP.TransferRange(aEParts[0],
 533       aEdgeEndParams[0], theParameter, Standard_False);
 534     aSATPP.TransferRange(aEParts[1],
 535       theParameter, aEdgeEndParams[1], Standard_False);
 536     aEParts[0].Orientation(aOrient);
 537     aEParts[1].Orientation(aOrient);
 538
 539     const int aFirstPI = (aOrient != TopAbs_REVERSED) ? 0 : 1;
 540     *aEIt = aEParts[aFirstPI];
 541     InsertAfter(aEIt, aEParts[1 - aFirstPI], myEdges);
 542     ++aEIt;
 543
 544     aResult = true;
 545   }
 546   else
 547   {
 548     TopoDS_Edge aNewEdge = ReplaceVertex(aEdge, (aParamI == 0) ? false : true);
 549     *aEIt = aNewEdge;
 550     if (aParamI > 0)
 551     {
 552       ++aEIt;
 553
 554       std::list<TopoDS_Edge>::iterator aEdgePosition = theEdgePosition;
 555       if (isClosed || ++aEdgePosition != myEdges.end())
 556       {
 557         aResult = true;
 558       }
 559     }
 560     else
 561     {
 562       if (isClosed || theEdgePosition != aFirstEIt)
 563       {
 564         aResult = true;
 565       }
 566     }
 567   }
 568
 569   // Calculate the curve parts.
 570   std::list<TopoDS_Edge>::iterator aLastEIt = myEdges.end();
 571   if (aEIt != aFirstEIt && aEIt != aLastEIt)
 572   {
 573     std::list<TopoDS_Edge>* aEdges = !isClosed ? &theCurve.myEdges : &myEdges;
 574     aEdges->splice(aEdges->begin(), myEdges, aEIt, aLastEIt);
 575   }
 576
 577   return aResult;
 578 }
 579
 580 void HYDROData_TopoCurve::Cut(
 581   const std::list<TopoDS_Edge>::const_iterator& theEdgePosition,
 582   const double theParameter,
 583   HYDROData_TopoCurve& theCurve1,
 584   HYDROData_TopoCurve& theCurve2) const
 585 {
 586   theCurve1 = *this;
 587   std::list<TopoDS_Edge>::const_iterator aEPos = myEdges.begin();
 588   std::list<TopoDS_Edge>::iterator aEPos1 = theCurve1.myEdges.begin();
 589   for (; aEPos != theEdgePosition; ++aEPos1, ++aEPos);
 590   theCurve1.Cut(aEPos1, theParameter, theCurve2);
 591 }
 592
 593 bool HYDROData_TopoCurve::Cut(
 594   const std::deque<std::list<double> >& theParameters,
 595   std::deque<HYDROData_TopoCurve>& theCurves) const
 596 {
 597   bool aResult = false;
 598   HYDROData_TopoCurve aCurves[2];
 599   aCurves[0] = *this;
 600   int aCI = 0;
 601   std::list<TopoDS_Edge>::iterator aEIt = aCurves[0].myEdges.begin();
 602   std::deque<std::list<double> >::const_iterator aPLIt = theParameters.begin();
 603   for (std::deque<std::list<double> >::const_iterator aLastPLIt =
 604     theParameters.end(); aPLIt != aLastPLIt; ++aPLIt)
 605   {
 606     TopoDS_Edge aNextEdge;
 607     {
 608       std::list<TopoDS_Edge>::iterator aNextEIt = aEIt;
 609       ++aNextEIt;
 610       if (aNextEIt != aCurves[aCI].myEdges.end())
 611       {
 612         aNextEdge = *aNextEIt;
 613       }
 614     }
 615
 616     for (Iterator<std::list<double>, std::list<double>::const_iterator> aPIt(
 617       *aPLIt, (aEIt->Orientation() != TopAbs_REVERSED)); aPIt.More(); ++aPIt)
 618     {
 619       const int aCI1 = 1 - aCI;
 620       aResult |= aCurves[aCI].Cut(aEIt, **aPIt, aCurves[aCI1]);
 621       if (!aCurves[aCI1].IsEmpty())
 622       {
 623         theCurves.push_back(HYDROData_TopoCurve());
 624         theCurves.back().append(aCurves[aCI]);
 625         aEIt = aCurves[aCI1].myEdges.begin();
 626         aCI = aCI1;
 627       }
 628       else
 629       {
 630         aEIt = aCurves[aCI].myEdges.begin();
 631       }
 632     }
 633
 634     if (!aNextEdge.IsNull() && !aEIt->IsEqual(aNextEdge))
 635     {
 636       ++aEIt;
 637     }
 638   }
 639   theCurves.push_back(aCurves[aCI]);
 640   return aResult;
 641 }
 642
 643 double HYDROData_TopoCurve::Project(
 644   const gp_XYZ& thePoint,
 645   std::list<TopoDS_Edge>::const_iterator& theEdgeIterator,
 646   double& theParameter) const
 647 {
 648   double aMinSqDist = DBL_MAX;
 649   std::list<TopoDS_Edge>::const_iterator aEIt = myEdges.begin();
 650   std::list<TopoDS_Edge>::const_iterator aLastEIt = myEdges.end();
 651   for (; aEIt != aLastEIt; ++aEIt)
 652   {
 653     double aParam;
 654     const double aSqDist = ProjectPointToEdge(thePoint, *aEIt, aParam);
 655     if (aMinSqDist > aSqDist)
 656     {
 657       aMinSqDist = aSqDist;
 658       theEdgeIterator = aEIt;
 659       theParameter = aParam;
 660     }
 661   }
 662   return aMinSqDist;
 663 }
 664
 665 int HYDROData_TopoCurve::Intersect(
 666   const TopoDS_Wire& theWire,
 667   std::deque<std::list<double> >& theParameters) const
 668 {
 669   std::string brepName = "theWireTool";
 670   brepName += ".brep";
 671   BRepTools::Write( theWire, brepName.c_str() );
 672
 673   int aIntCount = 0;
 674   theParameters.resize(myEdges.size());
 675   DEBTRACE("myEdges.size() " << myEdges.size());
 676   std::list<TopoDS_Edge>::const_iterator aEIt = myEdges.begin();
 677   std::list<TopoDS_Edge>::const_iterator aLastEIt = myEdges.end();
 678   std::deque<std::list<double> >::iterator aPIt = theParameters.begin();
 679   for (; aEIt != aLastEIt; ++aPIt, ++aEIt)
 680   {
 681     DEBTRACE("---");
 682     const TopoDS_Edge& aEdge = *aEIt;
 683     std::list<double>& aParams = *aPIt;
 684     TopExp_Explorer aEIt2(theWire, TopAbs_EDGE);
 685     for (; aEIt2.More(); aEIt2.Next())
 686     {
 687       aIntCount += IntersectEdge(aEdge,TopoDS::Edge(aEIt2.Current()), aParams);
 688       DEBTRACE("aParams.size() " << aParams.size());
 689     }
 690   }
 691   DEBTRACE("aIntCount " << aIntCount);
 692   return aIntCount;
 693 }
 694
 695 void HYDROData_TopoCurve::CloseCurve()
 696 {
 697   const TopoDS_Vertex aVertex = TopoDS::Vertex(TopExp::LastVertex(
 698     myEdges.back(), Standard_True).Oriented(TopAbs_FORWARD));
 699   TopoDS_Edge& aEdge = myEdges.front();
 700   const TopoDS_Edge aForwardEdge = TopoDS::Edge(aEdge.Oriented(TopAbs_FORWARD));
 701   aEdge = TopoDS::Edge(ShapeBuild_Edge().CopyReplaceVertices(
 702     aForwardEdge, aVertex, TopoDS_Vertex()).Oriented(aEdge.Orientation()));
 703 }
 704
 705 void HYDROData_TopoCurve::Merge(
 706   const int thePosition, HYDROData_TopoCurve& theCurve)
 707 {
 708   if (thePosition == 0)
 709   {
 710     const TopoDS_Vertex aVertex = TopoDS::Vertex(TopExp::LastVertex(
 711       theCurve.myEdges.back(), Standard_True).Oriented(TopAbs_FORWARD));
 712     TopoDS_Edge& aEdge = myEdges.front();
 713     aEdge = TopoDS::Edge(ShapeBuild_Edge().CopyReplaceVertices(
 714         TopoDS::Edge(aEdge.Oriented(TopAbs_FORWARD)), aVertex, TopoDS_Vertex()).
 715       Oriented(aEdge.Orientation()));
 716     prepend(theCurve);
 717   }
 718   else
 719   {
 720     const TopoDS_Vertex aVertex = TopoDS::Vertex(TopExp::FirstVertex(
 721       theCurve.myEdges.front(), Standard_True).Oriented(TopAbs_REVERSED));
 722     TopoDS_Edge& aEdge = myEdges.back();
 723     aEdge = TopoDS::Edge(ShapeBuild_Edge().CopyReplaceVertices(
 724         TopoDS::Edge(aEdge.Oriented(TopAbs_FORWARD)), TopoDS_Vertex(), aVertex).
 725       Oriented(aEdge.Orientation()));
 726     append(theCurve);
 727   }
 728 }
 729
 730 void HYDROData_TopoCurve::Merge(
 731   const double theTolerance, std::deque<HYDROData_TopoCurve>& theCurves)
 732 {
 733   // Process the curve closeness.
 734   const double aSqTol = theTolerance * theTolerance;
 735   const gp_XYZ aPs[] = {
 736     BRep_Tool::Pnt(TopExp::FirstVertex(myEdges.front(), Standard_True)).XYZ(),
 737     BRep_Tool::Pnt(TopExp::LastVertex(myEdges.back(), Standard_True)).XYZ()};
 738   bool isClosed = IsClosed();
 739   if (!isClosed && (aPs[0] - aPs[1]).SquareModulus() <= aSqTol)
 740   {
 741     CloseCurve();
 742     isClosed = true;
 743   }
 744
 745   // Find the merge places.
 746   HYDROData_TopoCurve* aCurves[] = {NULL, NULL};
 747   int aOrder = 0;
 748   if (!isClosed)
 749   {
 750     std::deque<HYDROData_TopoCurve>::iterator aCIt = theCurves.begin();
 751     std::deque<HYDROData_TopoCurve>::iterator aLastCIt = theCurves.end();
 752     for (; aCIt != aLastCIt; ++aCIt)
 753     {
 754       HYDROData_TopoCurve& aCurve = *aCIt;
 755       if (aCurve.IsEmpty() || aCurve.IsClosed())
 756       {
 757         continue;
 758       }
 759
 760       const gp_XYZ aP1 = BRep_Tool::Pnt(
 761         TopExp::FirstVertex(aCurve.myEdges.front(), Standard_True)).XYZ();
 762       if (aCurves[0] == NULL && (aPs[1] - aP1).SquareModulus() <= aSqTol)
 763       {
 764         aCurves[0] = &aCurve;
 765       }
 766
 767       const gp_XYZ aP2 = BRep_Tool::Pnt(
 768         TopExp::LastVertex(aCurve.myEdges.back(), Standard_True)).XYZ();
 769       if (aCurves[1] == NULL && (aPs[0] - aP2).SquareModulus() <= aSqTol)
 770       {
 771         aCurves[1] = &aCurve;
 772         aOrder = (aCurves[0] == NULL) ? 1 : 0;
 773       }
 774     }
 775   }
 776
 777   if (aCurves[0] == NULL && aCurves[1] == NULL)
 778   {
 779     theCurves.push_back(HYDROData_TopoCurve());
 780     theCurves.back().append(*this);
 781   }
 782   else if (aCurves[1] == NULL)
 783   {
 784     aCurves[0]->Merge(0, *this);
 785   }
 786   else if (aCurves[0] == NULL)
 787   {
 788     aCurves[1]->Merge(1, *this);
 789   }
 790   else
 791   {
 792     aCurves[aOrder]->Merge(aOrder, *this);
 793     if (aCurves[0] != aCurves[1])
 794     {
 795       aCurves[aOrder]->Merge(aOrder, *aCurves[1 - aOrder]);
 796     }
 797     else
 798     {
 799       aCurves[aOrder]->CloseCurve();
 800     }
 801   }
 802 }
 803
 804 bool HYDROData_TopoCurve::Connect(
 805   const Standard_Real theTolerance, std::deque<HYDROData_TopoCurve>& theCurves)
 806 {
 807   const double aSqTol = theTolerance * theTolerance;
 808   const gp_XYZ aPs[] = {
 809     BRep_Tool::Pnt(TopExp::FirstVertex(myEdges.front(), Standard_True)).XYZ(),
 810     BRep_Tool::Pnt(TopExp::LastVertex(myEdges.back(), Standard_True)).XYZ()};
 811   bool isClosed = IsClosed();
 812   if (!isClosed && (aPs[0] - aPs[1]).SquareModulus() <= aSqTol)
 813   {
 814     CloseCurve();
 815     isClosed = true;
 816   }
 817
 818   if (!isClosed)
 819   {
 820     std::deque<HYDROData_TopoCurve>::iterator aCIt = theCurves.begin();
 821     std::deque<HYDROData_TopoCurve>::iterator aLastCIt = theCurves.end();
 822     for (; aCIt != aLastCIt; ++aCIt)
 823     {
 824       HYDROData_TopoCurve& aCurve2 = *aCIt;
 825       if (aCurve2.IsEmpty() || aCurve2.IsClosed())
 826       {
 827         continue;
 828       }
 829
 830       const TopoDS_Edge* aEdges2[] =
 831         {&aCurve2.myEdges.front(), &aCurve2.myEdges.back()};
 832       const gp_XYZ aPs2[] = {
 833         BRep_Tool::Pnt(TopExp::FirstVertex(*aEdges2[0], Standard_True)).XYZ(),
 834         BRep_Tool::Pnt(TopExp::LastVertex(*aEdges2[1], Standard_True)).XYZ()};
 835       const double aSqDists[] =
 836         {(aPs[1] - aPs2[0]).SquareModulus(), (aPs[0] - aPs2[1]).SquareModulus()};
 837       const int aOrder = (aSqDists[0] <= aSqDists[1]) ? 0 : 1;
 838       if (aSqDists[aOrder] > aSqTol)
 839       {
 840         const TopoDS_Edge& aEdge =
 841           (aOrder == 0) ? myEdges.back() : myEdges.front();
 842         const gp_XYZ aPs3[] = {aPs[1 - aOrder], aPs2[aOrder]};
 843         const gp_XYZ aTs[] =
 844           {Tangent(aEdge, 1 - aOrder), Tangent(*aEdges2[aOrder], aOrder)};
 845         Handle(Geom_BSplineCurve) aBSpline;
 846         if (!Interpolate(aPs3[aOrder], aPs3[1 - aOrder],
 847           aTs[aOrder], aTs[1 - aOrder], aBSpline))
 848         {
 849           return false;
 850         }
 851
 852         HYDROData_TopoCurve aECurve = BRepBuilderAPI_MakeEdge(aBSpline).Edge();
 853         aCurve2.Merge(aOrder, aECurve);
 854       }
 855       aCurve2.Merge(aOrder, *this);
 856       if (aSqDists[1 - aOrder] <= aSqTol)
 857       {
 858         aCurve2.CloseCurve();
 859       }
 860       return true;
 861     }
 862   }
 863
 864   theCurves.push_back(HYDROData_TopoCurve());
 865   theCurves.back().append(*this);
 866   return true;
 867 }
 868
 869 int HYDROData_TopoCurve::BSplinePiecewiseCurve(
 870   const double theDeflection, HYDROData_TopoCurve& theCurve) const
 871 {
 872   int aPieceCount = 0;
 873   std::list<TopoDS_Edge>::const_iterator aLastEIt = myEdges.end();
 874   std::list<TopoDS_Edge>::const_iterator aEIt = myEdges.begin();
 875   TopoDS_Vertex aEndVertex;
 876   TopoDS_Edge aPrevEdge;
 877   for (; aEIt != aLastEIt; ++aEIt)
 878   {
 879     Handle(Geom_BSplineCurve) aBSpline =
 880       ::BSpline(BRepAdaptor_Curve(*aEIt), theDeflection);
 881     if (aBSpline.IsNull())
 882     {
 883       return 0;
 884     }
 885
 886     if (aEIt->Orientation() == TopAbs_REVERSED)
 887     {
 888       aBSpline->Reverse();
 889     }
 890
 891     TopoDS_Edge aEdge;
 892     BRep_Builder().MakeEdge(aEdge, aBSpline, Precision::Confusion());
 893     TopoDS_Vertex aVertex;
 894     BRep_Builder().MakeVertex(aVertex,
 895       aBSpline->Value(aBSpline->FirstParameter()), Precision::Confusion());
 896     if (!aPrevEdge.IsNull())
 897     {
 898       BRep_Builder().Add(aPrevEdge, aVertex.Oriented(TopAbs_REVERSED));
 899     }
 900     else
 901     {
 902       aEndVertex = aVertex;
 903     }
 904     BRep_Builder().Add(aEdge, aVertex);
 905     theCurve.myEdges.push_back(aEdge);
 906     aPieceCount += aBSpline->NbKnots() - 1;
 907     aPrevEdge = aEdge;
 908   }
 909
 910   if (!IsClosed())
 911   {
 912     BRepAdaptor_Curve aCurve(aPrevEdge);
 913     BRep_Builder().MakeVertex(aEndVertex,
 914       aCurve.Value(aCurve.LastParameter()), Precision::Confusion());
 915   }
 916   BRep_Builder().Add(aPrevEdge, aEndVertex.Oriented(TopAbs_REVERSED));
 917   return aPieceCount;
 918 }
 919
 920 bool HYDROData_TopoCurve::ValuesInKnots(std::list<gp_XYZ>& theValues) const
 921 {
 922   std::list<TopoDS_Edge>::const_iterator aLastEIt = myEdges.end();
 923   std::list<TopoDS_Edge>::const_iterator aEIt = myEdges.begin();
 924   for (; aEIt != aLastEIt; ++aEIt)
 925   {
 926     Handle(Geom_BSplineCurve) aCurve;
 927     {
 928       TopLoc_Location aLoc;
 929       double aParams[2];
 930       aCurve = Handle(Geom_BSplineCurve)::
 931         DownCast(BRep_Tool::Curve(*aEIt, aLoc, aParams[0], aParams[1]));
 932       if (!aLoc.IsIdentity() || aEIt->Orientation() != TopAbs_FORWARD ||
 933         aCurve.IsNull())
 934       {
 935         return false;
 936       }
 937     }
 938
 939     for (int aNbKnots = aCurve->NbKnots(), aKN = 1; aKN < aNbKnots; ++aKN)
 940     {
 941       theValues.push_back(aCurve->Value(aCurve->Knot(aKN)).XYZ());
 942     }
 943   }
 944
 945   if (!IsClosed())
 946   {
 947     TopLoc_Location aLoc;
 948     double aParams[2];
 949     Handle(Geom_Curve) aCurve =
 950       BRep_Tool::Curve(myEdges.back(), aLoc, aParams[0], aParams[1]);
 951     theValues.push_back(aCurve->Value(aParams[1]).XYZ());
 952   }
 953   return true;
 954 }