Sampling and Reconstruction

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell
 Reading
Required:
Watt, Section 14.1
Recommended:
Ron Bracewell, The Fourier Transform and Its
Applications, McGraw-Hill.
Don P. Mitchell and Arun N. Netravali,
“Reconstruction Filters in Computer Computer
Graphics ,” Computer Graphics, (Proceedings
of SIGGRAPH 88). 22 (4), pp. 221-228, 1988.

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 2
 What is an image?
We can think of an image as a function, f, from R2 to R:
f( x, y ) gives the intensity of a channel at position ( x, y )
Realistically, we expect the image only to be defined over a
rectangle, with a finite range:

f: [a,b]x[c,d]  [0,1]

A color image is just three functions pasted together. We

can write this as a “vector-valued” function:
! r ( x, y ) "
f ( x, y ) = # g ( x, y ) $
# $
#% b( x, y ) $&

We’ll focus in grayscale (scalar-valued) images for now.

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 3
 Images as functions

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 4
 Digital images
In computer graphics, we usually create or operate on
digital (discrete) images:
Sample the space on a regular grid
Quantize each sample (round to nearest integer)
If our samples are Δ apart, we can write this as:
f[i ,j] = Quantize{ f(i Δ, j Δ) }

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 5
 Motivation: filtering and resizing
What if we now want to:
smooth an image?
sharpen an image?
enlarge an image?
shrink an image?
Before we try these operations, it’s helpful
to think about images in a more
mathematical way…

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 6
 Fourier transforms
We can represent a function as a linear combination
(weighted sum) of sines and cosines.
We can think of a function in two complementary ways:
Spatially in the spatial domain
Spectrally in the frequency domain
The Fourier transform and its inverse convert between
these two domains:
$
F(s) = % f (x)e"i 2 # sx dx
Spatial "$ Frequency
$
domain i 2 " sx domain
f (x) = % F(s)e ds
#$
!

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 7
!
 Fourier transforms (cont’d)
$
F(s) = % f (x)e"i 2 # sx dx
Spatial "$ Frequency
$
domain i 2 " sx domain
f (x) = % F(s)e ds
#$

Where!do the sines and cosines come in?

!
f(x) is usually a real signal, but F(s) is generally
complex:
F(s) = Re(s) " iIm(s) = F(s) e"i 2 # s

If f(x) is symmetric, i.e., f(x) = f(-x)), then F(s) =

Re(s).
! Why?
University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 8
 1D Fourier examples

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 9
 2D Fourier transform
$ $
"i 2 # sy y
F(sx ,sy ) = %% f (x, y)e"i 2 # sx x e dxdy
Spatial $"$$"$
Frequency
domain f (x, y) = % % F(s ,s )e i 2 " sx x i 2 " sy y
e dsx dsy domain
x y
#$ #$

!
Spatial
! domain Frequency domain

f (x, y) F(sx ,sy )

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 10
! !
 2D Fourier examples

Spatial Frequency
domain domain
f (x, y) F(sx ,sy )

! !

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 11
 Convolution
One of the most common methods for filtering a
function is called convolution.
In 1D, convolution is defined as:

g(x) = f (x) * h(x)

$
= % f ( x ")h(x # x ")dx "
#$
$ )
= % f ( x ") h ( x " # x)dx "
#$

)
where h (x) = h("x)
!
University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 12
 Convolution properties
Convolution exhibits a number of basic, but
important properties.
Commutativity: a(x) " b(x) = b(x) " a(x)

Associativity: [a(x) " b(x)]" c(x) = a(x) "[b(x) " c(x)]

!
Linearity: a(x) "[k # b(x)] = k # [a(x) " b(x)]
!
a(x) " (b(x) + c(x)) = a(x) " b(x) + a(x) " c(x)

University of Texas at Austin

! CS384G - Computer Graphics Spring 2010 Don Fussell 13
 Convolution in 2D
In two dimensions, convolution becomes:
g(x, y) = f (x, y) " h(x, y)
% %
= && f ( x #, y #)h(x $ x #)(y $ y #)dx #dy #
$% $%
% % )
= && f ( x #, y #) h ( x # $ x)( y # $ y)dx #dy #
$% $%
)
where h (x, y) = h("x,"y)

! * =

f(x,y) h(x,y) g(x,y)

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 14
 Convolution theorems

Convolution theorem: Convolution in the

spatial domain is equivalent to
multiplication in the frequency domain.
f "h #F $H

Symmetric theorem: Convolution in the

frequency domain is equivalent to
multiplication
! in the spatial domain.
f "h #F $H
University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 15
 Convolution theorems
Theorem F ( f * g ) = F ( f )F ( g )
F ( fg ) = F ( f ) * F ( g )
F - 1( F * G ) = F - 1( F )F - 1(G )
F - 1( FG ) = F - 1( F ) * F - 1(G )
% %

Proof (1) F( f * g) = && f ( t ")g(t # t ")dt "e#i$ t dt

#% #%
% %
= && f ( t ")g(t # t ")e#i$ t "e#i$ (t# t " )dtdt "
#% #%
% %
= && f ( t ")g( t "")e#i$ t "e#i$ t "" dt ""dt "
#% #%
% %
#i$ t "
= & f ( t ")e dt " & g( t "")e#i$ t "" dt ""
#% #%

= F( f )F(g)
University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 16
 1D convolution theorem example

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 17
 2D convolution theorem example

f(x,y) |F(sx,sy)|

h(x,y) |H(sx,sy)|

g(x,y) |G(sx,sy)|

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 18
 The delta function
The Dirac delta function, δ(x), is a handy
tool for sampling theory.
It has zero width, infinite height, and unit
area.
It is usually drawn as:

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 19
 Sifting and shifting
For sampling, the delta function has two important
properties.
Sifting: f (x)" (x # a) = f (a)" (x # a)

Shifting: f (x) " # (x $ a) = f (x $ a)

!
University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 20
 The shah/comb function
A string of delta functions is the key to sampling.
The resulting function is called the shah or comb
function: III(x) = %"(x # nT) $

n= 0

which looks like:

!

Amazingly, the Fourier transform of the shah

function takes the same form:
$
III(s) = %" (s # ns0 )
n= 0

where so = 1/T.
University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 21
!
 Sampling
Now, we can talk about sampling.

The Fourier spectrum gets replicated by spatial sampling!

How do we recover the signal?
University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 22
 Sampling and reconstruction

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 23
 Sampling and reconstruction in 2D

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 24
 Sampling theorem
This result is known as the Sampling Theorem
and is generally attributed to Claude Shannon
(who discovered it in 1949) but was discovered
earlier, independently by at least 4 others:
A signal can be reconstructed from its samples without
loss of information, if the original signal has no energy in
frequencies at or above ½ the sampling frequency.
For a given bandlimited function, the minimum
rate at which it must be sampled is the Nyquist
frequency.

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 25
 Reconstruction filters
The sinc filter, while “ideal”, has two drawbacks:
It has large support (slow to compute)
It introduces ringing in practice
We can choose from many other filters…

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 26
 Cubic filters
Mitchell and Netravali (1988) experimented with cubic filters,
reducing them all to the following form:
#(12 " 9B " 6C) x 3 + ("18 + 12B + 6C) x 2 + (6 " 2B) x <1
1 % 3 2
r(x) = $(("B " 6C) x + (6B + 30C) x + ("12B " 48C) x + (8B + 24C) 1 ' x < 2
6%
&0 otherwise

The choice of B or C trades off between being too blurry or having too
much ringing. B=C=1/3 was their “visually best” choice.
!
The resulting reconstruction filter is often called the “Mitchell filter.”

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 27
 Reconstruction filters in 2D
We can also perform reconstruction in 2D…

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 28
 Aliasing

Sampling rate is too low

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 29
 Aliasing
What if we go below the Nyquist frequency?

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 30
 Anti-aliasing
Anti-aliasing is the process of removing the frequencies
before they alias.

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 31
Anti-aliasing by analytic prefiltering

We can fill the “magic” box with analytic pre-

filtering of the signal:

Why may this not generally be possible?

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 32
 Filtered downsampling
Alternatively, we can sample the image at a higher rate, and then filter that
signal:

We can now sample the signal at a lower rate. The whole process is called
filtered downsampling or supersampling and averaging down.

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 33
 Practical upsampling
When resampling a function (e.g., when resizing an image), you do not
need to reconstruct the complete continuous function.
For zooming in on a function, you need only use a reconstruction filter
and evaluate as needed for each new sample.
Here’s an example using a cubic filter:

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 34
 Practical upsampling
This can also be viewed as:
1.putting the reconstruction filter at the desired location
2.evaluating at the original sample positions
3.taking products with the sample values themselves
4. summing it up

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 35
 Practical downsampling
Downsampling is similar, but filter has larger support and smaller
amplitude.
Operationally:
1. Choose filter in downsampled space.
2. Compute the downsampling rate, d, ratio of new sampling rate to old
sampling rate
3. Stretch the filter by 1/d and scale it down by d
4. Follow upsampling procedure (previous slides) to compute new values

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 36
 2D resampling
We’ve been looking at separable filters:

r2D (x, y) = r1D (x)r1D (y)

!
How might you use this fact for efficient resampling in 2D?

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 37
 Next class: Image Processing
Reading:
Jain, Kasturi, Schunck, Machine Vision.
McGraw-Hill, 1995.
Sections 4.2-4.4, 4.5(intro), 4.5.5, 4.5.6, 5.1-5.4.
(from course reader)

Topics:
- Implementing discrete convolution
- Blurring and noise reduction
- Sharpening
- Edge detection

University of Texas at Austin CS384G - Computer Graphics Spring 2010 Don Fussell 38