function [seg, Am, area] = mergeregions(s1, s2, seg, Am, area)
    
    if s1==s2
        fprintf('s1 == s2!\n')
        return
    end
    
    % The area of s1 is now that of s1 and s2
    area(s1) = area(s1)+area(s2);
    area(s2) = NaN;
    
    % s1 inherits the adjacancy matrix entries of s2
    Am(s1,:) = Am(s1,:) | Am(s2,:);
    Am(:,s1) = Am(:,s1) | Am(:,s2);        
    
            Am(s1,s1) = 0;  % Ensure s1 is not connected to itself 
    
    % Relabel s2 with s1 in the segment image
    seg(seg==s2) = s1;
    
    