8
$\begingroup$

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

plots

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

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – xyz
    Feb 21 '16 at 4:45
12
+200
$\begingroup$

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}]}]

can you spot any difference?

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

$\endgroup$
10
  • $\begingroup$ Thanks a lot(^_^).Nice to use the BernsteinBasis[] listable attribute and Fold[#1.#2 &, pts, {row,col}] strategy. $\endgroup$
    – xyz
    Feb 20 '16 at 9:36
  • $\begingroup$ @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. $\endgroup$
    – J. M.'s torpor
    Feb 20 '16 at 9:55
  • $\begingroup$ @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. $\endgroup$
    – xyz
    Feb 20 '16 at 10:29
  • 1
    $\begingroup$ @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]}. $\endgroup$
    – J. M.'s torpor
    Feb 20 '16 at 10:46
  • $\begingroup$ @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$. $\endgroup$
    – xyz
    Feb 20 '16 at 11:04
11
+50
$\begingroup$

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} *)

Old Answer

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} *)
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – TheDoctor
    Feb 9 '16 at 3:51
  • 1
    $\begingroup$ As I show in my answer, you don't even need to transpose; you just have to perform the multiplications in a different order. $\endgroup$
    – J. M.'s torpor
    Feb 20 '16 at 6:56
  • $\begingroup$ @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. ) $\endgroup$
    – xzczd
    Feb 20 '16 at 7:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.