ByoungSooPark/cadquery
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #694 from CadQuery/splineApprox

Browse Source 
Implement makeSplineApprox for edges and faces
This commit is contained in:
Jeremy Wright
2021-03-31 15:06:41 -04:00
committed by GitHub
parent df550f071a 7a259bed61
commit 7546866159
5 changed files with 536 additions and 38 deletions
Download Patch File Download Diff File Expand all files Collapse all files

225
cadquery/cq.py
View File

@ -35,7 +35,7 @@ from typing import (
Dict,
)
from typing_extensions import Literal
from inspect import Parameter, Signature
from .occ_impl.geom import Vector, Plane, Location
from .occ_impl.shapes import (
@ -264,7 +264,15 @@ class Workplane(object):
return list(all.values())
@overload
def split(self: T, keepTop: bool = False, keepBottom: bool = False) -> T:
...
@overload
def split(self: T, splitter: Union[T, Shape]) -> T:
...
def split(self: T, *args, **kwargs) -> T:
"""
Splits a solid on the stack into two parts, optionally keeping the separate parts.
@ -284,28 +292,64 @@ class Workplane(object):
c = c.faces(">Y").workplane(-0.5).split(keepTop=True)
"""
if (not keepTop) and (not keepBottom):
raise ValueError("You have to keep at least one half")
# split using an object
if len(args) == 1 and isinstance(args[0], (Workplane, Shape)):
solid = self.findSolid()
arg = args[0]
maxDim = solid.BoundingBox().DiagonalLength * 10.0
topCutBox = self.rect(maxDim, maxDim)._extrude(maxDim)
bottomCutBox = self.rect(maxDim, maxDim)._extrude(-maxDim)
solid = self.findSolid()
tools = (
(arg,)
if isinstance(arg, Shape)
else [v for v in arg.vals() if isinstance(v, Shape)]
)
rv = [solid.split(*tools)]
top = solid.cut(bottomCutBox)
bottom = solid.cut(topCutBox)
if keepTop and keepBottom:
# Put both on the stack, leave original unchanged.
return self.newObject([top, bottom])
# split using the current wokrplane
else:
# Put the one we are keeping on the stack, and also update the
# context solidto the one we kept.
if keepTop:
return self.newObject([top])
# boilerplate for arg/kwarg parsing
sig = Signature(
(
Parameter(
"keepTop", Parameter.POSITIONAL_OR_KEYWORD, default=False
),
Parameter(
"keepBottom", Parameter.POSITIONAL_OR_KEYWORD, default=False
),
)
)
bound_args = sig.bind(*args, **kwargs)
bound_args.apply_defaults()
keepTop = bound_args.arguments["keepTop"]
keepBottom = bound_args.arguments["keepBottom"]
if (not keepTop) and (not keepBottom):
raise ValueError("You have to keep at least one half")
solid = self.findSolid()
maxDim = solid.BoundingBox().DiagonalLength * 10.0
topCutBox = self.rect(maxDim, maxDim)._extrude(maxDim)
bottomCutBox = self.rect(maxDim, maxDim)._extrude(-maxDim)
top = solid.cut(bottomCutBox)
bottom = solid.cut(topCutBox)
if keepTop and keepBottom:
# Put both on the stack, leave original unchanged.
rv = [top, bottom]
else:
return self.newObject([bottom])
# Put the one we are keeping on the stack, and also update the
# context solid to the one we kept.
if keepTop:
rv = [top]
else:
rv = [bottom]
return self.newObject(rv)
@deprecate()
def combineSolids(
@ -1729,6 +1773,20 @@ class Workplane(object):
return self.eachpoint(lambda loc: slot.moved(loc), True)
def _toVectors(
self, pts: Iterable[VectorLike], includeCurrent: bool
) -> List[Vector]:
vecs = [self.plane.toWorldCoords(p) for p in pts]
if includeCurrent:
gstartPoint = self._findFromPoint(False)
allPoints = [gstartPoint] + vecs
else:
allPoints = vecs
return allPoints
def spline(
self: T,
listOfXYTuple: Iterable[VectorLike],
@ -1742,10 +1800,9 @@ class Workplane(object):
makeWire: bool = False,
) -> T:
"""
Create a spline interpolated through the provided points.
Create a spline interpolated through the provided points (2D or 3D).
:param listOfXYTuple: points to interpolate through
:type listOfXYTuple: list of 2-tuple
:param tangents: vectors specifying the direction of the tangent to the
curve at each of the specified interpolation points.
@ -1807,13 +1864,7 @@ class Workplane(object):
that cannot be correctly interpreted as a spline.
"""
vecs = [self.plane.toWorldCoords(p) for p in listOfXYTuple]
if includeCurrent:
gstartPoint = self._findFromPoint(False)
allPoints = [gstartPoint] + vecs
else:
allPoints = vecs
allPoints = self._toVectors(listOfXYTuple, includeCurrent)
if tangents:
tangents_g: Optional[Sequence[Vector]] = [
@ -1844,31 +1895,141 @@ class Workplane(object):
return self.newObject([rv_w if makeWire else e])
def splineApprox(
self: T,
points: Iterable[VectorLike],
tol: Optional[float] = 1e-6,
minDeg: int = 1,
maxDeg: int = 6,
smoothing: Optional[Tuple[float, float, float]] = (1, 1, 1),
forConstruction: bool = False,
includeCurrent: bool = False,
makeWire: bool = False,
) -> T:
"""
Create a spline interpolated through the provided points (2D or 3D).
:param points: points to interpolate through
:param tol: tolerance of the algorithm (default: 1e-6)
:param minDeg: minimum spline degree (default: 1)
:param maxDeg: maximum spline degree (default: 6)
:param smoothing: optional parameters for the variational smoothing algorithm (default: (1,1,1))
:param includeCurrent: use current point as a starting point of the curve
:param makeWire: convert the resulting spline edge to a wire
:return: a Workplane object with the current point at the end of the spline
*WARNING* for advanced users.
"""
allPoints = self._toVectors(points, includeCurrent)
e = Edge.makeSplineApprox(
allPoints,
minDeg=minDeg,
maxDeg=maxDeg,
smoothing=smoothing,
**({"tol": tol} if tol else {}),
)
if makeWire:
rv_w = Wire.assembleEdges([e])
if not forConstruction:
self._addPendingWire(rv_w)
else:
if not forConstruction:
self._addPendingEdge(e)
return self.newObject([rv_w if makeWire else e])
def parametricCurve(
self: T,
func: Callable[[float], VectorLike],
N: int = 400,
start: float = 0,
stop: float = 1,
tol: float = 1e-6,
minDeg: int = 1,
maxDeg: int = 6,
smoothing: Optional[Tuple[float, float, float]] = (1, 1, 1),
makeWire: bool = True,
) -> T:
"""
Create a spline interpolated through the provided points.
Create a spline curve approximating the provided function.
:param func: function f(t) that will generate (x,y) pairs
:type func: float --> (float,float)
:param func: function f(t) that will generate (x,y,z) pairs
:type func: float --> (float,float,float)
:param N: number of points for discretization
:param start: starting value of the parameter t
:param stop: final value of the parameter t
:param tol: tolerance of the algorithm (default: 1e-3)
:param minDeg: minimum spline degree (default: 1)
:param maxDeg: maximum spline degree (default: 6)
:param smoothing: optional parameters for the variational smoothing algorithm (default: (1,1,1))
:param makeWire: convert the resulting spline edge to a wire
:return: a Workplane object with the current point unchanged
"""
diff = stop - start
allPoints = [func(start + diff * t / N) for t in range(N + 1)]
allPoints = self._toVectors(
(func(start + diff * t / N) for t in range(N + 1)), False
)
return self.spline(allPoints, includeCurrent=False, makeWire=makeWire)
e = Edge.makeSplineApprox(
allPoints, tol=tol, smoothing=smoothing, minDeg=minDeg, maxDeg=maxDeg
)
if makeWire:
rv_w = Wire.assembleEdges([e])
self._addPendingWire(rv_w)
else:
self._addPendingEdge(e)
return self.newObject([rv_w if makeWire else e])
def parametricSurface(
self: T,
func: Callable[[float, float], VectorLike],
N: int = 20,
start: float = 0,
stop: float = 1,
tol: float = 1e-2,
minDeg: int = 1,
maxDeg: int = 6,
smoothing: Optional[Tuple[float, float, float]] = (1, 1, 1),
) -> T:
"""
Create a spline surface approximating the provided function.
:param func: function f(u,v) that will generate (x,y,z) pairs
:type func: (float,float) --> (float,float,float)
:param N: number of points for discretization in one direction
:param start: starting value of the parameters u,v
:param stop: final value of the parameters u,v
:param tol: tolerance used by the approximation algorithm (default: 1e-3)
:param minDeg: minimum spline degree (default: 1)
:param maxDeg: maximum spline degree (default: 3)
:param smoothing: optional parameters for the variational smoothing algorithm (default: (1,1,1))
:return: a Workplane object with the current point unchanged
This method might be unstable and may require tuning of the tol parameter.
"""
diff = stop - start
allPoints = []
for i in range(N + 1):
generator = (
func(start + diff * i / N, start + diff * j / N) for j in range(N + 1)
)
allPoints.append(self._toVectors(generator, False))
f = Face.makeSplineApprox(
allPoints, tol=tol, smoothing=smoothing, minDeg=minDeg, maxDeg=maxDeg
)
return self.newObject([f])
def ellipseArc(
self: T,

130
cadquery/occ_impl/shapes.py
View File

@ -20,6 +20,8 @@ from .geom import Vector, BoundBox, Plane, Location, Matrix
import OCP.TopAbs as ta # Tolopolgy type enum
import OCP.GeomAbs as ga # Geometry type enum
from OCP.Precision import Precision
from OCP.gp import (
gp_Vec,
gp_Pnt,
@ -36,7 +38,7 @@ from OCP.gp import (
)
# Array of points (used for B-spline construction):
from OCP.TColgp import TColgp_HArray1OfPnt
from OCP.TColgp import TColgp_HArray1OfPnt, TColgp_HArray2OfPnt
# Array of vectors (used for B-spline interpolation):
from OCP.TColgp import TColgp_Array1OfVec
@ -109,7 +111,12 @@ from OCP.TopoDS import (
from OCP.GC import GC_MakeArcOfCircle, GC_MakeArcOfEllipse # geometry construction
from OCP.GCE2d import GCE2d_MakeSegment
from OCP.GeomAPI import GeomAPI_Interpolate, GeomAPI_ProjectPointOnSurf
from OCP.GeomAPI import (
GeomAPI_Interpolate,
GeomAPI_ProjectPointOnSurf,
GeomAPI_PointsToBSpline,
GeomAPI_PointsToBSplineSurface,
)
from OCP.BRepFill import BRepFill
@ -118,6 +125,8 @@ from OCP.BRepAlgoAPI import (
BRepAlgoAPI_Fuse,
BRepAlgoAPI_Cut,
BRepAlgoAPI_BooleanOperation,
BRepAlgoAPI_Splitter,
BRepAlgoAPI_BuilderAlgo,
)
from OCP.Geom import (
@ -938,7 +947,7 @@ class Shape(object):
self,
args: Iterable["Shape"],
tools: Iterable["Shape"],
op: BRepAlgoAPI_BooleanOperation,
op: Union[BRepAlgoAPI_BooleanOperation, BRepAlgoAPI_Splitter],
) -> "Shape":
"""
Generic boolean operation
@ -999,6 +1008,15 @@ class Shape(object):
return self._bool_op((self,), toIntersect, intersect_op)
def split(self, *splitters: "Shape") -> "Shape":
"""
Split this shape with the positional arguments.
"""
split_op = BRepAlgoAPI_Splitter()
return self._bool_op((self,), splitters, split_op)
def mesh(self, tolerance: float, angularTolerance: float = 0.1):
"""
Generate traingulation if none exists.
@ -1379,6 +1397,19 @@ class Edge(Shape, Mixin1D):
return curve, BRepAdaptor_HCurve(curve)
def close(self) -> Union["Edge", "Wire"]:
"""
Close an Edge
"""
rv: Union[Wire, Edge]
if not self.IsClosed():
rv = Wire.assembleEdges((self,)).close()
else:
rv = self
return rv
@classmethod
def makeCircle(
cls: Type["Edge"],
@ -1537,6 +1568,45 @@ class Edge(Shape, Mixin1D):
return cls(BRepBuilderAPI_MakeEdge(spline_geom).Edge())
@classmethod
def makeSplineApprox(
cls: Type["Edge"],
listOfVector: List[Vector],
tol: float = 1e-3,
smoothing: Optional[Tuple[float, float, float]] = None,
minDeg: int = 1,
maxDeg: int = 6,
) -> "Edge":
"""
Approximate a spline through the provided points.
:param listOfVector: a list of Vectors that represent the points
:param tol: tolerance of the algorithm (consult OCC documentation).
:param smoothing: optional tuple of 3 weights use for variational smoothing (default: None)
:param minDeg: minimum spline degree. Enforced only when smothing is None (default: 1)
:param maxDeg: maximum spline degree (default: 6)
:return: an Edge
"""
pnts = TColgp_HArray1OfPnt(1, len(listOfVector))
for ix, v in enumerate(listOfVector):
pnts.SetValue(ix + 1, v.toPnt())
if smoothing:
spline_builder = GeomAPI_PointsToBSpline(
pnts, *smoothing, DegMax=maxDeg, Tol3D=tol
)
else:
spline_builder = GeomAPI_PointsToBSpline(
pnts, DegMin=minDeg, DegMax=maxDeg, Tol3D=tol
)
if not spline_builder.IsDone():
raise ValueError("B-spline approximation failed")
spline_geom = spline_builder.Curve()
return cls(BRepBuilderAPI_MakeEdge(spline_geom).Edge())
@classmethod
def makeThreePointArc(
cls: Type["Edge"], v1: Vector, v2: Vector, v3: Vector
@ -1602,6 +1672,19 @@ class Wire(Shape, Mixin1D):
return curve, BRepAdaptor_HCompCurve(curve)
def close(self) -> "Wire":
"""
Close a Wire
"""
if not self.IsClosed():
e = Edge.makeLine(self.endPoint(), self.startPoint())
rv = Wire.combine((self, e))[0]
else:
rv = self
return rv
@classmethod
def combine(
cls: Type["Wire"], listOfWires: Iterable[Union["Wire", Edge]], tol: float = 1e-9
@ -2022,6 +2105,47 @@ class Face(Shape):
return cls(face).fix()
@classmethod
def makeSplineApprox(
cls: Type["Face"],
points: List[List[Vector]],
tol: float = 1e-2,
smoothing: Optional[Tuple[float, float, float]] = None,
minDeg: int = 1,
maxDeg: int = 3,
) -> "Face":
"""
Approximate a spline surface through the provided points.
:param points: a 2D list of Vectors that represent the points
:param tol: tolerance of the algorithm (consult OCC documentation).
:param smoothing: optional tuple of 3 weights use for variational smoothing (default: None)
:param minDeg: minimum spline degree. Enforced only when smothing is None (default: 1)
:param maxDeg: maximum spline degree (default: 6)
:return: an Face
"""
points_ = TColgp_HArray2OfPnt(1, len(points), 1, len(points[0]))
for i, vi in enumerate(points):
for j, v in enumerate(vi):
points_.SetValue(i + 1, j + 1, v.toPnt())
if smoothing:
spline_builder = GeomAPI_PointsToBSplineSurface(
points_, *smoothing, DegMax=maxDeg, Tol3D=tol
)
else:
spline_builder = GeomAPI_PointsToBSplineSurface(
points_, DegMin=minDeg, DegMax=maxDeg, Tol3D=tol
)
if not spline_builder.IsDone():
raise ValueError("B-spline approximation failed")
spline_geom = spline_builder.Surface()
return cls(BRepBuilderAPI_MakeFace(spline_geom, Precision.Confusion_s()).Face())
def fillet2D(self, radius: float, vertices: Iterable[Vertex]) -> "Face":
"""
Apply 2D fillet to a face

1
doc/apireference.rst
View File

@ -52,6 +52,7 @@ All 2-d operations require a **Workplane** object to be created.
Workplane.move
Workplane.spline
Workplane.parametricCurve
Workplane.parametricSurface
Workplane.threePointArc
Workplane.sagittaArc
Workplane.radiusArc

83
tests/naca.py Executable file
View File

@ -0,0 +1,83 @@
naca5305 = [
[1.000074, 0.000520],
[0.998545, 0.000829],
[0.993968, 0.001750],
[0.986369, 0.003267],
[0.975792, 0.005351],
[0.962301, 0.007964],
[0.945974, 0.011059],
[0.926909, 0.014580],
[0.905218, 0.018465],
[0.881032, 0.022648],
[0.854494, 0.027058],
[0.825764, 0.031621],
[0.795017, 0.036263],
[0.762437, 0.040911],
[0.728224, 0.045493],
[0.692585, 0.049938],
[0.655739, 0.054180],
[0.617912, 0.058158],
[0.579336, 0.061813],
[0.540252, 0.065095],
[0.500900, 0.067958],
[0.461525, 0.070367],
[0.422374, 0.072291],
[0.383692, 0.073709],
[0.345722, 0.074611],
[0.308702, 0.074992],
[0.272257, 0.074521],
[0.237079, 0.072481],
[0.203612, 0.069020],
[0.172090, 0.064350],
[0.142727, 0.058704],
[0.115719, 0.052328],
[0.091240, 0.045475],
[0.069442, 0.038397],
[0.050454, 0.031335],
[0.034383, 0.024516],
[0.021313, 0.018141],
[0.011308, 0.012382],
[0.004410, 0.007379],
[0.000639, 0.003232],
[0.000000, 0.000000],
[0.002443, -0.002207],
[0.007902, -0.003318],
[0.016322, -0.003385],
[0.027630, -0.002492],
[0.041737, -0.000752],
[0.058539, 0.001696],
[0.077918, 0.004691],
[0.099743, 0.008055],
[0.123875, 0.011591],
[0.150167, 0.015098],
[0.178462, 0.018365],
[0.208603, 0.021185],
[0.240422, 0.023351],
[0.273752, 0.024670],
[0.308614, 0.024992],
[0.345261, 0.024967],
[0.382862, 0.024875],
[0.421191, 0.024682],
[0.460016, 0.024358],
[0.499100, 0.023878],
[0.538207, 0.023226],
[0.577098, 0.022390],
[0.615534, 0.021370],
[0.653278, 0.020171],
[0.690099, 0.018807],
[0.725767, 0.017298],
[0.760061, 0.015670],
[0.792768, 0.013955],
[0.823684, 0.012188],
[0.852613, 0.010407],
[0.879374, 0.008651],
[0.903799, 0.006957],
[0.925731, 0.005364],
[0.945032, 0.003906],
[0.961578, 0.002615],
[0.975264, 0.001519],
[0.986001, 0.000641],
[0.993720, 0.000001],
[0.998372, -0.000389],
[0.999926, -0.000520],
]

135
tests/test_cadquery.py
View File

@ -2189,13 +2189,15 @@ class TestCadQuery(BaseTest):
self.assertEqual(8, result.solids().item(0).faces().size())
def testSplitError(self):
"""
Test split produces the correct error when called with no solid to split.
"""
# Test split produces the correct error when called with no solid to split.
w = Workplane().hLine(1).vLine(1).close()
with raises(ValueError):
w.split(keepTop=True)
# Split should raise ValueError when called with no side kept
with raises(ValueError):
w.split(keepTop=False, keepBottom=False)
def testBoxDefaults(self):
"""
Tests creating a single box
@ -3291,6 +3293,11 @@ class TestCadQuery(BaseTest):
self.assertTrue(res_closed.solids().val().isValid())
self.assertEqual(len(res_closed.faces().vals()), 3)
res_edge = Workplane("XY").parametricCurve(func, makeWire=False)
self.assertEqual(len(res_edge.ctx.pendingEdges), 1)
self.assertEqual(len(res_edge.ctx.pendingWires), 0)
def testMakeShellSolid(self):
c0 = math.sqrt(2) / 4
@ -4387,3 +4394,125 @@ class TestCadQuery(BaseTest):
with raises(ValueError):
r.chamfer2D(0.25, [vs[0]])
def testSplineApprox(self):
from .naca import naca5305
from math import pi, cos
pts = [Vector(e[0], e[1], 0) for e in naca5305]
# spline
e1 = Edge.makeSplineApprox(pts, 1e-6, maxDeg=6, smoothing=(1, 1, 1))
e2 = Edge.makeSplineApprox(pts, 1e-6, minDeg=2, maxDeg=6)
self.assertTrue(e1.isValid())
self.assertTrue(e2.isValid())
self.assertTrue(e1.Length() > e2.Length())
with raises(ValueError):
e4 = Edge.makeSplineApprox(pts, 1e-6, maxDeg=3, smoothing=(1, 1, 1.0))
pts_closed = pts + [pts[0]]
e3 = Edge.makeSplineApprox(pts_closed)
w = Edge.makeSplineApprox(pts).close()
self.assertTrue(e3.IsClosed())
self.assertTrue(w.IsClosed())
# Workplane method
w1 = Workplane().splineApprox(pts)
w2 = Workplane().splineApprox(pts, forConstruction=True)
w3 = Workplane().splineApprox(pts, makeWire=True)
w4 = Workplane().splineApprox(pts, makeWire=True, forConstruction=True)
self.assertEqual(w1.edges().size(), 1)
self.assertEqual(len(w1.ctx.pendingEdges), 1)
self.assertEqual(w2.edges().size(), 1)
self.assertEqual(len(w2.ctx.pendingEdges), 0)
self.assertEqual(w3.wires().size(), 1)
self.assertEqual(len(w3.ctx.pendingWires), 1)
self.assertEqual(w4.wires().size(), 1)
self.assertEqual(len(w4.ctx.pendingWires), 0)
# spline surface
N = 40
T = 20
A = 5
pts = [
[
Vector(i, j, A * cos(2 * pi * i / T) * cos(2 * pi * j / T))
for i in range(N + 1)
]
for j in range(N + 1)
]
f1 = Face.makeSplineApprox(pts, smoothing=(1, 1, 1), maxDeg=6)
f2 = Face.makeSplineApprox(pts)
self.assertTrue(f1.isValid())
self.assertTrue(f2.isValid())
with raises(ValueError):
f3 = Face.makeSplineApprox(pts, smoothing=(1, 1, 1), maxDeg=3)
def testParametricSurface(self):
from math import pi, cos
r1 = Workplane().parametricSurface(
lambda u, v: (u, v, cos(pi * u) * cos(pi * v)), start=-1, stop=1
)
self.assertTrue(r1.faces().val().isValid())
r2 = Workplane().box(1, 1, 3).split(r1)
self.assertTrue(r2.solids().val().isValid())
self.assertEqual(r2.solids().size(), 2)
def testEdgeWireClose(self):
# test with edge
e0 = Edge.makeThreePointArc(Vector(0, 0, 0), Vector(1, 1, 0), Vector(0, 2, 0))
self.assertFalse(e0.IsClosed())
w0 = e0.close()
self.assertTrue(w0.IsClosed())
# test with already closed edge
e1 = Edge.makeCircle(1)
self.assertTrue(e1.IsClosed())
e2 = e1.close()
self.assertTrue(e2.IsClosed())
self.assertEqual(type(e1), type(e2))
# test with already closed WIRE
w1 = Wire.makeCircle(1, Vector(), Vector(0, 0, 1))
self.assertTrue(w1.IsClosed())
w2 = w1.close()
self.assertTrue(w1 is w2)
def testSplitShape(self):
"""
Testing the Shape.split method.
"""
# split an edge with a vertex
e0 = Edge.makeCircle(1, (0, 0, 0), (0, 0, 1))
v0 = Vertex.makeVertex(0, 1, 0)
list_of_edges = e0.split(v0).Edges()
self.assertEqual(len(list_of_edges), 2)
self.assertTrue(Vector(0, 1, 0) in [e.endPoint() for e in list_of_edges])
# split a circle with multiple vertices
angles = [2 * math.pi * idx / 10 for idx in range(10)]
vecs = [Vector(math.sin(a), math.cos(a), 0) for a in angles]
vertices = [Vertex.makeVertex(*v.toTuple()) for v in vecs]
edges = e0.split(*vertices).Edges()
self.assertEqual(len(edges), len(vertices) + 1)
endpoints = [e.endPoint() for e in edges]
self.assertTrue(all([v in endpoints for v in vecs]))