Portofolio 3

MyPhoto.jpg

im = imread('myPhoto.jpg')
recon = freqcomp(im, 50, 0)


recon = freqcomp(im, 150, 0)


Gaussian Filter:







gaussfilter = gaussian('myPhoto.jpg')


Average Filter:








avfilt = Average('myPhoto.jpg')


pertukaran nilai fase:



im1



im2





portofolio3 = no4('me1.jpg','me2.jpg')






portofolio3 = no5('me1.jpg','me2.jpg')

Listing of freecomp.m:

% FREQCOMP - Demonstrates image reconstruction from Fourier components
%
% Usage: recon = freqcomp(im, Npts, delay)
%
% Arguments: im - Image to be reconstructed.
% Npts - Number of frequency components to consider
% (defaults to 50).
% delay - Optional time delay between animations of the
% reconstruction. If this is omitted the function
% waits for a key to be pressed before moving to
% the next component.
%
% Returns: recon - The image reconstructed from the specified
% number of components
%
% This program displays:
%
% * The image.
% * The Fourier transform (spectrum) of the image with a conjugate
% pair of Fourier components marked with red dots.
% * The sine wave basis function that corresponds to the Fourier
% transform pair marked in the image above.
% * The reconstruction of the image generated from the sum of the sine
% wave basis functions considered so far.

% Copyright (c) 2002-2003 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.

% March 2002 - Original version
% April 2003 - General cleanup of code

function recon = freqcomp(im, Npts, delay)
if ndims(im) == 3
im = rgb2gray(im);
warning('converting colour image to greyscale');
end

if nargin < 2
Npts = 50;
end

[rows,cols] = size(im);

% If necessary crop one row and/or column so that there are an even
% number of rows and columns, this makes things easier later.
if mod(rows,2) % odd
rows = rows-1;
end
if mod(cols,2) % odd
cols = cols-1;
end

im = im(1:rows,1:cols);
rc = fix(rows/2)+1; % Centre point
cc = fix(cols/2)+1;

% The following code constructs two arrays of coordinates within the FFT
% of the image that correspond to complex conjugate pairs of points that
% spiral out from the centre visiting every frequency component on
% the way.
p1 = zeros(Npts,2); % Path 1
p2 = zeros(Npts,2); % Path 2

m = zeros(rows,cols); % Matrix for marking visited points
m(rc,cc) = 1;
m(rc,cc-1) = 1;
m(rc,cc+1) = 1;

p1(1,:) = [rc cc-1];
p2(1,:) = [rc cc+1];
d1 = [0 -1]; % initial directions of the paths
d2 = [0 1];

% Mark out two symmetric spiral paths out from the centre (I wish I
% could think of a neater way of doing this)

for n = 2:Npts
l1 = [-d1(2) d1(1)]; % left direction
l2 = [-d2(2) d2(1)];

lp1 = p1(n-1,:) + l1; % coords of point in left direction
lp2 = p2(n-1,:) + l2;

if ~m(lp1(1), lp1(2)) % go left
p1(n,:) = lp1;
d1 = l1;
m(p1(n,1), p1(n,2)) = 1; % mark point as visited
else % go sraight ahead
p1(n,:) = p1(n-1,:) + d1;
m(p1(n,1), p1(n,2)) = 1; % mark point as visited
end

if ~m(lp2(1), lp2(2)) % go left
p2(n,:) = lp2;
d2 = l2;
m(p2(n,1), p2(n,2)) = 1; % mark point as visited
else % go sraight ahead
p2(n,:) = p2(n-1,:) + d2;
m(p2(n,1), p2(n,2)) = 1; % mark point as visited
end

end

% Having constructed the path of frequency components to be visited
% we take the FFT of the image and then enter a loop that
% incrementally reconstructs the image from its components.

IM = fftshift(fft2(im));
recon = zeros(rows,cols); % Initialise reconstruction matrix

if max(rows,cols) < 150
fontsze = 7;
else
fontsze = 10;
end

figure(1), clf
subplot(2,2,1),imagesc(im),colormap gray, axis image, axis off
title('Original Image','FontSize',fontsze);
subplot(2,2,2),imagesc(log(abs(IM))),colormap gray, axis image
axis off, title('Fourier Transform + frequency component pair','FontSize',fontsze);

warning off % Turn off warnings that might arise if the images cannot be
% displayed full size
truesize(1)

for n = 1:Npts

% Extract the pair of Fourier components
F = zeros(rows,cols);
F(p1(n,1), p1(n,2)) = IM(p1(n,1), p1(n,2));
F(p2(n,1), p2(n,2)) = IM(p2(n,1), p2(n,2));

% Invert and add to reconstruction
f = real(ifft2(fftshift(F)));
recon = recon+f;

% Display results
subplot(2,2,2),imagesc(log(abs(IM))),colormap gray, axis image
axis off, title('Fourier Transform + frequency component pair','FontSize',fontsze);
hold on, plot([p1(n,2), p2(n,2)], [p1(n,1), p2(n,1)],'r.'); hold off
subplot(2,2,3),imagesc(recon),colormap gray, axis image, axis off, title('Reconstruction','FontSize',fontsze);
subplot(2,2,4),imagesc(f),colormap gray, axis image, axis off
title('Basis function corresponding to frequency component pair','FontSize',fontsze);

if nargin == 3
pause(delay);
else
fprintf('Hit any key to continue \n'); pause
end

end

warning on % Restore warnings


listing of gaussian.m:

function gaussfilter = gaussian(img)
%gaussian filter
im = imread(img);
im = imresize(im,0.5)
im = rgb2gray(im);

gaussianfilter = fspecial('gaussian', size(im), 6);
imfft = fft2(im);
gaussianfilterfft = fft2(gaussianfilter);

newimagegafft = gaussianfilterfft.*imfft;

figure,imshow(im);
figure,imagesc(gaussianfilter),title('Gaussian filter'),colormap gray;
figure,imagesc(fftshift(log(abs(imfft)+eps))),title('FFT'),colormap gray;
figure,imagesc(fftshift(log(abs(gaussianfilterfft)+eps))),title('Gaussian filter FFT'),colormap gray;
figure,imagesc(fftshift(log(abs(newimagegafft)+eps))),title('Gaussian filter kali FFT'),colormap gray;
figure,imagesc(fftshift(real(ifft2(newimagegafft)))),title('smoothing'),colormap gray;


listing of average.m:

function avfilt = Average(img)
%Average filter
im = imread(img);
im = imresize(im,0.5);
im = rgb2gray(im);
averagefilt = averagefilter(im,[21 21]);

imfft = fft2(im);
averagefiltfft = fft2(averagefilt);
imaveragefft = averagefiltfft.*imfft;

figure,imshow(im);
figure,imagesc(averagefilt),title('Average filter'),colormap gray;
figure,imagesc(fftshift(log(abs(imfft)+eps))),title('FFT'),colormap gray;
figure,imagesc(fftshift(log(abs(averagefiltfft)+eps))),title('Average filter FFT'),colormap gray;
figure,imagesc(fftshift(log(abs(imaveragefft)+eps))),title('Average filter kali FFT'),colormap gray;
figure,imagesc(fftshift(real(ifft2(imaveragefft)))),title('smoothing'),colormap gray;


listing of averagefilter.m:

function av = averagefilter(img,n)

avefilt = fspecial('average',n);
[rows cols]=size (img);
av=zeros(rows,cols);
av((floor((rows+1)/2)-floor((n(1,1)+1)/2)):(floor((rows+1)/2)+floor((n(1,1)+1)/2)),(floor((cols+1)/2)-floor((n(1,2)+1)/2)):(floor((cols+1)/2)+floor((n(1,2)+1)/2)))=avefilt(1,1);


listing of no4.m:

function portofolio3 = no4(img1,img2)
im1 = imread(img1);
im2 = imread(img2);
% greyscale
im1 = rgb2gray(im1);
im2 = rgb2gray(im2);

im1fft = fft2(im1);
im2fft = fft2(im2);

im1phase = angle(im1fft);
im2phase = angle(im2fft);

im1mag = abs(im1fft);
im2mag = abs(im2fft);

imagefft1 = im1mag.*(cos(im2phase) + i*sin(im2phase));
imagefft2 = im2mag.*(cos(im1phase) + i*sin(im1phase));

figure,imagesc(real(ifft2(imagefft1))),title('image1'),colormap gray;
figure,imagesc(real(ifft2(imagefft2))),title('image2'),colormap gray;


listing of no5.m:

function potofolio3 = no5(img1,img2)
im1 = imread(img1);
im2 = imread(img2);
% greyscale
im1 = rgb2gray(im1);
im2 = rgb2gray(im2);

im1fft = fft2(im1);
im2fft = fft2(im2);

im1mag = abs(im1fft);
im2mag = abs(im2fft);

image1phase = im1fft./im1mag;
image2phase = im2fft./im2mag;

figure,imagesc(real(ifft2(image1phase))),title('image1'),colormap gray;
figure,imagesc(real(ifft2(image2phase))),title('image2'),colormap gray;

figure,imagesc(im1),title('image1 asli'),colormap gray;
figure,imagesc(im2),title('image2 asli'),colormap gray;