Tài liệu

# Fundamentals of Image Processing

Science and Technology

# Image Representation

## Continuous-domain, discrete-domain, and finite-size images

An image is a spatially varying signal $s(x,y)$ where $x$ and $y$ are two spatial coordinates. The signal value $s(x,y)$ at each spatial location $(x,y)$ can be either a scalar (e.g. light intensity for gray scale images) or a vector (e.g. 3 dimensional vector for RGB color images, or more general $P$-dimensional vector for multispectral images). In the latter case, we could treat each vector component separately as a scalar image (referred to as a channel), sometimes after a certain transformation in the $P$-dimensional vector space.

In digital image processing, images are discretized into samples at discrete spatial locations that are indexed by integer coordinates $[m,n]$. Typically, a discrete-domain image $s[m,n]$ is related to a continuous-domain image $s(x,y)$ through the sampling operation

$s [ m , n ] = s ( m Δ x , n Δ y ) ,$

where $Δx$ and $Δy$ are sampling intervals in $x$ and $y$ dimensions, respectively. More general, a discrete-space image is obtained through the generalized-sampling operation

$w [ m , n ] = ∫ - ∞ ∞ ∫ - ∞ ∞ s ( x , y ) φ m , n ( x , y ) d x d y ,$

where $φm,n(x,y)$ is the point-spread function of the image sensor (e.g. a photometric sensor in a digital camera) at the location indexed by $(m,n)$. Typically, point-spread functions at different locations are simply shifted versions of a single function as

$φ m , n ( x , y ) = φ ( x - m Δ x , y - n Δ y ) ,$

and $φm,n(x,y)$ is called the sampling kernel.

Furthermore, a discrete image $s[m,n]$ is often of finite size; for example $0≤m≤M-1,0≤n≤N-1$. Then $s[m,n]$ can also be treated as an $M×N$ matrix. The image sample $s[m,n]$ and the corresponding location $[m,n]$ is often called a pixel, or picture element.

## Fourier transforms and sampling theorem

It is often very effective, conceptually and computationally, to represent images in the frequency domain using the Fourier transform. For a continuous-domain image $s(x,y)$, its Fourier transform is defined as

$S ( u , v ) = ∫ - ∞ ∞ ∫ - ∞ ∞ s ( x , y ) e - j 2 π ( x u + y v ) d x d y .$

Here, $u$ and $v$ denote frequency variables and they have reciprocal unit with $x$ and $y$. For example, if the spatial coordinate $x$ has unit in $mm$, then the corresponding frequency variable $u$ has unit in $mm-1$. Under certain conditions, the image $s(x,y)$ can be exactly recovered from its frequency-domain $S(u,v)$ by the inverse Fourier transform

$s ( x , y ) = ∫ - ∞ ∞ ∫ - ∞ ∞ S ( u , v ) e j 2 π ( x u + y v ) d u d v .$

We denote this pair of signals related by the Fourier transform (FT) as

$s ( x , y ) ⟷ FT S ( u , v ) .$

For a discrete image $s[m,n]$ the discrete-space Fourier transform (DSFT) relation

$s [ m , n ] ⟷ DSFT S ( u , v ) .$

is defined as

$S d ( u , v ) = ∑ m = - ∞ ∞ ∑ n = - ∞ ∞ s [ m , n ] e - j 2 π ( m u + n v ) , s [ m , n ] = ∫ - 1 / 2 1 / 2 ∫ - 1 / 2 1 / 2 S d ( u , v ) e j 2 π ( m u + n v ) d u d v .$

It is easy to see that $Sd(u,v)$ is a periodic function

$S d ( u + k , v + l ) = S d ( u , v ) , for all k , l ∈ Z ,$

and thus we only need to consider the function in one period; e.g. $Sd(u,v)$ with $|u|≤1/2,|v|≤1/2$.

Theorem 1 (Sampling) Suppose that the discrete-domain image $s[m,n]$ is related to the continuous-domain image $s(x,y)$ through the sampling operation [link]. Then their Fourier transforms are related by

$S d ( u , v ) = 1 Δ x Δ y ∑ k ∈ Z ∑ l ∈ Z S u + k Δ x , v + l Δ y .$

(Sketch) One way to prove this is to express $s[m,n]$ using [link] by substituting $x=mΔx,y=nΔy$ and then “match” with the right-hand side of [link].

The summation on the right-hand side of [link] consists of $S(u/Δx,v/Δy)$ and its translated copies in frequency by $(k,l)$. These copies with $(k,l)≠(0,0)$ are called alias terms. If $s(x,y)$ is bandlimited such that

$S ( u , v ) = 0 , for | u | ≥ 1 / ( 2 Δ x ) , | v | ≥ 1 / ( 2 Δ v ) ,$

then these alias terms do not overlap with $S(u/Δx,v/Δy)$, and thus $S(u,v)$ can be exactly recovered from $Sd(u,v)$ simply by

$S ( u , v ) = Δ x Δ y rect ( Δ x u ) rect ( Δ y v ) S d ( Δ x u , Δ y v ) .$

Here the rectangular function is defined as

$rect ( x ) = 1 if | x | ≤ 1 / 2 0 else.$

We will show later that [link] in the spatial domain is equivalent to

$s ( x , y ) = ∑ k ∈ Z ∑ l ∈ Z s [ m , n ] sinc ( t / Δ x - m ) sinc ( t / Δ y - n ) ,$

where the sinc function is defined as

$sinc ( x ) = sin ( π x ) π x .$

For the discrete image $s[m,n]$ of finite size $M×N$ with $0≤m≤M-1,0≤n≤N-1$, we have the discrete Fourier transform (DFT) relation

$s [ m , n ] ⟷ DFT S [ k , l ] ,$

which is defined as

$S [ k , l ] = ∑ m = 0 M - 1 ∑ n = 0 N - 1 s [ m , n ] e - j 2 π ( m k / M + n l / N ) , s [ m , n ] = 1 M N ∑ k = 0 M - 1 ∑ l = 0 N - 1 S [ k , l ] e j 2 π ( m k / M + n l / N ) .$

Therefore the DFT maps an $M×N$ image in the spatial domain into an $M×N$ image in the frequency domain; both images can have complex values.

Relating [link] to [link], we see that if the $M×N$ image $s[m,n]$ is zero padded outside its support $[0,M-1]×[0,N-1]$ then $S[k,l]$ is a sampled image of $Sd(u,v)$,

$S [ k , l ] = S d ( k / M , l / N ) .$

In summary, we have seen the following three Fourier transforms

$continuous-domain ⟷ FT continuous-domain discrete-domain ⟷ DSFT continuous-domain discrete-domain ⟷ DFT discrete-domain$

Among these transforms, only the last one, the DFT, is computationally feasible (i.e. with summations of finite terms). Moreover, the DFT can be implemented efficiently with fast Fourier transform algorithms. In moving from the FT to the DSFT and then to the DFT, we first discretize the spatial domain and then the frequency domain. Therefore, it is important to understand [link] and [link] so that we can relate the computational results and images by the DFT to the frequency representation of the original image in the real world.

## Vector-space framework

In a more abstract framework, we can view each image as a vector in an appropriate vector space (i.e. for continuous-domain, discrete-domain, or discrete-domain of finite support). The associate Fourier transform is a linear mapping or linear operator that maps a vector in the spatial domain into a vector in the frequency domain. We can express the Fourier transform and its inverse using the matrix-vector multiplication notation

$Forward or Analysis: S = F s , Inverse or Synthesis: s = F - 1 S ,$

In the vector-space framework, it is particularly useful to view the inverse Fourier transform as a basis expansion. For example, the inverse DFT in [link] can be written as

$s = ∑ k = 0 M - 1 ∑ l = 0 N - 1 S [ k , l ] f k , l ,$

in which the image $s$ is expanded as a linear combination of basis images$fk,l$ where

$f k , l [ m , n ] = 1 M N e j 2 π ( m k / M + n l / N ) .$

Other transforms such as the wavelet transform simply provide other basis expansions.

Problems

1. Complete the proof of the sampling theorem.
2. Suppose that $s(x,y)$ is a bandlimited image with its Fourier transform $S(fx,fy)=0$ for $|fx|≥1/(2Δx),|fy|≥1/(2Δy)$. This image is captured by an CCD array where sensors are placed on a rectangular grid with spacing $Δx×Δy$ and each sensor measures the integral of light intensity falling in an area of size $Δx×Δy$ on this grid. Derive a formula to recover $s(x,y)$ exactly from these discrete CCD measurements.
3. Prove that the inverse DFT perfectly recovers the image. That is, suppose that $S[k,l]$ is given in [link] then show that the right-hand side of [link] indeed returns $s[m,n]$.
4. An 3 in by 4 in photo is discretized by a 200 dpi (dots-per-inch) scanner and results in a 600 by 800 digital image. Assume that alias is negligible. Suppose that we want to filter out all spatial high frequency above 50 $inch-1$ in both horizontal and vertical dimensions out of the image. One way is to zero-pad the image to the size of 1024 by 1024 so that we can use a fast Fourier transform (FFT) algorithm to compute the DFT of size 1024 by 1024 of the zero-padded image. Find the indexes of the DFT coefficients that need to be zero out before taking the inverse DFT to achieve the desired filtering effect.

# Image Filtering

## Convolution operations

The image (linear) filtering operation in the continuous-domain is defined by convolution

$r ( x , y ) = ( s * h ) ( x , y ) = ∫ - ∞ ∞ ∫ - ∞ ∞ h ( x ' , y ' ) s ( x - x ' , y - y ' ) d x ' d y ' .$

And similarly in the discrete-domain:

$r [ m , n ] = ( s * h ) [ m , n ] = ∑ m ' = - ∞ ∞ ∑ n ' = - ∞ ∞ h [ m ' , n ' ] s [ m - m ' , n - n ' ] .$

The two-dimensional signal $h(x,y)$ or $h[m,n]$ is called filter, mask, or point-spread function.

## Examples

Example 1 (First-order derivatives) The first-order derivatives in the $x$ and $y$ directions of a discrete image $s[m,n]$ can be approximated by finite differences

$∂ s ∂ x = s [ m + 1 , n ] - s [ m , n ] = ( s * h x ) [ m , n ] ∂ s ∂ y = s [ m , n + 1 ] - s [ m , n ] = ( s * h y ) [ m , n ] ,$

which are convolutions with the following filters

$h x ( 1 ) = 1 - 1 , h y ( 1 ) = 1 - 1 .$

Here in the matrix form, row and column indexes correspond to $x$ (first) and $y$ (second) dimensions, respectively; and the sample in the box corresponds to the original (i.e. $(m,n)=(0,0)$).

Example 2 (The Laplacian and image sharping) The Laplacian of a two-dimensional signal $s(x,y)$ is defined as

$Δ s = ∂ 2 s ∂ x 2 + ∂ 2 s ∂ x 2$

Extending the definition of the first-order derivatives above to the second-order derivatives, we have

$∂ 2 s ∂ x 2 = s [ m + 1 , n ] - 2 s [ m , n ] + s [ m - 1 , n ] ∂ 2 s ∂ y 2 = s [ m , n + 1 ] - 2 s [ m , n ] + s [ m , n - 1 ] .$

It follows that the Laplacian can be computed as

$Δ s [ m , n ] = s [ m + 1 , n ] + s [ m - 1 , n ] + s [ m , n + 1 ] + s [ m , n - 1 ] - 4 s [ m , n ] = ( s * h Δ ) [ m , n ] ,$

which is convolution with the following filter

$h Δ = 0 1 0 1 - 4 1 0 1 0 .$

Sometimes, the Laplacian is extended by adding two more terms for two diagonal directions leading to the following filter

$h Δ ' = 1 1 1 1 - 8 1 1 1 1 .$

Since the Laplacian is a derivative operator, it highlights intensity discontinuities (or edges) in an image. We can sharpen an image by adding the negative of the Laplacian image to the original

$r [ m , n ] = s [ m , n ] - Δ s [ m , n ] .$

Using the Laplacian filter [link] we can write the resulting sharpening operation as a convolution

$r [ m , n ] = ( s * h sharp ) [ m , n ] ,$

with the following filter

$h sharp = - 1 - 1 - 1 - 1 9 - 1 - 1 - 1 - 1 .$

Example 3 (Gausian smoothing filter) The two-dimensional Gaussian filter, which is often used for image smoothing, is defined as

$h Gauss (2D) ( x , y ) = 1 2 π σ 2 e - x 2 + y 2 2 σ 2 .$

The 2D Gaussian filter is separable, which means it is a product of 1D filters in each dimension

$h Gauss (2D) ( x , y ) = h Gauss (1D) ( x ) h Gauss (1D) ( y ) , where h Gauss (1D) ( x ) = 1 2 π σ e - x 2 2 σ 2 .$

Discrete-domain Gaussian filters used in practice are sampled and truncated versions of the above continuous-domain filters.

Example 4 (Sobel edge detector) The Sobel edge detector is obtained by smoothing the image in the perpendicular direction before computing the directional derivatives. The associate Sobel edge detector filters are given by

$h x (Sobel) = 1 2 1 0 0 0 - 1 - 2 - 1 , h y (Sobel) = 1 0 - 1 2 0 - 2 1 0 - 1 .$

Edges are detected as pixels $[m,n]$ where the magnitude of the gradient is above a certain threshold $T$; i.e.

$| ( s * h x (Sobel) ) [ m , n ] | + | ( s * h y (Sobel) ) [ m , n ] | ≥ T .$

## Frequency response of a filter

A key result in signal and image processing is that convolution in the space domain becomes multiplication in the frequency domain

$r ( x , y ) = ( s * h ) ( x , y ) ⟷ FT R ( u , v ) = S ( u , v ) H ( u , v ) , r [ m , n ] = ( s * h ) [ m , n ] ⟷ DSFT R d ( u , v ) = S d ( u , v ) H d ( u , v ) .$

Therefore, the Fourier transform $H(u,v)$ of the filter, called frequency response, indicates how certain frequency components of the input image $s(x,y)$ are amplified or attenuated in the resulting filtered image $r(x,y)$.

However, multiplication in the DFT domain corresponds to circular convolution in the space domain

$r [ m , n ] = ( s ⊛ M , N h ) [ m , n ] ⟷ DFT R [ k , l ] = S [ k , l ] H [ k , l ] .$

The circular convolution operation for images of size $M×N$ is defined as

$( s ⊛ M , N h ) [ m , n ] = ∑ m ' = 0 M - 1 ∑ n ' = 0 N - 1 h [ m ' , n ' ] s [ 〈 m - m ' 〉 M , 〈 n - n ' 〉 N ] ,$

where $〈n〉N$ denotes modulo $N$ of $n$.

Problems

1. Prove the property [link].
2. Prove the property [link].
3. Find and sketch the frequency response of the Laplacian filter given in [link]. Is this a lowpass or highpass filter?
4. Develop and write a simple Matlab function to detect edges at $45∘$ in an image using the
``conv2``
function.

# Image Denoising

Sometimes we obtain images that are contaminated by noise. Image denoising aims to remove or reduce noise present in the image. Some common image denoising methods are smooth filtering, median filtering, Wiener filtering, wavelet thresholding or shrinkage.

## Wiener filter

We now present the removal of additive noise using Wiener filtering; first in 1D and then extend to 2D. Let $s[n]$ denote the original clean signal and $z[n]$ denote the noise. The obtained noisy signal is then

$r [ n ] = s [ n ] + z [ n ] .$

The goal is to filter the obtained noisy signal $r[n]$ with a linear filter $g[n]$ so that the output $s^[n]=(r*g)[n]$ is a an estimate of the clean signal $s[n]$. The Wiener filter $g[n]$ minimizes the expectation of the squared error given by

$J = E { ( s ^ [ n ] - s [ n ] ) 2 } = E ∑ k g [ k ] r [ n - k ] - s [ n ] 2$

Differentiating the above with respect to $g[k]$ and setting the result to zero, we get

$∂ J ∂ g [ k ] = 2 E ( s ^ [ n ] - s [ n ] ) r [ n - k ] = 0 .$

Therefore, the Wiener filter $g[k]$ has to satisfies the following (called orthogonal condition)

$E { s ^ [ n ] r [ n - k ] } = E { s [ n ] r [ n - k ] } .$

If we assume that $s[n]$ and $z[n]$ are wide-sense stationary random processes, then so is $r[n]$ given by [link]. That means their auto-correlation and cross-correlation functions only depend on the difference of sample indexes. For example, we can write

$R s s [ k ] = E { s [ n ] s [ n - k ] } , R s r [ k ] = E { s [ n ] r [ n - k ] } .$

Since $s^[n]=(r*g)[n]$, we have

$E { s ^ [ n ] r [ n - k ] } = E { ∑ m g [ m ] r [ n - m ] r [ n - k ] } = ∑ m g [ m ] E { r [ n - m ] r [ n - k ] } = ∑ m g [ m ] R r r [ k - m ] = ( g * R r r ) [ k ]$

Therefore the orthogonal condition [link] for the Wiener filter $g[k]$ can be rewritten as

$( g * R r r ) [ k ] = R s r [ k ] .$

Taking the DTFT of the above we obtain a closed form formula for the Wiener filter

$G ( u ) = P s r ( u ) P r r ( u ) ,$

where function $P$ (called power spectral density) is defined as the DTFT of the corresponding correlation function $R$. For example, if $z[n]$ is a white Gaussian random process with zero mean and variance $σ2$, then $Rzz[k]=σ2δ[k]$ and $Pzz(u)=σ2$.

If furthermore, we assume that the clean signal $s[n]$ and the noise $z[n]$ are uncorrelated, and the noise $z[n]$ has zero mean. Then

$E { s [ n ] z [ k ] } = E { s [ n ] } E { z [ k ] } = 0 for all n , k .$

It follows that

$R s r [ k ] = R s s [ k ] , R r r [ k ] = R s s [ k ] + R z z [ k ] ,$

Therefore, we can express the Wiener filter for denoising solely in terms of the power spectral densities of the clean signal $s[n]$ and the noise $z[n]$ as

$G ( u ) = P s s ( u ) P s s ( u ) + P z z ( u )$

From this, it is easy to see that the Wiener filter for image denoising is given by

$G ( u , v ) = P s s ( u , v ) P s s ( u , v ) + P z z ( u , v )$

## Wavelet denoising

Wavelet-based methods have established as a state-of-the-art denoising approach for additive white Gaussian noise (AWGN)

Another common type of noise is the impulse noise, which is sometimes referred to as speckle or salt and pepper noise. For this type of noise, normally a median filtering can do a good job.
. The basic idea is that due to its excellent approximation property for typical signals that are piecewise smooth, in the wavelet domain most of the signal information is captured in a few significant wavelet coefficients. As a result, the wavelet coefficients of the original signal stand out from the noisy one (after an orthonormal transform, an AWGN becomes another AWGN of the same variance). Therefore, a simple thresholding in the wavelet domain can effectively remove the noise out of the signal. [link] illustrates this basic concept of wavelet thresholding.

Thresholding can be either “soft thresholding” or “hard thresholding”. A hard thresholding estimator is implemented with

$h T ( x ) = x if | x | > T , 0 if | x | ≤ T .$

A soft thresholding estimator is implemented with

$s T ( x ) = x - T if x > T , x + T if x < T , 0 if | x | ≤ T .$

What is remarkable is such a simple wavelet thresholding algorithm achieves optimal performance in certain sense. The asymptotically optimal threshold (also referred to as “the universal threshold”) is $T=σ2logeN$ for signals of length $N$. In practice, lower threshold such as $T=3σ$ improves the MSE significantly.

There are several variations of the basic wavelet thresholding scheme that offer performance gains.

• Translation invariant thresholding. An improved thresholding estimator called “cycle-spinning” is calculated by averaging estimators for translated versions of the signal. The algorithm is equivalent to thesholding the non-subsampled or `à trous' wavelet transform.
• Spatially adaptive thresholding.
• Denoising based on statistical modeling of wavelet coefficients.

Problems

2. Derive the denoising algorithm using Wiener filter when both $s[n]$ and $z[n]$ are white Gaussian random processes with zero mean and variances $σs2$ and $σz2$, respectively.
3. Suppose that you want to denoise an image that was contaminated by an additive white Gaussian noise. In the wavelet domain (via an orthornormal transform), the wavelet coefficient of the noisy image $y$ is related to the corresponding wavelet coefficient of the clean image $x$ by:
$y=x+n,$
where $n$ is Gaussian random variable with the following PDF:
$p(n)=const·e-n22σ2.$
For typical natural images, $x$ can be modeled as a generalized Gaussian random variable with the following PDF:
$p(x)=const·e-|x|αβ.$
The parameters $σ$, $α$ and $β$ can be estimated beforehand. Derive a simple closed form formula for the maximum a posteriori estimate $x^MAP$ from $y$:
$x^MAP=argmaxxp(x|y),$
using $σ$, $α$ and $β$.

# Image Deconvolution

Suppose that the obtained image is a noisy convoluted version of the original image

$r [ m , n ] = ( s * h ) [ m , n ] + z [ m , n ] .$

## Inverse filtering and Wiener filtering

Here we assume that the convolution filter $h$ is known and fixed. Without noise, a simple way to recover $s$ from $r$ is by inverse filtering

$S ( u , v ) = R ( u , v ) H ( u , v ) = R ( u , v ) 1 H ( u , v ) ︸ G ( u , v ) .$

However, in the presence of noise, then inverse filtering the obtained signal in [link] leads to

$R ( u , v ) H ( u , v ) = S ( u , v ) + Z ( u , v ) H ( u , v ) .$

Hence noise might be “blew up” in the second term of the above equation near the frequency such that $H(u,v)$ is close to zero.

Following the derivation of the Wiener filter for denoising in the last section, we obtain the Wiener filter for deconvolution as:

$G ( u , v ) = H * ( u , v ) P s s ( u , v ) | H ( u , v ) | 2 P s s ( u , v ) + P z z ( u , v ) = 1 H ( u , v ) | H ( u , v ) | 2 | H ( u , v ) | 2 + P z z ( u , v ) / P s s ( u , v ) .$

## Vector-space approach

In the vector-space framework, the deconvolution problem $(h*s)[m,n]=r[m,n]$ can be written in the familiar matrix-vector equation

$A x = b ,$

where $x$ is the column vector corresponding to the unknown image $s[m,n]$, $b$ is the column vector corresponding to the given filtering image $r[m,n]$, and $A$ is the matrix corresponding to convolution with the filter $h[m,n]$. For example, 1D filtering a signal ${s[n]}n=02$ by a filter ${h[n]}n=01$ can be expressed as

$h [ 0 ] 0 0 h [ 1 ] h [ 0 ] 0 0 h [ 1 ] h [ 0 ] 0 0 h [ 1 ] ︸ A s [ 0 ] s [ 1 ] s [ 2 ] ︸ x = r [ 0 ] r [ 1 ] r [ 2 ] r [ 3 ] ︸ b$

The matrix-vector formulation allows us to resort to a rich body of literature on solving linear inverse problems.

In the presence of additive noise, our problem becomes

$A x + z = b ,$

where $z$ is a random vector. Given data $b$, a common approach is to search for the maximum-likelihood (ML) solution:

$x ML * = max x Pr ( b | x ) .$

If $z$ is a white Gaussian noise of zero mean and variance $σ2$, then given $x$, $b=Ax+z$ is also a Gaussian random vector of mean $Ax$ and variance $σ2$. That is

$Pr ( b | x ) = 1 ( 2 π σ 2 ) d / 2 e - ∥ b - A x ∥ 2 2 2 σ 2 .$

In this case the ML solution of [link] is

$x ML * = arg min x ∥ A x - b ∥ 2 2 ,$

which is also the least-squares (LS) solution of [link].

If prior knowledge of $x$ is known in the form of $Pr(x)$, then we can search for the maximum a priori (MAP) solution

$x MAP * = arg max x Pr ( x | b ) = arg max x { Pr ( b | x ) Pr ( x ) } .$

The second term in [link] is often called penalty or regularization term.

If we again assume the noise $z$ is white Gaussian then the MAP solution becomes

$x MAP * = arg min x { - log Pr ( b | x ) - log Pr ( x ) } = arg min x 1 2 σ 2 ∥ A x - b ∥ 2 2 - log Pr ( x ) ,$

which is also called a regularized or penalized LS solution of [link].

Consider the ML or LS solution in [link]. From linear algebra, we know that this solution can be obtained by the pseudo-inverse $xML/LS*=A†b=(ATA)-1ATb$. However, in image processing applications, the size of $x$ is typically in the order of millions samples, and thus storing matrix $A$ of size $1˜06×106$ and computing its pseudo-inverse are impractical. On the other hand, although $A$ is a big matrix, it has compact description (convolution with a filter $h[m,n]$) and fast algorithms for computing multiplications by $A$ and $AT$ (convolution with short FIR filters or using FFT). Under these circumstances, iterative methods offer a practical solution.

## Iterative methods for linear inversion

Generally, in an iterative method, there is a sequence of iterations where each one would successively improve the solution and hopefully they would quickly converse to the true solution. To measure the improvement at each iteration, an objective function is used.

The ML/LS solution of [link] is the minimizer of the following objective function

$f ( x ) = ∥ A x - b ∥ 2 2 = ( A x - b ) T ( A x - b ) = x T A T A x - 2 b T A x + b T b .$

The gradient of this objective function is

$▽ f ( x ) = 2 A T A x - 2 A T b .$

A direct solver would set $▽f(x)=0$ which leads to the normal equation $ATAx=ATb$ and the solution using the pseudo-inverse.

The search for the minimizer of $f(x)$ is a multidimensional search problem. Therefore, at each iteration step we should limit the search space to a one-dimensional search. Specifically, assume that $xn$ the current estimate for the solution. Then for the next step, we will fix a search direction$pn$ and limit the search space along that direction:

$x n + 1 = x n + α n p n ,$

where $αn∈R$ is chosen as the optimal step size

$α n = arg min α f ( x n + α p n ) .$

This optimal step size $αn$ is the solution of the following equation

$d f d α ( x n + α p n ) = 0$

which can be easily derived to be

$α n = p n T A T ( b - A x n ) p n T A T A p n .$

Therefore, with an initial guess $x0$ and a sequence of search directions $p0,p1,...$, the iterative search process is completely determined.

A common choice for search directions is to use negative of the gradient of $f$ at each step, or

$p n = - ▽ f ( x n ) = 2 A T ( b - A x n ) .$

This choice of directions leads to the steepest descent procedure. Instead of using the optimal step size $αn$ at each iteration given in [link], which incurs some computation cost, a simple approach is to use a fixed step for all iteration. This leads to the following iteration which is known as Landweber

$x n + 1 = x n + λ A T ( b - A x n ) .$

We want to examine the convergence property of the Landweber iteration. Rewrite [link] as

$x n + 1 = ( I - λ A T A ) ︸ C x n + λ A T b ︸ d ,$

then we have

$x n = C n x 0 + ( I + C + C 2 + ... + C n - 1 ) d .$

Therefore, the sequence $xn$ converges if all the eigenvalues of $C=I-λATA$ has magnitude less than 1. This condition is equivalent to

$| 1 - λ λ i | < 1 for all eigenvalues λ i of A T A ⟺ 0 < λ < 2 λ max ( A T A ) .$

Under this condition we have

$lim n → ∞ C n = 0 lim n → ∞ ( I + C + C 2 + ... + C n - 1 ) = ( I - C ) - 1 ,$

and thus

$lim n → ∞ x n = ( A T A ) - 1 A T b ,$

which is, as expected, the same as the solution using the pseudo-inverse.

Problems

1. Derive the Wiener for deconvolution given in [link].
2. Adjoint of convolution Iterative methods often requires fast algorithm for computing multiplication with $AT$, or adjoint of $A$.
1. Suppose that $C˜$ is the circular convolution matrix of size $N×N$ for the circular convolution with a length $N$ real filter $h=[h0,h1,...,hN-1]$. Develop a MATLAB function for a fast implementation of $C˜Ty$ using
``fft``
and
``ifft``
commands.
2. Suppose that $C$ is the linear convolution matrix of size $(L+N-1)×L$ for the linear convolution of a length $N$ real FIR filter $h=[h0,h1,...,hN-1]$ with a length $L$ input signal. Suppose that $L$ is much greater than $N$. Develop a MATLAB function for a fast implementation of $CTy$ using the
``conv``
command, where
``````w = conv(u,v) convolves vectors u and v.
``````

# Image Interpolation

Interpolation can be considered as the inverse operation of sampling: it aims to reconstruct signal value at any locations given samples on a sampling grid. Image interpolation has numerous applications including general geometric transformations (image zooming, resizing, rotation, and warping) and digital-to-analog conversions.

We assume that the continuous-domain image $s(x,y)$ at resolution $(Δx,Δy)$ has the following representation

$s ( x , y ) = ∑ m = - ∞ ∞ ∑ m = - ∞ ∞ c [ m , n ] ψ ( x - m Δ x , y - n Δ y ) .$

Taking the Fourier transform of both sides of [link] we obtain the following equivalent representation in the frequency domain

$S ( u , v ) = C d ( Δ x u , Δ y v ) Ψ ( u , v ) .$

For a fixed $ψ$, the space of all possible functions $s(x,y)$ that can be represented in the form [link] for some discrete-domain signal ${c[m,n]}$ is a subspace, called shift-invariant subspace.

Three basic examples of shift-invariant subspaces are:

1. For
$ψ(x,y)=1if0≤x<Δx,0≤y<Δy0otherwise,$
we have a subspace of piecewise-constant functions with pieces are rectangular regions $(x,y)∈[mΔx,(m+1)Δx)×[nΔy,(n+1)Δy)$. This corresponds to a zero-order-hold (ZOH) digital-to-analog (D/A) converter.
2. For
$ψ(x,y)=sinc(x/Δx)sinc(y/Δy)$
we have a subspace of bandlimited functions to frequencies $|u|≤1/(2Δx),|v|≤1/(2Δy)$. This corresponds to an ideal D/A converter. In this case, applying [link]-[link] we obtain [link].
3. Between these two above extremes, if $ψ(x,y)$ is a B-spline function then we have a subspace of spline functions.

We consider the generalized sampling scheme where sampled image is given by

$w 1 [ m , n ] = ∫ - ∞ ∞ ∫ - ∞ ∞ s ( x , y ) φ 1 ( x - m Δ x , y - n Δ y ) d x d y = ( s * φ 1 ¯ ) ( m Δ x , n Δ y ) .$

Here we use the notation $φ¯(x,y)=φ(-x,-y)$ so that correlation can be written as a convolution.

If $s(x,y)$ has the form of [link] then substituting it to [link] we have

$w 1 [ m , n ] = ( c * b 1 ) [ m , n ] ,$

where

$b 1 [ k , l ] = ( ψ * φ 1 ¯ ) ( k Δ x , l Δ y ) .$

Note that in the last two convolutions, the first one is in the discrete domain while the second one is in the continuous domain.

Given $ψ$ and $φ1$ we can compute $b1$. Then solving the deconvolution problem [link] we can recover coefficients $c$ from sampled data $w$. Finally, using [link] we can compute the image $s(x,y)$ at any location$(x,y)$.

Now suppose that we are given sampled image at a lower resolution

$w 2 [ m , n ] = ∫ - ∞ ∞ ∫ - ∞ ∞ s ( x , y ) φ 2 ( x - 2 m Δ x , y - 2 n Δ y ) d x d y = ( s * φ 2 ¯ ) ( 2 m Δ x , 2 n Δ y ) .$

Then similarly, we have

$w 2 [ m , n ] = ( c * b 2 ) [ 2 m , 2 n ] ,$

where

$b 2 [ k , l ] = ( ψ * φ 2 ¯ ) ( k Δ x , l Δ y ) .$

In words, sampled image $w2$ is obtained from $c$ by filtering followed by downsampling.

Problems

3. Suppose that the input image $s(x,y)$ can be represented by a linear expansion
$s(x,y)=∑m∈Z∑n∈Zc[m,n]ψ(x-mΔ,y-nΔ),$
where $φ$ is a chosen compactly supported basis function (e.g. B-spline). Our imaging system over-samples the convolution of $s(x,y)$ with a known and compactly supported filter $φ(x,y)$ and return discrete image
$r[m,n]=(s*φ)(mΔ/2,nΔ/2),form,n∈Z.$
1. Derive a fast algorithm to implement the “forward” linear operator $A$ that maps discrete image ${c[m,n]}$ into discrete image ${r[m,n]}$.
2. Derive a fast algorithm to implement $AT$, the adjoint operator of $A$.

# Image Reconstruction

## Image reconstruction from projections

In several important applications such as computer tomography (CT), we can collect projection data of an object and would like to reconstruct the internal view of the object. The projection data of an image $s(x,y)$ is defined as

$p θ ( t ) = ∫ - ∞ ∞ s ( t cos θ - r sin θ , t sin θ + r cos θ ) d r = ∫ - ∞ ∞ ∫ - ∞ ∞ s ( x , y ) δ ( t - x cos θ - y sin θ ) d x d y .$

Theorem 2 (Projection-slice) Let $Pθ(f)$ be the one-dimensional Fourier transform of $pθ(t)$. Then

$P θ ( f ) = S ( f cos θ , f sin θ ) .$

From the definition of the Fourier transform and projection we have

$P θ ( f ) = ∫ - ∞ ∞ p θ ( t ) e - j 2 π f t d t = ∫ - ∞ ∞ ∫ - ∞ ∞ s ( t cos θ - r sin θ , t sin θ + r cos θ ) e - j 2 π f t d r d t$

Apply the following change of variables (in fact a rotation)

$x = t cos θ - r sin θ y = t sin θ + r cos θ ⟺ t = x cos θ + y sin θ r = - x sin θ + y cos θ$

to the last integral we have

$P θ ( f ) = ∫ - ∞ ∞ ∫ - ∞ ∞ s ( x , y ) e - j 2 π f ( x cos θ + y sin θ ) d x d y = S ( f cos θ , f sin θ ) .$

Given $S(fcosθ,fsinθ)$ take the inverse Fourier transform using polar coordinates to recover $s(x,y)$

$s ( x , y ) = ∫ 0 2 π ∫ - ∞ ∞ S ( f cos θ , f sin θ ) e j 2 π ( x f cos θ + y f sin θ ) f d f d θ = ∫ 0 π ∫ - ∞ ∞ S ( f cos θ , f sin θ ) e j 2 π f ( x cos θ + y sin θ ) | f | d f d θ = ∫ 0 π ∫ - ∞ ∞ P θ ( f ) | f | e j 2 π f ( x cos θ + y sin θ ) d f d θ .$

Let

$Q θ ( f ) = P θ ( f ) | f | ,$

then we have

$s ( x , y ) = ∫ 0 π q θ ( x cos θ + y sin θ ) d θ .$

From [link] we see that $qθ(t)$ is obtained by filtering $pθ(t)$. The operation [link] backproject to the spatial domain. Thus this called “filter-backprojection” algorithm for image reconstruction.

## General linear inversion

In more general setting, we have access to linear measurements of an unknown image and after discretization, image reconstruction amounts to solve the following linear inverse problem

$A x = b ,$

where $A$ is a matrix representing the linear image-data relationship, $x$ is a column vector representing the unknown image, and $b$ is a column vector representing the collected data. We can use tools for solving linear inverse problems as described before.

Problems

1. Suppose that you would like to reconstruct an unknown image from its given projection data. You planned to use the filtered back projection method, but forgot the ramp filtering step in [link] and instead back projected directly the projection data. Even worse, the acquired projection data were deleted! Study this problem in the continous-domain, show how can you recover the true image from the wrongly reconstructed image? What would be the major problem in this recovery and how can you overcome it?
2. Projection and back projection Consider the two by two image
$S=s11s12s21s22$
Given the “row projections” $pjR=∑isij$ and the “column projections” $piC=∑jsij$, can the elements of the image be recovered? Explain. Find the image $S^$ that is obtained by back projecting row and column projections, and compare $S^$ with $S$.
3. Regularization using SVD Consider the following regularized inverse problem
$x*=argminx{∥Ax-b∥22+λ∥x-m∥22}$
for solving $Ax=b$ while enforcing the solution $x$ to be closed to a known vector $m$. Using the SVD to diagonalize the above minimization problem and find the solution $x*$ explicitly.

### Tập tin đính kèm

Tải về
Đánh giá:
0 dựa trên 0 đánh giá
Nội dung cùng tác giả

Nội dung tương tự