Calculation of statistics of objects from the second central points

Question

I am currently working on writing a version of the RegionProps function for GNU Octave . I have implemented most of this, but I am still struggling with the implementation of several parts. Earlier, I asked about the second central points of the region.

It was theoretically useful, but I had problems with the implementation of the proposals. I get results wildly different from MATLAB (or common sense, for that matter), and really don't understand why.

Consider this test image:

We see that it is tilted 45 degrees along the X axis, with the small and large axes 30 and 100, respectively.

Running through the MATLAB RegionProps function confirms this:

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

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

My code so far:

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 so my code is trying to use them. I include comments for the values ​​that 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 these formulas correctly, but my results are stupid. I did not find any obvious errors in my instant functions, although, to be honest, my understanding of the moments is not the biggest start.

Can anyone see what I'm messing with? Thank you very much.

Answer 1

EDIT:

According to Wikipedia :

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

which is explained by:

 axisLength = 4 * sqrt(eigenValue) 

---

Below is my version of the code (I have vectorized the functions of the moments):

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 verify 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]); 

Answer 2

Are you sure the core of your raw_moments function? You can try

 amount = ((x-1) ** i) * ((y-1) ** j) * im(y,x); 

This does not seem to cause the problems you see, but it can be at least part.