球面谐波函数 （Spherical harmonic function）分析实际颗粒形状公式推导及数值实现 Part 2 颗粒形状的球面谐波分析，借助3D Slicer的SPHARM实现

一、3D Slicer软件简介及安装

       Slice 3D是一个免费的开源软件，用于医学、生物医学和其他3D图像和网格的可视化、处理、分割、匹配和分析，以及规划和导航图像引导的程序。支持Linux、MacOS、Windows操作环境，目前在医学影像学方面应用较多。

        官方网站：https://www.slicer.org/

        此外，为实现本研究中用到的球面谐波分析方法，还需要再3D Slicer中安装SPHARM-PDM插件，以下是安装方法（以5.4.0为例）：

（1）进入软件后，选择箭头所指的图标；

（2）搜索SPHARM-PDM，并选择安装

（3）安装Shape Population Viewer，以查看颗粒的形状。

（4）重启软件，即完成安装。

二、文件准备

       图像分割用到的是阈值分割，相互粘结的颗粒之间的分割，使用分水岭算法，给每个单独的颗粒设置编号，并提取出某个（某个区域中）颗粒的的所有单独颗粒的体素信息，并转化为3D slicer能够读取的数据（*.gipl.gz或*.nii）。相应代码如下：

其中需要提前下载*.nii文件生成的Matlab插件 Tools for NIfTI and ANALYZE image 

网址：https://ww2.mathworks.cn/matlabcentral/fileexchange/8797-tools-for-nifti-and-analyze-image

close all;clear all
% 读取.raw数据
file_path = '02.separate.raw';  % 你自己的原始数据
width = 501; % 数据宽度
height = 501; % 数据高度
depth = 501; % 数据深度

fid = fopen(file_path, 'rb');
data = fread(fid, width * height * depth, '*uint8'); % 注意.raw文件是8位的
fclose(fid);

% 将数据重新整形为三维矩阵
data = reshape(data, [width, height, depth]);
% data = data(100:400,100:400,100:400);
% 计算距离变换
D = -bwdist(~data);
D = imhmin(D, 1.2); % 可选：限制最小深度
D(~data) = Inf;

% 执行分水岭算法
L = watershed(D);
L(~data) = 0;
a = L(:,:,100);imshow(a,[]),colormap(jet)
figure;b = D(:,:,100);imshow(b,[]);colormap(jet)
figure;imshow(data(:,:,100),[])

% 显示编号/特定形状的颗粒
targetLabel = 620;
bwim = L == targetLabel;
% 使用 isosurface 函数提取表面网格
% isoValue = 0.3;  % isoValue 控制着提取的表面的阈值，可以根据数据特性调整
[faces, vertices] = isosurface(bwim);

% 显示表面网格
figure;
p = patch('Vertices', vertices, 'Faces', faces);
set(p, 'FaceColor', 'g', 'EdgeColor', 'none');
daspect([1 1 1]);
view(3);
camlight;
lighting gouraud;
title(['Isosurface of Label ', num2str(targetLabel), ' Particle']);
hold on
% 显示提取的表面点
X = vertices(:,1);
Y = vertices(:,2);
Z = vertices(:,3);
X = X';
Y = Y';
Z = Z';
scatter3(X,Y,Z,10,'filled');

% 获取某个编号的颗粒的所有体素
xScaleMin = min(Y);
xScaleMax = max(Y); 
yScaleMin = min(X); % 注意matlab在存储过程中会有一个转置
yScaleMax = max(X);
zScaleMin = min(Z);
zScaleMax = max(Z);
i_th_Particle_Voxel = bwim(xScaleMin:xScaleMax,yScaleMin:yScaleMax,zScaleMin:zScaleMax);

% 保存到指定的文件格式
i_th_Particle_Voxel  = uint8(i_th_Particle_Voxel);
nii = make_nii(i_th_Particle_Voxel);

save_nii(nii,'i_thParticle.nii');

三、操作步骤

1. 分成三个主要的步骤，在软件界面中也是如此

2. 具体流程

（1）在Modules模块选择Shape Analysis Module模块；

（2）设置输入输出

程序会将输入文件中的所有(.gipl, .gipl.gz, .nii, .nii.gz, .nrrd, .mgh,.mgh.gz, .mhd, or .hdr）格式的文件读入；

输出会在输出文件夹中生成三个文件夹：分别是三个步骤中的输出文件。

（3）随后在后续的步骤中设置不同的，并运行程序，详细的截图如下：

(4) 输出文件：

(5)可视化处理