function plotGamma(im1,im2,im3,im4,im5,im6)
% Function that takes 6 gray images scanned from the Macbeth
% color chart and plots the intensity versus the normalized RGB
% value

% the intensity values for the test patches from white to black:
y1 = 1.00;
y2 = 0.6255;
y3 = 0.4160;
y4 = 0.2295;
y5 = 0.0958;
y6 = 0.0519;

% Get the median RGB values for the test images
[r1,g1,b1] = getRGB(im1);
[r2,g2,b2] = getRGB(im2);
[r3,g3,b3] = getRGB(im3);
[r4,g4,b4] = getRGB(im4);
[r5,g5,b5] = getRGB(im5);
[r6,g6,b6] = getRGB(im6);

% Normalize the RGB values from [1,256] to [0,1]
i1 = (mean(uint8([r1,g1,b1]))-1)/255;
i2 = (mean(uint8([r2,g2,b2]))-1)/255;
i3 = (mean(uint8([r3,g3,b3]))-1)/255;
i4 = (mean(uint8([r4,g4,b4]))-1)/255;
i5 = (mean(uint8([r5,g5,b5]))-1)/255;
i6 = (mean(uint8([r6,g6,b6]))-1)/255;

% Will have the relative intensities in the form of Y values, one
% for each patch.
% With the 6 data patches we solve for gamma as follows:
% a) assume that rgb = (I/k)^(1/gamma)
% b) rewrite as: log(rgb) = (1/gamma)*logI - (1/gamma)*log (k)
% c) find the "best fit line" (poly order 1)
I_frame = [i6 i5 i4 i3 i2 i1];
I_intensity = [y6 y5 y4 y3 y2 y1];
p = polyfit(log(I_intensity(3:6)), log(I_frame(3:6)), 1)
gamma = 1/p(1)
k = exp(-p(2)*gamma)

% plot the gamma 
t = [0:0.1:1];
plot(t,(t/k).^(1/gamma));


% see that gamma = 1.8723







