function [] = cellsegmentation(varargin)
% CELLSEGMENTATION Cell segmentation of data stored on disk
%
%   CELLSEGMENTATION(FOLDER,STS,STE,PLS,PLE,MINVOL,MAXVOL,PRMFILE) 
%
%   Segmentation of cells in the folders specified in FOLDERS which can 
%   either be a single string or a cell array of folders. 
%   Reading all stacks within the numbering in VALSTACKS, 
%   assuming the name STS to STE. The planes are
%   specified as start and stop in PLS and PLE and the segmentation 
%   method in METHOD. If PLE = Inf all planes are used.
%
%   MINVOLFULL and MAXVOLFULL are the limits for the diameter of a cell, 
%   and SEGMCH is the channel to use for segmentation in the multidimensional 
%   (4D) image array. PRMFILE is the possibility to give a parameter file
%   for segmentation. This is very useful for reproducibility of settings
%   and results.
%
%   CELLSEGMENTATION(...,'prm',PRM)
%   adds a parameter struct with command-line defined parameters for the
%   segmentation process.
%
%   The segmentation method must be chosen by setting PRM.SEGMMETHOD.
%   Choices for PRM.SEGMMETHOD:
% 
%   PRM.SEGMMETHOD = 1 (default) : Calling SEGMSURFWAT for segmentation of surface stained 
%   cells. The marker method should be selected via MPRM.MINIMAMETOD.
%
%   PRM.SEGMMETHOD = 2: Calling SEGMCT for cytoplasmically stained cells. Can also
%   be used for segmentation of nuclei.
%  
%   The method for selecting markers should be specified.
%
%   MPRM.GETMINIMAMETHOD = 'MINIMACELL': Attempts to load manually defined markers
%   which must be of the format STACK1-MASK.MAT, STACK2-MASK.MAT and in the
%   same folder as the data. The manual markers must be binary images with
%   ones at the markers and zero elsewhere, and the image must have same
%   dimension as the segmentation image. There must be two variables, MINIMA 
%   with all markers and MINIMACELL with only cell markers. Default value 
%   MRPM.METHOD = 'PLANE'.
%
%   MPRM.GETMINIMAMETHOD = 'NUCLEUS': Using nucleus markers. This option
%   requires that PRM.NUCLEUSCH is defined.
%
%   MPRM.GETMINIMAMETHOD = 'AUTOMATED' (default): Automated detection of cell markers 
%   from the surface staining channel.
%   
%   Specific settings for CELLSEGMENTATION:
%
%   PRM.STARTPLANE 
%   Starts segmentation at a higher plane than the bottom
%   plane. This is useful is the signal is very weak at the bottom, which
%   can mislead the segmentation. Default value PRM.STARTPLANE = 1.
%
%   PRM.sEGMCH 
%   Defines the number of the segmentation channel. Default
%   PRM.SEGMCH = 1.
%
%   PRM.NUCLEUSCH 
%   Defines the number of the nucleus channel. Default
%   PRM.NUCLEUSCH = [] since there is no nucleus image available.
%   
%   PRM.SUBTRACT 
%   Defines whether the nucleus channel is subtracted (1) or not (0) from the
%   surface stain (in case of cross talk). Default value PRM.SUBTRACT = 1.
%
%   PRM.FORMAT
%   Defines the image format to be loaded and saved after processing. Can
%   be either 'mat' or 'tif'. Remember to also set
%   prm.nchannel when using option 'tif'. Default is 'mat'. 
%  
%   PRM.NCHANNEL
%   Defines the number of channels in your data. It is only used for reading the
%   multiple tif files (with prm.format = 'tif'). Default is 1 (the
%   segmentation channel and no nucleus channel).
%
%   The results from CELLSEGMENTATION can be visualized by using VIEWSEGM
%
%
%   Example
%   -------
%
%   segmcell('.',1:10,1,20,1,10,100,'prmfilesegmsurfwat')
%
%
%   Example
%   -------
%
%   folder = {'folder1','folder2','folder3'};    
%   segmcell(folder,1:10,1,20,1,10,100,'prmfilesegmsurfwat')
% 
%
% See also viewsegm

prm = [];

folder = varargin{1};
prm.stackstart = varargin{2};
prm.stackstop = varargin{3};
prm.planestart = varargin{4};
prm.planestop = varargin{5};
prm.minvolfull = varargin{6};
prm.maxvolfull = varargin{7};
prm.prmfile = varargin{8};
for i = 9 : 2 : nargin
    varname = varargin{i};
    switch varname
        case 'prm'
            prmincmd = varargin{i+1};
        otherwise
            error(['Wrong option to ' mfilename]);
    end;
end;

prm.stacks = prm.stackstart : prm.stackstop;
prm.nstacks = numel(prm.stacks);
prm.format = 'mat';

prmsegm.segmstart = 1;
prmsegm.subtract = 1;
prmsegm.segmch = 1;
prmsegm.nucleusch = [];
prmsegm.getminimamethod = 'automated';
prmsegm.h = [0.5 0.5 1.5];
prmsegm.method = 'segmsurfwat';
prmsegm.subtract = 1;
prmsegm.gpu = 0;


% the number of folders and make into cell array
if iscell(folder)
    nfolder = length(folder);
else
    nfolder = 1;
    folder = mat2celldirect(folder);
end;

%
% Load the parameter file if one is given
%

prmfile = prm.prmfile;
if ~isempty(prmfile)    
    isfolder = ~isempty(strfind(prmfile,'/'));
    [a b c] = fileparts(prmfile);
    if isfolder        
        A = pwd;
        cd(a);        
    end;    
    msg = sprintf('Reading from parameter file %s',prmfile);
    disp(msg);
    cmd = [b c];
    prmin = eval(cmd);

    if isfolder
        cd(A);
    end;

    % merge the structs by overriding the existing values
    prmsegm = mergestruct(prmsegm,prmin);

end;

msg = ['This is ' mfilename ' using settings'];
disp(msg);
printstructscreen(prm);

% from command line, these have priority
prmsegm = mergestruct(prmsegm,prmincmd);
if numel(prmsegm.segmstart) == 1
    prmsegm.segmstart = prmsegm.segmstart*ones(prm.stacks(end),nfolder);
end;

% in prmsegm we only want scalar settings, for this image only
prmsegm.segmstart(prmsegm.segmstart == 0) = NaN;
prm.segmstart = prmsegm.segmstart;
prmsegm = rmfield(prmsegm,'segmstart');

% % fill out the empty slots of NaN's
% d = size(prm.startplane);
% s = prm.startplane;
% prm.startplane = NaN(prm.stacks(end),nfolder);
% prm.startplane(1:d(1),1:d(2)) = s;
% d = size(prm.stopplane);
% s = prm.stopplane;
% prm.stopplane = NaN(prm.stacks(end),nfolder);
% prm.stopplane(1:d(1),1:d(2)) = s;

%
% Segmentation, loop over folders
%
for i = 1 : nfolder

    % segmentation in this folder
    segmentfolder(folder{i},prm,prmsegm);

end;%for


%-------------------------------------


function [eval] = segmentfolder(folderhere,prm,prmsegm)

nstacks = length(prm.stacks);
% loop over all images
for i = 1 : nstacks
    
    % the stack here
    stack = prm.stacks(i);
    name = [folderhere '/'  'stack' int2str(stack)];    

    % get the image
    [im,imsegm,imnucl,prmsegm] = getim(name,prm,prmsegm,i);

    if isempty(im)
        msg = ['Could not load ' name ', continuing'];
        disp(msg);
        continue;
    end;
        
    dim = size(im);
    if numel(dim) == 2
        dim = [dim 1];
    end;
    dim3 = dim(1:3);

    if dim(1) == 0 || dim(2) == 0 || dim(3) == 0        
        error('Wrong input of empty image');
        continue;
    end;
    if prmsegm.segmstart > dim(3) 
        msg = ['Segmentation plane is above stack, continuing'];
        disp(msg);
        continue;
    end;
    
    % information to return
    info = [];  
    

    %
    % Segmentation of surface stained cells
    %
    if isequal(prmsegm.method,'segmsurfwat')       
            
        % nucleus channel for markers
        if isequal(prmsegm.getminimamethod,'nucleus') && prmsegm.subtract  
            % subtract nucleus from WGA
            if dim(3) == 1
                filt = fspecial('gaussian',3,1);
                b1 = imfilter(imsegm,filt,'replicate');
                b2 = imfilter(imnucl,filt,'replicate');
            else
                b1 = smooth3(imsegm,'gaussian',[3 3 3]);
                b2 = smooth3(imnucl,'gaussian',[3 3 3]);
            end;
            
            % for watershed we use subtracted image                
            load ball2;se = getball(ball,2,1);            
            for j = 0.1:0.1:1.5
                % How much to subtract? Try iterating until values get
                % negative or loop is done
                imsegmtemp = b1 - j*b2;
                vol = imsegmtemp < 0;
                vol = imopen(vol,se);                    
                vol = sum(vol(:));
                % more than just noise?
                if vol > 0.10*numel(imsegm)                    
                    break;
                end;
                imsegm = imsegmtemp;                
            end;

            clear imsegmtemp
            imsegm(imsegm < 0) = 0;

        end;       

        % so we dont need to repeat plane for this image
        imsegminput = imsegm;
        
         % dont take the lowest planes
        
        if isequal(prmsegm.getminimamethod,'minimacell')
            % check for manual labels in the lower planes; we include them
            % if they exist! Thus we may change prmsegm.segmstart!!           
            for j = 1 : prmsegm.segmstart                    
               imhere = minima(:,:,j);
               if ~isempty(find(imhere,1));
                   prmsegm.segmstart = j;
                   break;
               end;
            end;
        end;

        msg = ['Segmentation from plane ' int2str(prmsegm.segmstart)];
        disp(msg);

%         if isequal(prmsegm.minimamethod,'minimacell')
%             minima = minima(:,:,prmsegm.segmstart:end);
%             minimacell = minimacell(:,:,prmsegm.segmstart:end);
%         elseif isequal(prmsegm.minimamethod,'nucleus')
%             imnucl = imnucl(:,:,prmsegm.segmstart:end);
%         end;         

        % cut the image
        imsegm = imsegm(:,:,prmsegm.segmstart:end);

        % manual minima given
        if isequal(prmsegm.getminimamethod,'minimacell')
            % load manual minima
            minimaload = [folderhere '/' 'stack' int2str(stack) '-mask.mat'];
            D = load(minimaload);            
            minima = D.minima;
            minimacell = D.minimacell;

            % cut the minima
            minima = minima(:,:,prmsegm.segmstart:end,:);
            minimacell = minimacell(:,:,prmsegm.segmstart:end,:);
            
            % check for dimension of loaded minima
            dimhere = size(minima);
            if numel(dimhere) == 2
                dimhere = [dimhere 1];
            end;
            if ~isequal(dimhere,dim3)
                error('Wrong dimension of images');
            end;
            
            % fill holes if some of the objects are not filled
            for j = 1 : dim(3)
                minima(:,:,j) = imfill(minima(:,:,j),'holes');
                minimacell(:,:,j) = imfill(minimacell(:,:,j),'holes');
            end;
        elseif isequal(prmsegm.getminimamethod,'nucleus')
            % cut nucleus image
            imnucl = imnucl(:,:,prmsegm.segmstart:end);
        end;

       
        %
        % Segmentation    
        %
        if isequal(prmsegm.getminimamethod,'minimacell')
            [cellbw,wat,imsegm,minima,minimacell,info] = ...
                segmsurfwat(imsegm,prm.minvolfull,prm.maxvolfull,'prm',prmsegm,'minima',minima,'minimacell',minimacell);                        
            clear temp;
        elseif isequal(prmsegm.getminimamethod,'nucleus')
            [cellbw,wat,imsegm,minima,minimacell,info] = ...
                segmsurfwat(imsegm,prm.minvolfull,prm.maxvolfull,'prm',prmsegm,'imnucleus',imnucl);                                        
        elseif isequal(prmsegm.getminimamethod,'automated')                  
            [cellbw,wat,imsegm,minima,minimacell,info] = ...
                segmsurfwat(imsegm,prm.minvolfull,prm.maxvolfull,'prm',prmsegm);                        
        else
            error('Wrong minima method: GETMINIMAMETHOD');
        end;
        
    elseif isequal(prmsegm.method,'segmct')
        
        
        % dont take the lowest planes
        if prm.segmplane > 1            
            msg = ['Segmentation from plane ' int2str(prm.segmstart)];
            disp(msg);
            imsegm = imsegm(:,:,prm.segmstart:end);
        end;

        % CTG staining        
        [cellbw,wat] = segmct(imsegm,prm.minvolfull,prm.maxvolfull);        
        minima = [];
    else
        error('Wrong option in METHOD');
    end;
            
    
    if prmsegm.segmstart > 1
        msg = ['Adding lower planes'];
        disp(msg);
        
        
        % fix the lower planes              
        cellbw = repeatplane(cellbw,prmsegm.segmstart);
        wat = repeatplane(wat,prmsegm.segmstart);                
        imsegm = imsegminput;
        m = zeros(dim(1:3));
        m(:,:,prmsegm.segmstart:end) = minima;
        minima = m;        
        if ~isempty(minimacell)
            m = zeros(dim(1:3));
            m(:,:,prmsegm.segmstart:end) = minimacell;
            minimacell = m;
        end;
    end;
    

    
    info.prm = prmsegm;
    wat = double(wat);    
  
    if isequal(prm.format,'mat') || isequal(prm.format,'all')
        % save variables
        name = [folderhere '/' 'stack' int2str(stack) '-segm.mat' ];
        save(name,'cellbw','wat','im','imsegm','minima','minimacell','info','-V6')    
        disp(sprintf('Saving %s\n',name))
    end;
    if isequal(prm.format,'tif')|| isequal(prm.format,'all')
        name = [folderhere '/' 'stack' int2str(stack) '-cellbw.tif' ];
        imwritemulttif(name,cellbw);
        name = [folderhere '/' 'stack' int2str(stack) '-wat.tif' ];
        imwritemulttif(name,wat);
        name = [folderhere '/' 'stack' int2str(stack) '-imsegm.tif' ];
        imwritemulttif(name,imsegm);
        name = [folderhere '/' 'stack' int2str(stack) '-minima.tif' ];
        imwritemulttif(name,minima);
        if ~isempty(minimacell)
            name = [folderhere '/' 'stack' int2str(stack) '-minimacell.tif' ];
            imwritemulttif(name,minimacell);
        end;
        
    end;

end%for


%---------------------------------

function [imnew] = repeatplane(im,start)

dim = size(im);
add = start-1;
dimnew = dim;
dimnew(3) = dimnew(3) + add;
imnew = zeros(dimnew);
imnew(:,:,start:end) = im;
for i = 1 : start - 1
    imnew(:,:,i) = im(:,:,1);
end;

%-------------------------------------------------

function [im,imsegm,imnucl,prmsegm,prm] = getim(name,prm,prmsegm,i)


try
    if isequal(prm.format,'mat')
        D = load(name);
    elseif isequal(prm.format,'tif')
        namehere = [name '.tif'];        
        D.im = imreadmulttif(namehere);
        D.im = reordermultipletif(D.im,prm.nchannel);
    else
        error('Wrong format');
    end;
    disp(sprintf('Segmentation of %s using channel %i',name,prmsegm.segmch));
catch ME1
    disp(sprintf('Unable to load %s',name));
    disp(ME1.message)
    im = [];
    imsegm = [];
    imnucl = [];
    return;
end;
im = D.im;

try
    % in microns
    prmsegm.h = D.h;
catch
    % use default
end;

clear D;
dim = size(im);

% start and stop planes in this stack
prmsegm.planestart = max(prm.planestart,1);
prmsegm.planestop = min(prm.planestop,dim(3));

% ALL planes in stack
prmsegm.plane = prmsegm.planestart:prmsegm.planestop;

% all channels
im = im(:,:,prmsegm.plane,:);

% the image to segment
imsegm = im(:,:,:,prmsegm.segmch);

% the nucleus image
imnucl = [];
if ~isempty(prmsegm.nucleusch)
    imnucl = im(:,:,:,prmsegm.nucleusch);
end;

% where to start and stop segmentation in this stack
prmsegm.segmstart = max(1,prm.segmstart(i));


%----------------------------------------
function [manualminima] = getmanualminima(folder,prm)

nfolder = size(folder,1);
for i = 1 : nfolder
    folderhere = deblank(folder(i,:));

    % only the first stack
    stackhere = prm.stacks(1);
    % get the image to segment
    name = [ folderhere '/' 'stack' int2str(stackhere)];               
    [im,imsegm] = getim(name,prm);

    [M N O] = size(im);
    mid = round(O/2);
    % draw minima on image
    bw = drawimage(imsegm(:,:,mid),1);
    
    %
    % Make background
    %
    bck = bckminima(imsegm,bw);

    % assign to minima
    bw(bck == 1) = 1;

    % assign
    manualminima(:,:,mid,i) = bw;                                                

end %for


%-------------------------------------------------------


function [bck] = bckminima(imsegm,bw)


[M N O] = size(imsegm);
% mid = round(O/2);

if isempty(find(bw,1))
    msg = 'Using default settings in BCKMINIMA in CELLSEGMENTATION';
    warning(msg);
    maxminima = 100;
    numpixminima = 1000;
    meanbw = mean(imsegm(:));
else
    [dim faser] = bwsize(bw);
    maxminima = max(dim);
    % background is probably lower intensities than MINIMA, use this here
    meanbw = mean(imsegm(bw == 1));
    numpixminima = length(find(bw == 1));
end;


th = 0.3;
% need to do iterativt to get a proper value, sometimes we have im,
% sometimes minkappa


while 1
    bck = lt(imsegm,th*meanbw);
    numpixbck = length(find(bck));
    % was at 4!!
    if numpixbck > 4*numpixminima || th >= 2
        break;
    end;
    th = th + 0.1;    
end;

% must close to get away all small holes
dim = 5;
se = strel('disk',dim);
bck = imclose(bck,se);

% must have no overlap to the given minima
load ball3;se = getball(ball,3,1);
bck(imdilate(bw,se) == 1) = 0;

% to remove small parts
se = strel('disk',dim+2);
bck = imopen(bck,se);


%     % fill the holes
%     minima = fillholes(minima,0,numpixlow,conn);

%     % erode
%     dim = 8;
%     se = strel('disk',dim);
%     bck = imerode(bck,se);

% dont want small background minima
bck = bwareaopen(bck,maxminima,4);
        
 
   