% Purpose: Given a set of low-resolution (LR) images and a desired
%          resolution factor, register and interpolate to create a
%          high-resolution (HR) image.
%
% Assuming K LR images:
%   pictures = N1xN2xK matrix with the pixel values of all the images
%            = [N1xN2 matrix of image 1's pixels; ... ;
%               N1xN2 matrix of image K's pixels]
% For each image in pictures:
%   N1 = # of pixels in the y-dimension
%   N2 = # of pixels in the x-dimension
% R = desired resolution enhancement factor
function HRimage = LRtoHR(pictures,K,R)

% Extract dimensions of images (N1xN2).
[N1,N2,K] = size(pictures);
totalPoints = K*N1*N2;

% Create transform matrices for all the images and apply these transforms
% to come up with new (x,y) locations for all the pixel values.
%   -regPtsX is an N1xN2 matrix of a picture's pixels' new x locations
%   -registeredPts combines all regPts matrices (for all pictures):
%       registeredPts = K*N1*N2x2 matrix,
%                       [xPix1, yPix1; ... ;
%                        xPixK*N1*N2; yPixK*N1*N2]
%                       where xPix1 = x location (column) of pixel 1
%   -pixelVals = pixel values matrix corresponding to registeredPts
% Use P1 as the base case and transform other images relative to it.
P1 = pictures(:,:,1);
regPtsX = zeros(N1,N2,K);
regPtsY = zeros(N1,N2,K); 
registeredPts = zeros(K*N1*N2,2);
pixelVals = zeros(K*N1*N2,1);
for row = 1:N1
    for col = 1:N2
        pixNum = (row-1)*N2 + col;
        pixelVals(pixNum,1) = pictures(row,col,1);
        registeredPts(pixNum,:) = [col,row];
    end
end
for imgNum = 2:K
    % 4 = number of levels, 2 = iterations per level
    [regPtsX,regPtsY] = ...
        registerImgs(proj_reg(P1,pictures(:,:,imgNum),4,2),N1,N2);
    for row = 1:N1
        for col = 1:N2
            pixNum = (imgNum-1)*N1*N2 + (row-1)*N2 + col;
            pixelVals(pixNum,1) = pictures(row,col,imgNum);
            registeredPts(pixNum,:) = [regPtsX(row,col),regPtsY(row,col)];
        end
    end
end

% Generate triangulation: delaunay(x,y)
% where (x,y) are pixel locations from vertices of registered LR images,
% and tri = Lx3 matrix, where each row defines 1 triangle
%           (the actual numbers in tri are indicies into x and y)
tri = delaunay(registeredPts(:,1)',registeredPts(:,2)',{'QJ'});
figure, triplot(tri,registeredPts(:,1)',registeredPts(:,2)')

% Estimate the gradient vector at each vertex
%vGradientVecs = getGradientVecs(tri,registeredPts,pixelVals);

% Approximate each triangle patch by a continuous function
%   1. Find the 9 c values for each patch (see a1 and b1 below)
%c = findBPolyCoefficients(tri,registeredPts,pixelVals,vGradientVecs);
%   2. For each point on the HR grid, find the associated triangle patch,
%      extract its c values, and use these values as the coefficients
%      in a bivariate polynomial to calculate the HR pixel value at
%      each grid point (x,y)
HRimage = estimateGridVals(tri,registeredPts,R,N1,N2,pixelVals);
imshow(HRimage,[]);

end