从第二个中心时刻计算对象统计 [英] Computing object statistics from the second central moments
问题描述
我目前正在编写一个MATLAB版本。
我的代码到目前为止：
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
我不知道我做错了什么。我似乎正确地遵循这些公式，但我的结果是废话。我没有发现任何明显的错误在我的时刻功能，虽然老实说，我对时刻的理解不是最好的开始。
任何人都可以看到我在哪里走错了？非常感谢。
编辑：
a href =http://en.wikipedia.org/wiki/Image_moment#Central_moments =nofollow noreferrer>维基百科：
eignevalues [...]是成比例
到特征向量轴的平方长度。
其解释如下：
axisLength = 4 * sqrt（eigenValue）
下面显示的是我的代码版本
my_regionprops.m
function props = my_regionprops b $ b 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从高到低
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） ％＃ 标志？
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 * ...
] -centroids（1））。
end
function outmom = raw_moments（im，i，j）
outmom = sum（sum（（1：size（im，1））' （1：size（im，2））。
end
...以及测试代码：
test.m
I = imread（'135deg100by30ell.png'）;
I = logical（I）;
>> p = regionprops（I，{'Eccentricity'MajorAxisLength''MinorAxisLength''Orientation'}）
p =
MajorAxisLength：101.34
MinorAxisLength：32.296
偏心率：0.94785
定向：-44.948
>> props = my_regionprops（I）
props =
MajorAxisLength：101.33
MinorAxisLength：32.275
偏心率：0.94792
方向：-44.948
％ ＃这些值是手动的;）
子图（121），imshow（I），imdistline（gca，[17 88]，[9 82]）;
subplot（122），imshow（I），imdistline（gca，[43 67]，[59 37]）;
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.
解决方案 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]);
这篇关于从第二个中心时刻计算对象统计的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
