Integral and Discrete Transforms in Image Processing
Sampling & Resampling
David Svoboda
email: svoboda@fi.muni.cz
Centre for Biomedical Image Analysis
Faculty of Informatics, Masaryk University, Brno, CZ
October 14, 2022
David Svoboda (CBIA@FI) Image Transforms autumn 2022 1 / 50
Outline
1 Sampling in 1D
Comb function
Nyquist-Shannon theorem
Aliasing
Reconstruction/Interpolation
2 Sampling in 2D
Comb function
Aliasing
Reconstruction/Interpolation
3 Resampling in 1D
Design of an ideal resampling filter
4 Resampling in 2D
Elliptical Weighted Average (EWA)
David Svoboda (CBIA@FI) Image Transforms autumn 2022 2 / 50
1 Sampling in 1D
Comb function
Nyquist-Shannon theorem
Aliasing
Reconstruction/Interpolation
2 Sampling in 2D
Comb function
Aliasing
Reconstruction/Interpolation
3 Resampling in 1D
Design of an ideal resampling filter
4 Resampling in 2D
Elliptical Weighted Average (EWA)
David Svoboda (CBIA@FI) Image Transforms autumn 2022 3 / 50
Comb function
In 1D space
III(x) =
∞
n=−∞
δ(x − n)
x
Notice: “III” is pronouced as shah (Cyrilic character).
David Svoboda (CBIA@FI) Image Transforms autumn 2022 4 / 50
Comb function
Some properties
III(−x) = III(x)
III(x + n) = III(x)
III(x − 1
2) = III(x + 1
2)
III(x) = 0 x ̸= n
III(ax) = 1
|a| δ(x − n
a )
III(x
τ ) ⊃ τIII(τω)
David Svoboda (CBIA@FI) Image Transforms autumn 2022 5 / 50
Sampling
Sampling in 1D:
a process of converting a continuous signal into a discrete sequence.
a multiplication with comb function:
III(x)f (x) =
∞
n=−∞
f (n)δ(x − n)
Question: What happens in frequency (Fourier) domain?
David Svoboda (CBIA@FI) Image Transforms autumn 2022 6 / 50
Sampling
Image domain versus Fourier domain
Image/Time domain:
multiplication of the function f and III
sampling
Fourier domain:
convolution of the function FT(f ) and FT(III)
convolution with Dirac impulses → replication of FT(f )
scaling property is also valid for III function
David Svoboda (CBIA@FI) Image Transforms autumn 2022 7 / 50
Sampling
x
f(x)
u
x u
III(x/t) t III(tu)
F(u)
1/t
max frequency
Notice: The comb function density must be high enough to guarantee
proper sampling
David Svoboda (CBIA@FI) Image Transforms autumn 2022 8 / 50
Sampling
III(x/t)f(x)
x
x u
u
u
t III(tu)*F(u)
x
David Svoboda (CBIA@FI) Image Transforms autumn 2022 9 / 50
Sampling
Nyquist-Shannon theorem
Exact reconstruction of a continuous signal from its samples is possible if
the signal is bandlimited and the sampling frequency is greater than twice
the signal maximal frequency
Harry Nyquist (1889 – 1976) & Claude Elwood Shannon (1916 – 2001)
Question: How to use N-S theorem, if the original signal is unlimited?
David Svoboda (CBIA@FI) Image Transforms autumn 2022 10 / 50
Sampling
Common problems – aliasing
The cause of aliasing:
when Nyquist-Shannon condition is broken, i.e.
sampling frequency is not high enough or
(time alias – wagon wheel effect)
the signal in not bandlimited
(PC games – far horizon)
David Svoboda (CBIA@FI) Image Transforms autumn 2022 11 / 50
Sampling
Common problems – aliasing
How to eliminate aliasing?
sampling at higher frequency
does it help if the signal is not band limited?
expensive for memory and time
OR
prefiltering
before sampling the input signal is “prefiltered” by lowpass filter
prefilter sampling
continuouscontinuous discrete
bandlimited
David Svoboda (CBIA@FI) Image Transforms autumn 2022 12 / 50
Sampling
Common problems – aliasing
Some lowpass filters
Gaussian filter
fσ(x) =
1
σ
√
2π
e− x2
2σ2
Sinc filter
f (x) =
sin(x)
x
B-spline filter
b1(x) =
1 |x| ≤ 1/2
0 |x| > 1/2
bn(x) = b1(x) ∗ b1(x) ∗ · · · ∗ b1(x) n-times
David Svoboda (CBIA@FI) Image Transforms autumn 2022 13 / 50
Reconstruction/Interpolation
Inverse process to sampling
The purpose: reconstruction of the original continuous signal from the
sampled sequence.
Reconstruction ≡ convolution with a low-pass filter.
Common reconstruction filters:
box (nearest neighbour)
tent (linear interpolation)
cubic B-spline (cubic polynomial interpolation)
Gaussian
sinc function
Lanczos (windowed sinc function)
Notice: The unit area under the curve.
David Svoboda (CBIA@FI) Image Transforms autumn 2022 14 / 50
Reconstruction/Interpolation
Lanczos filter
−4 −3 −2 −1 0 1 2 3 4
−0.2
0
0.2
0.4
0.6
0.8
1
1.2
Lanczos3
lanczos3(x) =
sin πx
πx
sin π x
3
π x
3
|x| < 3
0 |x| ≥ 3
David Svoboda (CBIA@FI) Image Transforms autumn 2022 15 / 50
Reconstruction/Interpolation
Example
Box filter
↔
David Svoboda (CBIA@FI) Image Transforms autumn 2022 16 / 50
Reconstruction/Interpolation
Example
Tent filter
↔
David Svoboda (CBIA@FI) Image Transforms autumn 2022 17 / 50
Reconstruction/Interpolation
Example
Cubic B-spline filter
↔
David Svoboda (CBIA@FI) Image Transforms autumn 2022 18 / 50
1 Sampling in 1D
Comb function
Nyquist-Shannon theorem
Aliasing
Reconstruction/Interpolation
2 Sampling in 2D
Comb function
Aliasing
Reconstruction/Interpolation
3 Resampling in 1D
Design of an ideal resampling filter
4 Resampling in 2D
Elliptical Weighted Average (EWA)
David Svoboda (CBIA@FI) Image Transforms autumn 2022 19 / 50
Comb function
In 2D space
2
III(x, y) =
∞
m=−∞
∞
n=−∞
δ(x − m, y − n)
Separability of delta function implies:
2
III(x, y) = III(x)III(y)
David Svoboda (CBIA@FI) Image Transforms autumn 2022 20 / 50
Sampling
Sampling in 2D:
a process of converting a continuous 2D function into a discrete grid
a multiplictation with comb function:
2
III(x, y)f (x, y) =
m n
f (m, n)δ(x − m, y − n)
Sampling in image domain
·
x
y
David Svoboda (CBIA@FI) Image Transforms autumn 2022 21 / 50
Sampling
Sampling in 2D:
a process of converting a continuous 2D function into a discrete grid
a multiplictation with comb function:
2
III(x, y)f (x, y) =
m n
f (m, n)δ(x − m, y − n)
Sampling in frequency domain
∗
ωx
ωy
David Svoboda (CBIA@FI) Image Transforms autumn 2022 21 / 50
Sampling
Common problems – aliasing
Insufficient sampling
reconstructed image
David Svoboda (CBIA@FI) Image Transforms autumn 2022 22 / 50
Sampling
Common problems – aliasing
Sufficient sampling
reconstructed image
David Svoboda (CBIA@FI) Image Transforms autumn 2022 23 / 50
Sampling
Common problems – aliasing
An example
David Svoboda (CBIA@FI) Image Transforms autumn 2022 24 / 50
Sampling
Common problems – aliasing
An example
David Svoboda (CBIA@FI) Image Transforms autumn 2022 25 / 50
Reconstruction/Interpolation
Task to solve
Problem of missing pixels when zooming into the digital image
3x magnification
David Svoboda (CBIA@FI) Image Transforms autumn 2022 26 / 50
Reconstruction/Interpolation
Solution
Nearest neighbour interpolation
Kernel3×3 =


1 1 1
1 1 1
1 1 1


0
0.5
-2
-1
1
0 211 0-12 -2
Note: Kernel3×3 = [1 1 1]T
[1 1 1]
David Svoboda (CBIA@FI) Image Transforms autumn 2022 27 / 50
Reconstruction/Interpolation
Solution
Completion of missing pixels (nearest neighbour)
3x magnification + nearest neighbour interpolation
David Svoboda (CBIA@FI) Image Transforms autumn 2022 28 / 50
Reconstruction/Interpolation
Solution
Bilinear interpolation
Kernel5×5 =
1
9






1 2 3 2 1
2 4 6 4 2
3 6 9 6 3
2 4 6 4 2
1 2 3 2 1





 0
-2
0.5
-1
1
20
1
1 0
-12 -2
Note: Kernel5×5 = 1
3 [1 2 3 2 1]T 1
3 [1 2 3 2 1]
David Svoboda (CBIA@FI) Image Transforms autumn 2022 29 / 50
Reconstruction/Interpolation
Solution
Completion of missing pixels (bilinear)
3x magnification + bilinear interpolation
David Svoboda (CBIA@FI) Image Transforms autumn 2022 30 / 50
1 Sampling in 1D
Comb function
Nyquist-Shannon theorem
Aliasing
Reconstruction/Interpolation
2 Sampling in 2D
Comb function
Aliasing
Reconstruction/Interpolation
3 Resampling in 1D
Design of an ideal resampling filter
4 Resampling in 2D
Elliptical Weighted Average (EWA)
David Svoboda (CBIA@FI) Image Transforms autumn 2022 31 / 50
Resampling in 1D
Let us design a 1D resampling filter
The filter should be easy to implement and fast for computation.
The filter should solve the alias problem.
David Svoboda (CBIA@FI) Image Transforms autumn 2022 32 / 50
Resampling in 1D
Design of the resampling filter
1 reconstruct the continuous signal from the discrete one
2 warp the domain of the continuous signal
3 prefilter the warped, continuous signal
4 sample this signal to produce the discrete output signal
m
f(m)
reconstruct
u
warp
reconstructed inputdiscrete input
h(u)
prefilter
warped input
h’(x)
x
continuous output
x
sample
g(n)h’’(x)
discrete output
n
f (m) . . . discrete input signal
r(u) . . . reconstruction filter
h(u) = (f ∗ r)(u) = k f (k)r(u − k) . . . reconstructed input signal
David Svoboda (CBIA@FI) Image Transforms autumn 2022 33 / 50
Resampling in 1D
Design of the resampling filter
1 reconstruct the continuous signal from the discrete one
2 warp the domain of the continuous signal
3 prefilter the warped, continuous signal
4 sample this signal to produce the discrete output signal
m
f(m)
reconstruct
u
warp
reconstructed inputdiscrete input
h(u)
prefilter
warped input
h’(x)
x
continuous output
x
sample
g(n)h’’(x)
discrete output
n
h′(x) = h(u) = h(γ−1(x)) . . . warped signal
x = γ(u) . . . mapping from one coordinate system to another one
David Svoboda (CBIA@FI) Image Transforms autumn 2022 33 / 50
Resampling in 1D
Design of the resampling filter
1 reconstruct the continuous signal from the discrete one
2 warp the domain of the continuous signal
3 prefilter the warped, continuous signal
4 sample this signal to produce the discrete output signal
m
f(m)
reconstruct
u
warp
reconstructed inputdiscrete input
h(u)
prefilter
warped input
h’(x)
x
continuous output
x
sample
g(n)h’’(x)
discrete output
n
h′′(x) = (h′ ∗ p)(x) = h′(t)p(x − t)dt . . . prefiltered signal
p(x) . . . (lowpass) prefilter
David Svoboda (CBIA@FI) Image Transforms autumn 2022 33 / 50
Resampling in 1D
Design of the resampling filter
1 reconstruct the continuous signal from the discrete one
2 warp the domain of the continuous signal
3 prefilter the warped, continuous signal
4 sample this signal to produce the discrete output signal
m
f(m)
reconstruct
u
warp
reconstructed inputdiscrete input
h(u)
prefilter
warped input
h’(x)
x
continuous output
x
sample
g(n)h’’(x)
discrete output
n
g(n) = h′′(x)III(x) . . . discrete output signal
David Svoboda (CBIA@FI) Image Transforms autumn 2022 33 / 50
Resampling in 1D
Important (implementation) notes
During the resampling process we actually never construct a
continuous signal h(u), h′(x) or h′′(x).
We pick up the individual positions in the resampled image g(n) and
look for their corresponding positions and their neighbourhood in the
original image f (m).
As the computation is inverted, we never use γ function.
We use only γ−1.
David Svoboda (CBIA@FI) Image Transforms autumn 2022 34 / 50
Resampling in 1D
Derivation of an ideal resampling filter
Computation of one sample point
g(n) = h′′
(n) = h′
(t)p(n − t)dt
= h(γ−1
(t))p(n − t)dt
= p(n − t)
k
f (k)r(γ−1
(t) − k)dt =
k
f (k)ρ(n, k)
where
ρ(n, k) = p(n − t)r(γ−1
(t) − k)dt
ρ(n, k) is called a resampling filter.
If γ is affine, we can derive: ρ(n, k) = p(γ−1(n) − k) ∗ r(γ−1(n) − k).
David Svoboda (CBIA@FI) Image Transforms autumn 2022 35 / 50
Resampling in 1D
Practical problems
If the mapping γ is not affine, the filter ρ(n, k) is space variant.
Solution (postfiltering/supersampling)
1 Reconstruct the continuous signal from the discrete input signal.
2 Warp the domain of the input signal.
3 Sample the warped signal at very high resolution to avoid alias.
4 Postfilter the signal to produce a lower resolution output signal.
Notice: The convolution is employed in the very end of this algorithm, i.e.
it is discrete and space invariant.
David Svoboda (CBIA@FI) Image Transforms autumn 2022 36 / 50
Resampling in 1D
Practical problems
If the mapping γ is not affine, the filter ρ(n, k) is space variant.
Solution (postfiltering/supersampling)
1 Reconstruct the continuous signal from the discrete input signal.
2 Warp the domain of the input signal.
3 Sample the warped signal at very high resolution to avoid alias.
4 Postfilter the signal to produce a lower resolution output signal.
Notice: The convolution is employed in the very end of this algorithm, i.e.
it is discrete and space invariant.
David Svoboda (CBIA@FI) Image Transforms autumn 2022 36 / 50
1 Sampling in 1D
Comb function
Nyquist-Shannon theorem
Aliasing
Reconstruction/Interpolation
2 Sampling in 2D
Comb function
Aliasing
Reconstruction/Interpolation
3 Resampling in 1D
Design of an ideal resampling filter
4 Resampling in 2D
Elliptical Weighted Average (EWA)
David Svoboda (CBIA@FI) Image Transforms autumn 2022 37 / 50
Resampling in 2D
Task
Design a 2D resampling filter
Obey the rules that are valid for 1D resampling filter.
The filter maps texels from texture space to screen space.
The filter might be anisotropic.
The filter should work fast.
David Svoboda (CBIA@FI) Image Transforms autumn 2022 38 / 50
Resampling in 2D
γ-mapping
γ−1 converts coordinates from screen space x = (x, y)
to texture space u = (u, v)
u
v
x
y
David Svoboda (CBIA@FI) Image Transforms autumn 2022 39 / 50
Resampling in 2D
γ-mapping
γ is projective
circular neighbourhood of one pixel is transformed into an elliptical
neighbourhood.
can be locally approximated by linear mapping – Jacobian matrix J:
J−1
=
∂u
∂x
∂u
∂y
∂v
∂x
∂v
∂y
=
Ux Uy
Vx Vy
γ : x = Ju
γ−1 : u = J−1x
David Svoboda (CBIA@FI) Image Transforms autumn 2022 40 / 50
Resampling in 2D
γ-mapping
Gaussian warped by γ gives Gaussian:
gI (u) = gI (J−1
x)
=
1
2π
e− 1
2 (J−1
x)T
J−1
x
=
=
1
2π
e
1
2 xT
(J−1T
J−1
)x
= /V −1
= J−1T
J−1
/
= |V |1/2
·
1
2π|V |1/2
e
1
2 xT
V −1
x
= /|V |1/2
= |J|/
= |J| · gV (x)
where
J−1T
J−1
. . . variance matrix (positive definite)
standard multivariate Gaussian:
gΣ(u) =
1
2π|Σ|1/2
e− 1
2 uT
Σ−1
u
David Svoboda (CBIA@FI) Image Transforms autumn 2022 41 / 50
Resampling in 2D
γ-mapping
Formation of Gaussian regions
Variance matrix V defines a conic:
V −1
=
A B/2
B/2 C
−1
= J−1T
J−1
from which the shape of zero-centered ellipse can be derived:
Q(u, v) = Au2
+ Buv + Cv2
< F
Here
A = V 2
x + V 2
y + 1
B = −2 (Ux Vx + Uy Vy )
C = U2
x + U2
y + 1
F = AC − B2
4
v
2
0u =(u ,v )0 0
u
Au + Buv + Cv = F2
F
David Svoboda (CBIA@FI) Image Transforms autumn 2022 42 / 50
Resampling in 2D
γ-mapping
Formation of Gaussian regions
Σ =
102 0
0 32 Σ =
102 20
20 32
David Svoboda (CBIA@FI) Image Transforms autumn 2022 43 / 50
Resampling in 2D
γ-mapping
BoundingBox for elliptical regions
u = ±2
CF
4AC − B2
v = ±2
AF
4AC − B2
Weight for a particular texture pixel
weight(u, v) = /ρ((x, y), (u, v))/
= e−αQ(u,v)
Notice: α . . . user defined constant (e.g. α = 2)
David Svoboda (CBIA@FI) Image Transforms autumn 2022 44 / 50
Resampling in 2D
γ-mapping
Image Pyramids (MIP map)
Size of ellipse determines level of detail in MIP map pyramid that
should be fetched from the memory.
MIP map pyramid is precomputed to avoid (pointless) repetitive
downsampling.
David Svoboda (CBIA@FI) Image Transforms autumn 2022 45 / 50
Resampling in 2D
Elliptical Weighted Average (EWA)
Implementation Notes
For each image pixel (x,y) from screen space:
1 Find corresponding point (u,v) in texture space.
2 Define the local affine transform γ.
3 Compute Jacobian J of this mapping.
4 Delineate the ellipse in texture space (find bounding box)
5 Using the ellipse size choose the two nearest MIP map levels
6 Perform the following steps in each MIP map level and combine the results
1 Build the Gaussian over the texture.
2 Evaluate direct convolution of texture image with Gaussian.
3 Get one value as a result.
7 Store the result (one value) in the screen pixel (x,y).
David Svoboda (CBIA@FI) Image Transforms autumn 2022 46 / 50
Resampling in 2D
EWA – An Example
Nearest neighbour Linear interpolation
Linear interpolation + mipMAP EWA + mipMAP
David Svoboda (CBIA@FI) Image Transforms autumn 2022 47 / 50
Resampling in 2D
EWA – An Example
Nearest neighbour Linear interpolation
Linear interpolation + mipMAP EWA + mipMAP
David Svoboda (CBIA@FI) Image Transforms autumn 2022 48 / 50
Bibliography
1 Paul Heckbert, Fundamentals of Texture Mapping and Image Warping,
Master’s thesis, UCB/CSD 89/516, CS Division, U.C. Berkeley, June 1989,
86 pp
2 McCormack, J., Perry, R.N., Farkas, K.I.; Jouppi, N.P. ”Feline: Fast
Elliptical Lines for Anisotropic Texture Mapping”, ACM SIGGRAPH, pp
243-250, August 1999
3 Pharr, M., Humphreys, G. Physically Based Rendering: From Theory to
Implementation. Morgan Kaufmann Publishers Inc., San Francisco, CA,
USA, 2004, https://www.pbr-book.org/
David Svoboda (CBIA@FI) Image Transforms autumn 2022 49 / 50
You should know the answers . . .
Show that the 2D DFT is separable transform.
Derive the complexity of 2D discrete FFT.
Explain the reciprocity of wide and narrow shapes in time and
frequency domain, respectively.
Derive (dot not formulate) the Nyquist-Shannon theorem for 2D
image data.
Show an example of the aliasing effect.
What is a prefilter?
What is the difference between a screen space and texture space?
Give an example of γ warping function both for 1D and 2D case.
What is the difference between projective and affine mappings?
Describe individual steps of EWA filter.
David Svoboda (CBIA@FI) Image Transforms autumn 2022 50 / 50