How do Bezier Patches work in the Utah Teapot?

Question

I prematurely posted a code golf challenge to draw the Utah Teapot using this dataset (just the teapot). (Revised and Posted teapot challenge) But when I looked deeper at the data in order to whip up a little example, I realized I have no idea what's going on with that data. I have a good understanding of Bezier curves in 2D, implemented deCasteljau. But for 3D does it work the same?

Yes! It does!

The data contains patches containing 16 vertices each. Is there a standard ordering for how these are laid out? And if they correspond to the 2D curves, then the four corner points actually touch the surface and the remaining 12 are controls, right?

Yes!

My "original plan" was to simplify the shape to rectangles, project them to the canvas, and draw filled shapes in a grey computed by the magnitude of the dot product of the patch normal to a light vector. If I simplify it that far, will it even look like a teapot? Does one have to use raytracing to achieve a recognizable image?

That's subjective. :-(

While this may look like several questions, but they are all aspects of this one: "Please, kindly Guru, school me on some Bezier Patches? What do I need to know to draw the teapot?"

---

Here's the code I've written so far. (uses this matrix library: mat.ps)

%!
%%BoundingBox: 115 243 493 487
%-115 -243 translate

(mat.ps)run  %include matrix library

/tok{ token pop exch pop }def
/s{(,){search{ tok 3 1 roll }{ tok exit }ifelse }loop }def
/f(teapot)(r)file def
/patch[ f token pop { [ f 100 string readline pop s ] } repeat ]def
/vert[ f token pop { [ f 100 string readline pop s ] } repeat ]def
%vert == patch == %test data input

/I3 3 ident def      % 3D identity matrix
/Cam [ 0 0 10 ] def  % world coords of camera center viewpoint
/Theta [ 0 0 0 ] def % y-rotation x-rotation z-rotation
/Eye [ 0 0 15 ] def  % eye relative to camera vp
/Rot I3 def          % initial rotation seq

/makerot {
    Theta 0 get roty        % pan
    Theta 1 get rotx matmul % tilt
    Theta 2 get rotz matmul % twist
} def

/proj {
    Cam {sub} vop  % translate to camera coords
    Rot matmul  % perform camera rotation
    0 get aload pop Eye aload pop  % extract dot x,y,z and eye xyz
    4 3 roll div exch neg          % perform perspective projection
    4 3 roll add 1 index mul
    4 1 roll 3 1 roll sub mul exch % (ez/dz)(dx-ex) (ez/dz)(dy-ey)
} def


/R 20 def
/H -3 def
/ang 0 def
{
    300 700 translate
    1 70 dup dup scale div setlinewidth

    /Cam [ ang sin R mul  H  ang cos R mul ] def % camera revolves around Y axis at height H, dist R
    /Theta [ ang  H R atan  0 ] def % rotate camera back to origin
    /Rot makerot def % squash rotation sequence into a matrix

    patch {
        % Four corners
        %[ exch dup 0 get exch dup 3 get exch dup 12 get exch 15 get ]

        % Boundary curves
        [ exch
            dup 8 get exch dup 4 get exch dup 0 get exch %curveto4

            dup 14 get exch dup 13 get exch dup 12 get exch %curveto3

            dup 7 get exch dup 11 get exch dup 15 get exch %curveto2

            dup 1 get exch dup 2 get exch dup 3 get exch %curveto1

            dup 0 get exch %moveto
        pop ]

        { 1 sub vert exch get proj } forall

        moveto
            curveto curveto curveto curveto

        stroke
        %flushpage flush (%lineedit)(r)file pop
    } forall
    pstack
    showpage

    %exit
    /ang ang 10 add def
} loop

Here's the original Newell Teapot dataset.

And here's my spectacularly bad image:

Update: bugfix. Maybe they are laid out 'normally' after all. Selecting the correct corners at least gives a symmetrical shape:

Update: boundary curves looks better.

Answer 1

A Bi-Cubic Bezier surface patch is a 4x4 array of 3D points. Yes, the four corners touch the surface; and the rows are Bezier curves, and columns are also Bezier curves. But the deCasteljau algorithm is based on calculating the median between two points, and is equally meaningful in 3D as in 2D.

The next step in completing the above code is subdividing the patches to cover smaller portions. Then the simple boundary curve extraction above becomes a suitable polygon mesh.

Start by flattening the patches, inserting the vertex data directly instead of using a separate cache. This code iterates through the patches, looking up points in the vertex array and constructs a new array of patches which is then redefined with the same name.

 /patch[ patch{ [exch { 1 sub vert exch get }forall ] }forall ]def

Then we need the deCasteljau algorithm to split Bezier curves. vop comes from the matrix library and applies a binary operation upon corresponding elements of a vector and produces a new vector as the result.

/median { % [x0 y0 z0] [x1 y1 z1]
    {add 2 div} vop % [ (x0+x1)/2 (y0+y1)/2 (z0+z1)/2 ]
} def   

/decasteljau { % [P0] P1 P2 P3  .  P0 P1' P2' P3'  P3' P4' P5' P3
    {p3 p2 p1 p0}{exch def}forall
    /p01 p0 p1 median def
    /p12 p1 p2 median def
    /p23 p2 p3 median def
    /p012 p01 p12 median def
    /p123 p12 p23 median def
    /p0123 p012 p123 median def
    p0 p01 p012 p0123  % first half-curve
    p0123 p123 p23 p3  % second half-curve
} def   

Then some stack manipulation to apply to each row of a patch and assemble the results into 2 new patches.

/splitrows { % [b0 .. b15]  .  [c0 .. c15] [d0 .. d15] 
    aload pop % b0 .. b15
    4 {                          % on each of 4 rows
        16 12 roll decasteljau   % roll the first 4 to the top
        8 4 roll                 % exch left and right halves (probably unnecessary)
        20 4 roll                % roll new curve to below the patch (pushing earlier ones lower)
    } repeat
    16 array astore              % pack the left patch
    17 1 roll 16 array astore    % roll, pack the right patch
} def

An ugly utility lets us reuse the row code for columns. The stack comments were necessary to write this procedure, so they're probably necessary to read it. n j roll rolls n elements (to the left), j times; == the top j elements above n-th element (counting from 1). So n steady decreases, selecting where to put the element, and j selects which element to put there (dragging everything else with it). If bind were applied, this procedure would be substantially faster than a dictionary-based procedure.

% [ 0  1  2  3
%   4  5  6  7
%   8  9 10 11
%  12 13 14 15 ]
/xpose {
    aload pop % 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
    15 12 roll % 0 4 5 6 7 8 9 10 11 12 13 14 15 1 2 3
    14 11 roll % 0 4 8 9 10 11 12 13 14 15 1 2 3 5 6 7
    13 10 roll % 0 4 8 12 13 14 15 1 2 3 5 6 7 9 10 11

    12 9 roll % 0 4 8 12 1 2 3 5 6 7 9 10 11 13 14 15
    11 9 roll % 0 4 8 12 1 5 6 7 9 10 11 13 14 15 2 3
    10 8 roll % 0 4 8 12 1 5 9 10 11 13 14 15 2 3 6 7
    9 7 roll % 0 4 8 12 1 5 9 13 14 15 2 3 6 7 10 11

    8 6 roll % 0 4 8 12 1 5 9 13 2 3 6 7 10 11 14 15
    7 6 roll % 0 4 8 12 1 5 9 13 2 6 7 10 11 14 15 3
    6 5 roll % 0 4 8 12 1 5 9 13 2 6 10 11 14 15 3 7
    5 4 roll % 0 4 8 12 1 5 9 13 2 6 10 14 15 3 7 11
    4 3 roll % 0 4 8 12 1 5 9 13 2 6 10 14 3 7 11 15
    16 array astore
} def
% [ 0 4  8 12
%   1 5  9 13
%   2 6 10 14
%   3 7 11 15 ]

/splitcols {
    xpose
    splitrows
    xpose
} def

Then apply these functions to the patch data. Again, redefining patch each time.

/patch[ patch{ splitrows }forall ]def
/patch[ patch{ splitrows }forall ]def
/patch[ patch{ splitcols }forall ]def
/patch[ patch{ splitcols }forall ]def

This gives the ability to deal with smaller fragments.

Add a visibility test.

/visible { % patch  .  patch boolean 
    dup % p p
    dup 3 get exch dup 0 get exch 12 get % p p3 p0 p12
    1 index {sub} vop % p p3 p0 v0->12
    3 1 roll {sub} vop % p v0->12 v0->3
    cross /normal exch def
    dup
    [ exch dup 0 get exch dup 3 get exch dup 12 get exch 15 get ]
    { Cam {sub} vop normal dot 0 ge } forall
    %add add add 4 div 0 lt
    or or or
} def

Producing

Update: test was backwards.

Update: Test is Useless! You can see from the image that the bottom piece is not oriented outward, and of course, backface-culling doesn't prevent the handle from showing through the pot. This calls for hidden surface removal. And since Postscript doesn't have support for a Z-buffer, I guess it'll have to be a Binary Space Partition. So it's back to the books for me.

Update: Add a Model->World transform to turn the thing upright.

/Model -90 rotx def % model->world transform

/proj {
    Model matmul 0 get             % perform model->world transform
    Cam {sub} vop                  % translate to camera coords
    ...

Producing this.

Complete program so far. (uses matrix library:mat.ps.) In ghostscript, you can view an animated rotation by holding [enter].

    %!
    %%BoundingBox: 109 246 492 487
    %-109 -246 translate

    (mat.ps)run  %include matrix library
    (det.ps)run  %supplementary determinant function

    /tok{ token pop exch pop }def
    /s{(,){search{ tok 3 1 roll }{ tok exit }ifelse }loop }def
    /f(teapot)(r)file def
    /patch[ f token pop { [ f 100 string readline pop s ] } repeat ]def
    /vert[ f token pop { [ f 100 string readline pop s ] } repeat ]def
    /patch[ patch{ [exch { 1 sub vert exch get }forall ] }forall ]def
    %vert == patch == %test data input flush quit

    /I3 3 ident def      % 3D identity matrix
    /Cam [ 0 0 10 ] def  % world coords of camera center viewpoint
    /Theta [ 0 0 0 ] def % y-rotation x-rotation z-rotation
    /Eye [ 0 0 15 ] def  % eye relative to camera vp
    /Rot I3 def          % initial rotation seq

    /Model -90 rotx def % model->world transform

    /makerot {
            Theta 0 get roty        % pan
            Theta 1 get rotx matmul % tilt
            Theta 2 get rotz matmul % twist
    } def

    /proj {
            Model matmul 0 get             % perform model->world transform
            Cam {sub} vop  % translate to camera coords
            Rot matmul  % perform camera rotation
            0 get aload pop Eye aload pop  % extract dot x,y,z and eye xyz
            4 3 roll div exch neg          % perform perspective projection
            4 3 roll add 1 index mul
            4 1 roll 3 1 roll sub mul exch % (ez/dz)(dx-ex) (ez/dz)(dy-ey)
    } def

    /median { % [x0 y0 z0] [x1 y1 z1]
            {add 2 div} vop % [ (x0+x1)/2 (y0+y1)/2 (z0+z1)/2 ]
    } def

    /decasteljau { % [P0] P1 P2 P3  .  P0 P1' P2' P3'  P3' P4' P5' P3
            {p3 p2 p1 p0}{exch def}forall
            /p01 p0 p1 median def
            /p12 p1 p2 median def
            /p23 p2 p3 median def
            /p012 p01 p12 median def
            /p123 p12 p23 median def
            /p0123 p012 p123 median def
            p0 p01 p012 p0123
            p0123 p123 p23 p3
    } def

    /splitrows { % [b0 .. b15]  .  [c0 .. c15] [d0 .. d15]
            aload pop % b0 .. b15
            4 {
                    16 12 roll decasteljau
                    %8 4 roll
                    20 4 roll
            } repeat
            16 array astore
            17 1 roll 16 array astore
    } def

    /xpose {
            aload pop % 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
            15 12 roll % 0 4 5 6 7 8 9 10 11 12 13 14 15 1 2 3
            14 11 roll % 0 4 8 9 10 11 12 13 14 15 1 2 3 5 6 7
            13 10 roll % 0 4 8 12 13 14 15 1 2 3 5 6 7 9 10 11

            12 9 roll % 0 4 8 12 1 2 3 5 6 7 9 10 11 13 14 15
            11 9 roll % 0 4 8 12 1 5 6 7 9 10 11 13 14 15 2 3
            10 8 roll % 0 4 8 12 1 5 9 10 11 13 14 15 2 3 6 7
            9 7 roll % 0 4 8 12 1 5 9 13 14 15 2 3 6 7 10 11

            8 6 roll % 0 4 8 12 1 5 9 13 2 3 6 7 10 11 14 15
            7 6 roll % 0 4 8 12 1 5 9 13 2 6 7 10 11 14 15 3
            6 5 roll % 0 4 8 12 1 5 9 13 2 6 10 11 14 15 3 7
            5 4 roll % 0 4 8 12 1 5 9 13 2 6 10 14 15 3 7 11
            4 3 roll % 0 4 8 12 1 5 9 13 2 6 10 14 3 7 11 14
            16 array astore
    } def

    /splitcols {
            xpose
            splitrows
            xpose exch xpose
    } def

    /patch[ patch{ splitrows }forall ]def
    /patch[ patch{ splitrows }forall ]def
    /patch[ patch{ splitrows }forall ]def
    /patch[ patch{ splitrows }forall ]def
    /patch[ patch{ splitcols }forall ]def
    /patch[ patch{ splitcols }forall ]def

    /color {normal light dot 1 add 4 div
%1 exch sub
setgray} def

    /visible { % patch  .  patch boolean
            dup % p p
            dup 3 get exch dup 0 get exch 12 get % p p3 p0 p12
            1 index {sub} vop % p p3 p0 v0->12
            3 1 roll {sub} vop % p v0->12 v0->3
            cross /normal exch def
            dup
            [ exch dup 0 get exch dup 3 get exch dup 12 get exch 15 get ]
            { Cam {sub} vop normal dot 0 ge } forall
            %add add add 4 div 0 lt
            or or or
    } def

    /drawpatch {

            % Four corners
            %[ exch dup 0 get exch dup 3 get exch dup 12 get exch 15 get ]

            visible {

                    [ exch
                            % control rows
                            %dup 4 get exch dup 5 get exch dup 6 get exch dup 7 get exch
                            %dup 11 get exch dup 10 get exch dup 9 get exch dup 8 get exch

                            % control columns
                            %dup 1 get exch dup 5 get exch dup 9 get exch dup 13 get exch
                            %dup 14 get exch dup 10 get exch dup 6 get exch dup 2 get exch

                            % Boundary curves
                            dup 8 get exch dup 4 get exch dup 0 get exch %curveto4
                            dup 14 get exch dup 13 get exch dup 12 get exch %curveto3
                            dup 7 get exch dup 11 get exch dup 15 get exch %curveto2
                            dup 1 get exch dup 2 get exch dup 3 get exch %curveto1
                            dup 0 get exch %moveto
                    pop ]

                    { proj } forall

                    moveto curveto curveto curveto curveto
                    %moveto lineto lineto lineto lineto lineto lineto lineto closepath
                    %moveto lineto lineto lineto lineto lineto lineto lineto closepath

                    stroke
                    %flushpage flush (%lineedit)(r)file pop

            }{
                    pop
            }ifelse

    } def

    /R 20 def
    /H -3 def
    /ang 10 def
    {
            300 700 translate
            1 70 dup dup scale div setlinewidth

            % camera revolves around Y axis at height H, dist R
            /Cam [ ang sin R mul  H  ang cos R mul ] def
            /Theta [ ang  H R atan  0 ] def   % rotate camera back to origin
            /Rot makerot def        % squash rotation sequence into a matrix

            patch {
                    drawpatch
            } forall
            pstack
            showpage

            %exit
            /ang ang 10 add def
    } loop

Answer 2

Based on assistance on math.StackExchange, I've been led to a sub-goal of supplementing the matrix library with a function to compute determinants.

So, this code passes some clumsy initial tests, but it's damned ugly, I must admit:

GS>[[1 0][0 1]] det
GS<1>=
1
GS>[[0 1][1 0]] det =
-1
GS>(mat.ps) run
GS>3 ident
GS<1>det =
1
GS>[[1 2 3][4 5 6][7 8 9]] det =
0
GS>

Update. A little more readable.

Update. A lot more readable using dot and cross. Thanks again, MvG.

(mat.ps) run % use dot and cross from matrix library

/elem { % M i j
    3 1 roll get exch get % M_i_j
} def

/det {
    dup length 1 index 0 get length ne { /det cvx /typecheck signalerror } if
    1 dict begin /M exch def
        M length 2 eq {
            M 0 0 elem
            M 1 1 elem mul
            M 0 1 elem
            M 1 0 elem mul sub
        }{
            M length 3 eq {
                M aload pop cross dot
            }{ /det cvx /rangecheck signalerror } ifelse
        } ifelse
    end
} def