function array_out = filter_candle_old(array_in, ...
                                         patchradius, ...
                                         searchradius, ...
                                         beta, ...
                                         fast_background)
% TODO - description
% TODO - validate inputs

% TODO - delete

if fast_background
    fprintf('Fast processing of the background\n');
else
    fprintf('Accurate processing of the whole image\n');
end

% pad
array_in = pad_array(array_in, patchradius);

% filtration
fact = 64 - searchradius;
dim = size(array_in, 3);
substack = ceil(dim/fact);
if substack == 2
    fact = fact*2;
    substack = 1;
end
if substack == 1
    fact = dim;
end
for ind = 1:substack
    % To avoid out of memory the stack is processed by substacks of 64 images
    s_init = 1 + (ind-1)*fact; % beguinning of the substack
    s_end = ind*fact; % end of the substack
    % Limits of the substack
    if s_end > dim
        s_end = dim;
    end
    if (s_end-s_init) < 12
        s_init = s_end-12;
    end
    if s_init<1
        s_init = 1;
    end
    % Padding between substacks
    if (s_init-searchradius) <= 1
        pad_init = 0;
    else
        pad_init = searchradius;
    end
    if (s_end+searchradius) > dim
        pad_end = 0;
    else
        pad_end = searchradius;
    end
    % Subimg to denoise
    img = single(array_in(:,:,s_init-pad_init:s_end+pad_end));
    fimg_med = extra.candle.median3D(img,1);
    % Background detection
    if fast_background == 1
        mini = min(fimg_med(:));
        maxi = max(fimg_med(:));
%         average_val = mean(fimgMed(:));
        delta = (maxi-mini)/10000;
        k=3;
        N=k^3;
        ConvMed=(convn(fimg_med,ones(k,k,k),'same')/N);
        bins = double(mini):delta:double(maxi);
%         [nb,histb] = histc(ConvMed(:),bins);
        [nb,~] = histcounts(ConvMed(:),bins);
        [~, loc] = max(nb);
        Threshold = bins(loc);
        size_mask = size(fimg_med);
        mask = ones(size_mask);
        mask( ConvMed<(Threshold+5*delta) ) = 0;
    else
        size_mask = size(fimg_med);
        mask = ones(size_mask);
    end
    % Anscombe transform to convert Poisson noise into Gaussian noise
    anscombe = @(img) 2 * sqrt(img + (3/8));
    
    img = anscombe(img);
    fimg_med = anscombe(fimg_med);
    
    [MAP, ~] = extra.candle.GaussianSTD_MAP(img, 2*searchradius);
    % Denoising
    fimg = extra.candle.adapt_onlm_collaborative(single(img), ...
                                                   searchradius, ...
                                                   patchradius, ...
                                                   single(MAP), ...
                                                   beta, ...
                                                   single(fimg_med), ...
                                                   single(mask));
    % Optimal Inverse Anscombe Transform
    fimg = extra.candle.OVST(fimg);
%     img =  extra.candle.OVST(img);
%     fimgMed = extra.candle.OVST(fimgMed);
    
    % Conversion
    array_out(:,:,s_init:s_end) = fimg(:,:,pad_init+1:end-pad_end);
    
    disp(['Dataset ', num2str(ind), ' of ' num2str(substack), ' processed.'])
end
array_out = trim_array(array_out, patchradius); % trim
end

