所需要的软件:Slicer或SlicerSALT(如果使用Slicer,需要单独安装SPHARM-PDM工具箱,SPHARM-PDM新版本已作为Slicer的扩展程序和SlicerSALT的一部分,不再单独提供)
整体流程:STL网格模型 → 二值label → 球面参数化 → SH展开 → PCA
1. 从stl网格模型得到二值label
(1)加载STL文件
- 打开 3D Slicer;
- File → Add Data → 选择STL文件;
(2)把STL变成 Segmentation
- 切到 Segmentations 模块;
- 点击“Add segment”,新建一个 segmentation;
- 在下方 Export/import models and labelmaps 区域:
- Operation选择Import
- Input type选择Models
- Input node选择导入的STL模型
- 点击Import

(3)导出labelmap
- 切到 Segmentations模块;
- 在下方 Export/import models and labelmaps 区域:
- Operation选择Export
- Output type选择Labelmap
- 点击Export

(4)修改体素分辨率
- 切到 Volumes 模块;
- 查看体素数量,体素大小,图像原点等信息;
- Image Dimensions体素数量,代表i, j, k三个方向上的体素数量;
- Image Spacing体素大小,数值越小,分辨率越高,可以类比位图中的像素,需要注意的是直接修改这个值,不会改变volume的体素大小,想要改变体素大小,需要对label进行重采样;
- Image Origin图像原点,表示体素索引的物理坐标,决定整个体积在世界坐标中的摆放位置;
- 点击Center Volume按钮可以自动把volume的中心放到坐标原点 (0,0,0) 位置;
- 重采样label,改变体素大小(可选,如果体素数量太多);
- 切到Resample Scalar Volume模块
- Spacing更改为想要的体素大小
- Interpolation插值选择nearestNeighbor
- Input Volume选择要重采样的volume
- Output Volume新建一个volume

(5)保存数据
File → Save Data,保存项目数据和二值label nrrd文件
2. 球面参数化
(1)拓扑修正与预处理
- SlicerSALT导入上面保存的二值label nrrd文件;
- 导入进来的volume,而不是labelmap,无法进行SegPostProcess,需要先进行转换。切到Volumes模块,找到Convert to labelmap,点击Convert;
- 切到SegPostProcess模块;
- Input Image选择导入的voulume
- Output Image新建voulume
- Extract this label,如果volume是二值图像(背景=0,模型=1),填1,表示只取值为1的那一类
- Scale to isotropic pixel dimensions,体素重采样成各向同性
- Enforced spacing in x, y and z directions,体素间距大小

(2)球面参数化
- 切到GenParaMesh模块;
- 设置参数:
- Input Volume选择要做参数化的volume
- Output Para Mesh,被映射到“参数空间”(通常是球面坐标)的网格
- Output Surface Mesh,原来物理空间中的表面,它的点的数量和顺序与 Para Mesh 一一对应
- Number of iterations,参数化优化的最大迭代次数
- Label id,如果volume是二值图像(背景=0,模型=1),填1

3. SH展开
(1)切到ParaToSPHARMMesh模块
- 切到ParaToSPHARMMesh模块
- 参数设置:
- Input Para mesh dataset,选GenParaMesh中生成的Para Mesh
- Input Surface mesh dataset,选GenParaMesh中生成的Surface Mesh
- Output Directory and base filename,选择文件输出路径
- Subdivision level for linear icosahedron subdivision,球面采样的分辨率,数值越大,球面网格越密。要注意的是这里采用的是基于频率K的正二十面体细分网格,而不是递归细分等级的正二十面体细分网格,两者计算方法不一样。具体如下:
- 频率 K 的 geodesic icosahedron:
- 递归细分level L的icosphere:
两者的对应关系:
- 频率 K 的 geodesic icosahedron:
两者的对应关系:
| 递归 level L | 对应 | ||
| 0 | 1 | 12 | 20 |
| 1 | 2 | 42 | 80 |
| 2 | 4 | 162 | 320 |
| 3 | 8 | 642 | 1280 |
| 4 | 16 | 2562 | 5120 |
| 5 | 32 | 10242 | 20480 |
- Degree of spherical harmonic expansion,球谐阶数
- Use the Fliptemplate for normalization,是否使用主轴对齐,强烈建议开启,特别是后续若要进行PCA的话
- Fliptemplate for normalization of axis-flips,输入用于主轴对齐的模板文件,当Use the Fliptemplate for normalization选项打开时生效,输入文件格式为.coef文件。
(3)输出文件说明:
_para.vtk:一个在单位球上的三角网格,每个顶点带有标量(通常是从原表面映射来的坐标或颜色),用来表示“原形状的哪个点被拉到了这个球面点上”;
_paraTheta.txt / _paraPhi.txt / _paraPhiHalf.txt / _paraMix.txt:球面坐标的采样信息,Theta:极角(0-π),Phi:方位角(0-2π);PhiHalf、paraMix 是内部用的一些组合/映射;
_SPHARM.vtk:用原始坐标系(未做椭球对齐)的 SPHARM 重建表面;
_SPHARM.coef:对应_SPHARM.vtk那一套坐标系下的球谐系数;
_SPHARM_ellalign.vtk:经过椭球对齐(ellipsoid alignment)后的 SPHARM 重建表面;
_SPHARM_ellalign.coef:统一对齐后的坐标系下的球谐系数;
_SPHARMMedialAxis.vtk:基于 SPHARM 得到的形状中轴(medial axis)网格;
_MedialAxisScalars.csv 和 _SPHARMMedialAxis.vtk:对应的标量信息。

4. PCA
(1)切到Population Analysis模块
(2)参数设置
- Creation of CSV File for Classification Groups
- 在Group中定义组的编号,之后选择要进行PCA的vtk文件所在文件夹,然后点击Export导出csv文件。
- 在Group中定义组的编号,之后选择要进行PCA的vtk文件所在文件夹,然后点击Export导出csv文件。
- Preview/Update Groups
在这里可以看到导入的需要进行PCA的文件,点击Preview可以预览模型。 - PCA
- Minimum explained variance,可以统计到的最小解释方差值
- Maximum number of eignvalues,最大的主成分个数
- 下面的滑块可以在1~3个标准差之间进行调节





