function [SimSeg,NumSeg,Sim,ChangeBackSeg]=sim_segmentation(Sim,varargin)
%--------------------------------------------------------------------------
% sim_segmentation function                                        ImBasic
% Description: Thresholding and image segmentation using the
%              Watershed Transform.
% Input  : - Set of images on which to perform segmentation.
%            Input can be either cell array of images, FITS images names
%            with wild cards, file lists, SIM and more. See images2sim.m
%            for options and help.
%          * Arbitrary number of pairs of arguments: ...,keyword,value,...
%            where keyword are one of the followings:
%            'Thresh' - The threshold value in units of the error.
%                       Default is 3.0    
%            'SubBack'- Subtract background prior to thresholding
%                       {true|false}. Default is true.
%            'BackIm' - Background image/s to subtract (see images2sim.m
%                       for options). This can be either a single image
%                       (to be used with all input images) or a list
%                       of images. If empty and if SubBack is true, then 
%                       will attempt to calculate the background using 
%                       sim_background.m. Default is empty.
%            'SimBackPar' - Cell array of parameters to pass to
%                       sim_background.m. Default is {}.
%            'CalcErr'- Calculate error prior to thresholding
%                       {true|false}. Default is true.
%            'ErrIm'  - Error image/s (see images2sim.m
%                       for options). This can be either a single image
%                       (to be used with all input images) or a list
%                       of images. If empty and if CalcErr is true, 
%                       then will attempt to calculate the error using 
%                       sim_std.m. Default is empty.
%            'SimStdPar' - Cell array of parameters to pass to
%                       sim_std.m. Default is {}.
%            'SegBackOne' - In case the background segment is not 1
%                       according to the watershed transform, change it               
%                       to one {true|false}. default is true.   
%            'FeedBackB1' - A massage appears in case the background segment
%                       is changed to 1 after the watershed transform. 
%            'ImageField' - The field name in the structure image (Sim)
%                      containing the image to segment. Default is 'Im'.
%            'BackImField' - Background field name in the SIM structure
%                       array. Default is 'BackIm'.
%            'ErrImField' -  Error field name in the SIM structure array.
%                       Default is 'ErrIm'.
%            'CopyHead' - Copy header from original image {'y' | 'n'}.
%                       Default is 'y'.
%            'AddHead' - Cell array with 3 columns containing additional
%                       keywords to be added to the header.
%                       See cell_fitshead_addkey.m for header structure
%                       information. Default is empty matrix.
%            'OutSIM' - Output is a SIM class (true) or a structure array
%                       (false). Default is true.
%            --- Additional parameters
%            Any additional key,val, that are recognized by one of the
%            following programs:
%            images2sim.m
% Output : - Structure array of segmented images (SIM).
%          - vector containing the number of segments in each image
%            (i.e.number of extracted objects).
%          - The Structure array of images (SIM) - the original / modified
%            one (after background substruction and Error calculation).
%          - In case a replacement is needed after the watershed transform
%            to change the background segment number to 1, this variable
%            is true, otherwise it is false. 
% Tested : Matlab R2013b
%     By : Tali Engel                      Dec 2014
%    URL : http://weizmann.ac.il/home/eofek/matlab/
% Example: [Sim_seg,num_seg]=sim_segmentation('lred0126.fits');
% Reliable: 
%--------------------------------------------------------------------------
FunName = mfilename;
ChangeBackSeg = false;

ImageField  = 'Im';
%HeaderField = 'Header';
BackImField = 'BackIm';
ErrImField  = 'ErrIm';

%--- default values ---
DefV.Thresh      = 3.0;
DefV.SubBack     = true;
DefV.BackIm      = [];
DefV.SimBackPar  = {};
DefV.CalcErr     = true;
DefV.ErrIm       = [];
DefV.SimStdPar   = {};
DefV.SegBackOne  = true;
DefV.FeedBackB1  = true;
DefV.ImageField  = ImageField;
DefV.BackImField = BackImField;
DefV.ErrImField  = ErrImField;
DefV.CopyHead    = 'y';
DefV.AddHead     = {};
DefV.OutSIM      = true;


InPar = set_varargin_keyval(DefV,'n','use',varargin{:});

ImageField  = InPar.ImageField;
BackImField = InPar.BackImField;
ErrImField  = InPar.ErrImField;


%--- read images ---
Sim = images2sim(Sim,varargin{:});
Nim = numel(Sim);    % number of images


% read image noise 
if (~isempty(InPar.ErrIm) && InPar.CalcErr)
    % read Error images
    SimErr = images2sim(InPar.ErrIm,varargin{:});
    Nstdim = numel(SimErr);
end

% read image background 
if (~isempty(InPar.BackIm) && InPar.SubBack)
    % read Background images
    SimBack = images2sim(InPar.BackIm,varargin{:});
    Nbckim  = numel(SimBack);
end


if (InPar.OutSIM),
    SimSeg = simdef(Nim); %SIM;    % output is of SIM class
else
    SimSeg = struct(ImageField,cell(Nim,1));
end


% initialization
NumSeg = zeros(Nim,1);

%--- Go over images, threshold and segment ---
for Iim=1:1:Nim,    
    % estimate image noise prior to background subtraction
    if (InPar.CalcErr),    
        if (isempty(InPar.ErrIm)), 
            % Error image is not available, 
            % calculate std image
            Sim(Iim) = sim_std(Sim(Iim),InPar.SimStdPar{:});
        else
            Sim(Iim).(ErrImField) = SimErr(min(Iim,Nstdim)).(ImageField);
        end
    end
    % background subtraction
    if (InPar.SubBack),
        % subtract background from SIM
        if (isempty(InPar.BackIm)),
            % background image is not available,
            % calculate background
            Sim(Iim) = sim_background(Sim(Iim),InPar.SimBackPar{:},'SubBack','true');
        else
            Sim(Iim).(BackImField) = SimBack(min(Iim,Nbckim)).(ImageField);
            Sim(Iim).(ImageField)  = Sim(Iim).(ImageField) - Sim(Iim).(BackImField);
        end
    end
                   
    % thresholding
    % (pixel > thresh_val => pixel=1,  pixel < thresh_val => pixel=0)    
    BinImage = zeros(size(Sim(Iim).(ImageField)));
    BinImage(Sim(Iim).(ImageField) > (InPar.Thresh.*Sim(Iim).(ErrImField))) = 1;
    
    % --- segmentation of the binary image ---
    % the distance transform of the complement of the binary image
    % complement the distance transform
    D = -bwdist(~BinImage);
    % force pixels that don't belong to the objects to be -Inf
    D(~BinImage) = -Inf;
    % apply the watershed transform
    SimSeg(Iim).(ImageField) = watershed(D);
    % number of segments in each image
    NumSeg(Iim) = max(SimSeg(Iim).(ImageField)(:)) - 1;   % subtract 1 - for the background segment 
    
    % check that the background is represented by segment #1, and if not,
    % fix it
    if (InPar.SegBackOne)
        Seg_Im = SimSeg(Iim).(ImageField);
        N = NumSeg(Iim)+1;             %number of segments, including the background
        % find the segment that includes maximum numbers of pixels
        Hist = histc(Seg_Im(:),1:1:N);
        [~,Ipix] = max(Hist);
        if (Ipix>1)
            Seg_Im(Seg_Im==1) = N+1;   
            Seg_Im(Seg_Im==Ipix) = 1;
            Seg_Im(Seg_Im==(N+1)) = Ipix;            
            SimSeg(Iim).(ImageField) = Seg_Im;
            ChangeBackSeg = true;
            if (InPar.FeedBackB1)
                disp('The background segment number was changed to 1 artificially (after the watershed transform)');
            end
        else
            %do nothing
        end
    end
        
        
    %--- Update header ---
    % Update the SimSeg header (not the Sim header)
    SimSeg(Iim) = sim_update_head(SimSeg(Iim),'CopyHeader',InPar.CopyHead,...
                                        'AddHead',InPar.AddHead,...
                                        'Comments',{sprintf('Created by %s.m written by Tali Engel',FunName)},...
                                        'History',{sprintf('Threshold value: %d',InPar.Thresh)});
                                    
end

    
