%
% CameraGamma.m
%
% Read in camera output of GrayPatch and calculate camera Gamma Correction
%
% Hui Li
% 3/7/2000
%	

% read in camera image
cim = imread('./image/GrayPatchOutput.jpg','jpg');
figure
imshow(cim);

im = double(cim)./255;
figure
imshow(im);


% Decide test area
p = zeros(210,300);
cim(107:316,105:404) = p+255;
cim(107:316,500:799) = p+255;
cim(107:316,880:1179) = p+255;

cim(392:601,105:404) = p+255;
cim(392:601,500:799) = p+255;
cim(392:601,880:1179) = p+255;

cim(692:901,110:409) = p+255;
cim(692:901,500:799) = p+255;
cim(692:901,880:1179) = p+255;

figure
imshow(cim);


% Retrieve test patch
im1 = im(107:316,105:404,:);
im2 = im(107:316,500:799,:);
im3 = im(107:316,880:1179,:);

im4 = im(392:601,105:404,:);
im5 = im(392:601,500:799,:);
im6 = im(392:601,880:1179,:);

im7 = im(692:901,110:409,:);
im8 = im(692:901,500:799,:);
im9 = im(692:901,880:1179,:);

imr = [im1 im2 im3; im4 im5 im6; im7 im8 im9];
figure
imshow(imr);


GPval = [0.0078125 0.125:0.125:1];

%%%% Compute camera RED gamma
cval(1) = mean(mean(im1(:,:,1)));
cval(2) = mean(mean(im2(:,:,1)));
cval(3) = mean(mean(im3(:,:,1)));
cval(4) = mean(mean(im4(:,:,1)));
cval(5) = mean(mean(im5(:,:,1)));
cval(6) = mean(mean(im6(:,:,1)));
cval(7) = mean(mean(im7(:,:,1)));
cval(8) = mean(mean(im8(:,:,1)));
cval(9) = mean(mean(im9(:,:,1)));

a = max(max(cval));
sval = cval./a;
save('sval.mat','sval');

figure
plot(GPval, sval, 'ro-');
grid on, hold on;

% approimate gamma
x = fmin('f', 0, 3, sqrt(2));
idealVal = GPval.^x;
gammaR = 1/x;

plot(GPval,idealVal,'k*-');
legend('Actual Data', 'Estimated Data', 2);
title('Camera Gamma Function of R Value');
xlabel('Input Values');
ylabel('Gamma Corrected R Values');
hold off;


%%%% Compute camera Green gamma
cval(1) = mean(mean(im1(:,:,2)));
cval(2) = mean(mean(im2(:,:,2)));
cval(3) = mean(mean(im3(:,:,2)));
cval(4) = mean(mean(im4(:,:,2)));
cval(5) = mean(mean(im5(:,:,2)));
cval(6) = mean(mean(im6(:,:,2)));
cval(7) = mean(mean(im7(:,:,2)));
cval(8) = mean(mean(im8(:,:,2)));
cval(9) = mean(mean(im9(:,:,2)));

a = max(max(cval));
sval = cval./a;
save('sval.mat','sval');

figure
plot(GPval, sval, 'go-');
grid on, hold on;

% approimate gamma
x = fmin('f', 0, 3, sqrt(2));
idealVal = GPval.^x;
gammaG = 1/x;

plot(GPval,idealVal,'k*-');
legend('Actual Data', 'Estimated Data', 2);
title('Camera Gamma Function of G Value');
xlabel('Input Values');
ylabel('Gamma Corrected G Values');
hold off;


%%%% Compute camera Blue gamma
cval(1) = mean(mean(im1(:,:,3)));
cval(2) = mean(mean(im2(:,:,3)));
cval(3) = mean(mean(im3(:,:,3)));
cval(4) = mean(mean(im4(:,:,3)));
cval(5) = mean(mean(im5(:,:,3)));
cval(6) = mean(mean(im6(:,:,3)));
cval(7) = mean(mean(im7(:,:,3)));
cval(8) = mean(mean(im8(:,:,3)));
cval(9) = mean(mean(im9(:,:,3)));

a = max(max(cval));
sval = cval./a;
save('sval.mat','sval');

figure
plot(GPval, sval, 'bo-');
grid on, hold on;

% approimate gamma
x = fmin('f', 0, 3, sqrt(2));
idealVal = GPval.^x;
gammaB = 1/x;

plot(GPval,idealVal,'k*-');
legend('Actual Data', 'Estimated Data', 2);
title('Camera Gamma Function of B Value');
xlabel('Input Values');
ylabel('Gamma B Corrected Values');
hold off;


%%%% Compute overall camera gamma
cval(1) = mean(mean(mean(im1)));
cval(2) = mean(mean(mean(im2)));
cval(3) = mean(mean(mean(im3)));
cval(4) = mean(mean(mean(im4)));
cval(5) = mean(mean(mean(im5)));
cval(6) = mean(mean(mean(im6)));
cval(7) = mean(mean(mean(im7)));
cval(8) = mean(mean(mean(im8)));
cval(9) = mean(mean(mean(im9)));

a = max(max(cval));
sval = cval./a;
save('sval.mat','sval');

figure
plot(GPval, sval, 'mo-');
grid on, hold on;

% approimate gamma
x = fmin('f', 0, 3, sqrt(2));
idealVal = GPval.^x;
gamma = 1/x;

plot(GPval,idealVal,'k*-');
legend('Actual Data', 'Estimated Data', 2);
title('Camera Gamma Function of Gray Values');
xlabel('Input Values');
ylabel('Gamma Corrected Gray Values');
hold off;

