%---------------------------------------------------------------------------------
    clear;
    rand('state',200);
    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);
%---------------------------------------------------------------------------------    
    % 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);
%    thirdImgsCount = imagesCount/3;
%    for i=1:thirdImgsCount;
%        ims1{i} = ims{i};
%        Ps1{i}  = Ps{i};
%    end
%    for i=thirdImgsCount+1:2*thirdImgsCount
%        ims2{i-thirdImgsCount} = ims{i};
%        Ps2{i-thirdImgsCount}  = Ps{i};
%    end
%    for i=2*thirdImgsCount+1:size(Ps,2)
%        ims3{i-2*thirdImgsCount} = ims{i};
%        Ps3{i-2*thirdImgsCount}  = Ps{i};
%    end
%    
%    for j = 1:size(planeZs,2)
%        avgIm1 = getAvgImgInPlane(ims1,Ps1,photoWidth,photoHeight,planeZs(j),sigmaForAvg);
%        avgIm2 = getAvgImgInPlane(ims2,Ps2,photoWidth,photoHeight,planeZs(j),sigmaForAvg);
%        avgIm3 = getAvgImgInPlane(ims3,Ps3,photoWidth,photoHeight,planeZs(j),sigmaForAvg);
%        
%        % Show the average photos:
%        figure(201); showIm(avgIm1);
%        figure(202); showIm(avgIm2);
%        figure(203); showIm(avgIm3);
%        
%        diff1 = abs(avgIm1 - avgIm2);
%        diff2 = abs(avgIm2 - avgIm3);
%        diff3 = abs(avgIm1 - avgIm3);
%        
%        similar = (diff1<=threshold & diff2<=threshold & diff3<=threshold);
%        
%        similar = blurrImage(similar,3);
%        
%        avgIm = (thirdImgsCount*avgIm1+thirdImgsCount*avgIm2+(imagesCount-2*thirdImgsCount)*avgIm3)/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 ROBUST SUPER RESOLUTION METHOD:
    threshold       = 1;
    estimatorStdDev = 10;
    
%    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)-5,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  = 300;
    reconsAreaHeight = 300;
    
    % 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);
        
%for i = 1:imagesCount
%    Fs{i} = getImgFormationMatrix(Ps{i},width,height,reconsAreaWidth,reconsAreaHeight,photoX(j),photoY(j),planeZs(j));
%end
%reconsPhoto = superResolution(ims,Fs,reconsAreaWidth,reconsAreaHeight);
%figure(300); showIm(reconsPhoto);
        
        Dx1 = (size(reconsPhoto1,2)-reconsAreaWidth)/2;
        Dy1 = (size(reconsPhoto1,1)-reconsAreaHeight)/2;
        initGuess1 = reconsPhoto1(Dy1+1:Dy1+reconsAreaHeight,Dx1+1:Dx1+reconsAreaWidth);
        
        Dx2 = (size(reconsPhoto2,2)-reconsAreaWidth)/2;
        Dy2 = (size(reconsPhoto2,1)-reconsAreaHeight)/2;
        initGuess2 = reconsPhoto2(Dy2+1:Dy2+reconsAreaHeight,Dx2+1:Dx2+reconsAreaWidth);
        
        [robSRPhoto1,Ws1] = robustSuperResolution(ims1,Fs1,reconsAreaWidth,reconsAreaHeight,initGuess1,estimatorStdDev,.001,3);
        [robSRPhoto2,Ws2] = robustSuperResolution(ims2,Fs2,reconsAreaWidth,reconsAreaHeight,initGuess2,estimatorStdDev,.001,3);
        
        % Show the Super Resolution reconstruction:
        figure(302); showIm(robSRPhoto1);
        figure(303); showIm(robSRPhoto2);
    
        rSRP1 = blurrImage(robSRPhoto1,4);
        rSRP2 = blurrImage(robSRPhoto2,4);
        diff = abs(rSRP1 - rSRP2);
        similar = (diff<=threshold);
        
        similar = blurrImage(similar,2);
        
        robSRPhoto = (robSRPhoto1+robSRPhoto2)/2; % OJO!!!*****
        similarIm = robSRPhoto;
        similarIm(~similar) = 0;
        
        % Show the region of similarity:
        figure(304); showIm(similarIm);
        
        for i = 1:halfImgsCount
%            aa = reshape(Fs1{i} * similarIm(:),200,200);
%            figure(305); showIm(aa);
            
            newIm = substractFromIm(ims{i},Fs1{i},similar,photoX(j),photoY(j),planeZs(j));
            
%            figure(306); showIm(ims1{i});
%            figure(307); showIm(newIm);
            
            ims1{i} = newIm;
        end
        
        for i = halfImgsCount+1:imagesCount;
%            aa = reshape(Fs2{i-halfImgsCount} * similarIm(:),200,200);
%            figure(305); showIm(aa);
            
            newIm = substractFromIm(ims{i},Fs2{i-halfImgsCount},similar,photoX(j),photoY(j),planeZs(j));
            
%            figure(306); showIm(ims{i});
%            figure(307); showIm(newIm);
            
            ims2{i-halfImgsCount} = newIm;
        end
    end