# How can I determine ideal NURBS surface control points for a multi-patch sphere?

I would like to know how to "parameterize" the control points for the 6 faces of the commonly used NURBS mutli-patch sphere ("Cobb sphere"), for an arbitrary-size control point grid. Ideally with degree 3 and uniform knot vectors.

For example, given a 6x6 grid, where would I place the control points and how would I weight them in order to achieve perfect curvature? Then, extend this to a 5x5 grid, 7x7, 10x10, etc. I have seen this done with 2D circles/arcs, and I know rational/non-one weights are probably needed, but I can't find a good method.

I asked a similar question in Mathematics Exchange, but it was perhaps too verbose or better suited for here, and did not receive any feedback. Thanks in advance for any advice or direction someone could provide.

Edit: I tried out the points suggested by user Reynolds from the Dedoncker et al paper, but the tiles don't perfectly fit a sphere. Adding an image showing some of my evaluations of those points in a few different platforms: The paper mentions the CVs were implemented in Octave (albeit, scaled or otherwise modified), and indeed I had actually seen them before in Octave during previous searches. They're called "NURBS", but now when I export the CVs (although, called "coefs" here...) and evaluate them with my Python NURBS functions, they do not match the Octave plot. Maybe they are actually Bézier? I believe I need to convert these from Bézier to NURBS but I don't totally understand that process, and ultimately I need to work with NURBS.

I was able to find control point locations for a biquartic rational B'ezier patch for a cube tiling by Dedoncker et al which is based on the paper by Cobb. The control point values can be found in table A.2. One important thing to notice is that the control points are homogeneous. The control positions of the other faces can be determined by symmetry.

To represent this patch as a NURBS surface remember that a Bezier curve is a special type of NURBS curve with a knot vector with multiple knots. So use the knot vectors for a quartic NURBS piece $$t_u$$ and $$t_v$$:

$$t_u = t_v = [0, 0, 0, 0, 0, 1, 1, 1, 1, 1]$$

Then, to obtain larger grid sizes you could use knot insertion, or plain subdivision of the Bezier patch.

Here I have a simple python script that creates a NURBS sphere from the control points from the paper of Dedoncker et al:

import matplotlib.pyplot as plt
from scipy.linalg import expm, norm
import numpy as np
import math

cps = np.zeros([5, 5, 4])
cps[0][0] = np.array([0.577350269189626, 0.577350269189626, -0.577350269189626, 1])
cps[1][0] = np.array([0.278838767912603, 0.632392158505876, -0.632392158505876, 0.891211203608397])
cps[2][0] = np.array([0, 0.647791890991355, -0.647791890991355, 0.859116756396542])
cps[3][0] = np.array([-0.278838767912603, 0.632392158505876, -0.632392158505876, 0.891211203608397])
cps[4][0] = np.array([-0.577350269189626, 0.577350269189626, -0.577350269189626, 1])

cps[0][1] = np.array([0.632392158505876, 0.278838767912603, -0.632392158505876, 0.891211203608397])
cps[1][1] = np.array([0.315090742770461, 0.315090742770461, -0.762259526419164, 0.762259526419164])
cps[2][1] = np.array([0, 0.328648516366383, -0.804938188574224, 0.71866517354005])
cps[3][1] = np.array([-0.315090742770461, 0.315090742770461, -0.762259526419164, 0.762259526419164])
cps[4][1] = np.array([-0.632392158505876, 0.278838767912603, -0.632392158505876, 0.891211203608397])

cps[0][2] = np.array([0.647791890991355, 0, -0.647791890991355, 0.859116756396542])
cps[1][2] = np.array([0.328648516366383, 0, -0.804938188574224, 0.71866517354005])
cps[2][2] = np.array([0, 0, -0.859116756396542, 0.671272431591931])
cps[3][2] = np.array([-0.328648516366383, 0, -0.804938188574224, 0.71866517354005])
cps[4][2] = np.array([-0.647791890991355, 0, -0.647791890991355, 0.859116756396542])

cps[0][3] = np.array([0.632392158505876, -0.278838767912603, -0.632392158505876, 0.891211203608397])
cps[1][3] = np.array([0.315090742770461, -0.315090742770461, -0.762259526419164, 0.762259526419164])
cps[2][3] = np.array([0, -0.328648516366383, -0.804938188574224, 0.71866517354005])
cps[3][3] = np.array([-0.315090742770461, -0.315090742770461, -0.762259526419164, 0.762259526419164])
cps[4][3] = np.array([-0.632392158505876, -0.278838767912603, -0.632392158505876, 0.891211203608397])

cps[0][4] = np.array([0.577350269189626, -0.577350269189626, -0.577350269189626, 1])
cps[1][4] = np.array([0.278838767912603, -0.632392158505876, -0.632392158505876, 0.891211203608397])
cps[2][4] = np.array([0, -0.647791890991355, -0.647791890991355, 0.859116756396542])
cps[3][4] = np.array([-0.278838767912603, -0.632392158505876, -0.632392158505876, 0.891211203608397])
cps[4][4] = np.array([-0.577350269189626, -0.577350269189626, -0.577350269189626, 1])

#rotation matrix to get other patch cps
def M(axis, theta):
M_3 = expm(np.cross(np.eye(3), axis/norm(axis)*theta))
M_4 = np.eye(4)
M_4[0:3, 0:3] = M_3
return M_4

Y_0 = M(np.array([0, 1, 0]), math.pi / 2)
Y_1 = M(np.array([0, 1, 0]), math.pi)
Y_2 = M(np.array([0, 1, 0]), 1.5 * math.pi)

Z_0 = M(np.array([0, 0, 1]), math.pi / 2)
Z_1 = M(np.array([0, 0, 1]), 1.5 * math.pi)

cps2 = np.dot(cps, Y_0)
cps3 = np.dot(cps, Y_1)
cps4 = np.dot(cps, Y_2)

cps5 = np.dot(cps2, Z_0)
cps6 = np.dot(cps2, Z_1)

#Bezier basis functions
def B_0(t):
return (1-t)**4

def B_1(t):
return 4*(1-t)**3 * t

def B_2(t):
return 6*(1-t)**2 * t**2

def B_3(t):
return 4*(1-t) * t**3

def B_4(t):
return t**4

def B(i, t):
if(i == 0):
return B_0(t)
elif(i == 1):
return B_1(t)
elif(i == 2):
return B_2(t)
elif(i == 3):
return B_3(t)
elif(i == 4):
return B_4(t)

def surface(u, v, cps):
S = np.array([0, 0, 0, 0])
W = 0
for i in range(0, 5):
for j in range(0, 5):
w = cps[i][j][3] * B(i, u) * B(j, v)
S = S + w * cps[i][j]
W = W + w

h = S / W
return h

def spheretile(cps, N):
U = np.linspace(0, 1, N)
V = np.linspace(0, 1, N)

spheretile = np.zeros([N, N, 4])
i = 0
j = 0
for u in U:
j = 0
for v in V:
spheretile[i][j] = surface(u, v, cps)
#division from homogeneous coordinates happens here
spheretile[i][j] = spheretile[i][j] / spheretile[i][j][3]

j = j+ 1
i = i + 1

return spheretile

s1 = spheretile(cps, 20)
s2 = spheretile(cps2, 20)
s3 = spheretile(cps3, 20)
s4 = spheretile(cps4, 20)
s5 = spheretile(cps5, 20)
s6 = spheretile(cps6, 20)

def plotSphereTile(s, ax):
Xs = s[:, :, 0]
Ys = s[:, :, 1]
Zs = s[:, :, 2]
ax.plot_surface(Xs, Ys, Zs)

fig = plt.figure()
plotSphereTile(s1, ax)
plotSphereTile(s2, ax)
plotSphereTile(s3, ax)
plotSphereTile(s4, ax)
plotSphereTile(s5, ax)
plotSphereTile(s6, ax)

plt.show()


Here is the resulting figure of the script:

• @ Reynolds, thanks for the reference to those papers. I had seen the Cobb paper before but couldn't work with the points as they were Bézier. I think I'm having a similar issue with the second paper, as seen in the image I added to my question in an edit just now. Do you have suggestions for converting these from Bézier to NURBS? Jan 30 at 20:48
• I added a few lines to explain how to represent this as a NURBS surface. Feb 1 at 8:18
• Thanks, I was unsure about the knot vector but it does seem that I was using that same multiple knots sequence (see plot titles in the main question post). In response to your other comment, I am multiplying by the weights in the numerator and dividing by the sum of the basis functions, although I'm not sure what you mean by "factored by their respective w"- I am following this formula: mathworld.wolfram.com/NURBSSurface.html Feb 1 at 23:27
• (and presumably Blender is using that general form too). However, those are the NURBS basis functions, not Bezier- do they use different basis functions, despite NURBS being a generalization of Bezier? I'm seeing details about Bernstein polynomials but I can't tell if that's equivalent to what I already have for NURBS. Feb 1 at 23:36
• Okay, I think what you might be missing here is that evaluating the NURBS surface with homogeneous points will leave you with a homogeneous surface. For each point generated S(u,v) = [x,y,z,w] do a perspective division [x/w,y/w,z/w] to obtain the 3D point. Feb 2 at 10:19