# Implementation of tensor product formula

For the Bézier surface, which owns the following matrix definition:

\begin{align*} \mathbf S(u,v)&=\sum_{i=0}^m \sum_{j=0}^n \mathbf P_{i,j} B_{m,i}(u) B_{n,j}(v)\\ &=\small \begin{pmatrix}B_{m,0}(u)&\cdots&B_{m,m}(u)\end{pmatrix}_{1\times(m+1)} \begin{pmatrix} \mathbf P_{0,0}&\cdots&\mathbf P_{0,n}\\ P_{1,0}&\cdots&\mathbf P_{1,n}\\\vdots&\vdots&\vdots\\ \mathbf P_{m,0}&\cdots&\mathbf P_{m,n} \end{pmatrix} \begin{pmatrix} B_{n,0}(v)\\B_{n,1}(v)\\ \vdots\\ B_{n,n}(v) \end{pmatrix}_{(n+1)\times 1} \end{align*}

where, $B_{n,i}(u)$ is Bernstein basis.

vec1 = {B[0],...,B[m]};
mat = {{P[0,0],...P[0,n]},...,{P[m,0],...P[m,n]}};
vec2 = {B[0],...,B[n]};

bez = vec1.mat1.vec2


However, the $P_{i,j}$ is the coordinate of a 3D point, which own this style:{x,y,z}. So I cannot use vec1.mat1.vec2 directly.

An alternative method is using Hold[] to unevaluate the coordinate {x,y,z}. Namely, Hold[{x,y,z}]. Lastly, with the help of ReleaseHold[] to evaluate the expression.

vec1 = {B[0],...,B[m]};
mat = Map[Hold,{{P[0,0],...P[0,n]},...,{P[m,0],...P[m,n]}},{2}];
vec2 = {B[0],...,B[n]};

bez = vec1.mat1.vec2 // ReleaseHold


Another way that I came up with is

vec1.mat[[All, All, #]].vec2 & /@ {1, 2, 3}


### Comparison

Bernstein[0, 0, u_?NumericQ] := 1
Bernstein[n_, i_, u_?NumericQ] := Binomial[n, i] u^i (1 - u)^(n - i)

BezierSurface2[pts_, u_?NumericQ, v_?NumericQ] :=
Module[{m, n, AllBasis},
{m, n} = Dimensions[pts, 2];
AllBasis =
Function[{deg, u0}, Bernstein[deg, #, u0] & /@ Range[0, deg]];
With[{row = AllBasis[m - 1, u], col = AllBasis[n - 1, v]},
row.Map[Hold, pts, {2}].col // ReleaseHold]
]

BezierSurface1[pts_, u_?NumericQ, v_?NumericQ] :=
Module[{m, n, AllBasis},
{m, n} = Dimensions[pts, 2];
AllBasis =
Function[{deg, u0}, Bernstein[deg, #, u0] & /@ Range[0, deg]];
With[{row = AllBasis[m - 1, u], col = AllBasis[n - 1, v]},
row.pts[[All, All, #]].col & /@ {1, 2, 3}]
]


cpts = Table[{i, j, RandomReal[{-1, 1}]}, {i, 5}, {j, 5}];
ParametricPlot3D[
BezierSurface1[cpts, u, v], {u, 0, 1}, {v, 0, 1}] // AbsoluteTiming

ParametricPlot3D[
BezierSurface2[cpts, u, v], {u, 0, 1}, {v, 0, 1}] // AbsoluteTiming

f = BezierFunction[cpts];
ParametricPlot3D[f[u, v], {u, 0, 1}, {v, 0, 1}] // AbsoluteTiming


So my question: is there a more efficient method to implement this formula?

• I accepted the J.M.'s answer because it is a general strategy to define CAGDBezierFunction[]. In addition, I will award the bounty to xzczd owns to that it is the first answer and his Transpose[] and Listable method is very useful for me.
– xyz
Commented Feb 21, 2016 at 4:45

n = 4;
BlockRandom[SeedRandom[42, Method -> "Legacy"]; (* for reproducibility *)
cpts = Table[{i, j, RandomReal[{-1, 1}]}, {i, n + 1}, {j, n + 1}]];

GraphicsRow[{ParametricPlot3D[BezierFunction[cpts][u, v], {u, 0, 1}, {v, 0, 1},
Evaluated -> True], (* built-in function *)
ParametricPlot3D[Evaluate[Fold[#2.#1 &, cpts, (* using dot-products *)
{BernsteinBasis[n, Range[0, n], u],
BernsteinBasis[n, Range[0, n], v]}]],
{u, 0, 1}, {v, 0, 1}]}]


You can of course replace BernsteinBasis[] with your own Bernstein[]; no need for Hold[] trickery!

• Thanks a lot(^_^).Nice to use the BernsteinBasis[] listable attribute and Fold[#1.#2 &, pts, {row,col}] strategy.
– xyz
Commented Feb 20, 2016 at 9:36
• @Shutao, in this simple case, you can of course just do col.(row.pts); the Fold[] approach becomes useful when you start considering $n$-dimensional Bézier hypersurfaces. Commented Feb 20, 2016 at 9:55
• @J.M. I have implemented the Bezier curve and surface in my CAGD package via The NURBS Book. Did you know the definition of $n$-dimentional Bezier hypersurface？The built-in BezierFunction[] owns the below usage：f=BezierFunction[RandomReal[1,{5,5,5,1}]].However,I didn't know its definition.
– xyz
Commented Feb 20, 2016 at 10:29
• @Shutao, well, if you have, say, the 4D case, with a set of control points represented by a $p\times q\times r\times 4$ list, you just plug this list into the second argument of Fold[], and then the third argument is something like {BernsteinBasis[p - 1, Range[0, p - 1], u], BernsteinBasis[q - 1, Range[0, q - 1], v], BernsteinBasis[r - 1, Range[0, r - 1], w]}. Commented Feb 20, 2016 at 10:46
• @J.M. Thanks a lot. Got it. For the curve case, the control-points owns dimention：$n \times i$, where $i=2,3,4$. For the surface case, the control nets owns dimention :$m \times n \times i$.
– xyz
Commented Feb 20, 2016 at 11:04

## Update

Fully compiling the code to C makes it as fast as the built-in:

cBernstein =
Compile @@ (Hold[{n, {j, _Real, 1}, u}, Table[expr u^i (1 - u)^(n - i), {i, j}],
CompilationTarget -> C] /. expr -> FunctionExpand[Binomial[n, i]])

BezierSurface4 =
With[{cBernstein = cBernstein},
With[{AllBasis = Function[{deg, u0}, cBernstein[deg, Range[0, deg], u0]]},
Compile[{{pts, _Real, 3}, u, v}, Module[{m, n}, {m, n} = Dimensions[pts, 2];

With[{row = AllBasis[m - 1, u], col = AllBasis[n - 1, v]},
row.Transpose[pts, {1, 3, 2}].col]],
RuntimeOptions -> {"EvaluateSymbolically" -> False},
CompilationOptions -> {"InlineCompiledFunctions" -> True},
CompilationTarget -> C]]];

f = BezierFunction[cpts];
ParametricPlot3D[f[u, v], {u, 0, 1}, {v, 0, 1}]; // AbsoluteTiming
(* {0.018014, Null} *)

ParametricPlot3D[BezierSurface4[cpts, u, v], {u, 0, 1}, {v, 0, 1}]; // AbsoluteTiming

(* {0.018006, Null} *)


Compile and Transpose and Listable will help:

cBernstein =
Compile[{n, {i, _Real, 1}, u},
Evaluate[FunctionExpand[Binomial[n, i]] u^i (1 - u)^(n - i)]]

BezierSurface3[pts_, u_?NumericQ, v_?NumericQ] :=
Module[{m, n, AllBasis}, {m, n} = Dimensions[pts, 2];
AllBasis = Function[{deg, u0}, cBernstein[deg, Range[0, deg], u0]];
With[{row = AllBasis[m - 1, u], col = AllBasis[n - 1, v]},
row.Transpose[pts, {1, 3, 2}].col]]

ParametricPlot3D[BezierSurface1[cpts, u, v], {u, 0, 1}, {v, 0, 1}]; // AbsoluteTiming
(* {0.219170, Null} *)
ParametricPlot3D[BezierSurface3[cpts, u, v], {u, 0, 1}, {v, 0, 1}]; // AbsoluteTiming
(* {0.098571, Null} *)

• Nice solution. However, the OP did not make it clear why he wanted to re-implement a built-in function. Yes, I can read what he asked for—but he did not explain why. Since your code is only as fast as the built-in code (not surprisingly—it would be really impressive if it was faster), the only reason he could have would be to implement something closely related that is not built-in. Commented Feb 9, 2016 at 3:51
• As I show in my answer, you don't even need to transpose; you just have to perform the multiplications in a different order. Commented Feb 20, 2016 at 6:56
• @J.M. Oh, you're right, row.Transpose[pts, {1, 3, 2}].col can be further simplified to col.(row.pts) in this case, the speed improvement is negligible though. (Transpose is so fast. ) Commented Feb 20, 2016 at 7:49