主轴长度测量是否等同于regionprops3输出中的物体直径？

Question

主轴长度测量是否等同于regionprops3输出中的物体直径？

3

我正在使用MATLAB中的regionprops3命令来计算物体的属性。我以为PrincipalAxisLength是我的物体的直径，但它并不是。

我创建了一个包含椭球形状的二进制图像I，其中半径分别为(40, 10, 10)。从regionprops3中得到以下数值：

sum(I(:)) = 16741
stats.volume = 16741
stats.EquivDiameter = 31.739  % 4/3*pi*(31.739/2)^3=16741
stats.PrincipalAxisLength = [71.535, 17.936, 17.908]

如果PrincipalAxisLength是一个椭球体的3个直径，则其体积为：

4/3*pi*71.535/2*17.936/2*17.908/2 = 12031

这与上述体积不等价。

那么，我该如何计算椭球的半径？

- ankit agrawal

2个回答

3

我曾经面临同样的问题，需要识别对应于二进制区域的椭球体。通过发展数学，可以检索出一个因子sqrt(5)，用于将特征值的平方根转换为椭球体半径。

Matlab的regionprops3函数使用因子4，因此因子sqrt(5/16)解决了差异。

数学涉及使用球坐标的积分和三角表达式的线性化。由于我找不到任何其他说明文档或资源，因此我建议在这里写下它（希望这是适当的...）。

发展

定义

设a>b>c为椭球体的三个半轴。我们可以假设椭球体与三个主轴居中并对齐，其中x轴沿最长半径，z轴沿最短半径。在这种情况下，特征值对应于三个二阶中心矩lambda_1 = mu200、lambda_2 = mu020和lambda_2 = mu002。

然后的想法是从椭球体的参数明确计算这些值。

对于mu200，定义如下：

mu200 = \int \int \int I(x,y,z) x^2 dx dy dz

当椭球内时，I（x，y，z）为1，否则为0。

使用球坐标进行积分

对于积分，使用球坐标 （rho，theta，phi） 更方便，对应半径，与竖直线的倾斜和方位角。它们对应：

x = rho * a * cos(phi) * sin(theta)
y = rho * b * sin(phi) * sin(theta)
z = rho * c * sin(theta)

改变坐标系需要使用变换的雅可比矩阵。在这种情况下，它的行列式为

a * b * c * rho^2 * sin(theta)

mu200的瞬时积分为：

mu200 = a * b *c \int_0^{2*pi} \int_0^pi \int_0^1 (a * rho * cos(phi) * sin(theta) )^2 * rho^2 * sin(theta) drho dtheta dphi

可以简化为

mu200 = ( a^3 * b * c / 5 ) \int_0^{2*pi} \int_0^pi (sin(theta))^3 dtheta (cos(phi))^2 dphi

使用正弦函数的三次方的线性化，可以对 theta 的积分进行求解:

sin^3(theta) = (3/4) * sin(theta) - (1/4) * sin(theta)

然后内部积分变为：

\int_0^pi (sin(theta))^3 dtheta = 4/3

将其纳入完整的积分中：

mu200 = ( a^3 * b * c / 5 ) * (4/3) \int_0^{2*pi} (cos(phi))^2 dphi
mu200 = ( a^3 * b * c / 5 ) * (4/3) * (1/2) 2*pi
mu200 = (4 * pi / 3) * (a^3 * b * c / 5)

鉴定

将椭球体的体积表示为m000 = (4 * pi / 3) * a * b * c，可以得到：

mu200 = (a^2 / 5) * V

在regionprops3函数中，协方差矩阵是通过除以体积进行归一化得到的。

因此我们有：

lambda_1 = a^2 / 5

和关系：

a = sqrt(5) * sqrt(lambda1)

- dlegland

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- gnovice · Accepted Answer

PrincipalAxisLength形状测量的文档描述为（重点是我）：

具有与区域相同归一化二阶中心矩的椭球体主轴长度（以体素为单位），返回为1x3向量。 regionprops3按值从高到低排序。

换句话说，它将多元概率分布拟合到该区域，并计算该分布的中心矩，返回PrincipalAxisLength测量的二阶中心矩（即方差）。

您看到的体积差异是由于这些中心矩不定义用于适应区域内数据的边界椭球体，而只是用于其内部分布数据的概率拟合的方差。对于包含椭球形状的区域，它会低估椭球的范围（因此总体积）。以下是一些可视化代码:

% Create 3D binary data containing ellipsoid:
[X, Y, Z] = meshgrid(-50:50, -20:20, -20:20);
R = [40 10 10];
I = ((X./R(1)).^2 + (Y./R(2)).^2 + (Z./R(3)).^2 <= 1);

% Calculate statistics with regionprops3:
stats = regionprops3(I, 'all');
L = stats.PrincipalAxisLength./2;

% Create ellipsoid surface from second central moments and plot:
[x, y, z] = ellipsoid(0, 0, 0, L(1), L(2), L(3), 100);
surf(x, y, z, 'FaceAlpha', 0.5, 'FaceColor', 'r', 'EdgeColor', 'none');
axis equal;
hold on;

% Create image of middle slice through ellipsoid and plot:
imageSlice = 256.*uint8(I(:, :, 21));
imageSlice = cat(3, imageSlice, imageSlice, imageSlice);
image(-50:50, -20:20, imageSlice);
view(0, 90);

这是最终的图形：

如果你想要将一个椭球面拟合到二进制对象上，可以使用 PrincipalAxisLength 的测量结果作为真实轴长的起始猜测值。在“File Exchange”上有这个函数可以进行椭球体拟合，但那个函数使用一组表面点而不是二进制体积。或许你可以调整那个代码以适应你的需求。

更新

在Cris的评论中提供的链接表明，第二中心矩和椭圆半径之间可能存在明确的数学关系。虽然我还没有时间处理这个数学问题，但我注意到仅通过缩放 PrincipalAxisLength 乘以sqrt(5/16)就可以得到非常接近R的值：

>> sqrt(5/16).*stats.PrincipalAxisLength

ans =

  40.008372204885049   9.970908565527971   9.970908565527971