%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Canon G1
%
%    1 2 3 4 5 6 
%  1 G M G M G M ...
%  2 Y C Y C Y C
%  3 G M G M G M ...
%  4 Y C Y C Y C
%    .
%    .
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clear; %close all;

% 15, 18, 21, 24, 27; 14, 26
fnum = 24;
fname = ['./spectrometer_data/CRW_03',num2str(fnum),'.ppm']

width = 2088; height = 1550; widthCr = 2088; heightCr = 1550;
fp = fopen(fname,'rb');
im = fread(fp,4*width*height,'uint16');
im2 = reshape(im,  width*4, height)';
imCrop = im2(1:heightCr, 1:widthCr*4);

Line = reshape(imCrop', 4, widthCr*heightCr);
gLine2 = Line(1, :); mLine2 = Line(2, :); 
cLine2 = Line(3, :); yLine2 = Line(4, :);
g = reshape(gLine2, widthCr, heightCr)'; 
m = reshape(mLine2, widthCr, heightCr)';
c = reshape(cLine2, widthCr, heightCr)'; 
y = reshape(yLine2, widthCr, heightCr)';

G = g; M = m; C = c; Y = y;
gamma = 2.2;

	[rsize,csize] = size(g);
	rG = G(1:2:rsize,1:2:csize);
	rM = M(1:2:rsize,2:2:csize);
	rC = C(2:2:rsize,2:2:csize);
	rY = Y(2:2:rsize,1:2:csize);
	G = rG; M = rM; C = rC; Y = rY; 
	[rsize,csize] = size(G);

%------------ Color CORRECTION -----------------%
rgb2gmcy =     [
			 0.11, 0.86, 0.08
			 0.50, 0.29, 0.51
			 0.11, 0.92, 0.75
			 0.81, 0.98, 0.08 ];
color_corr = pinv(rgb2gmcy)

color_corr(3,:) = color_corr(3,:)*2;



GMCY_arr = [G(:), M(:), C(:), Y(:)]';
RGB_arr = color_corr*GMCY_arr;
max_rgb = max(RGB_arr(:));
min_rgb = min(RGB_arr(:));
RGB_arr = RGB_arr.*(RGB_arr>0);
RGBarr = RGB_arr/max_rgb/2; % Now, RGBarr is from exactly 0 to 1.
gcRGBarr = (RGBarr).^(1/gamma); % Gamma Correction

imRGB = zeros(rsize,csize,3);
imRGB(:,:,1) = reshape(gcRGBarr(1,:),rsize,csize);
imRGB(:,:,2) = reshape(gcRGBarr(2,:),rsize,csize);
imRGB(:,:,3) = reshape(gcRGBarr(3,:),rsize,csize);
figure; image(imRGB); axis equal;
%------------ The END of Color CORRECTION -----------------%

if 1;
	dcm = zeros(24,4);
	for ind = 1:24;
		title([num2str(ind)]); BW = roipoly;
		dcm(ind,1) = mean(G(BW)); dcm(ind,2) = mean(M(BW)); 
		dcm(ind,3) = mean(C(BW)); dcm(ind,4) = mean(Y(BW));
	end;
else
	load dcm_crw_0308;
end;
close;
save(['dcm_crw',num2str(fnum)], 'dcm');
ls