function gamma = evalGamma

% Usage: gamma = evalGamma
%
% Calculates the Olympus camera's gamma value from our
% lab data using a least-squares fit to an exponential curve.
% This function also shows a plot of the results.
%
% KLY 98

load olympusColors;

graySeries = 19:24;
olympusGrays = olympusColors(1,:,graySeries);

% Convert to matrix form
olympusGrays = fliplr(reshape(olympusGrays,3,6));
olympusGrays = sum(olympusGrays);

% Surface reflectance values at 600nm
patchValue = [.0375 .0979 .2171 .4002 .6381 .9531];

normalizedPatch = scale(patchValue, 0, 1, .0375, .9531);
normalizedGrays = scale(olympusGrays, 0, 1, min(olympusGrays), ...
    max(olympusGrays));

figure
plot(normalizedPatch, normalizedGrays,'--');

% Fit data to exponential curve
% -----------------------------
% The data is of the form G = kP^x, 
% i.e. log(G) = log(k) + xlog(P), where x is the gamma
% coefficient. We then use polyfit to find the best line
% through those points. I am throwing out the value for
% the darkest gray because it seems a little off.

logP = log(normalizedPatch(3:6));
logG = log(normalizedGrays(3:6));
gamma = polyfit(logP, logG, 1);
gamma = gamma(1);

x=linspace(0,1,500);
y=x.^gamma;
hold on
plot(x,y)
grid on
xlabel('Normalized reflectance values')
ylabel('Normalized framebuffer values')
title('Gamma Curve for the Olympus Digital Camera')

text(.7, .78, ['gamma = ' num2str(gamma)])
legend('Original data', 'Least-squares fit',2)
hold off
