% Returns the image resulting from averaging the projection of all images in imageList onto the plane determined by the 3D points tl, tr, br, bl using pixXCount pixels in the X direction and pixYCount pixels in the Y direction.
function [avgIm] = getAvgImgInPlane(ims,Ps,pixXCount,pixYCount,planeZ,sigma)
    
    avgIm = zeros(pixYCount,pixXCount);
    count = zeros(pixYCount,pixXCount);
    
    [x,y] = meshgrid(1:pixXCount, 1:pixYCount);
    xyzw = [x(:)'; y(:)'; repmat(planeZ,1,pixXCount*pixYCount); repmat(1,1,pixXCount*pixYCount)];
    xyzw(1,:) = (xyzw(1,:)-0.5) - pixXCount/2;
    xyzw(2,:) = (xyzw(2,:)-0.5) - pixYCount/2;
    
    % Go over all images:
    for i = 1:size(ims,2)
        im = blurrImage(ims{i},sigma);
        P = Ps{i};
        
        width  = size(im,2);
        height = size(im,1);
        
        pix = P * xyzw;
        
        pix(1,:) = 1 + floor(pix(1,:)./pix(3,:)+width/2);
        pix(2,:) = 1 + floor(pix(2,:)./pix(3,:)+height/2);
        
        % Compute the index of the coords that fall inside of the image:
        insid = (1<=pix(1,:) & pix(1,:)<=width & 1<=pix(2,:) & pix(2,:)<=height);

        % Compute the image coords mapping for each plane (x,y) coord:
        idx = pix(2,:)+(pix(1,:)-1)*height;

        % We are only interested in the plane coords that map inside the image:
        idxInside = idx(insid);
        imInside  = im(idxInside);

        % Set the average image values to be the corresponding image pixel values:
        avgIm(insid) = avgIm(insid) + imInside;
        
        % Only count the pixels in the plane that are projected inside the image and that do not project to a null pixel:
%        nonZero = im(idx);
%        nonZero(~insid) = 0;
%        count(nonZero) = count(nonZero) + 1;
        count(insid) = count(insid) + 1;
    end
    
    count(count==0) = 1;
    avgIm = avgIm ./ count;
