matlab - Computing object statistics from the second central moments

Question

I'm currently working on writing a version of the MATLAB RegionProps function for GNU Octave. I have most of it implemented, but I'm still struggling with the implementation of a few parts. I had previously asked about the second central moments of a region.

This was helpful theoretically, but I'm having trouble actually implementing the suggestions. I get results wildly different from MATLAB's (or common sense for that matter) and really don't understand why.

Consider this test image:

We can see it slants at 45 degrees from the X axis, with minor and major axes of 30 and 100 respectively.

Running it through MATLAB's RegionProps function confirms this:

MajorAxisLength: 101.3362
MinorAxisLength: 32.2961
Eccentricity: 0.9479
Orientation: -44.9480

Meanwhile, I don't even get the axes right. I'm trying to use these formulas from Wikipedia.

My code so far is:

raw_moments.m:

function outmom = raw_moments(im,i,j)

  total = 0;
  total = int32(total);
  im = int32(im);

  [height,width] = size(im);

  for x = 1:width;
     for y = 1:height;
        amount = (x ** i) * (y ** j) * im(y,x);
        total = total + amount;
     end;
  end;

  outmom = total;

central_moments.m:

function cmom = central_moments(im,p,q);

  total = 0;
  total = double(total);
  im = int32(im);

  rawm00 = raw_moments(im,0,0);

  xbar = double(raw_moments(im,1,0)) / double(rawm00);
  ybar = double(raw_moments(im,0,1)) / double(rawm00);

  [height,width] = size(im);

  for x = 1:width;
    for y = 1:height;
      amount = ((x - xbar) ** p) *  ((y - ybar) ** q) * double(im(y,x));
      total = total + double(amount);
    end;
  end;

  cmom = double(total);

And here's my code attempting to use these. I include comments for the values I get at each step:

inim = logical(imread('135deg100by30ell.png'));

cm00 = central_moments(inim,0,0);          % 2567

up20 = central_moments(inim,2,0) / cm00;   % 353.94
up02 = central_moments(inim,0,2) / cm00;   % 352.89
up11 = central_moments(inim,1,1) / cm00;   % 288.31

covmat = [up20, up11; up11, up02];
%[ 353.94  288.31
%  288.31  352.89 ]

eigvals = eig(covmat);          % [65.106 641.730]

minoraxislength = eigvals(1);   % 65.106
majoraxislength = eigvals(2);   % 641.730

I'm not sure what I'm doing wrong. I seem to be following those formulas correctly, but my results are nonsense. I haven't found any obvious errors in my moment functions, although honestly my understanding of moments isn't the greatest to begin with.

Can anyone see where I'm going astray? Thank you very much.

Answer 1

EDIT:

According to Wikipedia:

the eignevalues [...] are proportional to the squared length of the eigenvector axes.

which is explained by:

axisLength = 4 * sqrt(eigenValue)

---

Shown below is my version of the code (I vectorized the moments functions):

my_regionprops.m

function props = my_regionprops(im)
    cm00 = central_moments(im, 0, 0);
    up20 = central_moments(im, 2, 0) / cm00;
    up02 = central_moments(im, 0, 2) / cm00;
    up11 = central_moments(im, 1, 1) / cm00;

    covMat = [up20 up11 ; up11 up02];
    [V,D] = eig( covMat );
    [D,order] = sort(diag(D), 'descend');        %# sort cols high to low
    V = V(:,order);

    %# D(1) = (up20+up02)/2 + sqrt(4*up11^2 + (up20-up02)^2)/2;
    %# D(2) = (up20+up02)/2 - sqrt(4*up11^2 + (up20-up02)^2)/2;

    props = struct();
    props.MajorAxisLength = 4*sqrt(D(1));
    props.MinorAxisLength = 4*sqrt(D(2));
    props.Eccentricity = sqrt(1 - D(2)/D(1));
    %# props.Orientation = -atan(V(2,1)/V(1,1)) * (180/pi);      %# sign?
    props.Orientation = -atan(2*up11/(up20-up02))/2 * (180/pi);
end

function cmom = central_moments(im,i,j)
    rawm00 = raw_moments(im,0,0);
    centroids = [raw_moments(im,1,0)/rawm00 , raw_moments(im,0,1)/rawm00];
    cmom = sum(sum( (([1:size(im,1)]-centroids(2))'.^j * ...
                     ([1:size(im,2)]-centroids(1)).^i) .* im ));
end

function outmom = raw_moments(im,i,j)
    outmom = sum(sum( ((1:size(im,1))'.^j * (1:size(im,2)).^i) .* im ));
end

... and the code to test it:

test.m

I = imread('135deg100by30ell.png');
I = logical(I);

>> p = regionprops(I, {'Eccentricity' 'MajorAxisLength' 'MinorAxisLength' 'Orientation'})
p = 
    MajorAxisLength: 101.34
    MinorAxisLength: 32.296
       Eccentricity: 0.94785
        Orientation: -44.948

>> props = my_regionprops(I)
props = 
    MajorAxisLength: 101.33
    MinorAxisLength: 32.275
       Eccentricity: 0.94792
        Orientation: -44.948

%# these values are by hand only ;)
subplot(121), imshow(I), imdistline(gca, [17 88],[9 82]);
subplot(122), imshow(I), imdistline(gca, [43 67],[59 37]);