Homec4science
Diffusion epfl (master) Browse

R2653/writing/sp4comm.multipub/130-image/10-improc.tex
f90ad791b9a3master

Last Change
Raw File

10-improc.tex
View Options

\section{Basic Signals and Operators}
The simplest two-dimensional signal is the 2D impulse:
\begin{equation}
\delta[n_1, n_2] = \begin{cases}
1 & \mbox{if $n_1= n_2 = 0$} \\
0 & \mbox{otherwise.}
\end{cases}
\end{equation}
The 2D equivalent of a zero-centered rectangular signal is
\begin{equation}
\rect\left(\frac{n_1}{2N_1}, \frac{n_2}{2N_2}\right) = \begin{cases}
1 & \mbox{if $|n_1| < N_1$ and $|n_2| < N_2$} \\
0 & \mbox{otherwise;}
\end{cases}
\end{equation}
%more generally, signals that are a nonzero constant over a given finite support are called {\em indicator functions} for the support.
these two signals are represented in Figure~\ref{basic2D} using support plots.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\begin{center}
\begin{tabular}{cc}
\begin{dspCP}[xticks=5,yticks=5,xout=true]{-6,6}{-6,6}
\dspPoints{0 0}
\end{dspCP}
&
\begin{dspCP}[xticks=5,yticks=5,xout=true]{-6,6}{-6,6}
\dspPoints{-3 2 -2 2 -1 2 0 2 1 2 2 2 3 2
-3 1 -2 1 -1 1 0 1 1 1 2 1 3 1
-3 0 -2 0 -1 0 0 0 1 0 2 0 3 0
-3 -1 -2 -1 -1 -1 0 -1 1 -1 2 -1 3 -1
-3 -2 -2 -2 -1 -2 0 -2 1 -2 2 -2 3 -2 }
\end{dspCP}
\\
{\sf (a)} & {\sf (b)}
\end{tabular}
\end{center}
\caption{Two dimensional impulse (a) and rectangular signal for $N_1=4$ and $N_2=3$.}\label{basic2D}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The impulse and the rect signals are examples of {\em separable} 2D signals, i.e. signals that can be written as the product of two independent 1D signals:
\[
x[n_1, n_2] = x_1[n_1]x_2[n_2].
\]
It is indeed obvious that
\[
\delta[n_1, n_2] = \delta[n_1]\delta[n_2]
\]
and that
\[
\rect\left(\frac{n_1}{2N_1}, \frac{n_2}{2N_2}\right) = \rect\left(\frac{n_1}{2N_1}\right) \rect\left(\frac{n_2}{2N_2}\right).
\]
Although separable signals constitute a rather ``artificial'' subset of 2D signals, they are very convenient in conjunction with filtering operations, as we will see shortly. To realize that separability is entirely dependent on the reference system consider for instance the diamond-shaped signal
\[
x[n_1, n_2] = \begin{cases}
1 & \mbox{if $|n_1| + |n_2| < N$} \\
0 & \mbox{otherwise}
\end{cases}
\]
which shown in Figure~\ref{nonsep2D}-(a) for $N=3$; clearly this is a 45-degree rotation of a square separable signal, yet the diamond signal is not separable. Similarly, the rectangular frame
\[
x[n_1, n_2] = \rect(n_1/2N, n_2/2N) - \rect(n_1/2M, n_2/2M)
\]
shown in Figure~\ref{nonsep2D}-(b) for $N = 4$ and $M=2$, although obtained as a linear combination of separable signal is not itself separable.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}\label{nonsep2D}
\begin{center}
\begin{tabular}{cc}
\begin{dspCP}[xticks=5,yticks=5,xout=true]{-6,6}{-6,6}
\dspPoints{0 2
-1 1 0 1 1 1
-2 0 -1 0 0 0 1 0 2 0
-1 -1 0 -1 1 -1
0 -2}
\end{dspCP}
&
\begin{dspCP}[xticks=5,yticks=5,xout=true]{-6,6}{-6,6}
\dspPoints{-3 3 -2 3 -1 3 0 3 1 3 2 3 3 3
-3 2 -2 2 -1 2 0 2 1 2 2 2 3 2
-3 1 -2 1 2 1 3 1
-3 0 -2 0 2 0 3 0
-3 -1 -2 -1 2 -1 3 -1
-3 -2 -2 -2 -1 -2 0 -2 1 -2 2 -2 3 -2
-3 -3 -2 -3 -1 -3 0 -3 1 -3 2 -3 3 -3 }
\end{dspCP}
\\
{\sf (a)} & {\sf (b)}
\end{tabular}
\end{center}
\caption{Diamond signal (a) and rectangular frame (b).}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The basic operations between 2D signals are the intuitive extension of their 1D counterparts. In particular, the convolution between two 2D signals, which defines the operation of processing an image with a linear, space-invariant system, is defined as
\begin{equation}\label{2dconv}
x[n_1, n_2] \ast h[n_1, n_2] = \sum_{k_1 = -\infty}^{\infty}\sum_{k_2 = -\infty}^{\infty} x[k_1,k_2]h[n_1-k_1,n_2-k_2].
\end{equation}
If $h[n_1, n_2] = h_1[n_1]h_2[n_2]$, i.e. if $h[n_1, n_2]$ is separable, we can rewrite~(\ref{2dconv}) as
\begin{align}\label{2dconv}
x[n_1, n_2] \ast h[n_1, n_2] &= \sum_{k_1 = -\infty}^{\infty}h_1[n_1-k_1] \sum_{k_2 = -\infty}^{\infty} x[k_1,k_2]h_2[n_2-k_2] \\
&= h_1[n_1] \ast (h_2[n_2] \ast x[n_1,n_2]).
\end{align}
From the algorithmic point of view, therefore, the convolution with a separable signal is particularly simple since it becomes the cascade of two one-dimensional convolutions; stated differently, if $h[k_1,k_2]$ is a separable impulse response, the filtering operation is equivalent to filtering each image row with the first component of the impulse response, followed by filtering each column by the second component. In the case of FIR filters, this reduces the number of operations {\em per output sample} from $M_1M_2$ to $M_1+M_2$, where $M_1\times M_2$ is the size of the filter's support.
\subsection{Translation, Scaling and Rotations}
An arbitrary function of two {\em real} variables $x(t_1, t_2)$ can be reshaped via an affine transform of the underlying coordinate system. An affine transform is defined by the formula
\begin{equation}
\begin{bmatrix} t'_1 \\ t'_2 \end{bmatrix} =
\begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix} \begin{bmatrix} t_1 \\ t_2 \end{bmatrix} - \begin{bmatrix} d_1 \\ d_2 \end{bmatrix} =
\mathbf{A} \begin{bmatrix} t_1 \\ t_2 \end{bmatrix} - \mathbf{d} \nonumber
\end{equation}
where $\mathbf{A}$ defines a linear mapping of the reference system while the vector $\mathbf{d}$ determines a translation of the origin. Particular forms for $\mathbf{A}$ and $\mathbf{d}$ give rise to special cases:
\begin{itemize}
\item {\bf Translation:} if $\mathbf{A} = \mathbf{I}$, the result is a simple translation of the origin to the point $(d_1, d_2)$.
\item {\bf Scaling:} if $\mathbf{A} = \mbox{diag}(a_1, a_2)$ the function is rescaled by $a_1$ and $a_2$ along the two dimensions. If $a_1 = a_2$ the scaling is uniform for both axes and the rescaling is said to preserve the {\em aspect ratio} of the original function.
\item {\bf Rotation:} if $\mathbf{A}$ is of the form
\[
\mathbf{A} = \begin{bmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{bmatrix}
\]
the result is a counterclockwise rotation around the origin by an angle $\theta$.
\item {\bf Shearing:} if $\mathbf{A}$ is of the form
\[
\mathbf{A} = \begin{bmatrix} 1 & s \\ 0 & 1 \end{bmatrix}
\]
the result is a horizontal shearing of the function, i.e. a horizontal slanting of the support aligned with the line $t_2 = st_1$. The transpose of the above matrix produces a vertical shear.
\item {\bf Flips:} the matrix
\[
\mathbf{A} = \begin{bmatrix} -1 & 0 \\ 0 & 1 \end{bmatrix}
\]
flips the function around the $t_2$ axis (horizontal flip) while the matrix
\[
\mathbf{A} = \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix}
\]
flips the function around the $t_1$ axis (vertical flip)
\end{itemize}
In the case of a digital image $x[n_1,n_2]$ the coordinate indices must be integers so that, in principle, we can only use affine transforms for which $(\mathbf{A}[n_1,n_2]^T - \mathbf{d}) \in \mathbb{N}^2$. This limits our toolbox to integer translations, subsampling (scaling) by an integer factor, flips and rotations by multiples of 90 degrees only. A way out of the impasse is to consider a continuous-space interpolation of the digital image; the resulting function of two real variables can be then arbitrarily transformed and the result can be sampled back on the integer grid. Interpolation of 2D signals is conceptually identical to the 1D case, leading to a global 2D sinc interpolation formula and an associated sampling theorem for bandlimited 2D functions. Instead of following that route, however, we will concentrate here on purely algorithmic (and therefore {\em local}) interpolation schemes that allow us to apply arbitrary affine transforms to a digital image.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}\label{bilinearInt}
\begin{center}
\begin{dspCP}[width=6cm,xticks=custom,yticks=custom,xout=true]{0.5,2.5}{0.5,2.5}
\dspPointValue{2,2}{$x[\eta_1+1,\eta_2+1]$}
\dspPointValue{1,1}{$x[\eta_1,\eta_2]$}
\dspPointValue{2,1}{}
\dspPointValue{1,2}{}
\dspPointValue[dotstyle=o]{1.3,1.7}{$(t_1,t_2)$}
\psline[linewidth=0.5pt]{<->}(1,1)(1,1.7)\dspText(0.9,1.3){$\tau_2$}
\psline[linewidth=0.5pt]{<->}(1,1.7)(1.3,1.7)\dspText(1.14,1.75){$\tau_1$}
\psline[linewidth=0.5pt,linestyle=dotted](1,1)(2,1)
\psline[linewidth=0.5pt,linestyle=dotted](1,2)(2,2)
\psline[linewidth=0.5pt,linestyle=dotted](1.3,1)(1.3,2)
\dspPointValue[dotstyle=triangle]{1.3,1}{}
\dspPointValue[dotstyle=triangle]{1.3,2}{}
\psline[linewidth=0.5pt]{<->}(1,1.7)(1.3,1.7)\dspText(1.14,1.75){$\tau_1$}
\dspCustomTicks[axis=x]{1 $\eta_1$ 2 $\eta_1+1$}
\dspCustomTicks[axis=y]{1 $\eta_2$ 2 $\eta_2+1$}
\end{dspCP}
\end{center}
\caption{Bilinear interpolation framework.}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Given an $N_1\times N_2$ image $x[n_1, n_2]$ and an $(\mathbf{A}, \mathbf{d})$ pair, call $y[m_1, m_2]$ the image we would obtain if we could apply the associated affine transform to $x$. The key idea is to consider the {\em inverse} transformation, which maps each pixel coordinate in $y$ back to the space where $x$ lives:
\[
\begin{bmatrix} t_1 \\ t_2 \end{bmatrix} = \mathbf{A}^{-1} \begin{bmatrix} m_1 + d_1 \\ m_2 + d_2 \end{bmatrix};
\]
here $t_1$ and $t_2$ are obviously real variables and, as such, they will not fall on the original integer grid. If either $t_1$ or $t_2$ fall outside the $[0,N_1-1]$ and $[0,N_2-1]$ respectively, then we simply set $y[m_1, m_2] = 0$ (or any other pre-arranged value) and we're done. Otherwise we can express the coordinate pair as
\[
(t_1, t_2) = (\eta_1 + \tau_1, \eta_2 + \tau_2)
\]
with $\eta_{1,2} \in \mathbb{N}$ and with $0 \le \tau_{1,2} < 1$. This notation highlights the fact that the coordinate pair points to a ``subpixel'' located at a distance $(t_1, t_2)$ from $x[\eta_1, \eta_2]$ and contained within a four-pixel square bounded by $x[\eta_1, \eta_2]$ and $x[\eta_1+1,\eta_2+1]$ (see Figure~\ref{bilinearInt}). A popular method to obtain an estimate of the value for this subpixel is the {\em bilinear interpolation} algorithm, a straightforward extension of the first order interpolator in the one-dimensional domain. Imagine applying the 1D triangular interpolator a row in the image; the result is a straight line interpolation between neighboring pixels in the line so that for each $\tau$ between zero and one we can write
\[
x_c(\eta_1 + \tau_1, n_2) = (1-\tau_1)x[\eta_1, n_2] + \tau_1 x[\eta_1 + 1, n_2].
\]
We can compute the above formula for $n_2 = \eta_2$ and $n_2 = \eta_2+1$ in order to obtain the interpolated points indicated by triangles in Figure~\ref{bilinearInt}; we can then apply the same first order interpolation column-wise, and interpolate between the ``triangles'' to obtain:
\[
x_c(\eta_1 + \tau_1, \eta_2 + \tau_2) = (1-\tau_2)x_c(\eta_1 + \tau_1, \eta_2) + \tau_2 x_c(\eta_1 + \tau_1, \eta_2).
\]
By combining the two expressions we finally obtain the bilinear estimate for the transformed image as
\begin{eqnarray}
y[m_1,m_2] =& (1-\tau_1)(1-\tau_2)x[\eta_1, \eta_2] + \tau_1(1-\tau_2)x[\eta_1+1, \eta_2] \nonumber \\
& + (1-\tau_1)\tau_2x[\eta_1, \eta_2+1] + \tau_1\tau_2x[\eta_1+1, \eta_2+2]
\end{eqnarray}
which is a linear combination of the values of the pixels surrounding the target coordinate obtained by inverting the affine transform. As for all local interpolation schemes, many other methods exist to estimate the value of the subpixel, each with interesting properties and increasing complexity; it should be clear from the more complete analysis conducted in one dimension that in all cases we are dealing with localized approximations to an ideal sinc interpolation, affine transform and resampling.
\section{Fourier Transforms}
\subsection{2D-DFT}
As we have seen, the space of $N_1\times N_2$ images is isomorphic to $\mathbb{C}^{N_1N_2}$ (although we will only deal with real-valued images, so in fact we could consider $\mathbb{R}^{N_1N_2} \subset \mathbb{C}^{N_1N_2}$). Amongst all possible bases for $\mathbb{C}^{N_1N_2}$, the Fourier basis holds a preeminent position and allows us to map an $N_1\times N_2$ signal into a complex $N_1\times N_2$ signal via the 2D-DFT as:
\begin{equation}
X[k_1, k_2] = \sum_{n_1 = 0}^{N_1-1}\sum_{n_2=0}^{N_2-1} x[n_1,n_2] e^{-j\frac{2\pi}{N_1} n_1 k_1} e^{-j\frac{2\pi}{N_2} n_2 k_2}.
\end{equation}
The set of $N_1N_2$ basis vectors
\begin{equation}
w_{N_1,N_2}^{k_1, k_2}[n_1,n_2] = e^{j\frac{2\pi}{N_1} n_1 k_1} e^{j\frac{2\pi}{N_2} n_2 k_2}.
\end{equation}
forms an orthogonal set in $\mathbb{C}^{N_1N_2}$:
\begin{equation}
\sum_{n_1 = 0}^{N_1-1}\sum_{n_2=0}^{N_2-1} w_{N_1,N_2}^{k_1, k_2}[n_1,n_2] w_{N_1,N_2}^{h_1, h_2}[n_1,n_2] = N_1N_2\delta[k_1-h_1, k_2-h_2]
\end{equation}
so that the 2D-DFT can be inverted as
\begin{equation}
x[n_1,n_2] = \frac{1}{N_1 N_2}\sum_{k_1 = 0}^{N_1-1}\sum_{k_2=0}^{N_2-1} X[k_1, k_2] e^{j\frac{2\pi}{N_1} n_1 k_1} e^{j\frac{2\pi}{N_2} n_2 k_2}.
\end{equation}
The basis vectors are planar oscillatory components whose frequency in the vertical and horizontal directions are determined by the indices $k_1$ and $k_2$. Some examples are displayed in Figures~\ref{2D-DFT-1} to~\ref{2D-DFT-3}; the figures plot the magnitude of the basis vectors by mapping the maximum value (1) to white and the minimum value (0) to black. Note that the Fourier basis is separable and that therefore the 2D-DFT can be computed efficiently as the cascade of a row-by-row 1D transform followed by a column-by-column transform. This algorithmic advantage is however somewhat offset by the fact that separability privileges the directions parallel to the vertical and horizontal axes which, in natural images, are not ``special'' directions at all; as a consequence, the basis is somewhat ill-equipped to encode naturally occurring patterns.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 2D-DFT basis vectors
%
% MATLAB code:
%%% Plot 2D-DFT basis vectors.
%%% List here all k1,k2 pairs to plot
%PAIRS = [[1 0]; [2 0]; [3 0]; [30 0]; [0 1]; [0 3]; [1 1]; [250 250]; [2 7]; [3 230]];
%%% image size (square)
%N = 256;
%for n=1:length(PAIRS)
% clf
% k=PAIRS(n,:);
% imagesc(real(exp(1j*(2*pi/N)*(0:N-1)'*k(1))*exp(1j*(2*pi/N)*(0:N-1)*k(2))));
% axis equal
% axis off
% set(gca, 'Position', get(gca, 'Position') - get(gca, 'TightInset') * [-1 0 1 0; 0 -1 0 1; 0 0 1 0; 0 0 0 1]);
% print('-depsc', '-r300', sprintf('2DDFT-%d-%d',k(1),k(2)))
%end
%
%
\begin{figure}\label{2D-DFT-1}
\begin{center}
\begin{tabular}{cc}
\basisImage{1}{0} & \basisImage{2}{0} \\
\basisImage{3}{0} & \basisImage{30}{0}
\end{tabular}
\end{center}
\caption{Some basis vectors for $\mathbb{C}^{256\times 256}$ (real part).}
\end{figure}
\begin{figure}\label{2D-DFT-2}
\begin{center}
\begin{tabular}{cc}
\basisImage{0}{1} & \basisImage{0}{3} \\
\basisImage{1}{1} & \basisImage{5}{250}
\end{tabular}
\end{center}
\caption{Some basis vectors for $\mathbb{C}^{256\times 256}$ (real part).}
\end{figure}
\begin{figure}\label{2D-DFT-3}
\begin{center}
\begin{tabular}{cc}
\basisImage{2}{7} & \basisImage{3}{230} \\
\end{tabular}
\end{center}
\caption{Some basis vectors for $\mathbb{C}^{256\times 256}$ (real part).}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Magnitude and Phase}
Just as the ear is largely insensitive to phase information in sound, the opposite happens with the eye. To understand why, consider that most of the semantic information in a scene can be almost always intelligibly summed up by an outline, i.e. by the edges of the objects contained in the picture. As we have seen several times by now, a sharp transition in a signal is encoded mostly in the phase domain, since the wavefronts of the underlying Fourier basis vectors must align properly in order for their weighted sum to transition sharply. This is immediately apparent if we perform the following experiment: take an image, compute its DFT and then try to reconstruct it either only from the magnitude spectrum or from the phase information. The results are shown in Figure~\ref{MagPhase} where the picture on the left depicts
\[
x_m[n_1,n_2] = \Re\{\mbox{IDFT} \{|X[k_1,k_2]|\}\}
\]
while the image on the right shows
\[
x_p[n_1,n_2] = \Re\{\mbox{IDFT} \{e^{j\angle X[k_1,k_2]}\}\}
\]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\begin{center}
\begin{tabular}{cc}
\showImage{magOnly}%
&
\showImage{phaseOnly}
\end{tabular}
\end{center}
\caption{Image reconstruction from only magnitude and phase}\label{MagPhase}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{2D-DTFT}
For theoretical derivations involving infinite-size images, the DTFT can be extended to the two-dimensional domain as
\begin{equation}
X(\omega_1, \omega_2) = \sum_{n_1 = -\infty}^{\infty}\sum_{n_2 = -\infty}^{\infty} x[n_1,n_2] e^{-j(\omega_1 n_1 + \omega_2 n_2)}
\end{equation}
The above expression, when it exists, maps an infinite two-dimensional discrete signal to a complex function of two real variables, $2\pi$-periodic in both arguments. It can be inverted via the formula
\begin{equation}
x[n_1,n_2] = \frac{1}{(2\pi)^2} \int_{-\pi}^{\pi} \int_{-\pi}^{\pi} X(\omega_1, \omega_2) e^{j(\omega_1 n_1 + \omega_2 n_2)} d\omega_1 d\omega_2.
\end{equation}
\section{Image Filtering}
To begin, a cautionary note: an unfortunate reality in image processing is that very little that can be accomplished by means of linear operators only; a filtering operation is usually just the preparatory step that preludes to some nonlinear and very often adaptive treatment. The reason for this is actually rather easy to understand: an image is a finite-size signal that depicts widely different types of information in different areas. Think of a photo of a person in front of a car in front of a building: the variability in types of surface, edge curvature, textures and so on is extremely large. A filter, by definition, is supposed to be linear and \emph{space-invariant} and therefore only very simple and very ``mild'' filters can be applied to the entirety of the data without negatively affecting one or the other of the several coexisting regions. As a consequence, the design of 2D image filters, although a theoretically rich and challenging subject, remains a much more specialized discipline than its one-dimensional counterpart.
\subsection{Filter Types}
Just as in the one-dimensional case, two-dimensional filters can have a finite or an infinite support, they can be causal or not, and they can be classified in terms of their frequency response. Generally speaking, by eliminating or attenuating the high frequency components, \textbf{lowpass} filters perform a smoothing operation; in the case of images, this correspond to \emph{blurring} the visual content: edges lose their sharpness and small details become less visible. \textbf{Highpass} filters, by contrast, reduce the low-frequency components; in images, this means that uniform areas are eliminated while the fast-varying regions are preserved. Since edges are points of discontinuity in an image and therefore possess a high frequency content, highpass filtering is a key first step in edge detection.
In image processing, in contrast to one-dimensional signal processing, causality is hardly ever a concern. Images are finite-size entities and, as such, we can rely on the fact that they will be available in full when processing begins; the convolution sum can therefore be computed in any chosen direction (row by row, or column by column, right to left or left to right).
In ``everyday'' image processing, IIR filters are not very common because of four main difficulties:
\begin{itemize}
\item \textbf{Linear phase.} In images linear phase is very important since, as we have seen, most of the visual information is contained in edges and a sharp edge require precise phase alignment. As for their 1D counterparts, IIR filter with exact linear phase do not exist and the design of approximately linear phase IIRs can be difficult.
\item \textbf{Border effects.} In one dimension, a finite-length signal has only two boundary points, ie its beginning and its end; as a consequence, even for moderately long signals, border effects can be safely neglected since they affect only a minute percentage of the whole data set. For finite-size images, on the other hand, the number of boundary points is proportional to the square root of the total number of pixels; in this case, border effects become much more significant; one way to combat the problem is to use FIR filters with a small support so that only a few pixels around the border are affected.
\item \textbf{Stability Issues.} Multidimensional CCDE's are more complex entities than their one-dimensional counterparts; multivariate polynomials do not factor into first-order terms so that a multidimensional $z$-transform representation does not lead to either a convenient stability test via root checking nor to simple inversion formulas via partial fraction expansion.
\item \textbf{Computability.} An arbitrary multidimensional CCDE may not be computable: consider for instance the filter defined by
\begin{eqnarray*}
y[n_1, n_2] = a_0y[n_1+1, n_2] + a_1y[n_1, n_2-1] + \\
a_2y[n_1-1, n_2] a_1y[n_1, n_2+1] + x[n_1, n_2];
\end{eqnarray*}
each output sample depends on four neighboring output samples (see Figure~\ref{nonComputable}) which, in turn, depend on the current output. Contrary to the one-dimensional case, in which a finite-order recursive CCDE can always be rewritten in causal form, here it is impossible to rearrange the equation to provide a computable algorithm.
\end{itemize}
For these reasons, most filters used in practice are small-support FIRs while the design and implementation of IIRs remains the province of more specialized papers and applications.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\begin{center}
\begin{dspCP}[xout=true]{-1.5,1.5}{-1.5,1.5}
\footnotesize
\dspPointValue[dotstyle=o]{0,0}{$y[0,0]$}
\dspPointValue{1,0}{$y[1,0]$}
\dspPointValue{0,1}{$y[0,1]$}
\dspPointValue{-1,0}{$y[-1,0]$}
\dspPointValue{0,-1}{$y[0,-1]$}
\end{dspCP}
\end{center}
\caption{Noncomputable 2D IIR: computing the value in $(0, 0)$ depends on four neighboring output values. \label{nonComputable}}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Lowpass Filtering}
Visually, lowpass filtering blurs the details in an image; while many lowpass characteristics can be used to perform this task, we will simply look at two alternatives, the two-dimensional moving average and the Gaussian blur. The latter is particularly interesting because its frequency response replicates the transfer function of a defocused camera lens.
\itempar{Moving Average.} The straightforward 2D extension of the moving average filter is the separable FIR
\[
h[n_1,n_2] = \frac{n_1}{(N+1)^2}\rect\left(\frac{n_1}{2N}, \frac{n_2}{2N}\right)
\]
which has a square support; again, in image processing causality is rarely a concern so that we can express the filter as centered around the origin. Each output pixel $y[n_1,n_2] = x[n_1,n_2] \ast h[n_1,n_2] = $ by averaging the original image over a square of side $N+1$ pixels centered around $(n_1, n_2)$. Examples of MA filtering are shown in Figure~\ref{MAfir}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\begin{center}
\begin{tabular}{cc}
\showImage{MA11}
&
\showImage{MA51}
\end{tabular}
\end{center}
\caption{Moving average blurs with a $11\times 11$ and a $51\times 51$ FIR.}\label{MAfir}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\itempar{Gaussian Filter.} \label{blurring} A preferred lowpass alternative to the moving average is the separable Gaussian filter
\[
h[n_1, n_2] = \frac{1}{2\pi\sigma^2}e^{-\frac{n_1^2+n_2^2}{2\sigma^2}}
\]
which has a smoother frequency characteristic (the Fourier transform of an Gaussian function is itself a Gaussian). Although a Gaussian response is nonzero over the entire real plane, a common approach is use a square-support FIR approximation with a width of approximately $3\sigma$ from the mean, which ensures that over 99\% of the total filter's energy is represented. For example, if $\sigma = 1$ we would use the $5\times 5$ FIR with values
\[
h[n_1, n_2] =
\begin{bmatrix}
1 & 4 & 7 & 4 & 1 \\
4 & 16 & 26 & 16 & 4 \\
7 & 26 & 41 & 26 & 7 \\
4 & 16 & 26 & 16 & 4 \\
1 & 4 & 7 & 4 & 1 \\
\end{bmatrix}
\cdot \frac{1}{273}
\]
where the taps have been normalized so that $\sum_{n_1}\sum_{n_2} h[n_1, n_2] = 1$. Note how, in the case of FIR filters, we can specify the impulse response as a $N_1\times N_2$ matrix with the assumption (unless otherwise noted) that the center value of the matrix coincides with the coordinate origin. Some examples of Gaussian blur are shown in Figures~\ref{GB1} and~\ref{GB2}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\begin{center}
\begin{tabular}{cc}
\showImage{GB1}
&
\showImage{GB2}
\end{tabular}
\end{center}
\caption{Gaussian blur: $\sigma = 1$ ($5\times 5$ FIR) and $\sigma = 1.8$ ($9\times 9$ FIR).}\label{GB1}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\begin{center}
\begin{tabular}{cc}
\begin{dspCP}[xout=true]{-14,14}{-14,14}%
\dspImageFile{Gauss5}%
\end{dspCP}
&
\begin{dspCP}[xout=true,xticks=custom,yticks=none]{0,28}{0,28}%
\dspImageFile{GB5}%
\end{dspCP}
\end{tabular}
\end{center}
\caption{Gaussian blur: filter's impulse response for $\sigma = 5$ ($29\times 29$ FIR) and filtering result.}\label{GB2}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Higpass Filtering and Edge Detection}
As we have seen, edges carry a significant portion of the semantic information contained in an image and therefore edge detection is an important type of operation in image processing. Mathematically, however, the concept of ``edge'' is rather difficult to define unambiguously; what appears to a human observer as an obvious boundary between objects may not correspond to a clearly identifiable pattern in the underlying intensity level of the corresponding pixels. A typical example is the horizon line between sea and sky, which we can easily identify thanks to accessory semantic cues, but which is hard to determine solely as a function of the intensity gradient. Broadly speaking, we can have two types of edges:
\begin{itemize}
\item a ``hard'' edges is a set of points across which the intensity changes abruptly, as in figure~\ref{fig:ip:edges}-(a); these edges are usually detected by searching for points of large magnitude for the gradient of the signa;
\item a ``soft'' edge is a set of points that correspond to a change of inflection in the intensity, as in figure~\ref{fig:ip:edges}-(b); these edges are usually detected by looking for the zero crossings of the Laplacian (that is, ponts where the sum of the second-order derivatives of the signal is zero).
\end{itemize}
Edge detection algorithms start by determining a set of candidate edge points using one of the above criteria and follow with a series of thresholding and heuristics to refine the set of edges. In both cases, as we will see, the initial filtering operation approximates a differentiation and, therefore, it implements a highpass filter.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\begin{center}
\begin{tabular}{cc}
\begin{dspCP}[xout=true]{-14,14}{-14,14}%
\dspImageFile{bw}%
\end{dspCP}
&
\begin{dspCP}[xout=true,xticks=custom,yticks=none]{0,28}{0,28}%
\dspImageFile{2DDFT-1-0}%
\end{dspCP}
\end{tabular}
\end{center}
\caption{Hard and soft edges.}\label{fig:ip:}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\itempar{Sobel Operators.} The Sobel filter provides a set of candidate edge points by finding the points where the derivative of the intensity is large. For a function of two variables $f(t_1, t_2)$ the gradient vector
\begin{equation}
\nabla f = \begin{bmatrix}
\displaystyle\frac{\partial f}{\partial t_1} &
\displaystyle\frac{\partial f}{\partial t_2}
\end{bmatrix}^T
\end{equation}
points in the direction of the maximum rate of change of the function around a given point. For a discrete-space image, the gradient can be approximated by taking the accumulated contributions of the first differences between neighboring points around a target location. The idea, due to Sobel, works as follows.
Consider the 8~closest neighbors to a target pixel $x$, labeled as in the following pattern:
\begin{equation}
\begin{bmatrix}
x_1 & x_2 & x_3 \\
x_8 & x & x_4 \\
x_7 & x_6 & x_5 \\
\end{bmatrix}
\end{equation}
We can approximate the magnitude of the first derivative around $x$ in eight possible directions as the difference between $x$ and each of its chosen neighbor, divided by the distance $d(x, x_i)$ between the two pixels. By summing all possible contributions, the gradient can be approximated as
\begin{equation} \label{eq:ip:gradapp}
\nabla x \approx \sum_{i=1}^{8} \frac{x-x_i}{d(x, x_i)}\mathbf{x}_i
\end{equation}
where $\mathbf{x}_i$ is a unit-norm vector that indicates the direction of each contribution, that is:
\begin{align*}
\mathbf{x}_1 &= \begin{bmatrix}-1 & 1\end{bmatrix}^T \\
\mathbf{x}_2 &= \begin{bmatrix}0 & 1\end{bmatrix}^T \\
\mathbf{x}_3 &= \begin{bmatrix}1 & 1\end{bmatrix}^T \\
\ldots \\
\mathbf{x}_8 &= \begin{bmatrix}-1 & 0\end{bmatrix}^T
\end{align*}
Without loss of generality, set the minimum distance between two pixels equal to two, so that
\[
d(x, x_i) = \begin{cases}
2 & \mbox{for $i$ even} \\
4 & \mbox{for $i$ odd}
\end{cases};
\]
this will just introduce a scaling factor of four in the estimation of the gradient's magnitude but will keep the computation simpler. By replacing the expressions for $d(x, x_i)$ and $\mathbf{x}_i$ in~(\ref{eq:ip:gradapp}) and collecting terms, we have
\begin{equation}
4\nabla x \approx \begin{bmatrix}
x_3 - x_1 +2(x_4-x_8) +x_5 - x_7 \\ x_7-x_1 +2(x_6-x_2) + x_5-x_3
\end{bmatrix}
\end{equation}
The first and second component of the gradient vector can thus be obtained as a linear combination of neighboring pixels using the \emph{horizontal} and \emph{vertical} Sobel FIR filters with impulse responses
\begin{equation} \label{eq:ip:horsobel}
s_h[n_1, n_2] =
\begin{bmatrix}
-1 & 0 & 1 \\
-2 & 0 & 2 \\
-1 & 0 & 1 \\
\end{bmatrix}.
\end{equation}
and
\begin{equation} \label{eq:ip:versobel}
s_v[n_1, n_2] =
\begin{bmatrix}
-1 & -2 & 1 \\
0 & 0 & 0 \\
1 & 2 & 1 \\
\end{bmatrix}.
\end{equation}
Interestingly, the impulse responses are separable; for the horizontal filter, for instance, we have
\[
s_h[n_1, n_2] = \begin{bmatrix}
1 \\
2 \\
1 \\
\end{bmatrix}
\begin{bmatrix} -1 & 0 & 1 \end{bmatrix};
\]
this decomposition of the impulse response shows that the filter can be interpreted as the cascade of a vertical three-point weighted average and a horizontal three-point truncation of the impulse response of an ideal discrete-time differentiator (see Secion 9.2). The vertical averaging prior to differentiation ensures that the operator is less sensitive to noise while still responsive to magnitude jumps.
Using the output of both filters, the magnitude of the gradient of an image $x[n_1, n_2]$ can be estimated via the nonlinear operator
\begin{equation}
% |\nabla x[n_1, n_2]|^2 = |s_h[n_1, n_2] \ast x[n_1, n_2]|^2 + |s_v[n_1, n_2]\ast x[n_1, n_2]|^2
|\nabla x|^2 = |s_h \ast x|^2 + |s_v\ast x|^2;
\end{equation}
since the gradient is large across jumps in signal magnitude, its local extrema can be used to discover edges in the image. An example of the gradient estimate obtained via the Sobel operator is shown on the left of Figure~\ref{SOBEL}; a common approach to edge detection would begin by thresholding this image as
\begin{equation}
y[n_1, n_2] = \begin{cases}
0 & \mbox{if $|\nabla x[n_1, n_2]|^2 > T$} \\
1 & \mbox{otherwise}
\end{cases}
\end{equation}
which yields the picture on the right.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\begin{center}
\begin{tabular}{cc}
\showImage{sobel}
&
\showImage{sobelTh}
\end{tabular}
\end{center}
\caption{Sobel gradient estimation: gradient image and thresholded gradient ($T=80$).}\label{SOBEL}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\itempar{Laplace Operator.} This operator provides a set of candidate edge points by finding the changes in inflection for the intensity values; these correspond to the values that minimize the Laplacian of the signal. For functions of two continuous variables, the Laplacian is defined as the sum of the second partial derivatives:
\[
\Delta f(t_1, t_2) = \frac{\partial^2 f}{\partial t_1^2} + \frac{\partial^2 f}{\partial t_2^2};
\]
the zero crossings of the Laplacian correspond to points of
In the discrete-space case, we can approximate the Laplacian operator as follows. First, in order to derive a compact approximation to the second partial derivative, consider the Taylor series expansion of a function of a real variable:
\[
f(t+\tau) = \sum_{n=0}^{\infty} \frac{f^{(n)}(t)}{n!}\tau^n
\]
If we write out the expansions for $f(t+\tau)$ and $f(t-\tau)$, stopping at the second order term, we have
\begin{eqnarray*}
f(x+\tau) &= f(t) + f'(t)\tau + \frac{1}{2}f''(t)\tau^2 \\
f(x-\tau) &= f(t) - f'(t)\tau + \frac{1}{2}f''(t)\tau^2
\end{eqnarray*}
the sum of which yields the approximation
\[
f''(t) = \frac{f(t-\tau) - 2f(t) + f(t+\tau)}{\tau^2}.
\]
On the discrete-space grid $\tau=1$ and we can therefore approximate the second derivative with the filter response $\begin{bmatrix}1 & -2 & 1 \end{bmatrix}$. By summing the horizontal and vertical approximations we finally obtain the nonseparable Laplacian filter:
\[
h[n_1, n_2] =
\begin{bmatrix}
0 & 1 & 0 \\
1 & -4 & 1 \\
0 & 1 & 0 \\
\end{bmatrix}.
\]
A common variation, which improves the rotational symmetry of the operator, includes the approximation of the derivatives along the 45-degree diagonals:
\[
h[n_1, n_2] =
\begin{bmatrix}
1 & 1 & 1 \\
1 & -8 & 1 \\
1 & 1 & 1 \\
\end{bmatrix}.
\]
The Laplacian is another useful tool in edge detection; indeed, since edges correspond to sharp amplitude transition, the ``curvature'' of the image signal will exhibit sharp positive and negative peaks around the edges. An example of image processing using the Laplacian filter is shown in Figure~\ref{LAPL}; the first figure shows the absolute value of the Laplacian, while the second figure shows the image obtained by thresholding the first. The figure clearly shows that the Laplacian operator is particularly sensitive to noise, which results in a lot of disconnected pixels in the thresholded image. A way to combat this problem is to precede the Laplacian operator by a Gaussian blur.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\begin{center}
\begin{tabular}{cc}
\showImage{lapl}
&
\showImage{laplTh}
\end{tabular}
\end{center}
\caption{Laplacian operator: magnitude and thresholded magnitude ($T=50$).}\label{LAPL}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c4science · Help