%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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
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
% Y = a1*exp(-((x-b1)/c1)^2)
% Fit the data set, constructing a cfit object containing the output
cFitTHLAdj = 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
FitTHLAdjParams=coeffvalues(cFitTHLAdj);
THLNoiseCentroid=FitTHLAdjParams(2)
THLNoiseCentroidWidth=FitTHLAdjParams(3)
THLNoiseCentroidSigma=THLNoiseCentroidWidth/sqrt(2)
% Round this to get the centre for use in later analysis
THLCentre=round(THLNoiseCentroid)
% Plot the data and the fit
figure;
plot(cFitTHLAdj,THLAdjBins,THLAdjEqualised)
% Adjust plot a bit
title('THL noise centroids after equalisation - 3D, day2','FontName','Arial','FontSize',16)
xlabel('THL value of centroid','FontName','Arial','FontSize',14)
ylabel('No of pixels','FontName','Arial','FontSize',14)
fig_handle=gcf;
set(fig_handle,'Units','centimeters')
set(fig_handle,'Position',[3,3,20,15])
% Position vector is [left, bottom, width, height]
axes_handle=gca;
set(axes_handle,'FontSize',12)
set(axes_handle,'FontName','Arial')
axis([380 460 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]
%}