%---------------------------------------------------------------------------------
    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;
    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);
%    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 Super Resolution Method:
    threshold       = 1;
    estimatorStdDev = 10;
    weightThreshold = .001;
    
%    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  = 150;
%    reconsAreaHeight = 150;
    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);

%        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,Ws] = robustSuperResolution(ims,Fs,reconsAreaWidth,reconsAreaHeight,initGuess,estimatorStdDev,.001,3);
        
        % Show the Super Resolution reconstruction:
        figure(301); showIm(robSRPhoto);
    
        % Truncate Ws up to a threshold:
        for k = 1:size(Ws,2)
            Ws2{k} = (Ws{k}>weightThreshold);
        end
        
        sigmaForIntersect = 1;
        intersecIm = getProjectionIntersection(Ws2,Ps,reconsAreaWidth,reconsAreaHeight,planeZs(j),sigmaForIntersect);
        avgIm      = getAvgImgInPlane(Ws,Ps,reconsAreaWidth,reconsAreaHeight,planeZs(j),sigmaForIntersect);
        
%        % Show the intersection image:
%        figure(302); showIm(intersecIm);
%        
%        % Show the average image:
%        figure(303); showIm(avgIm);
%        
%        for k = 1:size(Fs,2)
%            tt = .006;
%            av = avgIm>tt;
%            backProj = Fs{k} * av(:);
%            backProj = reshape(backProj,200,200);
%            figure(304); showIm(backProj);
%            
%            backProj2 = Fs{k} * intersecIm(:);
%            backProj2 = reshape(backProj2,200,200);
%            figure(305); showIm(backProj2);
%        end
    end
    