How to compute normal in quartic Walton-Meek's Gregory patch in tessellation shader?

Question

I am programing a tessellation shader in OpenGL which computes the quartic Walton-Meek's Gregory patch. I am searching for a local G1 method with good shading/visual results. So I am trying this patch. I didn't get good (visual/shading) results with PN-Triangles.

I've just had my first result with this WM patch, which isn't good as well. Perhaps the normal I am calculating is wrong. I am following these equations to compute it:

Then, I actually compute it through normalizing the cross product between the derivatives (last image).

If you want to take a look at the shader code, here it is (the first one is the tessellation control shader and the second one is the tessellation evaluation shader):

#version 430 core
layout (vertices = 3) out;

in VertOut
{
    vec3 normal;
} vertOut[];

out TescOut
{
    vec3 p0;
    vec3 p1;
    vec3 p2;
    vec3 p3;
    vec3 g0;
    vec3 g1;
    vec3 n;
} tescOut[];

void main()
{
    const float kTessLevel = 12.0;
    gl_TessLevelOuter[gl_InvocationID] = kTessLevel;
    gl_TessLevelInner[0] = kTessLevel;

    vec3 p0 = tescOut[gl_InvocationID].p0 = gl_in[gl_InvocationID].gl_Position.xyz;
    vec3 n = tescOut[gl_InvocationID].n = vertOut[gl_InvocationID].normal;

    const int nextInvID = gl_InvocationID < 2 ? gl_InvocationID + 1 : 0;
    vec3 edge = gl_in[nextInvID].gl_Position.xyz - p0;
    vec3 nNext = vertOut[nextInvID].normal;
    float d = length(edge), a = dot(n, nNext);
    vec3 gama = edge / d;
    float a0 = dot(n, gama), a1 = dot(nNext, gama);
    float ro = 6.0 * (2.0 * a0 + a * a1)/(4.0 - a * a);
    float sigma = 6.0 * (2.0 * a1 + a * a0)/(4.0 - a * a);

    vec3 v[4] = vec3[4]
    (
        p0,
        p0 + (d / 18.0) * (6.0 * gama - 2.0 * ro * n + sigma * nNext),
        gl_in[nextInvID].gl_Position.xyz - (d / 18.0) * (6.0 * gama + ro * n - 2.0 * sigma * nNext),
        edge = gl_in[nextInvID].gl_Position.xyz
    );

    vec3 w[3] = vec3[3]
    (
        v[1] - v[0],
        v[2] - v[1],
        v[3] - v[2]
    );

    vec3 A[3] = vec3[3]
    (
        cross(n, normalize(w[0])),
        vec3(0.0),
        cross(nNext, normalize(w[2]))
    );

    A[1] = normalize(A[0] + A[2]);

    vec3 l[5] = vec3[5]
    (
        v[0],
        0.25 * (v[0] + 3.0 * v[1]),
        0.25 * (2.0 * v[1] + 2.0 * v[2]),
        0.25 * (3.0 * v[2] + v[3]),
        v[3]
    );

    vec3 p1 = tescOut[gl_InvocationID].p1 = l[1];
    vec3 p2 = tescOut[gl_InvocationID].p2 = l[2];
    vec3 p3 = tescOut[gl_InvocationID].p3 = l[3];

    barrier();

    const int previousInvID = gl_InvocationID > 0 ? gl_InvocationID - 1 : 2;
    vec3 D[4] = vec3[4]
    (
        tescOut[previousInvID].p3 - 0.5 * (p0 + p1),
        vec3(0.0),
        vec3(0.0),
        tescOut[nextInvID].p1 - 0.5 * (p3 + tescOut[nextInvID].p0)
    );

    float mi[2] = float[2](dot(D[0], A[0]), dot(D[3], A[2]));
    float lambda[2] = float[2](dot(D[0], w[0])/dot(w[0], w[0]), dot(D[3], w[2])/dot(w[2], w[2]));

    tescOut[gl_InvocationID].g0 = 0.5 * (l[1] + l[2]) + (2.0/3.0) * (lambda[0] * w[1] + mi[0] * A[1]) + (1.0/3.0) * (lambda[1] * w[0] + mi[1] * A[0]);
    tescOut[gl_InvocationID].g1 = 0.5 * (l[2] + l[3]) + (1.0/3.0) * (lambda[0] * w[2] + mi[1] * A[2]) + (2.0/3.0) * (lambda[1] * w[1] + mi[1] * A[1]);
}

#version 430 core
layout(triangles, equal_spacing, ccw) in;

in TescOut
{
    vec3 p0;
    vec3 p1;
    vec3 p2;
    vec3 p3;
    vec3 g0;
    vec3 g1;
    vec3 n;
} tescOut[];

out TeseOut
{
    vec3 normal;
    vec3 viewPosition;
} teseOut;

uniform mat4 mvp;
uniform mat4 modelView;
uniform mat4 normalMatrix;
uniform bool isNormalLinearlyInterpolated;

#define uvw     gl_TessCoord

const float u = uvw.x, u2 = u * u, u3 = u2 * u, u4 = u2 * u2;
const float v = uvw.y, v2 = v * v, v3 = v2 * v, v4 = v2 * v2;
const float w = uvw.z, w2 = w * w, w3 = w2 * w, w4 = w2 * w2;

#define p400    tescOut[0].p0
#define p310    tescOut[0].p1
#define p220    tescOut[0].p2
#define p130    tescOut[0].p3
#define G01     tescOut[0].g0
#define G02     tescOut[0].g1
#define p040    tescOut[1].p0
#define p031    tescOut[1].p1
#define p022    tescOut[1].p2
#define p013    tescOut[1].p3
#define G11     tescOut[1].g0
#define G12     tescOut[1].g1
#define p004    tescOut[2].p0
#define p103    tescOut[2].p1
#define p202    tescOut[2].p2
#define p301    tescOut[2].p3
#define G21     tescOut[2].g0
#define G22     tescOut[2].g1

#define B400    u4
#define B040    v4
#define B004    w4
#define B310    (4.0 * u3 * v)
#define B031    (4.0 * v3 * w)
#define B103    (4.0 * u * w3)
#define B220    (6.0 * u2 * v2)
#define B022    (6.0 * v2 * w2)
#define B202    (6.0 * u2 * w2)
#define B130    (4.0 * u * v3)
#define B013    (4.0 * v * w3)
#define B301    (4.0 * u3 * w)
#define B211    (12.0 * u2 * v * w)
#define B121    (12.0 * u * v2 * w)
#define B112    (12.0 * u * v * w2)

#define B300    u3
#define B030    v3
#define B003    w3
#define B210    (3.0 * u2 * v)
#define B021    (3.0 * v2 * w)
#define B102    (3.0 * w2 * u)
#define B120    (3.0 * u * v2)
#define B012    (3.0 * v * w2)
#define B201    (3.0 * w * u2)
#define B111    (6.0 * u * v * w)

vec3 interpolate3D(vec3 p0, vec3 p1, vec3 p2)
{
    return gl_TessCoord.x * p0 + gl_TessCoord.y * p1 + gl_TessCoord.z * p2;
}

void main()
{
    vec4 pos = vec4(interpolate3D(tescOut[0].p0, tescOut[1].p0, tescOut[2].p0), 1.0);
    vec3 normal = normalize(interpolate3D(tescOut[0].n, tescOut[1].n, tescOut[2].n));

    if (u != 1.0 && v != 1.0 && w != 1.0)
    {
        vec3 p211 = (v * G12 + w * G21)/(v + w);
        vec3 p121 = (w * G02 + u * G11)/(w + u);
        vec3 p112 = (u * G22 + v * G01)/(u + v);

        vec3 barPos =   p400 * B400 +
                        p040 * B040 +
                        p004 * B004 +
                        p310 * B310 +
                        p031 * B031 +
                        p103 * B103 +
                        p220 * B220 +
                        p022 * B022 +
                        p202 * B202 +
                        p130 * B130 +
                        p013 * B013 +
                        p301 * B301 +
                        p211 * B211 +
                        p121 * B121 +
                        p112 * B112;

        pos = vec4(barPos, 1.0);

        vec3 dpdu = p400 * B300 +
                    p130 * B030 +
                    p103 * B003 +
                    p310 * B210 +
                    p121 * B021 +
                    p202 * B102 +
                    p220 * B120 +
                    p112 * B012 +
                    p301 * B201 +
                    p211 * B111 ;

        vec3 dpdv = p310 * B300 +
                    p040 * B030 +
                    p013 * B003 +
                    p220 * B210 +
                    p031 * B021 +
                    p112 * B102 +
                    p130 * B120 +
                    p022 * B012 +
                    p211 * B201 +
                    p121 * B111 ;

        normal = normalize(cross(dpdu, dpdv));
    }

    gl_Position = mvp * pos;
    pos = modelView * pos;
    teseOut.viewPosition = pos.xyz / pos.w;

    teseOut.normal = (normalMatrix * vec4(normal, 0.0)).xyz;
}

Some screenshots from my current results:

In the "good ones" images (first and third), the normals are computed by linear interpolation, while in the bad ones the normals are computed through the equations I said previously and are shown in the code.

So, how can I correctly compute the normals? Thanks in advance!

@edit

So, I tried @Reynolds (see comment below) derivatives, but still it didn't work properly, unfortunately. These are the results:

Results at the left (the dark ones) are using the code below and the ones at the right are flipping the normal (i.e., -normalize(cross(dpdu,dpdv)). The code:

vec3 dpdu = 4.0 * (B300 * (p301 - p400) + B030 * (p031 - p130) + B003 * (p004 - p103) + B210 * (p211 - p310) +
                    B120 * (p121 - p220) + B201 * (p202 - p301) + B102 * (p103 - p202) + B021 * (p022 - p121) +
                    B012 * (p013 - p112) + B111 * (p112 - p211));

vec3 dpdv = 4.0 * (B300 * (p310 - p400) + B030 * (p040 - p130) + B003 * (p013 - p103) + B210 * (p220 - p310) +
                    B120 * (p130 - p220) + B201 * (p211 - p301) + B102 * (p112 - p202) + B012 * (p022 - p112) +
                    B021 * (p031 - p121) + B111 * (p121 - p211));
normal = normalize(cross(dpdu, dpdv)); // -normalize(cross(dpdu, dpdv));

@edit2

Some more results:

See as the bottom part of the mesh disappeared. It seems to me a problem with the zero singularity. Probably I am not treating this correctly.

Here, the normals of the bottom part of the mesh have a negative component in Y (vertical axis). In the previous picture, Y = 0 (for normals).

Top-left: raw mesh with normals in white Top-right: raw mesh shading without tessellation (hence, without WM patch) Bottom-left: wireframe mode of mesh with tessellation (hence with WM patch) Bottom-right: mesh with tessellation (hence with WM patch)

Any ideas?? Thanks!

Answer 1

Computing exact derivatives of Gregory patches is hard due to the rational blending that occurs for the inner control points. Many people thus opt for an easier solution where the rational blending functions are first evaluated and then the resulting control points are interpreted as the control points of a standard polynomial Bézier triangle. This technique is used for instance in ACC-2 to approximate the normal. This will introduce singularities at the vertices due a zero denominator in the blending functions of the control points but these can easily be fixed as the normal is simply just the vertex normal at these instances.

Remember that the partial derivatives of a Bézier triangle of degree $p$ is a Bézier triangle of degree $p-1$ where the control points in that case are vectors between the control points. As seen from the figure below you can evaluate each partial derivative as a Bézier triangle. In your case you are using just control points, but you have to use differences between control points.

$$\frac{\partial P}{\partial u} = 4 \left(B_{300} (P_{400} - P_{301}) + B_{030} (P_{130} - P_{031}) + B_{003} (P_{103} - P_{004} ) + B_{210} ( P_{310} - P_{211} ) + B_{120} (P_{220} - P_{121}) + B_{201} (P_{301} - P_{202}) + B_{102} (P_{202} - P_{103}) + B_{021} (P_{121} - P_{022}) + B_{012} (P_{112} - P_{013}) + B_{111} (P_{211} - P_{112}) \right) $$

and

$$\frac{\partial P}{\partial v} = 4 \left(B_{300} (P_{310} - P_{301}) + B_{030} (P_{040} - P_{031}) + B_{003} (P_{013} - P_{004}) + B_{210} (P_{220} - P_{211}) + B_{120} (P_{130} - P_{121}) + B_{201} (P_{211} - P_{202}) + B_{102} (P_{112} - P_{103}) + B_{012} (P_{022} - P_{013}) + B_{021} (P_{031} - P_{022}) + B_{111} (P_{121} - P_{112}) \right). $$

Then proceed the same way in producing the normal by the cross product of the two partial derivatives.

To get some intuition about why this is this case we can look at the simpler case of a cubic Bézier curve with control points $[P_{30}, P_{21}, P_{12}, P_{03}]$

$$C(u) = (1-u)^3 P_{30} + 3 (1-u)^2 u P_{21} + 3 (1-u) u^2 P_{12} + u^3 P_{03},$$

if we take the derivative with respect to $u$ we obtain

$$\frac{\partial C(u)}{\partial u} = -3 (1-u)^2 P_{30} - 6 (1-u) u P_{21} + 3 (1-u)^2 P_{21} - 3 u^2 P_{12} + 6 (1-u) u P_{12} + 3 u^2 P_{03}.$$

If we then rearrange the terms we get

$$\frac{\partial C(u)}{\partial u} = 3 \left((1-u)^2 (P_{21} - P_{30}) + 2 (1-u) u (P_{12} - P_{21}) + u^2 ( P_{03} - P_{12})\right) .$$

As you can see the derivative is a quadratic Bézier where the control points are the vectors between the control points of the initial cubic Bézier curve. This same principle holds for Bézier triangles.