Creating Bezier surfaces using procedural graphics

Question

Are there any package for LaTeX that can create Bezier surfaces? I mean a actual existing function that takes matrices as arguments and generates the surface, not a function written by myself.

Answer 1

Yes: pgfplots can do it by means of its patchplots library.

pgfplots supports mesh, surf, and patch plots, each with different shader configurations and with different surface types (called patch type in pgfplots). The color at each vertex is typically determined from scalar "color data" combined with a color map (in this document, the "color data" uses the z value at each point and the colormap is the system default).

I suppose the patch plots are adequate here. These accept a series of patches, either by means of adjacency matrizes or simply by means of consecutive vertices which form the surfaces. The points can be read from tables or simply provided as coordinates in round braces.

Here is a survey over the supported patch types if you choose coordinates in round braces and consecutive vertices to define the surfaces. I wrote down just one surface per patch type, but any number is accepted.

\documentclass[a4paper]{article}

\usepackage{pgfplots}
\usepgfplotslibrary{patchplots}

\begin{document}
\thispagestyle{empty}

\pgfplotsset{
    nodes near coords=\coordindex,
    every node near coord/.style={circle,anchor=center,font=\tiny,fill=white,draw=black},
}

    \begin{tikzpicture}
    \begin{axis}[title=Rectangle from patch input]
    \addplot3[patch,shader=interp,patch type=rectangle] coordinates {
        (0,0,1) 
        (1,0,0) 
        (1,1,0) 
        (0,1,0) 
    };
    \end{axis}
    \end{tikzpicture}
    %
    \begin{tikzpicture}
    \begin{axis}[
        title=Bilinear from $4$--point patch input]
    \addplot3[patch,shader=interp,patch type=bilinear]
    coordinates {
        (0,0,1) 
        (1,0,0) 
        (1,1,0) 
        (0,1,0) 
    };
    \end{axis}
    \end{tikzpicture}

    \begin{tikzpicture}
    \begin{axis}[enlargelimits,
        title=Single Triangle patch]
    \addplot3[patch,shader=interp,] coordinates {
        (0,0,1) 
        (1,0,0) 
        (1,1,0) 
    };
    \end{axis}
    \end{tikzpicture}
    %
    \begin{tikzpicture}
    \begin{axis}[
        title=Quadratic Triangle]
    \addplot[patch,shader=interp,patch type=triangle quadr,
        ]
    coordinates {
        (0,0)  (5,4)  (0,7) 
        (2,3)  (3,6)  (-1,4) 
    };
    \end{axis}
    \end{tikzpicture}

    \begin{tikzpicture}
    \begin{axis}[
        title=Single Biquadratic Quadrilateral]
    \addplot[patch,shader=interp,patch type=biquadratic,
        ]
    coordinates {
        (0,0)  (6,1)  (5,5)  (-1,5) 
        (3,1)  (6,3)  (2,6)  (0,3) 
        (3,3.75) 
    };
    \end{axis}
    \end{tikzpicture}
    %
    %
    \begin{tikzpicture}
    \begin{axis}[
        title=Single Bicubic Quadrilateral]
    \addplot3[patch,shader=interp,patch type=bicubic,]
    coordinates {
        (0,0,1)  (1,0,0)  (2,0,0)  (3,0,0) 
        (0,1,0)  (1,1,0)  (2,1,0)  (3,1,0)
        (0,2,0)  (1,2,0)  (2,2,0)  (3,2,0)
        (0,3,0)  (1,3,0)  (2,3,0)  (3,3,0) 
    };
    \end{axis}
    \end{tikzpicture}

\end{document}

The nodes near coords instruction on top activates the small nodes with indices to illustrate in which sequence pgfplots expects the vertices. The sequence depends on the patch type.

There are also two "raw" patch types which allow to input coons patches and tensor bezier directly in their bezier basis representation (i.e. with control points rather than interpolation points). However, these are considered to be low-level and do not really support suitable z ordering.

Aside from explicitly defined patches, you can also sample a math expression using patch type sampling:

\documentclass[a4paper]{article}

\usepackage{pgfplots}
\usepgfplotslibrary{patchplots}

\begin{document}
\thispagestyle{empty}

\begin{tikzpicture}
    \begin{axis}
    \addplot3[
        patch,
        patch type sampling,
        patch type=bicubic,
        samples=5,
        shader=faceted interp,
        domain=-3:3,] 
        {exp(-x^2-y^2)};
    \end{axis}
\end{tikzpicture}
\end{document}

This samples a non-normalized normal distribution on a matrix of 5 x 5 bicubic patches. Each bicubic patch has 4x4 points (!).

Note that it also supports to sample non-bezier surfaces in which case you would typically use a surf plot, either from matrix input or using a math expression like

\documentclass[a4paper]{article}

\usepackage{pgfplots}

\begin{document}
\thispagestyle{empty}

\begin{tikzpicture}
    \begin{axis}
    \addplot3[
        surf,
        samples=15,
        shader=faceted interp,
        domain=-3:3,] 
        {exp(-x^2-y^2)};
    \end{axis}
\end{tikzpicture}
\end{document}

I chose a relatively small sampling density which could be increased by adopting the samples key (as you guessed).

DISCLAIMER NOTE: I am author of pgfplots.

You may want to inspect the manual at http://pgfplots.sourceforge.net/pgfplots.pdf to see if it suits your needs; its many example should make it simple to decide if this is what you had in mind.