%
% function [neighbours] = find_neighbours(R)
%
% this function takes in a NxMxP matrix that has P layers. The first layer
% is the binary image of the whole organ, and the other layers are defined
% by the grid. This function returns a matrix that is PxP in size,
% neighbours is an adjacancy graph where each row corresponds to a region
% and its neighbours, i.e. row 1 is region 1 and any 1s in this row
% correspond to neighbours.
%
% Dr. A. I. Hanna (2005)
%
function [neighbours] = find_neighbours(R)
% This function assumes that R>0 is the whole organ shape

waitbar_h = waitbar(0, 'Finding neighbours, please wait...');
num_regions = max(R(:));

neighbours = zeros(num_regions-1);
for i=2:num_regions
    Q = R==i;
    Q = imdilate(Q, ones(5));
    [x, y] = find(Q==1);
    Q = Q(min(x):max(x), min(y):max(y), :);
    subR = R(min(x):max(x), min(y):max(y));
    subR = subR.*(subR & Q);
    n = get_significant_neighbours(subR);
    neigh_default = zeros(1, num_regions-1);
    n(n==0) = [];
    if ~isempty(n)
        neigh_default(n)=1;
        neighbours(i-1,:) = neigh_default;
    end
    waitbar(i/num_regions, waitbar_h);
end
close(waitbar_h);
%%%
%
%%%
function n = get_significant_neighbours(subR)
n = unique(subR(:));
n(n==0) = [];
n(n==i) = [];
reg_to_remove = [];
% here we want to remove those neighbours that don't really contribute much
for ii=1:length(n)
    reg_im = subR==n(ii);
    if sum(reg_im(:))<100
        reg_to_remove = cat(1, reg_to_remove, ii);
    end
end
n(reg_to_remove) = [];
n = n - 1;
return
