从第二中心矩计算对象统计信息

Question

我目前正在编写一个 GNU Octave 版本的 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

我不确定我做错了什么。虽然我似乎正确地遵循了这些公式，但我的结果却是无意义的。虽然老实说我对矩的理解本来就不是很好，但我在我的矩函数中没有发现任何明显的错误。

有人能看出我哪里出错了吗？非常感谢。

- BigBeagle

Answer 1

编辑：

根据维基百科的介绍:

特征值[...]与特征向量轴的平方长度成比例。

这可以通过以下方式解释:

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

...并测试代码：

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

- Amro

---

是的，那个角度公式似乎有效。但我仍然不明白为什么轴的值不正确。 - BigBeagle

非常感谢你。我试过了，似乎运行得很好。而且我很感激你向我展示了更好的编写函数的方法。谢谢。 - BigBeagle

Answer 2

你确定你的 raw_moments 函数核心部分没问题吗？你可以尝试：

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

这似乎不足以引起您所看到的问题，但它可能至少是其中的一部分。

- mtrw

---

我尝试更改了那个，以及central_moments中的相同内容，但对最终结果没有任何影响。不过还是谢谢。 - BigBeagle

在MATLAB中，**运算符代表什么？ - glglgl

@glglgl - 在MATLAB中的^可以在Octave（以及Python和可能是Fortran）中写成**。 - mtrw

但是根据我所知，这不是在MATLAB中完成的...所以这是一种不兼容性（但很抱歉，我看到了“octave”标签和目标是Octave的句子有点晚...） - glglgl

此答案应在评论中。 - masad