clear;
aerial = imread('Fig3.09(a).jpg');
aerial = double(aerial);
[nx ny] = size(aerial);
nx
ny
u = aerial;

for i = 1:nx
  for j = 1:ny
      % This is aerial with periodic noise.
      u(i,j) = u(i,j) + ...
                    5.*(1+sin(2*pi*((i-1)/nx)*200))+...
                    5.*(1+sin(2*pi*((j-1)/ny)*200))+...
                    5.*(1+cos(2*pi*((i-1)/nx+(j-1)/ny)*141))+...
                    5.*(1+sin(2*pi*((i-1)/nx-(j-1)/ny)*141));;
  end
end
sineaerial = uint8(u);
imwrite(sineaerial, 'sineaerial.jpg');

c = 1.;      % Use the power transformation to darken.
gamma = 2;
f_fp =255*c*(u/255).^gamma;
u = f_fp;
fftu = fft2(u,2*nx-1,2*ny-1);
fftu = fftshift(fftu);
subplot(1,2,1)
mesh(log(1+(abs(fftu))));

filter = ones(2*nx-1,2*ny-1);
d0 = 400;      % Use Butterworth band reject filter.
n = 4;
w = 20;
for i = 1:2*nx-1
    for j =1:2*ny-1
        dist = ((i-(nx+1))^2 + (j-(ny+1))^2)^.5;
        if dist ~= d0
            filter(i,j)= 1/(1 + (dist*w/(dist^2 - d0^2))^(2*n));
        else
            filter(i,j) = 0;
        end
    end
end 
fil_aerial  = filter.*fftu;
subplot(1,2,2)
mesh(log(1+abs(fil_aerial)));

fil_aerial = ifftshift(fil_aerial);
fil_aerial = ifft2(fil_aerial,2*nx-1,2*ny-1);
fil_aerial = real(fil_aerial(1:nx,1:ny));
fil_aerial = uint8(fil_aerial);
imwrite(fil_aerial, 'sineaerial_fil.jpg');