GLSL Lens Distortion with Polynomial Radial Warps

This article explains how the Brown-Conrady polynomial radial distortion model is implemented in GLSL fragment shaders to evaluate and render optical barrel and pincushion lens distortions. By modifying normalized texture coordinates based on their radial distance from the optical center, developers can efficiently simulate lens imperfections, correct camera footage, or apply post-processing visual effects in real time.

The Brown-Conrady Algorithm

The primary algorithm used to model optical radial distortion is the Brown-Conrady model (often referred to simply as the polynomial radial distortion model). In physical lenses, light deviates symmetrically from the optical axis, causing straight lines to curve inward (pincushion distortion) or outward (barrel distortion).

The algorithm evaluates radial displacement using an even-powered polynomial series applied to the distance \(r\) from the optical center:

\[r_d = r \cdot (1 + k_1 r^2 + k_2 r^4 + k_3 r^6)\]

Where:

• \(r\) is the undistorted Euclidean distance from the lens center: \(r = \sqrt{x^2 + y^2}\).
• \(k_1, k_2, k_3\) are the radial distortion coefficients.
• \(r_d\) represents the resulting distorted radius.

When \(k_1 < 0\), the algorithm compresses coordinates at the perimeter, generating barrel distortion. When \(k_1 > 0\), outer coordinates are stretched outward, generating pincushion distortion.

Implementing the Algorithm in GLSL

In a GLSL fragment shader, applying this distortion requires remapping the input texture coordinate (uv) before sampling the scene texture.

1. Coordinate Normalization

Texture coordinates range from \([0.0, 1.0]\). The algorithm requires coordinates centered at the optical center (typically \((0.0, 0.0)\)) and adjusted for the viewport's aspect ratio:

vec2 uv = gl_FragCoord.xy / u_resolution.xy;
vec2 centeredCoord = uv - 0.5;
centeredCoord.x *= u_resolution.x / u_resolution.y; // Correct for non-square aspect ratio

2. Radius and Polynomial Evaluation

Compute the squared distance (\(r^2\)) to optimize performance, avoiding unnecessary square root operations where possible:

float r2 = dot(centeredCoord, centeredCoord);
float r4 = r2 * r2;

// Polynomial warp factor using coefficients k1 and k2
float distortionFactor = 1.0 + u_k1 * r2 + u_k2 * r4;

3. Remapping and Texture Sampling

Multiply the centered coordinates by the distortion factor, invert the aspect ratio correction, and shift the origin back to standard UV space:

vec2 distortedCoord = centeredCoord * distortionFactor;
distortedCoord.x /= u_resolution.x / u_resolution.y; // Restore aspect ratio
vec2 finalUV = distortedCoord + 0.5;

Complete Fragment Shader Example

#version 330 core

uniform sampler2D u_texture;
uniform vec2 u_resolution;
uniform float u_k1; // Primary distortion coefficient
uniform float u_k2; // Secondary distortion coefficient

in vec2 v_texCoord;
out vec4 fragColor;

void main() {
    // Center UV coordinates to [-0.5, 0.5]
    vec2 coord = v_texCoord - 0.5;
    
    // Account for aspect ratio
    float aspect = u_resolution.x / u_resolution.y;
    coord.x *= aspect;

    // Calculate radial distance squared
    float r2 = dot(coord, coord);

    // Evaluate polynomial radial warp
    float f = 1.0 + r2 * (u_k1 + u_k2 * r2);
    vec2 distorted = coord * f;

    // Return to [0, 1] UV space
    distorted.x /= aspect;
    vec2 finalUV = distorted + 0.5;

    // Discard fragments warped outside texture bounds
    if (finalUV.x < 0.0 || finalUV.x > 1.0 || finalUV.y < 0.0 || finalUV.y > 1.0) {
        fragColor = vec4(0.0, 0.0, 0.0, 1.0);
    } else {
        fragColor = texture(u_texture, finalUV);
    }
}

Forward vs. Inverse Mapping

• Forward Mapping (Distortion Simulation): Applying positive or negative coefficients directly to screen UVs calculates where undistorted source pixels map to on screen. This is computationally fast and ideal for post-processing effects such as CRT simulation, fish-eye lenses, or virtual reality (VR) headset compensation.
• Inverse Mapping (Undistortion/Correction): When correcting pre-existing physical camera distortion, an inverse polynomial (or a numerical solver like Newton-Raphson) must be used to calculate the undistorted coordinates from the distorted sensor inputs. In real-time GLSL applications, approximating the inverse with inverted polynomial coefficients is typically sufficient.