function StatArray = adjacency_stats2 (Im)
%
% Function to calculate the threshold adjacency statistics for a
% thresholded image.  Function takes a thresholded image (of a ny intiger,
% float or logic type, so long as 0 = background and non-0 = signal), and 
% for each pixel above threshold counts the number of neighbouring pixels 
% also above threshold.
%
% Based on: Hamilton et al. (2007).  BMC Bioinformatics Vol8 #110.
%
% Modifications from original algorithm:
%   -Works with images of any bit-depth by altering thresholds used:
%       -Threshold 1: mean intensity +/- 11.72% of the max intensity (Km)
%       -Threshold 2: (mean intensity - Km) to max
%       -Threshold 3: mean to max
%
% Inputs:
%   -Im = single-channel image
%
% Outputs:
%   -StatArray = a 1x27 array, of type double, arranged:
%       -[0n, 1n, 2n, 3n, 4n, 5n, 6n, 7n, 8n], where
%           -0n = fraction of pixels with no neighbours
%           -1n...8n = fraction of pixels with 1...8 neighbouring pixels 
%            above threshold
%       -This array is then repeated two more times for the additional
%       thresholding methods
%
% Examples:
%   -0 = below threshold   -X = pixel under analysis
%   -1 = above threshold
%
%                000  100  110  111  111  111  111  111  111
%                0X0  0X0  0X0  0X0  1X0  1X1  1X1  1X1  1X1
%                000  000  000  000  000  000  100  110  111 
%
% neighbours:     0    1    2    3    4    5    6    7    8
% StatArray = [0, 0.0278, 0.0556, 0.0833, 0.1111, 0.1389, 0.1667, 0.1944,
% 0.2222]
%
% Recommended usage:
% 1) Capture images of your samples with no or minimal alteration to your
% acquisition settings.
% 2) Crop images to single cells.
% 3) Import images into Matlab and pass to this function as a 2D array
% 4) Manually or via a script, analyze each image with this function
% 5) Use PCA or another clustering algorithm to identify clusters in the
% data.
%
%
%  --Bryan Heit--  
% V1: July 23, 2009 - first working copy
% V2: Sept 17, 2010  - changed Km to 0.5 STD


%% Check Input
if nargin ~= 1
    error('Incorrect number of input variables.  Input should be a single image')
end

if ndims(Im) ~= 2
    error ('Input is not a 2-dimentional image.');
end


%% Calculate Threshold Cutoffs

ImMean = mean(Im(:)); %mean image intensity
Kmax = max(Im(:)); %max image intensity
%Km = 0.1172*Kmax; %V1
Km = std(double(Im))/2; %V2, more versitile
StatArray = zeros(1,27);

%% Analyze image
for aa=1:3
    tArray = zeros (1,10);
    ImT = Im;
    % Generate Threshold
    if aa == 1
        ImT(ImT < (ImMean-Km)) = 0;
        ImT(ImT > (ImMean+Km)) = 0;
        ImT(ImT ~= 0) = 1;
    elseif aa == 2
        ImT(ImT < (ImMean-Km)) = 0;
        ImT(ImT ~= 0) = 1;
    else
        ImT(ImT < (ImMean)) = 0;
        ImT(ImT ~= 0) = 1;
    end
    for ii=2:max(size(Im(1,:)))-1
      for jj=2:max(size(Im(:,1)))-1
          if ImT(jj,ii)~=0
              Neighbours = 0; %reset neighbour counter
              tArray(10) = tArray(10)+1; %score pixel as above background
              if ImT(jj-1,ii)~=0
                  Neighbours = Neighbours + 1;
              end
              if ImT(jj-1,ii-1)~=0
                  Neighbours = Neighbours + 1;
              end
              if ImT(jj-1,ii+1)~=0
                  Neighbours = Neighbours + 1;
              end
              if ImT(jj,ii+1)~=0
                  Neighbours = Neighbours + 1;
              end
              if ImT(jj,ii-1)~=0
                  Neighbours = Neighbours + 1;
              end
              if ImT(jj+1,ii)~=0
                  Neighbours = Neighbours + 1;
              end
              if ImT(jj+1,ii-1)~=0
                  Neighbours = Neighbours + 1;
              end
              if ImT(jj+1,ii+1)~=0
                  Neighbours = Neighbours + 1;
              end
              Neighbours = Neighbours +1; % sets proper index point
              tArray(Neighbours) = tArray(Neighbours)+1; %score number of neighbours
          end % end if loop
      end %end internal (j) for loop
    end % end external (i) for loop

        % Convert stat array to fraction, and output data
    if tArray(10)>0 %this is required to avoid a computation error if a blank image is inputted
        tArray = tArray(1:9)./tArray(10);
    end

    rPos = ((aa-1)*9)+1;
    StatArray(rPos:rPos+8) = tArray;
    
end %end aa loop

end % end function