从第二个中心时刻计算对象统计信息

[poster]

我目前正在为GNU Octave编写一个MATLAB RegionProps函数版本。我已经实现了大部分，但是我仍在为几个部分的实现而苦苦挣扎。我以前曾问过一个地区的第二个中心时刻 。

从理论上讲这是有帮助的，但是我在实际执行建议时遇到了麻烦。我得到的结果与MATLAB的结果（或该问题的常识）完全不同，并且真的不明白为什么。

考虑以下测试图像：



我们可以看到它与X轴成45度倾斜，短轴和长轴分别为30和100。

通过MATLAB的RegionProps函数运行它可以确认这一点：

MajorAxisLength: 101.3362 MinorAxisLength: 32.2961 Eccentricity: 0.9479 Orientation: -44.9480 同时，我什至没有正确的轴。我正在尝试使用来自维基百科的这些公式 。

到目前为止，我的代码是：

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); 这是我的代码尝试使用这些代码。我包括对每个步骤得到的值的注释：

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 我不确定自己在做什么错。我似乎正确地遵循了这些公式，但是我的结果是胡说八道。我的矩函数没有发现任何明显的错误，尽管说实话，我对矩的理解并不是最重要的。

谁能看到我误入歧途？非常感谢你。


回答：
编辑：

根据维基百科 ：

eignevalues与特征向量轴的平方长度成正比 。

解释如下：

axisLength = 4 * sqrt(eigenValue) 下面显示的是我的代码版本（我将矩函数矢量化了）：

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 ...以及测试代码：

测试

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



更多&回答...