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;