% MATLAB Review Homework Solutions
% Edward Hsiao
% edhsiao@cmu.edu
% September 2, 2009

%% 1) MATLAB Vectorization
X1 = rand(1,1e6);
X2 = rand(1,1e6);

% vectorized
tic;
dist1 = sqrt(sum((X1-X2).^2)); 
toc;

% not vectorized
tic;
count = 0;
for i = 1:length(X1)
    count = count + (X1(i)-X2(i)).^2;
end
dist2 = sqrt(count);
toc;


%% 2) Basic Matrix Operations
A = [1 0 5 9
     2 1 3 8
     0 7 2 6
     4 7 1 5];
 
sum(A)      % 2.1
sum(A,2)    % 2.2
A'*A        % 2.3
det(A)      % 2.4
rank(A)     % 2.5
[e,v] = eig(A)      % 2.6
A.^2                % 2.7
sqrt(A)             % 2.8
[i,j] = ind2sub(size(A),find(A==0));    % 2.9
A(A==0) = -1;       % 2.10                       


%% 3.1) Distance of points to centroid
B = [1 2 5 13 34 1 1 2 3 5
     1 3 8 21 55 8 13 21 34 55];
 
dist = sqrt(sum((B - repmat(mean(B,2),[1,size(B,2)])).^2));


%% 3.2) Projection on Line Segment
A = [1;1]; B = [5;3]; C = [2;3];

b = (B-A);
b = b / norm(b);
c = C-A;

D = dot(b,c)*b + A

figure(1);
plot([A(1) B(1)],[A(2) B(2)], 'r'); hold on;
plot(C(1),C(2),'r.');
plot(D(1),D(2),'b.');
plot([C(1) D(1)],[C(2) D(2)], 'b-'); hold on;
axis equal;
xlim([0,6]); ylim([0,4]);


%% 3.3) Linear System of Equations
A = [2 1 1;7 3 -1;3 11 5]; b = [1;5;-2];

x = inv(A)*b
% or
x = A \ b       % equivalent to x = pinv(A) * b


%% 3.4) Matching
C = [1 5 3 9 7;5 3 2 0 4];
D = [6 1 0 2 5;9 3 2 8 1];

sqX1sqY1 = repmat(sum(C.^2)',[1,size(D,2)]);    % X1^2 + Y1^2
sqX2sqY2 = repmat(sum(D.^2),[size(C,2),1]);     % X2^2 + Y2^2
X1X2Y1Y2 = C'*D;                                % X1X2 + Y1Y2

% sqrt[ (X1^2 + Y1^2) + (X2^2 + Y2^2) - 2*(X1X2 + Y1Y2) ]
dist = sqrt(sqX1sqY1 + sqX2sqY2 - 2*X1X2Y1Y2);

% find the closest match for each row
[bestdist,bestmatch] = min(dist,[],2);


%% 4) Image Stretching
im = imread('apple.jpg');
H = [1 0 0; 
     0 2 0;
     0 0 1];
figure(1); imshow(im); axis image;
im = double(im);

% get size of matrix
[R,C,D] = size(im);

% compute the corner locations
p = [1 C C 1
     1 1 R R
     1 1 1 1];

% find the warped corners in the second image
wp = H * p;

% allocate the warp image by finding maximum rows and cols needed
maxx = max(wp(1,:));
maxy = max(wp(2,:));
warpim = zeros([maxy,maxx,D]);

% compute all points in the warped image
[wX,wY] = meshgrid(1:maxx,1:maxy);
allwp = [wX(:)';wY(:)';ones(1,maxx*maxy)];

% compute the inverse warped positions in the original image
alluwp = inv(H) * allwp;

% compute the image values for each dimension
for i = 1:D
    warpimi = interp2(im(:,:,i),alluwp(1,:),alluwp(2,:));
    warpim(:,:,i) = reshape(warpimi,maxy,maxx);
end

% show the image
warpim = uint8(warpim);
figure(2);imshow(warpim); axis image;


% EASY WAY TO DO THE SAME THING
im = imread('apple.jpg');
[R,C,D] = size(im);
easywarpim = imresize(im,[2*R C]);
figure(3); imshow(easywarpim);