CODEWITHSUNDEEP
imagematlabimage-processing
wrek
wrek

Reputation: 309

Object ellipticity based on the moment invariants - Matlab

Following a question about Hu moments and a question regarding object Ellipse Variance I used the code from this link to try to calculate ellipticity based on the moment invariants like in the image below (Based on A Beginner’s Guide to Image Shape Feature Extraction Techniques). I the code I used (only checking the first 2 objects) I do not get a value close to 1 (which should indicate an Ellipticity) when I have an ellipse object. Can someone please explain what I did wrong?

enter image description here

code:

clc;
clear;
close all;

% Image processing
I= rgb2gray(imread('https://az877327.vo.msecnd.net/~/media/images/products/2005/other/mipix%20ellipses%20full%20frame.jpg?v=1&h=550&w=800&crop=1'));
bw = imbinarize(I);
bw = imfill(bw,'holes');
bw = bwareaopen(bw, 100);
imshow(bw)
hold on;
[B,L] = bwboundaries(bw,'noholes');
stats = regionprops(L,'Centroid','Image');

%Loop - data acquisition
for i = 1 : 2 %numel(stats)
    b = B{i};
    c = stats(i).Centroid;
    f = stats(i).Image;
    y = b(:,1);
    x = b(:,2);
    plot( b(:,2),b(:,1),'Color','red','linewidth',1);
    text(c(1),c(2),num2str(i),'Color','red'); 

    %based on https://raw.githubusercontent.com/dipum/dipum-toolbox/master/dipum/invmoments.m
    [M,N] = size(f);
    [x,y] = meshgrid(1:N,1:M);

    % Turn x,y, and F into column vectors to make the summations a bit
    % easier to compute in the following.
    x = x(:);
    y = y(:);
    f = f(:);

    % DIPUM3E equation (13-27)
    m.m00 = sum(f);
    % Protect against divide-by-zero warnings.
    if (m.m00 == 0)
       m.m00 = eps;
    end
    % The other central moments:  
    m.m10 = sum(x .* f);
    m.m01 = sum(y .* f);
    m.m11 = sum(x .* y .* f);
    m.m20 = sum(x.^2 .* f);
    m.m02 = sum(y.^2 .* f);
    m.m30 = sum(x.^3 .* f);
    m.m03 = sum(y.^3 .* f);
    m.m12 = sum(x .* y.^2 .* f);
    m.m21 = sum(x.^2 .* y .* f);

    % DIPUM3E equations (13-28) through (13-30).
    xbar = m.m10 / m.m00;
    ybar = m.m01 / m.m00;

    e.eta11 = (m.m11 - ybar*m.m10) / m.m00^2;
    e.eta20 = (m.m20 - xbar*m.m10) / m.m00^2;
    e.eta02 = (m.m02 - ybar*m.m01) / m.m00^2;
    e.eta30 = (m.m30 - 3 * xbar * m.m20 + 2 * xbar^2 * m.m10) / ...
                m.m00^2.5;
    e.eta03 = (m.m03 - 3 * ybar * m.m02 + 2 * ybar^2 * m.m01) / ...
                m.m00^2.5;
    e.eta21 = (m.m21 - 2 * xbar * m.m11 - ybar * m.m20 + ...
               2 * xbar^2 * m.m01) / m.m00^2.5;
    e.eta12 = (m.m12 - 2 * ybar * m.m11 - xbar * m.m02 + ...
               2 * ybar^2 * m.m10) / m.m00^2.5;

    % DIPUM3E Table 13.8.
    phi(1) = e.eta20 + e.eta02;
    phi(2) = (e.eta20 - e.eta02)^2 + 4*e.eta11^2;
    phi(3) = (e.eta30 - 3*e.eta12)^2 + (3*e.eta21 - e.eta03)^2;
    phi(4) = (e.eta30 + e.eta12)^2 + (e.eta21 + e.eta03)^2;
    phi(5) = (e.eta30 - 3*e.eta12) * (e.eta30 + e.eta12) * ...
             ( (e.eta30 + e.eta12)^2 - 3*(e.eta21 + e.eta03)^2 ) + ...
             (3*e.eta21 - e.eta03) * (e.eta21 + e.eta03) * ...
             ( 3*(e.eta30 + e.eta12)^2 - (e.eta21 + e.eta03)^2 );
    phi(6) = (e.eta20 - e.eta02) * ( (e.eta30 + e.eta12)^2 - ...
                                     (e.eta21 + e.eta03)^2 ) + ...
             4 * e.eta11 * (e.eta30 + e.eta12) * (e.eta21 + e.eta03);
    phi(7) = (3*e.eta21 - e.eta03) * (e.eta30 + e.eta12) * ...
             ( (e.eta30 + e.eta12)^2 - 3*(e.eta21 + e.eta03)^2 ) + ...
             (3*e.eta12 - e.eta30) * (e.eta21 + e.eta03) * ...
             ( 3*(e.eta30 + e.eta12)^2 - (e.eta21 + e.eta03)^2 );


    disp('Object number:');
    disp(i);
    Is = (m.m20 * m.m02-m.m11^2)/m.m00^4;
    Em1= 16*pi^2*Is
    Em2= (16*pi^2*Is)^-1
    min(Em1, Em2)
end

Upvotes: 0

Views: 332

Answers (1)

beaker
beaker

Reputation: 16811

As I said in the comments, you need to use the central moments rather than raw moments as your code is doing. See Measuring Shape: Ellipticity, Rectangularity, and Triangularity

Removing the unnecessary calculations for eta and phi, the last part of your code should be:

% Calculate Raw Moments:  
m.m10 = sum(x .* f);
m.m01 = sum(y .* f);
m.m11 = sum(x .* y .* f);
m.m20 = sum(x.^2 .* f);
m.m02 = sum(y.^2 .* f);
m.m30 = sum(x.^3 .* f);
m.m03 = sum(y.^3 .* f);
m.m12 = sum(x .* y.^2 .* f);
m.m21 = sum(x.^2 .* y .* f);

% DIPUM3E equations (13-28) through (13-30).
xbar = m.m10 / m.m00;
ybar = m.m01 / m.m00;

% Calculate Central Moments:
m.mu11 = m.m11 - xbar*m.m01;
m.mu20 = m.m20 - xbar*m.m10;
m.mu02 = m.m02 - ybar*m.m01;

disp('Object number:');
disp(i);

Is = (m.mu20 * m.mu02 - m.mu11^2)/m.m00^4;

Em1= 16*pi^2*Is
Em2= (16*pi^2*Is)^-1
min(Em1, Em2)

With these changes, the results for the first two objects are:

Object number:
 1
Em1 =  1.000005232603003
Em2 =    9.999947674243773e-01
ans =    9.999947674243773e-01
Object number:
 2
Em1 =    9.999818710527637e-01
Em2 =  1.000018129275901
ans =    9.999818710527637e-01

Upvotes: 2

Related Questions