%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%       WORKING WITH THRESHOLD DISTRIBUTIONS
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% In Pixelman, after a threshold equalisation we can save the histogram of
% the threshold distributions with the THL adjust bits set low, high, and
% after equalisation. The results after equalisation are particularly
% useful for finding the noise centre.
%
% The data file consists of ASCII data, plus a header. The header includes
% the THL start and end points, and also the step size, so we need to grab
% these in order to understand the data.
%
% Then, we plot the THL centroids after equalisation, and make a Gaussian
% fit
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The only argument needed is 'ThreshEqFilename', the name of the file
% containing the histograms

% Open the data file. This function returns a "file ID" we use as an
% argument later

ThreshEqFilename = '3D1_THLAdjust_redo_redo_noisecentroid_histogram.txt';
ThreshEqFID = fopen(ThreshEqFilename); 

% Then, we use textscan to parse the file. If we use it repeatedly, then
% each time we restart where the previous one left off.
% Note that this function gives cell arrays - need to convert to suitable
% format, using DataRead as an intermediate variable

% Skip 2 header lines, and then read in 1 line and grab the THL start value
% The %*s tells Matlab to ignore the string at the start of the line
DataRead = textscan(ThreshEqFID, '%*s %d', 1, 'HeaderLines', 2);
THLStart=double(DataRead{1})

% Then, grab the end points and spacing
DataRead = textscan(ThreshEqFID, '%*s %d', 1);
THLEnd=double(DataRead{1})
DataRead = textscan(ThreshEqFID, '%*s %d', 1);
THLSpace=double(DataRead{1})

% Check no. of bins
THLNumSteps=(THLEnd-THLStart)/THLSpace

% Note that the results are a histogram - the first data entry will
% correspond to the bin between THLStart and THLStart+THLSpace. To deal
% with this, I will make the corresponding "X value" the middle of each
% bin. This is done as a column vector.

THLAdjBins=((THLStart+THLSpace/2):THLSpace:(THLEnd-THLSpace/2))';

% Skip 3 more lines to reach the data, and then just read in the 3 columns.
% We don't specify a number of lines, so we'll keep going to the end of the
% file
% The resulting data is stored as a cell array - we end up with 3 cell
% elements, each of which is a column vector.
DataRead = textscan(ThreshEqFID, '%d %d %d', 'HeaderLines', 3);

% Grab the appropriate histogram data, using double format
THLAdjLow=double(DataRead{1});
THLAdjHigh=double(DataRead{2});
THLAdjEqualised=double(DataRead{3});

% Plot the equalised distributions, if we want
%figure;
%plot(THLAdjBins,THLAdjEqualised)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FITTING A GAUSSIAN TO THE DATA

% We will fit all 3 of the distributions

% Y = a1*exp(-((x-b1)/c1)^2)

% Fit the data set, constructing a cfit object containing the output

[cFitTHLAdjEq,FitTHLAdjEqErrors] = fit(THLAdjBins,THLAdjEqualised,'gauss1');

% Code to return the names and values of paramters of fit
%coeffnames(cFitTHLAdj)
%coeffvalues(cFitTHLAdj)

% We specifically want to return the centre pos

FitTHLAdjEqParams=coeffvalues(cFitTHLAdjEq);
THLNoiseCentroidEq=FitTHLAdjEqParams(2)
THLNoiseCentroidWidthEq=FitTHLAdjEqParams(3);
THLNoiseCentroidSigmaEq=THLNoiseCentroidWidthEq/sqrt(2)



cFitTHLAdjLow = fit(THLAdjBins,THLAdjLow,'gauss1');
FitTHLAdjLowParams=coeffvalues(cFitTHLAdjLow);
THLNoiseCentroidLow=FitTHLAdjLowParams(2)
THLNoiseCentroidWidthLow=FitTHLAdjLowParams(3);
THLNoiseCentroidSigmaLow=THLNoiseCentroidWidthLow/sqrt(2)



cFitTHLAdjHigh = fit(THLAdjBins,THLAdjHigh,'gauss1');
FitTHLAdjHighParams=coeffvalues(cFitTHLAdjHigh);
THLNoiseCentroidHigh=FitTHLAdjHighParams(2)
THLNoiseCentroidWidthHigh=FitTHLAdjHighParams(3);
THLNoiseCentroidSigmaHigh=THLNoiseCentroidWidthHigh/sqrt(2)



% Plot the data and the fit

figure;
hold on;
plot(cFitTHLAdjEq,THLAdjBins,THLAdjEqualised)
plot(cFitTHLAdjLow,THLAdjBins,THLAdjLow)
plot(cFitTHLAdjHigh,THLAdjBins,THLAdjHigh)


%bar(ClusterCounts,ClusterHist,'c');
%linehandle=plot(cClusterPeakFit,'r');
%set(linehandle,'LineWidth',1);
hold off;



% Adjust plot a bit
xlabel('THL value of centroid','FontName','Arial','FontSize',12)
ylabel('No of pixels','FontName','Arial','FontSize',12)
fig_handle=gcf;
set(fig_handle,'Units','centimeters')
set(fig_handle,'Position',[3,3,12,10])
% Position vector is [left, bottom, width, height]
axes_handle=gca;
set(axes_handle,'FontSize',10)
set(axes_handle,'FontName','Arial')
axis([340 500 0 1])
axis 'auto y';



%{ 
% Extra code for convenient plotting of THL maps
title('THL adjust bits - Planar, day 2 minus day 1','FontName','Arial','FontSize',16)
fig_handle=gcf;
set(fig_handle,'Units','centimeters')
set(fig_handle,'Position',[3,3,18,15])
% Position vector is [left, bottom, width, height]
%}