]> SALOME platform Git repositories - modules/shaper.git/blob - src/GeomAlgoAPI/GeomAlgoAPI_CurveBuilder.cpp
Salome HOME
856d5d31ceb74f79d853e49a4a637089c9c50b1c
[modules/shaper.git] / 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 }