graphics - pierre - lineas curvas de bezier

¿Cómo funcionan los parches de Bezier en la Tetera de Utah?

Question

Publiqué prematuramente un desafío de golf de código para dibujar la Tetera de Utah con este conjunto de datos ( solo la tetera ). ( Desafío tetera revisado y publicado ) Pero cuando miré más profundamente en los datos con el fin de dar un pequeño ejemplo, me di cuenta de que no tengo idea de qué está pasando con esa información. Tengo una buena comprensión de las curvas de Bezier en 2D, implementado deCasteljau. Pero para 3D, ¿funciona igual?

¡Sí! ¡Lo hace!

Los datos contienen parches que contienen 16 vértices cada uno. ¿Hay un pedido estándar de cómo se presentan? Y si corresponden a las curvas 2D, entonces los cuatro puntos de esquina realmente tocan la superficie y los 12 restantes son controles, ¿verdad?

¡Sí!

Mi "plan original" fue simplificar la forma de los rectángulos, proyectarlos en el lienzo y dibujar formas rellenas en un gris calculado por la magnitud del producto de puntos del parche normal a un vector de luz. Si lo simplifico tan lejos, ¿se verá siquiera como una tetera? ¿Se debe usar raytracing para lograr una imagen reconocible?

Eso es subjetivo. :-(

Si bien esto puede parecer varias preguntas, pero todos son aspectos de este: "Por favor, gentil Gurú, enséñame algunos parches de Bezier? ¿Qué necesito saber para dibujar la tetera?"

Aquí está el código que he escrito hasta ahora. (usa esta biblioteca matriz: 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

Aquí está el conjunto de datos original de Newell Teapot .

Y aquí está mi imagen espectacularmente mala:

Actualización: corrección de errores. Quizás después de todo estén distribuidos "normalmente". Seleccionar las esquinas correctas al menos da una forma simétrica:

Actualización: las curvas de límite se ven mejor.

Answer 1

Un parche de superficie Bicibic Bezier es una matriz 4x4 de puntos 3D. Sí, las cuatro esquinas tocan la superficie; y las filas son curvas de Bezier, y las columnas también son curvas de Bezier. Pero el algoritmo deCasteljau se basa en el cálculo de la mediana entre dos puntos, y es igualmente significativo en 3D como en 2D.

El siguiente paso para completar el código anterior es subdividir los parches para cubrir porciones más pequeñas. Luego, la extracción de curva de límite simple anterior se convierte en una malla de polígono adecuada.

Comience aplanando los parches, insertando los datos de vértice directamente en lugar de usar un caché separado. Este código itera a través de los parches, busca puntos en la matriz de vértices y construye una nueva matriz de parches que luego se redefine con el mismo nombre.

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

Entonces necesitamos el algoritmo deCasteljau para dividir las curvas de Bezier. vop proviene de la biblioteca de matrices y aplica una operación binaria sobre los elementos correspondientes de un vector y produce un nuevo vector como resultado.

/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   

Luego, aplique un poco de manipulación a cada fila de un parche y ensamble los resultados en 2 nuevos parches.

/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

Una fea utilidad nos permite reutilizar el código de fila para las columnas. Los comentarios de la pila fueron necesarios para escribir este procedimiento, por lo que probablemente sean necesarios para leerlo. nj roll rolls n elementos (a la izquierda), j veces; == los elementos superiores j sobre el elemento n-ésimo (contando a partir de 1). Entonces n disminuye constantemente, seleccionando dónde colocar el elemento, yj selecciona qué elemento colocar allí (arrastrando todo lo demás con él). Si se aplicara bind , este procedimiento sería sustancialmente más rápido que un procedimiento basado en diccionario.

% [ 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

Luego, aplique estas funciones a los datos del parche. De nuevo, redefiniendo el parche cada vez.

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

Esto le da la capacidad de tratar con fragmentos más pequeños.

Agregue una prueba de visibilidad.

/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

Productor

Actualización: la prueba fue al revés.

Actualización: ¡la prueba es inútil! Puede ver por la imagen que la pieza inferior no está orientada hacia afuera, y por supuesto, la eliminación de cara posterior no impide que la manija se muestre a través del bote. Esto requiere eliminación de superficie oculta. Y dado que Postscript no tiene soporte para un Z-buffer, supongo que tendrá que ser una Partición espacial binaria. Así que volví a los libros para mí.

Actualización: agregue una transformación Model-> World para hacer que la cosa se levante.

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

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

Produciendo esto

Completa el programa hasta el momento. (usa la biblioteca de matriz: mat.ps. ) En ghostscript, puede ver una rotación animada manteniendo presionado [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

Sobre la base de la asistencia en math.StackExchange , me han llevado a una sub-meta de complementar la biblioteca de la matriz con una función para calcular los determinantes.

Entonces, este código pasa algunas pruebas iniciales torpes, pero es condenadamente feo , debo admitir:

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>

Actualizar. Un poco más legible

Actualizar. Mucho más legible usando punto y cruz. Gracias de nuevo, 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