from __future__ import print_function
from __future__ import absolute_import
import sys
from six.moves import range
sys.path.append('.')
import numpy as nm
from sfepy.base.base import Struct
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
nm_f64_eps = nm.finfo(nm.float64).eps
[docs]
def to_ndarray(a):
if a is None:
return None
else:
a = nm.asarray(a)
if len(a.shape) == 0:
a = a.reshape(1)
return a
[docs]
class BSpline(Struct):
"""
B-spline curve representation
"""
def __init__(self, degree=3, is_cyclic=False, ncp=0):
"""
Initialize B-spline class.
Parameters
----------
degree : int
The degree of the spline function.
is_cyclic : bool
Cyclic spline?.
ncp : int
The number of control points.
"""
self.degree = degree
self.knot_type = None
self.is_cyclic = is_cyclic
self.ncp = ncp
self.knots = None
self.basis = None
self.curve_coors = None
self.cp_coors = None
self.approx_coors = None
self.t = None
[docs]
def set_control_points(self, coors, cyclic_form=False):
"""
Set the B-spline control points.
Parameters
----------
coors : array
The coordinates of unique control points.
cyclic_form : bool
Are the control points in the cyclic form?
"""
coors = to_ndarray(coors)
if self.is_cyclic and not cyclic_form:
coors = nm.vstack((coors, coors[0:self.degree,:]))
self.cp_coors = coors
self.ncp = coors.shape[0]
[docs]
def get_control_points(self):
"""
Get the B-spline control points.
Returns
-------
coors : array
The coordinates of control points.
"""
if self.is_cyclic:
return self.cp_coors[:-self.degree,:]
else:
return self.cp_coors
[docs]
def set_param(self, t):
"""
Set the B-spline parametric vector.
Parameters
----------
t : array
The parameter vector of the B-spline.
"""
self.t = to_ndarray(t)
if self.knots is not None:
endval = self.knots[-(self.degree + 1)]
idxs = nm.where(self.t == endval)[0]
self.t[idxs] -= nm_f64_eps
[docs]
def set_param_n(self, n=100, knot_range=(0.0, 1.0)):
"""
Generate the B-spline parametric vector using the number of steps.
Parameters
----------
n : array
The number of steps in the B-spline parametric vector.
"""
self.t = nm.linspace(knot_range[0], knot_range[1], n)
self.t[-1] -= nm_f64_eps
[docs]
@staticmethod
def basis_function_dg0(t, knots, n):
"""
Basis function: degree = 0
Parameters
----------
t : array
The parametric vector.
knots : array
The knot vector.
n : int
The number of intervals.
Returns
-------
bfun : array
The spline basis function evaluated for given values.
"""
nt = len(t)
bfun = nm.zeros((nt,n), dtype=nm.float64)
for ii in range(n):
idxs = nm.where(nm.logical_and(knots[ii] <= t,
t < knots[ii + 1]))[0]
bfun[idxs,ii] = 1.0
return bfun
[docs]
@staticmethod
def basis_function_dg(degree, t, knots, n):
"""
B-spline basis functions.
Parameters
----------
degree : int
The degree of the spline function.
t : array
The parametric vector.
knots : array
The knot vector.
n : int
The number of intervals.
Returns
-------
bfun : array
The spline basis function evaluated for given values.
"""
if degree >= 1:
bfun_dgm1 = BSpline.basis_function_dg(degree - 1, t,
knots, n + 1)
nt = len(t)
bfun = nm.zeros((nt,n), dtype=nm.float64)
for ii in range(n):
c1 = t - knots[ii]
c2 = knots[ii + degree] - knots[ii]
if nm.abs(c2) > nm_f64_eps:
bfun[:,ii] = c1 / c2 * bfun_dgm1[:,ii]
c1 = knots[ii + degree + 1] - t
c2 = knots[ii + degree + 1] - knots[ii + 1]
if nm.abs(c2) > nm_f64_eps:
bfun[:,ii] += c1 / c2 * bfun_dgm1[:,ii + 1]
else:
bfun = BSpline.basis_function_dg0(t, knots, n)
return bfun
[docs]
def make_knot_vector(self, knot_type='clamped', knot_data=None,
knot_range=(0.0, 1.0)):
"""
Create a knot vector of the requested type.
Parameters
----------
knot_type : str
The knot vector type: clamped/cyclic/userdef.
knot_data :
The extra knot data.
"""
if self.is_cyclic and 'cyclic' not in knot_type:
knot_type = 'cyclic'
ncp = self.ncp
dg = self.degree
n_knots = dg + ncp + 1
n_inter = n_knots - 2 * dg
aux = nm.linspace(knot_range[0], knot_range[1], n_inter)
if knot_type == '' or knot_type == 'cyclic':
dd = aux[1]
self.knots = nm.hstack((nm.arange(-dg, 0) * dd,
aux,
nm.arange(1, dg + 1) * dd + 1))
elif knot_type == 'clamped':
self.knots = nm.array([aux[0]]* dg + list(aux) + [aux[-1]]* dg,
dtype=nm.float64)
else:
raise NotImplementedError
self.knot_type = knot_type
[docs]
def set_knot_vector(self, knots):
"""
Set the knot vector.
Parameters
----------
knots : array
The knot vector.
"""
self.knot_type = 'userdef'
self.knots = to_ndarray(knots)
[docs]
def get_knot_vector(self):
"""
Return the knot vector.
Returns
-------
knots : array
The knot vector.
"""
return self.knots
[docs]
def insert_knot(self, new):
"""
Insert a new knot into the knot vector.
Parameters
----------
new : float
The new knot value.
"""
kn = self.knots
dg = self.degree
ncp = self.ncp
cp = self.cp_coors
idx = nm.where(nm.logical_and(kn[:-1] <= new, new < kn[1:]))[0]
if len(idx) > 0:
multi = len(nm.where(kn == new)[0])
if multi < dg:
# new knot
newkn = nm.zeros((len(kn) + 1,), dtype=nm.float64)
newkn[:(idx + 1)] = kn[:(idx + 1)]
newkn[idx + 1] = new
newkn[(idx + 2):] = kn[(idx + 1):]
u1 = idx - dg + 1
u2 = idx + 1
# new control points
newcp = nm.zeros((ncp + 1, cp.shape[1]), dtype=nm.float64)
newcp[:u1,:] = cp[:u1,:]
newcp[u2:,:] = cp[(u2 - 1):,:]
for ii in range(u1, u2):
kn1 = kn[ii]
kn2 = kn[ii + dg]
dd = kn2 - kn1
newcp[ii,:] = (kn2 - new) / dd * cp[ii - 1] + \
(new - kn1) / dd * cp[ii]
self.knots = newkn
self.cp_coors = newcp
self.ncp = newcp.shape[0]
# evaluate B-spline base functions for new configuration
self.eval_basis()
else:
print('knot insertion failed: multiplicity = spline degree!')
else:
print('knot insertion failed: out of bounds!')
[docs]
def eval_basis(self, t=None, return_val=False):
"""
Evaluate the basis of the bpsline.
Parameters
----------
t : array
The parameter vector of the B-spline.
"""
if t is not None:
self.set_param(t)
if self.knots is None:
self.make_knot_vector()
if self.t is None:
self.set_param_n()
self.basis = self.basis_function_dg(self.degree, self.t,
self.knots, self.ncp)
if return_val:
return self.basis
[docs]
def eval(self, t=None, cp_coors=None):
"""
Evaluate the coordinates of the bpsline curve.
Parameters
----------
t : array
The parameter vector of the B-spline.
cp_coors : array
The coordinates of the control points.
"""
if cp_coors is not None:
self.set_control_points(cp_coors)
self.eval_basis(t)
self.curve_coors = nm.dot(self.basis, self.cp_coors)
return self.curve_coors
[docs]
def draw(self, ret_ax=False, ax=None, color='r', cp_id=True):
"""
Draw B-spline curve.
Parameters
----------
ret_ax : bool
Return an axes object?
ax : axes object
The axes to which will be drawn.
color : str
Line color.
cp_id : bool
If True, label control points.
"""
if self.curve_coors is None:
self.eval()
cc = self.curve_coors
cp = self.cp_coors
ci = self.approx_coors
if ci is not None and self.is_cyclic:
ci = nm.vstack((ci, ci[0,:]))
if cc.shape[1] == 3:
if ax is None:
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(cc[:,0], cc[:,1], cc[:,2], color + '-')
if cp_id:
ax.plot(cp[:,0], cp[:,1], cp[:,2], 'ko:', alpha=0.6)
if ci is not None:
ax.plot(ci[:,0], ci[:,1], ci[:,2], 'b--', alpha=0.6)
else:
if ax is None:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(cc[:,0], cc[:,1], color + '-')
if cp_id:
ax.plot(cp[:,0], cp[:,1], 'ko:', alpha=0.6)
for jj, icp in enumerate(self.cp_coors):
ax.text(icp[0], icp[1], 'N%d' % (jj + 1), fontsize=10)
if ci is not None:
ax.plot(ci[:,0], ci[:,1], 'b--', alpha=0.6)
ax.set_aspect('equal')
if ret_ax:
return ax
else:
plt.show()
[docs]
def draw_basis(self):
"""
Draw B-spline curve.
"""
plt.figure()
plt.plot(self.t,self.basis)
plt.legend(['b%d' % (ii + 1) for ii in range(self.basis.shape[1])])
plt.show()
[docs]
def approximate(self, coors, ncp=None, knot_type='clamped',
knots=None, alpha=0.5,
do_eval=False, do_param_correction=False):
"""
Approximate set of points by the B-spline curve.
Parameters
----------
coors : array
The coordinates of the approximated points.
ncp : int
The number of control points.
knot_type : str
The knot vector type.
knots : array
The knot vector.
alpha : float
The parameter vector distribution:
1.0 = chordal
0.5 = centripetal
do_eval : bool
Evaluate the curve coordinates?
do_param_correction : bool
Perform parametric corrections to improve the approximation?
"""
coors = to_ndarray(coors)
dg = self.degree
if ncp is not None:
if self.is_cyclic:
ncp += dg
self.ncp = ncp
self.make_knot_vector(knot_type)
if knots is not None:
self.knots = knots
self.knot_type = 'userdef'
ncp = len(knots) - dg - 1
self.ncp = ncp
# param vector
dist = nm.linalg.norm(coors[:-1,:] - coors[1:,:], axis=1)
dista = nm.power(dist, alpha)
self.t = nm.zeros((coors.shape[0],), dtype=nm.float64)
self.t[1:] += dista.cumsum() / dista.sum()
self.t[-1] -= nm_f64_eps
while True:
self.basis = self.basis_function_dg(dg, self.t,
self.knots, ncp)
A = nm.dot(self.basis.T, self.basis)
b = nm.dot(self.basis.T, coors)
# cyclic spline
if self.is_cyclic:
nred = ncp - dg
R = nm.zeros((ncp, nred), dtype=nm.float64)
for ii in range(nred):
R[ii,ii] = 1.0
for ii in range(self.degree):
R[nred + ii,ii] = 1.0
A = nm.dot(R.T, nm.dot(A, R))
b = nm.dot(R.T, b)
self.cp_coors = nm.dot(R, nm.dot(nm.linalg.inv(A), b))
else:
self.cp_coors = nm.dot(nm.linalg.inv(A), b)
self.approx_coors = coors
if not do_param_correction:
break
if do_eval:
self.curve_coors = nm.dot(self.basis, self.cp_coors)
[docs]
def set_approx_points(self, coors):
"""
Set the coordinates of approximated points.
Parameters
----------
coors : array
The coordinates of approximated points.
"""
self.approx_coors = to_ndarray(coors)
[docs]
class BSplineSurf(Struct):
"""
B-spline surface representation
"""
def __init__(self, degree=(3,3), is_cyclic=(False, False)):
"""
Initialize B-spline class.
Parameters
----------
degree : tuple of int
The degree of the spline functions.
is_cyclic : tuple of bool
Cyclic splines?.
"""
self.splines = [None, None]
for ii in range(2):
self.splines[ii] = BSpline(degree[ii], is_cyclic=is_cyclic[ii])
self.surf_coors = None
self.cp_coors = None
self.approx_coors = None
[docs]
def set_control_points(self, coors, cyclic_form=False):
"""
Set the B-spline control points.
Parameters
----------
coors : array
The coordinates of unique control points.
cyclic_form : bool
Are the control points in the cyclic form?
"""
coors = to_ndarray(coors)
if self.splines[0].is_cyclic and not cyclic_form:
coors = nm.vstack((coors, coors[0:self.splines[0].degree,:,:]))
if self.splines[1].is_cyclic and not cyclic_form:
coors = nm.hstack((coors, coors[:,0:self.splines[1].degree,:]))
self.cp_coors = coors
for ii in range(2):
self.splines[ii].ncp = coors.shape[ii]
[docs]
def get_control_points(self):
"""
Get the B-spline surface control points.
Returns
-------
coors : array
The coordinates of control points.
"""
aux = self.cp_coors
if self.splines[0].is_cyclic:
aux = aux[:-self.splines[0].degree,:,:]
if self.splines[1].is_cyclic:
aux = aux[:,:-self.splines[1].degree,:]
return aux
[docs]
def make_knot_vector(self, knot_type=('clamped', 'clamped'),
knot_data=(None, None)):
"""
Create a knot vector of the requested type.
Parameters
----------
knot_type : tuple of str
The knot vector types.
knot_data : tuple of ANY
The extra knot data.
"""
for ii in range(2):
self.splines[ii].make_knot_vector(knot_type[ii], knot_data[ii])
[docs]
def set_param_n(self, n=(100, 100)):
"""
Generate the B-spline parametric vector using the number of steps.
Parameters
----------
n : tuple of array
The number of steps in the B-spline parametric vectors.
"""
for ii in range(2):
self.splines[ii].set_param_n(n[ii])
[docs]
def set_approx_points(self, coors):
"""
Set the coordinates of approximated points.
Parameters
----------
coors : array
The coordinates of approximated points.
"""
self.approx_coors = to_ndarray(coors)
[docs]
def eval(self, t=(None, None), cp_coors=None):
"""
Evaluate the coordinates of the bpsline curve.
Parameters
----------
t : tuple of array
The parametric vector of the B-splines.
cp_coors : array
The coordinates of the control points.
"""
if cp_coors is not None:
self.set_control_points(cp_coors)
for ii in range(2):
self.splines[ii].eval_basis(t[ii])
nt = (len(self.splines[0].t), len(self.splines[1].t))
ncp = (self.splines[0].ncp, self.splines[1].ncp)
aux = nm.zeros((nt[0], ncp[1], 3), dtype=nm.float64)
for ii in range(ncp[1]):
aux[:,ii,:] = nm.dot(self.splines[0].basis, self.cp_coors[:,ii,:])
self.surf_coors = nm.zeros(nt + (3,), dtype=nm.float64)
for ii in range(nt[0]):
self.surf_coors[ii,:,:] = nm.dot(self.splines[1].basis, aux[ii,:,:])
return self.surf_coors
[docs]
def draw(self, ret_ax=False, ax=None):
"""
Draw B-spline surface.
Parameters
----------
ret_ax : bool
Return an axes object?
ax : axes object
The axes to which will be drawn.
"""
if self.surf_coors is None:
self.eval()
fig = plt.figure()
ax = Axes3D(fig)
coors = self.surf_coors
cs = coors.shape
for ii in range(cs[0] - 1):
for jj in range(cs[1] - 1):
verts = nm.array([coors[ii,jj,:],
coors[ii,jj + 1,:],
coors[ii + 1,jj + 1,:],
coors[ii + 1,jj,:]])
quad = Poly3DCollection([verts],
facecolors='gray', edgecolor='k',
linewidths=0.2, alpha=0.5)
ax.add_collection3d(quad)
cp = self.cp_coors
for ii in range(cp.shape[1]):
ax.plot(cp[:,ii,0], cp[:,ii,1], cp[:,ii,2], 'ro--', linewidth=2.0)
for ii in range(cp.shape[0]):
ax.plot(cp[ii,:,0], cp[ii,:,1], cp[ii,:,2], 'ro--', linewidth=2.0)
ax.set_aspect('equal')
plt.show()
[docs]
def approximate(self, coors, ncp, do_eval=False):
"""
Approximate set of points by the B-spline surface.
Parameters
----------
coors : array
The coordinates of the approximated points.
ncp : tuple of int
The number of control points.
"""
coors = to_ndarray(coors)
nsc = coors.shape[0:2]
aux = nm.zeros((nsc[0], ncp[1], 3), dtype=nm.float64)
spl1 = self.splines[1]
for ii in range(nsc[0]):
spl1.approximate(coors[ii,...], ncp[1])
aux[ii,...] = spl1.get_control_points()
self.cp_coors = nm.zeros((ncp[0], ncp[1], 3), dtype=nm.float64)
spl2 = self.splines[0]
for ii in range(ncp[1]):
spl2.approximate(aux[:,ii,:], ncp[0])
self.cp_coors[:,ii,:] = spl2.get_control_points()
self.approx_coors = coors
[docs]
def write_surface_vtk(self, filename, float_format='%.6f'):
"""
Write the spline surface to VTK file.
Parameters
----------
filename: str
Name of the VTK file.
float_format: str
Float formating.
"""
coors = self.surf_coors
cs0, cs1 = coors.shape[0:2]
nquads = (cs0 - 1) * (cs1 - 1)
quads = nm.zeros((nquads, 4), dtype=nm.int64)
kk = 0
for ii in range(cs0 - 1):
offs = ii * cs1
for jj in range(cs1 - 1):
quads[kk] = nm.array([jj + offs,
jj + offs + cs1,
jj + 1 + offs + cs1,
jj + 1 + offs])
kk += 1
f = open(filename, 'w')
f.write('# vtk DataFile Version 2.0\n')
f.write('spline surface\nASCII\nDATASET POLYDATA\n')
ff3 = ' '.join([float_format] * 3) + '\n'
f.write('POINTS %d float\n' % (cs0 * cs1))
for ii in range(cs0):
offs = ii * cs1
for jj in range(cs1):
f.write(ff3 % tuple(coors[ii,jj,:]))
f.write('POLYGONS %d %d\n' % (nquads, nquads * 5))
for ii in quads:
f.write('4 %s\n' % (' '.join([str(jj) for jj in ii])))
f.close()
[docs]
def write_control_polygon_vtk(self, filename, float_format='%.6f'):
"""
Write the control polygon to VTK file.
Parameters
----------
filename: str
Name of the VTK file.
float_format: str
Float formating.
"""
coors = self.cp_coors
cs0, cs1 = coors.shape[0:2]
lines = []
nlines = cs0 + cs1
nlpts = 0
for ii in range(cs0):
lines.append(nm.arange(cs1) + ii * cs1)
nlpts += cs1
for ii in range(cs1):
lines.append(nm.arange(cs0) * cs1 + ii)
nlpts += cs0
f = open(filename, 'w')
f.write('# vtk DataFile Version 2.0\n')
f.write('spline control polygon\nASCII\nDATASET POLYDATA\n')
ff3 = ' '.join([float_format] * 3) + '\n'
f.write('POINTS %d float\n' % (cs0 * cs1))
for ii in range(cs0):
for jj in range(cs1):
f.write(ff3 % tuple(coors[ii,jj,:]))
f.write('LINES %d %d\n' % (nlines, nlines + nlpts))
for ii in lines:
f.write('%d %s\n' % (len(ii), ' '.join([str(jj) for jj in ii])))
f.close()
[docs]
def get_2d_points(is3d=False):
"""
Returns the set of points.
Parameters
----------
is3d : bool
3D coordinates?
"""
out = nm.array([(-0.87, 0.15, 0),
(-0.70, 0.54, 0),
(-0.32, 0.80, 0),
(0.15, 0.70, 0),
(0.37, 0.26, 0),
(0.70, -0.07, 0),
(0.67, -0.49, 0),
(0.07, -0.81, 0),
(-0.44, -0.72, 0),
(-0.80, -0.34, 0)])
if is3d:
return out
else:
return out[:,:2]
[docs]
def approximation_example():
"""
The example of using BSplineSurf for approximation
of the surface given by the set of points.
"""
# define sample points in 2D
spts0 = get_2d_points(is3d=True)
# define sample points in 3D
spts = nm.array([spts0,
spts0 * 0.7 + nm.array([0.1,0,0.5]),
spts0 * 0.8 + nm.array([0.2,0,1.5]),
spts0 * 1.2 + nm.array([0.4,0,2.5])])
cyclic=(False, True)
spl1 = BSplineSurf((3, 3), is_cyclic=cyclic)
spl1.approximate(spts, (4,8))
cp = spl1.get_control_points()
spls2 = BSplineSurf((3, 3), is_cyclic=cyclic)
spls2.set_control_points(cp)
spls2.make_knot_vector()
spls2.set_param_n((12, 24))
spls2.eval()
spls2.draw()
[docs]
def simple_example():
"""
The example of using B-spline class.
"""
# define control points in 2D
cp = get_2d_points()
spl = BSpline(3, is_cyclic=True)
spl.set_control_points(cp)
spl.make_knot_vector()
spl.set_param_n(150)
spl.insert_knot(0.7)
spl.insert_knot(0.7)
spl.insert_knot(0.7)
spl.eval()
spl.draw()
[docs]
def main(argv):
# simple_example()
approximation_example()
if __name__ == '__main__':
main(sys.argv)