Blurred Blurring Model Figure 2- Kernel Convolution

Blurred Image Restoration using Blind Deconvolution, Wiener Filter and Lucy-Richardson algorithm Abdulrahman Ali Almarhabi Abstract: Image restoration have great importance nowadays. Image processing becomes difficult when image comes to blurred and noisy. In order to improve the quality of the image, there are variety of techniques and algorithms used for image recovering or de-blurring to make pictures sharp and useful such as Lucy-Richardson algorithm, blind deconvolution algorithm, Van Cittert algorithm Wiener filter, regularized filter and sharpening filter.Keywords: Image De-blurring, Blur, PSF, Deconvolution, Wiener, Lucy-Richardson, SSIM, PSNR. 1. INTRODUCTION Blurring of image is common problem while taking picture. It due to certain technical errors and shooting situations. Image restoration is the process of recovering the original image from the degraded image. In this study, I will use Matlab with Image Processing Toolbox to simulate the blurring of images using motion blur, Gaussian blur and box blur.The blurred images will restored using blind deconvolution algorithm, Lucy- Richardson algorithm and Wiener filter. The recovered image will compared based on types of blur, Peak Signal to Noise Ratio (PSNR), Mean Square Error (MSE) and Structure Similarity Index Map (SSIM). 1.1 Blurring 1.1.1Reason of Blurring There are two main causes for blurring of images: – Motionblur: which caused by either object movement or camera shake where the shutter speed dose not fast enough to freeze the action. – Out of focus or focusing on wrong object. 1.1.2 Blurring Model The blur defined as degradation of sharpness and contrast of the image.Mathematically it represents by convolve sharp or true image with certain kernel or point spread function (PSF) as shown in figure 1. ?? ( ?? , ?? ) = ?? ( ?? , ?? ) ? ?? ( ?? , ?? ) + ?? ( ?? , ?? ) . …. (1) Where: g(m,n)= blurred image, f(m,n)= sharp image, h(m,n)= convolution kernel and w(m,n)= noise.”Convolution is a common image processing technique that changes the intensities of a pixel to reflect the intensities of the surrounding pixels”. Figure 2 shows the convolution operation in images. Figure 1- Blurring Model Figure 2- Kernel Convolution 2 1.1.3 Type of Blurring – Box blur: also called “Average blur” smoothing the image by setting each pixel equal to the average pixel value of its specified box neighborhood.- Gaussian blur: is the result of blurring an image by Gaussian function. It used to approximate out of focus blur. – Motion blur: is blurring of image due to movement of the object or camera. The motion controlled by angle and direction. 1.2 De-blurring Concept De-blurring (deconvolution) is process used to reverse the effects of convolution on image. If we have the DFT of equation (1).Simply we can find solution for f(m,n) after doing the inverse DFT ?? ( ?? 1 , ?? 2 ) = ?? ( ?? 1 , ?? 2 ) × ?? ( ?? 1 , ?? 2 ) + ?? ( ?? 1 , ?? 2 ) ……. (2) ?? ? ( ?? 1 , ?? 2 ) = ?? ( ?? 1 , ?? 2 ) ? ?? ( ?? 1 , ?? 2 ) ?? ( ?? 1 , ?? 2 ) ……… ……… …. . … … (3) The noise can be increased if the ?? ( ?? 1 , ?? 2 ) have very small value. 1.2.1Type of De-blurring Image de-blurring (deconvolution) classified into two main techniques: – Non-blind deconvolution: where the PSF and noise level are either known or estimated. – Blind deconvolution: used effectively when no information about the blurring and noise. This is the most general case of image restoration. Figure 3- Debluring Techniques 1.2.2Blind Deconvolution Blind deconvolution is deconvolution algorithm that works without specific knowledge of the impulse response function used in the convolution. The simplest form the blind deconvolution algorithm is an iterative function, that takes an initial guess for PSF and blurred image, and returns de-blurred image and new guess for PSF. Good result will obtained when the size of initial PSF used in deconvolution is equal to the size of PSF used in blurringphase. 1.2.3Wiener Filter The Weiner filter is one of most general technique used for image de-blurring. The Wiener filter reduces the mean square error (MSE) between the estimated process and the desired process. The formula for the Wiener filter: ?? ? ( ?? 1 , ?? 2 ) = ?? ( ?? 1 , ?? 2 ) ?? ( ?? 1 , ?? 2 ) ?? | ?? ( ?? 1 , ?? 2 ) | 2 | ?? ( ?? 1 , ?? 2 ) | 2 + 1 ?????? ( ?? 1 , ?? 2 ) … . ( 4 ) The original image can be recovered perfectly when there is no noise and the point spread function is known. 1.3Image Quality To evaluate the recovered image, we will compare it with sharp image, which represents a reference. Image processing toolbox in MATLAB provides several function that used to measure the quality of images: 1.3.1 Mean squared error (MSE) It represents the cumulative squared error between the deblurred and the original image.The following equation used to compute the MSE: ???? = 1 ? ? ?? ( ?? , ?? ) ? ?? ( ?? , ?? ) 2 … … . ( 5 ) ?? ? 1 ?? = 0 ?? ? 1 ?? = 0 Where M x N represents the number of pixels in image. This measure is simple to calculate but sometimes does not align well with perceived quality by humans .Image Deblurring Techniques Non Blind Deconvolution Weiner Filter Lucy-Rechardson Method Blind Deconvolution 3 1.3.2 Peak Signal to Noise Ratio (PSNR) The PSNR computes the peak signal-to-noise ratio in dB between deblurred and the original image. The higher PSNR means better quality of recovered image . The following equation used to compute the PSNR: ???? = ?????? ?? 2 ?????? …………… … ……… …….(6) Where: L represents the max value of intensity of reference image. For example, if the input image has 16-bit depth L is 65535. 1.3.3 Structure Similarity Index Map (SSIM) SSIM is used to compare luminance, contrast and structure of two different images. Where: ? ?? , ? ?? represent mean intensity, ?? ?? , ?? ?? represent standard deviation, and ?? 1 , ?? 2 represent the constant to avoid instability when ? ?? 2 + ? ?? 2 is very close to zero and is defined as ?? ?? = ( ?? ?? ?? ) 2 in which ?? ?? ? 1 and L is the dynamic range of pixel values e.g.L=65535 for 16-bit gray scale image. 2. SIMULATION AND RESULTS The preparation for applying the de-blurring algorithm can be divided into the following steps: Step1: import sharp image using imread function in MATLAB. The grayscale image used with size 512 x 512 and 16-bit depth.Figure 4- Sharp (True) Image Step2: Simulate blur using fspecial function where the where the image is degraded by Gaussian, linear motion and box (average) blur. Step3: Simulate noise by adding Gaussian noise of mean=0 and variance=1/1000 to the blurred image using imnoise function. % Read Image (True Image) %****************************** Img = imread(‘MAKKAH.tif’); % ***** Gaussian Blur ***** PSF = fspecial(‘gaussian’,7,10); BLURRED_I=imfilter(Img,PSF,’symmetric’,’co nv’); BLURRED_DESCRIPTION =’Gaussian Blurred Image’; % ***** Linear Motion Blur ***** LEN = 30; THETA = 15; PSF = fspecial(‘motion’, LEN, THETA); BLURRED_I = imfilter(Img, PSF, ‘conv’, ‘circular’); BLURRED_DESCRIPTION =’Motion Blurred Image with LEN=30 and THETA=15′; % ***** Box (Linear)(Average) ***** PSF = fspecial(‘average’, 7); BLURRED_I= imfilter(Img, PSF,’conv’, ‘circular’); BLURRED_DESCRIPTION =’Average Blurred Image’; % ***** Gaussian Noise ***** NOISE_MEAN = 0; NOISE_VAR = 0.0001; BLURRED_NOISY = imnoise(BLURRED_I, ‘gaussian’, NOISE_MEAN, NOISE_VAR); NOISE_DESCRIPTION=’Blurred and Noisy Image’; 4 2.1Wiener Filter Figure 5- Restored of Blurred Image (NSR=0) Figure 6- Restored of Blurred Noisy Image (NSR=0) Figure 7- Restored of Blurred Noisy Image (Estimated NSR) Using Wiener Filter 2.2 Blind Deconvolution Figure 8- Restored of Blurred Image Using PSFs of various Size The ringing in the de-blurred image occurs along the areas of sharp intensity contrast in the image and along the image borders. weighting function used to reduce the ringing effect.%Restored Image, No Noise – NSR = 0 DEBLURRED_W_I1 = deconvwnr(BLURRED_I, PSF, 0); %Restoration of Blurred, Noisy Image – NSR = 0 DEBLURRED_W_I2 = deconvwnr(BLURRED_NOISY, PSF, 0); %Restoration of Blurred, Noisy Image – Estimated NSR Estimated_NSR = NOISE_VAR / var(im2double(Img(:))); DEBLURRED_W_I3 = deconvwnr(BLURRED_NOISY, PSF, Estimated_NSR); %Restore the Blurred Image Using Undersized PSF DEBLURRED_B_I1, P1 = deconvblind(BLURRED_I,ones(size(PSF)-4)); %Restore the Blurred Image Using Oversized PSF DEBLURRED_B_I2, P2 = deconvblind(BLURRED_I,ones(size(PSF)+4)); %Restore the Blurred Image Using Same Size of PSF DEBLURRED_B_I3, P3 = deconvblind(BLURRED_I,ones(size(PSF))); %to improve the ringing in the restored image WEIGHT = edge(BLURRED_I,’sobel’,.1); WEIGHT = 1- double(imdilate(WEIGHT,strel(‘disk’,2))); WEIGHT(1:3 end-(0:2),:) = 0; WEIGHT(:,1:3 end-(0:2)) = 0; DEBLURRED_B_I5, P5 = deconvblind(BLURRED_NOISY,ones(size(PSF)), 25,,WEIGHT); 5 Figure 9- Restored of Blurred Noisy Using Blind Deconvolution 2.3Lucy-Richardson Algorithm Figure 10- Restored of Blurred Using Lucy-Richardson Algorithm 2.4 Comparison Results Table-1 Blurred Noiseless Image PSF PSNR SSIM Qual ity Wiener filter Gaussian 7×7 7.689 0.024 L Average 7×7 49.420 0.991 H Motion LEN = 30 THETA = 15 54.313 0.997 H blind deconvolution Gaussian 7×7 31.459 0.938 M Average 7×7 31.272 0.937 M Motion LEN = 30 THETA = 15 20.433 0.658 L Gaussian 7×7 25.972 0.863 M Average 7×7 28.759 0.895 M Motion LEN = 30 THETA = 15 24.010 0.822 M Table2-Blurred Noisy Image PSF PSNR SSIM Qual ity Wiener filter Gaussian 7×7 20.499 0.304 L Average 7×7 23.215 0.347 M Motion LEN = 30 THETA = 15 21.296 0.280 M blind deconvolution Gaussian 7×7 24.769 0.444 M Average 7×7 24.776 0.448 M Motion LEN = 30 THETA = 15 20.049 0.475 L Gaussian 7×7 25.474 0.724 M Average 7×7 27.992 0.766 M Motion LEN = 30 THETA = 15 23.727 0.694 M %Noiseless DEBLURRED_L_I1 =deconvlucy(BLURRED_I,PSF); %Noisy DEBLURRED_L_I2 = deconvlucy(BLURRED_NOISY,PSF,,im2uint16( sqrt(NOISE_VAR))); 6 3.CONCLUSION The restoration or de-blurring images is not easy problem to resolve especially when the image is noisy and there is no information about PSF. From the above analysis and after applying some of image de-blurring techniques we can see that Wiener filter is more accurate and less complex than other approaches when the image is (noiseless) where the PSNR value is 49.420 and 54.313 for average and motion blur respectively.But it has very bad result in Gaussian blurred images. In addition, noise in Wiener filter has huge impact on the quality of the recovered image. Therefore, Wiener filter works well when the noise variance is small and PSF is known. The performance of blind deconvolution is efficient in both blurred and noisy image where (PSNR between 20 and 32) but need good estimation for the PSF.In most general case in image restoration, the PSF and noise is unknown. So in the future, I will apply the algorithm to obtain a good estimation for PSF before using de-blurring algorithm. 7 4. Appendix: Matlab Code close all clear clc %——————————————————————– % ***** Read Image (True Image) ***** Img = imread(‘MAKKAH.tif’); %——————————————————————– %Blurring of Image %——————————————————————– % ***** Gaussian Blur ***** PSF = fspecial(‘gaussian’,7,10); BLURRED_I=imfilter(Img,PSF,’symmetric’,’conv’); BLURRED_DESCRIPTION =’Gaussian Blurred Image’; % % ***** Linear Motion Blur ***** % LEN = 30; % THETA = 15; % PSF = fspecial(‘motion’, LEN, THETA); % BLURRED_I = imfilter(Img,PSF, ‘conv’, ‘circular’); % BLURRED_DESCRIPTION =’Motion Blurred Image with LEN=30 and THETA=15′; % % % ***** Box (Linear)(Average) ***** % PSF = fspecial(‘average’, 7); % BLURRED_I= imfilter(Img, PSF,’conv’, ‘circular’); % BLURRED_DESCRIPTION =’Average Blurred Image’; %——————————————————————– % Add Noise to Image %——————————————————————– % ***** Gaussian Noise ***** NOISE_MEAN = 0; NOISE_VAR = 0.0001; BLURRED_NOISY = imnoise(BLURRED_I, ‘gaussian’, NOISE_MEAN, NOISE_VAR); NOISE_DESCRIPTION=’Blurred and Noisy Image’; %——————————————————————– % Show images (Original – Blurred – Noisy) %——————————————————————– figure; subplot(131); imshow(Img); title(‘Original Image’); subplot(132); imshow(BLURRED_I); title(BLURRED_DESCRIPTION); subplot(133); imshow(BLURRED_NOISY); title(NOISE_DESCRIPTION); %——————————————————————– % Wiener Filter %——————————————————————– %Restored Image, No Noise – NSR = 0 DEBLURRED_W_I1 = deconvwnr(BLURRED_I, PSF, 0); %Restoration of Blurred, Noisy Image – NSR = 0 DEBLURRED_W_I2 = deconvwnr(BLURRED_NOISY, PSF, 0); 8 %Restoration of Blurred, Noisy Image – Estimated NSR Estimated_NSR = NOISE_VAR / var(im2double(Img(:))); DEBLURRED_W_I3 = deconvwnr(BLURRED_NOISY, PSF, Estimated_NSR); % Show images for wiener filter figure; subplot(221); imshow(DEBLURRED_W_I1); title(‘Restored Image, No Noise – NSR = 0 – Wiener -‘); subplot(222); imshow(DEBLURRED_W_I2) title(‘Restoration of Blurred, Noisy Image – NSR = 0 – Wiener -‘) subplot(223); imshow(DEBLURRED_W_I3) title(‘Restoration of Blurred, Noisy Image – Estimated NSR – Wiener -‘); % Quality of De-blurring fprintf(‘
Quality of De-blurring’); fprintf(‘
***** Wiener Filter (Noisless) *****’); fprintf(‘
MSE = %0.4f’, immse(DEBLURRED_W_I1, Img)); %MSE fprintf(‘
PSNR = %0.4f’, psnr(DEBLURRED_W_I1, Img)); %PSNR fprintf(‘
SSIM = %0.4f
‘,ssim(DEBLURRED_W_I1,Img)); %SSIM fprintf(‘
***** Wiener Filter (Noisy) *****’); fprintf(‘
MSE = %0.4f’, immse(DEBLURRED_W_I3, Img)); %MSE fprintf(‘
PSNR = %0.4f’, psnr(DEBLURRED_W_I3, Img)); %PSNR fprintf(‘
SSIM = %0.4f
‘,ssim(DEBLURRED_W_I3,Img)); %SSIM %——————————————————————– %Blind Deconvolution %——————————————————————– %Restore the Blurred Image Using Undersized PSF DEBLURRED_B_I1, P1 = deconvblind(BLURRED_I,ones(size(PSF)-4)); %Restore the Blurred Image Using Oversized PSF DEBLURRED_B_I2, P2 = deconvblind(BLURRED_I,ones(size(PSF)+4)); %Restore the Blurred Image UsingSame Size of PSF DEBLURRED_B_I3, P3 = deconvblind(BLURRED_I,ones(size(PSF))); %to improve the ringing in the restored image WEIGHT = edge(BLURRED_I,’sobel’,.1); WEIGHT = 1-double(imdilate(WEIGHT,strel(‘disk’,2))); WEIGHT(1:3 end-(0:2),:) = 0; WEIGHT(:,1:3 end-(0:2)) = 0; %(Niseless) DEBLURRED_B_I4, P4 = deconvblind(BLURRED_I,ones(size(PSF)),25,,WEIGHT); %(Noisy) DEBLURRED_B_I5, P5 = deconvblind(BLURRED_NOISY,ones(size(PSF)),25,,WEIGHT); % Show images for Blind Deconvolution figure; subplot(321); imshow(DEBLURRED_B_I1); title(‘Restoration of Blurred Image Using Undersized PSF – Blind Deconv.-‘); subplot(322); imshow(DEBLURRED_B_I2) title(‘Restoration of Blurred Image Using Oversized PSF – Blind Deconv. -‘) subplot(323); imshow(DEBLURRED_B_I3) title(‘Restoration of Blurred Image Using Same Size of PSF – Blind Deconv. -‘); 9 subplot(324); imshow(DEBLURRED_B_I4) title(‘Restoration of Blurred Image – Blind Deconv.-‘); subplot(325); imshow(DEBLURRED_B_I5) title(‘Restoration of Blurred, Noisy Image – Blind Deconv. -‘); % Show the Reconstructed PSF figure; subplot(321);imshow(PSF,,’InitialMagnification’,’fit’); title(‘True PSF’); subplot(322);imshow(P1,,’InitialMagnification’,’fit’); title(‘Reconstructed Undersized PSF’); subplot(323);imshow(P2,,’InitialMagnification’,’fit’); title(‘Reconstructed Oversized PSF’); subplot(324);imshow(P3,,’InitialMagnification’,’fit’); title(‘Reconstructed same size of true PSF’); subplot(325);imshow(P4,,’InitialMagnification’,’fit’); title(‘Reconstructed true PSF’); % Quality of De-blurring fprintf(‘
Quality of De-blurring’); fprintf(‘
***** Blind Deconvolution (Noisless) *****’); fprintf(‘
MSE = %0.4f’, immse(DEBLURRED_B_I4, Img)); %MSE fprintf(‘
PSNR = %0.4f’, psnr(DEBLURRED_B_I4, Img)); %PSNR fprintf(‘
SSIM = %0.4f
‘,ssim(DEBLURRED_B_I4,Img)); %SSIM fprintf(‘
***** Blind Deconvolution (Noisy) *****’); fprintf(‘
MSE = %0.4f’, immse(DEBLURRED_B_I5, Img)); %MSE fprintf(‘
PSNR = %0.4f’, psnr(DEBLURRED_B_I5, Img)); %PSNR fprintf(‘
SSIM = %0.4f
‘,ssim(DEBLURRED_B_I5,Img)); %SSIM %——————————————————————– %Lucy-Richardson %——————————————————————– %Noiseless DEBLURRED_L_I1 = deconvlucy(BLURRED_I,PSF); %Noisy DEBLURRED_L_I2 = deconvlucy(BLURRED_NOISY,PSF,,im2uint16(sqrt(NOISE_VAR))); % Show images for Lucy-Richardson figure; subplot(121); imshow(DEBLURRED_L_I1); title(‘Restoration of Blurred Image – Lucy-Richardson -‘); subplot(122);imshow(DEBLURRED_L_I2) title(‘Restoration of Blurred Noisy Image – Lucy-Richardson -‘) % Quality of De-blurring fprintf(‘
Quality of De-blurring’); fprintf(‘
***** Lucy-Richardson (Noisless) *****’); fprintf(‘
MSE = %0.4f’, immse(DEBLURRED_L_I1, Img)); %MSE fprintf(‘
PSNR = %0.4f’, psnr(DEBLURRED_L_I1, Img)); %PSNR fprintf(‘
SSIM = %0.4f
‘,ssim(DEBLURRED_L_I1,Img)); %SSIM % fprintf(‘
***** Lucy-Richardson (Noisy) *****’); fprintf(‘
MSE = %0.4f’, immse(DEBLURRED_L_I2, Img)); %MSE fprintf(‘
PSNR = %0.4f’, psnr(DEBLURRED_L_I2, Img)); %PSNR fprintf(‘
SSIM = %0.4f
‘,ssim(DEBLURRED_L_I2,Img)); %SSIM % ——————————————————————–INTERNET SOURCES: ——————————————————————————————