Lucas Parra CCNY City College of New York BME I5000 Biomedical Imaging Lecture 11 Point Spread Function Inverse Filtering Wiener Filtering Sharpening Lucas C Parra parra ccny cuny edu Blackboard http cityonline ccny cuny edu 1 Lucas Parra CCNY City College of New York Schedule 1 Introduction Spatial Resolution Intensity Resolution Noise 2 X Ray Imaging Mammography Angiography Fluoroscopy 3 Intensity manipulations Contrast Enhancement Histogram Equalisation 4 Computed Tomography 5 Image Reconstruction Radon Fourier Transform Filtered Back Projection 6 Nuclear Imaging PET and SPECT 7 Maximum Likelihood Reconstruction 8 Magnetic Resonance Imaging 9 Fourier reconstruction k space frequency and phase encoding 10 Optical imaging Fluorescence Microscopy Confocal Imaging 11 Enhancement Point Spread Function Filtering Sharpening Wiener filter 12 Segmentation Thresholding Matched filter Morphological operations 13 Pattern Recognition Feature extraction PCA Wavelets 14 Pattern Recognition Bayesian Inference Linear classification 2 Lucas Parra CCNY City College of New York Model Imaging System Model for a simple imaging system Linear Space Invariant system x Sampling and digitization h x source image image x digitized image often assumes linear shift invariance LSI The image of a point like source s x x on the detector combining all blurring effects of the imaging process is called the Points Spread Function h x h x LSI x 3 Lucas Parra CCNY City College of New York Model Imaging System A common simplified mathematical model for an imaging process is that of a linear system with additive noise s x y LSI h x y g x y n x y The source image s x y passes through a Linear Sift Invariant transformation h x y and sensing generates additive noise n i j The linear transformation is given by a convolution g x y h x y s x y n x y 4 Lucas Parra CCNY City College of New York Point Spread Function The image of an arbitrary source s x is then given by a convolution g x LSI s x LSI dx x x s x dx s x LSI x x dx s x h x x dx h x s x x h x s x For a 2D discrete array image we write the convolutions as N g x y M h x x y y s x y x 1 y 1 h x y s x y 5 Lucas Parra CCNY City College of New York Point Spread Function g conv2 h s 6 Lucas Parra CCNY City College of New York PSF Smoothing Sharpening There are a few simple choices of PSF we can apply to an image to improve image quality Smoothing Simple low pass filter to remover high frequency noise 1 2 1 hLP x y 1 16 2 4 2 1 2 1 Sharpening Simple high pass filter that enhances edges equivalent to a second derivative Laplacian filter 0 1 0 1 4 1 0 1 0 1 1 1 1 8 1 1 1 1 hHP x y hHP x y or Often hLP is implemented with a 2D Gaussian PSF and hHP with its second derivative This way it is easier to control the scale Un sharp masking h hp x y x y hlp x y 7 Lucas Parra CCNY City College of New York PSF Smoothing Sharpening 8 Lucas Parra CCNY City College of New York Assignment simple edge detection Use an image with a single item in the image that has a reasonable well defined border say of a cell in a dish Your goal is to measure the circumference Add 6dB noise to the image Low pass the image to reduce the noise Use a high pass filter second derivative to highlight the edge and threshold Select the connected component that corresponds to the boundary of the object Hopefully it is singly connected and you can determine the length of the edge The goal is not to get this perfectly but learn the effects of the different parameters low pass filter size threshold value region size etc Feel free to start without any noise and see if your algorithms still works reasonably well with the added noise Compare this to Canny and Sobel edge detection see matlab function edge 9 Lucas Parra CCNY City College of New York Point Spread Function K space In k space Fourier Domain this can be written as G k x k y H k x k y S k x k y N k x k y where G H S and N are the 2D Fourier Transforms FT of g h s and n respectively Often one considers the ideal noise free case N 0 G k x k y H k x k y S k x k y Demonstrate Low pass filter high pass filter band pass filter Laplacian un sharp masking etc in MATLAB 10 Lucas Parra CCNY City College of New York k space Effect of H on S in k space 11 Lucas Parra CCNY City College of New York Inverse Filtering In the case of zero noise N 0 and a known PSF we can undo the effect of H and recover the original image with an inverse filter 1 S k x k y G k x k y H k x k y That is by convolving with the inverse FT of H 1 s x y FT 1 1 g x y H k x k y 12 Lucas Parra CCNY City College of New York Inverse Filtering Problem H kx ky may be zero or very small for some frequencies At those frequencies even small noise will be stronger than the signal and 1 H kx ky will primarily enhance the noise Solution Wiener Filtering Assignment 10 1 Filter convolve an image with an PSF of your choice 2 Compute the inverse filter using the DFT 3 Recover the original image by filtering with this inverse filter 4 Show all three images your filter and its inverse filter and the difference image between original and recovered image 5 Corrupt the filtered image with additive noise 10dB SNR and repeat step 3 and 4 6 Now repeat the same convolution noise and recover the original image using a Wiener filter 13 Lucas Parra CCNY City College of New York Wiener Filtering Wiener filter considers the effect of noise N and the magnitude of H It give the optimal estimate S k x k y as S k x k y 2 S H G k x k y 2 2 2 H S N 2S E S 2 2N E N 2 are the power spectra of s and n Expected value of the absolute square of the FT H 2 is the magnitude response of h Absolute square of the FT Dependency on x y and kx ky is omitted here to simplify notation This estimate is optimal in that it represents the maximum aposteriory MAP estimate assuming zero mean Gaussian distributed spectra 14 Lucas Parra CCNY City College of New York Prob Estimation Maximum A Posteriori review Sometimes prior information on the model parameters is available in the form of a prior probability density p S S This …