How are 2D images formed on an image sensor from a real 3D scene? If you simply place an image sensor in front of the scene, the pixels on the sensor receive light from multiple points on the scene. The image will simply become a blurry shadow of the actual scene.

To produce a perfectly clear image, the light needs to be focused onto the image sensor such that the each pixel on the sensor corresponds to a single point point in 3D space. This can be done with a pinhole camera.

Let us simulate such a camera to see how it behaves.

In [1]:
from __future__ import annotations
import halide as hl
from halide_book.interactive import Slider, interact
from halide_book.dataclass import Ray, TypedFunc, Vector2, Vector3, getitem_func
from halide_book.variables import col, row, channel, img
from dataclasses import dataclass
import numpy as np

As in most ray or path tracers, we start at the end: Where the ray hits our camera. This is because we do not want to spend time on computing rays that do not hit the camera. So it is better to backtrace from the camera, through the scene and finally to the light source.

The pixel coordinate is given by default in Halide as the Var arguments to the left side. So the first thing we want to compute is the normalized camera pixel coordinate.

In the OpenCV calibration model, these are found using the intrinsic matrix of the camera $K$. This is typically defined as

$$ K = \begin{bmatrix} f_x & 0 & c_x \\ 0 & f_y & c_y \\ 0 & 0 & 1 \end{bmatrix} $$

The normalized device coordinates $u$ and $v$ are found as

$$u = \frac{x - c_x}{f_x}$$ $$v = \frac{y - c_y}{f_y}$$

We can implement this as a Func in Halide:

In [2]:
def create_ndc(
    camera_width: int,
    camera_height: int,
) -> hl.Func:

    c_x = camera_width / 2
    c_y = camera_height / 2

    f_x = camera_width
    f_y = camera_width

    u = (col - c_x) / f_x
    v = (row - c_y) / f_y

    ndc = hl.Func("ndc")
    ndc[col, row] = Vector2(x=u, y=v)
    return ndc

To visualize the coordinates, we can create a new function that stores the value of u in the red channel and the value of v in the blue channel:

In [3]:
if __name__ == "__main__":
    camera_width = 480
    camera_height = 320
    ndc = create_ndc(camera_width=camera_width, camera_height=camera_height)

    ndc_image = hl.Func("ndc_image")
    ndc_image[col, row, img] = 0.0
    ndc_image[col, row, 0] = 0.5 + 0.5 * ndc[col, row][0]
    ndc_image[col, row, 2] = 0.5 + 0.5 * ndc[col, row][1]

    interact(ndc_image, size=[camera_width, camera_height, 3])

With this in place, we are close to creating a ray in a pinhole camera model. But before we do so, we need something to look at. Let us define a basic checkerboard as

In [4]:
@dataclass
class Checkerboard:
    checkerboard: hl.Func
    light_checker_slider: Slider
    dark_checker_slider: Slider

    @staticmethod
    def create() -> Checkerboard:
        light_checker = hl.Param(hl.Float(32), name="light_checker", value=0.8)
        dark_checker = hl.Param(hl.Float(32), name="dark_checker", value=0.2)

        checkerboard = hl.Func("checkerboard")
        checkerboard[col, row, channel] = hl.select(
            ((col % 2) ^ (row % 2)) == 0,
            light_checker,
            dark_checker,
        )

        light_checker_slider = Slider(param=light_checker, min_value=0.0, max_value=1.0, step=0.01)
        dark_checker_slider = Slider(param=dark_checker, min_value=0.0, max_value=1.0, step=0.01)
        return Checkerboard(
            checkerboard=checkerboard,
            light_checker_slider=light_checker_slider,
            dark_checker_slider=dark_checker_slider,
        )


if __name__ == "__main__":
    checkerboard = Checkerboard.create()
    interact(
        checkerboard.checkerboard,
        size=[10, 10, 3],
        controls=[
            checkerboard.dark_checker_slider,
            checkerboard.light_checker_slider,
        ],
    )

dark_checker: 0.20000000298023224

light_checker: 0.800000011920929

Next, we define our ray. This can be done by assuming the camera center is in the origin and use the above coordinate and $z=1$ as the direction of the ray.

In [5]:
def create_ndc_ray(camera_width: int, camera_height: int) -> Ray:
    ndc = create_ndc(
        camera_width=camera_width,
        camera_height=camera_height,
    )
    ndc_vector = Vector2.from_tuple(ndc[col, row])
    direction = Vector3(
        x=ndc_vector.x,
        y=ndc_vector.y,
        z=1.0,
    )
    origin = Vector3(
        x=0.0,
        y=0.0,
        z=0.0,
    )
    return Ray(origin=origin, direction=direction)


if __name__ == "__main__":
    ray = create_ndc_ray(
        camera_height=camera_height,
        camera_width=camera_width,
    )

Next, we trace this ray and calculate the intersection with the plane:

In [6]:
@dataclass
class Plane:
    normal: Vector3
    right: Vector3
    origin: Vector3

    @property
    def up(self) -> Vector3:
        return self.normal.cross(self.right)


@dataclass
class Intersection:
    point: hl.Func
    valid: hl.Func

    def __getitem__(self, *coords) -> Intersection:
        return getitem_func(self, *coords)


def plane_intersection(
    ray: Ray,
    plane: Plane,
) -> Intersection:
    denom = plane.normal.normalized().dot(ray.direction)
    t = plane.normal.dot(plane.origin - ray.origin) / denom
    point = hl.Func("point")
    point[col, row] = ray.origin + t * ray.direction
    valid = hl.Func("valid")
    valid[col, row] = (denom > hl.f32(1e-6)) & (t > 0)

    return Intersection(
        point=point,
        valid=valid,
    )


if __name__ == "__main__":
    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),
    )
    intersection = plane_intersection(
        ray=ray,
        plane=plane,
    )

We can visualize this intersection as a color based on where in the plane-local coordinate system we hit. To do so, we divide the coordinate by some number, which we can vary with a slider.

In [7]:
if __name__ == "__main__":
    plane_extent = hl.Param(hl.Float(32), name="plane_extent", value=1000.0)
    plane_extent_slider = Slider(plane_extent, min_value=100, max_value=2000, step=1)

    point_image = hl.Func("point_image")
    point_image[col, row, img] = 0.0
    point_image[col, row, 0] = intersection.point[col, row][0] / plane_extent
    point_image[col, row, 1] = intersection.point[col, row][1] / plane_extent
    point_image[col, row, 2] = intersection.point[col, row][2] / plane_extent

    interact(point_image, size=[camera_width, camera_height, 3], controls=[plane_extent_slider])

plane_extent: 1000.0

Now we can find the texture coordinate on the checkerboard:

In [8]:
@dataclass
class CheckerboardTexture:
    checker_size: hl.Param
    checker_size_slider: Slider

    @staticmethod
    def create() -> CheckerboardTexture:
        checker_size = hl.Param(hl.Float(32), "checker_size", 100.0)

        checker_size_slider = Slider(checker_size, 10, 200, 1)

        return CheckerboardTexture(
            checker_size=checker_size,
            checker_size_slider=checker_size_slider,
        )

    def sample(
        self,
        checkerboard: Checkerboard,
        point: hl.Func,
        plane: Plane,
    ) -> hl.Func:
        intersection_vector = Vector3.from_tuple(point[col, row]) - plane.origin
        plane.pixel = Vector2(
            (intersection_vector).dot(plane.right),
            (intersection_vector).dot(plane.up),
        )

        texture_coordinate = Vector2(
            x=plane.pixel.x / self.checker_size,
            y=plane.pixel.y / self.checker_size,
        )

        texture = hl.Func("texture")
        texture[col, row, channel] = checkerboard.checkerboard[
            hl.clamp(hl.cast(hl.Int(32), hl.floor(texture_coordinate.x)), -100, 100),
            hl.clamp(hl.cast(hl.Int(32), hl.floor(texture_coordinate.y)), -100, 100),
            channel,
        ]
        return texture

Calling this function gives us a new Func that we can visualize with a slider for checker size. Let us make a function that we can reuse later that puts this all together:

In [9]:
@dataclass
class CheckerboardRenderer:
    checkerboard: Checkerboard
    texture: CheckerboardTexture

    def sample(self, ray: Ray, plane: Plane) -> Func:
        plane = Plane(
            normal=Vector3(
                x=0.25,
                y=0.0,
                z=0.35,
            ).normalized(),
            right=Vector3(
                x=0.35,
                y=0.25,
                z=0.0,
            ).normalized(),
            origin=Vector3(
                x=0.0,
                y=0.0,
                z=1000.0,
            ),
        )
        intersection = plane_intersection(ray=ray, plane=plane)
        return self.texture.sample(
            checkerboard=self.checkerboard,
            point=intersection.point,
            plane=plane,
        )


if __name__ == "__main__":
    ray = create_ndc_ray(
        camera_height=camera_height,
        camera_width=camera_width,
    )
    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=[
            renderer.texture.checker_size_slider,
            renderer.checkerboard.dark_checker_slider,
            renderer.checkerboard.light_checker_slider,
        ],
    )

checker_size: 100.0

dark_checker: 0.20000000298023224

light_checker: 0.800000011920929

So that was quite a bit of work just to render a checkerboard, but it sets us up for exploring quite a few more things in future sections. For instance, we can now start distorting the direction of the rays based on a lens model that is more advanced than the basic pinhole camera.

But before we move on, a small improvement we can make is to antialias by taking four samples for a given ray and averaging the result:

In [10]:
if __name__ == "__main__":
    ray = create_ndc_ray(
        camera_height=camera_height,
        camera_width=camera_width,
    )
    ray_func = TypedFunc("ray")
    ray_func[col, row] = ray
    ray_00 = ray_func[col, row]
    ray_01 = ray_func[col, row + 1]
    ray_10 = ray_func[col + 1, row]
    ray_11 = ray_func[col + 1, row + 1]
    checkerboard = Checkerboard.create()
    checkerboard_texture = CheckerboardTexture.create()
    renderer = CheckerboardRenderer(
        checkerboard=checkerboard,
        texture=checkerboard_texture,
    )
    output = hl.Func("output")
    output[col, row, img] = (
        renderer.sample(ray=ray_00, plane=plane)[col, row, img]
        + renderer.sample(ray=ray_01, plane=plane)[col, row, img]
        + renderer.sample(ray=ray_10, plane=plane)[col, row, img]
        + renderer.sample(ray=ray_11, plane=plane)[col, row, img]
    ) / 4.0
    interact(
        output,
        size=[camera_width, camera_height, 3],
        controls=[
            renderer.texture.checker_size_slider,
            renderer.checkerboard.dark_checker_slider,
            renderer.checkerboard.light_checker_slider,
        ],
    )

checker_size: 100.0

dark_checker: 0.20000000298023224

light_checker: 0.800000011920929