yudha weblog

Portofolio 7

2008-05-18T03:25:00.001-07:00

Portofolio 7:

>>portofolio 7

hasil:

pt3D =

-286.0283 159.4041 142.0436
-177.4378 163.8930 143.5515
-177.4575 17.6164 140.7850
-290.2918 12.1955 145.2484
-286.7534 160.7354 8.3899
-287.7561 12.6779 6.4816
-174.6022 14.6008 2.7149
-84.5654 -68.1487 0.4912
-68.3572 -146.0519 123.5506
18.9348 -178.8454 -4.1677
-136.9619 -207.2383 -1.7911
136.3903 -91.5529 62.7009
206.1682 -94.1534 62.1509
208.5392 -158.4030 57.5592
137.3153 -155.6600 59.8104
133.9942 -92.6254 -5.9321
133.3809 -160.9251 -8.6149
204.9606 -163.8064 -10.6342

Length matrix of face sides

slengths =

147.3052 108.6937 146.3028 113.0525
113.0525 138.1325 113.2328 138.7908
138.7908 148.0732 133.6623 147.3052
111.9686 133.6623 108.6937 138.7908
138.7908 146.3028 138.1325 149.7892
149.7892 113.2328 148.0732 111.9686

nface =

4 1 2 3 4
4 3 7 6 4
4 6 5 1 4
8 5 1 2 8
8 2 3 7 8
8 7 6 5 8

Length matrix of face sides

slengths =

146.5443 158.1369 151.6171
146.5443 155.4379 148.6489
151.6171 158.4790 148.6489
158.1369 158.4790 155.4379

nface =

1 2 3 1
1 2 4 1
1 3 4 1
2 3 4 2

Length matrix of face sides

slengths =

64.1789 69.8286 64.4571 71.3123
71.3123 68.5008 71.6662 68.7402
68.7402 68.3551 68.6831 64.1789
68.5777 68.6831 69.8286 68.7402
68.7402 64.4571 68.5008 64.5929
64.5929 71.6662 68.3551 68.5777

nface =

4 1 2 3 4
4 3 7 6 4
4 6 5 1 4
8 5 1 2 8
8 2 3 7 8
8 7 6 5 8

listing of portofolio7:

function portofolio7()
im1 = imread('stereo1.jpg' );
im2 = imread('stereo2.jpg' );

C1 = [0.6596 -0.7391 -0.0615 363.4235;
-0.1851 -0.1387 -0.9437 342.7417;
0.0005 0.0003 -0.0003 1.0000];
C2 = [0.9234 -0.2221 -0.0257 347.7796;
-0.0741 -0.2278 -0.9168 339.8960;
0.0002 0.0004 -0.0002 1.0000];

pt3D = stereo(im1, im2, C1, C2)

figure(3)
cube(pt3D(1:7,:));
tetrahedron(pt3D(8:11,:));
cube(pt3D(12:18,:));

% label coordinate axes
text(100,0,0,'x');
text(0,100,0,'y');
text(0,0,100,'z');

% draw in a set of coordinate axes
axislength = 100*eye(3);
for i=1:3
line([0, axislength(i,1)], [0, axislength(i,2)], [0, axislength(i,3)]);
end
axis equal; box on; rotate3D on; grid on;
end

function cube(cubepts3D)
% determine hidden vertex
cubepts3D(8,:) = - cubepts3D(4,:) + cubepts3D(6,:) + cubepts3D(2,:);
% define faces from standard numbering
cubefaces = [4 1 2 3
4 3 7 6
4 6 5 1
8 5 1 2
8 2 3 7
8 7 6 5];
% draw 'patcheds' from vertice and face matrix
patch('Faces',cubefaces,'Vertices',cubepts3D, 'FaceColor', 'none')
fprintf(1, 'Length matrix of face sides\n');
[slengths, nface] = sidelengths(cubepts3D, cubefaces)
end

function tetrahedron(tetrahedronpts3D)
% define faces from standard numbering
tetrahedronfaces = [1 2 3
1 2 4
1 3 4
2 3 4];
% draw 'patcheds' from vertice and face matrix
patch('Faces',tetrahedronfaces,'Vertices',tetrahedronpts3D, 'FaceColor', 'none')
fprintf(1, 'Length matrix of face sides\n');
[slengths, nface] = sidelengths(tetrahedronpts3D, tetrahedronfaces)
end

function [slengths, nface] = sidelengths(pt3D, face)
[rows, cols] = size(face);
nface = [face face(:,1)];
for i=1:cols
for j=1:rows
slengths(j,i) = norm(pt3D(nface(j,i),:)-pt3D(nface(j,i+1),:));
end
end
end

listing of digipits.m:

% DIGIPTS - digitise points in an image
%
% Function to digitise points in an image. Points are digitised by clicking
% with the left mouse button. Clicking any other button terminates the
% function. Each location digitised is marked with a red '+'.
%
% Usage: [u,v] = digipts
%
% where u and v are nx1 arrays of x and y coordinate values digitised in
% the image.
%
% This function uses the cross-hair cursor provided by GINPUT. This is
% much more useable than IMPIXEL

% Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% pk @ csse uwa edu au
% http://www.csse.uwa.edu.au/~pk
%
% May 2002

function [u,v] = digipts

hold on
u = []; v = [];
but = 1;
while but == 1
[x y but] = ginput(1);
if but == 1
u = [u;x];
v = [v;y];

plot(u,v,'r+');
end
end

hold off

listing of stereo.m:

% STEREO
%
% Usage: pt3D = stereo(im1, im2, C1, C2)
%
% Where: im1 and im2 are the two stereo images
% C1 and C2 are the calibration matrices
% for these two images respectively
%
% The function will prompt you to digitise some points in the first image
% (finishing by clicking the right button). The function then
% prompts you to digitise the equivalent points (which you must digitise
% in exactly the same sequence) in the second image.
% The function then solves the stereo equations
% and returns the 3D coordinates of the points in pt3D.
function [XYZ, uv1, uv2] = stereo(im1, im2, C1, C2)

fprintf(1, 'Digitise some points in figure 1\n');
figure(1)
imshow(im1);
[u1,v1] = digipts;
uv1 = [u1,v1]';
fprintf(1, 'Digitise some points in figure 2\n');
figure(2)
imshow(im2);
[u2,v2] = digipts;
uv2 = [u2,v2]';
% check if same number of points are selected
if length(u1) ~= length(u2)
fprintf(1, 'Same number of points not selected\n');
end

for i = 1:length(u1)
a = [C1(1:2,1:3) - [u1(i)*C1(3,1:3); v1(i)*C1(3,1:3)];
C2(1:2,1:3) - [u2(i)*C2(3,1:3); v2(i)*C2(3,1:3)]];
c = [u1(i) - C1(1,4);
v1(i) - C1(2,4);
u2(i) - C2(1,4);
v2(i) - C2(2,4)];
b(:, i) = a \ c;
end
XYZ = b';

Portofolio 6

2008-05-18T03:05:00.000-07:00

Portofolio 6:

stereo1.jpg:

im1 = imread('stereo1.jpg');
imshow(im1),uv = ginput(12);

uv = uv';
calibpts;
xyz = calibPts;
c = calibrate(im1,xyz,uv);

hasil:

C =

0.6789 -0.7263 -0.0627 362.5067
-0.1722 -0.1175 -0.9483 342.2501
0.0006 0.0003 -0.0003 1.0000

mean squared error is 2.986016e-001
error in satisfying the camera calibration matrix is 5.820335e-009

stereo2.jpg:

im2= imread('stereo1.jpg');
imshow(im2),uv = ginput(12);

uv = uv';
calibpts;
xyz = calibPts;
c = calibrate(1m2,xyz,uv);

hasil:

C =

0.9448 -0.1889 -0.0717 348.9422
-0.0682 -0.2070 -0.9236 338.6192
0.0002 0.0005 -0.0003 1.0000

mean squared error is 1.379448e-001
error in satisfying the camera calibration matrix is 6.008675e-010

listing of calibpts.m:

% 3D Coordinate data of the target calibration points
% All positions are expressed in millimetres
%
% If you run this script by simply typing
%
% >> calibpts
%
% you will end up with a local 12 x 3 matrix variable called
% calibPts containing these data points.

calibPts = [ 49 0 65 % x y z coords of point 1
129 0 65 % x y z coords of point 2
49 0 145 % etc
129 0 145
49 0 225
129 0 225
0 129 65
0 49 65
0 129 145
0 49 145
0 129 225
0 49 225 ];

listing of calibrate.m:

% CALIBRATE
%
% Function to perform camera calibration
%
% Usage: C = calibrate(im, XYZ, uv)
%
% Where: im - is the image of the calibration target.
% XYZ - is a n x 3 array of XYZ coordinates
% of the calibration target points.
% uv - is a 2 x n array of the image coordinates
% of the calibration target points.
% C - is the 3 x 4 camera calibration matrix.
%
% This function plots the uv coordinates onto the image of
% the calibration target. It also projects the XYZ coordinates
% back into image coordinates using the calibration matrix
% and plots these points too as a visual check on the accuracy of
% the calibration process.
% Lines from the origin to the vanishing points in the X, Y and
% Z directions are overlaid on the image.
% The mean squared error between the positions of the uv coodinates
% and the projected XYZ coordinates is also reported.
%
% The function should also report the error in satisfying the
% camera calibration matrix constraint - the magnitude of
% (q1 x q3).(q2 x q3)
%

function C = calibrate(im, XYZ, uv)
% obtain rows so arbitrary number of points can be used
[rows, cols] = size(XYZ);

XYZ1 = [XYZ, ones(rows, 1)];

% buat B matrix
for n = 1:rows
B(2*n-1, :) = [XYZ1(n, :) 0 0 0 0 -uv(1, n)*XYZ(n, :)];
B(2*n, :) = [0 0 0 0 XYZ1(n, :) -uv(2, n)*XYZ(n, :)];
end

c = B \ uv(:);
c(12) = 1;
C = reshape(c,4,3)'

XYZ1 = XYZ1';
for i = 1:rows
suv(:,i) = C*XYZ1(:,i);
suv(:,i) = suv(:,i)/suv(3,i);
end

% calculate the mean squared error between the positions of the uv coodinates and suv.
mse = mean(mean((uv - suv(1:2,:)).^2));
fprintf(1, 'mean squared error is %d\n', mse);

% calculate the error in satisfying the camera calibration matrix constraint
q1 = C(1,1:3)';
q2 = C(2,1:3)';
q3 = C(3,1:3)';
error = abs(dot(cross(q1, q3), cross(q2, q3)));
fprintf(1, 'error in satisfying the camera calibration matrix is %d\n', error);

% annotate image
imshow(im);
hold;
% plot uv coordinates
plot(uv(1, :), uv(2, :),'rx')

% plot XYZ coordinates
plot(suv(1,:), suv(2,:),'b+')

% draw vanishing lines
for i = 1:3
line([C(1,4), C(1,i)/C(3,i)],[C(2,4), C(2,i)/C(3,i)])
end
hold;

Portofolio 5

2008-05-18T00:47:00.000-07:00

portofolio 5:

portofolio5 = gaussian('myPhoto.jpg',3)

portofolio5 = gaussian('myPhoto.jpg',5)

portofolio5 = gaussian('myPhoto.jpg',9)

listing of gaussian.m:

function portofolio5 = gaussian(im,deviation)

im = imread(im)
g = rgb2gray(im);

image = filter2(fspecial('gaussian', size(g),deviation), g);
figure,imagesc(image),colormap gray;

cannyim = EDGE(image,'canny',[0.0500 0.1250]);
figure,imagesc(cannyim),colormap gray;

imhori = filter2([-1 0 1], image);
figure,imagesc(imhori),colormap gray;

imvert = filter2([-1; 0; 1], image);
figure,imagesc(imvert),colormap gray;

Portofolio 4

2008-05-18T00:40:00.001-07:00

Portofolio 4:

myPhoto.jpg

portofolio4 = porto4('myPhoto.jpg')

Highboost filter membuat image menjadi blur:

high order highpass filter membuat efek seperti riak air pada image:

low order pada highpass filter mempertegas batas image:

high order lowpass filter membuat efek riak air dan mempertegas batas pada image:

low order lowpass filter menyebabkan image menjadi blur:

listing of porto4.m:

function portofolio4 = porto4(img)
im = imread(img);
im = rgb2gray(im);

lowpassffthigh = lowpassfilter(size(im), 0.05, 50);
highpassffthigh = highpassfilter(size(im), 0.05, 50);
highboostfft = highboostfilter(size(im), 0.05, 50, 0.5);
lowpassfftlow = lowpassfilter(size(im), 0.05, 1);
highpassfftlow = highpassfilter(size(im), 0.05, 1);

imfft = fft2(im);

surfl(fftshift(lowpassffthigh)), shading interp;
print -dpng lowpassffthighfreq.png
surfl(fftshift(highpassffthigh)), shading interp;
print -dpng highpassffthighfreq.png
surfl(fftshift(highboostfft)), shading interp;
print -dpng highboostfftfreq.png
surfl(fftshift(lowpassfftlow)), shading interp;
print -dpng lowpassfftlowfreq.png
surfl(fftshift(highpassfftlow)), shading interp;
print -dpng highpassfftlowfreq.png

surfl(fftshift(real(ifft2(lowpassffthigh)))), shading interp;
print -dpng lowpassffthighspac.png
surfl(fftshift(real(ifft2(highpassffthigh)))), shading interp;
print -dpng highpassffthighspac.png
surfl(fftshift(real(ifft2(highboostfft)))), shading interp;
print -dpng highboostfftspac.png
surfl(fftshift(real(ifft2(lowpassfftlow)))), shading interp;
print -dpng lowpassfftlowspac.png
surfl(fftshift(real(ifft2(highpassfftlow)))), shading interp;
print -dpng highpassfftlowspac.png

newlowpassffthigh = lowpassffthigh.*imfft;
newhighpassffthigh = highpassffthigh.*imfft;
newhighboostfft = highboostfft.*imfft;
newlowpassfftlow = lowpassfftlow.*imfft;
newhighpassfftlow = highpassfftlow.*imfft;

figure,imagesc(real(ifft2(newlowpassffthigh))),title('lowpassffthigh'),colormap gray;
figure,imagesc(real(ifft2(newhighpassffthigh))),title('highpassffthigh'),colormap gray;
figure,imagesc(real(ifft2(newhighboostfft))),title('highboostfft'),colormap gray;
figure,imagesc(real(ifft2(newlowpassfftlow))),title('lowpassfftlow'),colormap gray;
figure,imagesc(real(ifft2(newhighpassfftlow))),title('highpassfftlow'),colormap gray;
figure,imagesc(im),title('im'),colormap gray;

Portofolio 3

2008-05-17T22:42:00.001-07:00

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;

Portofolio 2

2008-05-17T22:36:00.001-07:00

myPhoto.jpg

bw = myfunc('myPhoto.jpg', 0.4);

strel = circularstruct(2);
imclose(bw,strel);

imopen(bw, strel);

tugas3_1= openclose (bw, 2)

tugas3_2 = closeopen (bw, 2)

listing of myfunc.m:

function bw = myfunc(im, thresh)
image = imread(im);
g = rgb2gray(image);
bw = im2bw(image,thresh);
imshow(bw)
imwrite(bw,'image_blacknwhite.png');

listing of tugas3_1.m:

function tugas3_1 = openclose (image, radius)
strel = circularstruct(radius);
openclose = imclose(imopen(image, strel),strel);
imshow(openclose);

listing of tugas 3-2.m:

function tugas3_2 = closeopen (image, radius)
strel = circularstruct(radius);
closeopen = imopen(imclose(image, strel),strel);
imshow(closeopen);

listing of ciscularstruct.m:

% CIRCULARSTRUCT
%
% Function to construct a circular structuring element
% for morphological operations.
%
% function strel = circularstruct(radius)
%
% Note radius can be a floating point value though the resulting
% circle will be a discrete approximation
%
% Peter Kovesi March 2000

function strel = circularstruct(radius)

if radius <>
error('radius must be >= 1');
end

dia = ceil(2*radius); % Diameter of structuring element

if mod(dia,2) == 0 % If diameter is a odd value
dia = dia + 1; % add 1 to generate a `centre pixel'
end

r = fix(dia/2);
x = ones(dia,1)*[-r:r]; % Matrix with each pixel set to its x coordinate relative to the centre
y = [-r:r]'*ones(1,dia); % " " " " " " " y " " " " "
rad = sqrt(x.^2 + y.^2); % " " " " " " " radius from the middle
strel = rad<=radius

listing of locatelandmarks.m:

% LOCATELANDMARKS - locates landmarks on SURF form
%
% Usage: [tl, tr, bl, br] = locatelandmarks(im)
%
% Argument: im - Image to be processed, assumed binary.
%
% Returns: tl, tr, bl, br
% - Coordinates of the centroids of the top-left, top-right,
% bottom-left and bottom-right landmarks respectively.
% These coordinates are returned as column vectors in the
% form [row; col] for each landmark.
%
% The function should also display the image with the centroids of the
% landmarks overlayed.

function [tl, tr, bl, br] = locatelandmarks(im)

image = imread(im)
[rows,cols] = size(image);

% menginisialisasi koordinat landmark
tl=[rows;cols];
tr=[rows;0 ];
bl=[0 ;cols];
br=[0 ;0 ];

bw = ~image;

a = 8/784; %faktor skala

SE = ones(a*rows, a*cols);
%lakukan close-open
closeopen = imopen(imclose(bw,SE),SE);

% matrix yang setiap pixelnya di set untuk setiap x koordinat
x = ones(rows,1)*[1:cols];

% matrix yang setiap pixelnya di set untuk setiap x koordinat
y = [1:rows]'*ones(1,cols);

% melabeli semua blob
[L, num] = bwlabel(closeopen);

% mencari blob untuk menentukan titik tengah
for i = 1:num
img = L==i;
area(i) = sum(sum(img));
meanx = sum(sum(double(img).*x))/area(i);
meany = sum(sum(double(img).*y))/area(i);

% tentukan kiri atas
if tl(1)+tl(2) > meanx + meany
tl = [meanx; meany];
tlindex = i;
end
% tentukan kanan atas
if tr(1)-tr(2) > meanx - meany
tr = [meanx; meany];
trindex = i;
end
% tentukan kanan bawah
if br(1)+br(2) <>
br = [meanx; meany];
brindex = i;
end
% tentukan kiri bawah
if bl(1)-bl(2) <>
bl = [meanx; meany];
blindex = i;
end
end

imshow(image)
hold;
landmarks = [tl, tr, br, bl];
plot(landmarks(1,:), landmarks(2,:),'xr','MarkerSize',25,'LineWidth',2);
plot(landmarks(1,:), landmarks(2,:),'or','MarkerSize',15,'LineWidth',1);
hold;

end

Portofolio 1

2008-05-16T10:16:00.000-07:00

1. mengubah image jadi hitam putih
bw = myfunc('lego1.png', 0.6)

2. BWL = bwlabel(~bw)
ob1 = BWL==1
ob2 = BWL==2
ob3 = BWL==3
[area, centroid, theta, roundness] = moments(ob1)
[area, centroid, theta, roundness] = moments(ob2)
[area, centroid, theta, roundness] = moments(ob3)

Listing of myfunc.m :

function bw = myfunc(im, thresh)
image = imread(im);
g = rgb2gray(image);
bw = im2bw(image,thresh);
imshow(bw)
imwrite(bw,'image_blacknwhite.png');

Listing of moments.m :

% MOMENTS
%
% Fungsi menghitung momen sebuah citra biner dan mengembalikan
% area, centroid, sudut dari sumbu inertia minimum, dan sebuah
% ukuran ‘kebulatan'. Fungsi ini mengasumsikan bahwa hanya ada
% satu obyek pada citra biner.
%
% function [area, centroid, theta, roundness] = moments(im)
%
% Argument: im - citra biner berisi nilai 0 atau 1
%
% Returns: area - luasan citra biner
% centroid - vektor dengan 2 elemen
% theta - sudut sumbu kelembaman minimum (radian)
% roundness - ratio minimum inertia/maximum inertia.
%
% Fungsi ini juga menampilkan citra dan letak dari
% centroid dan sumbu inertia minimum. Area yang dihitung
% harus ditampilkan sebagai title di citra.

function [area, centroid, theta, roundness] = moments(im)
%mendapatkan ukuran matrix image
[rows,cols] = size(im);
x = ones(rows,1)*[1:cols]; % Matrix pada setiap pixel set pada koordinat x
y = [1:rows]'*ones(1,cols); % Matrix pada setiap pixel set pada koordinat y

area = sum(sum(im));
meanx = sum(sum(double(im).*x))/area;
meany = sum(sum(double(im).*y))/area;
centroid = [meanx, meany];

% mengubah koordinat ke pusat masa
x = x - meanx;
y = y - meany;

a = sum(sum(double(im).*x.^2));
b = sum(sum(double(im).*x.*y))*2;
c = sum(sum(double(im).*y.^2));

denom = b^2 + (a-c)^2;

if denom == 0
% tetha mengubah sudut
thetamin = 2*pi*rand;
thetamax = 2*pi*rand;
roundness = 1;
else
sin2thetamin = b/sqrt(denom); %solusi positif
sin2thetamax = -sin2thetamin; %solusi negatif
cos2thetamin = (a-c)/sqrt(denom); %solusi positif
cos2thetamax = -cos2thetamin; %solusi negatif

thetamin = atan2(sin2thetamin, cos2thetamin)/2;
thetamax = atan2(sin2thetamax, cos2thetamax)/2;
Imin = 0.5*(c+a) - 0.5*(a-c)*cos2thetamin - 0.5*b*sin2thetamin;
Imax = 0.5*(c+a) - 0.5*(a-c)*cos2thetamax - 0.5*b*sin2thetamax;
roundness = Imin/Imax;
end

% menggambar sumbu sesuai dengan objek
rho = sqrt(area)/(roundness + 0.5) + 5 ;
[X1,Y1] = pol2cart(thetamin, rho);
[X2,Y2] = pol2cart(thetamin + pi, rho);

imshow(im);
hold;
line([X1 + meanx, X2 + meanx],[Y1 + meany, Y2 + meany])
plot(meanx, meany,'r+')
hold;