function [cellbw,wat,imsegm,minima,minimacell,info] = segmsurfwat(varargin)
% SEGMSURFWAT 3D segmentation of surface stained cells
%
%   [CELLBW,WAT,IMSEGM,MINIMA,MINIMACELL] = SEGMSURFWAT(IM,MINVOLFULL,MAXVOLFULL)
%   Finding cells in a surface stained image IM. MINVOLFULL and MAXVOLFULL are the lower
%   and upper expected volume of a cell in thousand cubic microns.
%   Returning the watershed image WAT, the classified cells in CELLBW,
%   the ridge image IMRIDGE, the markers MINIMA, information of parameters
%   in INFO and the smoothed image IMSM.
%
%   [CELLBW,WAT,IMSEGM,MINIMA,MINIMACELL,INFO] = SEGMSURFWAT(IM,MINVOL,MAXVOL,'prmfile',PRMFILE)
%   specifies a parameter file that can be used for segmentation
%
%   [CELLBW,WAT,IMSEGM,MINIMA,MINIMACELL,INFO] = SEGMSURFWAT(IM,MINVOL,MAXVOL,OPTIONNAME,OPTION)
%   specifies an OPTIONNAME for segmentation that can be of one of
%   these options:
%   PRM         : General parameters for SEGMSURFWAT
%   MINIMA      : Markers, must be same size as image
%   SMPRM       : Smoothing parameters
%   CPRM        : Classification parameters
%   MPRM        : Marker parameters
%   IMRIDGE     : Ridge enhanced image
%   MINIMACELL  : Markers for cells, for instance from manual drawings
%   IMNUCLEUS   : Nucleus image for defining markers
%
%   Example
%   -------
%
%   A segmentation with cells with volume of between 3000 and 50000 mcm3
%   [cellbw,wat,imsegmout,minima,minimacell,info] = segmsurfwat(im,3,50);
%
%
%   Example
%   -------
%
%   Specifying a particular method for smoothing
%   smprm.method = 1;
%   [cellbw,wat,imsegmout,minima,minimacell,info] = segmsurfwat(im,5,50,'smprm',smprm);
%
%   Example
%   -------
%
%   Segmentation in 3D. Using nucleus image for markers and classification.
%
%   load surfstain-and-nucleus-3D.mat
%   % Smoothing
%   smprm.method = 2;
%   % No ridge filtering
%   prm.filterridges = 0;
%   % Subtract the nucleus channel from the surface staining to reduce the
%   % cross talk effect.
%   imsegm = imsegm - imnucl;
%   [cellbw,wat,imsegmout,minima,minimacell,info] = ...
%    segmsurfwat(imsegm,3,50,'imnucleus',imnucl,'smprm',smprm,'prm',prm);
%   % Show results: image, nucleus image, markers, watershed image, cells
%   showall(imsegm,imnucl,minima,wat,cellbw)
%
%   Example SURFSTAIN_AND_NUCLEUS_2D
%   -------
%
%   The same as above but in 2D. Using nucleus image for markers and classification.
%
%   % load the data
%   load surfstain_and_nucleus_3D.mat
%   imsegm = imsegm(:,:,5);imnucl = imnucl(:,:,5);
%
%   % Smoothing
%   smprm.method = 2;
%
%   % No ridge filtering
%   prm.filterridges = 0;
%
%   % Subtract the nucleus channel from the surface staining to reduce the
%   % cross talk effect.
%   imsegm1 = imsegm;
%   imsegm = imsegm1 - imnucl;
%   [cellbw,wat,imsegmout,minima,minimacell,info] = ...
%    segmsurfwat(imsegm,3,80,'imnucleus',imnucl,'smprm',smprm,'prm',prm);
%
%   % Show results: image, nucleus image, markers, watershed image, cells
%   showall(imsegm,imnucl,minima,wat,cellbw)
%
% 	=======================================================================================
% 	CELLSEGM: Cell segmentation, Version 2012
% 	Copyright (c): Erlend Hodneland
%   Jonas Lies Vei 91,
%   5009 Bergen
%   Email: erlend.hodneland@biomed.uib.no
% 	=======================================================================================
%

printmsg(['This is ' upper(mfilename) ' for segmentation of cells']);



% image for segmentation
imsegm = varargin{1};

% cell diameters, here after!
minvolfull = 1000*varargin{2};
maxvolfull = 1000*varargin{3};
imsegmini = imsegm;
prmfile = [];
imsm = [];
imridge = [];
imnucl = [];

% for visualilzation
test = 0;
testfinal = 0;

%
% Standard parameters
%

% voxel size
prm.h = [0.5 0.5 1.5];

% adjust for not fully circular shape of 3D cells
prm.just = 0.9;

% filter ridges
prm.filterridges = 1;

% % dilation of strong pixels
% prm.dilation = 0;

% classify cells
prm.classify = 1;

%
% Minima parameters
%
mprm.level = 0.30;

%
% Smoothing parameters
%

% 0 : no filtering
% 1 : coherence enhancing diffusion CED
% 2 : directional coherence enhancing diffusion DCED
smprm.method = 2;

% 2D or 3D smoothing
smprm.smoothdim = 2;

% number of iterations of smoothing in CED
smprm.numit = 3;

% dimension image
dim = size(imsegm);
if numel(dim) == 2
    dim = [dim 1];
end;

%
% Eempty
%
mprm = [];
cprm = [];

%
% Watershed parameters
%

% To test intensity of wat lines
wprm.merging = 0;

% The signicance of watershed lines, stronger if larger, done on ridgeim!!
% if the raw image, then use much smaller, around 1
% set to 1.05-1.10
% wprm.thint = 1.06;
wprm.thintmerging = 1.3;

% The conexity measure for the merging, we should not allow a very much
% smaller convexity after merging
% thconv = 1 if no change in convexity
% thconv > 1 if increase in convexity
% thconv < 1 if decrease in convexity
wprm.thconvmerging = 1.0;


if maxvolfull < minvolfull
    error('The upper cell diameter is lower than the lower, check the input')
end;

% read from parameter file?
for i = 4 : 2 : nargin
    namehere = varargin{4};
    varhere = varargin{5};
    switch(namehere)
        case('prmfile')
            prmfile = varhere;
        otherwise
            % nothing
    end;
end;

%
% Load the parameter file if one is given
%
if ~isempty(prmfile)
    isfolder = ~isempty(findstr(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,mprmin,cprmin,wprmin,smprmin] = eval(cmd);

    if isfolder
        cd(A);
    end;
    % merge the structs by overriding the existing values
    prm = mergestruct(prm,prmin);
    mprm = mergestruct(mprm,mprmin);
    cprm = mergestruct(cprm,cprmin);
    wprm = mergestruct(wprm,wprmin);
    smprm = mergestruct(smprm,smprmin);
end;


%
% The command line has the highest priority, merge the variables here
%
for i = 4 : 2 : nargin

    % this variable
    namehere = varargin{i};
    varhere = varargin{i+1};

    switch namehere
        case 'minima'
            minima = varhere;
        case('imridge')
            imridge = varhere;
        case('smprm')
            prmin = varhere;
            smprm = mergestruct(smprm,prmin);
        case('mprm')
            prmin = varhere;
            mprm = mergestruct(mprm,prmin);
        case('wprm')
            prmin = varhere;
            wprm = mergestruct(wprm,prmin);
        case('prm')
            prmin = varhere;
            prm = mergestruct(prm,prmin);
        case('cprm')
            prmin = varhere;
            cprm = mergestruct(cprm,prmin);
        case('minimacell')
            minimacell = varhere;
        case('imnucleus')
            imnucl = varhere;
        otherwise
            msg = sprintf('Wrong option: %s',namehere);
            error(msg)
    end;
end;


%
% Computed settings
%

% get the ADJUSTED 3D cell size
[prm.minvol,prm.minvolvox,prm.maxvol,prm.maxvolvox] = ...
    cellsize(minvolfull,maxvolfull,prm.h,prm.just,dim(3));
cprm.minvolfull = minvolfull/1000;
cprm.maxvolfull = maxvolfull/1000;
mprm.minvolfull = minvolfull/1000;
mprm.maxvolfull = maxvolfull/1000;
wprm.minvolfull = minvolfull/1000;
wprm.maxvolfull = maxvolfull/1000;
mprm.minvolfull = minvolfull/1000;
mprm.maxvolfull = maxvolfull/1000;


% pixel size
mprm.h = prm.h;
cprm.h = prm.h;
smprm.h = prm.h;
wprm.h = prm.h;


disp('Using settings:');
printstructscreen(prm);

% if you give in a nucleus image you want to use it!
if ~isempty(imnucl)
    mprm.method = 'nucleus';
end;

%--------------------------------------------%
%-----------Cell segmentation start----------%
%--------------------------------------------%


%
% Types of smoothing
%

% coherence enhancing diffusion
if smprm.method == 1
    msg = 'Coherence enhancing diffusion';
    disp(msg);
    if smprm.smoothdim == 2
        for i = 1 : dim(3)
            imhere = imsegm(:,:,i);
            imsegm(:,:,i) = cohenhdiff(imhere,0.1,smprm.numit,0.0001);
        end;
    elseif smprm.smoothdim == 3
        % true 3D
        imsegm = cohenhdiff(imsegm,0.1,smprm.numit,0.0001,'opt','num');
    else
        error('Wrong option to SMPRM.SMOOTHDIM')
    end;
% Directional coherence enhancing diffusion
elseif smprm.method == 2
    msg = 'Directional coherence enhancement filtering';
    disp(msg);
    if smprm.smoothdim == 2
        for i = 1 : dim(3)
            imhere = imsegm(:,:,i);
            imsegm(:,:,i) = dircohenh(imhere,13,1,1,3);
        end;
    elseif smprm.smoothdim == 3
        % true 3D
         imsegm = dircohenh(imsegm,13,1,1,3);
    else
        error('Wrong option to SMPRM.SMOOTHDIM')
    end;

else
    imsm = imsegm;
    msg = 'No filtering in SEGMSURFWAT, FILTIM is assigned as IM';
    disp(msg);
end;
% load test3D;
% imsegm = b;


%
% Ridge enhancement
%
if prm.filterridges == 1
    disp('Ridge enhancement')
    imsegm = ridgeenhhessian(imsegm,prm.h);
else
    msg = 'No ridge enhancement';
    disp(msg);
end;

%
% Finding minima
%

if ~exist('minima','var')
    % find the minima
    if isempty(imnucl)
        [minima,minimacell,mprm] = findminima(imsegm,mprm);
    else
        [minima,minimacell,mprm] = findminima(imsegm,mprm,imnucl);
    end;
else
    disp('Not creating minima, using manual minima');
end;

if test == 1
    showall(imsegm,minima)
end;

%
% 3D watershed segmentation
%

% segmentation
disp('3D watershed segmentation')
% save test imsegm minima minimacell
% pause
Iimpose = imimposemin(imsegm,minima);
wat = watershed(Iimpose);
clear Iimpose

if wprm.merging == 1
    % Test the significance of watershed lines, the new version
    % optlog = 1 : merge the strong lines (for CT)
    % optlog = 2 : merge the weak lines (for WGA)
    % optint = 1 : use the region based merging
    % optint = 2 : use the snake based merging
    optlog = 2;
    optint = 2;
    disp('  Merging watershed regions')
    wat = mergefragments(wat,imsegm,wprm.thintmerging,wprm.thconvmerging,optlog,optint);

    % re segment after fooling around with the image
    disp('  Resegment watershed')
    wat = resegmentwat(wat,imsegm);

end;


if test == 1
    show(minima,10)
    show(wat,11)
end;


%
% Classify cells
%

if prm.classify
    if ~isempty(minimacell)
        [cellbw,info.classifycells] = classifycells(wat,imsegmini,cprm,minimacell);
    else
        [cellbw,info.classifycells] = classifycells(wat,imsegmini,cprm);
    end;
else
    warning('No cell classification')
    cellbw = zeros(dim);
end;

% % dilation of strong regions
% if prm.dilation
%     disp('Dilation of strong voxels')
%     cellbw = dilatestrong(cellbw,imsegm);
% else
%     disp('No dilation of strong pixels')
% end;

if test == 1 || testfinal == 1
    if dim(3) == 1
        show(imsegm,1);title('Image');
        show(wat,10);title('Wat');
        show(minima,11);title('Minima')
        show(minimacell,12);title('Minimacell')
        show(cellbw,13);title('Final segmentation');
        pause
    else
        showall(imsegm,cellbw,wat,minima)
    end;

end;

% return variables
info.prm = prm;
info.mprm = mprm;
info.wprm = wprm;
info.cprm = cprm;


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

This is SEGMSURFWAT for segmentation of cells

Index exceeds matrix dimensions.

Error in ==> segmsurfwat at 97
imsegm = varargin{1};

Published with MATLAB® 7.12