function [segimg, top_slice, bot_slice, sumslices, zflip ] = preprocess( img, varargin )
%PREPROCESS segments and preprocesses a given image.
%
% Backwards compatibility mode
% Inputs
% ------
% img = image array to be processed
% imgFile = mat file containing image to be processed
% psfPath = path to psf to use
% downsample = [1xn] or single number defining the amount to downsample the image before preprocessing.
% adherent = boolean flag defining if the cell is adherent(forces largest slice to bottom)
% top_thresh = the fraction of fluorescence you wish to consider as the top slice
% bot_thresh = the fraction of fluorescence you wish to consider as the bottom slice
% display = boolean flag of whether to display progress on screen
% cellmask = binary valued image array masking the inside of the cell
%
% Outputs
% -------
% segimg = segmented image
% top_slice = integer value of the top slice where signal lies
% bot_slice = integer value of the top slice where signal lies
% sumslices = the cumsum of each slice
% zflip = boolean flag indicating if preprocess thinks the cell is upsidedown

% Copyright (C) 2006-2014 Murphy Lab
% Carnegie Mellon University
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published
% by the Free Software Foundation; either version 2 of the License,
% or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA.
%
% For additional information visit http://murphylab.web.cmu.edu or
% send email to murphy@cmu.edu

%Created by Devin Sullivan 3/21/13 moved from being a subfunction of
%seg_cell_and_nucleus.m
%
%G. Johnson 8/22/13 - support for 2D masks, even if img is 3D
%G. Johnson 8/26/13 - Re-mask the output of region_seg incase the snake moves
%                     out of bounds
%G. Johnson 1/15/14 - Remove psf as per Ivan's request
%D. Sullivan 4/24/14 Removed if statment that would not work if bot_slice and bot_seg were
%equal

%icaoberg 26/1/2014
%backward compatibility mode means that it will run as before but new call
%will simplify calling this method directly
if nargin == 9
    disp( 'Running method in backward compatibility mode' )
    imgFile = varargin{1};
    psfPath = varargin{2};
    downsample = varargin{3};
    adherent = varargin{4}; 
    top_thresh = varargin{5}; 
    bot_thresh = varargin{6};
    display = varargin{7};
    cellmask = varargin{8};
    param = [];
    param = ml_initparam( param, struct( 'display', false ) );
    param = ml_initparam( param, struct( 'debug', false ) );  
elseif nargin == 10
    disp( 'Running method in backward compatibility mode' )
    imgFile = varargin{1};
    psfPath = varargin{2};
    downsample = varargin{3};
    adherent = varargin{4}; 
    top_thresh = varargin{5}; 
    bot_thresh = varargin{6};
    display = varargin{7};
    cellmask = varargin{8};
    param = varargin{9};
    param = ml_initparam( param, struct( 'display', false ) );
    param = ml_initparam( param, struct( 'debug', false ) );  
elseif nargin == 2
    if ~isfield( param, 'image_file' )
        imgFile = '';
    else
        imgFile = param.image_file;
    end
    
    if ~isfield( param, 'psf_path' );
        psfPath = '';
    else
        psfPath = varargin{2};
    end
    
    if ~isfield( param, 'downsample' )
        downsample = [ 5, 5, 1 ];
    else
        downsample = param.downsample;
    end
    
    if ~isfield( param, 'adherent' )
        adherent = true;
    else
        adherent = param.adherent;
    end
    
    if ~isfield( param, 'top_threshold' )
        top_thresh = 0.98;
    else
        top_thresh = param.top_threshold;
    end
    
    if ~isfield( param, 'bottom_threshold' )
        bot_thresh = 0.02;
    else
        bot_thresh = param.bottom_threshold;
    end
    
    if ~isfield( param, 'display' )
        display = false;
    else
        display = param.display;
    end
    
    if ~isfield( param, 'cellmask' )
        cellmask = [];
    else
        cellmask = varargin{8};
    end 
    
    param = [];
    param = ml_initparam( param, struct( 'display', false ) );
    param = ml_initparam( param, struct( 'debug', false ) );  
else
    warning( 'Wrong number of input arguments. Exiting method.' );
    segimg = []; top_slice = []; bot_slice = []; sumslices = [];  zflip = [];
    return
end

if ~exist('cellmask', 'var') || isempty( cellmask )
    cellmask = ones(size(img));
end

imsize = size(img);
masksize = size(cellmask);

%icaoberg
%if the psfpath is empty then do not attempt to load psf
% Load PSFs
if ~isempty( psfPath )
    infPSF = imfinfo( psfPath);
    
    psf = zeros( infPSF(1).Height, infPSF(1).Width, length(infPSF));
    
    for I=1:length(infPSF)
        psf(:,:,I)=imread(psfPath,I);
    end
    
    psf = psf.^2; % Approximate confocal PSF
else
    psf = [];
end

%D. Sullivan 3/21/13 initialize the z flip to 0
zflip = 0;

%icaoberg 15/1/2014
%downsample PSFs and images as specified by user
%if the psf is not specified then the psf should be empty
if exist( imgFile, 'file' )
    load(imgFile)
else
    if ~isempty( psf )
        disp( 'Downsampling PSF and image'); tic;
    else
        disp( 'Downsampling image'); tic;
    end
    
    disp( [ 'Downsampling using vector [x,y,z]: [' ...
        num2str(downsample) ']' ] );
    
    if ~isempty( psf )
        psf = ml_downsize(psf,downsample,'linear');
    end
    
    resizeimg = ml_downsize( img, downsample, 'linear'); toc
    
%     if length(size(cellmask)) == 3 && masksize(3) ~= imsize(3)
%         [~,ind] = max(squeeze(sum(sum(cellmask,1),2)));
%         cellmask = cellmask(:,:,ind);
%     end
    
    if length(size(cellmask)) == 3
        %if the mask isnt the same size as the currently downsampled image
        if ~all(size(cellmask) == size(resizeimg))
            resizecellmask = ml_downsize(cellmask,downsample) > 0;
        else
            resizecellmask = cellmask;
        end
        
    else
        %D. Sullivan 5/5/14 - added check to see if downsampling is
        %necessary
         %if the mask isnt the same size as the currently downsampled image
        if (size(cellmask,1) ~= size(resizeimg,1)) || (size(cellmask,2)~=size(resizeimg,2))
            resizecellmask = ml_downsize(cellmask, downsample(1:2)) > 0;
        else
            resizecellmask = cellmask;
        end
        resizecellmask = repmat(resizecellmask, [1,1, size(resizeimg,3)]);
    end
    
    clear dnaim3
    
    %icaoberg 15/1/2014
    %only attempt to deconvolve image if psf is not empty
    if ~isempty( psf )
        fprintf(1,'%s\n','Deconvolving image'); tic;
        [image,psf2] = deconvblind(resizeimg,psf); toc
    else
        image = resizeimg;
    end
end


fprintf(1,'%s\n','Segmenting image'); tic;
sumslices = cumsum(squeeze(sum(sum(image))));
sumslices = sumslices/sumslices(end);
bot_slice=find(sumslices> bot_thresh); %0.05, 0.02
bot_slice=bot_slice(1);
top_slice=find(sumslices< top_thresh); %0.95, 0.98
%D. Sullivan 5/5/14 - if there are no slices less than the top_thresh, set
%top_slice to maximum value in image
if isempty(top_slice)
    top_slice = size(image,3);
else
    top_slice=top_slice(end);
end
%D.Sullivan 6/6/13, these two lines literally do nothing
% bot_slice=max(bot_slice,bot_slice);
% top_slice=min(top_slice,top_slice);

%fprintf(1,'Total slices=%i, Bottom slice=%i, Top slice=%i\n',length(sumslices),bot_slice,top_slice);
disp( ['Total number of slices: ' num2str(length(sumslices))] );
disp( ['Bottom slice index: ' num2str(bot_slice)] );
disp( ['Top slice index: ' num2str(top_slice)] );

mask = image > ml_rcthreshold(image(:));
for i = 1:size(mask,3)
    mask(:,:,i) = bwfill(mask(:,:,i), 'holes');
end
mask = repmat(sum(mask,3), [1,1,size(mask,3)]) > 0;

% bounds = ones(size(mask));
% bounds(3:end-3,3:end-3,:) = 0;
% mask(bounds >0) = 0;

disp( 'Cropping 3D image using mask' ); tic
mask = tz_maskimg_3d( mask, resizecellmask ); toc

disp( 'Region Based Active Contour Segmentation' ); tic
display = true;
maximum_iterations = 3000;
quit_tolerance = 0.00001;

%icaoberg 26/1/2014
%pass in parameter structure
segimg = region_seg(image, mask, ...
    maximum_iterations, 0.7, display, quit_tolerance, param ); toc

%icaoberg 26/1/2014
%this rearranges some files for reports
%needs cleanup
directory = [ pwd filesep 'showCurveAndPhi' ];
directory2 = [ param.preprocessingFolder filesep 'segmentation' ];
if ~exist( directory2 ); mkdir( directory2 ); end

try
    if exist( directory )
        files = dir( [ directory filesep 'iteration*' ] );
        img = [];
        for j=1:1:length(files)
            file = [ directory filesep files(j).name ];
            temp = imread( file );
            img = [ img; temp ];
        end
        
        if isfield( param, 'output_filename' )
            [ path, filename, ext ] = fileparts( param.output_filename );
            output_filename = [ param.preprocessingFolder filesep ...
                'segmentation' filesep filename '.png' ];
        else
            warning( 'Output filename not set. Inferring filename from number of files.' );
            number_of_files = length( dir( [directory2 filesep 'img*.png'] ) ) + 1;
            output_filename = [ param.preprocessingFolder filesep ...
                'segmentation' filesep 'img' num2str(number_of_files) ];
        end
        imwrite( img, output_filename );
    end
    
    if exist( directory )
        rmdir( directory, 's' );
    end
catch
    disp( 'Folder is empty. Removing directory' );
    rmdir( directory, 's' );
end

clear directory 
clear directory2
clear output_filename

segimg = ml_findmainobj(segimg);
segimg = tz_maskimg_3d(segimg, resizecellmask);

inds = find(sum(sum(segimg,1),2));
bot_seg = inds(1);
top_seg = inds(end);

bot_slice=max([bot_slice,bot_seg]);
top_slice=min([top_slice,top_seg]);

%devins 24/4/2014
%removed if statment that would not work if bot_slice and bot_seg were
%equal
if adherent %false, <param>
    %D. Sullivan  3/21/13 check if image is right side up.
    %since we expect the cell to generally get larger at the bottom for
    %adherent cell lines, if the slope of the areas is positive we should
    %flip the order
    
    %get total cell area per slice
    areas = squeeze(sum(sum(segimg)));
    %eliminate ones that are not in the cell
    cellareas = areas(areas~=0);
    %P = coefficients P(1)*X+P(2)
    P = polyfit([1:length(cellareas)],cellareas',1);
    if P(1)>0
        %set zflip flag
        zflip = 1;
        segimg = flipdim(segimg,3);
    end
    [bot_slice,top_slice]=fixadherent(segimg,bot_slice,top_slice);
end
clear image