The Unchained Blog

OpenGL Quad Interpolation

problem & solution



   
  
  by C.G. Sleuth
December, 2020

Quadrilaterals - four-sided polygons - might be more common in design than we realize. For some applications, like display of font characters or representation of ruled surfaces, surfaces of revolution, or patches, they are a more natural design element than triangles. You may know that OpenGL divides all quadrilaterals into triangles before rendering, but did you know there can be a problem inherent with this division?

The Problem

When OpenGL divides a quadrilateral into separately rasterized triangles, the results are often artifact-free, such as the following square with a grid texture.

 

Figure 1: Texture-mapped square rendered as two triangles
left: untransformed, right: oblique, perspective view

The geometry (the four points that define the square) and the texture domain (defined by the four corresponding points on the texture map) have the same, square shape. In perspective, the geometry appears as a trapezoid, but the texture-mapped image remains artifact-free. If, rather than a square, we start with a trapezoid as our geometry and divide it into two triangles that are rendered separately, the result contains a C1 discontinuity along the diagonal. Why does the trapezoid below have an artifact but the trapezoid above does not?

 

Figure 2: Texture-mapped trapezoid rendered as two triangles

To answer this question, consider a single, non-degenerate triangle defined by three screen-space points and an associated triangle of three texture-space points (i.e., uv-coordinates). A unique, affine transformation exists that maps the screen triangle to the texture triangle. This guarantees that the parallel rows of pixels produced when rasterizing the screen-space triangle correspond with parallel rows of samples along a texture map.

Figure 3: Rasterization
left: pixels in screen-space, right: samples in texture-space

Let us now add a fourth screen-space point and a fourth texture point, creating a second screen triangle and a second texture triangle. The rows of pixels in the two screen triangles obviously align with each other, but for the rows of samples to align between the two texture triangles, the affine transformation described above must be applied to the new screen location to determine the new texture location. If the new texture-space point is placed elsewhere, a discontinuity in the texture gradient occurs at the shared boundary of the triangles.


Figure 4: Effect on gradient of placement of fourth texture point
ideal location of new texture point shown as circle

The appendix provides code to test whether a quad's texture coordinates are affinely related to its corner locations.

Sometimes Unavoidable

It isn't always possible that the texture shape is an affine image of the geometric shape. For example, the longitude/latitude parameterization typical of a sphere implies that the quad geometry becomes increasingly trapezoidal the closer the quad is to a pole, whereas its texture domain remains rectangular. As the foreshortening increases, the texture discontinuity increases.


Figure 5: Quads on sphere
left: geometric space, middle: texture space, right: rendered as triangles

In the extreme, at the pole, Q2 is a quadrilateral whose degenerate upper edge is a point. When divided into triangles, geometrically, one triangle covers no pixels, so its corresponding texture is ignored.

 

Figure 6: Triangulated, degenerate quad
left: geometric space, right: texture space

Solution: Quad Bilinear Interpolation

The C1 texture discontinuity has been of concern since the late 1970s. In 2004, Hormann and Tarini proposed a quadrilateral rendering primitive that uses generalized barycentric coordinates to produce a bilinear interpolation of vertex attributes.

Their method can be implemented with the following OpenGL geometry shader. With GL_LINES_ADJACENCY as the first argument to glDrawArrays or glDrawElements, four vertices from the vertex shader are sent to the geometry shader as an array, for each quad (the draw subroutine does not accept GL_QUADS as an input primitive, but GL_LINES_ADJACENCY achieves the same result).

The perspective space coordinates of the vertices (gpos0-3) and the texture coordinates (guv0-3) are then sent to the rasterizer (as separate variables, not arrays) and then to the pixel shader for use in computing the proper uv-coordinates.

const char *geometryShader = R"(
    #version 330
    layout (lines_adjacency) in;                                       // 1 quad in
    layout (triangle_strip, max_vertices = 6) out;                     // 2 tris out
    in vec3 vpoint[4], vnormal[4];
    in vec2 vuv[4];
    out vec3 gpoint, gnormal;
    flat out vec4 gpos0, gpos1, gpos2, gpos3;                          // same each vrt
    flat out vec2 guv0, guv1, guv2, guv3;                              // same each vrt
    void EmitPoint(int i, bool provokeFlatOutput) {
        gpoint = vpoint[i];
        gnormal = vnormal[i];
        gl_Position = gl_in[i].gl_Position;
        if (provokeFlatOutput) {                                       // once per tri
            gpos0 = gl_in[0].gl_Position; gpos1 = gl_in[1].gl_Position;
            gpos2 = gl_in[2].gl_Position; gpos3 = gl_in[3].gl_Position;
            guv0 = vuv[0]; guv1 = vuv[1]; guv2 = vuv[2]; guv3 = vuv[3];
        }
        EmitVertex();                                                  // send vertex
    }
    void main() {
        EmitPoint(0, false); EmitPoint(1, false); EmitPoint(3, true);
        EndPrimitive();                                                // triangle 1
        EmitPoint(1, false); EmitPoint(2, false); EmitPoint(3, true);
        EndPrimitive();                                                // triangle 2
    }
)";

In UV(), below, the built-in gl_FragCoord input to the pixel shader is converted to perspective-space coordinates and used to compute the barycentric weights with respect to the quad corners, gpos0-3. The perspective-corrected weights are applied to guv0-3 to determine the correct, bilinearly interpolated uv-coordinates for the pixel. Double precision is used in BarycentricWeights() and UV().

    const char *pixelShader = R"(
    #version 430
    in vec3 gpoint, gnormal;
    flat in vec4 gpos0, gpos1, gpos2, gpos3;                           // quad points
    flat in vec2 guv0, guv1, guv2, guv3;                               // quad uvs
    uniform sampler2D image;
    uniform vec4 vp;                                                   // viewport
    uniform vec3 light;
    double Cross2d(dvec2 a, dvec2 b) { return(a.x*b.y)-(a.y*b.x); }
    dvec4 BarycentricWeights(dvec2 p) {
        double w[4] = { gpos0.w,gpos1.w,gpos2.w,gpos3.w }, r[4], t[4], u[4];
        dvec2 v[4] = { dvec2(gpos0)/w[0], dvec2(gpos1)/w[1], dvec2(gpos2)/w[2], dvec2(gpos3)/w[3] }, s[4];
        for (int i = 0; i < 4; i++) {
            s[i] = v[i]-p;
            r[i] = length(s[i])*sign(w[i]);
        }
        for (int i = 0; i < 4; i++) {
            double A = Cross2d(s[i], s[(i+1)%4]);
            double D = dot(s[i], s[(i+1)%4]);
            t[i] = (r[i]*r[(i+1)%4]-D)/A;
        }
        for (int i = 0; i < 4; i++)
            u[i] = (t[(i+3)%4]+t[i])/r[i];
        return dvec4(u[0],u[1],u[2],u[3])/(u[0]+u[1]+u[2]+u[3]);
    }
    vec2 UV() {
        // from A Quadrilateral Rendering Primitive by Hormann & Tarini
        // get barycentric weights, do persp div, set uv weighted-sum
        dvec2 p = dvec2(2*gl_FragCoord.x/vp[2]-1, 2*gl_FragCoord.y/vp[3]-1);
        vec4 wt = BarycentricWeights(p);
        float f[4] = { wt[0]/gpos0.w, wt[1]/gpos1.w, wt[2]/gpos2.w, wt[3]/gpos3.w };
        return (f[0]*guv0+f[1]*guv1+f[2]*guv2+f[3]*guv3)/(f[0]+f[1]+f[2]+f[3]);    
    }
    void main() {
        vec3 N = normalize(gnormal), L = normalize(light-gpoint);
        float d = max(0, dot(N, L)), i = clamp(.05+d, 0, 1);
        gl_FragColor = vec4(i*texture(image, UV()).rgb, 1);
    }
)";

 

Figure 7: Quads on sphere, rendered as triangles
left: with typical shaders, right: with above shaders

Perspective Correction

The diagonals in the texture above appear curved due to the curvature of the sphere. When viewed in perspective, however, lines can appear curved, even on a flat surface, if the barycentric weights are not corrected for perspective. The correction is performed by the division by gpos0-4.w when computing v[ ] in BarycentricWeights() and f[ ] in UV(). As a result, straight lines remain straight, as they should.

 

Figure 8: Textured square
left: quad rendered with bilinear interpolation without perspective correction, right: with correction
A blog post by Reed (Quadrilateral Interpolation) solves the texture discontinuity problem with bilinear interpolation, but does not account for perspective.

Et too, Shading?

C1 discontinuity also occurs when shading a quad without texture. A perspective-correct bilinear interpolation can be made by replacing guv0-3 in the geometry and pixel shaders with gcolor0-3.

 

Figure 9: A trapezoid with three dark corners and one bright corner
left: divided into triangles, right: bilinearly interpolated

The OpenGL SuperBible (6th ed., ch. 8, Rendering Quads using a Geometry Shader) and a blog post by Izdebski (How to Correctly Interpolate Vertex Attributes on a Parallelogram Using Modern GPUs) consider this shading issue, although the proposed solutions do not eliminate texture discontinuity.

Conclusion

An alternative method to reduce the triangulation artifact is to increase resolution. For example, for the sphere, the number of circles of latitude can be increased, especially near the poles. Even at twice the resolution, however, the artifact remains problematic.


Figure 10: Comparison
left: bilinear quads, right: two-fold resolution, divided into triangles

UV() requires 100+ arithmetic operations, so one good strategy is to invoke bilinear interpolation only for quads whose geometric and texture coordinates are not affinely related.

This post extends a discussion in Computer Graphics: Implementation and Explanation.

(Email the author)

Appendix: Testing a Quad's Texture / Geometry Correlation

If a quad's texture coordinates are an affine image of (i.e., affinely related to) its corners, the rate and direction of texture samples will match across the triangles and there will be no texture gradient discontinuity. If they are not an affine image, bilinear quad interpolation can be used to avoid the discontinuity.

Here is a simple test for affinity. For each corner of the quad, compute its barycentric coordinates with respect to the other three corners (this is a different use for barycentric coordinates than the use above for bilinear interpolation). Also compute the barycentric coordinates of the corresponding texture-space location with respect to the other three texture locations. If there is affinity, these two barycentric coordinates will match.


Figure 11: Barycentric computation for p inside or outside triangle
 
      vec3 cross(vec3 a, vec3 b) { return vec3(a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x); }
      float cross(vec2 a, vec2 b) { return vec2(a.x*b.y)-(a.y*b.x); }
      float area(vec2 p1, vec2 p2, vec2 p3) { return cross(p2-p1, p3-p2); }
      float area(vec3 p1, vec3 p2, vec3 p3) { (vec3 v = cross(p2-p1, p3-p2)); return (v.z > 0? 1 : -1)*length(v); }
      template <Type T> bool barycentric(T p, T p1, T p2, T p3, vec3 &bary) {
          float a1 = area(p, p2, p3), a2 = area(p, p3, p1), a3 = area(p, p1, p2), sum = a1+a2+a3;
          if (abs(sum) < FLT_EPSILON)
              return false;
          bary = vec3(a1, a2, a3)/sum;
          return true;
      }
bool affine_mapped(vec3 *pt, vec2 *uv, float tolerance = .01) {
    // is uv[] an affine image of pt[]?
    for (int i = 0; i < 4; i++) {
        int i1 = mod(i+1, 4), i2 = mod(i+2, 4), i3 = mod(i+3, 4);
        vec3 bpt, buv;
        bool ptOk = barycentric(pt[i], pt[i1], pt[i2], pt[i3], bpt);
        bool uvOk = barycentric(uv[i], uv[i1], uv[i2], uv[i3], buv);
        if (!ptOk || !uvOk || length(bpt-buv) > tolerance)
            return false;
    }
    return true;
}

If a fourth point is outside the triangle defined by the first three points, the weights will be of mixed sign. The barycentric coordinates are invariant under affine transformation (scale, shear, rotation, or translation) applied to the geometry or the texture coordinates. They are also invariant under perspective if they are inverse-scaled by the homogeneous coordinate.