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