% Code for the publication:
%  C.S. Smith, S. Stallinga, K.A. Lidke, B. Rieger, D. Grünwald, 
%  Probability-based particle detection that enables threshold-free and robust in vivo single molecule tracking,
%  Molecular Biology of the Cell, 2015 accepted. mbc.E15-06-0448
%
% Copyright 2015, University of Massachusetts, USA
%
% Carlas Smith, October 2015
%
% This code requires the toolbox DIPimage to be installed
% (http://www.diplib.org/download)

close all
clearvars
addpath(genpath('helperfunctions'))

%% Load data
% !!The necessary example data can be downloaded at: www.umassmed.edu/grunwaldlab/
fileStructure(1).Geller = '../data/grid_082015-1.tif';
fileStructure(1).Dark = '../data/dark-1.tif';
fileStructure(1).Data = '../data/merged3-1.tif';

%% camera calibration
%out(2) the 1/slope of the fit in [photons per ADU]: conversion ADU =
%detected #photons * slope !

a = bfopen( [ fileStructure(1).Geller]);
geller = double(cell2mat(permute(a{1}(1:min(100,end),1),[3 2 1])));

a = bfopen( [ fileStructure(1).Dark] );
bg = double(cell2mat(permute(a{1}(1:min(100,end),1),[3 2 1])));

out = cal_readnoise(geller, bg);

%%
PSFSigma=1.39; %size of PSF
nFrames = 20; %maximum number of frames to be processed

a = bfopen( [ fileStructure(1).Data] );
data = double(cell2mat(permute(a{1}(1:min(end,nFrames),1),[3 2 1])));
data=(data-repmat(mean(bg,3),[1 1 size(data,3)]))*out(2);

dataCam1 = data(1:512,1:512,1:min(end,nFrames));

%% GLRT Based Detection
[coordsCam1,detParCam1,cutProcessCam1] = LLRMapv2(dataCam1,PSFSigma);

% NOTE:
% % If you do not want to calculate cluster features you can specify: minPixels =[]
% [coordsCam1,detParCam1,cutProcessCam1] = LLRMapv2(dataCam1,PSFSigma,[]);
% 
% % If you want to compute the GLRT for all pixels without rejecting large areas with background compReduction=[]
% [coordsCam1,detParCam1,cutProcessCam1] = LLRMapv2(dataCam1,PSFSigma,0,[]);


%% Pre-Filter Detection Clusters Cam1
paramsPreFilterFits = getDefaultParamsPreFilterFits;
% paramsPreFilterFits = 
%     circularityMax: 3
%     circularityMin: 0.5000
%             pH1Max: 1
%             pH1Min: 0
%       minPixelDist: 0
%     clusterSizeMin: 0
%     clusterSizeMax: 50
[ maskPreFiltCam1 ] =  preFilterFits(coordsCam1,detParCam1,paramsPreFilterFits);

% Plot pre-filtered spots
optionsLLR = dipSubLoc2DSetOptions;
optionsLLR.BoxCenters = [coordsCam1(find(maskPreFiltCam1),1) coordsCam1(find(maskPreFiltCam1),2) coordsCam1(find(maskPreFiltCam1),3)];
optionsLLR.BoxSize = 2*(3*PSFSigma+1).*[1 1];
C=jet(256); 
index = double(floor(255.*(detParCam1.pH1(:,1)))+1);
optionsLLR.BoxColor = C(index,:);
optionsLLR.plotBoxes = 1;
optionsLLR.im = dip_image(cutProcessCam1,'uint8')';

h = dipSubLoc2D(optionsLLR);  

%% MLE Fit Intensities Cam1
paramsFit = getDefaultParamsFit;
% paramsFit       
%     FitSigma: 0
%     Iterations: 10
%     MaxCudaFits: 100000
%     PSFSigma: 1.3900
%     BoxSize: 11.3400
coodsUnCut=round(coordsCam1+(1.5*(2*PSFSigma+1)-0.5).*[ones(size(coordsCam1,1),1) ones(size(coordsCam1,1),1) zeros(size(coordsCam1,1),1)]);
paramsFit.FitSigma=true;

% Fit detections
[ rawFitResultsCam1 ] = fitBoxCenters( single(squeeze(dataCam1)),coodsUnCut,paramsFit);

%%
paramsFilterFits = getDefaultParamsFilterFits;
% paramsFilterFits
%       MinCRLBSTD: 0
%        MinPValue: 0
%       MinPhotons: 50
%            MinBg: 0
%       MaxCRLBSTD: 0.5000
%        MaxPValue: 1
%       MaxPhotons: Inf
%            MaxBg: Inf
%     MinPixelDist: 0
paramsFilterFits.MinPixelDist=2;
[ maskFilt1 ] =  filterFits(rawFitResultsCam1,paramsFilterFits);

% Plot filtered spots
optionsLLR = dipSubLoc2DSetOptions;
optionsLLR.BoxCenters = [rawFitResultsCam1.Coord(maskPreFiltCam1&maskFilt1,1) rawFitResultsCam1.Coord(maskPreFiltCam1&maskFilt1,2) rawFitResultsCam1.Frame(maskPreFiltCam1&maskFilt1)];
optionsLLR.BoxSize = 2*(3*PSFSigma+1).*[1 1];
C=jet(256); 
index = double(floor(255.*(detParCam1.pH1(maskPreFiltCam1&maskFilt1,1)))+1);
optionsLLR.BoxColor = C(index,:);
optionsLLR.plotBoxes = 1;
optionsLLR.im = dip_image(dataCam1);
h = dipSubLoc2D(optionsLLR);  