%---------------------------------------------------------------------------------
    clear;
    rand('state',1);
    SHOW_INPUT_IMAGES = 0;
%---------------------------------------------------------------------------------
    % GENERATE THE SCENE:
    
    % Generate scene layout:
    photo = double(imread('test.jpg'));
    photo = photo(:,:,1);

    w = 150;
    h = 150;
    xs = [336,264,276,276];
    ys = [180,638,401,284];
    
    ph1 = photo(xs(1)-w/2:xs(1)+w/2, ys(1)-h/2:ys(1)+h/2);
    ph2 = photo(xs(2)-w/2:xs(2)+w/2, ys(2)-h/2:ys(2)+h/2);
    ph3 = photo(xs(3)-w/2:xs(3)+w/2, ys(3)-h/2:ys(3)+h/2);
    ph4 = photo(xs(4)-w/2:xs(4)+w/2, ys(4)-h/2:ys(4)+h/2);
    
    pc1 = [-100, 50,100];
    pc2 = [   0, 40,250];
    pc3 = [-100,-30,400];
    pc4 = [ 100,  0,420];
    pcs = [pc1;pc2;pc3;pc4];
    
    % Scene center:
    sCx = (max(pcs(:,1))+min(pcs(:,1)))/2;
    sCy = (max(pcs(:,2))+min(pcs(:,2)))/2;
    sCz = (max(pcs(:,3))+min(pcs(:,3)))/2;
    sC = [sCx,sCy,sCz];
    
    photos      = {ph1,ph2,ph3,ph4};
    photoCoords = {pc1,pc2,pc3,pc4};
    scene = {sC,photos,photoCoords};
%---------------------------------------------------------------------------------    
    % GENERATE CAMERAS AT RANDOM POSITIONS:
    
    % Camera intrinsic parameters:
    f  = 100;
    Sx = 1;
    Sy = 1;
    Ox = 0;
    Oy = 0;
    
    % Camera center ranges:
    xRange = [-500, 500];
    yRange = [ -50,  50];
    zRange = [-300,   0];
    
    nrOfCameras = 12;
    [cs,Ps] = randGenerateCameras(nrOfCameras,f,Sx,Sy,Ox,Oy,xRange,yRange,zRange,sC);
    
%    c1 = [-200, -50, -150];
%    c2 = [-100, -25, -300];
%    c3 = [ 100,  25, -300];
%    c4 = [ 200,  50, -150];
%    cs = {c1,c2,c3,c4};
%    
%    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
%    
%    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);
%    Ps = {P1,P2,P3,P4};
%---------------------------------------------------------------------------------    
    % SHOW SCHEMATIC OF THE 3D SCENE AND CAMERAS:
    
    if(SHOW_INPUT_IMAGES)
        figure(50); clf; hold on;
        
        % Show the photo planes:
        photos = scene{2};
        photoCoords = scene{3};
        for i=1:size(photos,2)
            photoWidth = size(photos{i},2);
            coords = photoCoords{i};
            x = coords(1);
            z = coords(3);
            plot([x-photoWidth/2;x+photoWidth/2],[z,z],'r');
        end
        
        % Show the camera directions:
        for i=1:size(Ps,2)
            coord = cs(i,:);
            plot([coord(1),sC(1)],[coord(3),sC(3)],'b')
        end
        
        % Show the scene photos:
        for i = 1:size(photos,2)
            figure(i*100); showIm(photos{i});
        end
    end
%---------------------------------------------------------------------------------    
    % GET THE IMAGE FORMATION MATRICES:
    
    % Camera images size:
    width  = 200;
    height = 200;
    
    for i=1:size(Ps,2)
        % Define the gaussian std dev for the blurring of the photo in the plane prior to taking the picture:
        sigma = 1;
        
        % "Take" the picture of the scene:
        ims{i} = takePicture(scene,Ps{i},width,height,sigma);
        
        % Show each "taken" picture:
        if(SHOW_INPUT_IMAGES)
            figure(i); showIm(ims{i});
        end
    end
    
%---------------------------------------------------------------------------------    
    % Compute the average photo using all pictures:
    photoWidth  = 300;
    photoHeight = 300;
    
    planeZs = [pc1(3)/2,pc1(3),pc2(3),(pc1(3)+pc2(3))/2,pc3(3),pc4(3),pc4(3)+50];
    
    sigmaForAvg = .01;
    threshold   = 1;
    
    % Partition the images set in two disjoint subsets:
    imagesCount   = size(Ps,2);
    halfImgsCount = imagesCount/2;
    for i=1:halfImgsCount;
        imagesList1{i} = {ims{i};Ps{i}};
    end
    for i=halfImgsCount+1:size(Ps,2)
        imagesList2{i-halfImgsCount} = {ims{i};Ps{i}};
    end
    
    for j = 1:size(planeZs,2)
        avgIm1 = getAvgImgInPlane(imagesList1,photoWidth,photoHeight,planeZs(j),sigmaForAvg);
        avgIm2 = getAvgImgInPlane(imagesList2,photoWidth,photoHeight,planeZs(j),sigmaForAvg);
        
        % Show the average photos:
        figure(200+3*(j-1));   showIm(avgIm1);
        figure(200+3*(j-1)+1); showIm(avgIm2);
        
        diff = abs(avgIm1 - avgIm2);
        similar = (diff<=threshold);
        
        similar = blurrImage(similar,2);
        
        avgIm = (halfImgsCount*avgIm1+(imagesCount-halfImgsCount)*avgIm2)/imagesCount;
        
        similarIm = avgIm;
        similarIm(~similar) = 0;
        
        % Show the region of similarity:
        figure(200+3*(j-1)+2); showIm(similarIm);
        
if(j==2|j==3|j==5|j==6)
        for i = 1:size(Ps,2)
            ims{i} = substractFromIm(ims{i},Ps{i},similar,0,0,planeZs(j));
        end
    end
end
%---------------------------------------------------------------------------------    
    
    % Apply my Super Resolution Method:
    threshold   = 1;
    
    photoX  = [pc1(1)  ,pc1(1),pc2(1),pc2(1),           pc3(1),pc4(1),pc4(1)];
    photoY  = [pc1(2)  ,pc1(2),pc2(2),pc2(2),           pc3(2),pc4(2),pc4(2)];
    planeZs = [pc1(3)/2,pc1(3),pc2(3),(pc1(3)+pc2(3))/2,pc3(3),pc4(3),pc4(3)+50];
    
%    photoX  = [pc1(1),pc2(1),pc3(1),pc4(1)];
%    photoY  = [pc1(2),pc2(2),pc3(2),pc4(2)];
%    planeZs = [pc1(3),pc2(3),pc3(3),pc4(3)];
    
%    reconsAreaWidth  = 150;
%    reconsAreaHeight = 150;
    reconsAreaWidth  = 100;
    reconsAreaHeight = 100;
    
    % Partition the images set in two disjoint subsets:
    imagesCount   = size(Ps,2);
    halfImgsCount = imagesCount/2;
    
    ims1 = ims(1:halfImgsCount);
    ims2 = ims(halfImgsCount+1:imagesCount);
    
    for j = 1:size(planeZs,2)
        % Obtain the formation matrices of all images for the j-th plane:
        for i = 1:halfImgsCount;
            Fs1{i} = getImgFormationMatrix(Ps{i},width,height,reconsAreaWidth,reconsAreaHeight,photoX(j),photoY(j),planeZs(j));
        end
        
        for i = halfImgsCount+1:imagesCount;
            Fs2{i-halfImgsCount} = getImgFormationMatrix(Ps{i},width,height,reconsAreaWidth,reconsAreaHeight,photoX(j),photoY(j),planeZs(j));
        end
        
        % Obtain the SR reconstruction as initial guess for the Robust SR algorithm:
        reconsPhoto1 = superResolution(ims1,Fs1,reconsAreaWidth,reconsAreaHeight);
        reconsPhoto2 = superResolution(ims2,Fs2,reconsAreaWidth,reconsAreaHeight);
        
        figure(300); showIm(reconsPhoto1);
        figure(301); showIm(reconsPhoto2);
        
        diff = abs(reconsPhoto1 - reconsPhoto2);
        similar = (diff<=threshold);
        
        similar = blurrImage(similar,2);
        
        reconsPhoto = (reconsPhoto1+reconsPhoto2)/2; % OJO!!!*****
        
        similarIm = reconsPhoto;
        similarIm(~similar) = 0;
        
        % Show the region of similarity:
        figure(302); showIm(similarIm);
        
        
%        Dx = (size(reconsPhoto,2)-reconsAreaWidth)/2;
%        Dy = (size(reconsPhoto,1)-reconsAreaHeight)/2;
%        initGuess = reconsPhoto(Dy+1:Dy+reconsAreaHeight,Dx+1:Dx+reconsAreaWidth);
%        robSRPhoto = robustSuperResolution2(ims,Fs,reconsAreaWidth,reconsAreaHeight,initGuess,2,.001,3);
        
        % Show the Super Resolution reconstruction:
        figure(300+j); showIm(reconsPhoto);
    end
    