src/HYDROData/HYDROData_PolylineOperator.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_PolylineOperator.h>
  20 #include <HYDROData_Document.h>
  21 #include <BRepAdaptor_Curve.hxx>
  22 #include <BRep_Builder.hxx>
  23 #include <BRep_Tool.hxx>
  24 #include <BRepBuilderAPI_MakeEdge2d.hxx>
  25 #include <BRepBuilderAPI_MakeEdge.hxx>
  26 #include <BRepBuilderAPI_MakeWire.hxx>
  27 #include <Extrema_ExtCC.hxx>
  28 #include <Extrema_ExtPC.hxx>
  29 #include <GeomAPI_Interpolate.hxx>
  30 #include <NCollection_Vector.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 <TopoDS.hxx>
  38 #include <TopoDS_Edge.hxx>
  39 #include <TopoDS_Wire.hxx>
  40 #include <TopExp.hxx>
  41 #include <TopExp_Explorer.hxx>
  42 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
  43 #include <QString>
  44
  45 template<class T> void append( std::vector<T>& theList, const std::vector<T>& theList2 )
  46 {
  47   int aSize = theList.size();
  48   int aNewSize = aSize + theList2.size();
  49
  50   if( aSize==aNewSize )
  51     return;
  52
  53   theList.resize( aNewSize );
  54   for( int i=aSize, j=0; i<aNewSize; i++, j++ )
  55     theList[i] = theList2[j];
  56 }
  57
  58 static TopoDS_Edge ReplaceVertex(
  59   const TopoDS_Edge& theEdge, Standard_Integer theVertexIndex)
  60 {
  61   TopoDS_Vertex aVertices[] =
  62     {TopExp::FirstVertex(theEdge), TopExp::LastVertex(theEdge)};
  63   const TopAbs_Orientation aOrient = theEdge.Orientation();
  64   if (aOrient == TopAbs_REVERSED)
  65   {
  66     theVertexIndex ^= 1;
  67   }
  68   aVertices[theVertexIndex].EmptyCopy();
  69   aVertices[0].Orientation(TopAbs_FORWARD);
  70   aVertices[1].Orientation(TopAbs_REVERSED);
  71   TopoDS_Edge aFE = TopoDS::Edge(theEdge.Oriented(TopAbs_FORWARD));
  72   TopoDS_Edge aNewEdge =
  73     ShapeBuild_Edge().CopyReplaceVertices(aFE, aVertices[0], aVertices[1]);
  74   return TopoDS::Edge(aNewEdge.Oriented(aOrient));
  75 }
  76
  77 template<typename TCurveType>
  78 static Standard_Boolean WireToCurve(
  79   const TopoDS_Wire& theWire,
  80   TCurveType& theCurve,
  81   Standard_Boolean& theIsClosed)
  82 {
  83   TopTools_IndexedDataMapOfShapeListOfShape aVertexToEdges;
  84   TopExp::MapShapesAndAncestors(theWire,
  85     TopAbs_VERTEX, TopAbs_EDGE, aVertexToEdges);
  86   const Standard_Integer aVCount = aVertexToEdges.Extent();
  87   if (aVCount == 0)
  88   {
  89     return Standard_False;
  90   }
  91
  92   {
  93     Standard_Integer aEndCount = 0;
  94     for (Standard_Integer aVN = 1; aVN <= aVCount; ++aVN)
  95     {
  96       const Standard_Integer aEdgeCount =
  97         aVertexToEdges.FindFromIndex(aVN).Extent();
  98       if (aEdgeCount == 1)
  99       {
 100         ++aEndCount;
 101       }
 102       if (aEdgeCount > 2)
 103       {
 104         return Standard_False;
 105       }
 106     }
 107     theIsClosed = (aEndCount == 0);
 108     if (!theIsClosed && aEndCount != 2)
 109     {
 110       return Standard_False;
 111     }
 112   }
 113
 114   Standard_Integer aVN = 1;
 115   if (!theIsClosed)
 116   {
 117     while (aVN <= aVCount &&
 118       aVertexToEdges.FindFromIndex(aVN).Extent() == 2)
 119     {
 120       ++aVN;
 121     }
 122   }
 123
 124   TopTools_ListOfShape* aEdges = &aVertexToEdges.ChangeFromIndex(aVN);
 125   while (!aEdges->IsEmpty())
 126   {
 127     const TopoDS_Edge aEdge = TopoDS::Edge(aEdges->First());
 128     aEdges->RemoveFirst();
 129     theCurve.Append(aEdge);
 130     Standard_Integer aVN2 =
 131       aVertexToEdges.FindIndex(TopExp::FirstVertex(aEdge));
 132     if (aVN2 == aVN)
 133     {
 134       aVN2 = aVertexToEdges.FindIndex(TopExp::LastVertex(aEdge));
 135     }
 136     aVN = aVN2;
 137
 138     aEdges = &aVertexToEdges.ChangeFromIndex(aVN2);
 139     if (aEdges->First().IsEqual(aEdge))
 140     {
 141       aEdges->RemoveFirst();
 142     }
 143     else
 144     {
 145       const TopoDS_Edge aEdge = TopoDS::Edge(aEdges->First());
 146       aEdges->Clear();
 147       aEdges->Append(aEdge);
 148     }
 149   }
 150   return (!theIsClosed && theCurve.Size() == aVCount - 1) ||
 151     (theIsClosed && theCurve.Size() == aVCount);
 152 }
 153
 154 template<typename TCurveType>
 155 static void CurveToWire(const TCurveType& theCurve, TopoDS_Wire& theWire)
 156 {
 157   BRep_Builder aBulder;
 158   aBulder.MakeWire(theWire);
 159   for (typename TCurveType::Iterator aEIt(theCurve); aEIt.More(); aEIt.Next())
 160   {
 161     aBulder.Add(theWire, aEIt.Value());
 162   }
 163 }
 164
 165 static void CurveToWire(
 166   const NCollection_Vector<TopoDS_Edge>& theEdges1,
 167   const NCollection_Vector<TopoDS_Edge>& theEdges2,
 168   TopoDS_Wire& theWire)
 169 {
 170   BRep_Builder aBulder;
 171   aBulder.MakeWire(theWire);
 172   const NCollection_Vector<TopoDS_Edge>* aEdges[] = {&theEdges1, &theEdges2};
 173   for (Standard_Integer aEI = 0; aEI < 2; ++aEI)
 174   {
 175     NCollection_Vector<TopoDS_Edge>::Iterator aEIt(*aEdges[aEI]);
 176     for (; aEIt.More(); aEIt.Next())
 177     {
 178       aBulder.Add(theWire, aEIt.Value());
 179     }
 180   }
 181 }
 182
 183 static Standard_Real ProjectPointToCurve(
 184   const gp_Pnt& thePoint,
 185   const Adaptor3d_Curve& theCurve,
 186   Standard_Real& theParameter)
 187 {
 188   Extrema_ExtPC aAlgo(thePoint, theCurve);
 189   Standard_Integer aMinEN = -2;
 190   Standard_Real aMinSqDist = DBL_MAX;
 191   if (aAlgo.IsDone())
 192   {
 193     const Standard_Integer aECount = aAlgo.NbExt();
 194     for (Standard_Integer aEN = 1; aEN <= aECount; ++aEN)
 195     {
 196       const gp_Pnt& aP = aAlgo.Point(aEN).Value();
 197       const Standard_Real aSqDist = thePoint.SquareDistance(aP);
 198       if (aMinSqDist > aSqDist)
 199       {
 200         aMinSqDist = aSqDist;
 201         aMinEN = aEN;
 202       }
 203     }
 204   }
 205
 206   const Standard_Real aParams[] =
 207     {theCurve.FirstParameter(), theCurve.LastParameter()};
 208   const gp_Pnt aEnds[] =
 209     {theCurve.Value(aParams[0]), theCurve.Value(aParams[1])};
 210   const Standard_Real aSqDists[] =
 211     {thePoint.SquareDistance(aEnds[0]), thePoint.SquareDistance(aEnds[1])};
 212   for (Standard_Integer aEI = 0; aEI < 2; ++aEI)
 213   {
 214     if (aMinSqDist > aSqDists[aEI])
 215     {
 216       aMinSqDist = aSqDists[aEI];
 217       aMinEN = -aEI;
 218     }
 219   }
 220
 221   if (aMinEN <= 0)
 222   {
 223     theParameter = aParams[-aMinEN];
 224     return aMinSqDist;
 225   }
 226
 227   const Extrema_POnCurv& aPOnC = aAlgo.Point(aMinEN);
 228   const gp_Pnt& aPoint = aPOnC.Value();
 229   theParameter = aPOnC.Parameter();
 230   for (Standard_Integer aEI = 0; aEI < 2; ++aEI)
 231   {
 232     if (Abs(theParameter - aParams[aEI]) < Precision::PConfusion() ||
 233       aPoint.SquareDistance(aEnds[aEI]) < Precision::SquareConfusion())
 234     {
 235       theParameter = aParams[aEI];
 236     }
 237   }
 238   return aMinSqDist;
 239 }
 240
 241 static Standard_Real ProjectPointToEdge(
 242   const gp_Pnt& thePoint,
 243   const TopoDS_Edge& theEdge,
 244   Standard_Real& theParameter)
 245 {
 246   return ProjectPointToCurve(thePoint, BRepAdaptor_Curve(theEdge), theParameter);
 247 }
 248
 249 static void SplitCurveByPoint(
 250   const NCollection_Vector<TopoDS_Edge>& theEdges,
 251   const Standard_Integer theEdgeIndex,
 252   const Standard_Real theParameter,
 253   NCollection_Vector<TopoDS_Edge>& theEdges1,
 254   NCollection_Vector<TopoDS_Edge>& theEdges2)
 255 {
 256   for (Standard_Integer aEI = 0; aEI < theEdgeIndex; ++aEI)
 257   {
 258     theEdges1.Append(theEdges(aEI));
 259   }
 260
 261   const TopoDS_Edge& aEdge = theEdges(theEdgeIndex);
 262   BRepAdaptor_Curve aCurve(aEdge);
 263   Standard_Integer aParamI = -1;
 264   const Standard_Real aEdgeEndParams[] =
 265     {aCurve.FirstParameter(), aCurve.LastParameter()};
 266   if (Abs(theParameter - aEdgeEndParams[0]) < Precision::PConfusion())
 267   {
 268     aParamI = 0;
 269   }
 270   else if (Abs(theParameter - aEdgeEndParams[1]) < Precision::PConfusion())
 271   {
 272     aParamI = 1;
 273   }
 274
 275   const TopAbs_Orientation aOrient = aEdge.Orientation();
 276   if (aOrient == TopAbs_REVERSED)
 277   {
 278     aParamI ^= 1;
 279   }
 280
 281   NCollection_Vector<TopoDS_Edge>* aEdges = &theEdges2;
 282   const Standard_Integer aECount = theEdges.Size();
 283   if (aParamI == 0)
 284   {
 285     aEdges = (theEdgeIndex == 0) ? &theEdges1 : &theEdges2;
 286     aEdges->Append(ReplaceVertex(aEdge, 0));
 287   }
 288   else if (aParamI == 1)
 289   {
 290     theEdges1.Append(ReplaceVertex(aEdge, 1));
 291   }
 292   else
 293   {
 294     TopoDS_Edge aFE = TopoDS::Edge(aEdge.Oriented(TopAbs_FORWARD));
 295     ShapeAnalysis_TransferParametersProj aSATPP(aFE, TopoDS_Face());
 296     aSATPP.SetMaxTolerance(Precision::Confusion());
 297     TopoDS_Vertex aSplitV1, aSplitV2;
 298     BRep_Builder().MakeVertex(
 299       aSplitV1, aCurve.Value(theParameter), Precision::Confusion());
 300     BRep_Builder().MakeVertex(
 301       aSplitV2, aCurve.Value(theParameter), Precision::Confusion());
 302     TopoDS_Edge aEParts[] = {
 303       ShapeBuild_Edge().CopyReplaceVertices(aFE, TopoDS_Vertex(),
 304         TopoDS::Vertex(aSplitV1.Oriented(TopAbs_REVERSED))),
 305       ShapeBuild_Edge().CopyReplaceVertices(aFE, aSplitV2, TopoDS_Vertex())};
 306     ShapeBuild_Edge().CopyPCurves(aEParts[0], aFE);
 307     ShapeBuild_Edge().CopyPCurves(aEParts[1], aFE);
 308     BRep_Builder().SameRange(aEParts[0], Standard_False);
 309     BRep_Builder().SameRange(aEParts[1], Standard_False);
 310     BRep_Builder().SameParameter(aEParts[0], Standard_False);
 311     BRep_Builder().SameParameter(aEParts[1], Standard_False);
 312     aSATPP.TransferRange(aEParts[0],
 313       aEdgeEndParams[0], theParameter, Standard_False);
 314     aSATPP.TransferRange(aEParts[1],
 315       theParameter, aEdgeEndParams[1], Standard_False);
 316     aEParts[0].Orientation(aOrient);
 317     aEParts[1].Orientation(aOrient);
 318
 319     const Standard_Integer aFirstPI = (aOrient != TopAbs_REVERSED) ? 0 : 1;
 320     theEdges1.Append(aEParts[aFirstPI]);
 321     theEdges2.Append(aEParts[1 - aFirstPI]);
 322   }
 323   for (Standard_Integer aEI = theEdgeIndex + 1; aEI < aECount; ++aEI)
 324   {
 325     aEdges->Append(theEdges(aEI));
 326   }
 327 }
 328
 329 static void SplitCurveByPoints(
 330   const NCollection_Vector<TopoDS_Edge>& theCurve,
 331   const NCollection_Vector<NCollection_Sequence<Standard_Real> >& theParameters,
 332   NCollection_Vector<NCollection_Vector<TopoDS_Edge> >& theSplittedCurves)
 333 {
 334   NCollection_Vector<TopoDS_Edge> aCurves[3];
 335   aCurves[0] = theCurve;
 336   Standard_Integer aCI = 0, aEI = 0;
 337   NCollection_Vector<TopoDS_Edge>::Iterator aEIt(theCurve);
 338   for (NCollection_Vector<NCollection_Sequence<Standard_Real> >::Iterator
 339     aPLIt(theParameters); aPLIt.More(); ++aEI, aEIt.Next(), aPLIt.Next())
 340   {
 341     const Standard_Boolean isForward =
 342       (aEIt.Value().Orientation() != TopAbs_REVERSED);
 343     for (NCollection_Sequence<Standard_Real>::Iterator
 344       aPIt(aPLIt.Value(), isForward); aPIt.More(); aPIt.Next())
 345     {
 346       const Standard_Integer aCI1 = (aCI + 1) % 3, aCI2 = (aCI + 2) % 3;
 347       SplitCurveByPoint(aCurves[aCI], aEI,
 348         aPIt.Value(), aCurves[aCI1], aCurves[aCI2]);
 349       if (!aCurves[aCI2].IsEmpty())
 350       {
 351         theSplittedCurves.Append(aCurves[aCI1]);
 352         aCurves[aCI].Clear();
 353         aCI = aCI2;
 354         aEI = 0;
 355       }
 356       aCurves[aCI1].Clear();
 357     }
 358   }
 359   theSplittedCurves.Append(aCurves[aCI]);
 360 }
 361
 362 static Standard_Integer AppendIntersectionPoint(
 363   const Adaptor3d_Curve& theCurve,
 364   const Standard_Real theParameter,
 365   NCollection_Sequence<Standard_Real>& theParameters)
 366 {
 367   // Check the coincidence.
 368   NCollection_Sequence<Standard_Real> aEndParams;
 369   aEndParams.Append(theCurve.FirstParameter());
 370   aEndParams.Append(theCurve.LastParameter());
 371   NCollection_Sequence<Standard_Real>* aParams[] =
 372     {&theParameters, &aEndParams};
 373   const gp_Pnt aPoint = theCurve.Value(theParameter);
 374   for (Standard_Integer aLI = 0; aLI < 2; ++aLI)
 375   {
 376     NCollection_Sequence<Standard_Real>::Iterator aPIt(*aParams[aLI]);
 377     for (Standard_Integer aPI = 0; aPIt.More(); aPIt.Next(), ++aPI)
 378     {
 379       const Standard_Real aParam = aPIt.Value();
 380       if (Abs(theParameter - aParam) < Precision::PConfusion() ||
 381         aPoint.SquareDistance(theCurve.Value(aParam)) <=
 382           Precision::SquareConfusion())
 383       {
 384         Standard_Integer aIntCount = 0;
 385         if (aLI != 0)
 386         {
 387           if (aPI == 0)
 388           {
 389             theParameters.Prepend(aEndParams.First());
 390           }
 391           else
 392           {
 393             theParameters.Append(aEndParams.Last());
 394           }
 395           ++aIntCount;
 396         }
 397         return aIntCount;
 398       }
 399     }
 400   }
 401
 402   // Calculate the position to insert.
 403   NCollection_Sequence<Standard_Real>::Iterator aPIt(theParameters);
 404   if (aPIt.More() && aPIt.Value() < theParameter)
 405   {
 406     NCollection_Sequence<Standard_Real>::Iterator aPIt2 = aPIt;
 407     aPIt2.Next();
 408     for (; aPIt2.More() && aPIt2.Value() < theParameter;
 409       aPIt.Next(), aPIt2.Next());
 410     theParameters.InsertAfter(aPIt, theParameter);
 411   }
 412   else
 413   {
 414     theParameters.Prepend(theParameter);
 415   }
 416   return 1;
 417 }
 418
 419 static Standard_Integer IntersectEdge(
 420   const TopoDS_Edge& theEdge1,
 421   const TopoDS_Edge& theEdge2,
 422   NCollection_Sequence<Standard_Real>& theParameters)
 423 {
 424   // Process the ends.
 425   Standard_Integer aIntCount = 0;
 426   BRepAdaptor_Curve aCurve1 = BRepAdaptor_Curve(theEdge1);
 427   BRepAdaptor_Curve aCurve2 = BRepAdaptor_Curve(theEdge2);
 428   const gp_Pnt aEndPs[] = {aCurve2.Value(aCurve2.FirstParameter()),
 429     aCurve2.Value(aCurve2.LastParameter())};
 430   for (Standard_Integer aPI = 0; aPI < 2; ++aPI)
 431   {
 432     Standard_Real aParameter;
 433     if (ProjectPointToCurve(aEndPs[aPI], aCurve1, aParameter) <=
 434       Precision::SquareConfusion())
 435     {
 436       AppendIntersectionPoint(aCurve1, aParameter, theParameters);
 437     }
 438   }
 439
 440   // Process the internal extrema.
 441   Extrema_ExtCC aAlgo(aCurve1, aCurve2);
 442   if (aAlgo.IsDone())
 443   {
 444     const Standard_Integer aECount = aAlgo.NbExt();
 445     for (Standard_Integer aEN = 1; aEN <= aECount; ++aEN)
 446     {
 447       Extrema_POnCurv aP1, aP2;
 448       aAlgo.Points(aEN, aP1, aP2);
 449       if (aP1.Value().SquareDistance(aP2.Value()) <=
 450         Precision::SquareConfusion())
 451       {
 452         AppendIntersectionPoint(aCurve1, aP1.Parameter(), theParameters);
 453       }
 454     }
 455   }
 456   return aIntCount;
 457 }
 458
 459 static Standard_Integer IntersectCurve(
 460   const NCollection_Vector<TopoDS_Edge>& theEdges,
 461   const TopoDS_Wire& theWire,
 462   NCollection_Vector<NCollection_Sequence<Standard_Real> >& theParameters)
 463 {
 464   Standard_Integer aIntCount = 0;
 465   const Standard_Integer aECount1 = theEdges.Size();
 466   for (Standard_Integer aEI1 = 0; aEI1 < aECount1; ++aEI1)
 467   {
 468     const TopoDS_Edge& aEdge1 = theEdges(aEI1);
 469     TopExp_Explorer aEIt2(theWire, TopAbs_EDGE);
 470     for (; aEIt2.More(); aEIt2.Next())
 471     {
 472       aIntCount += IntersectEdge(aEdge1,
 473         TopoDS::Edge(aEIt2.Current()), theParameters(aEI1));
 474     }
 475   }
 476   return aIntCount;
 477 }
 478
 479 static void CloseCurve(NCollection_Sequence<TopoDS_Edge>& theCurve)
 480 {
 481   const TopoDS_Vertex aVertex = TopoDS::Vertex(TopExp::LastVertex(
 482     theCurve.Last(), Standard_True).Oriented(TopAbs_FORWARD));
 483   const TopoDS_Edge& aEdge = theCurve.First();
 484   const TopoDS_Edge aForwardEdge = TopoDS::Edge(aEdge.Oriented(TopAbs_FORWARD));
 485   theCurve.ChangeFirst() = TopoDS::Edge(ShapeBuild_Edge().CopyReplaceVertices(
 486     aForwardEdge, aVertex, TopoDS_Vertex()).Oriented(aEdge.Orientation()));
 487 }
 488
 489 static Standard_Boolean IsClosed(
 490   const NCollection_Sequence<TopoDS_Edge>& theCurve)
 491 {
 492   return TopExp::FirstVertex(theCurve.First(), Standard_True).
 493     IsSame(TopExp::LastVertex(theCurve.Last(), Standard_True));
 494 }
 495
 496 static void ExtendCurve(
 497   const Standard_Integer thePosition,
 498   NCollection_Sequence<TopoDS_Edge>& theCurve,
 499   NCollection_Sequence<TopoDS_Edge>& theExtension)
 500 {
 501   if (thePosition == 0)
 502   {
 503     const TopoDS_Vertex aVertex = TopoDS::Vertex(TopExp::LastVertex(
 504       theExtension.Last(), Standard_True).Oriented(TopAbs_FORWARD));
 505     TopoDS_Edge& aEdge = theCurve.ChangeFirst();
 506     aEdge = TopoDS::Edge(ShapeBuild_Edge().CopyReplaceVertices(
 507         TopoDS::Edge(aEdge.Oriented(TopAbs_FORWARD)), aVertex, TopoDS_Vertex()).
 508       Oriented(aEdge.Orientation()));
 509     theCurve.Prepend(theExtension);
 510   }
 511   else
 512   {
 513     const TopoDS_Vertex aVertex = TopoDS::Vertex(TopExp::FirstVertex(
 514       theExtension.First(), Standard_True).Oriented(TopAbs_REVERSED));
 515     TopoDS_Edge& aEdge = theCurve.ChangeLast();
 516     aEdge = TopoDS::Edge(ShapeBuild_Edge().CopyReplaceVertices(
 517         TopoDS::Edge(aEdge.Oriented(TopAbs_FORWARD)), TopoDS_Vertex(), aVertex).
 518       Oriented(aEdge.Orientation()));
 519     theCurve.Append(theExtension);
 520   }
 521 }
 522
 523 static void ExtendCurve(
 524   const Standard_Integer thePosition,
 525   const TopoDS_Edge& theExtension,
 526   NCollection_Sequence<TopoDS_Edge>& theCurve)
 527 {
 528   NCollection_Sequence<TopoDS_Edge> aExtension;
 529   aExtension.Append(theExtension);
 530   ExtendCurve(thePosition, theCurve, aExtension);
 531 }
 532
 533 static gp_XYZ Tangent(
 534   const Adaptor3d_Curve& theCurve, const Standard_Integer thePosition)
 535 {
 536   const Standard_Real aParam = (thePosition == 0) ?
 537     theCurve.FirstParameter() : theCurve.LastParameter();
 538   gp_Pnt aP;
 539   gp_Vec aV;
 540   theCurve.D1(aParam, aP, aV);
 541   Standard_Real aNorm = aV.Magnitude();
 542   aNorm = (aNorm >= Precision::PConfusion()) ? aNorm : 0;
 543   return ((1 / aNorm) * aV).XYZ();
 544 }
 545
 546 static gp_XYZ Tangent(
 547   const TopoDS_Edge& theEdge, const Standard_Integer thePosition)
 548 {
 549   BRepAdaptor_Curve aCurve(theEdge);
 550   return Tangent(BRepAdaptor_Curve(theEdge), thePosition);
 551 }
 552
 553 static Standard_Boolean Interpolate(
 554   const gp_Pnt thePoint1,
 555   const gp_Pnt thePoint2,
 556   const gp_Vec theTangent1,
 557   const gp_Vec theTangent2,
 558   Handle(Geom_BSplineCurve)& theBSpline)
 559 {
 560   Handle(TColgp_HArray1OfPnt) aPs = new TColgp_HArray1OfPnt(1, 2);
 561   TColgp_Array1OfVec aTs(1, 2);
 562   Handle(TColStd_HArray1OfBoolean) aTFs = new TColStd_HArray1OfBoolean(1, 2);
 563   aPs->SetValue(1, thePoint1);
 564   aPs->SetValue(2, thePoint2);
 565   aTs.SetValue(1, theTangent1);
 566   aTs.SetValue(2, theTangent2);
 567   aTFs->SetValue(1, Standard_True);
 568   aTFs->SetValue(2, Standard_True);
 569   GeomAPI_Interpolate aInterpolator(aPs, Standard_False, 0);
 570   aInterpolator.Load(aTs, aTFs, Standard_False);
 571   aInterpolator.Perform();
 572   const Standard_Boolean aResult = (aInterpolator.IsDone() == Standard_True);
 573   if (aResult)
 574   {
 575     theBSpline = aInterpolator.Curve();
 576   }
 577   return aResult;
 578 }
 579
 580 static Standard_Boolean Merge(
 581   const TopoDS_Wire& theWire,
 582   const Standard_Real theTolerance,
 583   NCollection_Vector<NCollection_Sequence<TopoDS_Edge> >& theMergedCurves)
 584 {
 585   NCollection_Sequence<TopoDS_Edge> aCurve;
 586   Standard_Boolean isClosed;
 587   if (!WireToCurve(theWire, aCurve, isClosed))
 588   {
 589     return Standard_False;
 590   }
 591
 592   if (isClosed)
 593   {
 594     theMergedCurves.Append(aCurve);
 595     return Standard_True;
 596   }
 597
 598   const Standard_Real aSqTol = theTolerance * theTolerance;
 599   const gp_Pnt aPs[] = {
 600     BRep_Tool::Pnt(TopExp::FirstVertex(aCurve.First(), Standard_True)),
 601     BRep_Tool::Pnt(TopExp::LastVertex(aCurve.Last(), Standard_True))};
 602   if (!isClosed && aPs[0].SquareDistance(aPs[1]) <= aSqTol)
 603   {
 604     CloseCurve(aCurve);
 605     theMergedCurves.Append(aCurve);
 606     return Standard_True;
 607   }
 608
 609   NCollection_Sequence<TopoDS_Edge>* aCurves[] = {NULL, NULL};
 610   Standard_Integer aOrder = 0;
 611   for (NCollection_Vector<NCollection_Sequence<TopoDS_Edge> >::Iterator
 612     aCIt(theMergedCurves); aCIt.More(); aCIt.Next())
 613   {
 614     NCollection_Sequence<TopoDS_Edge>& aC = aCIt.ChangeValue();
 615     if (aC.IsEmpty() || IsClosed(aC))
 616     {
 617       continue;
 618     }
 619
 620     const gp_Pnt aP1 =
 621       BRep_Tool::Pnt(TopExp::FirstVertex(aC.First(), Standard_True));
 622     if (aCurves[0] == NULL && aP1.SquareDistance(aPs[1]) <= aSqTol)
 623     {
 624       aCurves[0] = &aC;
 625     }
 626
 627     const gp_Pnt aP2 =
 628       BRep_Tool::Pnt(TopExp::LastVertex(aC.Last(), Standard_True));
 629     if (aCurves[1] == NULL && aP2.SquareDistance(aPs[0]) <= aSqTol)
 630     {
 631       aCurves[1] = &aC;
 632       aOrder = (aCurves[0] == NULL) ? 1 : 0;
 633     }
 634   }
 635
 636   if (aCurves[0] == NULL && aCurves[1] == NULL)
 637   {
 638     theMergedCurves.Append(aCurve);
 639   }
 640   else if (aCurves[1] == NULL)
 641   {
 642     ExtendCurve(0, *aCurves[0], aCurve);
 643   }
 644   else if (aCurves[0] == NULL)
 645   {
 646     ExtendCurve(1, *aCurves[1], aCurve);
 647   }
 648   else
 649   {
 650     ExtendCurve(aOrder, *aCurves[aOrder], aCurve);
 651     if (aCurves[0] != aCurves[1])
 652     {
 653       ExtendCurve(aOrder, *aCurves[aOrder], *aCurves[1 - aOrder]);
 654     }
 655     else
 656     {
 657       CloseCurve(*aCurves[aOrder]);
 658     }
 659   }
 660   return Standard_True;
 661 }
 662
 663 static Standard_Boolean Connect(
 664   const TopoDS_Wire& theWire,
 665   const Standard_Real theTolerance,
 666   NCollection_Vector<NCollection_Sequence<TopoDS_Edge> >& theMergedCurves)
 667 {
 668   NCollection_Sequence<TopoDS_Edge> aCurve;
 669   Standard_Boolean isClosed;
 670   if (!WireToCurve(theWire, aCurve, isClosed))
 671   {
 672     return Standard_False;
 673   }
 674
 675   if (isClosed)
 676   {
 677     theMergedCurves.Append(aCurve);
 678     return Standard_True;
 679   }
 680
 681   const Standard_Real aSqTol = theTolerance * theTolerance;
 682   const gp_Pnt aPs[] = {
 683     BRep_Tool::Pnt(TopExp::FirstVertex(aCurve.First(), Standard_True)),
 684     BRep_Tool::Pnt(TopExp::LastVertex(aCurve.Last(), Standard_True))};
 685   if (!isClosed && aPs[0].SquareDistance(aPs[1]) <= aSqTol)
 686   {
 687     CloseCurve(aCurve);
 688     theMergedCurves.Append(aCurve);
 689     return Standard_True;
 690   }
 691
 692   for (NCollection_Vector<NCollection_Sequence<TopoDS_Edge> >::Iterator
 693     aCIt(theMergedCurves); aCIt.More(); aCIt.Next())
 694   {
 695     NCollection_Sequence<TopoDS_Edge>& aCurve2 = aCIt.ChangeValue();
 696     if (aCurve2.IsEmpty() || IsClosed(aCurve2))
 697     {
 698       continue;
 699     }
 700
 701     const TopoDS_Edge* aEdges2[] = {&aCurve2.First(), &aCurve2.Last()};
 702     const gp_Pnt aPs2[] = {
 703       BRep_Tool::Pnt(TopExp::FirstVertex(*aEdges2[0], Standard_True)),
 704       BRep_Tool::Pnt(TopExp::LastVertex(*aEdges2[1], Standard_True))};
 705     const Standard_Real aSqDists[] =
 706       {aPs2[0].SquareDistance(aPs[1]), aPs2[1].SquareDistance(aPs[0])};
 707     const Standard_Integer aOrder = (aSqDists[0] <= aSqDists[1]) ? 0 : 1;
 708     if (aSqDists[aOrder] > aSqTol)
 709     {
 710       const TopoDS_Edge& aEdge = (aOrder == 0) ? aCurve.Last() : aCurve.First();
 711       const gp_Pnt aPs3[] = {aPs[1 - aOrder], aPs2[aOrder]};
 712       const gp_XYZ aTs[] =
 713         {Tangent(aEdge, 1 - aOrder), Tangent(*aEdges2[aOrder], aOrder)};
 714       Handle(Geom_BSplineCurve) aBSpline;
 715       if (!Interpolate(aPs3[aOrder], aPs3[1 - aOrder],
 716         aTs[aOrder], aTs[1 - aOrder], aBSpline))
 717       {
 718         return Standard_False;
 719       }
 720
 721       ExtendCurve(aOrder, BRepBuilderAPI_MakeEdge(aBSpline), aCurve2);
 722     }
 723     ExtendCurve(aOrder, aCurve2, aCurve);
 724     if (aSqDists[1 - aOrder] <= aSqTol)
 725     {
 726       CloseCurve(aCurve2);
 727     }
 728     return Standard_True;
 729   }
 730
 731   theMergedCurves.Append(aCurve);
 732   return Standard_True;
 733 }
 734
 735 bool HYDROData_PolylineOperator::Split( const Handle( HYDROData_Document )& theDoc,
 736                                         const Handle( HYDROData_PolylineXY )& thePolyline,
 737                                         const gp_Pnt2d& thePoint,
 738                                         double theTolerance ) const
 739 {
 740   std::vector<gp_Pnt2d> aPointsList( 1 );
 741   aPointsList[0] = thePoint;
 742   std::vector<TopoDS_Wire> aCurves = GetWires( thePolyline );
 743   bool isOK = true;
 744   for( int i=0, n=aCurves.size(); i<n; i++ )
 745   {
 746     std::vector<TopoDS_Shape> aCurvesList = Split( aCurves[i], thePoint, theTolerance );
 747     bool isLocalOK = CreatePolylines( theDoc, thePolyline->GetName(), aCurvesList, true );
 748     isOK = isOK && isLocalOK;
 749   }
 750   return isOK;
 751 }
 752
 753 bool HYDROData_PolylineOperator::Split( const Handle( HYDROData_Document )& theDoc,
 754               const Handle( HYDROData_PolylineXY )& thePolyline,
 755               const Handle( HYDROData_PolylineXY )& theTool,
 756               double theTolerance ) const
 757 {
 758   HYDROData_SequenceOfObjects aSeq;
 759   aSeq.Append( theTool );
 760   return split( theDoc, thePolyline, aSeq, theTolerance, -1 );
 761 }
 762
 763 bool HYDROData_PolylineOperator::Split( const Handle( HYDROData_Document )& theDoc,
 764                                         const HYDROData_SequenceOfObjects& thePolylines,
 765                                         double theTolerance )
 766 {
 767   int f = thePolylines.Lower(), l = thePolylines.Upper();
 768   for( int i=f; i<=l; i++ )
 769   {
 770     Handle( HYDROData_PolylineXY ) aPolyline = Handle( HYDROData_PolylineXY )::DownCast( thePolylines.Value( i ) );
 771     if( !split( theDoc, aPolyline, thePolylines, theTolerance, i ) )
 772       return false;
 773   }
 774   return true;
 775 }
 776
 777 bool HYDROData_PolylineOperator::Merge( const Handle( HYDROData_Document )& theDoc,
 778                                         const QString& theName,
 779                                         const HYDROData_SequenceOfObjects& thePolylines,
 780                                         bool isConnectByNewSegment,
 781                                         double theTolerance )
 782 {
 783   NCollection_Vector<NCollection_Sequence<TopoDS_Edge> > aMergedCurves;
 784   HYDROData_SequenceOfObjects::Iterator aPIt(thePolylines);
 785   for (; aPIt.More(); aPIt.Next())
 786   {
 787     Handle(HYDROData_PolylineXY) aPolyline =
 788       Handle(HYDROData_PolylineXY)::DownCast(aPIt.Value());
 789     std::vector<TopoDS_Wire> aWires = GetWires(aPolyline);
 790     for (std::vector<TopoDS_Wire>::const_iterator aWIt = aWires.begin(),
 791       aLastWIt = aWires.end(); aWIt != aLastWIt; ++aWIt)
 792     {
 793       const Standard_Boolean aResult = !isConnectByNewSegment ?
 794         ::Merge(*aWIt, theTolerance, aMergedCurves) :
 795         Connect(*aWIt, theTolerance, aMergedCurves);
 796       if (!aResult)
 797       {
 798         return false;
 799       }
 800     }
 801   }
 802
 803   TopoDS_Compound aWireSet;
 804   BRep_Builder aBuilder;
 805   aBuilder.MakeCompound(aWireSet);
 806   for (NCollection_Vector<NCollection_Sequence<TopoDS_Edge> >::Iterator
 807     aCIt(aMergedCurves); aCIt.More(); aCIt.Next())
 808   {
 809     if (!aCIt.Value().IsEmpty())
 810     {
 811       TopoDS_Wire aWire;
 812       CurveToWire(aCIt.Value(), aWire);
 813       aBuilder.Add(aWireSet, aWire);
 814     }
 815   }
 816
 817   std::vector<TopoDS_Shape> aPolylines(1);
 818   aPolylines[0] = aWireSet;
 819   CreatePolylines(theDoc, theName, aPolylines, false);
 820   return true;
 821 }
 822
 823 bool HYDROData_PolylineOperator::split( const Handle( HYDROData_Document )& theDoc,
 824                                         const Handle( HYDROData_PolylineXY )& thePolyline,
 825                                         const HYDROData_SequenceOfObjects& theTools,
 826                                         double theTolerance,
 827                                         int theIgnoreIndex ) const
 828 {
 829   std::vector<TopoDS_Wire> aCurves = GetWires( thePolyline );
 830   std::vector<TopoDS_Wire> aToolCurves;
 831   for( int i=theTools.Lower(), n=theTools.Upper(); i<=n; i++ )
 832     if( i!=theIgnoreIndex )
 833     {
 834       Handle( HYDROData_PolylineXY ) aToolPolyline =
 835         Handle( HYDROData_PolylineXY )::DownCast( theTools.Value( i ) );
 836       append( aToolCurves, GetWires( aToolPolyline ) );
 837     }
 838
 839   const int aPSCount = aCurves.size();
 840   const int aTSCount = aToolCurves.size();
 841   std::vector<TopoDS_Shape> aResult;
 842   for (int aPSI = 0; aPSI < aPSCount; ++aPSI)
 843   {
 844     NCollection_Vector<TopoDS_Edge> aCurve;
 845     Standard_Boolean aIsClosed;
 846     if (!WireToCurve(aCurves[aPSI], aCurve, aIsClosed))
 847     {
 848       continue;
 849     }
 850
 851     NCollection_Vector<NCollection_Sequence<Standard_Real> > aParams;
 852     aParams.SetValue(aCurve.Size() - 1, NCollection_Sequence<Standard_Real>());
 853     for (int aTSI = 0; aTSI < aTSCount; ++aTSI)
 854     {
 855       IntersectCurve(aCurve, aToolCurves[aTSI], aParams);
 856     }
 857
 858     NCollection_Vector<NCollection_Vector<TopoDS_Edge> > aSplittedCurves;
 859     SplitCurveByPoints(aCurve, aParams, aSplittedCurves);
 860
 861     Standard_Boolean aIsClosed2 = aIsClosed;
 862     if (aIsClosed2)
 863     {
 864       const NCollection_Sequence<Standard_Real>& aPs = aParams.First();
 865       if (!aPs.IsEmpty())
 866       {
 867         const TopoDS_Edge& aEdge = aCurve.First();
 868         const Standard_Boolean isForward =
 869           (aEdge.Orientation() != TopAbs_REVERSED);
 870         const Standard_Real aParam = isForward ? aPs.First() : aPs.Last();
 871         BRepAdaptor_Curve aCurve(aEdge);
 872         const Standard_Real aEndParam = isForward ?
 873           aCurve.FirstParameter() : aCurve.LastParameter();
 874         aIsClosed2 = (Abs(aParam - aEndParam) > Precision::PConfusion());
 875       }
 876
 877       if (aIsClosed2)
 878       {
 879         const NCollection_Sequence<Standard_Real>& aPs = aParams.Last();
 880         if (!aPs.IsEmpty())
 881         {
 882           const TopoDS_Edge& aEdge = aCurve.Last();
 883           const Standard_Boolean isForward =
 884             (aEdge.Orientation() != TopAbs_REVERSED);
 885           const Standard_Real aParam = isForward ? aPs.Last() : aPs.First();
 886           BRepAdaptor_Curve aCurve(aEdge);
 887           const Standard_Real aEndParam = isForward ?
 888             aCurve.LastParameter() : aCurve.FirstParameter();
 889           aIsClosed2 = (Abs(aParam - aEndParam) >= Precision::PConfusion());
 890         }
 891       }
 892     }
 893
 894     Standard_Integer aFSCI = 0, aLSCI = aSplittedCurves.Size() - 1;
 895     if (aIsClosed2 && aFSCI < aLSCI)
 896     {
 897       TopoDS_Wire aWire;
 898       CurveToWire(aSplittedCurves(aLSCI), aSplittedCurves(aFSCI), aWire);
 899       aResult.push_back(aWire);
 900       ++aFSCI;
 901       --aLSCI;
 902     }
 903
 904     for (Standard_Integer aSCI = aFSCI; aSCI <= aLSCI; ++aSCI)
 905     {
 906       TopoDS_Wire aWire;
 907       CurveToWire(aSplittedCurves(aSCI), aWire);
 908       aResult.push_back(aWire);
 909     }
 910   }
 911
 912   CreatePolylines(theDoc, thePolyline->GetName(), aResult, true);
 913   return true;
 914 }
 915
 916 std::vector<TopoDS_Wire> HYDROData_PolylineOperator::GetWires( const Handle( HYDROData_PolylineXY )& thePolyline )
 917 {
 918   std::vector<TopoDS_Wire> aResult;
 919
 920   TopoDS_Shape aShape = thePolyline->GetShape();
 921
 922   if( aShape.ShapeType()==TopAbs_WIRE )
 923   {
 924     aResult.push_back( TopoDS::Wire( aShape ) );
 925   }
 926   else
 927   {
 928     TopExp_Explorer anExp( aShape, TopAbs_WIRE );
 929     for( ; anExp.More(); anExp.Next() )
 930     {
 931       aResult.push_back( TopoDS::Wire( anExp.Current() ) );
 932     }
 933   }
 934   return aResult;
 935 }
 936
 937 std::vector<TopoDS_Shape> HYDROData_PolylineOperator::Split( const TopoDS_Wire& theWire,
 938                                                              const gp_Pnt2d& thePoint,
 939                                                              double theTolerance )
 940 {
 941   std::vector<TopoDS_Shape> aResult;
 942   NCollection_Vector<TopoDS_Edge> aEdges;
 943   Standard_Boolean isClosed;
 944   if (!WireToCurve(theWire, aEdges, isClosed))
 945   {
 946     aResult.push_back(theWire);
 947     return aResult;
 948   }
 949
 950   const gp_Pnt aP(thePoint.X(), thePoint.Y(), 0);
 951   Standard_Real aMinSqDist = DBL_MAX;
 952   int aSEI = -1;
 953   Standard_Real aSParam;
 954   for (int aECount = aEdges.Size(), aEI = 0; aEI < aECount; ++aEI)
 955   {
 956     Standard_Real aParam;
 957     const Standard_Real aSqDist =
 958       ProjectPointToEdge(aP, aEdges(aEI), aParam);
 959     if (aMinSqDist > aSqDist)
 960     {
 961       aMinSqDist = aSqDist;
 962       aSEI = aEI;
 963       aSParam = aParam;
 964     }
 965   }
 966
 967   NCollection_Vector<TopoDS_Edge> aEdges1, aEdges2;
 968   SplitCurveByPoint(aEdges, aSEI, aSParam, aEdges1, aEdges2);
 969   TopoDS_Wire aWire;
 970   if (!isClosed)
 971   {
 972     CurveToWire(aEdges1, aWire);
 973     if (!aEdges2.IsEmpty())
 974     {
 975       aResult.push_back(aWire);
 976       CurveToWire(aEdges2, aWire);
 977     }
 978   }
 979   else
 980   {
 981     if (!aEdges2.IsEmpty())
 982     {
 983       CurveToWire(aEdges2, aEdges1, aWire);
 984     }
 985     else
 986     {
 987       CurveToWire(aEdges1, aWire);
 988     }
 989   }
 990   aResult.push_back(aWire);
 991   return aResult;
 992 }
 993
 994 bool HYDROData_PolylineOperator::CreatePolylines( const Handle( HYDROData_Document )& theDoc,
 995                                                   const QString& theNamePrefix,
 996                                                   const std::vector<TopoDS_Shape>& theShapes,
 997                                                   bool isUseIndices )
 998 {
 999   if( theDoc.IsNull() )
1000     return false;
1001
1002   int n = theShapes.size();
1003   int anIndex = 1;
1004   for( int i=0; i<n; i++ )
1005   {
1006     Handle( HYDROData_PolylineXY ) aPolyline =
1007       Handle( HYDROData_PolylineXY )::DownCast( theDoc->CreateObject( KIND_POLYLINEXY ) );
1008     if( aPolyline.IsNull() )
1009       return false;
1010
1011     aPolyline->SetShape( theShapes[i] );
1012
1013     if( isUseIndices )
1014     {
1015       QString aNewName = theNamePrefix + "_" + QString::number( anIndex );
1016       if( theDoc->FindObjectByName( aNewName ).IsNull() )  // the object with such a name is not found
1017         aPolyline->SetName( aNewName );
1018       anIndex++;
1019     }
1020     else
1021     {
1022       aPolyline->SetName( theNamePrefix );
1023     }
1024   }
1025   return true;
1026 }