function [F] = getOcclusionImgFormationMatrix(Fold,similar,P,width,height,photoWidth,photoHeight,photoX,photoY,planeZ)
    
    sizeImage = width*height;
    sizePhoto = photoWidth*photoHeight;
    
    % This is the amount of pixels in the photo that map to a single pixel in the image.
    % It is related to the distance between the camera with matrix P and the planeZ.
    % The larger is such distance, the greater this value has to be.
    photoPixelsPerImagePixel = 50;
    
    [x,y] = find(similar(:));
    XYZW = [x(:)'; y(:)'];
    xyzw = [XYZW; repmat(planeZ,1,sizePhoto); repmat(1,1,sizePhoto)];
    xyzw(1,:) = (xyzw(1,:)-0.5) - photoWidth/2  + photoX;
    xyzw(2,:) = (xyzw(2,:)-0.5) - photoHeight/2 + photoY;
    
    xy = P * xyzw;
    
    xy(1,:) = floor(1 + xy(1,:)./xy(3,:)+width/2);
    xy(2,:) = floor(1 + xy(2,:)./xy(3,:)+height/2);
    
    % Compute the image coords mapping for each plane (x,y) coord:
    xyIdx   = xy(2,:)   + (xy(1,:)-1)*height;
    xyzwIdx = XYZW(2,:) + (XYZW(1,:)-1)*photoHeight;
    
    % Compute the index of the coords that fall inside of the image:
    insid = (1<=xy(1,:) & xy(1,:)<=width & 1<=xy(2,:) & xy(2,:)<=height);
    
    % Get only the coords of the pixels inside the image where the photo pixels project:
    xyIdxInside = xyIdx(insid);
    
    % Get only those coords of the photo that project inside the image
    xyzwIdxInside = xyzwIdx(insid);

    % Build the Image Formation Matrix:
    F = sparse(xyIdxInside,xyzwIdxInside,1,sizeImage,sizePhoto,photoPixelsPerImagePixel*sizeImage);
    
    % Normalize the rows of F to add up to 1 (by diving each row by the sum of its column values):
    rowSums = full(sum(F,2));
    rowSums(rowSums==0) = 1;
    F = sparse(1:sizeImage,1:sizeImage,1./rowSums) * F;
