%--------------------------------------------------------------------
% Function : rotate_image
%
% This function takes a specified image, and rotates it about the 
% x-axis, then y-axis, then z-axis based on the degrees given. It also 
% scales the image based on the scale argument by projecting it
% closer or further away
%
% Input:
%     pic           matrix of original image bmp (face)
%     x_deg         degrees to rotate about the x-axis
%     y_deg         degrees to rotate about the y-axis
%     z_deg         degrees to rotate about the z-axis
%     scale         Scale this image
%                   Value should be between 400 and 1100
%
% Output:
%     newimage      returns the source image rotated and scaled
%
% Example:
%     newimage = rotate_image('faces.bmp', 5, 40, -4, 800);
%
% Info:
%     Henry Lin
%     01/15/2008 for 600.161.01
%--------------------------------------------------------------------

function [newimage] = rotate_image(pic_path, x_deg, y_deg, z_deg, scale)

% read in image 
pic = imread(pic_path);

% convert degrees to radians
x_rad = (pi * x_deg)/180;
y_rad = (pi * y_deg)/180;
z_rad = (pi * z_deg)/180;

% get size of image
image_size = size(pic);

% set the origin of the image coord sys to the middle of the image
origin = [ceil(image_size(1)/2) ceil(image_size(2)/2)]

% create a 4-dim matrix that stores the coord(x,y,z,1) location values
coords = zeros(image_size(1), image_size(2), 4);
coords(:,:,4) = ones(image_size(1), image_size(2));

% fill in the location values in the coord matrix
for x=1:image_size(2)
  coords(:,x,1) = x - origin(2);
end

for y=1:image_size(1)
  coords(y,:,2) = y - origin(1);
end

x_coords = coords(:,:,1);
y_coords = coords(:,:,2);
z_coords = coords(:,:,3);
h_coords = coords(:,:,4);

x_coords = x_coords(:)';
y_coords = y_coords(:)';
z_coords = z_coords(:)';
h_coords = h_coords(:)';

locs = [x_coords; y_coords; z_coords; h_coords];

size_locs = size(locs);

% rotation matrices
rotate_x = [1 0 0 0; 0 cos(x_rad) sin(x_rad) 0; 0 -sin(x_rad) cos(x_rad) 0; 0 0 0 1]
rotate_y = [cos(y_rad) 0 -sin(y_rad) 0; 0 1 0 0 ;sin(y_rad) 0 cos(y_rad) 0; 0 0 0 1]
rotate_z = [cos(z_rad) sin(z_rad) 0 0; -sin(z_rad) cos(z_rad) 0 0; 0 0 1 0; 0 0 0 1]

% rotation matrix to multiply by
rotation = rotate_z * rotate_y * rotate_x 

% new coordinates
newlocs = rotation * locs;

newlocs_normx = [];
newlocs_normy = [];
newlocs_normz = [];

newlocs(3,:) = max(2* image_size) + newlocs(3, :);

size_newlocs = size(newlocs);

for rows=1:image_size(1)

  % project 3-d pixels onto 2d image plane. Where the image plane is
  % determines scale of image
  newlocs_normx = [newlocs_normx; (newlocs(1, rows:image_size(1):end) .* scale) ./ newlocs(3, rows:image_size(1):end)];
  newlocs_normy = [newlocs_normy; (newlocs(2, rows:image_size(1):end) .* scale) ./ newlocs(3, rows:image_size(1):end)];

end

newlocs_norm(:,:,1) = int16(newlocs_normx);
newlocs_norm(:,:,2) = int16(newlocs_normy);

min_y = int16(abs(min(newlocs_normy(:))) + 1);
min_x = int16(abs(min(newlocs_normx(:))) + 1);

y_min = min(newlocs_normy(:))
y_max = max(newlocs_normy(:))
x_min = min(newlocs_normx(:))
x_max = max(newlocs_normx(:))

image_size;
size(pic);

% calculated new pixel values for each pixel coord
for y=1:image_size(1)
  for x=1:image_size(2)

 	newimage(newlocs_norm(y,x,2)+min_y, newlocs_norm(y,x,1)+min_x,1) = pic(y,x,1);
 	newimage(newlocs_norm(y,x,2)+min_y, newlocs_norm(y,x,1)+min_x,2) = pic(y,x,2);
 	newimage(newlocs_norm(y,x,2)+min_y, newlocs_norm(y,x,1)+min_x,3) = pic(y,x,3);
 
  end
end

size(newimage);
imagesc(newimage)
axis image