B_matrix(C, N1, N2, p)

function B = B_matrix(C, N1, N2, p)
% B_matrix: B matrix for Self-Organizing Map after Cluster Reinforcement Phase (Manukyan et al, 2012)
% B = B_matrix(C)
% B = B_matrix(C, p)
%
% The B matrix is used for visualization of sharpened boundaries 
% after CR phase.
% It essentially stores 3 kinds of distances (horizontal, vertical and
% diagonal) between adjacent neurons. See reference below for more details.
%
% TO VISUALIZE THE B-MATRIX USE ONE OF THE FOLLOWING FUNCTIONS:
% a) Plot_B
% b) DrawBLines
% c) B_GUI
%
% PLEASE REFERENCE THE FOLLOWING IF YOU USE THIS CODE:
% N. Manukyan, M.J. Eppstein, D.M. Rizzo,"Data-driven cluster reinforcement 
% and visualization in sparsely-matched self-organizing maps,"
% Neural Networks and learning Systems, IEEE Transactions on, vol23, no 5, pp 846-852,may 2012.      
%
%
% INPUTS:
%         C:       N x M matrix of N SOM-trained neurons, each with M features
%                  (i.e., linear form of trained SOM with M component planes)
%         N1:      scalar representing the first dimension of the rectangular SOM map
%         N2:      scalar representing the second dimension of the rectangular SOM map
%               
%         p:       optional 1 x M vector of optional weighting for the M component planes (defaults to ones(1,M))
%
% OUTPUT:
%         B:       2*N1 x 2*N2 B-matrix (Manukyan et al, 2012), where N2 is
%                  the second dimension of the rectangular SOM map
%
%
% Authors and Contact information:
% Narine Manukyan 
% University e-mail: [e-mail]
% Lifetime email: [e-mail]
%
% Margaret J. Eppstein
% e-mail: [e-mail]
%
% Department of Computer Science
% University of Vermont
% Burlington Vermont, 05401
% 
% Code posted to Matlab File Exchange: 15-March-2012

if nargin < 1
    error('Requires at least one input argument.');
end

N = size(C,1); % number of neurons on the grid

M = size(C,2); % number of features

% Check if the features are weighted
if nargin ==4 && ~all(p==1)
    C = p(ones(N,1),:).*C; %apply the optional feature weights to each neuron 
end

C = reshape(C,N1,N2,M);% Convert the 2D matrix to 3D; 
% (first 2 dimensions are for the rectangular SOM, third dimension is the features)

% Define South and East neighbor indexes
South = [2:N1 1]; % allow for toroidal wrap (but still works if not toroidal)
East = [2:N2 1];

% Define odd element indexes on both dimensions of the map
oddElements1 = 1:2:2*N1-1;
oddElements2 = 1:2:2*N2-1;

% Define even element indexes on both dimensions of the map
evenElements1 = 2:2:2*N1;
evenElements2 = 2:2:2*N2;

% Initialize the B matrix
B = zeros(2*N1,2*N2);

% Calculate horizontal elements (distance between horizontal elements in C)
B(oddElements1,evenElements2) = eucdist(C,C(:,East,:));

% Calculate vertical elements (distance between vertical elements in C)
B(evenElements1,oddElements2) = eucdist(C,C(South,:,:));

% Calculate diagonal elements (mean of the distances between elements on two diagonals in C)
B(evenElements1,evenElements2) = (eucdist(C,C(South,East,:)) + eucdist(C(South,:,:),C(:,East,:)))./2;

end