Matlab: Matrix Neighbor

I have a large number of images that I have split into segments so that their matrices look like this:

img = [ 1 1 1 1 1 2 2 2 3 3 3 3  
        1 1 1 1 2 2 2 2 2 3 3 3  
        1 1 1 4 4 4 2 2 2 3 3 3  
        5 5 5 5 5 5 5 2 2 3 3 3 ];

      

        
        
        
      

    

where each number represents a different area, and each area is randomly formed. Thus, in this case, area 1 has neighbors 2, 4 and 5, area 2 has neighbors 1, 3 and 4, and so on.

I extracted all the regions into separate cells and got the statistics (mean, variance, etc.) that I plan to use to combine regions with statistics within a certain tolerance. I am struggling to find an efficient way to get each region's neighbors for this merger to take place.

I have a terrible solution that takes a very long time even for a single image:

referenceImage = [ 1 1 1 1 1 2 2 2 3 3 3 3;
                    1 1 1 1 2 2 2 2 2 3 3 3;
                    1 1 1 4 4 4 2 2 2 3 3 3;
                    5 5 5 5 5 5 5 2 2 3 3 3];

% Wish to extract each region into a separate cell
lastSP = 5;
sps = 1:lastSP;
% Could be a way to vectorise the below loop but it escapes me
superPixels(lastSP) = struct('Indices', 0, 'Neighbours', 0);
% Split data into separate cells
parfor a = 1 : lastSP
    inds = find(referenceImage == sps(a));
    superPixels(a).Indices = inds;
end

szs = size(referenceImage); % Sizes of RGB Image
for a = 1 : lastSP + 1
    mask = zeros(szs(1), szs(2)); % Just bin mask wanted
    mask(superPixels(a).Indices) = 1; % Mark the region pixels as one
    mask = xor(bwmorph(mask, 'thicken'), mask); % Obtain the outlying regions

    inds = find(mask ==1); % Fetch the external region indices

    neighbours = []; % Have to dynamically grow neighbours matrix
    neigh = 1;  

    for b = 1 : length(inds)
        found = false;
        if ~isempty(neighbours) % Check neighbours first
            for c = 1 : length(neighbours)
                if any(superPixels(neighbours(c)).Indices == inds(b))
                    found = true;
                    break;
                end
            end
        end

        if ~found
           for c = 1 : lastSP + 1 % Check every other region
               if any(superPixels(c).Indices == inds(b))
                    neighbours(neigh) = c;
                    neigh = neigh + 1;
                    break;
                end
           end
        end            
    end  
    superPixels(a).Neighbours = neighbours;
end

      

        
        
        
      

    

I am wondering if this is really the best way to approach this problem. I know the very last loop is the main problem, but I can't think of another way to intelligently write this unless I go into and check the neighbors of known neighbors.

Any help or nudge in the right direction would be greatly appreciated; thank!