function [seg,seg_arr,reg_info]=region_grow1(im,seg,criterion,connectivity)
  
% This function performs basic region growing on an image.
%
% The input arguments are the following:
% im: The image to be segmented
% seg: An initial labelling of the image, typically you would let
% this image be an image where every pixel has its own label
% criterion: The value of the similarity criterion that is used to decide
% between fusion of segments or not.
% connectivity: The connectivity of regions considered for fusion.
% 
% The output arguments are the following:
%
% seg: The final segmented image
% seg_arr:  A stack of all the segmented images from every
% iteration
% reg_info: A structure containing all the information about
% the different segments.
%
% Copyright Lars Aurdal, The Norwegian Computing Centre
  
% Get size of input image
  
rows=size(im,1);
cols=size(im,2);

% Here comes a slightly tricky part.  In order to keep the computation
% time down we need to use s data structure to keep track of information
% about the regions.  A MATLAB struct is what we need.  Each element in
% the reg_info array of structures holds the size in pixels of the
% element, the mean of the pixels assigned to the element as well as the
% row and column coordinates of all the pixels assigned to the region,
% this is nice to have in order to be able to quickly adress all the
% pixels belonging to a single element.

% Start by initialising

reg_info=struct('size',num2cell(zeros(1,max(max(seg)))),...  % Initially we don't know their size
                'crt',num2cell(zeros(1,max(max(seg)))),...   % ...nor their criterion value
                'row_coord',[],...                           % Coordinates go here
                'col_coord',[]);                             % and here

% Now loop over all the pixels in the initial segmented image and gather
% information about the existing regions

for m=1:rows
  for n=1:cols
    c_lab=seg(m,n);
    reg_info(c_lab).size      = reg_info(c_lab).size+1;
    
    % The following calculation of criterion corresponds to a calculation
    % of the segment means
    
    reg_info(c_lab).crt        = ((reg_info(c_lab).size-1)*reg_info(c_lab).crt+...
                                 im(m,n))/reg_info(c_lab).size;
    reg_info(c_lab).row_coord = [reg_info(c_lab).row_coord m];
    reg_info(c_lab).col_coord = [reg_info(c_lab).col_coord n];
  end
end

% Initially we assume that something will change in the first iteration

changed=1;

% Keep track of iterations in order to display this

iter=0;

% Now loop over the segmented image as long as something changes.

while(changed)
    
    % Give user some feedback
    
    disp(['Currently in iteration: ' num2str(iter+1)]);
    
    % Check if something changes in this iteration
            
    changed=0;
    
    % Now loop
    
    for m=1:rows
        for n=1:cols
        
            % Get the label of the current pixel considered
        
            c_lab=seg(m,n);
        
            % Get the coordinates of the current pixels neighbours.  This
            % can be done in 4 or 8 connectivity
                        
            if (connectivity==4)
                ngb=get4ngb(rows,cols,m,n);
            end
            if (connectivity==8)
                ngb=get8ngb(rows,cols,m,n);
            end
        
            % Loop over all these neighbours and check if they should
            % receive the same label as the pixel we are currently
            % examining
        
            for k=1:size(ngb,2)
                n_lab=seg(ngb(1,k),ngb(2,k));
                if(n_lab~=c_lab)
                    if(abs(reg_info(n_lab).crt-reg_info(c_lab).crt)<criterion)
                        
                        % First update the mean of the segment having the
                        % c_lab label
                        
                        reg_info(c_lab).crt=(reg_info(c_lab).size*reg_info(c_lab).crt+...
                                            reg_info(n_lab).size*reg_info(n_lab).crt)/...
                                            (reg_info(c_lab).size+reg_info(n_lab).size);
                                       
                        % Then we are ready to update the size
                        
                        reg_info(c_lab).size=reg_info(c_lab).size+reg_info(n_lab).size;
                        
                        % Then give all pixels with n_lab the c_lab label
                        
                        tmp_row_coord=reg_info(n_lab).row_coord;
                        tmp_col_coord=reg_info(n_lab).col_coord;
                        for p=1:reg_info(n_lab).size
                            seg(tmp_row_coord(p),tmp_col_coord(p))=c_lab;
                        end
                        
                        % Finally update the coordinates of the pixels
                        % belonging to this segment
                        
                        reg_info(c_lab).row_coord=[reg_info(c_lab).row_coord reg_info(n_lab).row_coord];
                        reg_info(c_lab).col_coord=[reg_info(c_lab).col_coord reg_info(n_lab).col_coord];
                        
                        % ...and mark it for deletion
                        
                        reg_info(n_lab).size=0;
                        
                        % If we have reached this point then we know that
                        % something changed in this iteration
                        
                        changed=1;
                    end
                end
            end
        end
    end
    iter=iter+1;
    
    % We will store the temporary results in seg_arr for demonstration
    % purposes, the first layer equals the initial segmentation

    seg_arr(:,:,iter)=seg;
    
end

% Clean up the reg_info structure

index=1;
for m=1:max(max(seg))
    if (reg_info(m).size>0)
        tmp_reg_info(index)=reg_info(m);
        index=index+1;
    end
end
reg_info=tmp_reg_info;