18、形状模型优化与评估：迭代更新与质量评测

形状模型优化与评估：迭代更新与质量评测

在形状分析领域，优化和评估统计模型是至关重要的环节。本文将介绍使用迭代更新形状图像的优化例程，以及如何利用地面真值数据评估这些模型的质量。

1. 迭代更新形状图像的优化例程

在优化过程中，常常需要在迭代方法中使用形状图像。下面将分别介绍开放表面和封闭表面的优化示例。

1.1 开放表面优化示例

这里使用的开放表面数据集是 10 个远端股骨的集合。以下是优化开放表面对应关系的算法伪代码：

Algorithm 8.4 : An Example of Optimising Correspondence on Open Surfaces
Using Shape Images.
procedure {shape vectori} ←optimisation example 3
({Shape Pointsi}, {Triangulationi})
Variables
• Shape Pointsi is a nP ×3 matrix whose rows contain the coordinates of nodes that define
the surface of the ith training example;
• Triangulationi is a nt × 3 matrix that defines the connectivity of the nodes: each row
defines a triangle by indexing the rows of Shape Pointsi – see glossary for description of
format;
• the function returns a set of 3n-dimensional shape vectors {shape vectori}, formed by
sampling each training set according to its optimal parameterisation.
Declarations
• y = transfer point(Points A, Points B, Triangulation, x)
transfers the point x from triangulation (Points A, Triangulation) to triangulation
(Points B, Triangulation) to form a new point y – the pseudocode is given in Al-
gorithm 7.3;
• b = bilinear interpolation(X, Y; a)
the n × 2 matrix X, holds values of a (vector-valued) function, evaluated at points whose
coordinates are stored in the n × 2 matrix Y; the function returns b, the value of the
function at point a, which is estimated by bilinear interpolation;
• {Warpi} = optimise fluid({Shape Imagei}, Grid Points)
given a set of shape images {Shape Imagei}, with pixels at positions Grid Points,
this function optimises correspondence by re-parameterisation using a fluid-based scheme;
the function returns the set of warps {Warpi}, which correspond to the optimal set of
re-parameterisations – pseudocode is given in Algorithm 8.3.
Initialization
1. let ref be the index of a reference shape, whose pose and parameterisation do not change
during optimisation;
2. parameterise and align the surfaces as described in step 2 of Algorithm 7.4; we end up
with a set of parameterisations, {Parameterisationi}, represented as a set of nP ×
2 matrices, whose rows contain the coordinates of nodes in parameter space (the unit
square), corresponding to the surface nodes stored in Shape Pointsi;
3. create a m × m grid over the unit square (m2 = n = 1600 was used here) and store the
coordinates in a n × 2 matrix Grid Points:
Grid Points =
[(0, 0)T ,
(0, 1/(m - 1))T , ... ,
(1/(m - 1) , 0)T ,
(1/(m - 1), 1/(m - 1))T , ... , (1, 1)T ]T ;
4. create a n × 2 matrix Sample Pointsi that holds the sample points for each training
example – the sample points are initially set to be regular grid points:
for i = 1 ... nS
Sample Pointsi ←Grid Points
Optimisation
for it = 1 ... niterations (we used niterations = 1500 here)
1. create a shape image for each example by sampling each shape according to the current
parameterisation:
for i = 1 ... nS, (if it > 1 , then i ≠ ref)
for j = 1 ... n
Shape Imagei(j, ·) ←transfer point(Parameterisationi,
Shape Pointsi, Triangulationi, Sample Pointsi(j, ·)).
2. optimise the set of shape images by re-parameterisation using fluid regularization; the
function returns a set of n × 2 matrices that represent the optimal re-parameterisation of
each shape;
{Reparameterisationi} ←fluid optimisation
({Shape Imagei}, Grid Points)
3. re-parameterise the sample points of each shape using the optimal re-parameterisation:
for i = 1 ... nS, i ≠ ref
for j = 1 ... n
Sample Pointsi(j, ·) ←bilinear interpolate(Grid Points,
Reparameterisationi, Sample Pointsi(j, ·))
Post-processing
1. use the optimal sample points to create a shape vector for each example:
for i = 1 ... nS
for j = 1 ... n
v ←transfer point(Parameterisationi,
Shape Pointsi, Triangulationi, Sample Pointsi(j, ·)),
shape vectori(j) ←v(1),
shape vectori(j + n) ←v(2),
shape vectori(j + 2n) ←v(3).
return {shape vectori}.

该算法的主要步骤如下：
1. 初始化 ：
- 选择一个参考形状的索引。
- 对表面进行参数化和对齐。
- 在单位正方形上创建网格。
- 初始化样本点为网格点。
2. 优化 ：
- 根据当前参数化对每个形状进行采样，创建形状图像。
- 使用流体正则化对形状图像集进行重新参数化优化。
- 使用优化后的重新参数化对样本点进行重新参数化。
3. 后处理 ：使用优化后的样本点为每个示例创建形状向量。

该算法与之前使用参数化表示（夹紧板样条）和三角化形状表面的方法产生的对应关系和模型几乎无法区分，但基于流体的算法收敛速度快约 100 倍。

1.2 封闭表面优化示例

封闭表面数据集是 82 个人类海马体的集合。为了处理封闭表面，需要扩展上述算法，并引入将平面上的点映射到单位球体的过程。

以下是将平面上的点映射到球体的算法伪代码：

Algorithm 8.5 : Transforming Points from the Plane onto a Sphere.
procedure sphere point ←plane to sphere(plane point)
Description
Transfers a point from the plane onto the surface of a unit sphere via an octahedron.
Variables
• plane point is a two-dimensional vector, representing the coordinates of the point on the
plane;
• sphere point is a three-dimensional vector, representing the coordinates of the mapped
point after projection onto the unit sphere
Declarations
• y = transfer point(Points A, Points B, Triangulation, x)
transfers the point x from triangulation (Points A, Triangulation) to triangulation
(Points B, Triangulation) to form a new point y – pseudocode is given in Algorithm 7.3.
Initialization
1. construct
a
flattened
triangulated
octahedron,
storing
the
coordinates
in
Flattened Octahedron Points
,
and
the
triangulation
in
a
matrix,
Triangulation (see Fig. 8.3);
Flattened Octahedron Points ←
[(0, 1)T , (0, 1/2 )T , (0, 0)T ,
( 1/2, 1)T , ( 1/2, 1/2)T , ( 1/2, 0)T , (1, 1)T , (1, 1/2 )T , (1, 0)T ]T ;
Triangulation ←
[(5, 8, 6)T , (5, 8, 4)T , (5, 2, 4), (5, 2, 6)T ,
(9, 8, 6)T , (7, 8, 4)T , (1, 2, 4)T , (3, 2, 6)T ]T ;
2. create an octahedron corresponding to the flattened version and store the nodes (with all
nodes of the octahedron at unit length from the origin) in a matrix Octahedron Points:
Octahedron Points ←((0, 0, −1), (0, 1, 0), (0, 0, −1), (1, 0, 0),
(0, 0, 1), (−1, 0, 0), (0, 0, −1), (0, −1, 0), (0, 0, −1)) ;
Projection
1. transfer the point from the flattened octahedron triangulation to the unflattened octahe-
dron:
sphere point ←transfer point(Flattened Octahedron Points,
Octahedron Points, Triangulation, plane point);
2. normalize the point to lie on the surface of the unit sphere:
sphere point ←
sphere point
||sphere point|| .
return sphere point;

以下是优化封闭表面对应关系的主要算法伪代码：

Algorithm 8.6 : An Example of Optimising Correspondence on Closed Surfaces
Using Shape Images.
procedure {shape vectori} ←optimisation example 4
({Shape Pointsi}, {Triangulationi})
Variables
• Shape Pointsi is a nP ×3 matrix whose rows contain the coordinates of nodes that define
the surface of the ith training example;
• Triangulationi is a nt × 3 matrix that defines the connectivity of the nodes: each row
defines a triangle by indexing the rows of Shape Pointsi
• the function returns a set of 3n-dimensional shape vectors {shape vectori}, formed by
sampling each training set according to its optimal parameterisation.
Declarations
• sphere point = plane to sphere(plane point)
takes a point on the plane, whose coordinates are held in the two-dimensional vector
plane point, and projects it onto the sphere via the octahedron to create a point whose
coordinates are held in the three-dimensional vector sphere point. See Algorithm 8.5;
• y = transfer point(Points A, Points B, Triangulation, x)
transfers the point x from triangulation (Points A, Triangulation) to triangulation
(Points B, Triangulation) to form a new point y – pseudocode is given in Algorithm 7.3;
• {Warpi} = optimise fluid({Shape Imagei}, Grid Points)
given a set of shape images {Shape Imagei}, with pixels at positions Grid Points,
this function optimises correspondence by re-parameterisation using a fluid-based scheme;
the function returns the set of warps {Warpi}, which correspond to the optimal set of
re-parameterisations – pseudocode is given in Algorithm 8.3.
Initialization
1. let ref be the index of a reference shape, whose pose and parameterisation do not change
during optimisation;
2. parameterise and align the surfaces as described in step 2 of Algorithm 7.4 – note, however,
that we are dealing here with closed surfaces, hence there is no need to extract and optimise
the boundary nodes; we end up with a set of parameterisations, {Parameterisationi},
represented as a set of nP × 3 matrices, whose rows contain the coordinates of nodes
in parameter space (the surface of the unit sphere), corresponding to the surface nodes
stored in Shape Pointsi;
3. create a m × m grid over the unit square (m2 = n = 1600 was used here) and store the
coordinates in a n × 2 matrix Grid Points:
Grid Points =
[(0, 0)T ,
(0, 1/(m - 1))T , ... ,
(1/(m - 1) , 0)T ,
(1/(m - 1), 1/(m - 1))T , ... , (1, 1)T ]T ;
4. project the gridpoints onto the sphere via an octahedron:
Sphere Grid Points ←plane to sphere(Grid Points);
5. create a ntg×3 matrix Grid Triangulation that defines the connectivity of the spherical
gridpoints: each row defines a triangle by indexing the rows of Sphere Grid Points;
6. create two copies of Sphere Grid Points: Sphere Grid Points A is an exact copy,
whereas the signs of the z-coordinates in Sphere Grid Points B are reversed:
Sphere Grid Points A ←Sphere Grid Points,
Sphere Grid Points B ←Sphere Grid Points,
Sphere Grid Points B(·, 3) ←−Sphere Grid Points B(·, 3);
what we end up with are two spherical projections of the gridpoints, with the (fixed) edges
of the grid at a different position on the sphere;
7. create two n×3 matrices, Sample Points Ai and Sample Points Bi for each shape:
their rows contain the coordinates of sample points; they are initialized to the spherical
projection of the gridpoints:
for i = 1 ... nS
Sample Points Ai ←Sphere Grid Points A,
Sample Points Bi ←Sphere Grid Points B.
Optimisation
for it = 1 ... niterations (we used niterations = 12500 here)
1. for each training example: create a shape image by using the current sample points to
sample the surface coordinates of the training example and store them in a n × 3 matrix
Shape Imagei; note that if the iteration number, it, is an odd number, then sample
points A are used otherwise sample points B are used:
for i = 1 ... nS, (if it > 1 , then i ≠ ref)
if ((2⌈it/2⌉−it) = 1)
then
for j = 1 ... n
Shape Imagei(j, ·) ←transfer point
(Parameterisationi, Shape Pointsi,
Triangulationi, Sample Points Ai(j, ·))
else
for j = 1 ... n
Shape Imagei(j, ·) ←transfer point
(Parameterisationi, Shape Pointsi,
Triangulationi, Sample Points Bi(j, ·))
2. optimise the set of shape images by re-parameterisation; the function returns a set of n×2
matrices, which represent the optimal re-parameterisation of each shape;
{Warpi} ←optimise fluid({Shape Imagei}, Grid Points)
3. project the re-parameterisation Warpi of each shape onto the sphere:
for i = 1 ... nS, i ≠ ref
for j = 1 ... n
Reparameterisationi(j, ·) ←plane to sphere(Warpi(j, ·))
4. now, apply the (spherical) re-parameterisation to both sets of sample points:
if ((2⌈it/2⌉−it) = 1)
then
for i = 1 ... nS, i ≠ ref
4.1. update both sets of sample points:
for j = 1 ... n
Sample Points Ai(j, ·) ←transfer point
(Sphere Grid Points A, Reparameterisationi,
Grid Triangulation, Sample Points Ai(j, ·))
Sample Points Bi(j, ·) ←transfer point
(Sphere Grid Points A, Reparameterisationi,
Grid Triangulation, Sample Points Bi(j, ·))
else
for i = 1 ... nS, i ≠ ref
4.1. since sample points B has been used, we need to reverse the sign of the
z-coordinates of each re-parameterisation point:
Reparameterisationi(·, 3) ←−Reparameterisationi(·, 3)
4.2. update both sets of sample points:
for j = 1 ... n
Sample Points Ai(j, ·) ←transfer point
(Sphere Grid Points B, Reparameterisationi,
Grid Triangulation, Sample Points Ai(j, ·))
Sample Points Bi(j, ·) ←transfer point
(Sphere Grid Points B, Reparameterisationi,
Grid Triangulation, Sample Points Bi(j, ·))
Post-processing
1. subdivide a unit icosahedron and project its coordinates onto the unit sphere (see
Sect. 7.3.3) to produce a nu × 3 matrix Uniform Sample Points;
2. for each shape: re-parameterise the uniform sample points with either sample points A or
B:
for i = 1 ... nS
for j = 1 ... nu
Reparameterised Sample Pointsi(j, ·) ←transfer point
(Sphere Grid Points A, Sample Points Ai,
Grid Triangulation, Uniform Sample Points(j, ·))
3. now use the reparameterised sample points to sample the original shapes and produce a
shape vector for each example:
for i = 1 ... nS
for j = 1 ... n
v ←transfer point
(Parameterisationi, Shape Pointsi,
Triangulationi, Reparameterised Sample Pointsi(j, ·)),
shape vectori(j) ←v(1),
shape vectori(j + n) ←v(2),
shape vectori(j + 2n) ←v(3).
return {shape vectori}.

该算法的主要步骤如下：
1. 初始化 ：
- 选择参考形状的索引。
- 对表面进行参数化和对齐。
- 创建单位正方形网格并将其投影到球体上。
- 创建网格三角剖分矩阵。
- 创建两份网格点的球形投影副本。
- 初始化每个形状的样本点。
2. 优化 ：
- 根据当前样本点为每个训练示例创建形状图像。
- 对形状图像集进行重新参数化优化。
- 将重新参数化投影到球体上。
- 应用球形重新参数化更新样本点。
3. 后处理 ：
- 细分单位二十面体并将其坐标投影到单位球体上。
- 对均匀样本点进行重新参数化。
- 使用重新参数化的样本点采样原始形状并生成形状向量。

通过该算法找到的海马体对应关系看起来合理，模型的第一变异模式反映了训练集中的形状变化。

2. 统计模型评估

在构建统计形状模型后，需要对其进行评估。评估问题不仅适用于形状模型，也适用于外观模型。评估主要涉及两个方面：表示的质量和基于该表示的模型的质量。

2.1 评估表示质量

假设我们有训练集和测试集，我们可以利用这些原始形状或图像作为地面真值数据来评估表示的质量。

对于每个训练示例，我们有原始形状向量 $x$ 和对应的参数向量 $b$，它们之间的映射由矩阵 $N$ 决定：
$x \rightarrow N^T (x - \bar{x}) = b(x)$，$b \rightarrow \bar{x} + Nb = \tilde{x}(b)$

如果保留所有模式，训练集的重建将是精确的。但在实际中，我们通常只保留足够的模式。为了检查训练形状是否得到充分表示，我们需要比较每个原始形状向量 $x$ 与其重建 $\tilde{x}$。

对于测试形状，即使保留所有模式，也不能保证精确重建，因为测试形状可能不在训练集所张成的子空间内。

比较两个形状时，最简单的定量度量是点到点距离的平方和，即形状空间中的欧几里得距离 $|x - \tilde{x}|$。但欧几里得点到点距离对垂直于形状轮廓或表面的分离和切向于局部形状的分离同等对待。为了考虑这一点，可以使用垂直的点到线或点到面的距离来评估表示的质量。

对于外观模型，通过比较原始形状和图像块与重建的形状和图像块来评估表示的质量，得到形状之间的欧几里得点到点或点到线距离，以及图像空间中两个图像块之间的欧几里得距离。

2.2 获取地面真值数据的方法

• 自动方法 ：一些自动方法可以提取形状信息，但可能无法提供地面真值对应信息。因此，对于这种来源的地面真值数据，我们可以使用点到线距离来验证形状表示，但无法用于详细评估模型。
• 手动或半自动注释 ：对于某些类别的图像和形状，手动或半自动注释可以提供可靠的形状轮廓地面真值信息，但不同注释者之间可能存在显著差异。对于二维图像数据，可重复的地标较难找到；对于三维数据，地标点的放置非常困难。

综上所述，评估形状表示质量所需的地面真值数据可以通过自动或人工注释获得。评估形状表示质量相当于检查我们在形状空间中选择的模型子空间是否合适。接下来，我们将重点关注评估在该模型子空间内构建的模型的质量。

下面是开放表面优化算法的 mermaid 流程图：

graph TD;
    A[初始化] --> B[选择参考形状索引];
    A --> C[参数化和对齐表面];
    A --> D[创建单位正方形网格];
    A --> E[初始化样本点];
    B --> F[优化];
    C --> F;
    D --> F;
    E --> F;
    F --> G[创建形状图像];
    F --> H[优化形状图像集];
    F --> I[重新参数化样本点];
    G --> J[后处理];
    H --> J;
    I --> J;
    J --> K[生成形状向量];

下面是封闭表面优化算法的主要步骤表格：
|步骤|操作|
|----|----|
|初始化|选择参考形状索引、参数化和对齐表面、创建网格并投影到球体、创建网格三角剖分、创建网格点副本、初始化样本点|
|优化|创建形状图像、优化形状图像集、投影重新参数化、更新样本点|
|后处理|细分单位二十面体、重新参数化均匀样本点、生成形状向量|

形状模型优化与评估：迭代更新与质量评测

3. 模型质量的评估

在评估模型质量之前，我们需要明确模型质量的定义，这在很大程度上取决于模型的使用目的。如果用于分析训练集本身，重要的是模型相对于训练集的质量；如果用于分析未见过的形状或图像示例，情况会更复杂，因为表现不佳可能是训练集本身不足，也可能是基于该训练集构建的模型不够好。

3.1 基于地面真值数据的评估

我们已经知道可以利用训练集和测试集的原始形状或图像作为地面真值数据来评估表示质量，同样也可以用这些数据评估模型质量。
- 训练集评估 ：构建模型后，对于每个训练示例，有原始形状向量 $x$ 和重建形状向量 $\tilde{x}$。为了检查模型对训练集的拟合程度，我们比较 $x$ 和 $\tilde{x}$。如前文所述，可使用点到点距离的平方和（即形状空间中的欧几里得距离 $|x - \tilde{x}|$ ）作为简单的定量度量。例如，对于一个由多个顶点构成的形状，计算每个对应顶点之间的距离的平方，然后求和。
- 测试集评估 ：对于测试集，即使保留所有模式，也不能保证精确重建，因为测试形状可能不在训练集所张成的子空间内。同样使用欧几里得距离 $|x - \tilde{x}|$ 来评估模型对测试集的泛化能力。如果测试集的欧几里得距离普遍较大，说明模型的泛化能力可能较差。

然而，欧几里得点到点距离对垂直于形状轮廓或表面的分离和切向于局部形状的分离同等对待。实际上，沿切线移动点对表示形状的有效性影响较小，而垂直于局部形状移动点影响较大。因此，为了更准确地评估，可使用垂直的点到线或点到面的距离。

对于外观模型，评估表示质量时，通过比较原始形状和图像块与重建的形状和图像块，得到形状之间的欧几里得点到点或点到线距离，以及图像空间中两个图像块之间的欧几里得距离。

3.2 地面真值数据获取的局限性

• 自动方法 ：一些自动方法可以提供良好的分割和连续形状轮廓或表面的提取，但缺乏地面真值对应信息。例如，某些图像分割算法可以将物体从背景中分离出来，但无法准确标记每个点的对应关系。所以，对于这种来源的地面真值数据，只能用点到线距离验证形状表示，无法用于详细评估模型。
• 手动或半自动注释 ：对于某些类别的图像和形状，手动或半自动注释可以提供可靠的形状轮廓地面真值信息。但不同注释者之间可能存在显著差异，这会影响评估的准确性。例如，在标注人脸时，不同人对某些特征点的位置判断可能会有偏差。对于二维图像数据，可重复的地标较难找到；对于三维数据，如 MR 图像体积，地标点的放置非常困难甚至不可能。虽然手动标注结构、表面或组织类别的体素级标签等地面真值数据对其他类型的评估有用，但不是模型评估理想的密集点到点对应地面真值数据。

3.3 不依赖地面真值的评估方法

除了基于地面真值数据的评估方法，还可以构建不依赖地面真值的评估方法。例如，可以通过分析模型的内部结构，如特征向量和特征值，来评估模型的质量。特征值反映了每个模式的方差，我们可以根据特征值的分布来选择保留的模式数量。如果某些特征值非常小，说明对应的模式对整体形状变化的贡献较小，可以考虑不保留这些模式。

还可以使用交叉验证的方法，将数据集划分为多个子集，轮流使用其中一个子集作为测试集，其余子集作为训练集，多次训练和测试模型，然后综合评估模型的性能。

4. 总结与展望

通过上述对形状模型优化和评估的介绍，我们了解了使用迭代更新形状图像的优化例程，包括开放表面和封闭表面的优化算法，以及如何评估统计模型的质量。

在优化方面，基于流体的算法在处理开放表面时收敛速度比之前的方法快约 100 倍，在处理封闭表面时也能找到合理的对应关系，模型的变异模式能反映训练集中的形状变化。

在评估方面，评估表示质量和模型质量是构建有效模型的重要环节。虽然利用地面真值数据可以进行一定程度的评估，但存在获取困难和局限性的问题，因此不依赖地面真值的评估方法具有重要的研究价值。

未来的研究方向可以包括进一步改进优化算法，提高收敛速度和优化效果；探索更有效的不依赖地面真值的评估方法，以更准确地评估模型质量；将形状模型应用到更多领域，如医学诊断、计算机视觉等，为实际问题提供更好的解决方案。

下面是模型评估流程的 mermaid 流程图：

graph TD;
    A[准备训练集和测试集] --> B[构建模型];
    B --> C[评估表示质量];
    C --> D[比较训练集形状向量];
    C --> E[比较测试集形状向量];
    D --> F[计算欧几里得距离];
    E --> F;
    F --> G[判断表示质量];
    B --> H[评估模型质量];
    H --> I[基于地面真值评估];
    H --> J[不依赖地面真值评估];
    I --> K[分析训练集拟合度];
    I --> L[分析测试集泛化能力];
    J --> M[分析特征向量和特征值];
    J --> N[交叉验证];
    K --> O[综合评估模型];
    L --> O;
    M --> O;
    N --> O;

下面是评估方法对比表格：
|评估方法|优点|缺点|适用场景|
|----|----|----|----|
|基于地面真值的欧几里得点到点距离|计算简单，直观反映形状差异|对垂直和切向分离同等对待，可能不准确|快速初步评估|
|基于地面真值的点到线或点到面距离|考虑了垂直和切向分离的不同影响，更准确|计算相对复杂|需要准确评估形状表示时|
|不依赖地面真值的特征分析|可根据模型内部结构选择模式，无需地面真值|缺乏直观的形状比较|地面真值数据难以获取时|
|不依赖地面真值的交叉验证|综合评估模型性能，减少数据划分的影响|计算量较大|需要全面评估模型泛化能力时|