list manipulation - 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?

Answer

This might help you get an idea:

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!