%% Main script for image analysis
% Choose folders to be analysed.
% 
% Execute this section only if you have one folder selected to be analysed 
% (not batch):

%global PARAM, PARAM.folder = uigetdir('S:\Jeuken Microscope Data\Ievgen Mazurenko\', 'Please select folders to be analysed');
%% 
% Defining global parameters to be used within the analysis 
%%
tic;
global PARAM;
PARAM.fitM = 6;        % half-height of the window for vesicle fitting (by rows)  
PARAM.fitN = 6;        % half-width of the window for vesicle fitting (by colomns) 
PARAM.pointErrCri = 15;% fit: maximal allowed error per point 
PARAM.r2Cri=0.8;    % fit: minimal allowed R^2 of the fit
PARAM.circularCri=0.8; % fit: minimal allowed circularity of the spot
PARAM.widthCriMin=0.5; % fit: minimal allowed width of the spot
PARAM.widthCriMax=2.5; % fit: maximal allowed width of the spot   
PARAM.amplCri=20;      % fit: minimal allowed amplitude of peak above baseline
%PARAM.baseCri=PARAM.threshold * 1.2;     % fit: maximal allowed baseline level of the peak 
PARAM.edgeCri=6;       % fit: minimal allowed distance to the image edges
PARAM.centreDisCri=2;  % fit: maximal allowed distance between observed and fitted maximum
PARAM.distanceCri=10;  % fit: minimal allowed distance between vesicles
PARAM.percentOfFramesCri=50; % frames fit: maximal allowed percentages of frames with fit unsatisfied by criterias
PARAM.atto = [0, 0.05, 1];% calibration curve to convert ATTO intensity to diameter. d = p1 + p2 * Int * p3
% p3 - correction factor for different concentrations of ATTO
PARAM.pHCalib = [7.68, 0.115, 1.729, 1]; % calibration data (pKa, Ra, Rb, b) to convert intensity ratio to pH 

fprintf('Analysing the folder %s \n', PARAM.folder);
fprintf('Analysis started at  %s \n', datetime('now'));
%% 
% Loading tif-images from the folder and averaging them
%%
[frames, meanImgCh1, meanImgCh2] = fTifLoad(PARAM.folder);
%% 
% Threshold calculation. Plotting two average images and binary mask image 
% for 410 nm to verify the threshold
%%
PARAM.threshold = fCalculateThreshold(meanImgCh2, 0.1);  % threshold to find vesicle maximum on channel 2 (HPTS_L)
PARAM.baseCri=PARAM.threshold * 1.2;                     % fit: maximal allowed baseline level of the peak
figure;
subplot(1,2,1), imshow(imadjust(meanImgCh1)), title('470 nm (HPTS_H)', 'FontSize', 14);
subplot(1,2,2), imshow(imadjust(meanImgCh2)), title('410 nm (HPTS_L)', 'FontSize', 14);
figure, imshow(meanImgCh2 > PARAM.threshold), title('Binary mask 410 nm (HPTS_L)', 'FontSize', 14);
%% 
% Finding initial maximums on th 410nm channel using <https://uk.mathworks.com/matlabcentral/fileexchange/37388-fast-2d-peak-finder 
% FastPeakFind> function. fspesial should correspond to expected peak size and 
% shape. The coordinates of the peaks are saved in the "vesicles" structure.
%%
[coorPeaks, binImage]=FastPeakFind(meanImgCh2, PARAM.threshold, fspecial('gaussian', 13,2));
vesicles = fCreateVesicles(coorPeaks);
%% 
% Fitting each vesicle to 2D-Gaussian signal
%%
[fittedImg, vesicles] = fGauss2DFit(meanImgCh2, vesicles);
%% 
% Filtering the fitted vesicles on the basis of defined in the parameters 
% criteria
%%
filteredVesicles = fVesicleFilter(vesicles, meanImgCh2);
fSeeRandomVesicles(meanImgCh2, fittedImg, filteredVesicles);
%% 
% Loading time stamps and fitting the vesicles on all time frames
%%
timeStamps = load(strcat(PARAM.folder, '\', 'DeltaTime.csv')); % replace with name of timestamps if different
[fittedFrames, vesicles3D] = fGauss2DFitFrames(frames, filteredVesicles, timeStamps);
%% 
% Filtering the vesicles on all timeframes. Saving the data about each vesicle 
% on the disk
%%
filteredVesicles3D = fVesicleFilterFrames(vesicles3D);
fSeeRandomTraces(filteredVesicles3D);
save(strcat(PARAM.folder, '\', 'filteredVesicles3D.mat'), 'filteredVesicles3D', '-v7');
%% 
% Additional functions to see and save results:
%%
%fWriteGif(frames, fittedFrames, filteredVesicles3D, [1,2]); % if you need to create gif-animation for a defined vesicle
fSaveVesicles(filteredVesicles3D); % saves vesicles ratio vs time traces in the ascii file for further use in external programs
medianTrace = fGetMedian(filteredVesicles3D); % calculates and save median of all vesicles pH-change
figure, plot(medianTrace(1,:), medianTrace(2,:)), title('Median trace of all vesicles', 'FontSize', 14);
%% 
% Loading Atto image, fitting and converting to the vesicles diameter
%%
filteredVesicles3D = fGetATTO(filteredVesicles3D);
%% 
% Converting vesicles structure for traces structure optimized for the work 
% with pH-traces
%%
traces = fConvertVesicles(filteredVesicles3D);
%% 
% Filtering vesicles traces for out-of-range points and outliers
%%
traces = fFilterTraces(traces);
%% 
% Converting ratios to pH
%%
traces = fCalculatepH(traces);
fSaveTraces(traces(~[traces(:).delFlag])); % saves filtered pH into ascii file
%%
toc;
fprintf('Analysis ended at  %s', datetime('now'));
%%
%end