% Gregory Ng (gregng at stanford)
% EE 362/Psych 221
% Stanford University
% March 2005

%%%%%%%%%%%%%%%%%%%%%%%%%
% Attempt to estimate how much noise there is at each
% digital intensity level

% This routine takes into account the gamma coefficient for the
% response of the camera.

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Configure
% Offsets of the square
% Size of square

do_luminance = 0; % else do_rgb
plot_snr_same_graph = 0; % Plot SNR on the same graph

offsets = [112,1026 ; 432,1021 ; 744,1019 ; 1057,1013 ; 1371,1011 ; 1679,1007 ];
offsetlen = size(offsets, 1);
folder = '\\thefootball\Documents\camera\tung_iso50';

sy = 256; sx = 256;
x1=0;
x2=0;
y1=0;
y2=0;


sprintf_pattern1 = 'iso50_tung2_%05d.JPG';
%sprintf_pattern1 = 'n1_iso50_tung2_%05d.PNG';
prefix = 'n1';
imcount = 10;
%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (cameratype && cameratype < 0) 
	cameratype = 6;
end
imbegin = 1;

gamma = 1;
if (cameratype==1) 
	% Nikon D70
	offsets = [112,1026 ; 432,1021 ; 744,1019 ; 1057,1013 ; 1371,1011 ; 1679,1007 ];
	offsetlen = size(offsets, 1);
	
	sy = 256; sx = 256;
	
	folder = 'G:\Labelled Files';
	

	x1=245;
	x2=2293;
	y1=287;
	y2=1825;
	
	sprintf_pattern1 = 'iso50_tung2_%05d.JPG';
	%sprintf_pattern1 = 'n1_iso50_tung2_%05d.PNG';
	prefix = 'n1';
elseif (cameratype==2) 
	% S1 IS
	offsetlen = size(offsets, 1);
	x1=254;
	x2=1935;
	y1=68;
	y2=1293;	
	sy = 256; sx = 256;
	
	folder = 'D:\MyDocs\My Pictures\colornoise\canon_s1is\gradient';
	
	sprintf_pattern1 = 'gradient_%05d.JPG';
	%sprintf_pattern1 = 'n1_iso50_tung2_%05d.PNG';
	prefix = 'n1';
	imbegin=1;
	imcount = 10;
	
	gamma = 0.6805;

elseif (cameratype==3) 
	% S30
	offsetlen = size(offsets, 1);
	x1=121;
	x2=2048;
	y1=1;
	y2=1536;	
	sy = 256; sx = 256;
	
	folder = 'D:\MyDocs\My Pictures\colornoise\canon_s30\gradient';
	
	sprintf_pattern1 = '%05d_RT8.png';
	%sprintf_pattern1 = 'n1_iso50_tung2_%05d.PNG';
	prefix = 'n1';
	imbegin=2;
	imcount = 14;
elseif (cameratype==4) 
	% S30 on the macbeth chart
	offsetlen = size(offsets, 1);
	x1=121;
	x2=2048;
	y1=1;
	y2=1536;	
	sy = 256; sx = 256;
	
	folder = 'D:\MyDocs\My Pictures\colornoise';
	
	sprintf_pattern1 = 's30_macbeth_%05d_RT8.tif';
	%sprintf_pattern1 = 'n1_iso50_tung2_%05d.PNG';
	prefix = 'n1';
	imbegin=1;
	imcount = 21;
elseif (cameratype==5)
	ztitle = 'S1 IS on the LCD (ISO 50)';
	% S30 on the macbeth chart
	offsetlen = size(offsets, 1);
	x1=288;
	x2=1783;
	y1=293;
	y2=1471;	
	sy = 512; sx = 512;
	
	folder = 'D:\MyDocs\My Pictures\colornoise\gradients\s1islcd';
	
	sprintf_pattern1 = 'gradiso050_%05d.jpg';
	%sprintf_pattern1 = 'n1_iso50_tung2_%05d.PNG';
	prefix = 'iso50';
	imbegin=1;
	imcount = 20;
elseif (cameratype==6)
	% S1 IS on the LCD
	% S30 on the macbeth chart
	ztitle = 'S1 IS on the LCD (ISO 100)';
	offsetlen = size(offsets, 1);
	x1=288;
	x2=1783;
	y1=293;
	y2=1471;	
	sy = 512; sx = 512;
	
	folder = 'D:\MyDocs\My Pictures\colornoise\gradients\s1islcd';
	
	sprintf_pattern1 = 'gradiso100_%05d.jpg';
	%sprintf_pattern1 = 'n1_iso50_tung2_%05d.PNG';
	prefix = 'iso100';
	imbegin=1;
	imcount = 20;
elseif (cameratype==7)
	% S1 IS on the LCD
	% S30 on the macbeth chart
	ztitle = 'S1 IS on the LCD (ISO 200)';
	offsetlen = size(offsets, 1);
	x1=288;
	x2=1783;
	y1=293;
	y2=1471;	
	sy = 512; sx = 512;
	
	folder = 'D:\MyDocs\My Pictures\colornoise\gradients\s1islcd';
	
	sprintf_pattern1 = 'gradiso200_%05d.jpg';
	%sprintf_pattern1 = 'n1_iso50_tung2_%05d.PNG';
	prefix = 'iso200';
	imbegin=1;
	imcount = 20;
	

elseif (cameratype==8)
	% S1 IS on the LCD
	% S30 on the macbeth chart
	ztitle = 'S1 IS on the LCD (ISO 400)';
	offsetlen = size(offsets, 1);
	x1=288;
	x2=1783;
	y1=293;
	y2=1471;	
	sy = 512; sx = 512;
	
	folder = 'D:\MyDocs\My Pictures\colornoise\gradients\s1islcd';
	
	sprintf_pattern1 = 'gradiso400_%05d.jpg';
	%sprintf_pattern1 = 'n1_iso50_tung2_%05d.PNG';
	prefix = 'iso400';
	imbegin=1;
	imcount = 20;
else
	return;
end




for idx = imbegin:imcount
	if (idx == imbegin) 
		imseries = sprintf( sprintf_pattern1 , idx);
	else 
		imseries = [imseries; sprintf(sprintf_pattern1 , idx); ];
	end
end
clear sprintf_pattern;

%Nikon D70 data.
if (cameratype==1) 
imseries = [
'640-1.JPG';
'640-2.JPG';
'640-3.JPG';
];
end


%% Find the sample mean from the set of images.
xmean = zeros(256, 1);
yvar  = zeros(256, 3);
ycount = zeros(256, 3); % Count of pixels with this intensity (in this channel)
if (do_luminance==1) 
	[imavg imvar] = average_images_ycbcr(folder, imseries, x1, x2, y1, y2);
else
	[imavg imvar] = average_images(folder, imseries, x1, x2, y1, y2, gamma);
end


% Bin the averages into discrete intensities
imavg = round(imavg) + 1; % In the range [1,256?]

for val = 1:256


	for chan = 1:3
	
		findmat = (imavg(:,:,chan) == val);
		count = sum(sum(findmat));
		
		
		
		% Mean of the variance
		if (count == 0) 
			varmean = 0;
		else
			[rows cols vals] = find(findmat .* imvar(:,:,chan));
			%varmean = sum(sum( imvar(:,:,chan) .* findmat )) / count;
			varmean = mean(vals);
		end
		
		%clear rows, cols, vals, varmx;	
		ycount(val, chan) = count;
		yvar(val, chan) = varmean;
		
	end
	%disp('.');
end
clear findmat;


figure; 

intensity_axis = 0:255;
%intensity_axis = ((intensity_axis./255) .^ (1/0.6805)) .* 255;



if (do_luminance) 
	subplot(211); 
	hold on;
	plot(intensity_axis, yvar(:,1), 'k-'); 
	plot(intensity_axis, yvar(:,2), 'r-'); 
	plot(intensity_axis, yvar(:,3), 'b-'); 
	title('Noise energy per pixel intensity.');
	axis tight;

	subplot(212); 
	plot(intensity_axis, ycount(:,1), 'k-'); 
	plot(intensity_axis, ycount(:,2), 'r-'); 
	plot(intensity_axis, ycount(:,3) , 'b-'); 
	title('Count of number of pixels at intensity');
	
	hold off;
	axis tight;
else
	subplot(211); 
	hold on;
	plot(intensity_axis, yvar(:,1), 'r-'); 
	plot(intensity_axis, yvar(:,2), 'g-'); 
	plot(intensity_axis, yvar(:,3), 'b-'); 
	title(sprintf('Noise energy per pixel intensity %s', ztitle));

	%title('Noise energy per pixel intensity.');
	axis tight;

	subplot(212); 
	hold on;
	plot(intensity_axis, ycount(:,1), 'r-'); 
	plot(intensity_axis, ycount(:,2), 'g-'); 
	plot(intensity_axis, ycount(:,3), 'b-'); 
	title(sprintf('Count of number of pixels at intensity %s', ztitle));
	%title('Count of number of pixels at intensity');
	axis tight;
	hold off;


	%%%%%%%%%%%%%5
	% Compute the SNR
	
	
	if (plot_snr_same_graph) 
		figure(25);
	else
		figure;
	end
	hold on;

	for chan = 1:3
		clear snraxis  snrval;
		snrcount =  0;
		for idx = 1:256
			if (yvar(idx,chan) > 0 && ycount(idx, chan) > 1000) 
				snrcount = snrcount + 1;
				snraxis(snrcount) = idx-1;
				snrval(snrcount) = 20 / log(10) * log((idx-1) / sqrt(yvar(idx, chan)));
			end
		end
		if (chan == 1)
			plottype = 'r-';
		elseif (chan == 2)
			plottype = 'g-';
		else 
			plottype = 'b-';
		end
		if (snrcount > 0)
			plot(snraxis, snrval, plottype);
		end
		
		% Save the SNR data to files.
		filename = sprintf('snr-%d-c%d-axis.txt', cameratype, chan);
		eval(sprintf('save %s snraxis -ASCII -TABS', filename) );
		filename = sprintf('snr-%d-c%d-data.txt', cameratype, chan);
		eval(sprintf('save %s snrval -ASCII -TABS', filename) );
	end
	
	
	
	if (plot_snr_same_graph)
		title('Signal to Noise Ratio (SNR)');
	else
		title(sprintf('SNR %s', ztitle));
	end
	xlabel('Intensity (ADU)');
	ylabel('SNR (dB)');
	axis ([0 255 0 50]);
	hold off;

end

cameratype = -1;

%figure; imshow(log(result), []);
%figure; imshow(log(result) / log(10), []);