function [blurred_image,h] = image_blurring (image, filter)

image = double(image); % option= 'DCT';

[blurred_image,h] = Image_distortion_periodic(image, filter);  

return;


%% embedded function 1
function [y_blurred,h] = Image_distortion_periodic(image,filter)

x = double(image); [M,N] = size(x);

switch filter{1}
    
    case 'gaussian'
        std = filter{2};
%     h = fspecial('gaussian', 25, 1.6);
    h = gaussian_kernel (255,255, std);
    centered_h = centered_kernel_2D(h,M,N); h = centered_h;
    H = fft2(centered_h);
    y_blurred = Image_degradation_periodic(H,x);

    case 'jinc'
        s = filter{2};   size_h = 55;      center = (size_h+1)/2;
        h = zeros(size_h, size_h); xx = h;
    for i=1:size_h
        for j=1:size_h
            xx(i,j) = sqrt( (i-center)^2+(j-center)^2 );   
        end
    end
    xx(center, center) = eps; xx =xx/s;
    h = besselj(1,xx)./xx; h = h.*h;     h = h./sum(sum(h));
    h = centered_kernel_2D(h,M,N);    H = fft2(h);
    y_blurred = Image_degradation_periodic(H,x);
    
    case 'non_symmetric'
    s = filter{2};   theta = filter{3};   sigma_1 = filter{4};    sigma_2 = filter{5}; 
    h = non_symmetric_gaussian_kernel ( 255, 255, theta,   sigma_1, sigma_2, s ) ;
    h = centered_kernel_2D(h,M,N);    H = fft2(h);
    y_blurred = Image_degradation_periodic(H,x);
    
end

return;

%% embedded function 2
function y_blurred = Image_degradation_periodic(H,x)

X = fft2(x); Y1 = H.*X; y_blurred = ifft2(Y1);

return;


%% embedded function 5
function h = gaussian_kernel(m,n,sigma)
h = ones(m,n);
center_x = (m+1)/2; center_y = (n+1)/2;
for i = 1:m
    for j = 1:n
        h(i,j) = exp(-((i-center_x)^2+(j-center_y)^2)/(2*sigma^2));
    end
end
h = h./sum(h(:));

return;


%% embedded function 5
function h = non_symmetric_gaussian_kernel ( m, n,  theta,    sigma_1, sigma_2, s )
h = ones(m,n);  theta = -theta;
center_x = (m+1)/2; center_y = (n+1)/2;
for i = 1:m
    for j = 1:n
        i_x = i-center_x; j_y = j-center_y; 
        h(i,j) = exp(- (i_x*cos(theta)-j_y*sin(theta) )^2/(1*s*sigma_1^2) ...
            - ( (i_x*sin(theta)+j_y*cos(theta) )^2/(1*s*sigma_2^2) ) );
    end
end
h = h./sum(h(:));

return;

%% embedded function 6
function centered_h = centered_kernel_2D(h,M,N)

[m,n] = size(h); size_h = m; centered_h = zeros(M,N);
clear m n;

if mod(size_h,2) == 1
% if size_h is odd ---> whole-point symmetry, often referring to convolution kernel 
% top-left
centered_h(1:(size_h+1)/2,1:(size_h+1)/2) = h((size_h+1)/2:size_h,(size_h+1)/2:size_h);
% bottom-right
centered_h(M-(size_h-1)/2+1:M,N-(size_h-1)/2+1:N) = h(1:(size_h-1)/2,1:(size_h-1)/2);
% top-right
centered_h(1:(size_h+1)/2,N-(size_h-1)/2+1:N) = h((size_h+1)/2:size_h,1:(size_h-1)/2);
% bottom-left
centered_h(M-(size_h-1)/2+1:M,1:(size_h+1)/2) = h(1:(size_h-1)/2,(size_h+1)/2:size_h);
end


if mod(size_h,2) == 0
% if size_h is even ---> half-point symmetry, often referring
% to Haar wavelet, for example.
% top-left of h
centered_h(1:(size_h/2+1),1:(size_h/2+1)) = h(size_h/2:end,size_h/2:end);
% bottom-right of h
centered_h(end-(size_h/2-1)+1:end,end-(size_h/2-1)+1:end) = h(1:size_h/2-1,1:size_h/2-1);
% top-right of h
centered_h(1:size_h/2+1,end-(size_h/2-1)+1:end) = h(size_h/2:end,1:size_h/2-1);
% bottom-left of h
centered_h(end-(size_h/2-1)+1:end,1:size_h/2+1) = h(1:size_h/2-1,end-(size_h/2+1)+1:end);
end

return;

