findind(A,B)

function [IB, IA] = findind(A,B)
% FINDIND - find indices of matching elements between two matrices
%
%   IB = FINDIND(A,B) returns a cell array IB which has the same size as A.
%   Each cell of IB contains a row vector with the indices of those elements
%   of B that are equal to the corresponding element of A.  In other words,
%   for each k, IB{k} equals "find(B == A(k))".
%
%   [IB, IA] = FINDIND(A,B) also return another cell array IA holding the
%   indices of those elements of A that are equal to the corresponding
%   element of B.  In other words, for each k, IA{k} equals "find(A == B(k))".
%
%   A and B may have repeated elements. If an element of A or B is not
%   present in B or A, the corresponding cell will be empty. A and B can be
%   cell array of strings.  
%
%   Example:
%      A = [1 4 2 ; 3 4 5] ;
%      B = [1 2 1 1 4 2 2 2 4 5 5] ;
%      x = findind(A,B)
%      % -> x = {[ 1 3 4 ], [ 5 9 ], [ 2 6 7 8 ] 
%      %         []       , [ 5 9 ], [10 11]     }
%
%   See also FIND, ISMEMBER, UNIQUE

% tested for Matlab R2008a (requires bsxfun)
% version 2.0 (jan 2009)
% (c) Jos van der Geest
% email: [e-mail]

% History:
% 1.0, jan 2009 - creation, inspired by Loren's blog
% 2.0, jan 2009 - 1) fixed bug when none of the elements matched
%                 2) fixed some lay-out issues, improved help
%                 3) added second return argument IA

error(nargchk(2,2,nargin)) ;

if ~((iscellstr(A) && iscellstr(B)) || (isnumeric(A) && isnumeric(B))) ,
    error('A and B should be of the same class (numeric or cell array of strings).') ;
end

% IA and IB refer to indices into A and into B resp. 

% first handle empty arrays
if numel(A)==0 || numel(B) == 0,    
    IB = cell(size(A)) ;
    if nargout > 1,
        IA = cell(size(B)) ;
    end
else    
    if iscellstr(A),        
        % Both A and B are cell strings. 
        % Convert to numbers using the UNION/ISMEMBER trick
        AB = union(A,B) ;
        [dum,A] = ismember(A,AB) ;
        [dum,B] = ismember(B,AB) ;
    end
    
    try
        M = bsxfun(@eq, A(:).',B(:)) ;
    catch
        % most probably the error is realted to an older matlab version,
        % which does not has BSXFUN. So, try repmat instead
        warning([mfilename ':BsxfunNotFound'],'BSXFUN error, using REPMAT instead.') ;
        try
            M = repmat(A(:).',numel(B),1) == repmat(B(:),1,numel(A)) ;
        catch
            % This error should not occur, except for rare installs
            error('Cannot apply BSXFUN or REPMAT.') ;
        end
    end
    
    % Each column of M corresponds to an element of A, and each row
    % corresponds to an element of B. 
    
    if any(M(:)),
        % if elements A(iA) and B(iB) match, M(iB,iA) will be true        
        [iB,iA] = find(M) ;
        % count the N(k) elements of B matching A(k),
        N = histc(iA,1:numel(A)) ;
        % the first N(1) indices iB correspond to A(1), etc. Use these numbers
        % to split the indices iB into cells. There will be numel(A) cells.
        IB = mat2cell(iB.',1,N) ;
        % and reshape into the same size as A
        IB = reshape(IB,size(A)) ;
        
        % This can be done in one step which is less clear
        % IB = reshape(mat2cell(iB.',1,histc(iA,1:numel(A))),size(A)) ; 
    
        % set empty cells to 0-by-0 empty, and not 1-by-0 to avoid confusion
        q = cellfun('isempty',IB) ;
        IB(q) = {[]} ;
        
        if nargout > 1,
            % and apply the same for the indices into A
            N = histc(iB,1:numel(B)) ;
            IA = mat2cell(iA.',1,N) ;
            IA = reshape(IA,size(B)) ;
            IA(cellfun('isempty',IA)) = {[]} ;
        end
        
    else % any(M)
        % none of elements matched so all cells will be empty
        IB = cell(size(A)) ;
        if nargout > 1,
            IA = cell(size(B)) ;
        end
    end
end