Slide 1Slide 2Slide 3Slide 4Slide 5Slide 6Slide 7Slide 8Slide 9Slide 10Slide 11Slide 12Slide 13Slide 14Slide 15Slide 16Slide 17Slide 18Slide 19Slide 20Slide 21Slide 221Lucas Parra, CCNYCity College of New YorkBME I5000: Biomedical ImagingLecture 11Point Spread Function, Inverse Filtering, Wiener Filtering, Sharpening, ... Lucas C. Parra, [email protected] Blackboard: http://cityonline.ccny.cuny.edu/2Lucas Parra, CCNYCity College of New YorkSchedule1. Introduction, Spatial Resolution, Intensity Resolution, Noise2. X-Ray Imaging, Mammography, Angiography, Fluoroscopy3. Intensity manipulations: Contrast Enhancement, Histogram Equalisation4. Computed Tomography5. Image Reconstruction, Radon & Fourier Transform, Filtered Back Projection6. Nuclear Imaging, PET and SPECT7. Maximum Likelihood Reconstruction8. Magnetic Resonance Imaging9. Fourier reconstruction, k-space, frequency and phase encoding 10. Optical imaging, Fluorescence, Microscopy, Confocal Imaging11. Enhancement: Point Spread Function, Filtering, Sharpening, Wiener filter12. Segmentation: Thresholding, Matched filter, Morphological operations13. Pattern Recognition: Feature extraction, PCA, Wavelets14. Pattern Recognition: Bayesian Inference, Linear classification3Lucas Parra, CCNYCity College of New YorkModel for a simple imaging system ...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): Model Imaging Systemsourceimageimagedigitized imageLinear Space Invariant systemSampling anddigitization h (x)=LSI [δ (x )]xδ( x)h (x)4Lucas Parra, CCNYCity College of New YorkA common simplified mathematical model for an imaging process is that of a linear system with additive noise: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)LSIh(x,y)s(x, y)+n(x, y)g(x, y)Model Imaging System5Lucas Parra, CCNYCity College of New YorkPoint Spread FunctionThe image of an arbitrary source s(x) is then given by a convolutionFor a 2D discrete array (image) we write the convolutions asg ( 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 )g ( x , y)=∑x '=1N∑y '=1Mh(x− x' , y− y' )s ( x' , y ')=h( x , y)∗s( x , y)6Lucas Parra, CCNYCity College of New YorkPoint Spread Function>> g = conv2(h,s);7Lucas Parra, CCNYCity College of New YorkPSF – Smoothing, Sharpening1 2 12 4 21 2 1hLP(x,y) = 1/16 ×0 1 01 -4 10 1 0hHP(x,y) =-1 -1 -1-1 8 -1-1 -1 -1hHP(x,y) =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 noiseSharpening: Simple high pass filter that enhances edges, equivalent to a second derivative (Laplacian filter) 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: hhp(x , y )=δ( x , y)−hlp( x , y)8Lucas Parra, CCNYCity College of New YorkPSF – Smoothing, Sharpening9Lucas Parra, CCNYCity College of New YorkAssignment – 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().10Lucas Parra, CCNYCity College of New YorkIn k-space (Fourier Domain) this can be written as 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:Demonstrate: Low-pass filter, high-pass filter, band-pass filter, Laplacian, un-sharp masking, etc. in MATLAB.Point Spread Function - K-spaceG(kx, ky)= H (kx, ky) S (kx, ky)+ N (kx, ky)G(kx, ky)= H (kx, ky) S (kx, ky)11Lucas Parra, CCNYCity College of New YorkEffect of H on S in k-space:k-space12Lucas Parra, CCNYCity College of New YorkIn 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:That is by convolving with the inverse FT of H-1: Inverse FilteringS (kx, ky)=1H (kx, ky)G(kx, ky)s( x , y)=FT−1[1H (kx, ky)]∗g (x , y )13Lucas Parra, CCNYCity College of New YorkProblem: 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 FilteringAssignment 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. Inverse Filtering14Lucas Parra, CCNYCity College of New YorkWiener filter considers the effect of noise N and the magnitude of H . It give the optimal estimate as:σ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.Wiener FilteringS (kx, ky)=σS2H*∣H∣2σS2+σN2G(kx, ky)S (kx, ky)15Lucas Parra,