<p>I was able to find control point locations for a biquartic rational B'ezier patch for a cube tiling by <a href="https://lirias.kuleuven.be/bitstream/123456789/607807/1/SphereTiling.pdf" rel="nofollow noreferrer">Dedoncker et al</a> which is based on the paper by <a href="https://collections.lib.utah.edu/dl_files/4e/77/4e7746dd53c79f8557272b92b47d2d407da4931a.pdf" rel="nofollow noreferrer">Cobb</a>. 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.</p>
<p>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 <span class="math-container">$t_u$</span> and <span class="math-container">$t_v$</span>:</p>
<p><span class="math-container">$$t_u = t_v = [0, 0, 0, 0, 0, 1, 1, 1, 1, 1]$$</span></p>
<p>Then, to obtain larger grid sizes you could use knot insertion, or plain subdivision of the Bezier patch.</p>
<p>Here I have a simple python script that creates a NURBS sphere from the control points from the paper of Dedoncker et al:</p>
<pre><code>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()
ax = fig.add_subplot(111, projection='3d')
plotSphereTile(s1, ax)
plotSphereTile(s2, ax)
plotSphereTile(s3, ax)
plotSphereTile(s4, ax)
plotSphereTile(s5, ax)
plotSphereTile(s6, ax)

plt.show()
</code></pre>
<p>Here is the resulting figure of the script:</p>
<p><a href="https://i.sstatic.net/PHWZl.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/PHWZl.png" alt="enter image description here" /></a></p>

CC BY-SA 4.0
2024-01-30T20:48:24.647Z
@ 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?
Chaosoahc
CC BY-SA 4.0
2024-02-01T08:18:15.793Z
I added a few lines to explain how to represent this as a NURBS surface.
Reynolds
CC BY-SA 4.0
2024-02-01T23:27:27.780Z
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: https://mathworld.wolfram.com/NURBSSurface.html
Chaosoahc
CC BY-SA 4.0
2024-02-01T23:36:10.807Z
(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.
Chaosoahc
CC BY-SA 4.0
2024-02-02T10:19:03.243Z
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.
Reynolds
CC BY-SA 4.0
2024-02-02T11:12:33.057Z
I added a simple python script that uses the control points from the paper and creates a NURBS sphere.
Reynolds
CC BY-SA 4.0
2024-02-14T20:05:26.340Z
Thank you for the code demonstration. I can see how your B() simplifies the NURBS basis function by assuming Bezier-style knots. However, I'm still looking for a solution using uniform knot vectors (0, 1, 2, 3, ...), not Bezier-style (0, 0, 0,... 1, 1, 1). The Bezier method does look good, but it's not ideal, as I was hoping to keep uniform knot vectors to have easier continuity when editing between adjacent patches, when I attach one of these tiles to other shapes. In the meantime I'm trying optimization algorithms to fit uniform NURBS to the surface I want, but that isn't perfect...
Chaosoahc
CC BY-SA 4.0
2024-01-29T08:03:19.833Z
2024-02-02T12:11:54.693Z
2024-02-02T12:11:54.693Z
Reynolds
<p>I would like to know how to &quot;parameterize&quot; the control points for the 6 faces of the commonly used NURBS mutli-patch sphere (&quot;Cobb sphere&quot;), for an arbitrary-size control point grid. Ideally with degree 3 and uniform knot vectors.</p>
<p>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.</p>
<p>I asked a <a href="https://math.stackexchange.com/questions/4842413/defining-a-quad-spherical-cube-tile-as-a-uniform-nurbs-surface">similar question</a> 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.</p>
<hr />
<p>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: <a href="https://i.sstatic.net/UcPkd.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/UcPkd.png" alt="enter image description here" /></a>
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 &quot;NURBS&quot;, but now when I export the CVs (although, called &quot;coefs&quot; 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.</p>

<mathematics><matrix><nurbs><spline>
How can I determine ideal NURBS surface control points for a multi-patch sphere?
Chaosoahc
# How can I determine ideal NURBS surface control points for a multi-patch sphere?

(Asked by: username_0 on 2024-01-25)

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](https://math.stackexchange.com/questions/4842413/defining-a-quad-spherical-cube-tile-as-a-uniform-nurbs-surface) 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 username_1 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: [![enter image description here](https://i.sstatic.net/UcPkd.png)](https://i.sstatic.net/UcPkd.png)
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.

**Tags:** mathematics, matrix, nurbs, spline


---

## 1 answers

### Answer by username_1 on 2024-01-29. Score: 0

I was able to find control point locations for a biquartic rational B'ezier patch for a cube tiling by [Dedoncker et al](https://lirias.kuleuven.be/bitstream/123456789/607807/1/SphereTiling.pdf) which is based on the paper by [Cobb](https://collections.lib.utah.edu/dl_files/4e/77/4e7746dd53c79f8557272b92b47d2d407da4931a.pdf). 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()
ax = fig.add_subplot(111, projection='3d')
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:

[![enter image description here](https://i.sstatic.net/PHWZl.png)](https://i.sstatic.net/PHWZl.png)

> **Comments:**
> **username_0**: @ username_1, 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?  
> **username_1**: I added a few lines to explain how to represent this as a NURBS surface.  
> **username_0**: 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: https://mathworld.wolfram.com/NURBSSurface.html  
> **username_0**: (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.  
> **username_1**: 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.  
> **username_1**: I added a simple python script that uses the control points from the paper and creates a NURBS sphere.  
> **username_0**: Thank you for the code demonstration. I can see how your B() simplifies the NURBS basis function by assuming Bezier-style knots. However, I'm still looking for a solution using uniform knot vectors (0, 1, 2, 3, ...), not Bezier-style (0, 0, 0,... 1, 1, 1). The Bezier method does look good, but it's not ideal, as I was hoping to keep uniform knot vectors to have easier continuity when editing between adjacent patches, when I attach one of these tiles to other shapes. In the meantime I'm trying optimization algorithms to fit uniform NURBS to the surface I want, but that isn't perfect...  

---

