从第二个中心矩计算对象统计量 [英] Computing object statistics from the second central moments
问题描述
我目前正在编写一个版本的 MATLAB .
到目前为止我的代码是:
raw_moments.m:
function outmom = raw_moments(im,i,j)总计 = 0；总计 = int32(总计)；im = int32(im);[高度，宽度] = 尺寸(im)；对于 x = 1:宽度；对于 y = 1:高度；金额 = (x ** i) * (y ** j) * im(y,x);总计 = 总计 + 金额；结尾;结尾;outmom = 总计；central_moments.m:
function cmom = central_moments(im,p,q);总计 = 0；总计 = 双倍(总计)；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);[高度，宽度] = 尺寸(im)；对于 x = 1:宽度；对于 y = 1:高度；金额 = ((x - xbar) ** p) * ((y - ybar) ** q) * double(im(y,x));总计 = 总计 + 双倍(金额)；结尾;结尾;cmom = 双倍(总计)；这是我尝试使用这些的代码.我包括对我得到的值的评论在每一步:
inim = logical(imread('135deg100by30ell.png'));cm00 = central_moments(inim,0,0);% 2567up20 = central_moments(inim,2,0)/cm00;% 353.94up02 = central_moments(inim,0,2)/cm00;% 352.89up11 = central_moments(inim,1,1)/cm00;% 288.31covmat = [up20, up11;上11，上02]；%[ 353.94 288.31% 288.31 352.89 ]eigvals = eig(covmat);% [65.106 641.730]小轴长度 = eigvals(1);% 65.106长轴长度 = eigvals(2);% 641.730我不确定我做错了什么.我似乎正确地遵循了这些公式，但我的结果是无稽之谈.我在矩函数中没有发现任何明显的错误，尽管老实说，我对矩的理解并不是最好的.
任何人都可以看到我误入歧途的地方吗?非常感谢.
根据维基百科:
<块引用>特征值 [...] 成比例为特征向量轴的平方长度.
解释如下:
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');%# 排序 cols 从高到低V = V(:,订单);%# 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.MajorAxisLength = 4*sqrt(D(1));props.MinorAxisLength = 4*sqrt(D(2));props.离心率 = 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);结尾函数 cmom = central_moments(im,i,j)rawm00 = raw_moments(im,0,0);质心 = [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 ));结尾函数 outmom = raw_moments(im,i,j)outmom = sum(sum( ((1:size(im,1))'.^j * (1:size(im,2)).^i) .* im ));结尾...以及测试它的代码:
test.m
I = imread('135deg100by30ell.png');I = 逻辑(I)；>>p = regionprops(I, {'偏心' 'MajorAxisLength' 'MinorAxisLength' '方向'})p =长轴长度:101.34小轴长度:32.296偏心率:0.94785方向:-44.948>>道具 = my_regionprops(I)道具=长轴长度:101.33小轴长度:32.275偏心率:0.94792方向:-44.948%# 这些值只能手动计算 ;)子图(121)，imshow(I)，imdistline(gca，[17 88]，[9 82])；子图(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屋！
