Purpose

The purpose of this article is to propose a new method that preserves C^1 continuity over the common edge of two generated triangles from convex quadrilaterals. This new method is based on an algebraic solution for the Bilinear interpolation coefficient expressed in Barycentric coordinates. Bilinear interpolation has the advantage of being the simplest interpolation from an algebraic perspective. Consequently, the computational overhead is negligibly small. Additionally, linear interpolation allows for the easy construction of other types of interpolation, such as polynomial interpolation. The algebraic solution will then be implemented and tested using the three available hardware-accelerated pipelines supported by GPU hardware.

State of the Art

Currently, efforts to address the problem outlined in this article can be categorized into the following approaches:

Proposition of the New Method

The method presented in the article consists of an equation that can express Bilinear interpolation coefficients u, v using Barycentric coordinates \alpha, \beta, \gamma of one of the two triangles into which each convex quadrilateral can be divided.

The Euclidean coordinates of every point P in a given quadrilateral on the plane can be expressed as a Bilinear interpolation of the Euclidean coordinates of its vertices.

The line segment AB is parameterized as E(u)=A + u(B-A), , u \in[0,1],

the line segment DC is parameterized as F(u)=D + u(C-D), , u \in[0,1],

the line segment DA is parameterized as G(v)=A + v(D-A), , v \in[0,1],

the line segment BC is parameterized as H(v)=G + v(H-G), , v \in[0,1].

The point P is the intersection of the two segments E(u)F(u) and G(v)H(v).

Alternatively, the Euclidean coordinates of point P can be expressed as a Barycentric coordinates of

a triangle ABD \ P = \alpha A + \beta B + \gamma D

where \alpha,\beta,\gamma are positive real numbers such that \alpha + \beta + \gamma = 1..

This can be shown geometrically in two-dimensional space.

The Euclidean coordinates of point P are the same whether it is inside the quadrilateral ABCD or the triangle ABC.

Depending on the Euclidean coordinates of point P, it may lie inside triangle ABC or triangle ACD. There is a system of equations that describe the coordinates of point P using the vertices of quadrilateral ABCD and the Bilinear interpolation coefficients u,v, as well as the vertices of triangles ABC or ACD and the Barycentric coordinates for the same point P.

The coordinates of point P in \mathbb{R}^2 can be calculated using Barycentric coordinates \alpha, \beta, \gamma and Bilinear interpolation coefficients u and v, which can be expressed by the following two independent set of equations:

The coordinates of point P in \mathbb{R}^2 can be calculated using Barycentric coordinates \alpha, \beta, \gamma and Bilinear interpolation coefficients u and v, which can be expressed by the following two independent sets of equations:

\begin{cases} A_x + v(D_x - A_x) = G_x \\ A_y + v(D_y - A_y) = G_y \\ B_x + v(C_x - B_x) = H_x \\ B_y + v(C_y - B_y) = H_y \\ G_x + s(H_x - G_x) = P_x \\ G_y + s(H_y - G_y) = P_y \\ \alpha A_x + \beta B_x + \gamma D_x = P_x \\ \alpha A_y + \beta B_y + \gamma D_y = P_y \\ \end{cases} \newline \newline \begin{cases} A_x + u(B_x - A_x) = E_x \\ A_y + u(B_y - A_y) = E_y \\ D_x + u(C_x - D_x) = F_x \\ D_y + u(C_y - D_y) = F_y \\ E_x + t(F_x - E_x) = P_x \\ E_y + t(F_x - E_y) = P_y \\ \alpha A_x + \beta B_x + \gamma D_x = P_x \\ \alpha A_y + \beta B_y + \gamma D_y = P_y \\ \end{cases}

The equations can be reduced to the form:

\begin{cases} A_x + u(B_x - A_x) + s(D_x + u(C_x - D_x) - A_x - u(B_x - A_x)) = \alpha A_x + \beta B_x + \gamma D_x \\ A_y + u(B_y - A_y) + s(D_y + u(C_y - D_y) - A_y - u(B_y - A_y)) = \alpha A_y + \beta B_y + \gamma D_y \\ \end{cases} \newline \newline \begin{cases} A_x + v(D_x - A_x) + t(B_x + v(C_x - B_x) - A_x - v(D_x - A_x)) = \alpha A_x + \beta B_x + \gamma D_x \\ A_y + v(D_y - A_y) + t(B_y + v(C_y - B_y) - A_y - v(D_y - A_y)) = \alpha A_y + \beta B_y + \gamma D_y \\ \end{cases}

We can further simplify using:

\begin{cases} \vec{P} = D-A \\ \vec{R} = B-A \\ \vec{Q} = C-B-D+A \\ \alpha = (1-\beta-\gamma) \\ \end{cases}

to a two-system of equations:

\begin{cases} u\vec{R}_x + s(\vec{P}_x + u\vec{Q}_x) = \beta \vec{R}_x + \gamma \vec{P}_x \\ u\vec{R}_y + s(\vec{P}_y + u\vec{Q}_y) = \beta \vec{R}_y + \gamma \vec{P}_y \\ \end{cases} \newline \newline \begin{cases} v\vec{P}_x + t(\vec{R}_x + v\vec{Q}_x) = \beta \vec{R}_x + \gamma \vec{P}_x \\ v\vec{P}_y + t(\vec{R}_y + v\vec{Q}_y) = \beta \vec{R}_y + \gamma \vec{P}_y \\ \end{cases}

Now, solving the two systems of equations for u and v leads to the problem of solving two quadratic equations:

\begin{cases} u^2(\vec{R}_x\vec{Q}_y - \vec{Q}_x\vec{R}_y) + u(\vec{R}_x\vec{P}_y - \vec{P}_x\vec{R}_y + \beta(\vec{Q}_x\vec{R}_y-\vec{R}_x\vec{Q}_y) + \gamma (\vec{Q}_x\vec{P}_y-\vec{P}_x\vec{Q}_y)) - \\ \beta(\vec{P}_x\vec{R}_y-\vec{R}_x\vec{P}_y) = 0\\ \end{cases} \newline \newline \begin{cases} v^2(\vec{P}_x\vec{Q}_y - \vec{Q}_x\vec{P}_y) + v(\vec{P}_x\vec{R}_y - \vec{R}_x\vec{P}_y + \beta(\vec{Q}_x\vec{R}_y-\vec{R}_x\vec{Q}_y) + \gamma(\vec{Q}_x\vec{P}_y-\vec{P}_x\vec{Q}_y)) - \\ \gamma(\vec{R}_x\vec{P}_y-\vec{P}_x\vec{R}_y) = 0\\ \end{cases}

Embedding this solution in \mathbb{R}^3 at the plane z=0 with the surface normal vector for a triangle \vec{N} = (0,0,1), the equation can be written as:

u^2 \vec{N}\cdot(\vec{R}\times \vec{Q}) + u(\vec{N}\cdot(\vec{R}\times \vec{P}) + \beta \vec{N}\cdot(\vec{Q}\times \vec{R}) + \gamma \vec{N}\cdot(\vec{Q}\times \vec{P})) + \beta \vec{N}\cdot(\vec{P}\times \vec{R}) = 0 \newline \newline v^2 \vec{N}\cdot(\vec{P}\times \vec{Q}) + v(\vec{N}\cdot(\vec{P}\times \vec{R}) + \beta \vec{N}\cdot(\vec{Q}\times \vec{R}) + \gamma \vec{N}\cdot(\vec{Q}\times \vec{P})) - \gamma \vec{N}\cdot(\vec{R}\times \vec{P}) = 0

where \times is the cross product in \mathbb{R}^3 and \cdot is the dot product in \mathbb{R}^3.

Substituting the triple dot product into the equation:

\begin{cases} A_u = \vec{N}\cdot(\vec{R}\times \vec{Q}) \\ A_v = \vec{N}\cdot(\vec{P}\times \vec{Q}) \\ B = \vec{N}\cdot(\vec{P}\times \vec{R}) \\ \end{cases}

the final two equations are as follows:

u^2 A_u - u(\beta A_u + \gamma A_v + B) + \beta B = 0 \\ v^2 A_v - v(\beta A_u + \gamma A_v - B) - \gamma B = 0 \\

The last step is to solve these two quadratic equations for the parameters u and v.

The equation can be used to calculate bilinear interpolation coefficients u and v for any quadrilateral in three-dimensional space. This is because if a quadrilateral is embedded in any plane in \mathbb{R}^3, it can be transformed isometrically to the plane z=0, and the premise still holds after isometric mapping.

It is worth noting that the above equations show that the bilinear interpolation coefficients u and v depend on all four Euclidean coordinates of the vertices in the quadrilateral. If we are dealing with static meshes, these coefficients can be calculated in advance (this is not necessary at runtime). However, if we are working with dynamic (animated) meshes, it is necessary to calculate these coefficients at runtime. Therefore, we need to use a Graphics Pipeline stage where we still have access to all four positions of the vertices of the quadrilateral.

Hardware-accelerated implementation

The hardware-accelerated implementation of the algebraic solution presented in the previous section can be achieved using three different Graphics Pipeline stages: Geometry Shader, Tessellation Shader, and Mesh Shader. This three-way implementation demonstrates that the technique is feasible on GPU hardware released as far back as 2008. The following source code uses the GLSL shading language.

The table below describes what API versions and extensions are necessary to implement the presented technique:

The implementation of the presented method requires Barycentric coordinates of the currently processed fragment in the Fragment Shader. This is possible when the following extensions are available:

• Vulkan® VK_KHR_fragment_shader_barycentric extension.

• DirectX® 12 Shader Model 6.1.

• OpenGL® 4 GL_NV_fragment_shader_barycentric or GL_AMD_shader_explicit_vertex_parameter.

If the GPU does not support the above extension, there is a way to generate them in the Geometry, Tessellation, or Mesh Shader stage.

The implementation of the new method presented in this article can be divided into the following steps:

1. Delivering Quadrilateral Primitives to the Graphics Pipeline.

   Quadrilateral primitives are not available in modern graphics APIs (excluding the legacy OpenGL GL_QUADS and GL_QUAD_STRIP primitive types). To support quadrilateral-like primitives, use the following methods:

   • Geometry Shader: With the Geometry Shader stage, adjacency primitives can be used to mimic quadrilateral primitives:

     • OpenGL® GL_LINES_ADJACENCY

     • Vulkan® VK_PRIMITIVE_TOPOLOGY_LINE_LIST_WITH_ADJACENCY

     • DirectX® 12/11 D3D_PRIMITIVE_TOPOLOGY_LINELIST_ADJ

   • With Tessellation stages, a new Patch primitive was introduced. In this topology, the number of vertices defined in the meshes is not fixed, but the so-called Control Points range from 1 to 32 per Patch. These Control Points can be treated as Control Points of Bézier patches, but they can also be treated as vertices of a polygon. The implementation of the solution presented in this article uses a Patch with 4 Control Points and treats them as vertices of a quadrilateral.

   • Mesh Shader: Using Mesh Shaders, there are many ways to introduce a quadrilateral primitive into the Graphics Pipeline. This stage was developed to provide more flexibility to the user. Therefore, it is up to the user to generate meshlets that will contain quadrilaterals in some way (this is not within the scope of this article).

2. The next step is to calculate the constants A_u, A_v, B (using one of the steps above), which will then be used in the Fragment Shader. We need to calculate these constants in one of the above pipeline stages because, in the Fragment Shader, we do not have access to the four attributes of the quadrilateral vertices (the Barycentric coordinate extension, if present, only gives us access to the three attributes of the triangle vertices). By analyzing the algebraic solution of the presented method, we can see that these constants depend on the positions of the vertices of the quadrilateral.

The following function needs to be executed in the Geometry Shader, Tessellation Evaluation Shader, or Mesh Shader.

    vec3 calculateConstants(vec3 v0, vec3 v1, vec3 v2, vec3 v3)
    {
        const vec3 v10 = v1 - v0;
        const vec3 v30 = v3 - v0;
        const vec3 crossv10v30 = cross(v10, v30);
        const vec3 n = normalize(crossv10v30);
        const vec3 d = v0 - v1 + v2 - v3;
        float A_u = dot(n, cross(d, v30));
        float A_v = dot(n, cross(d, v10));
        const float B =  dot(n, crossv10v30);
        A_u = abs(A_u) < 0.00001 ? 0.0 : A_u;
        A_v = abs(A_v) < 0.00001 ? 0.0 : A_v;
        return vec3(A_u, A_v, B);
    }

Where v0, v1, v2, v3 are the Euclidean coordinates of the quadrilateral vertices in the following order:

The calculated values of A_u, A_v and B are then passed to the Fragment Shader.

2. Passing vertex attributes to the Fragment Shader.

For Bilinear interpolation to work in a Fragment Shader, we need to be able to fetch all four vertex attributes (such as texture coordinates, vertex normal, and vertex color) in it.

• If the Barycentric coordinate extension is available, we have access to three of these attributes, so another stage must export the missing vertex attribute of the quadrilateral. The following code can be used in the Fragment Shader:

layout (location = 0) pervertexEXT in block
{
	vec3 Normal;
	flat vec3 NormalExtra;
	vec2 Texcoord;
	flat vec2 TexcoordExtra;
	vec4 Color;
	flat vec4 ColorExtra;
	flat vec3 BaryConstants;
} In[];

• If the Barycentric coordinate extension is not available, all four vertex attributes need to be exported as a four-element array, with a variable calculated to imitate the Barycentric Coordinate. The following code can be used in the Fragment Shader:

layout(location = 0) in block
{
    flat vec3 Normal[4];
    flat vec2 Texcoord[4];
	flat vec4 Color[4];
	flat vec3 BaryConstants;
	vec3 BarycentricCoordinate;	
} In;

3. Orders of Barycentric coordinates in the triangles formed when dividing the quadrilateral.

When using Barycentric coordinates for interpolation of vertex attributes, the order in which they are defined in the triangle doesn’t matter. However, in the implementation method described in this article, the specific order is important. This order can be visually represented as follows:

Where red color represents the influence of the \alpha Barycentric coordinate, green represents \beta, and blue represents \gamma.

This order of Barycentric coordinates can be supplied by hardware when the proper extension is available, or these data can be generated in the Tessellation or Geometry Shader stage.

4. Calculate u, v Bilinear coefficients.

In the Fragment Shader, the calculation of Bilinear interpolation coefficients u and v can be done by calling the following function.

    vec2 calculateBilinear(float beta, float gamma, float A_u, float A_v, float B)
    {
      const float B_u = B;
      const float B_v = -B;
      const float B2 = beta * A_v + gamma * A_u;
      const float t_u = B_u - B2;
      const float t_v = B_v - B2;
      const float b_u = t_u / A_u * 0.5;
      const float b_v = t_v / A_v * 0.5;
      const float c_u = gamma * B_v;
      const float c_v = beta * B_u;
      const float d_u = b_u * b_u - c_u / A_u;
      const float d_v = b_v * b_v - c_v / A_v;
      const float s_u = A_u >= 0.0 ? -1.0 : 1.0;
      const float s_v = A_v >= 0.0 ? 1.0 : -1.0;
      const float u = A_u != 0.0 ? (-b_u - s_u * sqrt(d_u)) : (-c_u / t_u);
      const float v = A_v != 0.0 ? (-b_v - s_v * sqrt(d_v)) : (-c_v / t_v);
      return vec2(u, v);
    }

Where beta is \beta Barycentric coordinate, gamma is \gamma Barycentric coordinate, A_u, A_v, B are constants calculated in previous Graphics Pipeline Stage.

5. Bilinear interpolation of vertices attributes.

When Bilinear interpolation coefficients u and v are calculated, the interpolation of vertex attributes can be achieved using functions.

vec4 interpolateBilinear(vec4 v0, vec4 v1, vec4 v2, vec4 v3, vec2 uv)
{
    return mix(mix(v0, v1, uv.y), mix(v2, v3, uv.y), uv.x); 
}
vec3 interpolateBilinear(vec3 v0, vec3 v1, vec3 v2, vec3 v3, vec2 uv)
{
    return mix(mix(v0, v1, uv.y), mix(v2, v3, uv.y), uv.x); 
}
vec2 interpolateBilinear(vec2 v0, vec2 v1, vec2 v2, vec2 v3, vec2 uv)
{
    return mix(mix(v0, v1, uv.y), mix(v2, v3, uv.y), uv.x); 
}
float interpolateBilinear(float v0, float v1, float v2, float v3, vec2 uv)
{
    return mix(mix(v0, v1, uv.y), mix(v2, v3, uv.y), uv.x); 
}

Where v0, v1, v2, v3 are the quadrilateral vertex attributes.

vec4 calculateConstants(vec4 A, vec4 B, vec4 C, vec4 D)
{
    return A - B + C - D;
}
vec3 calculateConstants(vec3 A, vec3 B, vec3 C, vec3 D)
{
    return A - B + C - D;
}
vec2 calculateConstants(vec2 A, vec2 B, vec2 C, vec2 D)
{
    return A - B + C - D;
}
float calculateConstants(float A, float B, float C, float D)
{
    return A - B + C - D;
}

vec4 interpolateBilinear(vec4 attribute, vec2 barycentric, vec4 extraVal)
{
    return attribute + barycentric.x * barycentric.y * extraVal;
}
vec3 interpolateBilinear(vec3 attribute, vec2 barycentric, vec3 extraVal)
{
    return attribute + barycentric.x * barycentric.y * extraVal;
}
vec2 interpolateBilinear(vec2 attribute, vec2 barycentric, vec2 extraVal)
{
    return attribute + barycentric.x * barycentric.y * extraVal;
}
float interpolateBilinear(float attribute, vec2 barycentric, float extraVal)
{
    return attribute + barycentric.x * barycentric.y * extraVal;
}

Results

The output of the algebraic solution presented in this article remains consistent across different implementations. Depending on the type of vertex attributes, the results can be summarized in the following table.