This preview shows page 1. Sign up to view the full content.
Unformatted text preview: Medical Image Processing 1 (MIP1)
Image Signal Analysis and Processing (ISAP)
SS’09
Linear Filters
Dr. Pierre Elbischger
School of Medical Information Technology
Carinthia University of Applied Sciences Linear shift invariant system (LSI)
Most imaging systems have for the most part three significant properties
H
LSISystem f(x) H(f(x)) Linearity
Superpostion: The response to the
sum of stimuli is the sum of the
individual responses.
Scaling: The response to the scaled
stimuli is the response scaled by the
same sacling factor to the unscaled
stimuli .
Shift invariance: The response to a translated stimulus is just a translation of the response
to the stimulus. A device that is linear and shift invariant is known as a shift invariant linear system, or often
just as a system. The response of a shift invariant linear system to a stimulus is obtained by a
convolution of the input signal with the systems’ impulse response function.
12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 2 Test whether a system is a LSI system
System under test Linearity
The input is the superposition of two different and scaled functions Scaling an superposition is verified. The system is linear.
Shift invariance
The shifted output signal is given by
This is not the same as the system applied to the shifted input
The system is not shift invariant
12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 3 The Dirac delta function
The Dirac delta function can be loosely thought of as a function on the real line
which is zero everywhere except at the origin, where it is infinite. The area below the
function is constrained to be one. The delta function can be viewed as the limit of a so called nascent delta function Two nascent delta functions are: Rectangular function
12.05.2009 Normal distribution
Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 4 2D Dirac delta function
The ideal impulse in the image plane is defined by Sifting property
The convolution of the Dirac delta function with the twodimension function f(x,y)
provides the value of the function f(x,y) at the point λ,µ. 12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 5 2D convolution
Convolution is a very useful linear, translationinvariant operation. The convolution expresses a
linear filtering process using the linear shift invariant system defined by its impulse function
(filter) h(x,y) that is the output of the system given a Dirac impulse as the input signal. The animations graphically illustrate the 1D convolution of two rectangle functions (left) and
two Gaussians (right). In the plots, the green curve shows the convolution of the blue and red
curves as a function of t. The gray region indicates the product of both functions as a function
of t, so its area as a function of t is precisely the convolution. 12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 6 2D Convolution 12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 7 Discrete 2D LSISystem – Discrete convolution
In the general case the output signal of a linear system is the weighted superposition of
the system’s impulse response function h(n,m) with weights representing the input
signal x(k,l). In shift variant systems the impulse response is a function of translation. x(n,m) h(k,l) y(n,m) In a shift invariant system a translation of the input results in a translation of the output. The output can now be described by the discrete convolution of x(k,l) with h(n,m). 12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 8 Point spread function (PSF)
Point spread function (PSF)
Each infinitesimal point of an object
causes a blurred point in the image. If the
optical system can be assumed to be
linear, then the blurred point is represented
by the impulse response function h(x,y)
of the system
– the PSF. Optical transfer function (OTF)
The Fourier transform of the PSF is the so called OTF. 12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 9 Discrete convolution (filtering)  Illustration
y origin h(1,1) h(1,0) h(1,1) h(0,1) h(0,0) h(0,1) h(1,1) mask h(1,0) h(1,1) f(x1,y1) f(x1,y) f(x1,y+1) Mask coefficients x
f(x,y1) f(x,y) f(x,y+1) f(x+1,y1) f(x+1,y) f(x+1,y+1)
Image section
under mask
12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 10 Considerations about the border
The convolution of an NxN image with an LxL impulse response function results in no
meaningful data at the (L1)/2 border pixels of the (full size MxM) output image.
Problem
At the borders the filter
kernel is multiplied with no
well defined image pixels.
(same size) Solutions
1. Reject the border pixels →
The image becomes smaller
2. User zeropadding
3. Continue the image at its
borders using reflected
copies of the image (full size) (valid size) Matlab: imfilter, conv2
12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 11 Frequency response of a 2D LSISystem
The complex exponential is an eigenfunction of the LSISystem. Thus, it passes
the system without changing its basic shape but only its amplitude and phase is
affected System’s frequency response. The frequency response is a periodic function of ωx and ωy with period 2π. Matlab: freqz2
12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 12 Separable filter kernels
A filter kernel is called separable if The 2D convolution of an image with a separable digital filter can be implemented as the
cascade of two 1D convolutions. This significantly reduces the
computational complexity.
The table lists the number of
mathematical operations that are
required to compute the 2D
convolution at a singe location
(pixel) using a filter kernel of size
LxL. In order to obtain the full
convolution result, the kernel has
to be evaluated at NxN (output
image has size of the input image)
locations.
12.05.2009 # mathematical
operations 2D convolution 2 x 1Dconvolution L2 2L Additions (L1)2 2(L1) Locations N2 N2 Multiplications Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 13 Filtering in the frequency domain
• Make use of the FFT that is a fast implementation of the DFT. • If the dimension of the filter kernel L is sufficiently large with respect to the
dimension of the input array N, Fourier domain convolution will be more efficient
than direct convolution, perhaps by an order of magnitude or more.
In many signal processing applications, the same impulse response is used on
different data. Hence, the Fourier transform of the impulse response can be precomputed and, therefore, actually does not cause computational costs.
H(u,v) must be of the same size as F(u,v) • • 12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 14 Zeropadding (circular convolution)
• The periodicity of the image f(n,m) that is implicitly introduced by the DFT may
lead to a wrong convolution result at the boundaries due to circular convolution. • Use zeropadding of the image to obtain correct convolution results Matlab: padarray
12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 15 Illustration of zeropadding effects
Proper zeropadding Improper zeropadding Filter Input image Output image 12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 16 T10 Filtering in the frequency domain
Frequency domain convolution
Operation
Forward and backward FFT
Scalar multiplication
Sum
Efficiency of the FFT in the 1Dcase Direct convolution
requires on the order of
J2L2 multiplications Order
2 J2 ld(J)
J2
J2 (1 + 2 ld(J))
perform the
computation
consecutively along
the column and
row dimension do for all vectors
in the dimension
(e.g. along all
columns) 1DFFT
along a single
dimension 2⋅ J ⋅ J⋅ld(J)
The 2D Fourier Transforma
can be computed by
consecutively applying the
onedimensional
FFT along the rows and
columns of the image, see
section ‘Fourier Transform’.
12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters Matlab: imfilter, fft2 17 Filters properties (1) – Linear phase
•
• Use filters of linear phase (symmetric or skew symmetric h(m,n)) leads to a full
sample shift of the output image that can be corrected no shift of image strucutres!
An odd number of filter taps lead to a zero crossing of the magnitude of H(π,π) 12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 18 Filters properties (2)
Preservation of the mean value for lowpass filtering Monotonically decreasing transfer function for lowpass filtering
This is especially of interest for scale spaces, where small structures should vanish before
larger ones. Isotropy
The smoothing characteristic of a typical smoothing filters should not depend on the direction,
thus the impulse function and the transfer function are circularly symmetric. Matlab: fspecial … various standard filter kernels
12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 19 Ideal lowpass filter •
•
•
•
• Sharp cutoff can physically not be
realized
Causes strong ringing effects
Isotropic
Not separable
Very long impulse response (large
filter kernel) 12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 20 Smoothing filter – Uniform kernel (rectangular) Frequency response function. The ringing
depends on the number of used filter taps.
anisotropic
12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 21 Butterworth filter Increasing n
•
•
•
• Smoothness of the cutoff depends on the
filter order n
Ringing increases with filter order
Isotropic
Not separable 12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 22 Butterworth filter – Isotropic 12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 23 Gaussian filter Increasing σ •
•
•
• Isotropic
Separable
A stronger smoothing does not lead to a stronger ringing
A Gaussian convolved with a Gaussian results in another Gaussian 12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 24 Scale space – Gaussian pyramid 1. Set image to the finest scale
2. For each layer from the finest to coarsest 1. Obtain this layer by smoothing the next finest
2. Then subsample the smoothed image Matlab: mA(1:2:end, 1:2:end) … subsampling of array mA by a factor of 2
12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 25 Bandreject filter 12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 26 Gradient filters g(x,y) edge direction For faster but less accurate computation: Matlab: atan2, abs
12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 27 Discrete approximation of the gradient T11 E1 The ideal gradient filter The first order approximation of the gradient filter The left and right discrete derivative
impulse response function of the
derivative operator The symmetric derivative 12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 28 Sobel filter E2 • The derivative operator amplifies the high frequency content and, therefore, is
sensitive to noise. • Smooth in both dimensions and compute the derivative in one dimension • The filter is symmetric and, therefore, has linear (zero) phase. symmetric derivative
Normalization of the smoothing
kernels (preserve the mean value,
see a couple of slides before) lead
to a scaling factor of 1/8, that have
to be used, if the correct value of
the derivative is of interest –
assuming a sample spacing of one. First order difference:
Smoothing kernels: 12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 29 Gradient images – Symmetric derivatives Original image 12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 30 Derivative of Gaussian
The Gaussian kernel is defined by: The partial derivative in direction x is computed on the image f(x,y) that has
been smoothed using a Gaussian kernel Kσ(x,y). This operation is equivalent to
the convolution of the image with the derivative of the Gaussian kernel. 12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters Gaussian derivative 31 Isotropy of the symmetric difference filter
+90° phase shift (differentiator) Error of the gradient as a function of the spatial
frequency and the orientation respectively. The maximum error
occurs at 22.5° Error in magnitude 12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters Error in the angle 32 Optimized filter kernels
The separable kernel vectors d0 and d1 of the derivative kernel can be obtained by
solving an optimization problem, where the objective function is the difference
between the ideal transfer function of the derivative and the transfer function of the
kernel. Column vectors for optimized 2D separable derivative kernels 12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 33 Angle error of different gradient kernels E3, E4 Magnitude of the angle error in degree for the direction of the largest error (22.5°).
Upper line: 3 x 3 filters (a symmetric difference, b Sobel, c optimized).
Lower line: 5 x 5 filters (d [1 8 0 8 1]/16, e derivative of Gauss, f optimized) 12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 34 Laplacian of Gaussian (LoG) T12 Use the Laplacian (2. derivative) of the analytic Gaussian function G(x,y) to compute
the second derivative. A discrete LoG filter for some particular σ.
0
0 hLoG (n, m) = 1 0
0 01
12
2 −16
12
01 0
1
2
1
0 0
0 1 0
0 1 1 1
hLoG (n, = 1 −8 1
m) 1 1 1 • The coefficients sum to zero
• In homogeneous regions the filter response is small: it becomes large at locations of a large contrast.
• Sometimes the negative of the filter mask is used.
LoG  Mexican hat
12.05.2009 Dr. Pierre Elbischger  MIP1'WS08  Linear Filters 35 ...
View Full
Document
 Spring '09
 PIERREELSCHBINGER

Click to edit the document details