function im_seed = CERW_seeds(im_input,options)

im_input_0 = CERW_preprocess_image(im_input,options);

phan_update = false;

phan_radius = [10,10];
if isfield(options,'phan_radius')
    phan_radius(1) = options.phan_radius;
    phan_radius(2) = options.phan_radius;
end
if isfield(options.BG_options,'phan_radius')
    phan_radius(1) = options.BG_options.phan_radius;
end
if isfield(options.FG_options,'phan_radius')
    phan_radius(2) = options.FG_options.phan_radius;
end
phan_update = (phan_update || ~isequal(phan_radius(1),phan_radius(2)));

phan_p = [2,2];
if isfield(options,'phan_p')
    phan_p(1) = options.phan_p;
    phan_p(2) = options.phan_p;
end
if isfield(options.BG_options,'phan_p')
    phan_p(1) = options.BG_options.phan_p;
end
if isfield(options.FG_options,'phan_p')
    phan_p(2) = options.FG_options.phan_p;
end
phan_update = (phan_update || ~isequal(phan_p(1),phan_p(2)));


phan_r = [1,1];
if isfield(options,'phan_r')
    phan_r(1) = options.phan_r;
    phan_r(2) = options.phan_r;
end
if isfield(options.BG_options,'phan_r')
    phan_r(1) = options.BG_options.phan_r;
end
if isfield(options.FG_options,'phan_r')
    phan_r(2) = options.FG_options.phan_r;
end
phan_update = (phan_update || ~isequal(phan_r(1),phan_r(2)));


phan_restrict = [false,false];
if isfield(options,'phan_restrict')
    phan_restrict(1) = options.phan_restrict;
    phan_restrict(2) = options.phan_restrict;
end
if isfield(options.BG_options,'phan_restrict')
    phan_restrict(1) = options.BG_options.phan_restrict;
end
if isfield(options.FG_options,'phan_restrict')
    phan_restrict(2) = options.FG_options.phan_restrict;
end
phan_update = (phan_update || ~isequal(phan_restrict(1),phan_restrict(2)));


% generate background seeds
if(phan_restrict(1))
    im_phan = imphansalkar(im_input_0,phan_radius(1),phan_p(1),phan_r(1),'otsu');
else
    im_phan = imphansalkar(im_input_0,phan_radius(1),phan_p(1),phan_r(1));
end
cc = bwconncomp(im_phan,6);
num_pixels = cellfun(@numel,cc.PixelIdxList);
[~,max_idx] = max(num_pixels);
for i = 1:cc.NumObjects
    if i ~= max_idx
        im_phan(cc.PixelIdxList{i}) = 0;
    end
end
im_bg = dil_fil_ero(im_phan,options.BG_options.r_dil,options.BG_options.r_ero);


if isfield(options.FG_options,'method')
    m_no = options.FG_options.method;
    if m_no == 1 || m_no == 3
        im_fg = dil_fil_ero(im_phan,options.FG_options.r_dil,options.FG_options.r_ero);
    end
    if m_no == 2 || m_no == 3
        im_input = CERW_preprocess_image(im_input,options.FG_options);
        
            if(phan_restrict(2))
                im_phan = imphansalkar(im_input,phan_radius(2),phan_p(2),phan_r(2),'otsu');
            else
                im_phan = imphansalkar(im_input,phan_radius(2),phan_p(2),phan_r(2));
            end
            cc = bwconncomp(im_phan,6);
            num_pixels = cellfun(@numel,cc.PixelIdxList);
            [~,max_idx] = max(num_pixels);
            for i = 1:cc.NumObjects
                if i ~= max_idx
                    im_phan(cc.PixelIdxList{i}) = 0;
                end
            end
        if m_no == 2
            im_fg = local_max(im_input,im_phan);
        else
            im_fg = max(im_fg,local_max(im_input,im_phan));
        end
    end
else
    im_fg = dil_fil_ero(im_phan,options.FG_options.r_dil,options.FG_options.r_ero);
end

im_seed = ((im_fg>0) + (im_bg>0)) / 2;

end

%
% Dilate and erode

function im_seed = dil_fil_ero(im_phan,r_dil,r_ero)
% faster than inbuilt dilate and erode
im_phan = (bwdist(im_phan>0) < r_dil);
im_phan = imfill(im_phan,'holes');
im_phan = (bwdist(im_phan==0) > r_ero);
im_seed = zeros(size(im_phan));
%image edges are always background and non-foreground
im_seed(2:(end-1),2:(end-1),2:(end-1)) = im_phan(2:(end-1),2:(end-1),2:(end-1));
end

% Maximal surfaces method

function im_seed = local_max(im_input,im_phan)
% save original size for output
sz0 = size(im_input);

% resize for faster processing
h_sum = squeeze(sum(im_phan,[2,3]));
w_sum = squeeze(sum(im_phan,[1,3]));
d_sum = squeeze(sum(im_phan,[1,2]));
h_min = find(h_sum>0,1);
w_min = find(w_sum>0,1);
d_min = find(d_sum>0,1);
h_max = find(h_sum>0,1,'last');
w_max = find(w_sum>0,1,'last');
d_max = find(d_sum>0,1,'last');
im_input = im_input(h_min:h_max,w_min:w_max,d_min:d_max);
im_phan = im_phan(h_min:h_max,w_min:w_max,d_min:d_max);
[h1,w1,d1] = size(im_input);
vol1 = h1*w1*d1;

% Create gradient image

im_grad = zeros(h1, w1, d1, 4);
[im_grad(:,:,:,2), im_grad(:,:,:,1), im_grad(:,:,:,3)] = imgradientxyz(im_input);
grad2 = im_grad(:,:,:,1:3).*im_grad(:,:,:,1:3);
im_grad(:,:,:,4) = sqrt(sum(grad2,4));
clear grad2;

% Find second derivatives for turning point classification

[gxx,~,~] = imgradientxyz(im_grad(:,:,:,2));
[~,gyy,~] = imgradientxyz(im_grad(:,:,:,1));
[~,~,gzz] = imgradientxyz(im_grad(:,:,:,3));
im_grad2 = (gxx + gyy + gzz) / 44^2;
for d = 1:3
    im_grad2 = movmean(im_grad2,3,d);
end

clear gxx gyy gzz;

% Set gradient direction to the locally maximal direction

im_grad_max = im_grad;
for d = 1:3
    for i = -1:2:1
        im_grad_shift = circshift(im_grad_max,i,d);
        new_max_idx = find(im_grad_shift(:,:,:,4) > im_grad_max(:,:,:,4));
        for j = 0:3
            im_grad_max(new_max_idx + j * vol1) = im_grad_shift(new_max_idx + j * vol1);
        end
    end
end
clear im_grad_shift new_max_idx

% Non-minimum gradient suppression

% TODO: vectorize this
im_seed_temp = zeros(h1, w1, d1);
for r = 2:h1-1
    for c = 2:w1-1
        for z = 2:d1-1
            % restrict by phansalkar threshold and ensure maxima
            if im_phan(r,c,z) == 1 && im_grad2(r,c,z) < 0
                [g_max, i_max] = max(abs(im_grad_max(r,c,z,1:3)));
                if g_max > 0 % only consider points with assigned direction
                    % scale gradient vector so that the longest axis is
                    % equal to 1.
                    g_loc = [0, 0, 0];
                    for i = 1:3
                        g_loc(i) = im_grad_max(r,c,z,i)/g_max;
                    end
                    % list axes with longest first
                    i_list = [i_max, mod(i_max,3) + 1, mod(i_max + 1,3) + 1];
                    % compute average gradient intensity along g_loc in
                    % each direction and save the minimum
                    v_min = im_grad(r,c,z,4) + 1; % only care if the minimum is less than grad
                    for d = -1:2:1
                        p_cen = [r, c, z] + d * g_loc;
                        p = zeros(4,3);
                        for k = 1:2
                            % all corners have same i_max coordinate as
                            % p_cen
                            p(k,i_list(1)) = floor(p_cen(i_list(1)));
                            p(k + 2,i_list(1)) = ceil(p_cen(i_list(1)));
                            % first two corners share a lower coordinate
                            p(k,i_list(2)) = floor(p_cen(i_list(2)));
                            % first and third corners share a lower
                            % coordinate
                            p(2 * k - 1,i_list(3)) = floor(p_cen(i_list(3)));
                            % last two corners share an upper
                            % coordinate
                            p(k + 2,i_list(2)) = ceil(p_cen(i_list(2)));
                            % second and fourth corners share an upper
                            % coordinate
                            p(2 * k,i_list(3)) = ceil(p_cen(i_list(3)));
                        end
                        %bilinear interpolation
                        v_sum = 0;
                        for k = 1:4
                            p_diff = abs(p(k,:) - p_cen);
                            v_sum = v_sum + (1 - p_diff(i_list(2))) * (1 - p_diff(i_list(3))) * im_grad(p(k,1),p(k,2),p(k,3),4);
                        end
                        v_min = min([v_min, v_sum]);
                    end
                    if(im_grad(r,c,z,4) < v_min)
                        im_seed_temp(r,c,z) = 1;
                    end
                end
            end
        end
    end
end
clear im_grad_max im_grad im_phan_FG

cc = bwconncomp(im_seed_temp,6);
num_pixels = cellfun(@numel,cc.PixelIdxList);
[~,max_idx] = max(num_pixels);
max_mean = mean(im_input(cc.PixelIdxList{max_idx}));
for i = 1:cc.NumObjects
    if (mean(im_input(cc.PixelIdxList{i})) < max_mean) && i ~= max_idx
        im_seed_temp(cc.PixelIdxList{i}) = 0;
    end
end

im_seed = zeros(sz0);
im_seed(h_min:h_max,w_min:w_max,d_min:d_max) = im_seed_temp;

end