cell_refinement.m

function seg = cell_refinement(seg_raw,img,basic_t,area_min,area_max)

% dx = [0  0  0  -1  1 -1 -1 1  1];
% dy = [0 -1  1   0  0 -1  1 1 -1];

dx = [0  0  0  -1  1 ];
dy = [0 -1  1   0  0 ];

img = double(img);
s = size(seg_raw);
cell_N = max(seg_raw(:));

I_peaks = zeros(1,cell_N);
p_peaks = zeros(cell_N,2);

for current_cell = 1:cell_N   
   
    pixels = find(seg_raw == current_cell);    
    if isempty(pixels)
        continue
    end
    [I_m,I_m_ind] = max(img(pixels));
    [I_m_r,I_m_c] = ind2sub(s,pixels(I_m_ind));
    
    I_peaks(current_cell) = I_m;
    p_peaks(current_cell,:) = [I_m_r,I_m_c];

end

[~, ord] = sort(I_peaks,'ascend');

seg = zeros(size(img));

for p = 1:cell_N   
    
    current_cell = ord(p);

    if I_peaks(current_cell) == 0
        continue
    end
    
    edge_I = max(I_peaks(current_cell) * 0.5,basic_t);
    
    list = zeros(1000,2);
    
    pixels = find(seg_raw == current_cell);  
%     seg(pixels) = current_cell;
    [rr, cc] = ind2sub(s,pixels);
    
    l = length(rr);
    list(1:l,:) =  horzcat(rr(:),cc(:));
    
    
    % the flood-fill procedure that marks the voxels with particle number
    while l > 0
        
        % take the voxel from the tail of the stack
        pixel = list(l,:);
        x = pixel(1);
        y = pixel(2);
        l = l - 1;
        
        % look around the current voxel in the neighbors
        for search_voxel = 1 : size(dx,2)
            
            % the indices of the surrounding voxel
            xx = x + dx(search_voxel);
            yy = y + dy(search_voxel);
            
            if (xx < 1) || (xx > s(1)) || (yy < 1) || (yy > s(2))
                continue
            end;
            
            % if this voxel is brighter than the basic threshold
            % and dimmer than the current voxel
            if (img(xx,yy) >= edge_I) && (img(xx,yy) <= img(x,y)) %I_peaks(current_cell)) %<= img(x,y)) %
                
                % if this voxel has not been visited before
                if (seg(xx,yy) == 0)
                    
                    % mark this voxel with current particle number
                    seg(xx,yy) = current_cell;
                    
                    % place this voxel into the stack
                    l = l + 1;
                    list(l,:) = [xx,yy];
                    
                end
                
                % if this voxel has been visited from other cores
                if  (seg(xx,yy) > 0) && (seg(xx,yy) ~= current_cell)
                    
                    % this voxel doesn't belong to a particle core
                    % since it is under the influece of more than one particle
                    seg(xx,yy) = -1;
                    
                end
            end
        end
    end

end

% mark all the non-core voxels with zero
seg(seg == -1) = 0;
for current_cell = 1 : cell_N
    
    BW = imfill((seg == current_cell),'holes');
    seg(find(BW)) = 0;
    
    L = bwlabel(BW);
    areas = cell2mat(struct2cell(regionprops(BW,'area')));
    [~,max_ind] = max(areas);

    if isempty(max_ind)
        continue
    end
    BW = (L == max_ind(1));
    
    ind = find(BW);
    pixels_count = length(ind);
    
    if (pixels_count >= area_min) && (pixels_count <= area_max)
        seg(find(BW)) = current_cell;
    end   
      
end