4 Images as Functions on Domains

A digital image is a finite array, but many of the most important ideas in computer vision become clean only after we view the array as samples of an underlying function. This chapter develops the bridge between continuous images, discrete samples, interpolation, finite differences, aliasing, random fields, and differentiable image transformations. The goal is to give applied mathematics and statistics students a precise language for moving between calculus, probability, numerical analysis, and modern AI vision models.

4.1 Learning Goals

After reading this chapter, students should be able to:

model a grayscale or color image as a function $I: \Omega \to \mathbb{R}^C$ on a spatial domain;
distinguish the continuous image, the sampled image, and the reconstructed/interpolated image;
write sampling as an operator and understand discretization error;
explain nearest-neighbor, bilinear, and smoother interpolation from a mathematical viewpoint;
compute finite-difference approximations of image gradients;
explain aliasing using sinusoidal examples and the Nyquist idea;
describe geometric transformations as function composition and pullback operations;
connect classical image models with differentiable preprocessing and AI vision pipelines.

4.2 3.1 Why Images-as-Functions Matters

Chapter 2 introduced pixels, coordinates, and color spaces. That chapter emphasized that a stored image is usually an array

\[ I \in \mathbb{R}^{H \times W \times C}. \]

This array notation is necessary, but it is not always sufficient. When we resize an image, rotate it, blur it, compute a gradient, warp it under a camera model, or feed it to a neural network with differentiable data augmentation, we are implicitly treating the image as a function on a spatial domain.

A clean mathematical description is

\[ I: \Omega \subset \mathbb{R}^2 \longrightarrow \mathbb{R}^C, \]

where $\Omega$ is an image domain and $C$ is the number of channels. For a grayscale image, $C=1$. For an RGB image, $C=3$. For a feature map inside a convolutional neural network, $C$ may be 64, 256, 768, or larger.

The central idea of this chapter is the distinction

\[ \text{continuous image} \quad \longrightarrow \quad \text{sampled image} \quad \longrightarrow \quad \text{reconstructed image}. \]

A digital image is not just a function and not just an array. It is usually a sampled and quantized representation of an underlying signal, together with assumptions about how to reconstruct or interpolate it.

This viewpoint is useful across the whole book:

Fourier analysis studies the frequency content of image functions.
Filtering applies linear operators to image functions.
Edges are locations where derivatives or finite differences are large.
Camera geometry maps scene functions to image functions.
Segmentation partitions the image domain $\Omega$.
Deep learning repeatedly transforms sampled functions into new feature functions.
Diffusion models learn probability distributions over images as high-dimensional random functions.

The mathematical language of functions also prevents many practical errors. For example, a rotation by $13^\circ$ does not map integer pixel centers to integer pixel centers. A model must decide how to assign intensity values at non-integer coordinates. That decision is interpolation, and it has mathematical consequences.

4.3 3.2 Continuous Images and Spatial Domains

Let $\Omega \subset \mathbb{R}^2$ be a rectangular domain,

\[ \Omega = [0,W] \times [0,H]. \]

A grayscale image is a scalar-valued function

\[ I: \Omega \to \mathbb{R}. \]

A color image is a vector-valued function

\[ I: \Omega \to \mathbb{R}^3, \qquad I(x,y) = \begin{bmatrix} R(x,y) \\ G(x,y) \\ B(x,y) \end{bmatrix}. \]

More generally,

\[ I(x,y) \in \mathbb{R}^C. \]

Here $C$ is the number of channels. In a raw camera image, the channels may represent physical sensor responses. In a deep network, the channels may represent learned feature responses.

4.3.1 Function spaces

To do analysis, we often place images in function spaces. Common choices include:

\[ L^2(\Omega) = \left\{ I: \Omega \to \mathbb{R} : \int_{\Omega} |I(x,y)|^2\,dx\,dy < \infty \right\}, \]

and, for differentiable images,

\[ C^1(\Omega) = \left\{ I: \Omega \to \mathbb{R} : I_x \text{ and } I_y \text{ exist and are continuous} \right\}. \]

In practice, real images may contain discontinuities at object boundaries. Therefore $C^1(\Omega)$ is often too strong. A more realistic mathematical model says that images are smooth in regions but may have jumps along curves.

Piecewise-smooth image model. A useful idealization is that an image is smooth inside objects and has discontinuities across object boundaries. Edges are then not mysterious artifacts; they are the visible traces of discontinuities or rapid changes in the image function.

This model helps explain why gradients and edges are central in computer vision. If $I$ is differentiable at $(x,y)$, then its first-order local approximation is

\[ I(x+\Delta x, y+\Delta y) \approx I(x,y) + I_x(x,y)\Delta x + I_y(x,y)\Delta y. \]

The vector

\[ \nabla I(x,y) = \begin{bmatrix} I_x(x,y) \\ I_y(x,y) \end{bmatrix} \]

points in the direction of greatest local increase of intensity. Its magnitude $\|\nabla I(x,y)\|$ measures local contrast.

4.4 3.3 Sampling Operators

A camera does not store the whole function $I$. It samples it on a finite grid.

Let

\[ \Lambda_{H,W} = \{0,1,\ldots,H-1\} \times \{0,1,\ldots,W-1\} \]

be the set of array indices. Under the pixel-center convention, the sample location for row $i$ and column $j$ is

\[ (x_j,y_i)=\left(j+\frac12,\ i+\frac12\right). \]

The sampling operator $S$ maps a continuous image to a discrete image:

\[ S(I)[i,j] = I\left(j+\frac12, i+\frac12\right). \]

If the image has channels, then

\[ S(I)[i,j,c] = I_c\left(j+\frac12, i+\frac12\right). \]

This formula is simple, but it hides an important assumption: the pixel value equals the value of the image function at the pixel center. A physically more realistic camera model may average light over a pixel cell:

\[ S_{\text{avg}}(I)[i,j] = \frac{1}{|P_{ij}|} \int_{P_{ij}} I(x,y)\,dx\,dy, \]

where

\[ P_{ij} = [j,j+1] \times [i,i+1]. \]

The point-sample model and the cell-average model lead to different behavior, especially when the image has high-frequency texture.

Sampling is not a harmless bookkeeping step. If a continuous image changes too rapidly relative to the pixel grid, the sampled image may contain misleading low-frequency patterns. This is aliasing.

The interactive figure above shows a continuous image surface sampled on a finite grid. Increasing the spatial frequency of the underlying function would require a finer sampling grid to preserve the visible structure.

4.5 3.4 Discrete Images as Functions on a Grid

Once sampled, an image is a function on a finite set:

\[ I_h: \Lambda_{H,W} \to \mathbb{R}^C. \]

Here $h$ represents grid spacing. In many image arrays, the grid spacing is normalized to $h=1$, but for numerical analysis it is useful to write

\[ x_j = x_0 + jh_x, \qquad y_i = y_0 + ih_y. \]

Then the discrete image is

\[ I_h[i,j] = I(x_j,y_i). \]

The continuous function $I$ and the discrete function $I_h$ are mathematically different objects. The continuous function has values at every point in $\Omega$. The discrete image only has values at finitely many grid locations.

A common source of confusion is to write $I(x,y)$ for both. In this book we will often use:

\[ I(x,y) \quad \text{for a continuous image}, \]

and

\[ I[i,j] \quad \text{for a sampled image array}. \]

When we move between the two notations, we are making a modeling choice.

4.6 3.5 Reconstruction and Interpolation

Suppose we know only the discrete samples $I[i,j]$. What value should the image have at a non-grid point $(x,y)$?

This is a reconstruction problem. A reconstruction operator $R$ maps a discrete image back to a continuous or piecewise-continuous image:

\[ R(I_h)(x,y) \approx I(x,y). \]

Different reconstruction rules produce different images.

4.6.1 Nearest-neighbor interpolation

Nearest-neighbor interpolation chooses the value of the closest pixel center:

\[ R_{\text{nn}}(I_h)(x,y) = I[i^*,j^*], \]

where

\[ (i^*,j^*) = \arg\min_{(i,j) \in \Lambda_{H,W}} \left\|(x,y)-(x_j,y_i)\right\|_2. \]

This method is simple and preserves original pixel values, but it creates blocky artifacts.

4.6.2 Bilinear interpolation

Bilinear interpolation uses the four neighboring pixel values. Let

\[ x \in [j,j+1], \qquad y \in [i,i+1], \]

and define

\[ a = x-j, \qquad b = y-i, \qquad 0 \le a,b \le 1. \]

For a scalar image, the bilinear interpolant is

\[ \begin{aligned} R_{\text{bilin}}(I_h)(x,y) =&(1-a)(1-b)I[i,j] + a(1-b)I[i,j+1] \\ &+ (1-a)b I[i+1,j] + ab I[i+1,j+1]. \end{aligned} \]

This formula is a convex combination of the four corner values. The weights are nonnegative and add to one:

\[ (1-a)(1-b) + a(1-b) + (1-a)b + ab = 1. \]

Therefore, if the four corner values lie in $[0,1]$, the interpolated value also lies in $[0,1]$.

4.6.3 Interpolation as basis expansion

Interpolation can also be written as a basis expansion:

\[ R(I_h)(x,y) = \sum_{i,j} I[i,j] \varphi_{ij}(x,y), \]

where $\varphi_{ij}$ is a basis or kernel centered near pixel $(i,j)$.

For nearest-neighbor interpolation, the basis functions are piecewise constant indicator functions. For bilinear interpolation, the basis functions are tent functions. For smoother interpolation, one may use cubic kernels or sinc-based reconstructions.

This perspective is important because many image transformations in deep learning are differentiable only through the chosen interpolation formula. Bilinear interpolation is popular because it is simple, local, and piecewise differentiable with respect to both image values and coordinates.

4.7 3.6 Discretization Error

Let $I$ be a continuous image and $I_h=S(I)$ its sampled version. After reconstruction, we obtain $R(I_h)$. A natural error is

\[ \|I - R(S(I))\|. \]

The size of this error depends on:

the smoothness of $I$;
the grid spacing $h$;
the reconstruction method $R$;
the frequency content of $I$;
whether prefiltering was used before sampling.

For a smooth one-dimensional function $f$, linear interpolation has local error of order $h^2$ under suitable assumptions:

\[ |f(x) - R_{\text{linear}}(S(f))(x)| \le C h^2 \max_{\xi} |f''(\xi)|. \]

The two-dimensional bilinear case has analogous behavior when the second derivatives are controlled. The important qualitative lesson is:

Smoother images can be sampled and reconstructed more accurately than highly oscillatory images. High-frequency content requires finer sampling or appropriate low-pass filtering.

This principle appears again in Fourier analysis, scale-space theory, anti-aliased image resizing, and neural network preprocessing.

4.8 3.7 Aliasing and the Nyquist Idea

Aliasing occurs when different continuous signals produce the same samples on a grid.

Consider the one-dimensional sinusoid

\[ f(x) = \sin(2\pi \omega x), \]

sampled at integer points $x=n$. The samples are

\[ f[n] = \sin(2\pi \omega n). \]

If $\omega$ is replaced by $\omega + k$, where $k$ is an integer, then

\[ \sin(2\pi(\omega+k)n) = \sin(2\pi \omega n + 2\pi kn) = \sin(2\pi \omega n). \]

Thus many different continuous frequencies have exactly the same discrete samples.

The Nyquist idea says that, for uniform sampling with spacing $h$, frequencies above approximately

\[ \frac{1}{2h} \]

cannot be represented unambiguously without additional assumptions.

In images, aliasing appears as jagged edges, moire patterns, false textures, and strange artifacts after downsampling. A good image-resizing pipeline usually applies a low-pass filter before reducing resolution.

Anti-aliasing principle. Before sampling more coarsely, remove frequencies that the coarser grid cannot represent. In practice, this means blur before downsampling.

This principle is not only classical signal processing. It also affects modern AI systems. A classifier may change its prediction under small shifts or resizes partly because the preprocessing pipeline introduces aliasing.

4.9 3.8 Finite Differences and Image Derivatives

For a differentiable image $I(x,y)$, the gradient is

\[ \nabla I(x,y) = \begin{bmatrix} I_x(x,y) \\ I_y(x,y) \end{bmatrix}. \]

For a sampled image $I[i,j]$, derivatives must be approximated. Common finite-difference formulas include:

4.9.1 Forward differences

\[ D_x^+ I[i,j] = I[i,j+1] - I[i,j], \]

\[ D_y^+ I[i,j] = I[i+1,j] - I[i,j]. \]

4.9.2 Backward differences

\[ D_x^- I[i,j] = I[i,j] - I[i,j-1], \]

\[ D_y^- I[i,j] = I[i,j] - I[i-1,j]. \]

4.9.3 Central differences

\[ D_x^0 I[i,j] = \frac{I[i,j+1]-I[i,j-1]}{2}, \]

\[ D_y^0 I[i,j] = \frac{I[i+1,j]-I[i-1,j]}{2}. \]

With grid spacings $h_x$ and $h_y$, the denominators become $h_x$, $h_y$, $2h_x$, or $2h_y$ accordingly.

The discrete gradient magnitude is often approximated by

\[ \|\nabla_h I[i,j]\|_2 = \sqrt{(D_x I[i,j])^2 + (D_y I[i,j])^2}. \]

This discrete derivative viewpoint is foundational for edge detection, optical flow, variational segmentation, total variation denoising, and numerical PDE methods in imaging.

4.10 3.9 Difference Operators as Matrices

Finite differences can also be written as linear algebra. Suppose a one-dimensional signal is

\[ f = \begin{bmatrix} f_0 & f_1 & \cdots & f_{n-1} \end{bmatrix}^T. \]

A forward-difference operator can be represented by the matrix

\[ D = \begin{bmatrix} -1 & 1 & 0 & \cdots & 0 \\ 0 & -1 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 \end{bmatrix}. \]

Then

\[ Df = \begin{bmatrix} f_1-f_0 \\ f_2-f_1 \\ \vdots \\ f_{n-1}-f_{n-2} \end{bmatrix}. \]

For two-dimensional images, difference operators act along rows and columns. If the image is vectorized into a vector $u \in \mathbb{R}^{HW}$, then horizontal and vertical difference operators may be written as matrices

\[ D_x u, \qquad D_y u. \]

The discrete gradient operator is

\[ D u = \begin{bmatrix} D_x u \\ D_y u \end{bmatrix}. \]

This notation leads naturally to optimization problems such as

\[ \min_u \frac12\|u-y\|_2^2 + \lambda \|Du\|_1, \]

which is a simplified form of total-variation denoising. The first term keeps $u$ close to the observed image $y$, while the second term encourages piecewise-smooth structure.

4.11 3.10 Geometric Transformations as Function Composition

Let $T: \Omega \to \Omega$ be a spatial transformation. For example:

translation;
rotation;
scaling;
affine transformation;
projective transformation;
nonlinear warp.

If we want to transform an image, we need to define the transformed function. A common convention is the pullback:

\[ I_T(x,y) = I(T^{-1}(x,y)). \]

This means that the output value at location $(x,y)$ is obtained by looking backward to the corresponding location in the input image.

Why use $T^{-1}$ instead of $T$? Because we want to assign a value to every output pixel. For each output location, we ask: where did this point come from in the input image?

Backward warping. In image processing and differentiable vision layers, warping is often implemented by looping over output coordinates, mapping each one backward to an input coordinate, and then interpolating the input image at that coordinate.

For example, if $T$ is a translation by $(a,b)$,

\[ T(x,y)=(x+a,y+b), \]

then

\[ I_T(x,y) = I(x-a,y-b). \]

If $(x-a,y-b)$ is not a pixel center, interpolation is required.

This is the mathematical foundation behind image augmentation, spatial transformer networks, differentiable rendering, optical flow warping, and many registration algorithms.

4.12 3.11 Images as Random Functions

For statistics students, it is useful to view images as random objects.

A random image may be written as

\[ X: \Omega \times \Xi \to \mathbb{R}^C, \]

where $\Xi$ is a probability space. For each outcome $\xi \in \Xi$, the function

\[ X(\cdot,\cdot,\xi): \Omega \to \mathbb{R}^C \]

is one image.

A simple observation model is

\[ Y(x,y) = I(x,y) + \varepsilon(x,y), \]

where $\varepsilon$ is noise. In a discrete model,

\[ Y[i,j] = I[i,j] + \varepsilon[i,j]. \]

If the noise is independent Gaussian noise,

\[ \varepsilon[i,j] \sim N(0,\sigma^2), \]

then the negative log-likelihood is proportional to the squared error

\[ \sum_{i,j} (Y[i,j]-I[i,j])^2. \]

This is why least-squares losses appear naturally in denoising, reconstruction, registration, and many image-to-image prediction tasks.

More structured statistical models treat images as random fields with covariance

\[ K\bigl((x,y),(x',y')\bigr) = \operatorname{Cov}\bigl(X(x,y), X(x',y')\bigr). \]

A stationary covariance depends only on displacement:

\[ K\bigl((x,y),(x',y')\bigr) = k(x-x',y-y'). \]

This connects image modeling to Gaussian processes, Markov random fields, texture models, and modern generative models.

4.13 3.12 Function Viewpoint Inside Deep Networks

A convolutional neural network does not only process the original image. It builds a sequence of feature maps:

\[ F_0, F_1, F_2, \ldots, F_L. \]

The input might be

\[ F_0: \Omega_0 \to \mathbb{R}^3, \]

while an intermediate layer might be

\[ F_\ell: \Omega_\ell \to \mathbb{R}^{C_\ell}. \]

Each spatial location has a feature vector. Therefore a feature map is a vector-valued function on a discrete spatial domain.

A convolutional layer is approximately a local operator:

\[ F_{\ell+1}(p) = \sigma\left(\sum_{q \in N(p)} W(q-p)F_\ell(q)+b\right), \]

where:

$p$ is a spatial location;
$N(p)$ is a local neighborhood;
$W$ is a learned kernel;
$\sigma$ is a nonlinear activation.

A vision transformer uses a different representation. It partitions the image into patches and converts each patch into a token vector. This changes the spatial function into a sequence:

\[ I: \Omega \to \mathbb{R}^3 \quad \longrightarrow \quad (z_1,z_2,\ldots,z_N), \qquad z_i \in \mathbb{R}^d. \]

Even then, the construction begins with a sampled image function and a choice of patch geometry.

4.14 3.13 AI Component: Using Models to Reason About Image Functions

Large language models and multimodal AI systems can help students reason about computer vision, but only if the prompts are mathematically specific.

4.14.1 Prompt pattern 1: separate objects carefully

A useful study prompt is:

I have a continuous image $I:\Omega\to\mathbb{R}$, a sampled image $I_h[i,j]$, and a reconstructed image $R(I_h)(x,y)$. Explain the difference among these three objects, and give one example where confusing them causes an error in computer vision.

A weak answer will say only that one is continuous and one is discrete. A strong answer will mention sampling, interpolation, coordinate conventions, and possible aliasing.

4.14.2 Prompt pattern 2: ask for formula verification

Derive the bilinear interpolation formula from the requirement that the interpolated value is linear in $x$ for fixed $y$ and linear in $y$ for fixed $x$. Check that the four weights sum to one.

This prompt forces the AI tool to work through the mathematics instead of giving only implementation advice.

4.14.3 Prompt pattern 3: connect theory to code

Given a NumPy image array of shape (H, W), write a function for bilinear interpolation at non-integer coordinates (x, y). Explain how the formula corresponds to the mathematical weights $(1-a)(1-b)$, $a(1-b)$, $(1-a)b$, and $ab$.

This type of prompt is useful because it asks for both computation and mathematical interpretation.

4.14.4 Prompt pattern 4: audit preprocessing

A vision model resizes all input images from $1024\times 1024$ to $224\times 224$. Explain why this is a sampling problem. What can go wrong if no low-pass filtering is used before downsampling?

This prompt connects classical sampling theory with modern deep learning practice.

When using AI for study, do not ask only “explain interpolation.” Ask the model to specify the domain, the grid, the coordinates, the interpolation weights, and the assumptions. Computer vision errors often hide in these details.

4.15 3.14 Python Example: Sampling a Continuous Image

The following code constructs a synthetic continuous image function and samples it on a grid.

import numpy as np
import matplotlib.pyplot as plt

def continuous_image(x, y):
    """A smooth synthetic image function on [0, 1] x [0, 1]."""
    low_freq = 0.55 + 0.25 * np.sin(2 * np.pi * x) * np.cos(2 * np.pi * y)
    blob = 0.45 * np.exp(-45 * ((x - 0.65)**2 + (y - 0.35)**2))
    edge = 0.25 * (x + y > 1.05)
    return np.clip(low_freq + blob + edge, 0, 1)

H, W = 80, 120
ys = (np.arange(H) + 0.5) / H
xs = (np.arange(W) + 0.5) / W
X, Y = np.meshgrid(xs, ys)
I = continuous_image(X, Y)

plt.figure(figsize=(7, 4))
plt.imshow(I, cmap="gray", origin="upper", extent=[0, 1, 1, 0])
plt.colorbar(label="intensity")
plt.title("Sampled synthetic image")
plt.xlabel("x")
plt.ylabel("y")
plt.show()

Notice the use of pixel centers:

\[ x_j = \frac{j+1/2}{W}, \qquad y_i = \frac{i+1/2}{H}. \]

Changing this convention slightly may shift the image and affect geometric alignment.

4.16 3.15 Python Example: Bilinear Interpolation

The next example implements bilinear interpolation at a non-integer coordinate.

import numpy as np

def bilinear_interpolate(image, x, y):
    """
    Bilinear interpolation for a grayscale image.

    Parameters
    ----------
    image : ndarray of shape (H, W)
        Input image.
    x, y : float
        Continuous coordinates in array coordinate convention:
        x is column coordinate, y is row coordinate.

    Returns
    -------
    value : float
        Interpolated intensity.
    """
    H, W = image.shape

    # Clamp to valid interpolation range.
    x = np.clip(x, 0, W - 1)
    y = np.clip(y, 0, H - 1)

    j0 = int(np.floor(x))
    i0 = int(np.floor(y))
    j1 = min(j0 + 1, W - 1)
    i1 = min(i0 + 1, H - 1)

    a = x - j0
    b = y - i0

    return (
        (1 - a) * (1 - b) * image[i0, j0]
        + a * (1 - b) * image[i0, j1]
        + (1 - a) * b * image[i1, j0]
        + a * b * image[i1, j1]
    )

# Example
small = np.array([
    [0.0, 0.2, 0.4],
    [0.3, 0.6, 0.9],
    [0.5, 0.7, 1.0]
])

print(bilinear_interpolate(small, x=1.25, y=0.75))

A useful debugging check is to verify that if $(x,y)$ lands exactly on a pixel coordinate, the function returns that pixel value.

4.17 3.16 Python Example: Aliasing Under Downsampling

This example shows why high-frequency patterns need anti-aliasing before downsampling.

import numpy as np
import matplotlib.pyplot as plt

n_high = 1024
x_high = np.linspace(0, 1, n_high, endpoint=False)

freq = 37
signal_high = np.sin(2 * np.pi * freq * x_high)

# Downsample by taking every 16th point.
factor = 16
x_low = x_high[::factor]
signal_low = signal_high[::factor]

plt.figure(figsize=(8, 4))
plt.plot(x_high, signal_high, alpha=0.5, label="continuous-looking high-resolution signal")
plt.scatter(x_low, signal_low, label="samples after downsampling")
plt.legend()
plt.title("Aliasing from undersampling")
plt.xlabel("x")
plt.ylabel("signal")
plt.show()

A better pipeline would apply a low-pass filter before keeping every 16th sample.

4.18 3.17 Python Example: Finite-Difference Gradients

The following code computes a finite-difference gradient magnitude image.

import numpy as np
import matplotlib.pyplot as plt

H, W = 120, 160
y = np.linspace(0, 1, H)
x = np.linspace(0, 1, W)
X, Y = np.meshgrid(x, y)

I = 0.4 + 0.3 * X + 0.2 * np.exp(-60 * ((X - 0.65)**2 + (Y - 0.45)**2))
I += 0.35 * (X + Y > 1.15)
I = np.clip(I, 0, 1)

Dx = np.zeros_like(I)
Dy = np.zeros_like(I)

Dx[:, 1:-1] = (I[:, 2:] - I[:, :-2]) / 2
Dy[1:-1, :] = (I[2:, :] - I[:-2, :]) / 2

grad_mag = np.sqrt(Dx**2 + Dy**2)

plt.figure(figsize=(7, 4))
plt.imshow(grad_mag, cmap="magma")
plt.colorbar(label="gradient magnitude")
plt.title("Finite-difference gradient magnitude")
plt.axis("off")
plt.show()

The large values occur near sharp transitions and high-contrast structures.

4.19 3.18 Python Example: Differentiable Resizing in TensorFlow

Modern AI vision pipelines often resize images before feeding them to a model. The operation looks simple, but it is an interpolation and sampling operation.

import tensorflow as tf

# A batch of RGB images with shape (batch, height, width, channels).
images = tf.random.uniform(shape=(8, 256, 256, 3))

# Resize to the input size expected by many image models.
resized = tf.image.resize(
    images,
    size=(224, 224),
    method="bilinear",
    antialias=True
)

print(resized.shape)

The argument antialias=True matters when the output resolution is smaller than the input resolution. From the viewpoint of this chapter, it tells TensorFlow to include a low-pass filtering step during downsampling.

4.20 3.19 Worked Mathematical Example: Bilinear Weights

Suppose the four pixel values around a point are

\[ I[i,j]=0.1, \qquad I[i,j+1]=0.5, \qquad I[i+1,j]=0.4, \qquad I[i+1,j+1]=0.9. \]

Let $(x,y)$ be such that

\[ a = x-j = 0.25, \qquad b = y-i = 0.60. \]

The bilinear weights are

\[ \begin{aligned} w_{00} &= (1-a)(1-b) = 0.75 \cdot 0.40 = 0.30, \\ w_{10} &= a(1-b) = 0.25 \cdot 0.40 = 0.10, \\ w_{01} &= (1-a)b = 0.75 \cdot 0.60 = 0.45, \\ w_{11} &= ab = 0.25 \cdot 0.60 = 0.15. \end{aligned} \]

They sum to one:

\[ 0.30+0.10+0.45+0.15=1. \]

Thus

\[ \begin{aligned} R_{\text{bilin}}(I_h)(x,y) &=0.30(0.1)+0.10(0.5)+0.45(0.4)+0.15(0.9)\\ &=0.395. \end{aligned} \]

This calculation is simple, but it is exactly the local operation used inside many resizing, warping, grid-sampling, and differentiable augmentation routines.

4.21 3.20 Worked Mathematical Example: Two Frequencies with the Same Samples

Consider

\[ f_1(x)=\sin(2\pi(0.2)x) \]

and

\[ f_2(x)=\sin(2\pi(1.2)x). \]

At integer sample locations $x=n$,

\[ f_2(n)=\sin(2\pi(1.2)n)=\sin(2\pi(0.2)n+2\pi n)=\sin(2\pi(0.2)n)=f_1(n). \]

Therefore $f_1$ and $f_2$ produce the same samples even though they are different continuous functions. This is aliasing.

In image processing, this means two different continuous scenes or textures can produce the same sampled image array. Any reconstruction algorithm must therefore rely on assumptions.

4.22 3.21 Summary

This chapter introduced a fundamental bridge:

\[ I: \Omega \to \mathbb{R}^C \quad \xrightarrow{\text{sampling}} \quad I_h: \Lambda_{H,W} \to \mathbb{R}^C \quad \xrightarrow{\text{reconstruction}} \quad R(I_h): \Omega \to \mathbb{R}^C. \]

The main lessons are:

A continuous image and a sampled image are different mathematical objects.
Sampling loses information unless the image is sufficiently band-limited relative to the grid.
Interpolation is a modeling choice, not a purely technical detail.
Finite differences approximate derivatives and lead to edge, motion, and variational methods.
Geometric image transformations are naturally written using function composition.
Modern AI vision pipelines still depend on classical sampling, interpolation, and anti-aliasing principles.

4.23 Exercises

4.23.1 Conceptual Exercises

Explain the difference between an image as a matrix and an image as a function.
Why does rotating a digital image require interpolation?
Give an example where nearest-neighbor interpolation might be preferable to bilinear interpolation.
Explain why downsampling without smoothing can create artifacts.
What is the difference between a point-sample camera model and a pixel-average camera model?

4.23.2 Mathematical Exercises

Let $I(x,y)=x^2+y^2$. Sample it on the grid $x_j=jh$, $y_i=ih$. Compute $D_x^+I[i,j]$ and compare it with $I_x(x_j,y_i)$.
Prove that the bilinear interpolation weights are nonnegative and sum to one when $0\le a,b\le 1$.
Show that $\sin(2\pi\omega n)$ and $\sin(2\pi(\omega+k)n)$ have the same integer samples for every integer $k$.
Let $T(x,y)=(2x,2y)$. Write the transformed image $I_T$ using the pullback convention.
Suppose $Y[i,j]=I[i,j]+\varepsilon[i,j]$ with independent $N(0,\sigma^2)$ noise. Derive the negative log-likelihood up to an additive constant.

4.23.3 Python Exercises

Implement nearest-neighbor interpolation and compare it visually with bilinear interpolation.
Create a checkerboard image and downsample it by factors of 2, 4, and 8. Observe aliasing.
Implement central differences and compute gradient magnitude for a synthetic image.
Write a function that rotates an image around its center using backward warping and bilinear interpolation.
Compare resizing with and without anti-aliasing using TensorFlow or another image library.

4.23.4 AI-Assisted Study Exercises

Ask an AI tool to explain bilinear interpolation. Then revise the prompt to force it to state the coordinate convention and the four weights.
Ask an AI tool to generate NumPy code for rotating an image. Check whether it uses forward warping or backward warping.
Ask an AI tool why aliasing matters for CNN preprocessing. Evaluate whether the answer mentions low-pass filtering before downsampling.
Ask an AI tool to prove the aliasing identity for sinusoids. Verify every algebraic step yourself.
Give an AI tool a faulty bilinear interpolation function and ask it to find indexing errors. Then test the code on a $2\times2$ image where the answer can be computed by hand.

4.24 Notes and Further Reading

The function-based view of images is standard in classical image processing and computer vision; see Gonzalez and Woods (2018) and Szeliski (2022). The statistical view of images as observations generated from latent scenes is developed in model-based computer vision references such as Prince (2012). Later chapters will connect this foundation to Fourier analysis, convolution, edge detection, geometric warping, neural networks, and diffusion models.

Gonzalez, Rafael C., and Richard E. Woods. 2018. Digital Image Processing. 4th ed. Pearson.

Prince, Simon J. D. 2012. Computer Vision: Models, Learning, and Inference. Cambridge University Press. https://www.computervisionmodels.com/.

Szeliski, Richard. 2022. Computer Vision: Algorithms and Applications. 2nd ed. Springer. https://szeliski.org/Book/.

# Images as Functions on Domains ::: {.hero-box} A digital image is a finite array, but many of the most important ideas in computer vision become clean only after we view the array as samples of an underlying function. This chapter develops the bridge between continuous images, discrete samples, interpolation, finite differences, aliasing, random fields, and differentiable image transformations. The goal is to give applied mathematics and statistics students a precise language for moving between calculus, probability, numerical analysis, and modern AI vision models. ::: ## Learning Goals After reading this chapter, students should be able to: 1. model a grayscale or color image as a function $I: \Omega \to \mathbb{R}^C$ on a spatial domain; 2. distinguish the continuous image, the sampled image, and the reconstructed/interpolated image; 3. write sampling as an operator and understand discretization error; 4. explain nearest-neighbor, bilinear, and smoother interpolation from a mathematical viewpoint; 5. compute finite-difference approximations of image gradients; 6. explain aliasing using sinusoidal examples and the Nyquist idea; 7. describe geometric transformations as function composition and pullback operations; 8. connect classical image models with differentiable preprocessing and AI vision pipelines. ## 3.1 Why Images-as-Functions Matters Chapter 2 introduced pixels, coordinates, and color spaces. That chapter emphasized that a stored image is usually an array $$ I \in \mathbb{R}^{H \times W \times C}. $$ This array notation is necessary, but it is not always sufficient. When we resize an image, rotate it, blur it, compute a gradient, warp it under a camera model, or feed it to a neural network with differentiable data augmentation, we are implicitly treating the image as a function on a spatial domain. A clean mathematical description is $$ I: \Omega \subset \mathbb{R}^2 \longrightarrow \mathbb{R}^C, $$ where $\Omega$ is an image domain and $C$ is the number of channels. For a grayscale image, $C=1$. For an RGB image, $C=3$. For a feature map inside a convolutional neural network, $C$ may be 64, 256, 768, or larger. ::: {.math-spine} The central idea of this chapter is the distinction $$ \text{continuous image} \quad \longrightarrow \quad \text{sampled image} \quad \longrightarrow \quad \text{reconstructed image}. $$ A digital image is not just a function and not just an array. It is usually a sampled and quantized representation of an underlying signal, together with assumptions about how to reconstruct or interpolate it. ::: This viewpoint is useful across the whole book: - **Fourier analysis** studies the frequency content of image functions. - **Filtering** applies linear operators to image functions. - **Edges** are locations where derivatives or finite differences are large. - **Camera geometry** maps scene functions to image functions. - **Segmentation** partitions the image domain $\Omega$. - **Deep learning** repeatedly transforms sampled functions into new feature functions. - **Diffusion models** learn probability distributions over images as high-dimensional random functions. The mathematical language of functions also prevents many practical errors. For example, a rotation by $13^\circ$ does not map integer pixel centers to integer pixel centers. A model must decide how to assign intensity values at non-integer coordinates. That decision is interpolation, and it has mathematical consequences. ## 3.2 Continuous Images and Spatial Domains Let $\Omega \subset \mathbb{R}^2$ be a rectangular domain, $$ \Omega = [0,W] \times [0,H]. $$ A grayscale image is a scalar-valued function $$ I: \Omega \to \mathbb{R}. $$ A color image is a vector-valued function $$ I: \Omega \to \mathbb{R}^3, \qquad I(x,y) = \begin{bmatrix} R(x,y) \\ G(x,y) \\ B(x,y) \end{bmatrix}. $$ More generally, $$ I(x,y) \in \mathbb{R}^C. $$ Here $C$ is the number of channels. In a raw camera image, the channels may represent physical sensor responses. In a deep network, the channels may represent learned feature responses. ### Function spaces To do analysis, we often place images in function spaces. Common choices include: $$ L^2(\Omega) = \left\{ I: \Omega \to \mathbb{R} : \int_{\Omega} |I(x,y)|^2\,dx\,dy < \infty \right\}, $$ and, for differentiable images, $$ C^1(\Omega) = \left\{ I: \Omega \to \mathbb{R} : I_x \text{ and } I_y \text{ exist and are continuous} \right\}. $$ In practice, real images may contain discontinuities at object boundaries. Therefore $C^1(\Omega)$ is often too strong. A more realistic mathematical model says that images are smooth in regions but may have jumps along curves. ::: {.cv-card} **Piecewise-smooth image model.** A useful idealization is that an image is smooth inside objects and has discontinuities across object boundaries. Edges are then not mysterious artifacts; they are the visible traces of discontinuities or rapid changes in the image function. ::: This model helps explain why gradients and edges are central in computer vision. If $I$ is differentiable at $(x,y)$, then its first-order local approximation is $$ I(x+\Delta x, y+\Delta y) \approx I(x,y) + I_x(x,y)\Delta x + I_y(x,y)\Delta y. $$ The vector $$ \nabla I(x,y) = \begin{bmatrix} I_x(x,y) \\ I_y(x,y) \end{bmatrix} $$ points in the direction of greatest local increase of intensity. Its magnitude $\|\nabla I(x,y)\|$ measures local contrast. ## 3.3 Sampling Operators A camera does not store the whole function $I$. It samples it on a finite grid. Let $$ \Lambda_{H,W} = \{0,1,\ldots,H-1\} \times \{0,1,\ldots,W-1\} $$ be the set of array indices. Under the pixel-center convention, the sample location for row $i$ and column $j$ is $$ (x_j,y_i)=\left(j+\frac12,\ i+\frac12\right). $$ The sampling operator $S$ maps a continuous image to a discrete image: $$ S(I)[i,j] = I\left(j+\frac12, i+\frac12\right). $$ If the image has channels, then $$ S(I)[i,j,c] = I_c\left(j+\frac12, i+\frac12\right). $$ This formula is simple, but it hides an important assumption: the pixel value equals the value of the image function at the pixel center. A physically more realistic camera model may average light over a pixel cell: $$ S_{\text{avg}}(I)[i,j] = \frac{1}{|P_{ij}|} \int_{P_{ij}} I(x,y)\,dx\,dy, $$ where $$ P_{ij} = [j,j+1] \times [i,i+1]. $$ The point-sample model and the cell-average model lead to different behavior, especially when the image has high-frequency texture. ::: {.warning-note} Sampling is not a harmless bookkeeping step. If a continuous image changes too rapidly relative to the pixel grid, the sampled image may contain misleading low-frequency patterns. This is aliasing. ::: <div id="cv-ch03-sampling-demo" class="plotly-demo"></div> The interactive figure above shows a continuous image surface sampled on a finite grid. Increasing the spatial frequency of the underlying function would require a finer sampling grid to preserve the visible structure. ## 3.4 Discrete Images as Functions on a Grid Once sampled, an image is a function on a finite set: $$ I_h: \Lambda_{H,W} \to \mathbb{R}^C. $$ Here $h$ represents grid spacing. In many image arrays, the grid spacing is normalized to $h=1$, but for numerical analysis it is useful to write $$ x_j = x_0 + jh_x, \qquad y_i = y_0 + ih_y. $$ Then the discrete image is $$ I_h[i,j] = I(x_j,y_i). $$ The continuous function $I$ and the discrete function $I_h$ are mathematically different objects. The continuous function has values at every point in $\Omega$. The discrete image only has values at finitely many grid locations. A common source of confusion is to write $I(x,y)$ for both. In this book we will often use: $$ I(x,y) \quad \text{for a continuous image}, $$ and $$ I[i,j] \quad \text{for a sampled image array}. $$ When we move between the two notations, we are making a modeling choice. ## 3.5 Reconstruction and Interpolation Suppose we know only the discrete samples $I[i,j]$. What value should the image have at a non-grid point $(x,y)$? This is a reconstruction problem. A reconstruction operator $R$ maps a discrete image back to a continuous or piecewise-continuous image: $$ R(I_h)(x,y) \approx I(x,y). $$ Different reconstruction rules produce different images. ### Nearest-neighbor interpolation Nearest-neighbor interpolation chooses the value of the closest pixel center: $$ R_{\text{nn}}(I_h)(x,y) = I[i^*,j^*], $$ where $$ (i^*,j^*) = \arg\min_{(i,j) \in \Lambda_{H,W}} \left\|(x,y)-(x_j,y_i)\right\|_2. $$ This method is simple and preserves original pixel values, but it creates blocky artifacts. ### Bilinear interpolation Bilinear interpolation uses the four neighboring pixel values. Let $$ x \in [j,j+1], \qquad y \in [i,i+1], $$ and define $$ a = x-j, \qquad b = y-i, \qquad 0 \le a,b \le 1. $$ For a scalar image, the bilinear interpolant is $$ \begin{aligned} R_{\text{bilin}}(I_h)(x,y) =&(1-a)(1-b)I[i,j] + a(1-b)I[i,j+1] \\ &+ (1-a)b I[i+1,j] + ab I[i+1,j+1]. \end{aligned} $$ This formula is a convex combination of the four corner values. The weights are nonnegative and add to one: $$ (1-a)(1-b) + a(1-b) + (1-a)b + ab = 1. $$ Therefore, if the four corner values lie in $[0,1]$, the interpolated value also lies in $[0,1]$. <div id="cv-ch03-bilinear-demo" class="plotly-demo"></div> ### Interpolation as basis expansion Interpolation can also be written as a basis expansion: $$ R(I_h)(x,y) = \sum_{i,j} I[i,j] \varphi_{ij}(x,y), $$ where $\varphi_{ij}$ is a basis or kernel centered near pixel $(i,j)$. For nearest-neighbor interpolation, the basis functions are piecewise constant indicator functions. For bilinear interpolation, the basis functions are tent functions. For smoother interpolation, one may use cubic kernels or sinc-based reconstructions. This perspective is important because many image transformations in deep learning are differentiable only through the chosen interpolation formula. Bilinear interpolation is popular because it is simple, local, and piecewise differentiable with respect to both image values and coordinates. ## 3.6 Discretization Error Let $I$ be a continuous image and $I_h=S(I)$ its sampled version. After reconstruction, we obtain $R(I_h)$. A natural error is $$ \|I - R(S(I))\|. $$ The size of this error depends on: 1. the smoothness of $I$; 2. the grid spacing $h$; 3. the reconstruction method $R$; 4. the frequency content of $I$; 5. whether prefiltering was used before sampling. For a smooth one-dimensional function $f$, linear interpolation has local error of order $h^2$ under suitable assumptions: $$ |f(x) - R_{\text{linear}}(S(f))(x)| \le C h^2 \max_{\xi} |f''(\xi)|. $$ The two-dimensional bilinear case has analogous behavior when the second derivatives are controlled. The important qualitative lesson is: ::: {.math-spine} Smoother images can be sampled and reconstructed more accurately than highly oscillatory images. High-frequency content requires finer sampling or appropriate low-pass filtering. ::: This principle appears again in Fourier analysis, scale-space theory, anti-aliased image resizing, and neural network preprocessing. ## 3.7 Aliasing and the Nyquist Idea Aliasing occurs when different continuous signals produce the same samples on a grid. Consider the one-dimensional sinusoid $$ f(x) = \sin(2\pi \omega x), $$ sampled at integer points $x=n$. The samples are $$ f[n] = \sin(2\pi \omega n). $$ If $\omega$ is replaced by $\omega + k$, where $k$ is an integer, then $$ \sin(2\pi(\omega+k)n) = \sin(2\pi \omega n + 2\pi kn) = \sin(2\pi \omega n). $$ Thus many different continuous frequencies have exactly the same discrete samples. The Nyquist idea says that, for uniform sampling with spacing $h$, frequencies above approximately $$ \frac{1}{2h} $$ cannot be represented unambiguously without additional assumptions. <div id="cv-ch03-aliasing-demo" class="plotly-demo"></div> In images, aliasing appears as jagged edges, moire patterns, false textures, and strange artifacts after downsampling. A good image-resizing pipeline usually applies a low-pass filter before reducing resolution. ::: {.cv-card} **Anti-aliasing principle.** Before sampling more coarsely, remove frequencies that the coarser grid cannot represent. In practice, this means blur before downsampling. ::: This principle is not only classical signal processing. It also affects modern AI systems. A classifier may change its prediction under small shifts or resizes partly because the preprocessing pipeline introduces aliasing. ## 3.8 Finite Differences and Image Derivatives For a differentiable image $I(x,y)$, the gradient is $$ \nabla I(x,y) = \begin{bmatrix} I_x(x,y) \\ I_y(x,y) \end{bmatrix}. $$ For a sampled image $I[i,j]$, derivatives must be approximated. Common finite-difference formulas include: ### Forward differences $$ D_x^+ I[i,j] = I[i,j+1] - I[i,j], $$ $$ D_y^+ I[i,j] = I[i+1,j] - I[i,j]. $$ ### Backward differences $$ D_x^- I[i,j] = I[i,j] - I[i,j-1], $$ $$ D_y^- I[i,j] = I[i,j] - I[i-1,j]. $$ ### Central differences $$ D_x^0 I[i,j] = \frac{I[i,j+1]-I[i,j-1]}{2}, $$ $$ D_y^0 I[i,j] = \frac{I[i+1,j]-I[i-1,j]}{2}. $$ With grid spacings $h_x$ and $h_y$, the denominators become $h_x$, $h_y$, $2h_x$, or $2h_y$ accordingly. The discrete gradient magnitude is often approximated by $$ \|\nabla_h I[i,j]\|_2 = \sqrt{(D_x I[i,j])^2 + (D_y I[i,j])^2}. $$ <div id="cv-ch03-gradient-demo" class="plotly-demo"></div> This discrete derivative viewpoint is foundational for edge detection, optical flow, variational segmentation, total variation denoising, and numerical PDE methods in imaging. ## 3.9 Difference Operators as Matrices Finite differences can also be written as linear algebra. Suppose a one-dimensional signal is $$ f = \begin{bmatrix} f_0 & f_1 & \cdots & f_{n-1} \end{bmatrix}^T. $$ A forward-difference operator can be represented by the matrix $$ D = \begin{bmatrix} -1 & 1 & 0 & \cdots & 0 \\ 0 & -1 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 \end{bmatrix}. $$ Then $$ Df = \begin{bmatrix} f_1-f_0 \\ f_2-f_1 \\ \vdots \\ f_{n-1}-f_{n-2} \end{bmatrix}. $$ For two-dimensional images, difference operators act along rows and columns. If the image is vectorized into a vector $u \in \mathbb{R}^{HW}$, then horizontal and vertical difference operators may be written as matrices $$ D_x u, \qquad D_y u. $$ The discrete gradient operator is $$ D u = \begin{bmatrix} D_x u \\ D_y u \end{bmatrix}. $$ This notation leads naturally to optimization problems such as $$ \min_u \frac12\|u-y\|_2^2 + \lambda \|Du\|_1, $$ which is a simplified form of total-variation denoising. The first term keeps $u$ close to the observed image $y$, while the second term encourages piecewise-smooth structure. ## 3.10 Geometric Transformations as Function Composition Let $T: \Omega \to \Omega$ be a spatial transformation. For example: - translation; - rotation; - scaling; - affine transformation; - projective transformation; - nonlinear warp. If we want to transform an image, we need to define the transformed function. A common convention is the pullback: $$ I_T(x,y) = I(T^{-1}(x,y)). $$ This means that the output value at location $(x,y)$ is obtained by looking backward to the corresponding location in the input image. Why use $T^{-1}$ instead of $T$? Because we want to assign a value to every output pixel. For each output location, we ask: where did this point come from in the input image? ::: {.cv-card} **Backward warping.** In image processing and differentiable vision layers, warping is often implemented by looping over output coordinates, mapping each one backward to an input coordinate, and then interpolating the input image at that coordinate. ::: For example, if $T$ is a translation by $(a,b)$, $$ T(x,y)=(x+a,y+b), $$ then $$ I_T(x,y) = I(x-a,y-b). $$ If $(x-a,y-b)$ is not a pixel center, interpolation is required. This is the mathematical foundation behind image augmentation, spatial transformer networks, differentiable rendering, optical flow warping, and many registration algorithms. ## 3.11 Images as Random Functions For statistics students, it is useful to view images as random objects. A random image may be written as $$ X: \Omega \times \Xi \to \mathbb{R}^C, $$ where $\Xi$ is a probability space. For each outcome $\xi \in \Xi$, the function $$ X(\cdot,\cdot,\xi): \Omega \to \mathbb{R}^C $$ is one image. A simple observation model is $$ Y(x,y) = I(x,y) + \varepsilon(x,y), $$ where $\varepsilon$ is noise. In a discrete model, $$ Y[i,j] = I[i,j] + \varepsilon[i,j]. $$ If the noise is independent Gaussian noise, $$ \varepsilon[i,j] \sim N(0,\sigma^2), $$ then the negative log-likelihood is proportional to the squared error $$ \sum_{i,j} (Y[i,j]-I[i,j])^2. $$ This is why least-squares losses appear naturally in denoising, reconstruction, registration, and many image-to-image prediction tasks. More structured statistical models treat images as random fields with covariance $$ K\bigl((x,y),(x',y')\bigr) = \operatorname{Cov}\bigl(X(x,y), X(x',y')\bigr). $$ A stationary covariance depends only on displacement: $$ K\bigl((x,y),(x',y')\bigr) = k(x-x',y-y'). $$ This connects image modeling to Gaussian processes, Markov random fields, texture models, and modern generative models. ## 3.12 Function Viewpoint Inside Deep Networks A convolutional neural network does not only process the original image. It builds a sequence of feature maps: $$ F_0, F_1, F_2, \ldots, F_L. $$ The input might be $$ F_0: \Omega_0 \to \mathbb{R}^3, $$ while an intermediate layer might be $$ F_\ell: \Omega_\ell \to \mathbb{R}^{C_\ell}. $$ Each spatial location has a feature vector. Therefore a feature map is a vector-valued function on a discrete spatial domain. A convolutional layer is approximately a local operator: $$ F_{\ell+1}(p) = \sigma\left(\sum_{q \in N(p)} W(q-p)F_\ell(q)+b\right), $$ where: - $p$ is a spatial location; - $N(p)$ is a local neighborhood; - $W$ is a learned kernel; - $\sigma$ is a nonlinear activation. A vision transformer uses a different representation. It partitions the image into patches and converts each patch into a token vector. This changes the spatial function into a sequence: $$ I: \Omega \to \mathbb{R}^3 \quad \longrightarrow \quad (z_1,z_2,\ldots,z_N), \qquad z_i \in \mathbb{R}^d. $$ Even then, the construction begins with a sampled image function and a choice of patch geometry. ## 3.13 AI Component: Using Models to Reason About Image Functions Large language models and multimodal AI systems can help students reason about computer vision, but only if the prompts are mathematically specific. ### Prompt pattern 1: separate objects carefully A useful study prompt is: > I have a continuous image $I:\Omega\to\mathbb{R}$, a sampled image $I_h[i,j]$, and a reconstructed image $R(I_h)(x,y)$. Explain the difference among these three objects, and give one example where confusing them causes an error in computer vision. A weak answer will say only that one is continuous and one is discrete. A strong answer will mention sampling, interpolation, coordinate conventions, and possible aliasing. ### Prompt pattern 2: ask for formula verification > Derive the bilinear interpolation formula from the requirement that the interpolated value is linear in $x$ for fixed $y$ and linear in $y$ for fixed $x$. Check that the four weights sum to one. This prompt forces the AI tool to work through the mathematics instead of giving only implementation advice. ### Prompt pattern 3: connect theory to code > Given a NumPy image array of shape `(H, W)`, write a function for bilinear interpolation at non-integer coordinates `(x, y)`. Explain how the formula corresponds to the mathematical weights $(1-a)(1-b)$, $a(1-b)$, $(1-a)b$, and $ab$. This type of prompt is useful because it asks for both computation and mathematical interpretation. ### Prompt pattern 4: audit preprocessing > A vision model resizes all input images from $1024\times 1024$ to $224\times 224$. Explain why this is a sampling problem. What can go wrong if no low-pass filtering is used before downsampling? This prompt connects classical sampling theory with modern deep learning practice. ::: {.warning-note} When using AI for study, do not ask only “explain interpolation.” Ask the model to specify the domain, the grid, the coordinates, the interpolation weights, and the assumptions. Computer vision errors often hide in these details. ::: ## 3.14 Python Example: Sampling a Continuous Image The following code constructs a synthetic continuous image function and samples it on a grid. ```{python} #| label: ch03-sample-continuous-image #| eval: false import numpy as np import matplotlib.pyplot as plt def continuous_image(x, y): """A smooth synthetic image function on [0, 1] x [0, 1].""" low_freq = 0.55 + 0.25 * np.sin(2 * np.pi * x) * np.cos(2 * np.pi * y) blob = 0.45 * np.exp(-45 * ((x - 0.65)**2 + (y - 0.35)**2)) edge = 0.25 * (x + y > 1.05) return np.clip(low_freq + blob + edge, 0, 1) H, W = 80, 120 ys = (np.arange(H) + 0.5) / H xs = (np.arange(W) + 0.5) / W X, Y = np.meshgrid(xs, ys) I = continuous_image(X, Y) plt.figure(figsize=(7, 4)) plt.imshow(I, cmap="gray", origin="upper", extent=[0, 1, 1, 0]) plt.colorbar(label="intensity") plt.title("Sampled synthetic image") plt.xlabel("x") plt.ylabel("y") plt.show() ``` Notice the use of pixel centers: $$ x_j = \frac{j+1/2}{W}, \qquad y_i = \frac{i+1/2}{H}. $$ Changing this convention slightly may shift the image and affect geometric alignment. ## 3.15 Python Example: Bilinear Interpolation The next example implements bilinear interpolation at a non-integer coordinate. ```{python} #| label: ch03-bilinear-interpolation #| eval: false import numpy as np def bilinear_interpolate(image, x, y): """ Bilinear interpolation for a grayscale image. Parameters ---------- image : ndarray of shape (H, W) Input image. x, y : float Continuous coordinates in array coordinate convention: x is column coordinate, y is row coordinate. Returns ------- value : float Interpolated intensity. """ H, W = image.shape # Clamp to valid interpolation range. x = np.clip(x, 0, W - 1) y = np.clip(y, 0, H - 1) j0 = int(np.floor(x)) i0 = int(np.floor(y)) j1 = min(j0 + 1, W - 1) i1 = min(i0 + 1, H - 1) a = x - j0 b = y - i0 return ( (1 - a) * (1 - b) * image[i0, j0] + a * (1 - b) * image[i0, j1] + (1 - a) * b * image[i1, j0] + a * b * image[i1, j1] ) # Example small = np.array([ [0.0, 0.2, 0.4], [0.3, 0.6, 0.9], [0.5, 0.7, 1.0] ]) print(bilinear_interpolate(small, x=1.25, y=0.75)) ``` A useful debugging check is to verify that if $(x,y)$ lands exactly on a pixel coordinate, the function returns that pixel value. ## 3.16 Python Example: Aliasing Under Downsampling This example shows why high-frequency patterns need anti-aliasing before downsampling. ```{python} #| label: ch03-aliasing-python #| eval: false import numpy as np import matplotlib.pyplot as plt n_high = 1024 x_high = np.linspace(0, 1, n_high, endpoint=False) freq = 37 signal_high = np.sin(2 * np.pi * freq * x_high) # Downsample by taking every 16th point. factor = 16 x_low = x_high[::factor] signal_low = signal_high[::factor] plt.figure(figsize=(8, 4)) plt.plot(x_high, signal_high, alpha=0.5, label="continuous-looking high-resolution signal") plt.scatter(x_low, signal_low, label="samples after downsampling") plt.legend() plt.title("Aliasing from undersampling") plt.xlabel("x") plt.ylabel("signal") plt.show() ``` A better pipeline would apply a low-pass filter before keeping every 16th sample. ## 3.17 Python Example: Finite-Difference Gradients The following code computes a finite-difference gradient magnitude image. ```{python} #| label: ch03-finite-difference-gradient #| eval: false import numpy as np import matplotlib.pyplot as plt H, W = 120, 160 y = np.linspace(0, 1, H) x = np.linspace(0, 1, W) X, Y = np.meshgrid(x, y) I = 0.4 + 0.3 * X + 0.2 * np.exp(-60 * ((X - 0.65)**2 + (Y - 0.45)**2)) I += 0.35 * (X + Y > 1.15) I = np.clip(I, 0, 1) Dx = np.zeros_like(I) Dy = np.zeros_like(I) Dx[:, 1:-1] = (I[:, 2:] - I[:, :-2]) / 2 Dy[1:-1, :] = (I[2:, :] - I[:-2, :]) / 2 grad_mag = np.sqrt(Dx**2 + Dy**2) plt.figure(figsize=(7, 4)) plt.imshow(grad_mag, cmap="magma") plt.colorbar(label="gradient magnitude") plt.title("Finite-difference gradient magnitude") plt.axis("off") plt.show() ``` The large values occur near sharp transitions and high-contrast structures. ## 3.18 Python Example: Differentiable Resizing in TensorFlow Modern AI vision pipelines often resize images before feeding them to a model. The operation looks simple, but it is an interpolation and sampling operation. ```{python} #| label: ch03-tensorflow-resize #| eval: false import tensorflow as tf # A batch of RGB images with shape (batch, height, width, channels). images = tf.random.uniform(shape=(8, 256, 256, 3)) # Resize to the input size expected by many image models. resized = tf.image.resize( images, size=(224, 224), method="bilinear", antialias=True ) print(resized.shape) ``` The argument `antialias=True` matters when the output resolution is smaller than the input resolution. From the viewpoint of this chapter, it tells TensorFlow to include a low-pass filtering step during downsampling. ## 3.19 Worked Mathematical Example: Bilinear Weights Suppose the four pixel values around a point are $$ I[i,j]=0.1, \qquad I[i,j+1]=0.5, \qquad I[i+1,j]=0.4, \qquad I[i+1,j+1]=0.9. $$ Let $(x,y)$ be such that $$ a = x-j = 0.25, \qquad b = y-i = 0.60. $$ The bilinear weights are $$ \begin{aligned} w_{00} &= (1-a)(1-b) = 0.75 \cdot 0.40 = 0.30, \\ w_{10} &= a(1-b) = 0.25 \cdot 0.40 = 0.10, \\ w_{01} &= (1-a)b = 0.75 \cdot 0.60 = 0.45, \\ w_{11} &= ab = 0.25 \cdot 0.60 = 0.15. \end{aligned} $$ They sum to one: $$ 0.30+0.10+0.45+0.15=1. $$ Thus $$ \begin{aligned} R_{\text{bilin}}(I_h)(x,y) &=0.30(0.1)+0.10(0.5)+0.45(0.4)+0.15(0.9)\\ &=0.395. \end{aligned} $$ This calculation is simple, but it is exactly the local operation used inside many resizing, warping, grid-sampling, and differentiable augmentation routines. ## 3.20 Worked Mathematical Example: Two Frequencies with the Same Samples Consider $$ f_1(x)=\sin(2\pi(0.2)x) $$ and $$ f_2(x)=\sin(2\pi(1.2)x). $$ At integer sample locations $x=n$, $$ f_2(n)=\sin(2\pi(1.2)n)=\sin(2\pi(0.2)n+2\pi n)=\sin(2\pi(0.2)n)=f_1(n). $$ Therefore $f_1$ and $f_2$ produce the same samples even though they are different continuous functions. This is aliasing. In image processing, this means two different continuous scenes or textures can produce the same sampled image array. Any reconstruction algorithm must therefore rely on assumptions. ## 3.21 Summary This chapter introduced a fundamental bridge: $$ I: \Omega \to \mathbb{R}^C \quad \xrightarrow{\text{sampling}} \quad I_h: \Lambda_{H,W} \to \mathbb{R}^C \quad \xrightarrow{\text{reconstruction}} \quad R(I_h): \Omega \to \mathbb{R}^C. $$ The main lessons are: 1. A continuous image and a sampled image are different mathematical objects. 2. Sampling loses information unless the image is sufficiently band-limited relative to the grid. 3. Interpolation is a modeling choice, not a purely technical detail. 4. Finite differences approximate derivatives and lead to edge, motion, and variational methods. 5. Geometric image transformations are naturally written using function composition. 6. Modern AI vision pipelines still depend on classical sampling, interpolation, and anti-aliasing principles. ## Exercises ### Conceptual Exercises 1. Explain the difference between an image as a matrix and an image as a function. 2. Why does rotating a digital image require interpolation? 3. Give an example where nearest-neighbor interpolation might be preferable to bilinear interpolation. 4. Explain why downsampling without smoothing can create artifacts. 5. What is the difference between a point-sample camera model and a pixel-average camera model? ### Mathematical Exercises 1. Let $I(x,y)=x^2+y^2$. Sample it on the grid $x_j=jh$, $y_i=ih$. Compute $D_x^+I[i,j]$ and compare it with $I_x(x_j,y_i)$. 2. Prove that the bilinear interpolation weights are nonnegative and sum to one when $0\le a,b\le 1$. 3. Show that $\sin(2\pi\omega n)$ and $\sin(2\pi(\omega+k)n)$ have the same integer samples for every integer $k$. 4. Let $T(x,y)=(2x,2y)$. Write the transformed image $I_T$ using the pullback convention. 5. Suppose $Y[i,j]=I[i,j]+\varepsilon[i,j]$ with independent $N(0,\sigma^2)$ noise. Derive the negative log-likelihood up to an additive constant. ### Python Exercises 1. Implement nearest-neighbor interpolation and compare it visually with bilinear interpolation. 2. Create a checkerboard image and downsample it by factors of 2, 4, and 8. Observe aliasing. 3. Implement central differences and compute gradient magnitude for a synthetic image. 4. Write a function that rotates an image around its center using backward warping and bilinear interpolation. 5. Compare resizing with and without anti-aliasing using TensorFlow or another image library. ### AI-Assisted Study Exercises 1. Ask an AI tool to explain bilinear interpolation. Then revise the prompt to force it to state the coordinate convention and the four weights. 2. Ask an AI tool to generate NumPy code for rotating an image. Check whether it uses forward warping or backward warping. 3. Ask an AI tool why aliasing matters for CNN preprocessing. Evaluate whether the answer mentions low-pass filtering before downsampling. 4. Ask an AI tool to prove the aliasing identity for sinusoids. Verify every algebraic step yourself. 5. Give an AI tool a faulty bilinear interpolation function and ask it to find indexing errors. Then test the code on a $2\times2$ image where the answer can be computed by hand. ## Notes and Further Reading The function-based view of images is standard in classical image processing and computer vision; see @gonzalez2018digital and @szeliski2022computer. The statistical view of images as observations generated from latent scenes is developed in model-based computer vision references such as @prince2012computer. Later chapters will connect this foundation to Fourier analysis, convolution, edge detection, geometric warping, neural networks, and diffusion models. <script> function cvCh03MakeRange(start, stop, n) { const arr = []; if (n === 1) return [start]; for (let i = 0; i < n; i++) arr.push(start + (stop - start) * i / (n - 1)); return arr; } function cvCh03ContinuousImage(x, y) { const low = 0.55 + 0.22 * Math.sin(2 * Math.PI * x) * Math.cos(2 * Math.PI * y); const blob = 0.38 * Math.exp(-45 * ((x - 0.65) ** 2 + (y - 0.35) ** 2)); const edge = (x + y > 1.08) ? 0.25 : 0.0; return Math.max(0, Math.min(1, low + blob + edge)); } function cvCh03RenderSamplingDemo() { const el = document.getElementById('cv-ch03-sampling-demo'); if (!el || typeof Plotly === 'undefined') return; const xFine = cvCh03MakeRange(0, 1, 70); const yFine = cvCh03MakeRange(0, 1, 50); const zFine = yFine.map(y => xFine.map(x => cvCh03ContinuousImage(x, y))); const xs = []; const ys = []; const zs = []; const H = 10; const W = 14; for (let i = 0; i < H; i++) { const y = (i + 0.5) / H; for (let j = 0; j < W; j++) { const x = (j + 0.5) / W; xs.push(x); ys.push(y); zs.push(cvCh03ContinuousImage(x, y)); } } const surface = { x: xFine, y: yFine, z: zFine, type: 'surface', opacity: 0.82, colorscale: 'Viridis', showscale: false, name: 'continuous image function', hovertemplate: 'x=%{x:.2f}<br>y=%{y:.2f}<br>I=%{z:.3f}<extra></extra>' }; const samples = { x: xs, y: ys, z: zs, type: 'scatter3d', mode: 'markers', marker: {size: 4, color: zs, colorscale: 'Viridis', line: {width: 1}}, name: 'samples', hovertemplate: 'sample x=%{x:.2f}<br>sample y=%{y:.2f}<br>I=%{z:.3f}<extra></extra>' }; Plotly.newPlot(el, [surface, samples], { title: 'A continuous image function sampled on a finite grid', margin: {t: 55, r: 20, b: 35, l: 20}, scene: { xaxis: {title: 'x'}, yaxis: {title: 'y'}, zaxis: {title: 'intensity', range: [0, 1.2]}, camera: {eye: {x: 1.55, y: -1.75, z: 1.25}} } }, {responsive: true, displaylogo: false}); } function cvCh03RenderBilinearDemo() { const el = document.getElementById('cv-ch03-bilinear-demo'); if (!el || typeof Plotly === 'undefined') return; const I00 = 0.1, I10 = 0.8, I01 = 0.45, I11 = 0.95; const xs = cvCh03MakeRange(0, 1, 35); const ys = cvCh03MakeRange(0, 1, 35); const z = ys.map(b => xs.map(a => (1 - a) * (1 - b) * I00 + a * (1 - b) * I10 + (1 - a) * b * I01 + a * b * I11 )); const corners = { x: [0, 1, 0, 1], y: [0, 0, 1, 1], z: [I00, I10, I01, I11], type: 'scatter3d', mode: 'markers+text', marker: {size: 6}, text: ['I[i,j]', 'I[i,j+1]', 'I[i+1,j]', 'I[i+1,j+1]'], textposition: 'top center', name: 'corner samples', hovertemplate: '%{text}<br>value=%{z:.2f}<extra></extra>' }; Plotly.newPlot(el, [ { x: xs, y: ys, z: z, type: 'surface', colorscale: 'Viridis', opacity: 0.9, name: 'bilinear interpolant', hovertemplate: 'a=%{x:.2f}<br>b=%{y:.2f}<br>value=%{z:.3f}<extra></extra>' }, corners ], { title: 'Bilinear interpolation over one pixel cell', margin: {t: 55, r: 20, b: 35, l: 20}, scene: { xaxis: {title: 'a = x - j'}, yaxis: {title: 'b = y - i'}, zaxis: {title: 'interpolated intensity', range: [0, 1]}, camera: {eye: {x: 1.45, y: -1.65, z: 1.1}} } }, {responsive: true, displaylogo: false}); } function cvCh03RenderAliasingDemo() { const el = document.getElementById('cv-ch03-aliasing-demo'); if (!el || typeof Plotly === 'undefined') return; const xFine = cvCh03MakeRange(0, 10, 1000); const sampleX = []; for (let n = 0; n <= 10; n++) sampleX.push(n); const fLow = 0.2; const fHigh = 1.2; const lowY = xFine.map(x => Math.sin(2 * Math.PI * fLow * x)); const highY = xFine.map(x => Math.sin(2 * Math.PI * fHigh * x)); const sampleY = sampleX.map(x => Math.sin(2 * Math.PI * fLow * x)); Plotly.newPlot(el, [ { x: xFine, y: lowY, type: 'scatter', mode: 'lines', name: 'frequency 0.2', hovertemplate: 'x=%{x:.2f}<br>f=%{y:.3f}<extra></extra>' }, { x: xFine, y: highY, type: 'scatter', mode: 'lines', name: 'frequency 1.2', line: {dash: 'dash'}, hovertemplate: 'x=%{x:.2f}<br>f=%{y:.3f}<extra></extra>' }, { x: sampleX, y: sampleY, type: 'scatter', mode: 'markers', marker: {size: 9}, name: 'integer samples', hovertemplate: 'n=%{x}<br>sample=%{y:.3f}<extra></extra>' } ], { title: 'Aliasing: two different continuous sinusoids with the same integer samples', margin: {t: 55, r: 25, b: 55, l: 65}, xaxis: {title: 'x'}, yaxis: {title: 'signal value', range: [-1.15, 1.15]}, legend: {orientation: 'h'} }, {responsive: true, displaylogo: false}); } function cvCh03RenderGradientDemo() { const el = document.getElementById('cv-ch03-gradient-demo'); if (!el || typeof Plotly === 'undefined') return; const H = 45; const W = 60; const z = []; for (let i = 0; i < H; i++) { const row = []; const y = i / (H - 1); for (let j = 0; j < W; j++) { const x = j / (W - 1); const ramp = 0.35 + 0.25 * x; const blob = 0.35 * Math.exp(-55 * ((x - 0.65) ** 2 + (y - 0.45) ** 2)); const jump = (x + y > 1.15) ? 0.25 : 0.0; row.push(Math.max(0, Math.min(1, ramp + blob + jump))); } z.push(row); } const g = []; for (let i = 0; i < H; i++) { const row = []; for (let j = 0; j < W; j++) { const jm = Math.max(j - 1, 0), jp = Math.min(j + 1, W - 1); const im = Math.max(i - 1, 0), ip = Math.min(i + 1, H - 1); const dx = (z[i][jp] - z[i][jm]) / 2; const dy = (z[ip][j] - z[im][j]) / 2; row.push(Math.sqrt(dx * dx + dy * dy)); } g.push(row); } Plotly.newPlot(el, [ { z: z, type: 'heatmap', colorscale: 'Gray', showscale: false, name: 'image', xaxis: 'x', yaxis: 'y', hovertemplate: 'col=%{x}<br>row=%{y}<br>I=%{z:.3f}<extra></extra>' }, { z: g, type: 'heatmap', colorscale: 'Viridis', name: 'gradient magnitude', xaxis: 'x2', yaxis: 'y2', hovertemplate: 'col=%{x}<br>row=%{y}<br>|grad|=%{z:.3f}<extra></extra>' } ], { title: 'Finite differences: image intensity and gradient magnitude', margin: {t: 55, r: 30, b: 45, l: 45}, grid: {rows: 1, columns: 2, pattern: 'independent'}, xaxis: {title: 'column'}, yaxis: {title: 'row', autorange: 'reversed'}, xaxis2: {title: 'column'}, yaxis2: {title: 'row', autorange: 'reversed'} }, {responsive: true, displaylogo: false}); } document.addEventListener('DOMContentLoaded', () => { cvCh03RenderSamplingDemo(); cvCh03RenderBilinearDemo(); cvCh03RenderAliasingDemo(); cvCh03RenderGradientDemo(); }); </script>