% Given a surface triangulated by tri and mapped to pixel values by 
% polynomials associated with the coefficients in c, this function 
% evaluates these polynomials at uniformly spaced locations on an
% R*N1xR*N2 grid.  (That is, it calculates the HR pixel value at each
% grid point (x,y).)
function HRimage = estimateGridVals(tri,registeredPts,c,R,N1,N2,pixelVals)

[totalTris,three] = size(tri);

M = N1*R-1;
N = N2*R-1;
HRimage = zeros(M,N);
HRtriangles = zeros(M,N);
[X,Y] = meshgrid(1:1/R:N2,1:1/R:N1);

% Check all the triangles in order noting in which triangle each HR
% grid point occurs.
for triNum = 1:totalTris
    pts = registeredPts(tri(triNum,:),:);
    IN = inpolygon(X,Y,pts(:,1),pts(:,2));  % NxM
    HRtriangles(ind2sub(size(IN),find(IN==1))) = triNum;
end

for y = 1:M   % row
    for x = 1:N   % col
        
        % For testing, just average the pixels from the vertices of the
        % triangle the HR point is in.
        pix = pixelVals(tri(HRtriangles(y,x),:));
        HRimage(y,x) = (pix(1) + pix(2) + pix(3))/3;
        
        % Extract appropriate set of 9 c values
        %HRptC = c(HRtriangles(y,x),:);
        % Bivariate polynomial
        %HRimage(y,x) = sum(HRptC.*[1,x,y,x^2,y^2,x^3,(x^2)*y,x*(y^2),y^3]);
    end
end

end