如何改进isosurface方法，实现与regionprops3相当的3D表面积计算精度？

我需要在MATLAB中计算3D物体（二值掩码）的球形度（sphericity），尤其是针对**非各向同性体素间距（non-isotropic voxel spacing）的场景，因此必须获得精确的表面积（surface area）**测量结果。

我在各向同性体素的同一二值球体上对比了两种计算方法（代码见下文）。原本预期如果第二种方法在各向同性体素图像上能达到与regionprops3一致的结果，那么它在非各向同性体素图像上也能给出正确结果，但实际对比结果显示两者精度差异很大：

对比结果 #

regionprops3法 #

• 体积 = 4169
• 表面积 = 1262.5112
• 球形度 = 0.99221

isosurface + 叉积法 #

• 体积 = 4169
• 表面积 = 1366.8365
• 球形度 = 0.91648

可以看到，第二种方法的表面积精度远低于regionprops3，对应的球形度精度也更低。

问题 #

能否在不修改输入3D二值掩码的前提下改进第二种方法，使其精度接近regionprops3？

clc; clear; close all;

% ----------------------
% 创建二值球体
% ----------------------
N = 50;                  % 体素尺寸 (50x50x50)
mask = zeros(N, N, N);   % 3D二值矩阵
radius = 10;             % 球体半径
center = [25, 25, 25];   % 球体中心
% 创建3D网格
[x, y, z] = ndgrid(1:N, 1:N, 1:N);
mask(((x-center(1)).^2 + (y-center(2)).^2 + (z-center(3)).^2) <= radius^2) = 1;

% ----------------------
% regionprops3方法
% ----------------------
stats = regionprops3(mask, 'Volume', 'SurfaceArea');
V_rp = stats.Volume;
A_rp = stats.SurfaceArea;
psi_rp = (pi^(1/3) * (6 * V_rp)^(2/3)) / A_rp;

disp('--- regionprops3 method ---');
disp(['Volume = ', num2str(V_rp)]);
disp(['Surface Area = ', num2str(A_rp)]);
disp(['Sphericity = ', num2str(psi_rp)]);

% ----------------------
% isosurface + 叉积法
% ----------------------
dx = 1; dy = 1; dz = 1;  % 各向同性体素间距 
fv = isosurface(mask, 0.5);

% 将顶点缩放至真实世界单位
fv.vertices(:,1) = fv.vertices(:,1) * dx;
fv.vertices(:,2) = fv.vertices(:,2) * dy;
fv.vertices(:,3) = fv.vertices(:,3) * dz;

% 计算体积
V_iso = nnz(mask) * dx * dy * dz;

% 计算表面积
A_iso = 0;
for i = 1:size(fv.faces,1)
    v1 = fv.vertices(fv.faces(i,1),:);
    v2 = fv.vertices(fv.faces(i,2),:);
    v3 = fv.vertices(fv.faces(i,3),:);
    A_iso = A_iso + 0.5 * norm(cross(v2-v1, v3-v1));
end

psi_iso = (pi^(1/3) * (6 * V_iso)^(2/3)) / A_iso;

disp('--- isosurface + cross method ---');
disp(['Volume = ', num2str(V_iso)]);
disp(['Surface Area = ', num2str(A_iso)]);
disp(['Sphericity = ', num2str(psi_iso)]);

内容的提问来源于stack exchange，提问作者Ellie