function output = panorama(imagepath)
% Creates a panorama image from a sequence of images in `imagepath' directory.
% Assumes images are in order and are the only files in this directory. 
% Note:  This skeleton code does NOT include all steps or all details required.

% Process input imagepath
files = dir(fullfile(imagepath, '*.JPG'));

% Eliminate . and .. and assume everything else in directory is an input image. 
% Also assume that all images are color. 
imagelist = files(1:end);

for i=1:length(imagelist)-1
    image1 = strcat(imagepath,'/',imagelist(i).name);
    image2 = strcat(imagepath,'/',imagelist(i+1).name);
    
    % Find matching feature points between current two images using SIFT features and SIFT algorithm for matching. 
    % Note that this code does NOT include the use of RANSAC to find and use only good matches.
    [~, matchIndex, loc1, loc2] = match(image1, image2);
    im1_ftr_pts = loc1(find(matchIndex > 0), 1:2);
    im2_ftr_pts = loc2(matchIndex(find(matchIndex > 0)), 1:2);

    % Calculate 3x3 homography matrix, H, mapping coordinates in image2 into coordinates in image1
    H = calcH(im1_ftr_pts, im2_ftr_pts);
    H_list(i) = {H};
    % Need this done backwards for H_map later
    H = calcH(im2_ftr_pts, im1_ftr_pts);
    H_list2(i) = {H};
  
end


% Select one input image as the reference image, ideally one in the middle of the sequence of input images so that there is less distortion in the output. 
% Compute new homographies that map every other image *directly* to the reference image by composing H matrices in H_list
% Save these new homographies in a list called H_map
% Hence H_map is a list of 3 x 3 homography matrices that map each image into the reference image's coordinate system. 
% The homography in H_map that is associated with the reference image should be the identity matrix, created using eye(3)
% The homographies in H_map for the other images (both before and after the reference image) are computed by using already defined matrices in H_map and H_list as described in the homework.  Note:  Composing A with B is not the same as composing B with A.
%Notice H_map and H_list are cell arrays, which are general containers in Matlab. See http://www.mathworks.com/help/matlab/matlab_prog/what-is-a-cell-array.html for more info on using cell arrays.

H_map = calcHMap(H_list, H_list2, length(imagelist));

% Compute size of output panorama image
min_row = 1;
min_col = 1;
max_row = 0;
max_col = 0;

% for each input image
for i=1:length(H_map)
    cur_image = imread(strcat(imagepath,'/',imagelist(i).name));
    [rows,cols,~] = size(cur_image);
    
    % create a matrix with the coordinates of the four corners of the current image
    pt_matrix = cat(3, [1,1,1]', [1,cols,1]', [rows, 1,1]', [rows,cols,1]');
    
    % Map each of the 4 corner's coordinates into the coordinate system of the reference image
    for j=1:4
        result = H_map{i}*pt_matrix(:,:,j);
    
        min_row = floor(min(min_row, result(1)));
        min_col = floor(min(min_col, result(2)));
        max_row = ceil(max(max_row, result(1)));
        max_col = ceil(max(max_col, result(2))); 
    end
    
end

% Calculate output image size
im_rows = max_row - min_row + 1;
im_cols = max_col - min_col + 1;

% Calculate offset of the upper-left corner of the reference image relative to the upper-left corner of the output image
row_offset = 1 - min_row;
col_offset = 1 - min_col;

% Initialize output image & alpha to black (0)
pan_image = zeros(im_rows, im_cols, 3);
pan_alpha = zeros(im_rows, im_cols, 1);

% Perform inverse mapping for each input image
for i=1:length(H_map)
    
    cur_image = im2double(imread(strcat(imagepath,'/',imagelist(i).name)));
    
    % Initialize alpha values
    [rows, cols, z] = size(cur_image);
    cur_alpha = zeros(rows, cols);
    
    % Calculate alpha values for feathering
    for row=1:rows
        if row <= rows/2
            cur_alpha(row, :) = row/(rows/2);
        else
            cur_alpha(row, :) = 2 - row/(rows/2);
        end
    end
    
    % Multiply images by alpha values prior to mapping
    cur_image(:,:,1) = cur_image(:,:,1).*cur_alpha(:,:);
    cur_image(:,:,2) = cur_image(:,:,2).*cur_alpha(:,:);
    cur_image(:,:,3) = cur_image(:,:,3).*cur_alpha(:,:);
    
    % Create a list of all pixels' coordinates in output image
    [x,y] = meshgrid(1:im_cols, 1:im_rows);
    % Create list of all row coordinates and column coordinates in separate vectors, x and y, including offset
    x = reshape(x,1,[]) - col_offset;
    y = reshape(y,1,[]) - row_offset;
    
    % Create homogeneous coordinates for each pixel in output image
    pan_pts(1,:) = y;
    pan_pts(2,:) = x;
    pan_pts(3,:) = ones(1,size(pan_pts,2));
    
    % Perform inverse warp to compute coordinates in current input image
    image_coords = H_map{i}\pan_pts;
    row_coords = reshape(image_coords(1,:),im_rows, im_cols);
    col_coords = reshape(image_coords(2,:),im_rows, im_cols);
    % Note:  Some values will return as NaN ("not a number") because they map to points outside the domain of the input image
    
    % Bilinear interpolate color values
    pixel_color_r = interp2(cur_image(:,:,1), col_coords, row_coords, 'linear', 0);
    pixel_color_g = interp2(cur_image(:,:,2), col_coords, row_coords, 'linear', 0);
    pixel_color_b = interp2(cur_image(:,:,3), col_coords, row_coords, 'linear', 0);
    
    % Map alpha to panorama
    pixel_alpha = interp2(cur_alpha(:,:), col_coords, row_coords, 'linear', 0);
    
    % Accumulate the color values
    pan_image(:,:,1) = pan_image(:,:,1) + pixel_color_r;
    pan_image(:,:,2) = pan_image(:,:,2) + pixel_color_g;
    pan_image(:,:,3) = pan_image(:,:,3) + pixel_color_b;
    
    % Accumulate alpha values
    pan_alpha = pan_alpha + pixel_alpha;

end
    %figure();
    %imshow(pan_image);
    %figure();
    %imshow(pan_alpha);
    [rows,cols] = size(pixel_alpha);
    
    for row=1:rows
        for col=1:cols
            if pixel_alpha(row, col) > 0
                pan_image(row, col, 1) = pan_image(row, col, 1) ./ pan_alpha(row, col);
                pan_image(row, col, 2) = pan_image(row, col, 2) ./ pan_alpha(row, col);
                pan_image(row, col, 3) = pan_image(row, col, 3) ./ pan_alpha(row, col);
            end
        end
    end
    output = pan_image;
    figure();
    imshow(output);
    
end