Distortion
In the previous section, we introduced a pinhole camera model. However, most real cameras are made with lenses that introduce some distortion.
To add distortion to the pinhole model, we can extend it such that the normalized device coordinates $u$ and $v$ are defined by
$$ \begin{bmatrix} u \\ v \end{bmatrix} = \begin{bmatrix} f_x x'' + c_x \\ f_y y'' + c_y \end{bmatrix} $$
where $x''$ and $y''$ are the distortion coefficients. The distortion coefficients are functions of the $x'$ and $y'$ coordinates and determined by terms that we will introduce in this and following sections: Radial distortion, tangential distortion, and prism distortion.
We can extend the existing pinhole model in code by importing the relevant definitions from the previous section.
from __future__ import annotations
from dataclasses import dataclass
import halide as hl
from halide_book.interactive import Slider, interact
from halide_book.dataclass import Ray, Vector2, Vector3
from halide_book.pinhole import (
CheckerboardRenderer,
create_ndc_ray, Checkerboard, CheckerboardTexture, Plane,
)
from halide_book.variables import col, row, channel, img
First, let us visualize the pinhole model again without any distortion:
if __name__ == "__main__":
camera_width = 640
camera_height = 480
plane = Plane(
normal=Vector3(x=0.55, y=0.0, z=0.35).normalized(),
right=Vector3(x=0.55, y=0.35, z=0.0).normalized(),
origin=Vector3(x=0.0, y=0.0, z=1000.0),
)
ray = create_ndc_ray(camera_width=camera_width, camera_height=camera_height)
checkerboard = Checkerboard.create()
checkerboard_texture = CheckerboardTexture.create()
renderer = CheckerboardRenderer(
checkerboard=checkerboard,
texture=checkerboard_texture,
)
output = renderer.sample(ray=ray, plane=plane)
interact(
output,
size=[camera_width, camera_height, 3],
controls=[
checkerboard_texture.checker_size_slider,
renderer.checkerboard.dark_checker_slider,
renderer.checkerboard.light_checker_slider,
],
)
Next, we can start introducing distortion effects.
Commonly, a group of parameters define the radial distortion. This type of distortion makes straight lines look curved. A basic radial distortion definition is defined from three parameters, $k_1$, $k_2$, and $k_3$.
$$ \begin{bmatrix} x'' \\ y'' \end{bmatrix} = \begin{bmatrix} x' \left ( 1 + k_1 r^2 + k_2 r^4 + k_3 r^6 \right) \\ y' \left ( 1 + k_1 r^2 + k_2 r^4 + k_3 r^6 \right) \\ \end{bmatrix} $$
def make_slider(name: str):
param = hl.Param(hl.Float(32), name, 0.0)
return Slider(param, min_value=-1.0, max_value=1.0, step=0.01)
@dataclass
class RadialDistortionSimple:
k1: Slider
k2: Slider
k3: Slider
def apply(self, vector: Vector2) -> Vector2:
x = vector.x
y = vector.y
r2 = x**2 + y**2
r4 = r2 * r2
r6 = r2 * r4
radial_coefficient = 1.0 + self.k1.param * r2 + self.k2.param * r4 + self.k3.param * r6
return Vector2(
x * radial_coefficient,
y * radial_coefficient,
)
@staticmethod
def create() -> RadialDistortionSimple:
return RadialDistortionSimple(
k1=make_slider("k1"),
k2=make_slider("k2"),
k3=make_slider("k3"),
)
We can now visualize the effect of such a distortion. To do so, we need to take into account that the distortion applies to the undistorted pinhole image after it has been created. Distortion would take us from an undistorted pixel in an image to a distorted pixel on the sensor. But we already know the pixel on the sensor. That is the pixel we are starting from. Instead, we need the inverse of the distortion.
Since the distortion is non-linear, its inverse is not trivial. It can, however, be found with an iterative method.
As long as the parameters are within reasonable bounds, we can use something called fixed-point iteration. This method is not guaranteed to converge, but is pretty fast to compute when it does work.
The method works by assuming an initial guess for the undistorted pixel. If the distortion is not large, we can set this to be the same as the distorted.
$$x ~= x'$$
Next, we assume a fixed-point function on the form
$$x = x + x' - \mathrm{distort}(x)$$
If this converges, it means that the difference between the distorted $x'$ is approaching the estimate $\mathrm{distort}(x)$. And if that is true, $x$ must be the undistorted value.
In code, we can apply this with a small number of iterations:
current = ndc_vector
for i in range(4):
current = current + ndc_vector - full_distortion(current)
distorted = Vector3.from_vector2(current, z=1.0)
Note again that the above will not converge if the distortion is large. If that is the case, we might need to apply something like gradient descent or Newton's method instead, with a small step size. Whether the above method works or not can be checked by a round-trip test for a given set of parameters. The result of running distortion and then undistortion should result in the same value. I will skip that for now, but if you are using this in production, I would strongly suggest adding such a test for any given set of parameters.
By implementing undistortion, we can render an image that corresponds to the distortion model. This is what we do below.
Try to vary the parameters. Note that if you keep $k_2 = k_3 = 0$, then a value of $k_1 < 0.0$ will give a so-called barrel distortion, while a value of $k_1 > 0.0$ gives a pincushion distortion.
if __name__ == "__main__":
ray = create_ndc_ray(camera_height=camera_height, camera_width=camera_width)
pixel = Vector2(x=ray.direction.x, y=ray.direction.y)
undistorted = Vector2(x=ray.direction.x, y=ray.direction.y)
distortion = RadialDistortionSimple.create()
for i in range(6):
undistorted = pixel + undistorted - distortion.apply(undistorted)
undistorted = Vector3.from_vector2(undistorted, z=1.0)
undistorted_ray = Ray(
origin=ray.origin,
direction=undistorted,
)
output = renderer.sample(ray=undistorted_ray, plane=plane)
interact(
output,
size=[camera_width, camera_height, 3],
controls=[
checkerboard_texture.checker_size_slider,
renderer.checkerboard.dark_checker_slider,
renderer.checkerboard.light_checker_slider,
distortion.k1,
distortion.k2,
distortion.k3,
],
)