mip1_04_linear_filters_090512_1455219

mip1_04_linear_filters_090512_1455219 - Medical Image...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

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 LSI-System 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 two-dimension 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, translation-invariant 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 LSI-System – 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(x-1,y-1) f(x-1,y) f(x-1,y+1) Mask coefficients x f(x,y-1) f(x,y) f(x,y+1) f(x+1,y-1) 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 (L-1)/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 zero-padding 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 LSI-System The complex exponential is an eigenfunction of the LSI-System. 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 1D-convolution L2 2L Additions (L-1)2 2(L-1) 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 Zero-padding (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 zero-padding of the image to obtain correct convolution results Matlab: padarray 12.05.2009 Dr. Pierre Elbischger - MIP1'WS08 - Linear Filters 15 Illustration of zero-padding effects Proper zero-padding Improper zero-padding 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 1D-case 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) 1D-FFT along a single dimension 2⋅ J ⋅ J⋅ld(J) The 2D Fourier Transforma can be computed by consecutively applying the one-dimensional 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 low-pass filtering Monotonically decreasing transfer function for low-pass 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 low-pass filter • • • • • Sharp cut-off 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 cut-off 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 Band-reject 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

Ask a homework question - tutors are online