% Description:
%
%       radmapbeach.m generates the final LDR beach image in our project
%       and produces the response curves.
%       For scenes with less than 7 exposures, repeat an image used
%       earlier.
%       The sample pixels were hand-picked.
%       A smoothness factor of 2000 was used.
%
% Input:
%
%       l - smoothness factor
%       exposures - array of exposure times (NOT LOG)
%       image* - image at a given exposure
%
% Output:
%
%       final - final LDR image
%       gRopt - R response curve
%       gGopt - G response curve
%       gBopt - B response curve
%

function [final gRopt gGopt gBopt] = radmapbeach(l,exposures,image1,image2,image3,image4,image5,image6,image7)

fix(clock)

n = 7; % n = number of images
[row col rgb] = size(image1);
B = log(exposures);
w = weight();

% Select the pixel points

prow = 50; % prow = number of row pixels
pcol = 50; % pcol = number of column pixels
p = prow*pcol; % p = number of pixels used to generate g

im1 = double(image1(201:250,881:930,:));
im1 = reshape(im1,prow*pcol,rgb);

im2 = double(image2(201:250,881:930,:));
im2 = reshape(im2,prow*pcol,rgb);

im3 = double(image3(201:250,881:930,:));
im3 = reshape(im3,prow*pcol,rgb);

im4 = double(image4(201:250,881:930,:));
im4 = reshape(im4,prow*pcol,rgb);

im5 = double(image5(201:250,881:930,:));
im5 = reshape(im5,prow*pcol,rgb);

im6 = double(image6(201:250,881:930,:));
im6 = reshape(im6,prow*pcol,rgb);

im7 = double(image7(201:250,881:930,:));
im7 = reshape(im7,prow*pcol,rgb);

% Compute g for each color

im1R = im1(:,1);
im1G = im1(:,2);
im1B = im1(:,3);

im2R = im2(:,1);
im2G = im2(:,2);
im2B = im2(:,3);

im3R = im3(:,1);
im3G = im3(:,2);
im3B = im3(:,3);

im4R = im4(:,1);
im4G = im4(:,2);
im4B = im4(:,3);

im5R = im5(:,1);
im5G = im5(:,2);
im5B = im5(:,3);

im6R = im6(:,1);
im6G = im6(:,2);
im6B = im6(:,3);

im7R = im7(:,1);
im7G = im7(:,2);
im7B = im7(:,3);

ZR = zeros(p,n);
ZR(:,1) = im1R;
ZR(:,2) = im2R;
ZR(:,3) = im3R;
ZR(:,4) = im4R;
ZR(:,5) = im5R;
ZR(:,6) = im6R;
ZR(:,7) = im7R;
[gR lER] = gsolve(ZR,B,l,w);

ZG = zeros(p,n);
ZG(:,1) = im1G;
ZG(:,2) = im2G;
ZG(:,3) = im3G;
ZG(:,4) = im4G;
ZG(:,5) = im5G;
ZG(:,6) = im6G;
ZG(:,7) = im7G;
[gG lEG] = gsolve(ZG,B,l,w);

ZB = zeros(p,n);
ZB(:,1) = im1B;
ZB(:,2) = im2B;
ZB(:,3) = im3B;
ZB(:,4) = im4B;
ZB(:,5) = im5B;
ZB(:,6) = im6B;
ZB(:,7) = im7B;
[gB lEB] = gsolve(ZB,B,l,w);

% Compute the radiance map for each color

im1  = double(reshape(image1,row*col,rgb));
im1R = im1(:,1);
im1G = im1(:,2);
im1B = im1(:,3);

im2  = double(reshape(image2,row*col,rgb));
im2R = im2(:,1);
im2G = im2(:,2);
im2B = im2(:,3);

im3  = double(reshape(image3,row*col,rgb));
im3R = im3(:,1);
im3G = im3(:,2);
im3B = im3(:,3);

im4  = double(reshape(image4,row*col,rgb));
im4R = im4(:,1);
im4G = im4(:,2);
im4B = im4(:,3);

im5  = double(reshape(image5,row*col,rgb));
im5R = im5(:,1);
im5G = im5(:,2);
im5B = im5(:,3);

im6  = double(reshape(image6,row*col,rgb));
im6R = im6(:,1);
im6G = im6(:,2);
im6B = im6(:,3);

im7  = double(reshape(image7,row*col,rgb));
im7R = im7(:,1);
im7G = im7(:,2);
im7B = im7(:,3);

ZR = zeros(row*col,n);
ZR(:,1) = im1R;
ZR(:,2) = im2R;
ZR(:,3) = im3R;
ZR(:,4) = im4R;
ZR(:,5) = im5R;
ZR(:,6) = im6R;
ZR(:,7) = im7R;
             
ZG = zeros(row*col,n);
ZG(:,1) = im1G;
ZG(:,2) = im2G;
ZG(:,3) = im3G;
ZG(:,4) = im4G;
ZG(:,5) = im5G;
ZG(:,6) = im6G;
ZG(:,7) = im7G;

ZB = zeros(row*col,n);
ZB(:,1) = im1B;
ZB(:,2) = im2B;
ZB(:,3) = im3B;
ZB(:,4) = im4B;
ZB(:,5) = im5B;
ZB(:,6) = im6B;
ZB(:,7) = im7B;

Ecur = 0;

for i = 1:.2:1,
    gRtemp = i*gR;
    for j = 1:.2:1,
        gGtemp = j*gG;
        for k = 1:.2:1,
            gBtemp = k*gB;

            lER = map(ZR,B,w,gRtemp);
            lEG = map(ZG,B,w,gGtemp);
            lEB = map(ZB,B,w,gBtemp);
            lE = zeros(row,col,rgb); % Note the form lE takes for this function

            for x = 1:1:row,
                for y = 1:1:col,
                    lE(x,y,1) = lER(x+(y-1)*row);
                    lE(x,y,2) = lEG(x+(y-1)*row);
                    lE(x,y,3) = lEB(x+(y-1)*row);
                end
            end
            
            Max = max(lE(:));
            Min = min(lE(:));
            temp = uint8(round((lE-Min)*255/(Max-Min)));
            
            if Ecur < entropyC(temp)
                gRopt = gRtemp;
                gGopt = gGtemp;
                gBopt = gBtemp;
                final = temp;
                Ecur = entropyC(temp);
            end

        end
    end
end

fix(clock)