src/GeomAlgoAPI/GeomAlgoAPI_CurveBuilder.cpp

   1 // Copyright (C) 2014-2024  CEA, EDF
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #include "GeomAlgoAPI_CurveBuilder.h"
  21
  22 #include <GeomAPI_Edge.h>
  23 #include <GeomAPI_Pnt.h>
  24 #include <GeomAPI_Vertex.h>
  25 #include <GeomAPI_ShapeExplorer.h>
  26
  27 #include <Approx_ParametrizationType.hxx>
  28 #include <BRepBuilderAPI_MakeEdge.hxx>
  29 #include <BSplSLib.hxx>
  30 #include <Geom_BSplineCurve.hxx>
  31 #include <Geom_BSplineSurface.hxx>
  32 #include <GeomAPI_Interpolate.hxx>
  33 #include <GeomAPI_PointsToBSplineSurface.hxx>
  34 #include <gp_Pnt.hxx>
  35 #include <Precision.hxx>
  36 #include <TColgp_Array2OfPnt.hxx>
  37 #include <TColgp_HArray1OfPnt.hxx>
  38 #include <TopExp_Explorer.hxx>
  39 #include <TopoDS.hxx>
  40 #include <TopoDS_Edge.hxx>
  41
  42
  43 GeomEdgePtr GeomAlgoAPI_CurveBuilder::edge(const std::list<GeomPointPtr>& thePoints,
  44                                            const bool thePeriodic,
  45                                            const bool theIsToReorder,
  46                                            const GeomDirPtr& theStartTangent,
  47                                            const GeomDirPtr& theEndTangent)
  48 {
  49   std::list<GeomPointPtr> aPointsCopy = thePoints;
  50
  51   // If the curve to be closed - remove last point if it is too close to the first one
  52   bool isClose = aPointsCopy.front()->distance(aPointsCopy.back()) <= gp::Resolution();
  53   if (isClose && thePeriodic) {
  54     aPointsCopy.pop_back();
  55   }
  56
  57   // Reorder points if required
  58   if (theIsToReorder) {
  59     reorderPoints(aPointsCopy);
  60   }
  61
  62   // Prepare points array
  63   Handle(TColgp_HArray1OfPnt) aPoints = new TColgp_HArray1OfPnt(1, (int)aPointsCopy.size());
  64   std::list<GeomPointPtr>::const_iterator anIt = aPointsCopy.begin();
  65   for (int i = 1; anIt != aPointsCopy.end(); anIt++, i++) {
  66     GeomPointPtr aPoint = *anIt;
  67     aPoints->SetValue(i, aPoint->impl<gp_Pnt>());
  68   }
  69
  70   // Initialize interpolator
  71   GeomAPI_Interpolate anInterp(aPoints, thePeriodic, gp::Resolution());
  72
  73   // Assign tangents if defined
  74   if (theStartTangent && theEndTangent) {
  75     gp_Dir aDir = theStartTangent->impl<gp_Dir>();
  76     gp_Vec anInitialTangent(aDir.XYZ());
  77     aDir = theEndTangent->impl<gp_Dir>();
  78     gp_Vec aFinalTangent(aDir.XYZ());
  79
  80     anInterp.Load(anInitialTangent, aFinalTangent);
  81   }
  82
  83   // Compute
  84   if (aPoints->Length() > 1) {
  85     anInterp.Perform();
  86   }
  87
  88   // Set result in form of edge
  89   TopoDS_Edge anEdge;
  90   if (anInterp.IsDone()) {
  91     anEdge = BRepBuilderAPI_MakeEdge(anInterp.Curve()).Edge();
  92   }
  93
  94   GeomEdgePtr aResultShape(new GeomAPI_Edge);
  95   aResultShape->setImpl(new TopoDS_Shape(anEdge));
  96
  97   return aResultShape;
  98 }
  99
 100 GeomEdgePtr GeomAlgoAPI_CurveBuilder::approximate(const std::list<GeomPointPtr>& thePoints,
 101                                                   const bool thePeriodic,
 102                                                   const double thePrecision)
 103 {
 104   // Prepare points array to be able to build a surface.
 105   // This surface is based on two sets of points: the first is an original and
 106   // the second is shifted along orthogonal direction.
 107   // This is a workaround, because GeomAPI_PointsToBSpline algorithm cannot produce
 108   // the periodic curve, but GeomAPI_PointsToBSplineSurface can.
 109   TColgp_Array2OfPnt aPoints(1, (int)thePoints.size(), 1, 2);
 110   gp_Pnt aPlaneBase[3]; // base points to calculate the normal direction
 111   int aNbPlanePoints = 0;
 112   gp_Dir aNormal;
 113   std::list<GeomPointPtr>::const_iterator anIt = thePoints.begin();
 114   for (int i = 1; anIt != thePoints.end(); anIt++, i++) {
 115     const gp_Pnt& aPoint = (*anIt)->impl<gp_Pnt>();
 116     aPoints.SetValue(i, 1, aPoint);
 117     aPoints.SetValue(i, 2, aPoint);
 118     if (aNbPlanePoints < 3) {
 119       if (aNbPlanePoints == 0 ||
 120           aPoint.SquareDistance(aPlaneBase[0]) > Precision::SquareConfusion())
 121         aPlaneBase[aNbPlanePoints++] = aPoint;
 122       if (aNbPlanePoints == 3) {
 123         gp_Vec aVec12(aPlaneBase[0], aPlaneBase[1]);
 124         gp_Vec aVec13(aPlaneBase[0], aPlaneBase[2]);
 125         if (aVec12.CrossSquareMagnitude(aVec13) > Precision::SquareConfusion())
 126           aNormal = gp_Dir(aVec12 ^ aVec13);
 127         else
 128           --aNbPlanePoints;
 129       }
 130     }
 131   }
 132   if (aNbPlanePoints < 3)
 133     aNormal = gp::DZ();
 134   // shifted points
 135   for (int i = aPoints.LowerRow(); i <= aPoints.UpperRow(); i++)
 136     aPoints.ChangeValue(i, 2).ChangeCoord() += aNormal.XYZ();
 137
 138   // If the curve to be closed - remove last point if it is too close to the first one
 139   bool isClose = aPoints.Value(aPoints.LowerRow(), 1).Distance(
 140                  aPoints.Value(aPoints.UpperRow(), 1)) <= gp::Resolution();
 141   if (isClose && thePeriodic) {
 142     aPoints.Resize(aPoints.LowerRow(), aPoints.UpperRow() - 1,
 143                    aPoints.LowerCol(), aPoints.UpperCol(), Standard_True);
 144   }
 145
 146   // Initialize and perform approximator
 147   static const Standard_Integer DEGREE_MIN = 3;
 148   static const Standard_Integer DEGREE_MAX = 8;
 149   GeomAPI_PointsToBSplineSurface anApprox;
 150   anApprox.Init(aPoints, Approx_ChordLength, DEGREE_MIN, DEGREE_MAX,
 151                 GeomAbs_C2, thePrecision, thePeriodic);
 152
 153   // Set result in form of edge
 154   TopoDS_Edge anEdge;
 155   if (anApprox.IsDone()) {
 156     // build a curve along U-direction of the surface
 157     Handle(Geom_BSplineSurface) aSurface = anApprox.Surface();
 158     Handle(Geom_Curve) aCurve = aSurface->VIso(aSurface->VKnots().First());
 159
 160     anEdge = BRepBuilderAPI_MakeEdge(aCurve).Edge();
 161   }
 162
 163   GeomEdgePtr aResultShape(new GeomAPI_Edge);
 164   aResultShape->setImpl(new TopoDS_Shape(anEdge));
 165
 166   return aResultShape;
 167 }
 168
 169 void GeomAlgoAPI_CurveBuilder::reorderPoints(std::list<GeomPointPtr>& thePoints)
 170 {
 171   if (thePoints.size() < 3) {
 172     return;
 173   }
 174
 175   std::list<GeomPointPtr>::iterator aPIt = thePoints.begin();
 176   GeomPointPtr aPrevPnt = *aPIt;
 177   for (; aPIt != thePoints.end(); ++aPIt) {
 178     GeomPointPtr aPnt = *aPIt;
 179     std::list<GeomPointPtr>::iterator aNextIt = aPIt;
 180     std::list<GeomPointPtr>::iterator aNearestIt = ++aNextIt;
 181     double aMinDist = RealLast();
 182     while (aNextIt != thePoints.end()) {
 183       double aDist = aPnt->distance(*aNextIt);
 184       if (aDist < Precision::Confusion()) {
 185         // remove duplicates
 186         std::list<GeomPointPtr>::iterator aRemoveIt = aNextIt++;
 187         thePoints.erase(aRemoveIt);
 188         // update iterator showing the nearest point, because it may become invalid
 189         aNearestIt = aPIt; ++aNearestIt;
 190         continue;
 191       }
 192       if (aDist < aMinDist && (aMinDist - aDist) > Precision::Confusion()) {
 193         aNearestIt = aNextIt;
 194         aMinDist = aDist;
 195       }
 196       ++aNextIt;
 197     }
 198     aNextIt = aPIt; ++aNextIt;
 199     if (aNearestIt != aNextIt) {
 200       // Keep given order of points to use it in case of equidistant candidates
 201       //              .--<---<--.
 202       //             |           |
 203       // o  o  o  c  o->o->o->o->n  o  o
 204       //          |  |           |
 205       //          i i+1       nearest
 206       GeomPointPtr aPointToMove = *aNearestIt;
 207       thePoints.erase(aNearestIt);
 208       thePoints.insert(aNextIt, aPointToMove);
 209     }
 210   }
 211 }