function [SNR_ResIm , ResIm] = wiener_laplacian(OrigImage, Noised_Blurred_Image, Blur, Noise, type,motion)
[x,y] = size(Noised_Blurred_Image);    %get size of the image
K = 0.5;
gamma = 25;

G_fft=fft2(Noised_Blurred_Image);
Blur_fft = fft2(Blur,x,y);
Blur_fft_conj = conj(Blur_fft);
Blur_fft_sq = Blur_fft .* Blur_fft_conj;

OrigImagefft = fft2(OrigImage);
Noisefft = fft2(Noise,x,y);

Laplacian = [0 1 0; 1 -4 1; 0 1 0]; %Laplacian
Laplacian2 = [1 1 1; 1 -8 1; 1 1 1]; %Laplacian

Laplacian2_fft = fft2(Laplacian2, x , y);
Laplacian2_fft_conj = conj(Laplacian2_fft);
Laplacian2_fft_sq = Laplacian2_fft .* Laplacian2_fft_conj;

Laplacian_fft = fft2(Laplacian, x , y);
Laplacian_fft_conj = conj(Laplacian_fft);
Laplacian_fft_sq = Laplacian_fft .* Laplacian_fft_conj;

switch type
    
    case 1 % Constrained least squares restoration method
        DeBlurred_fft = (Blur_fft_conj .* G_fft) ./ (Blur_fft_sq + (gamma .* Laplacian_fft_sq));
        
    case 2 % Original Wiener
        Noise_sigma = fspecial('gaussian',[x,y],Noise);
        Noise_sigma_fft = fft2(Noise_sigma);
        
        Image_sigma_fft = fft2(OrigImage);
        N_sq = real(Noise_sigma_fft).^2 + imag(Noise_sigma_fft).^2;
        F_sq = real(Image_sigma_fft).^2 + imag(Image_sigma_fft).^2;
        K = N_sq ./ F_sq;
        DeBlurred_fft = (Blur_fft_conj .* G_fft) ./ (Blur_fft_sq + K);
        
  
    case 3 % Original Wiener with contsant
        
        DeBlurred_fft = (Blur_fft_conj .* G_fft) ./ (Blur_fft_sq + K);
        
%     case 4 % Recursive Wiener % Implemented in another file
%         Noise_Im = (im2double(Noised_Blurred_Image) - im2double(PrevIm));
%         Noise_Im_fft = fft2(Noise_Im);
%         ResIm_fft = fft2(PrevIm);
%         ratio = (Noise_Im_fft ./ ResIm_fft);
%         DeBlurred_fft = (Blur_fft_conj .* G_fft) ./ (Blur_fft_sq + ratio);
%         
    case 5 % Power Spectrum Equalization method
        %Noise_Im = (im2double(Noised_Blurred_Image) - im2double(OrigImage));
        Noise_sigma = fspecial('gaussian',[x,y],Noise);
        Noise_im_fft = fft2(Noise_sigma);
        
        Image_im_fft = fft2(OrigImage);
        
        L = sqrt(((Image_im_fft) ./ (Blur_fft_sq .* Image_im_fft + Noise_im_fft)));
        DeBlurred_fft = G_fft .* L;
end     

DeBlurred = ifft2(DeBlurred_fft);
DeBlurred_normolized = DeBlurred ./ max(max(DeBlurred));
ResIm = DeBlurred_normolized;

        
%figure;
A = deconvwnr(Noised_Blurred_Image,Blur);
%imshow(A);title('Matlab filtered Image: deconvwnr');

B = wiener2(Noised_Blurred_Image);
figure;
%imshow(B);title('Matlab filtered Image: wiener2');


subplot(2,2,1);imshow(Noised_Blurred_Image);title('Blurred & Noises image');
tmp = sprintf('Motion = %d pixel , Noise = %0.3g', motion , Noise);
subplot(2,2,2);imshow(ResIm);title(tmp);
subplot(2,2,3);imshow(uint8(A));title('Matlab filtered Image: deconvwnr');
subplot(2,2,4);imshow(uint8(B));title('Matlab filtered Image: wiener2');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
OrigImage = double(OrigImage);
OrigImage_Normolized = OrigImage ./ max(max(OrigImage));
%%%%%%%%%%%%%% RES %%%%%%%%%%%%%%%%%%%%%%%
diff1 = abs(im2double(DeBlurred_normolized) - im2double(OrigImage_Normolized)); %Noised_Blurred_Image
diff_sq1 = diff1' * diff1;
a = sum(sum(diff_sq1)) / (x*y)
im_mean1 = mean(mean(diff_sq1));
SNR_ResIm = 10*log10(max(max(im2double(OrigImage_Normolized)))/im_mean1) %Noised_Blurred_Image

%%%%%%%%%%%%%% A %%%%%%%%%%%%%%%%%%%%%%%
diff2 = abs(im2double(A) - im2double(OrigImage_Normolized));
diff_sq2 = diff2' * diff2;
a = sum(sum(diff_sq2)) / (x*y)
im_mean2 = mean(mean(diff_sq2));
SNR_A = 10*log10(max(max(im2double(OrigImage_Normolized)))/im_mean2)
%%%%%%%%%%%%% B %%%%%%%%%%%%%%%%%%%%%%%

diff3 = abs(im2double(B) - im2double(OrigImage_Normolized));
diff_sq3 = diff3' * diff3;
a = sum(sum(diff_sq3)) / (x*y)
im_mean3 = mean(mean(diff_sq3));
SNR_B = 10*log10(max(max(im2double(OrigImage_Normolized)))/im_mean3)