Shape matching tactics for shapes and complex images

Hi everyone.

First up, I must declare that I am a student of neuroscience who studied biology at the bachelor's level and therefore my mathematical knowledge is pretty basic. I just started using MATLAB in October and therefore my programming fluency is pretty basic too.

So here goes...

1) The need

Our lab has this set of image pairs. They need to be matched for shape similarity and a metric needs to be generated, i.e. shape A is x% similar to shape B.

2) My tactics till now -

a) MSE - I extracted the edges of the images using edge and it gave me a matrix with the dimensions of the image and appropriate places filled with 1s and the rest 0s. I then subtracted IM1 from IM2, calculated the MSE. Immediately, we see that this is not scale-invariant.

b) Pixel-Comparison - A little better than MSE but still is scale-invariant.

c) Distance search - I made a 'window' of 5px x 5px run over the image matrix. Wherever there was a 1 in this tiny window, I calculated its distance to the top-right hand corner point on the window. Thus I compared both images and said wherever the distance is same, give me a 1 and rest 0. Then I summed up the 1s, divided by the length of the matrix and got my %age. Much better but still would foul up many times when I scaled down the image considerably keeping the canvas size same.

d) Then I tried the same window distance search algorithm but I calculated distances from each point to each point on the vertical and horizontal sides of the window, asked for each minimum and generted an array of the minimums from each image and took their ratio. I said, if the ratio lies between 0.7 and 1.3, give me a 1 else 0, summed up the 1s and calculated %age. This worked for a circle and an ellipse, two different ellipses etc but surprisingly didn't work for a perfect square and a square with rounded edges of the same size...it gave me a similarity percentage of 0 because all the ratios are either 1.526 or greater.

I'm now at the end of my tether! No idea how to proceed except one. But I don't know how to implement it.

I want to calculate the first derivative to each point on the outline of the image. If the derivate at two corresponding pixels in two of the images are the same, then the shape there is the same, otherwise it isn't. This would be truly scale-invariant. But since the boundary matrix is logical (1s and 0s), I don't know how to calculate this first derivative test. I tried asking for the tangents & slopes but of course it is a dumb idea and it didn't work.

If anyone could provide me any help, guidance, pointers, tips; I'd be indeed very grateful.

Attached below is my commented code from the window-search paradigm.

clear all

% This is the minimum distance ratio based shape matcher. Please note that
% this version requires that the images be BMP in format and have
% dimensions divisible by 100, e.g. 200x400 or 500x500 etc. In the next
% versions (if this is successful), I will generalise it to accept jpg and
% png and svg images too, and it will work even if the dimensions are not
% divisible by 100. I have the echo on so that you can see the computations
% happening.

% Abhilash Dwarakanath, GSN-LMU (Munich), Berlin, 25.1.2011.

echo on

% Allowing user to dynamically choose the images he/she wants to compare

shape1_name = input('Please enter the name of the 1st file without the .bmp within single quotes : ');
shape2_name = input('Please enter the name of the 2nd file without the .bmp within single quotes : ');

%Reading the images into a matrix of values with 1s on the boundary and 0s
%elsewhere

shape1 = imread(shape1_name, 'bmp');
shape2 = imread(shape2_name, 'bmp');

% Checking if the image is an 8-bit coloured or 1-bit BW bmp and converting
% to 1-bit BW values of image is coloured.

l1 = islogical(shape1);
l2 = islogical(shape2);

if l1 ~= 1 && l2 ~= 1

      shape1_rgb = im2double(shape1);
      shape2_rgb = im2double(shape2);
      shape1_int = logical(rgb2gray(shape1_rgb));
      shape2_int = logical(rgb2gray(shape2_rgb));
else

      shape1_int = shape1;
      shape2_int = shape2;
end

% Detecting the edge using the Canny method.

shape1_bound = 0+(edge(shape1_int, 'canny'));
shape2_bound = 0+(edge(shape2_int, 'canny'));

% Taking all matrices to nearest 100

shape1_bound = shape1_bound(1:400,1:400);

shape2_bound = shape2_bound(1:400,1:400);

% Calculating similarity - This similarity calculation is based on a ratio
% based test. If one were to just calculate similar pixels, it would result
% in errors if the shapes were absolutely the same, but in different
% corners of the canvas.

% Therefore, I use a square window which is 1/100th the dimensions of the
% canvas and move it over the canvas. At every step, it will sample all
% those indices in which the value of boundary matrix is 1. And then it
% will calculate the distance from that point to each point on the two
% sides of the window, i.e. to the vertical and the horizontal. Any zero
% elements are squeezed out by parsing, because if the distance is 0, then
% the point is a tangent to the shape and therefore is of no consequence in
% estimating which is the minimum distance. Then, the minimum distance is
% added to an array.

dist1 = [];
dist2 = [];

[r,c] = size(shape1_bound);

for i = 1:r/100:r-((r/100)-1) % Running row index from 1 to 1/100th of row

      for j = 1:c/100:c-((c/100)-1) % Same for column index.
          bridge1 = shape1_bound(i:i+((r/100)-1),j:j+((c/100)-1));
          % This bridge is now my window which has a dimension 1/100th of R x
          % 1/100th of C. The loop is setup in such a way that at every
          % iteration, it will sample the next r/100 and c/100 points.
          [rows,cols] = size(bridge1);
          % Now the size of this window is referenced for looping
          for k = 1:rows
              for l = 1:cols
                  if bridge1(k,l) == 1
                      % In that window, wherever in the shape matrix there is
                      % a 1, i.e. an edge is detected, I calculate the
                      % distance from that point to all the points on the
                      % edge of the window.
                      for k1 = i:i+4
                          d1x(k1) = sqrt((k1 - k)^2 + (1 - l)^2); %#ok<*SAGROW>
                      end
                      for l1 = j:j+4
                          d1y(l1) = sqrt((1 - k)^2 + (l1 - k)^2);
                      end
                  % I now squeeze out all the zeros to eliminate getting 0 
                  % when I ask for the min    
                      d1cx = 0;
                      for m = 1:length(d1x)
                          if d1x(m) == 0
                              d1cx = d1cx+1;
                          end
                      end
                      d1xx = d1x(d1cx+1:length(d1x));
                      d1cy = 0;
                      for n = 1:length(d1y)
                          if d1y(n) == 0
                              d1cy = d1cy+1;
                          end
                      end
                      d1yy = d1y(d1cy+1:length(d1y));
                      % I now concatenate the x and y distances and ask for
                      % the min of them
                      d1 = [d1xx,d1yy];
                      d1min = min(d1);
                      dist1 = [dist1 d1min]; %#ok<*AGROW>
                  end
              end
          end
      end
end

% The above steps are repeated for the 2nd image

for i = 1:r/100:r-((r/100)-1) % Running row index from 1 to 1/100th of row

      for j = 1:c/100:c-((c/100)-1)
          bridge2 = shape2_bound(i:i+((r/100)-1),j:j+((c/100)-1));
          [rows,cols] = size(bridge2);
          for k = 1:rows
              for l = 1:cols
                  if bridge2(k,l) == 1
                      for k1 = i:i+4
                          d2x(k1) = sqrt((k1 - k)^2 + (1 - l)^2); %#ok<*SAGROW>
                      end
                      for l1 = j:j+4
                          d2y(l1) = sqrt((1 - k)^2 + (l1 - k)^2);
                      end
                      d2cx = 0;
                      for m = 1:length(d2x)
                          if d2x(m) == 0
                              d2cx = d2cx+1;
                          end
                      end
                      d2xx = d2x(d2cx+1:length(d2x));
                      d2cy = 0;
                      for n = 1:length(d2y)
                          if d2y(n) == 0
                              d2cy = d2cy+1;
                          end
                      end
                      d2yy = d2y(d2cy+1:length(d2y));
                      d2 = [d2xx,d2yy];
                      d2min = min(d2);
                      dist2 = [dist2 d2min]; %#ok<*AGROW>
                  end
              end
          end
      end
end

% There is a chance that one of the distance matrices might be smaller than
% the other. If that happens, we have to make them both the same size for
% the matrix divison to be possible

if length(dist1) < length(dist2)

      dist2 = dist2(1:length(dist1));
else

      dist1 = dist1(1:length(dist2));
end

% Now, we have two arrays, each containin the minimum distances from any
% given point on the edge to the nearest perpendicular point. I now
% calculate the ratio of these distances. The logic I use is this. Suppose
% I have two unit circles. The minimum distance from a point on this
% circle, when divided by the minimum distance from a point on the other
% unit circle to a similar point on the edge of our standard square will be
% the same, and their ratios will be 1. However, if we now take an ellipse,
% the ratio will differ as the ratio of the semi-minor axis to the
% semi-major axis. Now, since a circle is an ellipse with its semi-minor
% axis = semi-major axis, this ratio will be constant for every tangential
% projection, i.e. for a corresponding pixel, the ratios of orientation
% will be constant! This will then give us a good measure of similarity or
% dissimilarity.

% Therefore, I ask the program for ratio values that lie between certain
% numbers and assign them different weights. Then I add them up, sum and
% calculate the %age. NaNs are counted as 1.

ratio = dist1./dist2;

ratio = abs(ratio - 1);

for i = 1:length(ratio)

      if ratio(i) >= 0.9 && ratio(i) <= 1.1
          hit(i) = 1;
      elseif ratio(i) >= 0.5 && ratio(i) <= 0.6
          hit(i) = 0.55;
      elseif ratio(i) >= 0.6 && ratio(i) <= 0.7
          hit(i) = 0.65;
      elseif ratio(i) >= 0.7 && ratio(i) <= 0.8
          hit(i) = 0.75;
      elseif ratio(i) >= 0.8 && ratio(i) <= 0.9
          hit(i) = 0.85;
       elseif ratio(i) >= 1.1 && ratio(i) <= 1.2
          hit(i) = 0.85;
      elseif ratio(i) >= 1.2 && ratio(i) <= 1.3
          hit(i) = 0.75;
      elseif ratio(i) >= 1.3 && ratio(i) <= 1.4
          hit(i) = 0.65;
      elseif ratio(i) >= 1.4 && ratio(i) <= 1.5
          hit(i) = 0.55;
      else
          hit(i) = 0;
      end
end

nans = isnan(ratio);

n = 0;

for i = 1:length(nans)

      if nans(i) == 1
          n = n+1;
      else
          n = n+0;
      end
end

total = sum(hit) + n;

percent_similar = (total/length(ratio))*100 %#ok<NOPTS>

echo off

% I have tried here to compute a p measure (not to be confused with the p
% measure of an ANOVA) that dynamically is updated as and when similarity
% is calculated. Somehow it doesn't work completely. After a number of
% iterations, it gives me indexing error even though I haven't set the
% number of elements of P at all. 

% In the future, I plan to turn this into a Bayesian probability
% calculation for similarity. For that, our parameters must be discussed
% and defined.

% p(1) = 0;
% 
% for i = 1:length(ratio)
%     
%     if ratio(i) >= 0.6 && ratio(i) <= 1.4
%         
%         p(i+1) = p(i) + abs(1 - ratio(i));
%         
%     elseif ratio(i) < 0.6
%         
%         p(i+1) = p(i) - (1 - ratio(i));
%         
%     elseif ratio(i) > 1.4
%         
%         p(i+1) = p(i) + (1 - ratio(i));
%         
%     end
%     
% end
% 
% similarity = p

I can generate a contour with the contour plot, but can I incorporate that into my program?

[Doug Eastman]

I'm not an image processing expert, but if your original image is simply a white object on a black background could you use REGIONPROPS to extract various statistics such as eccentricity, perimeter/area, etc.? That could help neglect translations and scaling, but depending on the shapes it may not give you enough granularity.

[Abhilash]

Hi,

Thanks! No that wouldn't work because the images are pretty complex. For instance, a banana is being compared to an industrial crane whose neck is oriented in the same direction as the banana. Etc.

Answer 1

[Matt Wolinsky]

This is a difficult problem!

However the Computer Vision community has developed some pretty good solutions.

I'm not a Computer Vision expert, but some suggestions:

1) Any sort of local similarity/difference metric, e.g. norm(A(:)-B(:)), will be sensitive to noise/scale/rotation/illumination, etc.

So my first suggestion would be to use a robust similarity metric. A simple and robust one would be the "Pyramid Match Kernel":

Grauman and Darrell (2005)

The pyramid match kernel: Discriminative classification with sets of image features

http://scholar.google.com/scholar?cluster=17504420260500492874&hl=en&as_sdt=0,44

This is technically a metric on point sets (or histograms thereof) . . . but can work pretty well on images directly if you take the image as the "level 0 histogram" (see ref above).

2) If that doesn't work by itself, you can use it in its original framework: As a distance metric between "feature vectors" in a high-dimensional but sparse "feature space" derived from the image. There is a large literature on feature detectors, and I haven't ever used them . . . but this looks like a good bet:

OpenSURF (including Image Warp)

by Dirk-Jan Kroon (July 2010) SURF (Speeded Up Robust Features) image feature point detection / matching, as in SIFT

http://www.mathworks.com/matlabcentral/fileexchange/28300-opensurf-including-image-warp

The presentation here gives a brief overview of these approaches:

Scalable Image Recognition and Retrieval

http://www.krellinst.org/doecsgf/conf/2007/pres/grauman.shtml

Hope that helps!

--Matt Wolinsky

[Abhilash]

Thanks Matt!

Answer 2

[Abhilash]

I'd really appreciate any help on this matter.

Thanks