Bezier Surface Patch

Ross Beveridge, December 4, 2018

This notebook sets up the linear algebra to create a bicubic Bezier Surface Patch.  Keep in mind that a Bezier Surface patch is created by essentially combining Bezier curves in one parameter, the first dimension, with a second Bezier curve in a separate parameter, the second dimension.  Always be careful to distinguish the dimensionality of an object versus the space into which it is embedded.  A Bezier Surface Patche is 2D - that is why two scalars are needed to specify location on the surface.  However, because of our formulation, it lives in 3D space (is embedded in); and the tip off is that for every pair of parameters there are three distinct coordinates. 

In [10]:
%display latex
latex.matrix_delimiters(left='|', right='|')
latex.vector_delimiters(left='[', right=']')
var('u,v')
TV = Matrix(SR, 4,1, ((u**3,u**2,u,1)))
SV = Matrix(SR, 4,1, ((v**3,v**2,v,1)))
MB = Matrix(ZZ, 4,4, ((-1,3,-3,1),(3,-6,3,0),(-3,3,0,0),(1,0,0,0)))

The Symbolic Version

Symbolic math pagackes are a blessing and curse; you can quickly build and display very large constructs including polynomials.  

First, we see here how the X, Y, and Z coordinates of the control points segragate nicely; only the X coordinates influece the bicubic term for X and similarlly for Y and Z.  Also, see that expressed as a row vector times a matrix times a matrix times a matrix times a column vector the entire equation may be written in a moderatly compact form.  

In [11]:
var('x11,x12,x13,x14,x21,x22,x23,x24,x31,x32,x33,x34,x41,x42,x43,x44')
var('y11,y12,y13,y14,y21,y22,y23,y24,y31,y32,y33,y34,y41,y42,y43,y44')
var('z11,z12,z13,z14,z21,z22,z23,z24,z31,z32,z33,z34,z41,z42,z43,z44')
GX = Matrix(SR, 4,4, (x11,x12,x13,x14,x21,x22,x23,x24,x31,x32,x33,x34,x41,x42,x43,x44))
GY = Matrix(SR, 4,4, (y11,y12,y13,y14,y21,y22,y23,y24,y31,y32,y33,y34,y41,y42,y43,y44))
GZ = Matrix(SR, 4,4, (z11,z12,z13,z14,z21,z22,z23,z24,z31,z32,z33,z34,z41,z42,z43,z44))
pretty_print(LatexExpr('Q_{x}(u,v) = '),SV.transpose(), MB.transpose(), GX, MB, TV)
pretty_print(LatexExpr('Q_{y}(u,v) = '),SV.transpose(), MB.transpose(), GY, MB, TV)
pretty_print(LatexExpr('Q_{z}(u,v) = '),SV.transpose(), MB.transpose(), GZ, MB, TV)

The next three equations are correct by nature of how SageMath is constructing them.  Quickly can them to gain a general impression for what is taking place.  Then, recognized that while symbolic math packages allow us to see equations in the full glory (gory detail) there is too much going on to develop and intuition from the algebra. 

In [12]:
QX = (SV.transpose() * MB.transpose() * GX * MB * TV)[0][0]
QY = (SV.transpose() * MB.transpose() * GY * MB * TV)[0][0]
QZ = (SV.transpose() * MB.transpose() * GZ * MB * TV)[0][0]
QX = QX.expand()
QY = QY.expand()
QZ = QZ.expand()
pretty_print(LatexExpr('Q_{x}(u,v) = '), QX)
pretty_print(LatexExpr('Q_{y}(u,v) = '), QY)
pretty_print(LatexExpr('Q_{z}(u,v) = '), QZ)

 Example 1: Surface over regular X and Y

This first example plays by very constrained rules for the X and Y values.  Essentially, but setting up regularly space samples between 0 and 3 on the X and Y axes the surface becomes a height - Z - over X and Y.  Note for this example the 16 control points are explicitly displayed in a 4 by 4 matrix in order to allow the relative position of the control points to be seen.  Then the actual 3d surface is plottted. 

In [13]:
GX = Matrix(SR, 4,4, ((0,1,2,3),(0,1,2,3),(0,1,2,3),(0,1,2,3)))
GY = Matrix(SR, 4,4, ((0,0,0,0),(1,1,1,1),(2,2,2,2),(3,3,3,3)))
GZ = Matrix(SR, 4,4, ((5,6,3,5),(3,5,6,3),(4,3,5,6),(5,4,3,5)))
QX = SV.transpose() * MB.transpose() * GX * MB * TV
QY = SV.transpose() * MB.transpose() * GY * MB * TV
QZ = SV.transpose() * MB.transpose() * GZ * MB * TV
PP = map(lambda x,y,z : Matrix(ZZ, 3,1, (x,y,z)), GX.list(), GY.list(), GZ.list())
MP = Matrix(SR, 4,4, PP); MP
Out[13]:
In [14]:
gf = (QX[0][0], QY[0][0], QZ[0][0])
def cf(u,v): return (u + v) * 0.5
fig1  = sum(map(lambda x,y,z,c : point((x,y,z),color=c,size=16), GX.list(), GY.list(), GZ.list(),rainbow(16)))
fig1 += parametric_plot3d(gf,(u, 0.0,1.0),(v, 0.0, 1.0),color=(cf,colormaps.PiYG))
fig1.show()

Example 2: Pillow

Here we begin to see the flexibility that comes with a Bezier Patch.  Just for example, realize that what you see in the next example cannot be expressed as a function of x an y.  In other words, for certain pairs (x,y) there are multiple values of z.

In [15]:
GX = Matrix(SR, 4,4, ((2,1,4,3),(2,1,4,3),(2,1,4,3),(2,1,4,3)))
GY = Matrix(SR, 4,4, ((2,1,1,2),(2,1,1,2),(3,4,4,3),(3,4,4,3)))
GZ = Matrix(SR, 4,4, ((0,1,1,0),(1,2,2,1),(1,2,2,1),(0,1,1,0)))
QX = SV.transpose() * MB.transpose() * GX * MB * TV
QY = SV.transpose() * MB.transpose() * GY * MB * TV
QZ = SV.transpose() * MB.transpose() * GZ * MB * TV
PP = map(lambda x,y,z : Matrix(SR, 3,1, (Rational(x),Rational(y),Rational(z))), GX.list(), GY.list(), GZ.list())
MP = Matrix(SR, 4,4, PP); MP
Out[15]:
In [16]:
gf = (QX[0][0], QY[0][0], QZ[0][0])
def cf(u,v): return (u + v) * 0.5
fig2  = sum(map(lambda x,y,z,c : point((x,y,z),color=c,size=16), GX.list(), GY.list(), GZ.list(),rainbow(16)))
fig2 += parametric_plot3d(gf,(u, 0.0,1.0),(v, 0.0, 1.0),color='green',opacity=0.5)
fig2.show()

Example 3: Tie the ends together

Extending the previous example, this example proves a simple point: it is possible to bring all the ends of the curve to a single point in space.

In [17]:
GX = Matrix(SR, 4,4, ((2.5,1, 4, 2.5),(2, 1, 4, 3),(2, 1, 4, 3),(2.5, 1, 4, 2.5)))
GY = Matrix(SR, 4,4, ((2.5,1, 1, 2.5),(2, 1, 1, 2),(3, 4, 4, 3),(2.5, 4, 4, 2.5)))
GZ = Matrix(SR, 4,4, ((0,  1, 1,   0),(1, 2, 2, 1),(1, 2, 2, 1),(0,   1, 1,   0)))
QX = SV.transpose() * MB.transpose() * GX * MB * TV
QY = SV.transpose() * MB.transpose() * GY * MB * TV
QZ = SV.transpose() * MB.transpose() * GZ * MB * TV
PP = map(lambda x,y,z : Matrix(SR, 3,1, (Rational(x),Rational(y),Rational(z))), GX.list(), GY.list(), GZ.list())
MP = Matrix(SR, 4,4, PP); MP
Out[17]:
In [18]:
gf = (QX[0][0], QY[0][0], QZ[0][0])
def cf(u,v): return (u + v) * 0.5
fig3  = sum(map(lambda x,y,z,c : point((x,y,z),color=c,size=16), GX.list(), GY.list(), GZ.list(),rainbow(16)))
fig3 += parametric_plot3d(gf,(u, 0.0,1.0),(v, 0.0, 1.0),color=(cf,colormaps.PiYG))
fig3.show()