function [roi_idx, match_qlt] = findsubs(A, B, max_match_qlt)
% suchen einer Submatrix B in Hauptmatrix A
%
%  input   : A             - Hauptmatrix
%            B             - Submatrix
%            max_match_qlt - Minimale Überlappungsqualität
%  output  : roi_idx       - Indizien zu Matrix A
%            match_qlt     - Überlappungsqualität
%  example :
%         A     = magic(12);
%         B     = [71 70;86 87];
%         [idx,qlt] = findsubs(A,B,100)
%         isequal(A(idx{1}),B)

if nargin == 0 || nargin == 1 
  return;
end

a_size = size(A);
b_size = size(B);

if nargin == 2
  max_match_qlt = 100;
end


pos_y = zeros(1,(a_size(1) - b_size(1)+1)*(a_size(2) - b_size(2)+1)); 
pos_x = zeros(1,(a_size(1) - b_size(1)+1)*(a_size(2) - b_size(2)+1));   
match = zeros(1,(a_size(1) - b_size(1)+1)*(a_size(2) - b_size(2)+1));  

for y    = 1 : a_size(1) - b_size(1)+1;   
   for x = 1 : a_size(2) - b_size(2)+1;
      pos_y(sub2ind([a_size(1) - b_size(1)+1,a_size(2) - b_size(2)+1],y,x))= y;
      pos_x(sub2ind([a_size(1) - b_size(1)+1,a_size(2) - b_size(2)+1],y,x))= x;
      mat_match = bsxfun(@eq, A(y:y+b_size(1)-1, x:x + b_size(2)-1), B);
      match(sub2ind([a_size(1) - b_size(1)+1,a_size(2) - b_size(2)+1],y,x))= sum(mat_match(:));
   end   
end

match_qlt = (match/(b_size(1) * b_size(2)))*100;
pos_qlt = find(match_qlt >= max_match_qlt);
roi_idx = cell(1,length(pos_qlt));

for k = 1:length(roi_idx)
   pos   = pos_qlt(k);
   [x,y] = meshgrid(pos_y(pos):pos_y(pos)+b_size(1)-1,pos_x(pos):pos_x(pos) + b_size(2)-1);
   roi_idx{k}= sub2ind(a_size,x,y)';
end

match_qlt = match_qlt(pos_qlt);
[match_qlt,sort_idx] = sort(match_qlt,'descend');
roi_idx = roi_idx(sort_idx);