function [coords,dectectionPar,cutProcess] = LLRMapv2(process,PSFSigma,minPixels,compReduction,significance,iterations,split,maxCudaFits)
%LLRMapv2 Detects single mocules according to the GLRT framework on CUDA 6.5 enables device.
% SYNOPSIS:
%  [coords,dectectionPar,cutProcess] = LLRMapv2(process,PSFSigma,minPixels,
%                               compReduction,iterations,split,maxCudaFits)
% 
% PARAMETERS:
%     process: single molcule data (corrected for gain and dark counts)
% 
%     PSFSigma: sigma of diffraction limited PSF in pixels
% 
%     minPixels: minimal size of detection cluster. minPixels = [] obbits
%     calculation of detection features dectectionPar =[]
% 
%     compReduction: computational complexity reduction level in standard 
%     deviations from mean. Rejectregions with small change of single 
%     molecules, based on wavelet filtering. compReduction = [] results in 
%     computation of GLRT on all pixels.
% 
%     significance: false discovery rate 
% 
%     iterations: number of iterations for MLE algorithm.
% 
%     split: split clusters using watershed algorithm
% 
%     maxCudaFits: maximum numbers of cuda fits on single GPU device. If
%     n CUDA enables devices are ready 1xn vector is expected.
% 
% 
% DEFAULTS:
%     minPixels = 0; 
%     compReduction = 2;
%     significance = 0.05;
%     split = false;
%     maxCudaFits = 1000000;
% 
% 
% OUTPUTS:
%   coords: matrix with detection coordinates.
%   dectectionPar
%         .circularity: circularity of cluster using P2A feature.
%                 .pH1: detection probability of cluster.
%         .clusterSize: cluster size
%                  .ll: labeled clusters
%                  .hh: binary image containing raw detections
%             .pfa_adj: false discovery rate corrected false positive 
%                       probailty
%             
%   cutProcess: cropped image with the actuall test pixels.

% COMMENT
% This code requires the NVIDA CUDA 6.5 Driver and Toolkit to be installed. 
%  
% LITERATURE:
% Probability-based particle detection that enables threshold-free and 
% robust in vivo single molecule tracking, by Carlas .S. Smith, 
% Sjoerd Stallinga, Keith .A. Lidke, Bernd Rieger, David Grünwald
    
    if nargin < 3
        minPixels = 0;
    end
    if nargin < 4
        compReduction = 2;
    end
	if nargin < 5
        significance = 0.05;
    end
    if nargin < 6
        iterations = 8;
    end
    if nargin < 7
        split = false;
    end
    if nargin < 8
        maxCudaFits = 1000000;
    end

    dectectionPar=[];
    szx = 1.5*(2*PSFSigma+1)-0.5:1:size(process,1)-1.5*(2*PSFSigma+1)-0.5;
    szy = 1.5*(2*PSFSigma+1)-0.5:1:size(process,2)-1.5*(2*PSFSigma+1)-0.5;
    szx=szx(1:floor(size(szx,2)/4)*4);
    szy=szy(1:floor(size(szy,2)/4)*4);
    
    [Ym,Xm,Zm] = meshgrid(szx,szy,0:size(process,3)-1);
    XX = reshape(Xm,[ numel(Xm) 1]);
    YY = reshape(Ym,[ numel(Ym) 1]); 
    ZZ = reshape(Zm,[ numel(Zm) 1]);

    H2 = 1/16;
    H1 = 1/4;
    H0 = 3/8;
    g{1} = [H2,H1,H0,H1,H2];
    g{2} = [H2,0,H1,0,H0,0,H1,0,H2];

    n_smooth=2;

    cutProcess = double(cut(process,[size(Ym,1) size(Ym,2) size(Ym,3)]));
    img =max(cutProcess-double(smooth(cutProcess,[1 1 0]*sqrt(n_smooth)*9)),0);
    
    locIm = ones(size(cutProcess));
    if ~isempty(compReduction)
        kernel(1).filter=g{1};
        kernel(1).origin = 3;
        kernel(2).filter=g{1};
        kernel(2).origin = 3;
        kernel(3).filter=1;
        kernel(2).origin =1;
        V{1} = dip_separableconvolution(cutProcess,kernel,[1 1 0]);
        kernel(1).filter=g{2};
        kernel(1).origin = 5;
        kernel(2).filter=g{2};
        kernel(2).origin = 5;
        kernel(3).filter=1;
        kernel(2).origin =1;
        V{2} =dip_separableconvolution(V{1},kernel,[1 1 0]);
        W{2} = double(V{1}-V{2});
        locIm =  W{2}	>mean(dip_image(W{2}),[],[1 2])+compReduction*std(dip_image(W{2}),[],[1 2]);
    end
    idxIm = find(locIm);

    [ ROIStack ] = cMakeSubregions(XX(idxIm),YY(idxIm),ZZ(idxIm),3*(2*PSFSigma+1),single(process));
    [~, ~, LLr] = gpuGaussMLEv3(permute(single(ROIStack),[2 1 3 4]),single(PSFSigma),iterations,0,true,maxCudaFits); 

    pfa =ones(size(Ym));
    pfa(idxIm)  =LLr(3,:);

    [ h,~,pfa_adj ]=fdr_bh(reshape(pfa,[numel(pfa) 1]),significance,'dep','no');

    hh=reshape(h,size(pfa));
    pfa_adj=reshape(pfa_adj,size(pfa));
    coords=[];
    circularity=[];
    clusterSize=[];
    pH1=[];
    ll=newim(size(hh),'sint32');

    if ~isempty(minPixels)
        for i=1:size(hh,3)
            if split
                h=squeeze(bclosing(hh(:,:,i)));
                D = -dt(h);
                D(~h) = +Inf;
                ll(:,:,i-1) = label(~dipwatershed(double(D),1).*squeeze(hh(:,:,i)),1);
            else
                ll(:,:,i-1) = label(squeeze(hh(:,:,i)),1);

                msr = measure(squeeze(ll(:,:,i-1)), squeeze(cutProcess(:,:,i)), ({'Gravity','P2A','Size'}));
                coords =cat(1,coords, [msr(msr.Size> minPixels).Gravity'  (i-1).*ones(size( msr(msr.Size> minPixels).Gravity,2),1)]);
                circularity = cat(1,circularity,msr(msr.Size> minPixels).P2A');
                clusterSize = cat(1,clusterSize,msr(msr.Size> minPixels).Size');
                pH1 = cat(1,pH1,(min(max(2.*msr(msr.Size> minPixels).Size'./((2*(PSFSigma))^2*pi),0),1)));  
            end
        end
        dectectionPar.circularity = circularity;
        dectectionPar.pH1 = pH1;
        dectectionPar.clusterSize=clusterSize;
        dectectionPar.ll=ll;
        dectectionPar.hh=hh;
        dectectionPar.pfa_adj=pfa_adj;    
    else
            im_max = (cutProcess>=.999*maxf(cutProcess,[(2*PSFSigma+1) (2*PSFSigma+1) 0],'rectangular'))&hh;
            coords=findcoord(im_max);
    end
end