%% Main script for image analysis
% This live-script is the starting point of image analysis. Alternatively, a ".m"-file with the same name can be used
% 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 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