Iss14

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

splipy/SplineObject.py

@@ -8,6 +8,7 @@
 from splipy import BSplineBasis
 from splipy.utils import reshape, rotation_matrix, is_singleton, ensure_listlike,\
                          check_direction, ensure_flatlist, check_section, sections
+import splipy.state as state
 
 __all__ = ['SplineObject']
 
@@ -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,40 @@ def raise_order(self, *raises):
             raise ValueError("Cannot lower order using raise_order")
         if all(r == 0 for r in raises):
             return
+
+        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
+        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
@@ -1394,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] 
+

splipy/utils/__init__.py

@@ -171,5 +171,5 @@ def uniquify(iterator):
     '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'
 ]