<p>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 <a href="https://dl.acm.org/citation.cfm?id=1618497" rel="nofollow noreferrer">ACC-2</a> 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.</p>

<p>Remember that the partial derivatives of a Bézier triangle  of degree <span class="math-container">$p$</span> is a Bézier triangle of degree <span class="math-container">$p-1$</span> 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. </p>

<p><a href="https://i.sstatic.net/4dxSV.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/4dxSV.png" alt="cubic derivatives"></a></p>

<p><span class="math-container">$$\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) $$</span></p>

<p>and</p>

<p><span class="math-container">$$\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). $$</span></p>

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

<p>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 <span class="math-container">$[P_{30}, P_{21}, P_{12}, P_{03}]$</span></p>

<p><span class="math-container">$$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},$$</span></p>

<p>if we take the derivative with respect to <span class="math-container">$u$</span> we obtain</p>

<p><span class="math-container">$$\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}.$$</span></p>

<p>If we then rearrange the terms we get</p>

<p><span class="math-container">$$\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) .$$</span></p>

<p>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.</p>

CC BY-SA 4.0
2019-11-08T16:35:51.420Z
First of all, thank you very much for your answer, Reynolds!

About what you said that many people opt to evaluate and then interpret as standar polynomial Bézier triangle, isn't that exactly what I am doing? The equations I put as images are from Bézier-Bernstein triangle according to this [article](https://www.gamasutra.com/view/feature/131389/b%C3%A9zier_triangles_and_npatches.php?print=1)

 Which singularities do you mean? I only know the zero denominator singularity. By the way, how exactly can they be removed?
YardenJ2R
CC BY-SA 4.0
2019-11-08T16:42:38.323Z
Are your images and equations correct? Why isn't there P040 and why P310 is repeated? Why should I use the differences instead of only "downgrading" it to the Bézier triangles and finding control points between control points? Can you please explain to me in more detail?

I used the two partial derivatives equations you wrote. Unfortunately, the results aren't good yet.

Again, thank you very much for your help.
YardenJ2R
CC BY-SA 4.0
2019-11-08T17:00:18.123Z
Sure, no problem! I really appreciate your help.
YardenJ2R
CC BY-SA 4.0
2019-11-12T15:02:47.670Z
Reynolds, unfortunately, it didn't work yet. I've updated the question with the new results. Any ideas? Thanks again for your help!
YardenJ2R
CC BY-SA 4.0
2019-11-12T16:09:27.470Z
Unfortunately I got the direction of the derivatives wrong, again. These should be the right derivatives.
Reynolds
CC BY-SA 4.0
2019-11-12T16:15:29.100Z
At the second derivative, should it be $B_{300}(P_{310} - P_{301})$? You repeated the control point.
YardenJ2R
CC BY-SA 4.0
2019-11-12T16:17:48.987Z
Yes you are right. Please also consult the figures if you encounter more errors.
Reynolds
CC BY-SA 4.0
2019-11-12T16:24:30.020Z
Ok! I tried it out. Unfortunately, the results are the same as the last ones with normal flipping (the last pictures on the right). Even so, thanks for your help!
YardenJ2R
CC BY-SA 4.0
2019-11-19T15:12:17.160Z
Hey there, Reynolds! Would you mind elaborating on how should I treat the vertices and its normals considering the zero singularity, please? Perhaps this is where I am doing wrong.
YardenJ2R
CC BY-SA 4.0
2019-11-08T12:52:42.527Z
2019-11-12T16:16:38.137Z
2019-11-12T16:16:38.137Z
Reynolds
<p>I am programing a tessellation shader in OpenGL which computes the <a href="https://reader.elsevier.com/reader/sd/pii/0010448595000461?token=4190E80EF50B0C8E2FB4EADFFEF26DF93ACEDF3AED59AA6B97D59B2C1FC52B0B965B255B27749ECC11F49B44C1610688" rel="nofollow noreferrer">quartic Walton-Meek's Gregory patch</a>. 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.</p>

<p>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:</p>

<p><a href="https://i.sstatic.net/AQq0v.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/AQq0v.png" alt="Bézier-Bernstein Triangle and Derivatives Equations"></a></p>

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

<p>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):</p>

<pre><code>#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 &lt; 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 &gt; 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]);
}
</code></pre>

<pre><code>#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 &amp;&amp; v != 1.0 &amp;&amp; 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;
}
</code></pre>

<p>Some screenshots from my current results:</p>

<p><a href="https://i.sstatic.net/YIW42.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/YIW42.png" alt="enter image description here"></a>
<a href="https://i.sstatic.net/7bfza.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/7bfza.png" alt="enter image description here"></a>
<a href="https://i.sstatic.net/5rw2D.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/5rw2D.png" alt="enter image description here"></a>
<a href="https://i.sstatic.net/a1MBK.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/a1MBK.png" alt="enter image description here"></a></p>

<p>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.</p>

<p>So, how can I correctly compute the normals? Thanks in advance!</p>

<p>@edit</p>

<p>So, I tried @Reynolds (see comment below) derivatives, but still it didn't work properly, unfortunately. These are the results: <a href="https://i.sstatic.net/d26TV.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/d26TV.png" alt="New results"></a></p>

<p>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:</p>

<pre class="lang-cpp prettyprint-override"><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));
</code></pre>

<p>@edit2</p>

<p>Some more results:</p>

<p><a href="https://i.sstatic.net/AVd7g.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/AVd7g.png" alt="Example of some errors"></a>
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.</p>

<p><a href="https://i.sstatic.net/NaNdJ.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/NaNdJ.png" alt="Another example"></a>
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).</p>

<p>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)</p>

<p>Any ideas?? Thanks!</p>

<opengl><glsl><computational-geometry><tesselation>
How to compute normal in quartic Walton-Meek's Gregory patch in tessellation shader?
YardenJ2R
# How to compute normal in quartic Walton-Meek's Gregory patch in tessellation shader?

(Asked by: username_0 on 2019-11-07)

I am programing a tessellation shader in OpenGL which computes the [quartic Walton-Meek's Gregory patch](https://reader.elsevier.com/reader/sd/pii/0010448595000461?token=4190E80EF50B0C8E2FB4EADFFEF26DF93ACEDF3AED59AA6B97D59B2C1FC52B0B965B255B27749ECC11F49B44C1610688). 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:

[![Bézier-Bernstein Triangle and Derivatives Equations](https://i.sstatic.net/AQq0v.png)](https://i.sstatic.net/AQq0v.png)

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:

[![enter image description here](https://i.sstatic.net/YIW42.png)](https://i.sstatic.net/YIW42.png)
[![enter image description here](https://i.sstatic.net/7bfza.png)](https://i.sstatic.net/7bfza.png)
[![enter image description here](https://i.sstatic.net/5rw2D.png)](https://i.sstatic.net/5rw2D.png)
[![enter image description here](https://i.sstatic.net/a1MBK.png)](https://i.sstatic.net/a1MBK.png)

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 @username_1 (see comment below) derivatives, but still it didn't work properly, unfortunately. These are the results: [![New results](https://i.sstatic.net/d26TV.png)](https://i.sstatic.net/d26TV.png)

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:

```cpp
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:

[![Example of some errors](https://i.sstatic.net/AVd7g.png)](https://i.sstatic.net/AVd7g.png)
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.

[![Another example](https://i.sstatic.net/NaNdJ.png)](https://i.sstatic.net/NaNdJ.png)
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!

**Tags:** opengl, glsl, computational-geometry, tesselation


---

## 1 answers

### Answer by username_1 on 2019-11-08. Score: 4

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](https://dl.acm.org/citation.cfm?id=1618497) 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.

[![cubic derivatives](https://i.sstatic.net/4dxSV.png)](https://i.sstatic.net/4dxSV.png)

$$\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.

> **Comments:**
> **username_0**: First of all, thank you very much for your answer, username_1!  About what you said that many people opt to evaluate and then interpret as standar polynomial Bézier triangle, isn't that exactly what I am doing? The equations I put as images are from Bézier-Bernstein triangle according to this [article](https://www.gamasutra.com/view/feature/131389/b%C3%A9zier_triangles_and_npatches.php?print=1)   Which singularities do you mean? I only know the zero denominator singularity. By the way, how exactly can they be removed?  
> **username_0**: Are your images and equations correct? Why isn't there P040 and why P310 is repeated? Why should I use the differences instead of only "downgrading" it to the Bézier triangles and finding control points between control points? Can you please explain to me in more detail?  I used the two partial derivatives equations you wrote. Unfortunately, the results aren't good yet.  Again, thank you very much for your help.  
> **username_0**: Sure, no problem! I really appreciate your help.  
> **username_0**: username_1, unfortunately, it didn't work yet. I've updated the question with the new results. Any ideas? Thanks again for your help!  
> **username_1**: Unfortunately I got the direction of the derivatives wrong, again. These should be the right derivatives.  
> **username_0**: At the second derivative, should it be $B_{300}(P_{310} - P_{301})$? You repeated the control point.  
> **username_1**: Yes you are right. Please also consult the figures if you encounter more errors.  
> **username_0**: Ok! I tried it out. Unfortunately, the results are the same as the last ones with normal flipping (the last pictures on the right). Even so, thanks for your help!  
> **username_0**: Hey there, username_1! Would you mind elaborating on how should I treat the vertices and its normals considering the zero singularity, please? Perhaps this is where I am doing wrong.  

---

