How to locate the the position of a 3D subimage in larger image in Matlab?

I have a 3D subimage saved as X.mat , which is a subimage of the 3D image Y.mat. How can I find the position of this subimage in the image Y?

Can anyone help me?

Thanks,

Answer 1

[David Young]

Assuming you have used load to read the mat-files and the images are stored as X and Y, you can look for an exact match by checking each subimage of Y.

found = false;
for r = 1:size(Y,1)-size(X,1)+1
  for c = 1:size(Y,2)-size(X,2)+1
    for s = 1:size(Y,3)-size(X,3)+1
      if isequal(Y(r:r+size(X,1)-1, c:c+size(X,2)-1, s:s+size(X,3)-1), X)
        found = true;
        break
      end
    end
    if found
      break
    end
  end
  if found
    break
  end
end
if found
  fprintf('Match found at: %d %d %d\n', r, c, s);
else
  fprintf('No match found\n');
end

This reports the coordinates of the low-index corner of the match.

If the subimage is not an exact subimage, you need to use some other technique, such as normalised cross-correlation.

[margol]

Dear David,

Thank you very much for your help and time. The algorithm you suggested found just one coordinate, but one friend suggested the following code which essentially scans through the large array A and compare each element against the first one in the small array B. If an equal is found, a part of A, of which the size is same as of B, is picked out and compared against B.

clear;clc
% reproduce your scenario
A = randi(100, [30, 20, 10]);
B = A(20:30, 1:18, 4:end);
% counter-case verification
% B(end)=200;
% speed up
lenA = numel(A);
[sa1,sa2,sa3] = size(A);
[sb1,sb2,sb3] = size(B);
% a cumbersome method
eqflag = 0;
counterA = 1;
while (counterA <= lenA)
    if A(counterA) == B(1)
        [subA1,subA2,subA3] = ind2sub([sa1,sa2,sa3],counterA);
        if ( (subA1+sb1-1)<=sa1 ) && ( (subA2+sb2-1)<=sa2 ) ...
                && ( (subA3+sb3-1)<=sa3 ) && isequal( B, ...
            A(subA1+(1:sb1)-1,subA2+(1:sb2)-1,subA3+(1:sb3)-1) )
            eqflag = 1;
            break;
        end
    end
    counterA = counterA + 1;
end
if eqflag
    fprintf('found matching starting at A(%d, %d, %d).\n', ...
        subA1, subA2, subA3);
    fprintf('matching region A(%d:%d, %d:%d, %d:%d).\n', ...
        subA1, subA1+sb1-1, subA2, subA2+sb2-1, subA3, subA3+sb3-1);
else
    fprintf('no matching found.\n');
end
clearvars sa* lenA counterA
% --------------
% a parallel way
[sa1,sa2,sa3] = size(A);
match_first = find(A==B(1));
[m1,m2,m3] = ind2sub([sa1,sa2,sa3],match_first);
region_first_ind = intersect( intersect(find(m1+sb1-1<=sa1), ...
    find(m2+sb2-1<=sa2)),find(m3+sb3-1<=sa3)); % array size issue
region_first = num2cell( [m1(region_first_ind),m2(region_first_ind),...
    m3(region_first_ind)], 2);
region = cellfun(@(v) [v;v+[sb1,sb2,sb3]-1], region_first, ...
    'UniformOutput', false);
region_match = cellfun(@(v) isequal(A(v(1):v(2), v(3):v(4), v(5):v(6)),...
    B), region, 'UniformOutput', false);
match = cell2mat(region([region_match{:}]));
if ~isempty(match)
    fprintf('found matching starting at A(%d, %d, %d).\n', ...
        match(1), match(3), match(5));
    fprintf('matching region A(%d:%d, %d:%d, %d:%d).\n', ...
        match(1), match(2), match(3), match(4), match(5), match(6));
else
    fprintf('no matching found.\n');
end

Margol

[David Young]

I'm not sure you've understood my answer. What does it not do that you want to do?

[Image Analyst]

It should find only one coordinate if it occurs just one time in the larger volume. (I didn't look over the code in detail but assumes it finds the upper left corner or something like that.) If it occurs twice it will find two . If you wanted a listing of all coordinates in the rectangular block at the one "found" location then you can to that by adding the width, length, and height to the one found coordinate (upper left corner). That will get you a whole "block" of coordinates.

Answer 2

[Ahmet Cecen]

Depending on the size and type of data, it might be preferable to use a convolution or an erosion. I can't give a specific code snippet without further information.

[Alex Taylor]

I just wanted to chain off of Ahmet that this is a 3-D template matching problem. David Young's solution is very helpful, and will work fine if you are truly looking for one and only one EXACT match in the larger volume. More typically when doing template matching, either the template or the larger image are not exact matches of one another (there are noise, illumination differences, etc.). In these cases, people usually turn to algorithms like normalized cross correlation and determine the position in the correlation matrix with the highest peak correlation. This method will work even when there isn't an exact match.

You might want to look at

normxcorr2

in IPT as an introduction to the 2-D template matching workflow. You could try to implement your own 3-D extension to normxcorr2 using

convn

or

fftn

(depending on your template image size either the spatial domain or the frequency domain will be faster).