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

% TODO - go through the code, study and try to make it more efficient

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

if nargin < 6
   stack_size = size(array_in,3);
  % stack_size = 64;
end

cell_array = cell_partitioning(array_in, stack_size, patchradius);
anscombe = @(img) 2 * sqrt(img + (3/8)); % Anscombe transform to convert Poisson noise into Gaussian noise
parfor (m = 1:length(cell_array), 6)
    cell_array{m} = pad_array(cell_array{m}, patchradius); % pad array
    % Subimg to denoise
    fimg_med = extra.candle.median3D(cell_array{m}, 1);
    
    % Background detection
    if fast_background
        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
    
    cell_array{m} = anscombe(cell_array{m});
    fimg_med = anscombe(fimg_med);
    
    [MAP, ~] = extra.candle.GaussianSTD_MAP(cell_array{m}, 2*searchradius);
    % Denoising
    cell_array{m} = extra.candle.adapt_onlm_collaborative(single(cell_array{m}), ...
                                                   searchradius, ...
                                                   patchradius, ...
                                                   single(MAP), ...
                                                   beta, ...
                                                   single(fimg_med), ...
                                                   single(mask));
    fimg_med = [];
    % Optimal Inverse Anscombe Transform
    cell_array{m} = extra.candle.OVST(cell_array{m});
%     fimgMed = extra.candle.OVST(fimgMed);
    
    % Conversion
    cell_array{m} = trim_array(cell_array{m}, patchradius); % trim
    
    disp(['Done with batch ', num2str(m)])
end

array_out = single([]);
for k = 1:length(cell_array)
    if k == 1 % make sure the part of the output array is correctly indexed
        start = 1;
    else
        start = size(array_out, 3) + 1;
    end
    num_pics = size(cell_array{k}, 3);
    array_out(:, :, start:(start+num_pics-1)) = cell_array{k};
    cell_array{k} = [];
end

end

