Automatic interior edges in bezier surfaces

Question 1

Looks like 3ds Max constructs a parallelogram to get interior points.

patch[5] = patch[0] + (patch[1] - patch[0]) + (patch[4] - patch[0]);
patch[5] = patch[1] + patch[4] - patch[0];

And the whole function:

void CreatePatchVertices(Vector3[] s, int tessellation, bool isMirrored)
{
    Debug.Assert(s.Length == 16);

    for (int i = 0; i <= tessellation; i++)
    {
        float ti = (float)i / tessellation;

        for (int j = 0; j <= tessellation; j++)
        {
            float tj = (float)j / tessellation;

            // Compute automatic interior edges.
            s[5] = s[1] + s[4] - s[0];
            s[6] = s[2] + s[7] - s[3];
            s[9] = s[8] + s[13] - s[12];
            s[10] = s[11] + s[14] - s[15];

            // Perform four horizontal bezier interpolations
            // between the control points of this patch.
            Vector3 p1 = Bezier(s[0], s[1], s[2], s[3], ti);
            Vector3 p2 = Bezier(s[4], s[5], s[6], s[7], ti);
            Vector3 p3 = Bezier(s[8], s[9], s[10], s[11], ti);
            Vector3 p4 = Bezier(s[12], s[13], s[14], s[15], ti);

            // Perform a vertical interpolation between the results of the
            // previous horizontal interpolations, to compute the position.
            Vector3 position = Bezier(p1, p2, p3, p4, tj);

            // Perform another four bezier interpolations between the control
            // points, but this time vertically rather than horizontally.
            Vector3 q1 = Bezier(s[0], s[4], s[8], s[12], tj);
            Vector3 q2 = Bezier(s[1], s[5], s[9], s[13], tj);
            Vector3 q3 = Bezier(s[2], s[6], s[10], s[14], tj);
            Vector3 q4 = Bezier(s[3], s[7], s[11], s[15], tj);

            // Compute vertical and horizontal tangent vectors.
            Vector3 tangentA = BezierTangent(p1, p2, p3, p4, tj);
            Vector3 tangentB = BezierTangent(q1, q2, q3, q4, ti);

            // Cross the two tangent vectors to compute the normal.
            Vector3 normal = Vector3.Cross(tangentA, tangentB);
            normal.Normalize();

            // Create the vertex.
            AddVertex(position, normal);
        }
    }
}

This is what a teapot primitive with automatic interior points looks like both in 3ds Max and Primitives3D sample:

Question 2

This answer must be redundant but I would like to point your attention to another way of constructing bezier surfaces called Coons patches. Coons patches explicitly use only boundary curves, the conversion to bicubic patches is quite straight forward.

The reason I would recommend this, is that, it was developed for extrapolating boundary feature lines typical of automotive designs rather precisely and to calculate the interior points (for bicubic bezier patches) or tessellation automatically.

Given four boundary curves bu0, bu1, bv0, bv1;

bu0, bu1  // boundaries in u direction u <- [0, 1]
bv0, bv1  // in v direction  v <- [0, 1]
// bu0, bu1, bv0, bv1 are arrays of bezier coefficients

Interior points and points on surface of a coons patch can be evaluated by first evaluating ruled surfaces in u and v directions independently, and a bilinear interpolation in u and v both.

// In matrix form, so that you can reduce it all together

                      ⎡        ⎤
ruled_u = [ (1-u) u ] ⎜ bv0(v) ⎟
                      ⎜ bv1(v) ⎟
                      ⎣        ⎦ 

// Similarly
                      ⎡        ⎤
ruled_v = [ (1-v) v ] ⎜ bu0(u) ⎟
                      ⎜ bu1(u) ⎟
                      ⎣        ⎦ 

// Assuming cubic bezier curves as boundary, if you select appropriate control point 
// instead of bv0(u) etc, i.e. bv0[1] for u = 1/3 and bv0[2] for u = 2/3, we can convert
// a coons patch into a bicubic patch, instead of evaluating the point on surface directly

                          ⎡               ⎤ ⎡         ⎤
bilinear_uv = [ (1-u) u ] ⎜ bu0(0) bu0(1) ⎟ ⎜ (1 - v) ⎟
                          ⎜ bu1(0) bu1(1) ⎟ ⎜    v    ⎟
                          ⎣               ⎦ ⎣         ⎦

// Then the 'interior' point at u and v is
Coons (u, v) = ruled_u + ruled_v - bilinear_uv

On a coons patch the interior control points can be calculated as

Coons(1/3, 1/3)
Coons(1/3, 2/3)
Coons(2/3, 1/3)
Coons(2/3, 2/3)

while just selecting bv0[1] for 1/3 instead of bv0(1/3) etc. in above equations.

In Coons patches, the interpolation feels more natural than just linearly interpolating boundary control points. Also is more 'correct' in my experience. Not just aesthetically but parametrically, after tessellating as well.

Recently, I had to generate bezier surfaces approximating sRGB gamut in LCH colour space. It is just a transformed and distorted RGB cube. All I had was boundaries of the patches. The interior values were wrong while using just linear interpolation. Finally coons patches gave the least error by far.

Here is a sample of the generated bicubic surface from a set coons patches.

enter image description here