“Linear algebra is a fantastic subject. On the one hand it is clean and beautiful.” – Gilbert Strang This wonderful branch of mathematics is both beautiful and useful. It is the cornerstone upon which signal and image processing is built. This short chapter can not be a comprehensive survey of linear algebra; it is meant only as a brief introduction and review. The ideas and presentation order are modeled after Strang’s highly recommended Linear Algebra and its Applications. At the heart of linear algebra is machinery for solving linear equations. In the simplest case, the number of unknowns equals the number of equations. For example, here are a two equations in two unknowns: 2x − y = 1

0.4 Basis 0.5 Inner Products and Projections 0.6 Linear Transforms

y

2x−y=1 (x,y)=(2,3)

x + y = 5.

(1)
x+y=5

x

There are at least two ways in which we can think of solving these equations for x and y. The ﬁrst is to consider each equation as describing a line, with the solution being at the intersection of the lines: in this case the point (2, 3), Figure 0.1. This solution is termed a “row” solution because the equations are considered in isolation of one another. This is in contrast to a “column” solution in which the equations are rewritten in vector form: 2 1 x+ −1 y = 1 1 . 5 (2)

Each equation now corresponds to a plane, and the row solution corresponds to the intersection of the planes (i.e., the intersection of two planes is a line, and that line intersects the third plane at a point: in this case, the point u = 1, v = 1, w = 2). In vector form, the equations take the form: 2 1 1 5  4  u +  −6  v +  0  w =  −2  . −2 7 2 9
       

(4)

The solution again amounts to ﬁnding values for u, v, and w that scale the vectors on the left so that their sum is equal to the vector on the right, Figure 0.3.
(5,−2,9)

Figure 0.3 “Column” solution

In the context of solving linear equations we have introduced the notion of a vector, scalar multiplication of a vector, and vector sum. In its most general form, a n-dimensional column vector is represented as: x1  x2    x =  . , .  . 
 

and on the left are the three unknowns that can also be written as a column vector: u  v . w
 

(11)

The set of nine coeﬃcients (3 rows, 3 columns) can be written in matrix form: 2 1  4 −6 −2 7


1 0 2


(12)

Matrices, like vectors, can be added and scalar multiplied. Not surprising, since we may think of a vector as a skinny matrix: a matrix with only one column. Consider the following 3 × 3 matrix: a1 A =  a4 a7


a2 a5 a8

a3 a6  . a9


(13)

The matrix cA, where c is a scalar value, is given by: ca1  ca4 cA = ca7


ca2 ca5 ca8

ca3 ca6  . ca9


(14)

And the sum of two matrices, A = B + C, is given by: a1  a4 a7


a2 a5 a8

b1 + c1 a3 a6  =  b4 + c4 b7 + c7 a9
 

b2 + c2 b5 + c5 b8 + c8

b3 + c3 b6 + c6  . b9 + c9


(15)

5

With the vector and matrix notation we can rewrite the three equations in the more compact form of Ax = b: 2  4 −2


Where the multiplication of the matrix A with vector x must be such that the three original equations are reproduced. The ﬁrst component of the product comes from “multiplying” the ﬁrst row of A (a row vector) with the column vector x as follows: u ( 2 1 1 )  v  = ( 2u + 1v + 1w ) . w
 

1 1 u 5 −6 0   v  =  −2  . 7 2 w 9
   

(16)

(17)

This quantity is equal to 5, the ﬁrst component of b, and is simply the ﬁrst of the three original equations. The full product is computed by multiplying each row of the matrix A with the vector x as follows: 2  4 −2


In its most general form the product of a m × n matrix with a n dimensional column vector is a m dimensional column vector whose ith component is:
n

aij xj ,
j=1

(19)

where aij is the matrix component in the ith row and j th column. The sum along the ith row of the matrix is referred to as the inner product or dot product between the matrix row (itself a vector) and the column vector x. Inner products are another central idea in linear algebra (more on this later). The computation for multiplying two matrices extends naturally from that of multiplying a matrix and a vector. Consider for example the following 3 × 4 and 4 × 2 matrices: a11 A =  a21 a31


That is, the i, j component of the product C is computed from an inner product of the ith row of matrix A and the j th column of matrix B. Notice that this deﬁnition is completely consistent with the product of a matrix and vector. In order to multiply two matrices A and B (or a matrix and a vector), the column dimension of A must equal the row dimension of B. In other words if A is of size m × n, then B must be of size n × p (the product is of size m × p). This constraint immediately suggests that matrix multiplication is not commutative: usually AB = BA. However matrix multiplication is both associative (AB)C = A(BC) and distributive A(B + C) = AB + AC. The identity matrix I is a special matrix with 1 on the diagonal and zero elsewhere: 1 0 ... 0 0 0 1 ... 0 0   I = . . . .. . . . . . 0 0 ... 0 1
 

(22)

Given the deﬁnition of matrix multiplication, it is easily seen that for any vector x, Ix = x, and for any suitably sized matrix, IA = A and BI = B. In the context of solving linear equations we have introduced the notion of a vector and a matrix. The result is a compact notation for representing linear equations, Ax = b. Multiplying both sides by the matrix inverse A−1 yields the desired solution to the linear equations: A−1 Ax = A−1 b Ix = A−1 b x = A−1 b (23)

A matrix A is invertible if there exists 1 a matrix B such that BA = I and AB = I, where I is the identity matrix. The matrix B is the inverse of A and is denoted as A−1 . Note that this commutative property limits the discussion of matrix inverses to square matrices. Not all matrices have inverses. Let’s consider some simple examples. The inverse of a 1 × 1 matrix A = ( a ) is A−1 = ( 1/a ); but the inverse does not exist when a = 0. The inverse of a 2 × 2
The inverse of a matrix is unique: assume that B and C are both the inverse of matrix A, then by deﬁnition B = B(AC) = (BA)C = C, so that B must equal C.
1

7

matrix can be calculated as: a b c d
−1

=

1 ad − bc

d −b −c a

,

(24)

but does not exist when ad − bc = 0. Any diagonal matrix is invertible: A=
 

a1

..

. an

  

and

A−1 = 




1/a1

..

. 1/an



  , (25)

as long as all the diagonal components are non-zero. The inverse of a product of matrices AB is (AB)−1 = B −1 A−1 . This is easily proved using the associativity of matrix multiplication. 2 The inverse of an arbitrary matrix, if it exists, can itself be calculated by solving a collection of linear equations. Consider for example a 3 × 3 matrix A whose inverse we know must satisfy the constraint that AA−1 = I: 2  4 −2


1 1 −6 0   x1 7 2


x2

x3  =  e1





e2

1 0 0  =  0 1 0  .(26) e3 0 0 1
  

This matrix equation can be considered “a column at a time” yielding a system of three equations Ax1 = e1 , Ax2 = e2 , and Ax3 = e3 . These can be solved independently for the columns of the inverse matrix, or simultaneously using the Gauss-Jordan method. A system of linear equations Ax = b can be solved by simply left multiplying with the matrix inverse A−1 (if it exists). We must naturally wonder the fate of our solution if the matrix is not invertible. The answer to this question is explored in the next section. But before moving forward we need one last deﬁnition. The transpose of a matrix A, denoted as At , is constructed by placing the ith row of A into the ith column of At . For example: A= 1 2 4 −6 1 0 1 4 t  2 −6  and A = 1 0
 

the transposes: (A + B)t = At + B t . The transpose of a product of two matrices has the familiar form (AB)t = B t At . And the transpose of the inverse is the inverse of the transpose: (A−1 )t = (At)−1 . Of particular interest will be the class of symmetric matrices that are equal to their own transpose At = A. Symmetric matrices are necessarily square, here is a 3 × 3 symmetric matrix: 2 1 A =  1 −6 4 0


Vector spaces need not be ﬁnite dimensional, R∞ is a vector space. Matrices can also make up a vector space. For example the space of 3 × 3 matrices can be thought of as R9 (imagine stringing out the nine components of the matrix into a column vector). A subspace of a vector space is a non-empty subset of vectors that is closed under vector addition and scalar multiplication. That is, the following constraints are satisﬁed: (1) the sum of any two vectors in the subspace remains in the subspace; (2) multiplication of any vector by a scalar yields a vector in the subspace. With the closure property veriﬁed, the eight properties of a vector space automatically hold for the subspace.
Example 0.1 Consider the set of all vectors in R2 whose com-

ponents are greater than or equal to zero. The sum of any two

9

vectors in this space remains in the space, but multiplication of, 1 −1 for example, the vector by −1 yields the vector 2 −2 which is no longer in the space. Therefore, this collection of vectors does not form a vector space.

Vector subspaces play a critical role in understanding systems of linear equations of the form Ax = b. Consider for example the following system: u1  u2 u3


Unlike the earlier system of equations, this system is over-constrained, there are more equations (three) than unknowns (two). A solution to this system exists if the vector b lies in the subspace of the columns of matrix A. To see why this is so, we rewrite the above system according to the rules of matrix multiplication yielding an equivalent form: v1 u1 b1 x1  u2  + x2  v2  =  b2  . v3 u3 b3
     

v1 v2  v3


x1 x2

b1  b2  = b3
 

(29)

(30)

In this form, we see that a solution exists when the scaled columns of the matrix sum to the vector b. This is simply the closure property necessary for a vector subspace. The vector subspace spanned by the columns of the matrix A is called the column space of A. It is said that a solution to Ax = b exists if and only if the vector b lies in the column space of A.
Example 0.2 Consider the following over-constrained system: Ax 1 5 2 0 4 4 u v = = b b1 b2 b3

The column space of A is the plane spanned by the vectors ( 1 5 2 )t and ( 0 4 4 )t . Therefore, the solution b can not be an arbitrary vector in R3 , but is constrained to lie in the plane spanned by these two vectors.

At this point we have seen three seemingly diﬀerent classes of linear equations of the form Ax = b, where the matrix A is either: 1. square and invertible (non-singular), 10

2. square but not invertible (singular), 3. over-constrained. In each case solutions to the system exist if the vector b lies in the column space of the matrix A. At one extreme is the invertible n×n square matrix whose solutions may be any vector in the whole of Rn . At the other extreme is the zero matrix A = 0 with only the zero vector in it’s column space, and hence the only possible solution. In between are the singular and over-constrained cases, where solutions lie in a subspace of the full vector space. The second important vector space is the nullspace of a matrix. The vectors that lie in the nullspace of a matrix consist of all solutions to the system Ax = 0. The zero vector is always in the nullspace.
Example 0.3 Consider the following system: 1 5 2 0 4 4 1 9 6 Ax u v w = = 0 0 0 0

The nullspace of A contains the zero vector ( u v w )t = ( 0 0 0 )t . Notice also that the third column of A is the sum of the ﬁrst two columns, therefore the nullspace of A also contains all vectors of the form ( u v w )t = ( c c −c )t (i.e., all vectors lying on a one-dimensional line in R3 ).

0.4 Basis Recall that if the matrix A in the system Ax = b is invertible, then left multiplying with A−1 yields the desired solution: x = A−1 b. In general it is said that a n × n matrix is invertible if it has rank n or is full rank, where the rank of a matrix is the number of linearly independent rows in the matrix. Formally, a set of vectors u1 , u2, ..., un are linearly independent if: c1u1 + c2 u2 + ... + cn un = 0 (31)

(2,2)

(−1,−1)

(2,2)

(−2,0)

is true only when c1 = c2 = ... = cn = 0. Otherwise, the vectors are linearly dependent. In other words, a set of vectors are linearly dependent if at least one of the vectors can be expressed as a sum of scaled copies of the remaining vectors. Linear independence is easy to visualize in lower-dimensional subspaces. In 2-D, two vectors are linearly dependent if they lie along a line, Figure 0.4. That is, there is a non-trivial combination of the 11

(−1,2) (2,2)

(−2,0)

Figure 0.4 Linearly dependent (top/bottom) and independent (middle).

vectors that yields the zero vector. In 2-D, any three vectors are guaranteed to be linearly dependent. For example, in Figure 0.4, the vector ( −1 2 ) can be expressed as a sum of the remaining 3 linearly independent vectors: 2 ( −2 0 ) + ( 2 2 ). In 3-D, three vectors are linearly dependent if they lie in the same plane. Also in 3-D, any four vectors are guaranteed to be linearly dependent. Linear independence is directly related to the nullspace of a matrix. Speciﬁcally, the columns of a matrix are linearly independent (i.e., the matrix is full rank) if the matrix nullspace contains only the zero vector. For example, consider the following system of linear equations: u1  u2 u3


v1 v2 v3

0 x1 w1 w2   x2  =  0  . x3 w3 0
   

(32)

Recall that the nullspace contains all vectors x such that x1 u + x2 v + x3w = 0. Notice that this is also the condition for linear independence. If the only solution is the zero vector then the vectors are linearly independent and the matrix is full rank and invertible. Linear independence is also related to the column space of a matrix. If the column space of a n × n matrix is all of Rn , then the matrix is full rank. For example, consider the following system of linear equations: u1  u2 u3


v1 v2 v3

w1 x1 b1   x2  =  b2  . w2 w3 x3 b3
   

(33)

If the the matrix is full rank, then the solution b can be any vector in R3. In such cases, the vectors u, v, w are said to span the space. Now, a linear basis of a vector space is a set of linearly independent vectors that span the space. Both conditions are important. Given an n dimensional vector space with n basis vectors v1 , ..., vn, any vector u in the space can be written as a linear combination of these n vectors: u = a1 v1 + ... + an vn . (34)

In addition, the linear independence guarantees that this linear combination is unique. If there is another combination such that: u = b1v1 + ... + bn vn , 12 (35)

which would violate the linear independence condition. While the representation is unique, the basis is not. A vector space has inﬁnitely many diﬀerent bases. For example in R2 any two vectors that do not lie on a line form a basis, and in R3 any three vectors that do not lie in a common plane or line form a basis.
Example 0.4 The vectors ( 1 0 ) and ( 0 1 ) form the canonical 2 basis for R . These vectors are both linearly independent and span the entire vector space.

Example 0.5

The vectors ( 1
3

0

0 ), ( 0

1

0 ) and ( −1

0

0)

do not form a basis for R . These vectors lie in a 2-D plane and do not span the entire vector space.

Example 0.6

The vectors ( 1

0

0 ), ( 0

−1

0 ), ( 0

0

2 ),

and ( 1 −1 0 ) do not form a basis. Although these vectors span the vector space, the fourth vector is linearly dependent on the ﬁrst two. Removing the fourth vector leaves a basis for R3 .

For notational convenience, we will often drop the cumbersome notation of Equation (1.1), and refer to the entire sequence simply as f [x]. Discrete-time signals often arise from the periodic sampling of continuous-time (analog) signals, a process that we will cover fully in later chapters. 1.2 Discrete-Time Systems In its most general form, a discrete-time system is a transformation that maps a discrete-time signal, f [x], onto a unique g[x], and is denoted as: g[x] = T {f [x]}
Example 1.1 Consider the following system: g[x] = 1 2N + 1
N

f[x]

T

g[x]

Figure 1.2 Discrete-time system

(1.3)

f [x + k].
k=−N

In this system, the kth number in the output sequence is comf[x]

puted as the average of 2N + 1 elements centered around the kth input element. As shown in Figure 1.3, with N = 2, the output value at k = 5 is computed as 1/5 times the sum of the ﬁve input
x 3 5 7

Figure 1.3 Moving Average

elements between the dotted lines. Subsequent output values are computed by sliding these lines to the right.

Although in the above example, the output at each position k depended on only a small number of input values, in general, this may not be the case, and the output may be a function of all input values.

14

1.3 Linear Time-Invariant Systems Of particular interest to us will be a class of discrete-time systems that are both linear and time-invariant. A system is said to be linear if it obeys the rules of superposition, namely: T {af1[x] + bf2[x]} = aT {f1[x]} + bT {f2[x]}, (1.4)

for any constants a and b. A system, T {·} that maps f [x] onto g[x] is shift- or time-invariant if a shift in the input causes a similar shift in the output: g[x] = T {f [x]}
Example 1.2 g[x]

The precise reason why we are particularly interested in linear time-invariant systems will become clear in subsequent chapters. But before pressing on, the concept of discrete-time systems is reformulated within a linear-algebraic framework. In order to accomplish this, it is necessary to ﬁrst restrict ourselves to consider input signals of ﬁnite length. Then, any discrete-time linear system can be represented as a matrix operation of the form: g = Mf, (1.6)

where, f is the input signal, g is the output signal, and the matrix M embodies the discrete-time linear system.

Note that according to the initial deﬁnition of the system, the output signal at x = 1 is undeﬁned (i.e., g[1] = f [0]). In the above matrix formulation we have adopted a common solution to this problem by considering the signal as wrapping around itself and setting g[1] = f [8].

Any system expressed in the matrix notation of Equation (1.6) is a discrete-time linear system, but not necessarily a time-invariant system. But, if we constrain ourselves to Toeplitz matrices, then the resulting expression will be a linear time-invariant system. A Toeplitz matrix is one in which each row contains a shifted copy of the previous row. For example, a 5 × 5 Toeplitz matrix is of the form m1 m  5  =  m4   m3 m2


M

m2 m1 m5 m4 m3

m3 m2 m1 m5 m4

m4 m3 m2 m1 m5

It is important to feel comfortable with this formulation because later concepts will be built upon this linear algebraic framework. 16

m5 m4    m3   m2  m1



(1.7)

2. Linear Time-Invariant Systems
Our interest in the class of linear time-invariant systems (LTI) is motivated by the fact that these systems have a particularly convenient and elegant representation, and this representation leads us to several fundamental tools in signal and image processing. 2.1 Space: Convolution Sum In the previous section, a discrete-time signal was represented as a sequence of numbers. More formally, this representation is in terms of the discrete-time unit impulse deﬁned as: δ[x] = 1, x = 0 0, x = 0. (2.1)

2.1 Space: Convolution Sum 2.2 Frequency: Fourier Transform

1

x

Figure 2.1 Unit impulse

Any discrete-time signal can be represented as a sum of scaled and shifted unit-impulses:
∞

Let’s now consider what happens when we present a linear timeinvariant system with this new representation of a discrete-time signal: g[x] = T {f [x]} = T
  
∞ k=−∞

f [k]δ[x − k] .


 

(2.3)

17

By the property of linearity, Equation (1.4), the above expression may be rewritten as:
∞

g[x] =
k=−∞

f [k]T {δ[x − k]}.

(2.4)

Imposing the property of time-invariance, Equation (1.5), if h[x] is the response to the unit-impulse, δ[x], then the response to δ[x−k] is simply h[x−k]. And now, the above expression can be rewritten as:
∞

g[x] =
k=−∞

f [k]h[x − k].

(2.5)

Consider for a moment the implications of the above equation. The unit-impulse response, h[x] = T {δ[x]}, of a linear time-invariant system fully characterizes that system. More precisely, given the unit-impulse response, h[x], the output, g[x], can be determined for any input, f [x]. The sum in Equation (2.5) is commonly called the convolution sum and may be expressed more compactly as: g[x] = f [x] h[x]. (2.6)

A more mathematically correct notation is (f h)[x], but for later notational considerations, we will adopt the above notation.
Example 2.2 response: Consider the following ﬁnite-length unit-impulse h[x]
h[x] 0

The next output sum at x = −1, is computed by “sliding” the unit-impulse response along the input signal and computing a similar sum.

Since linear time-invariant systems are fully characterized by convolution with the unit-impulse response, properties of such systems can be determined by considering properties of the convolution operator. For example, convolution is commutative:
∞

f [x] h[x] =
k=−∞

f [k]h[x − k], 18

let j = x − k

∞

∞

=
j=−∞

f [x − j]h[j] =

j=−∞

h[j]f [x − j] (2.7)

= h[x] f [x]. Convolution also distributes over addition:
∞

f [x] (h1[x] + h2 [x]) =
k=−∞ ∞

f [k](h1[x − k] + h2 [x − k]) f [k]h1 [x − k] + f [k]h2 [x − k]
∞

=
k=−∞ ∞

=
k=−∞

f [k]h1 [x − k] +

k=−∞

f [k]h2 [x − k] (2.8)

= f [x] h1 [x] + f [x] h2 [x].

A ﬁnal useful property of linear time-invariant systems is that a cascade of systems can be combined into a single system with impulse response equal to the convolution of the individual impulse responses. For example, for a cascade of two systems: (f [x] h1 [x]) h2 [x] = f [x] (h1 [x] h2 [x]). (2.9)

f[x]

h1[x]

h2[x]

g[x]

f[x]

h1[x] * h2[x]

g[x]

Figure 2.3 Identical LTIs

This property is fairly straight-forward to prove, and oﬀers a good exercise in manipulating the convolution sum: g1[x] = f [x] h1 [x]
∞

=
k=−∞

f [k]h1[x − k]

and, (2.10)

g2[x] = (f [x] h1 [x]) h2 [x] = g1[x] h2 [x]
∞

=
j=−∞ ∞

g1[j]h2[x − j] substituting for g1[x],
 
∞ k=−∞

=
j=−∞ ∞

=
k=−∞ ∞

f [k] 



f [k]h1 [j − k] h2 [x − j]
∞



j=−∞ ∞

=
k=−∞

= f [x] (h1 [x] h2 [x]).

f [k] 



h1 [j − k]h2[x − j] h1 [i]h2[x − i − k]




let i = j − k,

i=−∞

(2.11)

Let’s consider now how these concepts ﬁt into the linear-algebraic framework. First, a length N signal can be thought of as a point 19

where each row contains a shifted and time-reversed copy of the unit-impulse response, h[x]. The convolution matrix can then be thought of as simply transforming the basis set. As expected, this matrix is a Toeplitz matrix of the form given earlier in Equation (1.7). The reason for the time-reversal can be seen directly from the convolution sum of Equation (2.5). More speciﬁcally, the output signal g[x] at a ﬁxed x, is determined by summing the products of f [k]h[x − k] for all k. Note that the signal h[x − k] can be equivalently written as h[−k + x], which is a shifted (by x) and time-reversed (because of the −k) copy of the impulse response. Note also that when expressed in matrix form, it becomes immediately clear that the convolution sum is invertible, when h is not identically zero: g = M f and f = M −1 g.
f[x]

    , (2.12)   0   h−1 h0

Before pressing on, let’s try to combine the main ideas seen so far into a single example. We will begin by deﬁning a simple discretetime system, show that it is both linear and time-invariant, and compute its unit-impulse response
Example 2.3 Deﬁne the discrete-time system, T {·} as: g[x] = f [x + 1] − f [x − 1].

g[x]

f[x]

This system is linear because it obeys the rule of superposition: T {af1 [x] + bf2 [x]}
g[x]

In the previous section the expression of a discrete-time signal as a sum of scaled and shifted unit-impulses led us to the convolution sum. In a similar manner, we will see shortly how expressing a signal in terms of sinusoids will lead us to the Fourier transform, and then to a new way to think about discrete-time systems. The basic signal used as the building block is the sinusoid: A cos[ωx + φ], −∞ < x < ∞, (2.13)

0

−1 1

0

where A is the amplitude, ω is the frequency, and φ is the phase of the sinusoid. Shown in Figure 2.6, from top to bottom, are cos[x], cos[2x], and cos[x + π/2]. Consider next the following, seemingly unrelated, complex exponential eiωx with frequency ω, and i the √ complex value −1. This function is simply a compact notation for describing a sum of the sine and cosine function: Ae
iωx

−1 1

0

−1

= A cos(ωx) + iA sin(ωx). 21

(2.14)

Figure 2.6 f [x] = A cos[ωx + φ]

The complex exponential has a special relationship with linear time-invariant systems - the output of a linear time-invariant system with unit-impulse response h[x] and a complex exponential as input is: g[x] = eiωx h[x]
∞

=
k=−∞

h[k]eiω(x−k)
∞ iωx k=−∞

= e

h[k]e−iωk

(2.15)

Deﬁning H[ω] to be the summation component, g[x] can be expressed as: g[x] = H[ω]eiωx, (2.16)

Imaginary

H=R+I |H|

that is, given a complex exponential as input, the output of a linear time-invariant system is again a complex exponential of the same frequency scaled by some amount. 3 The scaling of the complex exponential, H[w], is called the frequency response and is generally a complex-valued function expressed in terms of its real and imaginary components: H[ω] = HR [ω] + iHI [ω], (2.17)

Intuitively, this should make perfect sense. This system simply takes an input signal and outputs a delayed copy, therefore, there is no change in the magnitude of each sinusoid, while there is a phase shift proportional to the delay, x0 .

So, why the interest in sinusoids and complex exponentials? As we will show next, a broad class of signals can be expressed as a linear combination of complex exponentials, and analogous to the impulse response, the frequency response completely characterizes the system. Let’s begin by taking a step back to the more familiar sinusoids, and then work our way to the complex exponentials. Any periodic discrete-time signal, f [x], can be expressed as a sum of scaled, phase-shifted sinusoids of varying frequencies: f [x] = 1 2π
π

In the language of linear algebra, the sinusoids are said to form a basis for the set of periodic signals, that is, any periodic signal can be written as a linear combination of the sinusoids. Recall 23

that in deriving the convolution sum, the basis consisted of shifted copies of the unit-impulse. But note now that this new basis is not ﬁxed because of the phase term, φk . It is, however, possible to rewrite the Fourier series with respect to a ﬁxed basis of zerophase sinusoids. With the trigonometric identity cos(A + B) = cos(A) cos(B)−sin(A) sin(B), the Fourier series of Equation (2.18) may be rewritten as: f [x] = = = 1 2π 1 2π 1 2π
π

ck cos[kx + φk ]
k=−π π

ck cos[φk ] cos[kx] + ck sin[φk ] sin[kx]
k=−π π

ak cos[kx] + bk sin[kx]
k=−π

(2.19)

In this expression, the constants ak and bk are the Fourier coefﬁcients and are determined by the Fourier transform. In other words, the Fourier transform simply provides a means for expressing a signal in terms of the sinusoids. The Fourier coeﬃcients are given by:
∞ ∞

ak =
j=−∞

f [j] cos[kj] and bk =
j=−∞

f [j] sin[kj]

(2.20)

Notice that the Fourier coeﬃcients are determined by projecting the signal onto each of the sinusoidal basis. That is, consider both the signal f [x] and each of the sinusoids as T -dimensional vectors, f and b, respectively. Then, the projection of f onto b is: f0 b0 + f1 b1 + ... =
j

fj bj ,

(2.21)

where the subscript denotes the jth entry of the vector. Often, a more compact notation is used to represent the Fourier series and Fourier transform which exploits the complex exponential and its relationship to the sinusoids: eiωx = cos(ωx) + i sin(ωx), (2.22)

√ where i is the complex value −1. Under the complex exponential notation, the Fourier series and transform take the form: f [x] = 1 π ck eikx 2π k=−π
∞

Comparing the Fourier transform (Equation (2.24)) with the frequency response (Equation (2.16)) we see now that the frequency response of a linear time-invariant system is simply the Fourier transform of the unit-impulse response:
∞

Both the impulse response and the magnitude of the frequency response are shown below, where for clarity, the frequency response was drawn as a continuous function.

26

Space (h[x])
0.5

Frequency (|H[ω]|)
1

0

−0.5 −1 0 1

0 −pi

0

pi

This system is an (approximate) diﬀerentiator, and can be seen from the deﬁnition of diﬀerentiation: df (x) dx =
ε→0

lim

f (x + ε) − f (x − ε) , ε

where, in the case of the system T {·}, ε is given by the distance between samples of the discrete-time signal f [x]. Let’s see now if we can better understand the frequency response of this system, recall that the magnitude was given by | sin(ω)| and the phase by π . Consider now the derivative of a ﬁxed fre2 quency sinusoid sin(ωx), diﬀerentiating with respect to x gives ω cos(ωx) = ω sin(ωx − π/2). Note that diﬀerentiation causes a phase shift of π/2 and a scaling by the frequency of the sinusoid. Notice that this is roughly in-line with the Fourier transform, the diﬀerence being that the amplitude is given by | sin(ω)| instead of ω. Note though that for small ω, | sin(ω)| ≈ ω. This discrepancy makes the system only an approximate, not a perfect, diﬀerentiator.

Linear time-invariant systems can be fully characterized by their impulse, h[x], or frequency responses, H[ω], both of which may be used to determine the output of the system to any input signal, f [x]: g[x] = f [x] h[x] and G[ω] = F [ω]H[ω], (2.26)

where the output signal g[x] can be determined from its Fourier transform G[ω], by simply applying the inverse Fourier transform. This equivalence illustrates an important relationship between the space and frequency domains. Namely, convolution in the space domain is equivalent to multiplication in the frequency domain. This is fairly straight-forward to prove: g[x] = f [x] h[x]
∞ ∞

Notice that this equation embodies both the Fourier transform and the Fourier series of Equation (2.24). The above form is the Fourier transform, and the Fourier series is gotten by left-multiplying with the inverse of the matrix, M −1 F = f .

28

3. Sampling: Continuous to Discrete (and back)
It is often more convenient to process a continuous-time signal with a discrete-time system. Such a system may consist of three distinct stages: (1) the conversion of a continuous-time signal to a discrete-time signal (C/D converter); (2) the processing through a discrete-time system; and (3) the conversion of the output discretetime signal back to a continuous-time signal (D/C converter). Earlier we focused on the discrete-time processing, and now we will concentrate on the conversions between discrete- and continuoustime signals. Of particular interest is the somewhat remarkable fact that under certain conditions, a continuous-time signal can be fully represented by a discrete-time signal! 3.1 Continuous to Discrete: Space
f[x]

for integer values x. In this expression, the quantity T is the sampling period. In general, continuous-time signals will be denoted with rounded parenthesis (e.g., f (·)), and discrete-time signals with square parenthesis (e.g., f [·]). This sampling operation may be considered as a multiplication of the continuous time signal with an impulse train, Figure 3.2. The impulse train is deﬁned as:
∞

D/C

g(x)

Figure 3.1 Processing block diagram

s(x) =
k=−∞

δ(x − kT ),

(3.2)

f(x)

where δ(·) is the unit-impulse, and T is the sampling period - note that the impulse train is a continuous-time signal. Multiplying the impulse train with a continuous-time signal gives a sampled signal: fs (x) = f (x)s(x), (3.3)

x f[x]

Figure 3.2 Sampling: space

Note that the sampled signal, fs (x), is indexed on the continuous variable x, while the ﬁnal discrete-time signal, f [x] is indexed on the integer variable x. It will prove to be mathematically convenient to work with this intermediate sampled signal, fs (x).

29

3.2 Continuous to Discrete: Frequency
F(w)

w −w n S(w) wn

In the space domain, sampling was described as a product between the impulse train and the continuous-time signal (Equation (3.3)). In the frequency domain, this operation amounts to a convolution between the Fourier transform of these two signals: Fs (ω) = F (ω) S(ω) (3.4)

w −w s 0 Fs(w) ws

w w s− w n

Figure 3.3 Sampling: no aliasing

For example, shown in Figure 3.3 (from top to bottom) are the Fourier transforms of the continuous-time function, F (ω), the impulse train, S(ω), itself an impulse train, and the results of convolving these two signals, Fs (ω). Notice that the Fourier transform of the sampled signal contains multiple (yet exact) copies of the Fourier transform of the original continuous signal. Note however the conditions under which an exact replica is preserved depends on the maximum frequency response ωn of the original continuous-time signal, and the sampling interval of the impulse train, ωs which, not surprisingly, is related to the sampling period T as ωs = 2π/T . More precisely, the copies of the frequency response will not overlap if: ωn < ω s − ω n or (3.5)

F(w)

w −w n S(w) wn

ωs > 2ωn ,

w −w s 0 Fs(w) ws

The frequency ωn is called the Nyquist frequency and 2ωn is called the Nyquist rate. Shown in Figure 3.4 is another example of this sampling process in the frequency domain, but this time, the Nyquist rate is not met, and the copies of the frequency response overlap. In such a case, the signal is said to be aliased. Not surprisingly, the Nyquist rate depends on both the characteristics of the continuous-time signal, and the sampling rate. More precisely, as the maximum frequency, ωn , of the continuous-time signal increases, the sampling period, T must be made smaller (i.e., denser sampling), which in turn increases ωs , preventing overlap of the frequency responses. In other words, a signal that changes slowly and smoothly can be sampled fairly coarsely, while a signal that changes quickly requires more dense sampling. 3.3 Discrete to Continuous If the Nyquist rate is met, then a discrete-time signal fully characterizes the continuous-time signal from which it was sampled. On the other hand, if the Nyquist rate is not met, then the sampling leads to aliasing, and the discrete-time signal does not accurately represent its continuous-time counterpart. In the former case, it 30

w

Figure 3.4 Sampling: aliasing

is possible to reconstruct the original continuous-time signal, from the discrete-time signal. In particular since the frequency response of the discrete-time signal contains exact copies of the original continuous-time signals frequency response, we need only extract one of these copies, and inverse transform the result. The result will be identical to the original signal. In order to extract a single copy, the Fourier transform of the sampled signal is multiplied by an ideal reconstruction ﬁlter as shown in Figure 3.5. This ﬁlter has unit value between the frequencies −π/T to π/T and is zero elsewhere. This frequency band is guaranteed to be greater than the Nyquist frequency, ωn (i.e., ωs = 2π/T > 2ωn , so that π/T > ωn ). In the space domain, this ideal reconstruction ﬁlter has the form: sin(πx/T ) , πx/T

Fs(w)

w pi/T

Figure 3.5 Reconstruction

h(x) =

(3.6)
1

and is often referred to as the ideal sync function. Since reconstruction in the frequency domain is accomplished by multiplication with the ideal reconstruction ﬁlter, we could equivalently reconstruct the signal by convolving with the ideal sync in the space domain.

a sinusoid with frequency ω0 . We will eventually be interested in sampling this function and seeing how the eﬀects of aliasing are manifested. But ﬁrst, let’s compute the Fourier transform of this signal:
∞

F (ω)

=
k=−∞ ∞

f (k)e−iωk

=
k=−∞ ∞

cos(ω0 k)(cos(ωk) − i sin(ωk))

=
k=−∞

cos(ω0 k) cos(ωk) − i cos(ω0 k) sin(ωk)

First let’s consider the product of two cosines. It is easy to show from basic trigonometric identities that cos(A) cos(B) = 0 when A = B, and is equal to π when |A| = |B|. Similarly, one can show that cos(A) sin(B) = 0 for all A and B. So, the Fourier transform of cos(ω0 x) = π for |ω| = ω0 , and is 0 otherwise (see below). If the sampling rate is greater than 2ω0 , then there will be no aliasing, but if the sampling rate is less than 2ω0 , then the reconstructed signal will be of the form cos((ωs − ω0 )x), that is, the reconstructed signal will be appear as a lower frequency sinusoid - it will be aliased.

31

F(w)

w −w 0 w0

Sampling

No Aliasing

Aliasing Fs(w)

Fs(w)

w −w 0 w0 −w 0 w0

w

We will close this chapter by drawing on the linear algebraic framework for additional intuition on the sampling and reconstruction process. First we will need to restrict ourselves to the sampling of an already sampled signal, that is, consider a m-dimensional signal sub-sampled to a n-dimensional signal. We may express this operation in matrix form as follows: 1 0 0 0 0 1 g1  .  .   .  = . . .  0 0 0 gn 1 0 0
  

where the subscripts denote the vector and matrix dimensions, and in this example n = m/2. Our goal now is to determine when it is possible to reconstruct the signal f , from the sub-sampled signal g. The Nyquist sampling theory tells us that if a signal is band-limited (i.e., can be written as a sum of a ﬁnite number of sinusoids), then we can sample it without loss of information. We can express this constraint in matrix notation: fm = Bm×n wn , (3.8)

where the columns of the matrix B contains the basis set of sinusoids - in this case the ﬁrst n sinusoids. Substituting into the above sampling equation gives: gn = Sn×m Bm×n wn = Mn×n wn .

(3.9)

If the matrix M is invertible, then the original weights (i.e., the representation of the original signal) can be determined by simply left-multiplying the sub-sampled signal g by M −1 . In other 32

words, Nyquist sampling theory can be thought of as simply a matrix inversion problem. This should not be at all surprising, the trick to sampling and perfect reconstruction is to simply limit the dimensionality of the signal to at most twice the number of samples.

33

4. Digital Filter Design
Recall that the class of linear time-invariant systems are fully characterized by their impulse response. More speciﬁcally, the output of a linear time-invariant system to any input f [x] can be determined via a convolution with the impulse response h[x]: g[x] = f [x] h[x]. (4.1)

Therefore the ﬁlter h[x] and the linear-time invariant system are synonymous. In the frequency domain, this expression takes on the form: G[ω] = F [ω]H[ω]. (4.2)

In other words, a ﬁlter modiﬁes the frequencies of the input signal. It is often the case that such ﬁlters pass certain frequencies and attenuate others (e.g., a lowpass, bandpass, or highpass ﬁlters). The design of such ﬁlters consists of four basic steps: 1. 2. 3. 4. choose the desired frequency response choose the length of the ﬁlter deﬁne an error function to be minimized choose a minimization technique and solve

The choice of frequency response depends, of course, on the designers particular application, and its selection is left to their discretion. We will however provide some general guidelines for choosing a frequency response that is amenable to a successful design. In choosing a ﬁlter size there are two conﬂicting goals, a large ﬁlter allows for a more accurate match to the desired frequency response, however a small ﬁlter is desirable in order to minimize computational demands 4 . The designer should experiment with varying size ﬁlters until an equitable balance is found. With a frequency response and ﬁlter size in hand, this chapter will provide the computational framework for realizing a ﬁnite length ﬁlter that “best” approximates the speciﬁed frequency response. Although there are numerous techniques for the design of digital ﬁlters we will cover only two such techniques chosen for their simplicity and generally good performance (see Digital Filter Design by T.W. Parks and C.S. Burrus for a full coverage of many other approaches).

4

In multi-dimensional ﬁlter design, separability is also a desirable property.

34

4.1 Choosing a Frequency Response A common class of ﬁlters are bandpass in nature, that is, they pass certain frequencies and attenuate others. An ideal lowpass, bandpass, and highpass ﬁlter are illustrated in Figure 4.1 Shown is the magnitude of the frequency response in the range [0, π], since we are typically interested in designing real-valued, linear-phase ﬁlters, we need only specify one-half of the magnitude spectrum (the response is symmetric about the origin). The responses shown in Figure 4.1 are often referred to as brick wall ﬁlters because of their abrupt fall-oﬀ. A ﬁnite-length realization of such a ﬁlter produces undesirable “ringing” known as Gibbs phenomena as shown below in the magnitude of the frequency response of ideal lowpass ﬁlters of length 64, 32, and 16 (commonly referred to as the ﬁlter tap size). 64 taps
1 1
| H(w) |

Passband Stopband w 0 pi

32 taps
1

16 taps

Figure 4.1 Ideal lowpass, bandpass, and highpass

0.5

0.5

0.5

0 0

pi/2

pi

0 0

pi/2

pi

0 0

pi/2

pi

| H(w) |

These eﬀects are particularly undesirable because neighboring frequencies may be passed or attenuated by wildly varying amounts, leading to general instabilities. To counter this problem, the designer is resigned to sacriﬁcing the ideal response for a “softer” frequency response, Figure 4.2. Such a frequency response is amenable to a small ﬁnite-length ﬁlter free of signiﬁcant ringing artifacts. The speciﬁc functional form of the soft falloﬀ is somewhat arbitrary, however one popular form is a raised cosine. In its most general form, the frequency response takes the following form, where the bandpass nature of the response is controlled through ω0 , ω1 , ∆ω0 , and ∆ω1 .
  0,   1  [1 − cos(π(ω − ω )/∆ω )] ,  2 0 0       

In general, a larger transition width (i.e., ∆ω0 , and ∆ω1 ) allows for smaller ﬁlters with minimal ringing artifacts. The tradeoﬀ, of course, is that a larger transition width moves the frequency response further from the ideal brick-wall response. In specifying only half the frequency response (from [0, π]) we are implicitly imposing a symmetric frequency response and thus assuming that the desired ﬁlter is symmetric. 35

0.5

ω0

ω1

0 0

pi

Figure 4.3 Frequency response

4.2 Frequency Sampling Once the desired frequency response has been chosen, the more diﬃcult task of designing a ﬁnite-length ﬁlter that closely approximates this response begins. In general this problem is hard because by restricting ourselves to a small ﬁnite-length ﬁlter we are in eﬀect asking to ﬁt an arbitrarily complex function (the desired frequency response) with the sum of a small number of sinusoids. We begin with the most straight-forward design method - frequency sampling. Our goal is to design a ﬁlter h whose Fourier transform best approximates the speciﬁed response H, that is: F (h) = H (4.3)

In other words, the design simply involves inverse Fourier transforming the speciﬁed frequency response. This series of steps can be made more explicit and practical for computer implementation by expressing the initial constraint (Equation (4.3)) in matrix notation: Mh = H (4.5)

where H is the n-dimensional sampled frequency response, M is the n × n Fourier matrix (Equation (2.28)), and n is the chosen ﬁlter size. The ﬁlter h can be solved for by left multiplying both sides of Equation (4.5) by the inverse of the matrix F : h = M −1 H. (4.6)

Since the matrix M is square, this design is equivalent to solving for n unknowns (the ﬁlter taps) from n linearly independent equations. This fact illustrates the shortcomings of this approach, namely, that this method produces a ﬁlter with a frequency response that exactly matches the sampled response, but places no constraints on the response between sampled points. This restriction can often lead to poor results that are partially alleviated by the least-squares design presented next.

36

4.3 Least-Squares Our goal is to design a ﬁlter h that “best” approximates a speciﬁed frequency response. As before this constraint can be expressed as: M h = H, (4.7)

where M is the N × n Fourier matrix (Equation (2.28)), H is the N sampled frequency response, and the ﬁlter size is n. Note that unlike before, this equation is over constrained, having n unknowns in N > n equations. We can solve this system of equations in a least-squares sense by ﬁrst writing a squared error function to be minimized: E(h) = | M h − H |2 (4.8)

Shown in Figure 4.4 are a set of lowpass, bandpass, and highpass 16-tap ﬁlters designed using this technique. In this design, the frequency response was of the form given in Equation (4.3), with start and stop bands of [ω0 , ω1] = [0, π/2], [π/4, 3π/4], and [π/2, π], and a transition width of ∆ω0 = ∆ω1 = π/4. The frequency response was sampled at a rate of N = 512. Finally, note that the previous frequency sampling design is equivalent to the least-squares design when the sampling of the Fourier basis is the same as the ﬁlter size (i.e., N = n). 4.4 Weighted Least-Squares One drawback of the least-squares method is that the error function (Equation (4.8)) is uniform across all frequencies. This is easily rectiﬁed by introducing a weighting on the least-squares error function: E(h) = W | M h − H |2
5

0.5

0 0 1

pi/2

pi

0.5

0 0 1

pi/2

pi

0.5

0 0

pi/2

pi

(4.11)

Figure 4.4 Leastsquares: lowpass, bandpass, and highpass

Because of Parseval’s theorem (which amounts to the orthonormality of the Fourier transform), the minimal error in the frequency domain equates to a minimal error in the space domain.

37

where W is a diagonal weighting matrix. That is, the diagonal of the matrix contains the desired weighting of the error across frequency. As before, we minimize by diﬀerentiating with respect to h: dE(h) dh = 2M tW | M h − H | = 2M tW M h − 2W M t H, then set equal to zero, and solve:
1

(4.12)

h = (M t W M )−1 M t W H.

(4.13)

0.5

Note that this solution will be equivalent to the original leastsquares solution (Equation (4.10)) when W is the identity matrix (i.e., uniform weighting).
pi/2 pi

0 0 1

0.5

0 0

Shown in Figure 4.5 is a comparison of an 8-tap lowpass ﬁlter designed with a uniform weighting, and with a weighting that 1 emphasizes the errors in the low frequency range, W (ω) = (|ω|+1)8 . Note that in the case of the later, the errors in the low frequencies are smaller, while the errors in the high frequencies have increased.
pi/2 pi

Figure 4.5 Least-squares and weighted least squares

38

5. Photons to Pixels

5.1 Pinhole Camera The history of the pinhole camera (or camera obscura) dates back as early as the ﬁfth century B.C., and continues to be popular today among students, artists, and scientists. The Chinese philosopher Mo Ti is believed to be the ﬁrst to notice that objects reﬂect light in all directions and that the light rays that pass through a small hole produce an inverted image. In its simplest form a pinhole camera is a light-tight box with a tiny hole in one end and a photo-sensitive material on the other. Remarkably, this simple device is capable of producing a photograph. However, the pinhole camera is not a particularly eﬃcient imaging system (often requiring exposure times as long as several hours) and is more popular for its artistic value than for its practical value. Nevertheless, the pinhole camera is convenient because it aﬀords a simple model of more complex imaging systems. That is, with a pinhole camera model, the projection of points from the three-dimensional world onto the two-dimensional sensor takes on a particularly simple form. Denote a point in the three-dimensional world as a column vector, P = ( X Y Z )t and the projection of this point onto the two dimensional image plane as p = ( x y )t . Note that the world and image points are expressed with respect to their own coordinate systems, and for convenience, the image coordinate system is chosen to be orthogonal to the Z-axis, i.e., the origins of the two systems are related by a one-dimensional translation along the Z−axis or optical axis. It is straight-forward to show from a similar triangles argument that the relationship between the world and image point is: dX x=− Z and dY y=− , Z (5.1)

5.1 Pinhole Camera 5.2 Lenses 5.3 CCD

Figure 5.1 Pinhole image formation

Y P Z

y

p x X

where d is the displacement of the image plane along the Z-axis 6 These equations are frequently referred to as the perspective projection equations. Although non-linear in their nature, the perspective projection equations may be expressed in matrix form
The value d in Equation (5.1) is often referred to as the focal length. We do not adopt this convention primarily because it is a misnomer, under the pinhole model all points are imaged in perfect focus.
6

Figure 5.2 Perspective projection

39

using the homogeneous equations:


xs s

 ys 



−d 0 =  0 −d 0 0


X 0 0   Y 0 0 , Z 1 0 1
  

(5.2)

where the ﬁnal image coordinates are given by ( x

y )t = ( x s s

ys t s ) .

An approximation to the above perspective projection equations is orthographic projection, where light rays are assumed to travel from a point in the world parallel to the optical axis until they intersect the image plane. Unlike the pinhole camera and perspective projection equations, this model is not physically realizable and is used primarily because the projection equations take on a particularly simple linear form:
Y P p Z x X y

x=X And in matrix form: x y =

and

y = Y.

(5.3)

Figure 5.3 Orthographic projection

1 0

X 0 0   Y 1 0 Z
 

(5.4)

Orthographic projection is a reasonable approximation to perspective projection when the diﬀerence in depth between points in the world is small relative to their distance to the image plane. In the special case when all the points lie on a single frontal-parallel d surface relative to the image plane (i.e., Z is constant in Equation (5.1)), the diﬀerence between perspective and orthographic is only a scale factor. 5.2 Lenses It is important to remember that both the perspective and orthographic projection equations are only approximations of more complex imaging systems. Commercial cameras are constructed with a variety of lenses that collect and focus light onto the image plane. That is, light emanates from a point in the world in all directions and, whereas a pinhole camera captures a single light ray, a lens collects a multitude of light rays and focuses the light to a small region on the image plane. Such complex imaging systems are often described with the simpler thin-lens model. Under the thin-lens model the projection of the central or principal ray obeys the rules of perspective projection, Equation (5.1): the point P = ( X Y Z )t is projected onto the image plane centered about the point ( x y )t = ( −dX −dY )t . If the point P Z Z 40

Y P Z
y

Figure 5.4 Thin lens

is in perfect focus, then the remaining light rays captured by the lens also strike the image plane at the point p. A point is imaged in perfect focus if its distance from the lens satisﬁes the following thin-lens equation: 1 1 + Z d = 1 , f (5.5)

where d is the distance between the lens and image plane along the optical axis, and f is the focal length of the lens. The focal length is deﬁned to be the distance from the lens to the image plane such that the image of an object that is inﬁnitely far away is imaged in perfect focus. Points at a depth of Zo = Z are imaged onto a small region on the image plane, often modeled as a blurred circle with radius r: r =
1 f

R 1 − Zo

1 1 − f Zo

−

1 , d

(5.6)

where R is the radius of the lens. Note that when the depth of a point satisﬁes Equation (5.5), the blur radius is zero. Note also that as the lens radius R approaches 0 (i.e., a pinhole camera), the blur radius also approaches zero for all points independent of its depth (referred to as an inﬁnite depth of ﬁeld). Alternatively, the projection of each light ray can be described in the following more compact matrix notation: l2 α2 = 1
1 −R n2 −n1 n2

0
n1 n2

l1 α1

,

(5.7)

where R is the radius of the lens, n1 and n2 are the index of refraction for air and the lens material, respectively. l1 and l2 are the height at which a light ray enters and exits the lens (the thin lens idealization ensures that l1 = l2 ). α1 is the angle between the entering light ray and the optical axis, and α2 is the angle between the exiting light ray and the optical axis. This formulation is particularly convenient because a variety of lenses can be described in matrix form so that a complex lens train can then be modeled as a simple product of matrices.
Y

Image formation, independent of the particular model, is a threedimensional to two-dimensional transformation. Inherent to such a transformation is a loss of information, in this case depth information. Speciﬁcally, all points of the form Pc = ( cX cY cZ )t , for any c ∈ R, are projected to the same point ( x y )t - the projection is not one-to-one and thus not invertible. In addition to this geometric argument for the non-invertibility of image formation, a similarly straight-forward linear algebraic argument holds. 41

y

P 1 Z

P 2

P 3 p x X

Figure 5.5 Non-invertible projection

In particular, we have seen that the image formation equations may be written in matrix form as, p = Mn×m P , where n < m (e.g., Equation (5.2)). Since the projection is from a higher dimensional space to a lower dimensional space, the matrix M is not invertible and thus the projection is not invertible. 5.3 CCD To this point we have described the geometry of image formation, how light travels through an imaging system. To complete the image formation process we need to discuss how the light that strikes the image plane is recorded and converted into a digital image. The core technology used by most digital cameras is the chargecoupled device (CCD), ﬁrst introduced in 1969. A basic CCD consists of a series of closely spaced metal-oxide-semiconductor capacitors (MOS), each one corresponding to a single image pixel. In its most basic form a CCD is a charge storage and transport device: charge is stored on the MOS capacitors and then transported across these capacitors for readout and subsequent transformation to a digital image. More speciﬁcally, when a positive voltage, V , is applied to the surface of a P-type MOS capacitor, positive charge migrates toward ground. The region depleted of positive charge is called the depletion region. When photons (i.e., light) enter the depletion region, the electrons released are stored in this region. The value of the stored charge is proportional to the intensity of the light striking the capacitor. A digital image is subsequently formed by transferring the stored charge from one depletion region to the next. The stored charge is transferred across a series of MOS capacitors (e.g., a row or column of the CCD array) by sequentially applying voltage to each MOS capacitor. As charge passes through the last capacitor in the series, an ampliﬁer converts the charge into a voltage. An analog-to-digital converter then translates this voltage into a number (i.e., the intensity of an image pixel).

Light

V

Depletion region

MOS Ground

Figure 5.6 MOS capacitor

42

6. Point-Wise Operations

6.1 Lookup Table The internal representation of a digital image is simply a matrix of numbers representing grayscale or color values. But when an image is displayed on a computer monitor we typically do not see a direct mapping of the image. An image is ﬁrst passed through a lookup table (LUT) that maps the image intensity values to brightness values, Figure 6.1. If the lookup table is linear with unit slope and zero intercept then the image is directly mapped to the display, otherwise, the displayed image will not be an exact representation of the underlying image. For example, most computer monitors intentionally impose a non-linear LUT of the general form D = I α (i.e., gamma correction), where α is a real number, and D and I are the displayed and image values. A variety of interesting visual eﬀects can be achieved by simple manipulations of the functional form of the LUT. Keep in mind though that in manipulating the lookup table, the underlying image is left untouched, it is only the mapping from pixel value to display brightness that is being eﬀected. 6.2 Brightness/Contrast Perhaps the most common and familiar example of a LUT manipulation is to control the brightness or darkness of an image as shown in Figure 6.3. The bright and dark images of Einstein were created by passing the middle image through the LUTs shown in the same ﬁgure. The functional form of the LUT is a unit-slope line with varying intercepts: g(u) = u + b, with the image intensity values u ∈ [0, 1]. A value of b > 0 results in a brightening of the image and b < 0 a darkening of the image. Another common manipulation is that of controlling the contrast of an image as shown in Figure 6.4. The top image is said to be high contrast and the bottom image low contrast, with the corresponding LUTs shown in the same ﬁgure. The functional form of these LUTs is linear: g(u) = mu + b, where the relationship between the slope and intercept is b = 1/2(1 − m). The image contrast is increased with a large slope and negative intercept (in the limit, m → ∞ and b → −∞), and the contrast is reduced with a small slope and positive intercept (m → 0 and b → 1/2). The image is inverted with m = −1 and b = 1, that is white is mapped to black, and black is mapped to white. Of course, an image can be both contrast enhanced and brightened or darkened by simply passing the image through two (or more) LUTs. 43

Autoscaling is a special case of contrast enhancement where the minimum image intensity value is mapped to black and the maximum value is mapped to white, Figure 6.2. Autoscaling maximizes the contrast without saturating at black or white. The problem with this sort of autoscaling is that a few stray pixels can dictate the contrast resulting in a low-contrast image. A less sensitive approach is to sort the intensity values in increasing order and map the 1% and 99% intensity values to black and white, respectively. Although this will lead to a small amount of saturation, it rarely fails to produce a high-contrast image. 6.3 Gamma Correction Typically, high contrast images are visually more appealing. However a drawback of linear contrast enhancement described above is that it leads to saturation at both the low and high end of the intensity range. This may be avoided by employing a non-linear contrast adjustment scheme, also realizable as a lookup table (LUT) manipulation. The most standard approach is gamma correction, where the LUT takes the form: g(u) = uα ,
Dark

Bright

(6.1)

Figure 6.3 Brightness

where α > 1 increases contrast, and α < 1 reduces contrast. Shown in Figure 6.5 are contrast enhanced (top: α = 2) and contrast reduced (bottom: α = 1/2) images. Note that with the intensity values u scaled into the range [0, 1], black (0) and white (1) are mapped to themselves. That is, there is no saturation at the low or high end. Gamma correction is widely used in a number of devices because it yields reasonable results and is easily parameterized. One drawback to this scheme is that the gray values are mapped in an asymmetric fashion with respect to midlevel gray (0.5). This may be alleviated by employing a sigmoidal non-linearity of the form g(u) = 1 . 1 + eαu+β (6.2)

High

In order that g(u) be bound by the interval [0, 1], it must be scaled as follows: c2(g(u) − c1), where c1 = 1/(1 + eβ ) and c2 = (1 + e−α+β ) − c1. This non-linearity, with its two degrees of freedom, is more versatile and can produce a more balanced contrast enhancement. Shown in Figure 6.6 is a contrast enhanced image with α = 12 and β = 6.
Low

Figure 6.4 Contrast

44

6.4 Quantize/Threshold A digital image, by its very nature, is quantized to a discrete number of intensity values. For example an image quantized to 8-bits contains 28 = 256 possible intensity values, typically in the range [0, 255]. An image can be further quantized to a lower number of bits (b) or intensity values (2b ). Quantization can be accomplished by passing the image through a LUT containing a step function, Figure 6.7, where the number of steps governs the number of intensity values. Shown in Figure 6.7 is an image of Einstein quantized to ﬁve intensity values, notice that all the subtle variations in the curtain and in his face and jacket have been eliminated. In the limit, when an image is quantized to one bit or two intensity values, the image is said to be thresholded. Shown in Figure 6.8 is a thresholded image of Einstein and the corresponding LUT, a twostep function. The point at which the step function transitions from zero to one is called the threshold and can of course be made to be any value (i.e., slid left or right). 6.5 Histogram Equalize The intensity values of a typical image are often distributed unevenly across the full range of 0 to 255 (for an 8-bit image), with most the mass near mid-gray (128) and falling oﬀ on either side, Figure 6.9. An image can be transformed so that the distribution of intensity values is ﬂat, that is, each intensity value is equally represented in the image. This process is known has histogram equalization 7. Although it may not be immediately obvious an image is histogram equalized by passing it through a LUT with the functional form of the cumulative distribution function. More speciﬁcally, deﬁne N (u) as the number of pixels with intensity value u, this is the image histogram and a discrete approximation to the probability distribution function. Then, the cumulative distribution function is deﬁned as:
u

Figure 6.5 Contrast: Gamma

Figure 6.6 Contrast: Sigmoid

C(y) =
i=0

N (i),

(6.3)

that is, C(u) is the number of pixels with intensity value less than or equal to u. Histogram equalization then amounts to simply inserting this function into the LUT. Shown in Figure 6.9 is Einstein before and after histogram equalization. Notice that the eﬀect is similar to contrast enhancement, which intuitively should make sense since we increased the number of black and white pixels.
Why anyone would want to histogram equalize an image is a mystery to me, but here it is in case you do.
7

45

In all of these examples the appearance of an image was altered by simply manipulating the LUT, the mapping from image intensity value to display brightness value. Such a manipulation leaves the image content intact, it is a non-destructive operation and thus completely invertible. These operations can be made destructive by applying the LUT operation directly to the image. For example an image can be brightened by adding a constant to each pixel, and then displaying with a linear LUT. Since such an operation is destructive it may not be inveritble, for example when brightening an 8-bit image, all pixels that exceed the value 255 will be truncated to 255.
Figure 6.7 Quantize

Figure 6.8 Threshold

15000 10000 5000 0 0

100

200

15000 10000 5000 0 0

100

200

46

Figure 6.9 Histogram equalize

7. Linear Filtering

7.1 Convolution The one-dimensional convolution sum, Equation (2.5), formed the basis for much of our discussion on discrete-time signals and systems. Similarly the two-dimensional convolution sum will form the basis from which we begin our discussion on image processing and computer vision. The 1-D convolution sum extends naturally to higher dimensions. Consider an image f [x, y] and a two-dimensional ﬁlter h[x, y]. The 2-D convolution sum is then given by:
∞ ∞

In 1-D, the intuition for a convolution is that of computing inner products between the ﬁlter and signal as the ﬁlter “slides” across the signal. The same intuition holds in 2-D. Inner products are computed between the 2-D ﬁlter and underlying image as the ﬁlter slides from left-to-right/top-to-bottom. In the Fourier domain, this convolution is equivalent to multiplying the, now 2-D, Fourier transforms of the ﬁlter and image, where the 2-D Fourier transform is given by:
∞ ∞

ωx

Figure 7.1 2-D Frequency

F [ωx , ωy ] =
k=−∞ l=−∞

f [k, l]e−i(ωxk+ωy l) .

(7.2)

The notion of low-pass, band-pass, and high-pass ﬁltering extends naturally to two-dimensional images. Shown in Figure 7.1 is a simpliﬁed decomposition of the 2-D Fourier domain parameterized by ωx and ωy ∈ [−π, π]. The inner disc corresponds to the lowest frequencies, the center annulus to the middle (band) frequencies, and the outer dark area to the highest frequencies. Two of the most common (and opposing) linear ﬁltering operations are blurring and sharpening. Both of these operations can be accomplished with a 2-D ﬁlter and 2-D convolution, or more eﬃciently with a 1-D ﬁlter and a pair of 1-D horizontal and vertical convolutions. For example, a 2-D convolution with the blur ﬁlter: 0.0625 0.1250 0.0625  0.1250 0.2500 0.1250  0.0625 0.1250 0.0625
 

47
Figure 7.2 Low-, Band-, High-pass

can be realized by convolving in the horizontal and vertical directions with the 1-D ﬁlter: blur = ( 0.25 0.50 0.25 ) . (7.3)

That is, an outer-product of the 1-D ﬁlter with itself yields the 2-D ﬁlter - the ﬁlters are xy-separable. The separability of 2-D ﬁlters is attractive for two reasons: (1) it is computationally more eﬃcient and (2) it simpliﬁes the ﬁlter design. A generic blur ﬁlter may be constructed from any row of the binomial coeﬃcients: 1 1 2 1 1

1 3 3 1 1 4 6 4 1 1 5 10 10 5 1 1 6 15 20 15 6 1 where each row (ﬁlter) should be normalized by it’s sum (i.e., blur ﬁlters should always be unit-sum so as not to increase or decrease the mean image intensity). The amount of blur is then directly proportional to the size of the ﬁlter. Blurring simply reduces the high-frequency content in an image. The opposing operation, sharpening, is meant to enhance the high-frequencies. A generic separable sharpening ﬁlter is of the form: sharp = ( 0.08 −1.00 0.08 ) . (7.4)

This ﬁlter leaves the low-frequencies intact while enhancing the contribution of the high-frequencies. Shown in Figure 7.3 are results from blurring and sharpening. 7.2 Derivative Filters
Figure 7.3 Blur Sharpen and

Discrete diﬀerentiation forms the foundation for many applications in image processing and computer vision. We are all familiar with the deﬁnition of the derivative of a continuous signal f (x): D{f (x)} = lim f (x + ε) − f (x) . ε→0 ε (7.5)

This deﬁnition requires that the signal f (x) be well deﬁned for all x ∈ R. So, does it make sense to diﬀerentiate a discretely sampled signal, D{f [x]}, which is only deﬁned over an integer sampling lattice? Strictly speaking, no. But our intuition may be that this is not such an unreasonable request. After all, we know how to diﬀerentiate f (x), from which the sampled signal f [x] was derived, so why not just diﬀerentiate the continuous signal f (x) and then sample the result? Surely this is what we have in mind when we 48

ask for the derivative of a sampled signal. But one should not be fooled by the seeming simplicity of our intuition, as we will soon discover the design of an accurate and eﬃcient discrete derivative operator will prove to be suﬃciently challenging. Recall from earlier chapters that under certain conditions (Nyquist theory), the relationship between the continuous and sampled signals can be expressed precisely as: f (x) = f [x] h(x), (7.6)

On the left-hand side of the above equation is the desired quantity, the derivative of the sampled signal. On the right-hand side is a discrete convolution between two known quantities, the sampled derivative of the sync and the original sampled signal. The derivative of the sync can be expressed analytically by simply differentiating the sync function: h (x) = π 2 x/T 2 cos(πx/T ) − π/T sin(πx/T ) , (πx/T )2 (7.10)

0

0

where T is the sampling period at which f (x) was sampled. So, if the signal f (x) is sampled above the Nyquist rate and if it is in 49

Figure 7.4 Ideal and its derivative

sync

fact diﬀerentiable, then Equation (7.9) tells us that we can exactly compute the derivative of the sampled signal f [x], an altogether happy ending. If you are feeling a bit uneasy it is for a good reason. Although mathematically correct, we have a solution for diﬀerentiating a discretely sampled signal that is physically unrealizable. In particular the derivative of the sync, h (x), is spatially inﬁnite in extent, meaning that it cannot be implemented on a ﬁnite machine. And even worse, h (x) falls oﬀ slowly from the origin so that truncation will cause signiﬁcant inaccuracies. So we are going to have to part with mathematical perfection and design a ﬁnite-length ﬁlter. To begin we need to compute the frequency response of the ideal derivative ﬁlter. We can compute the response indirectly by ﬁrst expressing f [x] in terms of its Fourier series: f [x] = 1 2π
π

Diﬀerentiation in the space domain is then seen to be equivalent to multiplying the Fourier transform F [ω] by an imaginary ramp iω. And since multiplication in the frequency domain is equivalent to convolution in the space domain, an imaginary ramp is the frequency response of the ideal derivative ﬁlter. Trying to directly design a ﬁnite length ﬁlter to this response is futile because of the discontinuity at −π/π, which of course accounts for the spatially inﬁnite extent of h (x). So we are resigned to designing a ﬁlter with a periodic frequency response that “best” approximates a ramp. The simplest such approximation is that of a sinusoid where, at least in the low-frequency range, the match is reasonably good (i.e., sin(ω) = ω, for small ω). Employing the least-squares ﬁlter design technique (Equation (4.8)) we formulate a quadratic error function to be minimized: E(h) = | M h − H |2, (7.13)

To minimize we diﬀerentiate, set equal to zero and solve for the minimal solution: h = (M tM )−1 M t H (7.14)

Since the desired frequency response, a sinusoid, has only two degrees of freedom, amplitude and phase, a 2-tap ﬁlter will suﬃce (i.e., n = 2). The resulting ﬁlter is of the form h = ( 0.5 −0.5 ). Intuitively this is exactly what we should have expected - for example, applying this ﬁlter via a convolution and evaluating at, arbitrarily, n = 0 yields: f [x] = h[x] f [x]
∞

=
k=−∞

h[x − k]f [k] (7.15)

f [0] = h[1]f [−1] + h[0]f [0] = 0.5f [0] − 0.5f [−1].

Note that the derivative is being approximated with a simple twopoint diﬀerence, that is, a discrete approximation to the continuous deﬁnition in Equation (7.5). We could of course greatly improve on this ﬁlter design. But since we are really interested in multi-dimensional diﬀerentiation, let’s put aside further analysis of the one-dimensional case and move on to the two-dimensional case. It has been the tendency to blindly extend the one-dimensional design to higher-dimensions, but, as we will see shortly, in higherdimensions the story becomes slightly more complicated. In the context of higher-dimensional signals we ﬁrst need to consider partial derivatives. For example the partial derivative of a two dimensional signal f (x, y) in it’s ﬁrst argument is deﬁned as: fx (x, y) ≡ ∂f (x, y) ∂x f (x + ε, y) − f (x, y) = lim . ε→0 ε

(7.16)

According to the Nyquist theory, the continuous and discrete signals (if properly sampled) are related by the following equality: f (x, y) = f [x, y] h(x, y), (7.17)

Notice that calculating the partial derivative requires a pair of one-dimensional convolutions: a derivative ﬁlter, h [x], in the dimension of diﬀerentiation, and an interpolation ﬁlter, h[y], in the other dimension (for multi-dimensional signals, all remaining dimensions would be convolved with the interpolation ﬁlter). Since two-dimensional diﬀerentiation reduces to a pair of onedimensional convolutions it is tempting to simply employ the same diﬀerentiation ﬁlter used in the one-dimensional case. But since a pair of ﬁlters are now required perhaps we should give this some additional thought. In some ways the choice of ﬁlters seems trivial: chose an interpolation function h(x), diﬀerentiate it to get the derivative function h (x), and sample these functions to get the ﬁnal digital ﬁlters h[x] and h [x]. So how is this diﬀerent from the one-dimensional case? In the one-dimensional case only the derivative ﬁlter is employed, whereas in the two-dimensional case we require the pair of ﬁlters. And by our formulation we know that the pair of ﬁlters should satisfy the relationship that one is the derivative of the other h (x) = D(h(x)). And in fact this constraint is automatically enforced by the very nature in which the continuous functions are chosen, but in the ﬁnal step, these functions are sampled to produce discrete ﬁlters. This sampling step typically destroys the required derivative relationship, and, although a seemingly subtle point, has dramatic eﬀects on the accuracy of the resulting derivative operator. For example √ consider the often √ used Sobel derivative ﬁlters with h[x] = ( 1 2 1 ) /(2 + 2) and h [x] = ( 1 0 −1 ) /3. Shown in Figure 7.7 are the magnitudes of the Fourier transform of the derivative ﬁlter (solid line) and the interpolation ﬁlter times iω (i.e., frequency domain diﬀerentiation). If the ﬁlters obeyed the required derivative relationship than these curves would be exactly matched, which they clearly 52

f[y]
Figure 7.6 Horizontal partial diﬀerentiation

−pi

0

pi

Figure 7.7 Sobel frequency response

are not. The mismatching of the ﬁlters results in gross inaccuracies in derivative measurements. Let’s see then if we can design a better set of ﬁlters. We begin by writing down the desired relationship between the derivative and interpolation ﬁlters, most conveniently expressed in the frequency domain: H (ω) = iωH(ω), (7.21)

where the columns of the matrix Fm×n contain the ﬁrst n Fourier basis functions (i.e., a discrete-time Fourier transform), the matrix Fm×n = iωFm×n , and Wm×m is a diagonal frequency weighting matrix. Note that the dimension n is determined by the ﬁlter size and the dimension m is the sampling rate of the continuous Fourier basis functions, which should be chosen to be suﬃciently large to avoid sampling artifacts. This error function can be expressed more concisely as: E(u) = |M u|2 , (7.24)

The minimal unit vector u is then simply the minimal-eigenvalue eigenvector of the matrix M tM . After solving for u, the derivative and interpolation ﬁlters can be “unpacked” and normalized so that the interpolation ﬁlter is unit sum. Below are the resulting ﬁlter values for a 3-tap and 5-tap ﬁlter pair.

53

Shown in Figure 7.8 are the matching of these ﬁlters in the frequency domain. Notice that the 5-tap ﬁlters are nearly perfectly matched. h h h h
−pi 0 pi

Higher-order derivative ﬁlters can be designed by replacing the initial constraint in Equation (7.21) with H (ω) = (iω)k H(ω) for a kth order derivative. A peculiar aspect of this ﬁlter design is that nowhere did we explicitly try to model a speciﬁed frequency response. Rather, the design fell naturally from the relationship between the continuousand discrete-time signals and the application of the continuous derivative operator, and in this way is quite distinct from the onedimensional case. The proper choice of derivative ﬁlters can have a dramatic impact on the applications which utilize them. For example, a common application of diﬀerential measurements is in measuring motion from a movie sequence f (x, y, t). The standard formulation for motion estimation is: fx (x, y, t)vx(x, y) + fy (x, y, t)vy (x, y) + ft (x, y, t) = 0, (7.26)

−pi

0

pi

Figure 7.8 Frequency response of matched derivative ﬁlters

where the motion vector is v = ( vx vy )t , and fx (·), fy (·), and ft (·) are the partial derivatives with respect to space and time. Shown in Figure 7.9 are the resulting motion ﬁelds for a simple translational motion with the Sobel (top panel) and matched (bottom) derivative ﬁlters used to compute the various derivatives. Although these ﬁlters are the same size, the diﬀerence in accuracy is signiﬁcant.
Figure 7.9 Diﬀerential motion estimation

7.3 Steerable Filters In the previous section we showed how to compute horizontal and vertical partial derivatives of images. One may naturally wonder how to compute a derivative in an arbitrary direction. Quite remarkably it turns out that we need not design a new set of ﬁlters for each possible direction because the derivative in any direction can be synthesized from a linear combination of the horizontal and vertical derivatives. This property of derivatives has been termed steerability. There are several formulations of this property, we chose to work in the frequency domain where diﬀerentiation takes on a particularly simple form. 54

To begin, we express a two-dimensional image with respect to its Fourier series:
π π

Now, diﬀerentiation in the vertical direction v is equivalent to multiplying the Fourier transform by an imaginary vertical ramp. This trend generalizes to arbitrary directions, that is, the partial derivative in any direction α can be computed by multiplying the Fourier transform by an imaginary oriented ramp −jωα :
π π

fα (x, y) =
ωx =−π ωy =−π

−jωα F (ωx , ωy )e−j(ωx x+ωy y) ,(7.30)

where the oriented ramp can be expressed in terms of the horizontal and vertical ramps: ωα = cos(α)ωx + sin(α)ωy . (7.31)

Substituting this deﬁnition back into the partial derivative in α, Equation (7.30), gives:
π π

fα (x, y) =
ωx =−π ωy =−π π

−j[cos(α)ωx + sin(α)ωy ]F (ωx , ωy )e−j(ωx x+ωy y)
π

= cos(α)
ωx =−π ωy =−π π

−jωx F (ωx , ωy )e−j(ωx x+ωy y)
π

+ sin(α)
ωx =−π ωy =−π

−jωy F (ωx , ωy )e−j(ωx x+ωy y) (7.32)

= cos(α)fx (x, y) + sin(α)fy (x, y).
Recall that the derivative of an exponential is an exponential, so that according to the chain rule, Dx {eax } = aeax .
8

55

Notice that we obtain the horizontal and vertical derivatives when α = 0 and α = 90. This equation embodies the principle of steerability - the derivative in any direction α can be synthesized from a linear combination of the partial horizontal and vertical derivatives, fx (x, y) and fy (x, y). Pause to appreciate how remarkable this is, a pair of directional derivatives is suﬃcient to represent an inﬁnite number of other directional derivatives, i.e., α can take on any real-valued number. From the previous section we know how to compute the horizontal and vertical derivatives via convolutions with an interpolation and derivative ﬁlter. To compute any other directional derivative no more convolutions are required, simply take the appropriate linear combinations of the horizontal and vertical derivatives as speciﬁed in Equation (7.32). Shown in Figure 7.10 from top to bottom is a disc f (x, y), its horizontal derivative fx (x, y), its vertical derivative fy (x, y), and its steered derivative f45 (x, y), where the steered derivative was synthesized from the appropriate linear combinations of the horizontal and vertical derivatives. The obvious beneﬁt of steerability is that the derivative in any direction can be synthesized with minimal computational costs. Steerability is not limited to ﬁrst-order derivatives. Higher-order derivatives are also steerable; the N th-order derivative is steerable with a basis set of size N + 1. For example, the second-order derivative in an arbitrary direction can be synthesized as follows: fαα = cos2 (α)fxx + 2 cos(α) sin(α)fxy + sin2(α)fyy , (7.33)

where for notational simplicity the spatial arguments (x, y) have been dropped, and the multiple subscripts denote higher-order diﬀerentiation. Note that three partial derivatives are now needed to steer the second-order derivative. Similarly, the third-order derivative can be steered with a basis of size four: fααα = cos3 (α)fxxx + 3 cos2 (α) sin(α)fxxy +3 cos(α) sin2 (α)fxyy + sin3 (α)fyyy . (7.34) You may have noticed that the coeﬃcients needed to steer the basis set look familiar, they are the binomial coeﬃcients that come from a polynomial expansion. More speciﬁcally, as in Equation(7.30) the N th-order derivative in the frequency domain is computed by multiplying the Fourier transform by an imaginary oriented ramp raised to the N th power, (−jωα )N . Expressing this oriented ramp in terms of the horizontal and vertical ramps provides the basis and coeﬃcients needed to steer derivatives of arbitrary order: (ωα)N = (cos(α)ωx + sin(α)ωy )N . 56 (7.35)

Figure 7.10 Steerability

Although presented in the context of derivatives, the principle of steerability is not limited to derivatives. In the most general case, a two-dimensional ﬁlter f (x, y) is steerable in orientation if it can be expressed as a polar-separable function, g(r)h(θ), where h(θ) is band-limited. More speciﬁcally, for an arbitrary radial component g(r), and for h(θ) expressed as:
N

h(θ) =
n=1

an cos(nθ) + bn sin(nθ)

(7.36)

then the ﬁlter is steerable with a basis size of 2N . 7.4 Edge Detection Discrete diﬀerentiation forms the foundation for many applications in computer vision. One such example is edge detection - a topic that has received an excessive amount of attention, but is only brieﬂy touched upon here. An edge is loosely deﬁned as an extended region in the image that undergoes a rapid directional change in intensity. Diﬀerential techniques are the obvious choice for measuring such changes. A basic edge detector begins by computing ﬁrst-order spatial derivatives of an image f [x, y]: fx [x, y] = (f [x, y] h [x]) h[y] fy [x, y] = (f [x, y] h[x]) h [y], (7.37) (7.38)

where h [·] and h[·] are the derivative and preﬁlter deﬁned in Section 7.2. The “strength” of an edge at each spatial location is deﬁned to be the magnitude of the gradient vector [x, y] = ( fx [x, y] fy [x, y] ), deﬁned as:
Figure 7.11 Edges

|

[x, y]| =

2 2 fx [x, y] + fy [x, y].

(7.39)

As shown in Figure 7.11, the gradient magnitude is only the beginning of a more involved process (not discussed here) of extracting and localizing the salient and relevant edges. 7.5 Wiener Filter For any of a number of reasons a digital signal may become corrupted with noise. The introduction of noise into a signal is often modeled as an additive process, s = s+n. The goal of de-noising is ˆ to recover the original signal s from the corrupted signal s. Given ˆ a single constraint in two unknowns this problem is equivalent to my asking you “37 is the sum of two numbers, what are they?” Lacking clairvoyant powers or knowledge of how the individual numbers were selected we have little hope of a solution. But by 57
n

s

+

s

Figure 7.12 Additive noise

making assumptions regarding the signal and noise characteristics and limiting ourselves to a linear approach, a solution can be formulated known as the Wiener ﬁlter, of famed Mathematician Norbert Wiener (1894-1964). Having restricting ourselves to a linear solution, our goal is to design a ﬁlter h[x] such that: s[x] = h[x] s[x] ˆ = h[x] (s[x] + n[x]),
0 −pi

1

(7.40)

0

pi

that is, when the ﬁlter is convolved with the corrupted signal the original signal is recovered. With this as our goal, we reformulate this constraint in the frequency domain and construct a quadratic error functional to be minimized: E(H(ω)) = dω [H(ω)(S(ω) + N (ω)) − S(ω)]2. (7.41)

At an intuitive level this frequency response makes sense - when the signal is strong and the noise is weak the response is close to 1 (i.e., frequencies are passed), and when the signal is weak and the noise is strong the response is close to 0 (i.e., frequencies are stopped). So we now have an optimal (in the least-squares sense) 58

frequency response in terms of the signal and noise characteristics, but of course we don’t typically know what those are. But we can instantiate them by making assumptions about the general statistical nature of the signal and noise, for example a common choice is to assume white noise, N (ω) is constant for all ω, and, for natural images, to assume that S(ω) = 1/ω p. The frequency response in the top panel of Figure 7.13 was constructed under these assumptions. Shown in the bottom panel is a 7-tap ﬁlter derived from a least-squares design. This one-dimensional formulation can easily be extended to two or more dimensions. Shown in Figure 7.14 from top to bottom, is Einstein, Einstein plus noise, and the results of applying a 7 × 7 Wiener ﬁlter. Note that the noise levels are reduced but that much of the sharp image structure has also been lost, which is an unfortunate but expected side eﬀect given that the Wiener ﬁlter is low-pass in nature.

59

8. Non-Linear Filtering

8-1 Median Filter 8-2 Dithering

8.1 Median Filter Noise may be introduced into an image in a number of diﬀerent ways. In the previous section we talked about how to remove noise that has been introduced in an additive fashion. Here we look at a diﬀerent noise model, one where a small number of pixels are corrupted due to, for example, a faulty transmission line. The corrupted pixels randomly take on a value of white or black, hence the name salt and pepper used to describe such noise patterns (Figure 8.1). Shown in the middle panel of Figure 8.1 is the disastrous result of applying the solution from the additive noise model (Wiener ﬁlter) to the salt and pepper noise image in the top panel. Trying to average out the noise in this fashion is equivalent to asking for the average salary of a group of eight graduate students and Bill Gates. As the income of Gates will skew the average salary so does each noise pixel when its value is so disparate from its neighbors. In such cases, the mean is best replaced with the median, computed by sorting the set of numbers and reporting on the value that lies midway. Shown in the bottom panel of Figure 8.1 is the much improved result of applying a 3 × 3 median ﬁlter to the salt and pepper noise image. More speciﬁcally, the center pixel of each 3 × 3 neighborhood is replaced with the median of the nine pixel values in that neighborhood. Depending on the density of the noise the median ﬁlter may need to be computed over a larger neighborhood. The tradeoﬀ being that a larger neighborhood leads to a loss of detail, however this loss of detail is quite distinct from that of averaging. For example, shown in Figure 8.2 is the result of applying a 15×15 median ﬁlter. Notice that although many of the internal details have been lost the boundary contours (edges) have been retained, this is often referred to as posterization. This eﬀect could never be achieved with an averaging ﬁlter which would indiscriminately smooth over all image structures. Because of the non-linear sorting step of a median ﬁlter it cannot be implemented via a simple convolution and is thus often more costly to implement. Outside the scope of this presentation there are a number of tricks for reducing the computational demands of a median ﬁlter.

Salt & Pepper Noise

Wiener

Median
Figure 8.1 Median ﬁlter

Figure 8.2 15 × 15 median ﬁlter

60

8.2 Dithering Dithering is a process by which a digital image with a ﬁnite number of gray levels is made to appear as a continuous-tone image. For example, shown at the top of Figure 8.3 is an 8-bit (i.e., 256 gray values) grayscale image of Richard P. Feynman. Shown below are two 1 bit (i.e., 2 gray values) images. The ﬁrst was produced by thresholding, and the second by dithering. Although in both images each pixel takes on only one of two gray values (black or white), it is clear that the ﬁnal image quality is critically dependent on the way in which pixels are quantized. There are numerous dithering algorithms (and even more variants within each algorithm). Sadly there are few quantitative metrics for measuring the performance of these algorithms. A standard and reasonably eﬀective algorithm is a stochastic error diﬀusion algorithm based on the Floyd/Steinberg algorithm. The basic Floyd/Steinberg error diﬀusion dithering algorithm tries to exploit local image structure to reduce the eﬀects of quantization. For simplicity, a 1-bit version of this algorithm is described here, the algorithm extends naturally to an arbitrary number of gray levels. This algorithm operates by scanning through the image, left to right and top to bottom. At each pixel, the gray value is ﬁrst thresholded into “black” or “white”, the diﬀerence between the new pixel value and the original value is then computed and distributed in a weighted fashion to its neighbors. Typically, the error is distributed to four neighbors with the following weighting:
1 16

×

• 7 3 5 1

,

where the • represents the thresholded pixel, and the position of the weights represent spatial position on a rectangular sampling lattice. Since this algorithm makes only a single pass through the image, the neighbors receiving a portion of the error must consist only of those pixels not already visited (i.e., the algorithm is casual). Note also that since the weights have unit sum, the error is neither ampliﬁed nor reduced. As an example, consider the quantization of an 8-bit image to a 1-bit image. An 8-bit image has 255 gray values, so all pixels less than 128 (mid-level gray) are thresholded to 0 (black), and all values greater than 128 are thresholded to 255 (white). A pixel at position (x, y) with intensity value 120 is thresholded to 0. The error, 120-0 = 120, is distributed to four neighbors as follows: (7/16)120 is added to the pixel at position (x + 1, y), (3/16)120 is added to the pixel at 61

Figure 8.3 Thresholding and Dithering

position (x − 1, y + 1), (5/16)120 to pixel (x, y + 1) and (1/16)120 to pixel (x + 1, y + 1). The intuition behind this algorithm is that when the pixel value of 120 is thresholded to 0 that pixel is made darker. By propagating this error the surrounding pixels are brightened making it more likely that they will be thresholded to something brighter. As a result the local neighborhood maintains it’s average gray value. Qualitatively, the error diﬀusion algorithm reduces the eﬀects of quantization via simple thresholding. However this algorithm does introduce correlated artifacts due to the deterministic nature of the algorithm and scanning order. These problems may be partially alleviated by introducing a stochastic process in a variety of places. Two possibilities include randomizing the error before distributing it to its neighbors (e.g., randomly scaling the error in the range 0.9 to 1.1), and alternating the scanning direction (e.g., odd lines are scanned left to right, and even lines from right to left).

62

9. Multi-Scale Transforms

63

10. Motion Estimation

10-1 Diﬀerential Motion 10-2 Diﬀerential Stereo

10.1 Diﬀerential Motion Our visual world is inherently dynamic. People, cars, dogs, etc. are (usually) moving. These may be gross motions, walking across the room, or smaller motions, scratching behind your ear. Our task is to estimate such image motions from two or more images taken at diﬀerent instances in time. With respect to notation, an image is denoted as f (x, y) and an image sequence is denoted as f (x(t), y(t), t), where x(t) and y(t) are the spatial parameters and t is the temporal parameter. For example, a sequence of N images taken in rapid succession may be represented as f (x(t), y(t), t + i∆t) with i ∈ [0, N − 1], and ∆t representing the amount of time between image capture (typically on the order of 1/30th of a second). Given such an image sequence, our task is to estimate the amount of motion at each point in the image. For a given instant in space and time, we require an estimate of motion (velocity) v = ( vx vy ), where vx and vy denote the horizontal and vertical components of the velocity vector v. Shown in Figure 10.1 are a pair of images taken at two moments in time as a textured square is translating uniformly across the image. Also shown is the corresponding estimate of motion often referred to as a ﬂow ﬁeld. The ﬂow ﬁeld consists of a velocity vector at each point in the image (shown of course are only a subset of these vectors). In order to estimate motion, an assumption of brightness constancy is made. That is, it is assumed that as a small surface patch is moving, its brightness value remains unchanged. This constraint can be expressed with the following partial diﬀerential equation: ∂f (x(t), y(t), t) = 0. ∂t (10.1)

where the partials of the spatial parameters x and y with respect to time correspond to the velocity components: fx vx + fy vy + ft = 0. 64 (10.3)

The subscripts on the function f denote partial derivatives. Note again that this constraint holds for each point in space and time but that for notational simplicity the spatial/temporal parameters are dropped. This transformed brightness constancy constraint is rewritten by packing together the partial derivatives and velocity components into row and column vectors. ( fx fy ) vx vy + ft = 0 (10.4)

t fs v + ft = 0.

The space/time derivatives fs and ft are measured quantities, leaving us with a single constraint in two unknowns (the two components of the velocity vector, v). The constraint can be solved by assuming that the motion is locally similar, and integrating this constraint over a local image neighborhood. A least-squares error function takes the form:
2

If the matrix M is invertible (full rank), then the velocity can be estimated by simply left multiplying by the inverse matrix: v = −M −1 b (10.8)

The critical question then is, when is the matrix M invertible? Generally speaking the matrix is rank deﬁcient, and hence not invertible, when the intensity variation in a local image neighborhood varies only one-dimensionally (e.g., fx = 0 or fy = 0) or zero-dimensionally (fx = 0 and fy = 0). These singularities are sometimes referred to as the aperture and blank wall problem. The motion at such points simply can not be estimated. 65

Motion estimation then reduces to computing, for each point in space and time, the spatial/temporal derivatives fx , fy , and ft . Of course the temporal derivative requires a minimum of two images, and is typically estimated from between two and seven images. The spatial/temporal derivatives are computed as follows. Given a temporal sequence of N images, the spatial derivatives are computed by ﬁrst creating a temporally preﬁltered image. The spatial derivative in the horizontal direction fx is estimated by preﬁltering this image in the vertical y direction and diﬀerentiating in x. Similarly, the spatial derivative in the vertical direction fy is estimated by preﬁltering in the horizontal x direction and diﬀerentiating in y. Finally, the temporal derivative is estimated by temporally differentiating the original N images, and preﬁltering the result in both the x and y directions. The choice of ﬁlters depends on the image sequence length: an N tap pre/derivative ﬁlter pair is used for an image sequence of length N (See Section 7). 10.2 Diﬀerential Stereo Motion estimation involves determining, from a single stationary camera, how much an object moves over time (its velocity). Stereo estimation involves determining the displacement disparity of a stationary object as it is imaged onto a pair of spatially oﬀset cameras. As illustrated in Figure 10.2, these problems are virtually identical: velocity (v) ≡ disparity (∆). Motion and stereo estimation are often considered as separate problems. Motion is thought of in a continuous (diﬀerential) framework, while stereo, with its discrete pair of images, is thought of in terms of a discrete matching problem. This dichotomy is unnecessary: stereo estimation can be cast within a diﬀerential framework. Stereo estimation typically involves a pair of cameras spatially oﬀset in the horizontal direction such that their optical axis remain parallel (Figure 10.2). Denoting an image as f (x, y), the image that is formed by translating the camera in a purely horizontal direction is given by f (x + ∆(x, y), y). If a point in the world (X, Y, Z) is imaged to the image position (x, y), then the shift ∆(x, y) is inversely proportional to the distance Z (i.e., nearby objects have large disparities, relative to distant objects). Given this, a stereo pair of images is denoted as: fL (x + δ(x, y), y) and fR (x − δ(x, y), y), (10.9)

V

D

Figure 10.2 Motion and Stereo

where the disparity ∆ = 2δ. Our task is to determine, for each point in the image, the disparity (δ) between the left and right images. That is, to ﬁnd the shift that brings the stereo pair back 66

where for notational convenience, the spatial parameters have been dropped. Diﬀerentiating, setting the result equal to zero and solving for δ yields: dE(δ) = 2(fL + fR )x [(fL − fR ) + δ(fL + fR )x] dδ = 0 fL − fR δ = − (10.13) (fL + fR )x Stereo estimation then reduces to computing, for each point in the image, spatial derivatives and the diﬀerence between the left and right stereo pair (a crude derivative with respect to viewpoint). Why, if motion and stereo estimation are similar, do the mathematical formulations look so diﬀerent? Upon closer inspection they are in fact quite similar. The above formulation amounts to a constrained version of motion estimation. In particular, because of the strictly horizontal shift of the camera pair, the disparity was constrained along the horizontal direction. If we reconsider the motion estimation formulation assuming motion only along the horizontal direction, then the similarity of the formulations becomes evident. Recall that in motion estimation the brightness constancy assumption led to the following constraint: ∂f ∂x ∂f ∂y ∂f + + ∂x ∂t ∂y ∂t ∂t = 0, (10.14)

where the partial derivative of the spatial parameter x with respect to time correspond to the motion (speed) in the horizontal direction: fx vx + ft = 0. (10.16)

Unlike before, this leads to a single constraint with a single unknown which can be solved for directly: vx = − ft . fx (10.17)

This solution now looks very similar to the solution for diﬀerential stereo in Equation 10.13. In both solutions the numerator is a derivative, in one case with respect to time (motion) and in the other with respect to viewpoint (stereo). Also in both solutions, the denominator is a spatial derivative. In the stereo case, the denominator consists of the spatial derivative of the sum of the left and right image pair. This may seem odd, but recall that diﬀerentiation of a multi-dimensional function requires diﬀerentiating along the desired dimension and preﬁltering along all other dimensions (in this case the viewpoint dimension). In both the diﬀerential motion and stereo formulations there exists singularities when the denominator (spatial derivative) is zero. As with the earlier motion estimation this can be partially alleviated by integrating the disparities over a local image neighborhood. However, if the spatial derivative is zero over a large area, corresponding to a surface in the world with no texture, then disparities at these points simply can not be estimated.

68

11. Useful Tools

11.1 Expectation/Maximization The Expectation/Maximization (EM) algorithm simultaneously segments and ﬁts data generated from multiple parametric models. For example, shown in Figure 11.1 are a collection of data points (x, y) generated from one of two linear models of the form: y(i) = a1 x(i) + b1 + n1 (i) or y(i) = a2 x(i) + b2 + n2 (i),(11.1) where the model parameters are a1 , b1 and a2 , b2, and the system is modeled with additive noise n1 (i) and n2 (i). If we are told the model parameters, then determining which data point was generated by which model would be a simple matter of choosing, for each data point i, the model k that minimizes the error between the data and the model prediction: rk (i) = |ak x(i) + bk − y(i))|, (11.2)

for k = 1, 2 in our current example. On the other hand, if we are told which data points were generated by which model, then estimating the model parameters reduces to solving, for each model k, an over-constrained set of linear equations: xk (1)  xk (2)   .  . .


Figure 11.1 Data two models

from

xk (n) 1

1 1  . . .


ak bk

yk (1)  yk (2)    =  . ,  .  .
 

(11.3)

yk (n)

where the xk (i) and yk (i) all belong to model k. In either case, knowing one piece of information (the model assignment or parameters) makes determining the other relatively easy. But, lacking either piece of information makes this a considerably more diﬃcult estimation problem. The EM algorithm is an iterative two step algorithm that estimates both the model assignment and parameters. The “E-step” of EM assumes that the model parameters are known (initially, the model parameters can be assigned random values) and calculates the likelihood of each data point belonging to each model. In so doing the model assignment is made in a “soft” probabilistic fashion. That is, each data point is not explicitly assigned a single model, instead each data point i is assigned a 69

where, σN is proportional to the amount of noise in the data, and for each data point i, k wk (i) = 1. The “M-step” of EM takes the likelihood of each data point belonging to each model, and re-estimates the model parameters using weighted least-squares. That is, the following weighted error function on the model parameters is minimized: Ek (ak , bk ) =
i

wk (i)[ak x(i) + bk − y(i)]2.

(11.7)

The intuition here is that each data point contributes to the estimation of each model’s parameters in proportion to the belief that it belongs to that particular model. This quadratic error function is minimized by computing the partial derivatives with respect to the model parameters, setting the result equal to zero and solving for the model parameters. Diﬀerentiating: ∂Ek (ak , bk ) ∂ak ∂Ek (ak , bk ) ∂bk =
i

2wk (i)x(i)[akx(i) + bk − y(i)] 2wk (i)[ak x(i) + bk − y(i)], (11.8)

=
i

setting both equal to zero yields the following set of linear equations: ak
i

wk (i)x(i)2 + bk
i

wk (i)x(i) =
i

wk (i)x(i)y(i)(11.9) wk (i)y(i). (11.10)
i

ak
i

wk (i)x(i) + bk
i

wk (i) = 70

Rewriting in matrix form:
i

wk (i)x(i)2 i wk (i)x(i)

i wk (i)x(i) i wk (i)

ak bk

=

i wk (i)x(i)y(i) i

w(i)y(i) (11.11)

Axk = b xk = A−1 b,

yields a weighted least squares solution for the model parameters. Note that this solution is identical to solving the set of linear equations in Equation (11.3) using weighted least-squares. The EM algorithm iteratively executes the “E” and “M” step, repeatedly estimating and reﬁning the model assignments and parameters. Shown in Figure 11.2 are several iterations of EM applied to ﬁtting data generated from two linear models. Initially, the model parameters are randomly assigned, and after six iterations, the algorithm converges to a solution. Beyond the current scope are proofs on the convergence and rate of convergence of EM. At the end of this chapter is a Matlab implementation of EM for ﬁtting multiple linear models to two-dimensional data. 11.2 Principal Component Analysis 11.3 Independent Component Analysis