
    % Scene layout:
    photo = double(imread('test.jpg'));
    photo = photo(:,:,1);
    photoWidth  = size(photo,2);
    photoHeight = size(photo,1);
    
    % Scene center:
    planeZ = 100;
    sC = [0,0, planeZ];
    
    % Camera intrinsic parameters:
    f  = 100;
    Sx = 1;
    Sy = 1;
    Ox = 0;
    Oy = 0;
    
    % Camera images size:
    width  = 200;
    height = 200;
    
    c1 = [-200, -50, -150];
    c2 = [-100, -25, -300];
    c3 = [ 100,  25, -300];
    c4 = [ 200,  50, -150];
    
    d1 = sC([1,3])-c1([1,3]); d1 = d1/norm(d1,2); theta1 = acos(d1(2)); if(c1(1)<0) theta1 = -theta1; end
    d2 = sC([1,3])-c2([1,3]); d2 = d2/norm(d2,2); theta2 = acos(d2(2)); if(c2(1)<0) theta2 = -theta2; end
    d3 = sC([1,3])-c3([1,3]); d3 = d3/norm(d3,2); theta3 = acos(d3(2)); if(c3(1)<0) theta3 = -theta3; end
    d4 = sC([1,3])-c4([1,3]); d4 = d4/norm(d4,2); theta4 = acos(d4(2)); if(c4(1)<0) theta4 = -theta4; end
    
    % Show schematic of the 3D scene and cameras:
    clf; figure(50); hold on;
    plot([-photoWidth/2;photoWidth/2],[planeZ,planeZ],'r');
    plot([c1(1),sC(1)],[c1(3),sC(3)],'b')
    plot([c2(1),sC(1)],[c2(3),sC(3)],'b')
    plot([c3(1),sC(1)],[c3(3),sC(3)],'b')
    plot([c4(1),sC(1)],[c4(3),sC(3)],'b')
    
    R1 = build3DYRotation(theta1);
    R2 = build3DYRotation(theta2);
    R3 = build3DYRotation(theta3);
    R4 = build3DYRotation(theta4);

    P1 = buildCameraMatrix(R1,c1,f,Sx,Sy,Ox,Oy);
    P2 = buildCameraMatrix(R2,c2,f,Sx,Sy,Ox,Oy);
    P3 = buildCameraMatrix(R3,c3,f,Sx,Sy,Ox,Oy);
    P4 = buildCameraMatrix(R4,c4,f,Sx,Sy,Ox,Oy);
    
    % Compute the image formation matrices: 
    F1 = getImgFormationMatrix(P1,width,height,photoWidth,photoHeight,0,0,planeZ);
    F2 = getImgFormationMatrix(P2,width,height,photoWidth,photoHeight,0,0,planeZ);
    F3 = getImgFormationMatrix(P3,width,height,photoWidth,photoHeight,0,0,planeZ);
    F4 = getImgFormationMatrix(P4,width,height,photoWidth,photoHeight,0,0,planeZ);
    
    % Define the gaussian std dev for the blurring of the photo in the plane prior to picture taking:
    sigma = 1;
    
    % "Take" the pictures:
    im1 = takePictureUsingImgFormationMat(F1,width,height,photo,sigma);
    im2 = takePictureUsingImgFormationMat(F2,width,height,photo,sigma);
    im3 = takePictureUsingImgFormationMat(F3,width,height,photo,sigma);
    im4 = takePictureUsingImgFormationMat(F4,width,height,photo,sigma);
    
    % Show the original photo in planeZ:
    figure(100); showIm(photo);
    
    % Show each "taken" picture:
    figure(1); showIm(im1);
    figure(2); showIm(im2);
    figure(3); showIm(im3);
    figure(4); showIm(im4);
    
%    % Compute the average photo using all pictures:
%    sigmaForAvg = .01;
%    img1 = {im1;P1}; img2 = {im2;P2}; img3 = {im3;P3}; img4 = {im4;P4};
%    imagesList=[img1,img2,img3,img4];
%    avgIm = getAvgImgInPlane(imagesList,photoWidth,photoHeight,planeZ,sigmaForAvg);
%    
%    % Show the average photo:
%    figure(200); showIm(avgIm);
    
    % Apply my Super Resolution Method:
    
%    F1 = getImgFormationMatrix(P1,width,height,photoWidth,photoHeight,0,0,planeZ+5);
%    F2 = getImgFormationMatrix(P2,width,height,photoWidth,photoHeight,0,0,planeZ+5);
%    F3 = getImgFormationMatrix(P3,width,height,photoWidth,photoHeight,0,0,planeZ+5);
%    F4 = getImgFormationMatrix(P4,width,height,photoWidth,photoHeight,0,0,planeZ+5);
    Fs  = {F1,F2,F3,F4};
    
    ims = {im1,im2,im3,im4};
    reconstPhoto = superResolution(ims,Fs,photoWidth,photoHeight);
    
    reconsAreaWidth  = 100;
    reconsAreaHeight = 100;
    
    Dx = (size(reconstPhoto,2)-reconsAreaWidth)/2;
    Dy = (size(reconstPhoto,1)-reconsAreaHeight)/2;
    initGuess = reconstPhoto(Dy+1:Dy+reconsAreaHeight-1,Dx+1:Dx+reconsAreaWidth-1);
    
%    robSRPhoto = robustSuperResolution(ims,Fs,photoWidth,photoHeight,reconstPhoto,2,.001,3);
    robSRPhoto = robustSuperResolution(ims,Fs,reconsAreaWidth,reconsAreaHeight,initGuess,2,.001,3);
    
    % Show the Super Resolution reconstruction:
    figure(300); showIm(reconstPhoto);
    