频域低通：理想、巴特沃斯、高斯

im = imread('circlesBrightDark.png');
im = double(im)/255;

figure;
subplot(121);
imshow(im);
title('Original Image');
        
IM = fftshift(fft2(im));
subplot(122);
imshow(mat2gray(log(abs(IM) + 1e-10)));
title('Frequency Domain ')

imsize = size(im);
D = 1:imsize(2);
V = 1:imsize(1);
[D, V] = meshgrid(D, V);
R = sqrt((D-imsize(2)/2).^2 + (V-imsize(1)/2).^2);

H1 = zeros(imsize);
r = imsize(1)/20;
H1(R <= r) = 1;

myfilter(im, H1, 'Ideal Low-pass Filter')

N = 2;
H2 = 1./(1 + (R/r).^(2*N));
myfilter(im, H2, 'Butterworth Low-pass Filter');

sigma = 20;
H3 = exp(-R.^2./sigma.^2/2);
myfilter(im, H3, 'Gaussian Low-pass Filter');

频域高通：理想、巴特沃斯、高斯、同态滤波

H4 = 1 - H1;
myfilter(im, H4, 'Ideal High-pass Filter');

N = 2;
H5 = 1./(1 + (r./R).^(2*N));
myfilter(im, H5, 'Butterworth High-pass Filter');

H6 = 1 - H3;
myfilter(im, H6, 'Gaussian High-pass Filter');

im = imread('mandi.tif');
im = double(im)/255;
figure;
subplot(221);
imshow(im);
title('Original Image');
        
IM = fftshift(fft2(im));
subplot(222);
imshow(mat2gray(log(abs(IM) + 1e-10)));
title('Frequency Domain ')

imsize = size(im);
D = 1:imsize(2);
V = 1:imsize(1);
[D, V] = meshgrid(D, V);
R = sqrt((D-imsize(2)/2).^2 + (V-imsize(1)/2).^2);

r_H = 1;
r_L = 0.4;
c = 0.2;
r = imsize(1) / 20;

H7 = (r_H - r_L)*(1 - exp(-c*((R/r).^2))) + r_L;

im_In = log(im + 1);
IM_In = fftshift(fft2(im_In));
subplot(223);
imshow(H7);
title('Filter')

IM_In_ = IM_In.*H7;

im_In_ = ifft2(fftshift(IM_In_));
im_ = exp(abs(im_In_)) - 1;
im_ = mat2gray(im_);

subplot(224);
imshow(im_);
title('Result');

suptitle('Homomorphic Filter')

function myfilter(im, H, imtitle)
    
    IM = fftshift(fft2(im));
    
    figure;
    subplot(121);
    imshow(H);
    title('Filter');
    
    IM_ = IM.*H;
    im_ = ifft2(fftshift(IM_));
    subplot(122);
    imshow(abs(im_));
    title('Result');
    
    suptitle(imtitle);
end