From 8d73ee1225880357ea251f8026b659ae58447be1 Mon Sep 17 00:00:00 2001 From: Zetison Date: Thu, 6 Feb 2020 14:40:56 +0100 Subject: [PATCH 1/3] Implemented explicit version of raise_order for open knot vectors related to issue #14 --- splipy/SplineObject.py | 107 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 106 insertions(+), 1 deletion(-) diff --git a/splipy/SplineObject.py b/splipy/SplineObject.py index 63e437e2..f46049ca 100644 --- a/splipy/SplineObject.py +++ b/splipy/SplineObject.py @@ -1,5 +1,6 @@ # -*- coding: utf-8 -*- +import splipy.state as state import numpy as np import copy from operator import attrgetter, methodcaller @@ -426,7 +427,7 @@ def set_order(self, *order): def raise_order(self, *raises): """ Raise the polynomial order of the object. If only one argument is - given, the order is raised equally over all directions. + given, the order is raised equally over all directions. The explicit version is only implemented on open knot vectors. The function raise_order_implicit is used otherwise :param int u,v,...: Number of times to raise the order in a given direction. @@ -438,7 +439,111 @@ def raise_order(self, *raises): raise ValueError("Cannot lower order using raise_order") if all(r == 0 for r in raises): return + + if any(abs(b.knots[0] - b.knots[b.order-1]) > state.knot_tolerance or b.periodic > -1 for b in self.bases): + self.raise_order_implicit(*raises) + return + + new_bases = [b.raise_order(r) for b, r in zip(self.bases, raises)] + + def raise_order_1D(n, k, T, P, m, periodic): + """ Implementation of method in "Efficient Degree Elevation and Knot Insertion + for B-spline Curves using Derivatives" by Qi-Xing Huang a Shi-Min Hu, Ralph R Martin. Only the case of open knot vector is fully implemented + :param int n: (n+1) is the number of initial basis functions + :param int k: spline order + :param T: knot vector + :param P: weighted NURBS coefficients + :param int m: number of degree elevations + :param int periodic: Number of continuous derivatives at start and end. -1 is not periodic, 0 is continuous, etc. + :return Q: new control points + """ + + u = np.unique(T[k-1:-k+1]) + S = u.size - 1 + d = P.shape[0] # dimension of spline + + # Find multiplicity of the knot vector T + z = [1] * (S+1) + l = 1 + for i in range(k,T.size-k-1): + if abs(T[i] - T[i+1]) < state.knot_tolerance: + z[l] += 1 + else: + l += 1 + + # Step 1: Find Pt_i^j + Pt = np.zeros((d,n+1,k)) + Pt[:,:,0] = P + Pt = np.concatenate((Pt,Pt[:,0:periodic+1,:]),axis=1) + n += periodic+1 + + beta = np.cumsum(z[1:-1],dtype=int) + beta = np.insert(beta,0,0) # include the empty sum (=0) + for l in range(1,k): + for i in range(0,n+1-l): + if T[i+l] < T[i+k]: + Pt[:,i,l] = (Pt[:,i+1,l-1] - Pt[:,i,l-1])/(T[i+k]-T[i+l]) + + # Step 2: Create new knot vector Tb + nb = n + S*m + Tb = np.zeros(nb+m+k+1) + Tb[:k-1] = T[:k-1] + Tb[-k+1:] = T[-k+1:] + j = k-1 + for i in range(0,len(z)): + Tb[j:j+z[i]+m] = u[i] + j = j+z[i]+m + + # Step 3: Find boundary values of Qt_i^j + Qt = np.zeros((d,nb+1,k)) + l_arr = np.array(range(1,k)) + alpha = np.cumprod((k-l_arr)/(k+m-l_arr)) + alpha = np.insert(alpha,0,1) # include the empty product (=1) + indices = range(0,k) + Qt[:,0,indices] = np.multiply(Pt[:,0,indices],np.reshape(alpha[indices],(1,1,k))) # (21) + for p in range(0,S): + indices = range(k-z[p],k) + Qt[:,beta[p]+p*m,indices] = np.multiply(Pt[:,beta[p],indices],np.reshape(alpha[indices],(1,1,z[p]))) # (22) + + for p in range(0,S): + idx = beta[p]+p*m + Qt[:,idx+1:m+idx+1,k-1] = np.repeat(Qt[:,idx:idx+1,k-1],m,axis=1) # (23) + + # Step 4: Find remaining values of Qt_i^j + for j in range(k-1,0,-1): + for i in range(0,nb): + if Tb[i+k+m] > Tb[i+j]: + Qt[:,i+1,j-1] = Qt[:,i,j-1] + (Tb[i+k+m] - Tb[i+j])*Qt[:,i,j] #(20) with Qt replacing Pt + + return Qt[:,:,0] + + d_p = self.pardim + + controlpoints = self.controlpoints + for i in range(0,d_p): + dimensions = np.array(controlpoints.shape) + indices = np.array(range(0,d_p+1)) + indices[i],indices[d_p] = d_p,i + controlpoints = np.transpose(controlpoints,indices) + controlpoints = np.reshape(controlpoints,(np.prod(dimensions[indices[:-1]]),dimensions[i])) + controlpoints = raise_order_1D(controlpoints.shape[1]-1,self.order(i), + self.bases[i].knots,controlpoints,raises[i],self.bases[i].periodic) + controlpoints = np.reshape(controlpoints,np.append(dimensions[indices[:-1]],controlpoints.shape[1])) + controlpoints = np.transpose(controlpoints,indices) + + self.controlpoints = controlpoints + self.bases = new_bases + + return self + def raise_order_implicit(self, *raises): + """ Raise the polynomial order of the object. If only one argument is + given, the order is raised equally over all directions. + :param int u,v,...: Number of times to raise the order in a given + direction. + :return: self + """ + new_bases = [b.raise_order(r) for b, r in zip(self.bases, raises)] # Set up an interpolation problem From 4b98c43aa664faebc00c233c19ef6ad4fc30df4e Mon Sep 17 00:00:00 2001 From: Zetison Date: Fri, 7 Feb 2020 14:17:32 +0100 Subject: [PATCH 2/3] Added suggestions --- splipy/SplineObject.py | 79 ++-------------------------------------- splipy/utils/__init__.py | 68 +++++++++++++++++++++++++++++++++- 2 files changed, 71 insertions(+), 76 deletions(-) diff --git a/splipy/SplineObject.py b/splipy/SplineObject.py index f46049ca..720fe5f9 100644 --- a/splipy/SplineObject.py +++ b/splipy/SplineObject.py @@ -1,6 +1,5 @@ # -*- coding: utf-8 -*- -import splipy.state as state import numpy as np import copy from operator import attrgetter, methodcaller @@ -8,7 +7,8 @@ from bisect import bisect_left from splipy import BSplineBasis from splipy.utils import reshape, rotation_matrix, is_singleton, ensure_listlike,\ - check_direction, ensure_flatlist, check_section, sections + check_direction, ensure_flatlist, check_section, sections,\ + raise_order_1D __all__ = ['SplineObject'] @@ -440,83 +440,12 @@ def raise_order(self, *raises): if all(r == 0 for r in raises): return - if any(abs(b.knots[0] - b.knots[b.order-1]) > state.knot_tolerance or b.periodic > -1 for b in self.bases): + if any(b.continuity(b.knots[0]) < b.order or b.periodic > -1 for b in self.bases): self.raise_order_implicit(*raises) return new_bases = [b.raise_order(r) for b, r in zip(self.bases, raises)] - - def raise_order_1D(n, k, T, P, m, periodic): - """ Implementation of method in "Efficient Degree Elevation and Knot Insertion - for B-spline Curves using Derivatives" by Qi-Xing Huang a Shi-Min Hu, Ralph R Martin. Only the case of open knot vector is fully implemented - :param int n: (n+1) is the number of initial basis functions - :param int k: spline order - :param T: knot vector - :param P: weighted NURBS coefficients - :param int m: number of degree elevations - :param int periodic: Number of continuous derivatives at start and end. -1 is not periodic, 0 is continuous, etc. - :return Q: new control points - """ - - u = np.unique(T[k-1:-k+1]) - S = u.size - 1 - d = P.shape[0] # dimension of spline - - # Find multiplicity of the knot vector T - z = [1] * (S+1) - l = 1 - for i in range(k,T.size-k-1): - if abs(T[i] - T[i+1]) < state.knot_tolerance: - z[l] += 1 - else: - l += 1 - - # Step 1: Find Pt_i^j - Pt = np.zeros((d,n+1,k)) - Pt[:,:,0] = P - Pt = np.concatenate((Pt,Pt[:,0:periodic+1,:]),axis=1) - n += periodic+1 - - beta = np.cumsum(z[1:-1],dtype=int) - beta = np.insert(beta,0,0) # include the empty sum (=0) - for l in range(1,k): - for i in range(0,n+1-l): - if T[i+l] < T[i+k]: - Pt[:,i,l] = (Pt[:,i+1,l-1] - Pt[:,i,l-1])/(T[i+k]-T[i+l]) - - # Step 2: Create new knot vector Tb - nb = n + S*m - Tb = np.zeros(nb+m+k+1) - Tb[:k-1] = T[:k-1] - Tb[-k+1:] = T[-k+1:] - j = k-1 - for i in range(0,len(z)): - Tb[j:j+z[i]+m] = u[i] - j = j+z[i]+m - - # Step 3: Find boundary values of Qt_i^j - Qt = np.zeros((d,nb+1,k)) - l_arr = np.array(range(1,k)) - alpha = np.cumprod((k-l_arr)/(k+m-l_arr)) - alpha = np.insert(alpha,0,1) # include the empty product (=1) - indices = range(0,k) - Qt[:,0,indices] = np.multiply(Pt[:,0,indices],np.reshape(alpha[indices],(1,1,k))) # (21) - for p in range(0,S): - indices = range(k-z[p],k) - Qt[:,beta[p]+p*m,indices] = np.multiply(Pt[:,beta[p],indices],np.reshape(alpha[indices],(1,1,z[p]))) # (22) - - for p in range(0,S): - idx = beta[p]+p*m - Qt[:,idx+1:m+idx+1,k-1] = np.repeat(Qt[:,idx:idx+1,k-1],m,axis=1) # (23) - - # Step 4: Find remaining values of Qt_i^j - for j in range(k-1,0,-1): - for i in range(0,nb): - if Tb[i+k+m] > Tb[i+j]: - Qt[:,i+1,j-1] = Qt[:,i,j-1] + (Tb[i+k+m] - Tb[i+j])*Qt[:,i,j] #(20) with Qt replacing Pt - - return Qt[:,:,0] - + d_p = self.pardim controlpoints = self.controlpoints diff --git a/splipy/utils/__init__.py b/splipy/utils/__init__.py index 41d7d599..c4ab142b 100644 --- a/splipy/utils/__init__.py +++ b/splipy/utils/__init__.py @@ -166,10 +166,76 @@ def uniquify(iterator): seen.add(i) yield i +def raise_order_1D(n, k, T, P, m, periodic): + """ Implementation of method in "Efficient Degree Elevation and Knot Insertion + for B-spline Curves using Derivatives" by Qi-Xing Huang a Shi-Min Hu, Ralph R Martin. Only the case of open knot vector is fully implemented + :param int n: (n+1) is the number of initial basis functions + :param int k: spline order + :param T: knot vector + :param P: weighted NURBS coefficients + :param int m: number of degree elevations + :param int periodic: Number of continuous derivatives at start and end. -1 is not periodic, 0 is continuous, etc. + :return Q: new control points + """ + + u = np.unique(T[k-1:-k+1]) + S = u.size - 1 + d = P.shape[0] # dimension of spline + + # Find multiplicity of the knot vector T + b = BSplineBasis(k, T) + z = [k-1-b.continuity(t0) for t0 in b.knot_spans()] + + # Step 1: Find Pt_i^j + Pt = np.zeros((d,n+1,k)) + Pt[:,:,0] = P + Pt = np.concatenate((Pt,Pt[:,0:periodic+1,:]),axis=1) + n += periodic+1 + + beta = np.cumsum(z[1:-1],dtype=int) + beta = np.insert(beta,0,0) # include the empty sum (=0) + for l in range(1,k): + for i in range(0,n+1-l): + if T[i+l] < T[i+k]: + Pt[:,i,l] = (Pt[:,i+1,l-1] - Pt[:,i,l-1])/(T[i+k]-T[i+l]) + + # Step 2: Create new knot vector Tb + nb = n + S*m + Tb = np.zeros(nb+m+k+1) + Tb[:k-1] = T[:k-1] + Tb[-k+1:] = T[-k+1:] + j = k-1 + for i in range(0,len(z)): + Tb[j:j+z[i]+m] = u[i] + j = j+z[i]+m + + # Step 3: Find boundary values of Qt_i^j + Qt = np.zeros((d,nb+1,k)) + l_arr = np.array(range(1,k)) + alpha = np.cumprod((k-l_arr)/(k+m-l_arr)) + alpha = np.insert(alpha,0,1) # include the empty product (=1) + indices = range(0,k) + Qt[:,0,indices] = np.multiply(Pt[:,0,indices],np.reshape(alpha[indices],(1,1,k))) # (21) + for p in range(0,S): + indices = range(k-z[p],k) + Qt[:,beta[p]+p*m,indices] = np.multiply(Pt[:,beta[p],indices],np.reshape(alpha[indices],(1,1,z[p]))) # (22) + + for p in range(0,S): + idx = beta[p]+p*m + Qt[:,idx+1:m+idx+1,k-1] = np.repeat(Qt[:,idx:idx+1,k-1],m,axis=1) # (23) + + # Step 4: Find remaining values of Qt_i^j + for j in range(k-1,0,-1): + for i in range(0,nb): + if Tb[i+k+m] > Tb[i+j]: + Qt[:,i+1,j-1] = Qt[:,i,j-1] + (Tb[i+k+m] - Tb[i+j])*Qt[:,i,j] #(20) with Qt replacing Pt + + return Qt[:,:,0] + __all__ = [ 'nutils', 'refinement', 'image', 'NACA', 'curve', 'smooth', 'rotation_matrix', 'sections', 'section_from_index', 'section_to_index', 'check_section', 'check_direction', 'ensure_flatlist', 'is_singleton', 'ensure_listlike', 'rotate_local_x_axis', 'flip_and_move_plane_geometry', - 'reshape', + 'reshape','raise_order_1D' ] From f2003ad606bbac6f427c273a2bbe7f24f7f6e49e Mon Sep 17 00:00:00 2001 From: Zetison Date: Wed, 19 Feb 2020 11:45:22 +0100 Subject: [PATCH 3/3] Moved raise_order_1D to the end of SplineObject.py --- splipy/Curve.py | 2 +- splipy/SplineObject.py | 77 +++++++++++++++++++++++++++++++++++++--- splipy/utils/__init__.py | 66 ---------------------------------- 3 files changed, 74 insertions(+), 71 deletions(-) diff --git a/splipy/Curve.py b/splipy/Curve.py index e3253c0f..123e15ef 100644 --- a/splipy/Curve.py +++ b/splipy/Curve.py @@ -263,7 +263,7 @@ def torsion(self, t, above=True): return nominator / np.power(magnitude, 2) - def raise_order(self, amount): + def raise_order_implicit(self, amount): """ Raise the polynomial order of the curve. :param int amount: Number of times to raise the order diff --git a/splipy/SplineObject.py b/splipy/SplineObject.py index 720fe5f9..4c36caf7 100644 --- a/splipy/SplineObject.py +++ b/splipy/SplineObject.py @@ -7,8 +7,8 @@ from bisect import bisect_left from splipy import BSplineBasis from splipy.utils import reshape, rotation_matrix, is_singleton, ensure_listlike,\ - check_direction, ensure_flatlist, check_section, sections,\ - raise_order_1D + check_direction, ensure_flatlist, check_section, sections +import splipy.state as state __all__ = ['SplineObject'] @@ -440,12 +440,12 @@ def raise_order(self, *raises): if all(r == 0 for r in raises): return - if any(b.continuity(b.knots[0]) < b.order or b.periodic > -1 for b in self.bases): + if any(not(b.continuity(b.knots[0]) == -1) or b.periodic > -1 for b in self.bases): self.raise_order_implicit(*raises) return new_bases = [b.raise_order(r) for b, r in zip(self.bases, raises)] - + d_p = self.pardim controlpoints = self.controlpoints @@ -1428,3 +1428,72 @@ def make_splines_identical(cls, spline1, spline2): m = min(c1-c2, p[i]-1-c2) # c1=np.inf if knot does not exist spline1.insert_knot([k]*m, direction=i) +def raise_order_1D(n, k, T, P, m, periodic): + """ Implementation of method in "Efficient Degree Elevation and Knot Insertion + for B-spline Curves using Derivatives" by Qi-Xing Huang a Shi-Min Hu, Ralph R Martin. Only the case of open knot vector is fully implemented + :param int n: (n+1) is the number of initial basis functions + :param int k: spline order + :param T: knot vector + :param P: weighted NURBS coefficients + :param int m: number of degree elevations + :param int periodic: Number of continuous derivatives at start and end. -1 is not periodic, 0 is continuous, etc. + :return Q: new control points + """ + + d = P.shape[0] # dimension of spline + + # Find multiplicity of the knot vector T + b = BSplineBasis(k, T, periodic) + z = [k-1-b.continuity(t0) for t0 in b.knot_spans()] + z[0] = 1 + z[-1] = 1 + + u = b.knot_spans(include_ghost_knots=False) + S = len(u) - 1 + + # Step 1: Find Pt_i^j + Pt = np.zeros((d,n+1,k)) + Pt[:,:,0] = P + Pt = np.concatenate((Pt,Pt[:,0:periodic+1,:]),axis=1) + n += periodic+1 + + beta = np.cumsum(z[1:-1],dtype=int) + beta = np.insert(beta,0,0) # include the empty sum (=0) + for l in range(1,k): + for i in range(0,n+1-l): + if T[i+l] < T[i+k]: + Pt[:,i,l] = (Pt[:,i+1,l-1] - Pt[:,i,l-1])/(T[i+k]-T[i+l]) + + # Step 2: Create new knot vector Tb + nb = n + S*m + Tb = np.zeros(nb+m+k+1) + Tb[:k-1] = T[:k-1] + Tb[-k+1:] = T[-k+1:] + j = k-1 + for i in range(0,len(z)): + Tb[j:j+z[i]+m] = u[i] + j = j+z[i]+m + + # Step 3: Find boundary values of Qt_i^j + Qt = np.zeros((d,nb+1,k)) + l_arr = np.array(range(1,k)) + alpha = np.cumprod((k-l_arr)/(k+m-l_arr)) + alpha = np.insert(alpha,0,1) # include the empty product (=1) + indices = range(0,k) + Qt[:,0,indices] = np.multiply(Pt[:,0,indices],np.reshape(alpha[indices],(1,1,k))) # (21) + for p in range(0,S): + indices = range(k-z[p],k) + Qt[:,beta[p]+p*m,indices] = np.multiply(Pt[:,beta[p],indices],np.reshape(alpha[indices],(1,1,z[p]))) # (22) + + for p in range(0,S): + idx = beta[p]+p*m + Qt[:,idx+1:m+idx+1,k-1] = np.repeat(Qt[:,idx:idx+1,k-1],m,axis=1) # (23) + + # Step 4: Find remaining values of Qt_i^j + for j in range(k-1,0,-1): + for i in range(0,nb): + if Tb[i+k+m] > Tb[i+j]: + Qt[:,i+1,j-1] = Qt[:,i,j-1] + (Tb[i+k+m] - Tb[i+j])*Qt[:,i,j] #(20) with Qt replacing Pt + + return Qt[:,:,0] + diff --git a/splipy/utils/__init__.py b/splipy/utils/__init__.py index c4ab142b..5e845321 100644 --- a/splipy/utils/__init__.py +++ b/splipy/utils/__init__.py @@ -166,72 +166,6 @@ def uniquify(iterator): seen.add(i) yield i -def raise_order_1D(n, k, T, P, m, periodic): - """ Implementation of method in "Efficient Degree Elevation and Knot Insertion - for B-spline Curves using Derivatives" by Qi-Xing Huang a Shi-Min Hu, Ralph R Martin. Only the case of open knot vector is fully implemented - :param int n: (n+1) is the number of initial basis functions - :param int k: spline order - :param T: knot vector - :param P: weighted NURBS coefficients - :param int m: number of degree elevations - :param int periodic: Number of continuous derivatives at start and end. -1 is not periodic, 0 is continuous, etc. - :return Q: new control points - """ - - u = np.unique(T[k-1:-k+1]) - S = u.size - 1 - d = P.shape[0] # dimension of spline - - # Find multiplicity of the knot vector T - b = BSplineBasis(k, T) - z = [k-1-b.continuity(t0) for t0 in b.knot_spans()] - - # Step 1: Find Pt_i^j - Pt = np.zeros((d,n+1,k)) - Pt[:,:,0] = P - Pt = np.concatenate((Pt,Pt[:,0:periodic+1,:]),axis=1) - n += periodic+1 - - beta = np.cumsum(z[1:-1],dtype=int) - beta = np.insert(beta,0,0) # include the empty sum (=0) - for l in range(1,k): - for i in range(0,n+1-l): - if T[i+l] < T[i+k]: - Pt[:,i,l] = (Pt[:,i+1,l-1] - Pt[:,i,l-1])/(T[i+k]-T[i+l]) - - # Step 2: Create new knot vector Tb - nb = n + S*m - Tb = np.zeros(nb+m+k+1) - Tb[:k-1] = T[:k-1] - Tb[-k+1:] = T[-k+1:] - j = k-1 - for i in range(0,len(z)): - Tb[j:j+z[i]+m] = u[i] - j = j+z[i]+m - - # Step 3: Find boundary values of Qt_i^j - Qt = np.zeros((d,nb+1,k)) - l_arr = np.array(range(1,k)) - alpha = np.cumprod((k-l_arr)/(k+m-l_arr)) - alpha = np.insert(alpha,0,1) # include the empty product (=1) - indices = range(0,k) - Qt[:,0,indices] = np.multiply(Pt[:,0,indices],np.reshape(alpha[indices],(1,1,k))) # (21) - for p in range(0,S): - indices = range(k-z[p],k) - Qt[:,beta[p]+p*m,indices] = np.multiply(Pt[:,beta[p],indices],np.reshape(alpha[indices],(1,1,z[p]))) # (22) - - for p in range(0,S): - idx = beta[p]+p*m - Qt[:,idx+1:m+idx+1,k-1] = np.repeat(Qt[:,idx:idx+1,k-1],m,axis=1) # (23) - - # Step 4: Find remaining values of Qt_i^j - for j in range(k-1,0,-1): - for i in range(0,nb): - if Tb[i+k+m] > Tb[i+j]: - Qt[:,i+1,j-1] = Qt[:,i,j-1] + (Tb[i+k+m] - Tb[i+j])*Qt[:,i,j] #(20) with Qt replacing Pt - - return Qt[:,:,0] - __all__ = [ 'nutils', 'refinement', 'image', 'NACA', 'curve', 'smooth', 'rotation_matrix', 'sections', 'section_from_index', 'section_to_index',