%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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;

if 0;
	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);

	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; % 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;
end;
	%------------ The END of Color CORRECTION -----------------%

	if 0;
		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
		% 15, 18, 21, 24, 27; 14, 26
		load dcm_crw15;dcm1 = dcm
		load dcm_crw18;dcm2 = dcm
		load dcm_crw21;dcm3 = dcm
		load dcm_crw24;dcm4 = dcm
		load dcm_crw27;dcm5 = dcm
		dcm = dcm1 + dcm2 +dcm3 +dcm4 + dcm5;
	end;


if 1;
	load spec1.mat; S1 = spect';
	load spec2.mat; S2 = spect';
	load spec3.mat; S3 = spect';
	load spec4.mat; S4 = spect';
	So = S1;% + S2 + S3 + S4;
	%
	%
	% dcm: Digital Camera Measurements of 24 color patches. 24x4 (G,M,C,Y)
	% S: Reflected spectrum of 24 color patches measured by spectrometer. 24x101
	% Resp: Responsivity of CMYG sensors. 101x4
	%
	% dcm = S*Resp;
	%

	S = So(:,1:5:101);
	size(S)
	cond(S)

	Resp =	pinv(S)*dcm;
	%figure; plot(wavelength,Resp); 
	figure; plot(Resp/max(Resp(:))); 
	title('G,M,C,Y responsivity function');

	if 0;
		ns = null(S);
		dd = Resp(:,1)*ones(1,77) + 1e6*ns;

		figure;
		subplot(411); plot(dd(:,1:20));
		title('possible Green sensor solutions');
		subplot(412); plot(dd(:,21:40));
		subplot(413); plot(dd(:,41:60));
		subplot(414); plot(dd(:,61:77));
	end;
end;


if 1;
load var_jong;
figure;
R_u % my spectral responsivities
S_u = S_u + 0*1e-6*rand(24,34); % my spectra
for ind = 2:4;
	y = S_u*R_u + 10^(ind)*0.0007*randn(24,4); 
	% my responses (in electrons) with noise


		%estimate the responsivities
		S_uu = S_u(:,1:2:34);
		size(S_uu)
		cond(S_uu)

		R_hat = pinv(S_uu)*y; 
		%R_hat = pinv(S_uu)*dcm; 



	% Thoughts:
	y_2 = S_uu*R_hat;
	%Even if the R_hat are wrong they can give the right response? ->
	% Oh yeah, even completely screwed up sensors can give the same responses
	% as the right sensors (nullspace!)
	% Is the reason that the nullspace is so big?
	% How do we get a smaller nullspace? The dim of the nullspace is: 34(wavelengths)-24(surfaces) 
	%  -> More surfaces
	%  -> less wavelength
	% Best is: Measure more (independent) surfaces then wavelength, then we'll have a least
	% square problem as opposed to a least norm!

	% A dirty solution would be to reduce the number of wavelength by applying
	% a blur filter in the wavelength domain, and then only to solve for say 5-15 wavelengths

	if 1;
	subplot(221); 
	plot(R_u(:,1)/max(R_u(:)),'g'); hold on;
	plot(R_u(:,2)/max(R_u(:)),'c')
	plot(R_u(:,3)/max(R_u(:)),'m')
	plot(R_u(:,4)/max(R_u(:)),'y')
	axis tight;
	title('6*10^4');
	end;

	subplot(2,2,ind); 
	plot(R_hat(:,1)/max(R_u(:)),'g'); hold on;
	plot(R_hat(:,2)/max(R_u(:)),'c')
	plot(R_hat(:,3)/max(R_u(:)),'m')
	plot(R_hat(:,4)/max(R_u(:)),'y')
	axis tight;
end;
%figure; plot(y_2)
end;