超级会员免费看
非参数正则化的示例实现
1. 避免奇异性的重网格化方法
在处理连续变换时，尽管连续变换本身可能是同胚的，但经过有限采样并使用分段平面函数插值后的版本可能并非如此。对于基于流体的变换，同样会出现类似问题，因为我们使用的是对实际连续流体演化方程的有限采样近似。
基于流体的变换是在规则网格上定义的，因此可以很容易地计算变换在每个点的雅可比矩阵。监测该雅可比矩阵行列式的值，为预测和避免变换中的奇异性（即折叠）提供了一种便捷的方法。当雅可比矩阵的行列式在任何点为零或小于零时，就会出现奇异性。为避免这种情况，我们可以只允许雅可比矩阵在每个点都大于某个预定义值的变换。然而，这种方法只能产生有限的变换集合，在许多实际应用中，这不足以建立对应关系。
一种替代方法是重网格化方法。我们不考虑单个可能在变形较大时变得奇异的变形场，而是将变换分解为一系列较小的、非奇异的变换。当变换在任何点接近奇异性时，就会执行重网格化。具体过程如下：
1. 保存当前变换并将其应用于形状图像，创建一个新的传播形状图像。
2. 将工作变换重置为恒等变换，并使用新的形状图像重新开始优化。
3. 优化完成后，将这些变换组合成所需的单个变换。
2. 非参数正则化的高效算法
接下来将介绍一种使用基于流体的正则化器来优化形状图像训练集之间对应关系的高效算法。如果形状图像具有足够的分辨率，该算法可以单独使用，但在许多计算机系统上，这并不实际。因此，我们采用迭代方法。
2.1 实现选项
- 目标函数 ：基于行列式的目标函数（公式(4.9)），正则化常数 Δ = 0.001。
-
迭代次数
：
niterations设置为 25。 - 流体方程 ：使用有限差分方案和 LU 分解来求解流体方程。
- 姿态 ：不优化姿态变换。
对于形状图像，如果表示的是开放表面（有单个边界边缘），则需要在形状图像上设置边界条件，以防止点移出边缘。对于封闭表面，使用八面体映射，只需固定形状图像的四个奇异点。为简化实现，我们选择固定整个边界，这样该算法就可以在迭代方案中用于开放或封闭表面。
2.2 算法流程
以下是使用流体正则化进行优化的算法：
procedure {Reparameterisationi} ←optimise fluid
({Shape Imagei}, Grid Points)
Description
Uses fluid regularization to find the set of re-parameterisation functions that minimises a group-
wise objective function.
Variables
•
Shape Imagei is a n × 3 shape image representing the ith training shape; the pixel
coordinates are stored in the rows of a n × 2 matrix Grid Points;
•
{Reparameterisationi} is a set of n×2 matrices that represent the re-parameterisation
functions that can be applied to the corresponding shape images to produce the optimal
value of the objective function.
Declarations
•
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;
•
amax = max(a)
returns the element of the vector a with maximum value.
2.3 初始化步骤
-
确定参考形状的索引
ref。 - 将所有形状的重网格化计数器初始化为零：
for i = 1 ... nS
regrid counteri ←0;
- 复制原始形状图像：
for i = 1 ... nS
Original Shape Imagei ←Shape Imagei;
-
将规则网格
Grid Points的坐标连接成一个 2n 维向量：
grid points ←
[
Grid Points(·, 1)T , Grid Points(·, 2)T
]
- 初始化每个形状的位移场向量为零：
for i = 1 ... nS
displacementi ←(0, 0, ..., 0) ;
- 为每个训练示例创建一个 3n 维的形状向量：
for i = 1 ... nS
for j = 1 ... n
sample point ←
[
grid points(j) - displacementi(j),
grid points(j + n) - displacementi(j + n)
]
v ←bilinear interpolate (Grid Points, Shape Imagei, sample point)
shape vectori(j) ←v(1)
shape vectori(j + n) ←v(2)
shape vectori(j + 2n) ←v(3)
-
创建一个 n 维的布尔向量
free node，表示节点是否可以移动。 -
预计算 2n × 2n 的稀疏矩阵
D，用于构建二阶导数集合：
8.1. initialize by setting all elements of D to zero;
8.2. let:
a ←(2μ + λ) / h2
b ←μ / h2
c ←(μ + λ) / (4h2)
d ←(λ + μ) / h2
e ←(λ + 3μ) / (4h2)
where the grid-spacing h = 1/m; fluid viscosity values of μ = 1 and λ = 1 were used here;
8.3. create a 3n-dimensional vector to store the shape function and set it to the mean
shape:
shape function ←
1 / nS * sum(shape vectori)
8.4. populate the matrix using a finite difference scheme:
for k = 1 ... n
if (free node(k))
calculate the gradient of the shape function in the i and j directions and
store them in three-dimensional vectors:
grad i ←
[
shape function(k) - shape function(k + 1),
shape function(k + n) - shape function(k + 1 + n),
shape function(k + 2n) - shape function(k + 1 + 2n)
]
grad j ←
[
shape function(k) - shape function(k + m),
shape function(k + n) - shape function(k + m + n),
shape function(k + 2n) - shape function(k + m + 2n)
]
calculate 2 × 2 the matrix Metric – this corresponds to g in (8.32):
Metric ←1 / h3 *
[
grad i • grad i, grad i • grad j
grad j • grad i, grad j • grad j
]
calculate the inverse of the metric matrix:
Inverse Metric ←Metric−1
fill in the finite difference operations for the i component:
let:
˜a ←a * Inverse Metric(1, 1)
˜b ←b * Inverse Metric(2, 2)
˜c ←c * Inverse Metric(1, 1)
˜d ←d * Inverse Metric(1, 2)
˜e ←e * Inverse Metric(1, 2)
D(k, k) ←−2˜a - 2˜b % vi(i, j)
D(k, k - 1) ←˜a % vi(i + 1, j)
D(k, k + 1) ←˜a % vi(i + 1, j)
D(k, k - m) ←˜b % vi(i, j - 1)
D(k, k + m) ←˜b % vi(i, j + 1)
D(k, k + 1 + m) ←˜e % vi(i + 1, j + 1)
D(k, k - 1 + m) ←−˜e % vi(i - 1, j + 1)
D(k, k + 1 - m) ←−˜e % vi(i + 1, j - 1)
D(k, k - 1 - m) ←˜e % vi(i - 1, j - 1)
D(k, k + n) ←−2 ˜d % vj(i, j)
D(k, k + m + n) ←˜d % vj(i, j + 1)
D(k, k - m + n) ←˜d % vj(i, j - 1)
D(k, k + 1 + m + n) ←˜c % vj(i + 1, j + 1)
D(k, k - 1 + m + n) ←−˜c % vj(i - 1, j + 1)
D(k, k + 1 - m + n) ←−˜c % vj(i + 1, j - 1)
D(k, k - 1 - m + n) ←˜c % vj(i - 1, j - 1)
now do the same for the j component:
let:
˜a ←a * Inverse Metric(2, 2)
˜b ←b * Inverse Metric(1, 1)
˜c ←c * Inverse Metric(2, 2)
˜d ←d * Inverse Metric(2, 1)
˜e ←e * Inverse Metric(2, 1)
D(k + n, k + n) ←−2˜a - 2˜b % vj(i, j)
D(k + n, k - m + n) ←˜a % vj(i, j - 1)
D(k + n, k + m + n) ←˜a % vj(i, j + 1)
D(k + n, k - 1 + n) ←˜b % vj(i + 1, j)
D(k + n, k + 1 + n) ←˜b % vj(i + 1, j)
D(k + n, k + 1 + m + n) ←˜e % vj(i + 1, j + 1)
D(k + n, k - 1 + m + n) ←−˜e % vj(i - 1, j + 1)
D(k + n, k + 1 - m + n) ←−˜e % vj(i + 1, j - 1)
D(k + n, k - 1 - m + n) ←˜e % vj(i - 1, j - 1)
D(k + n, k) ←−2 ˜d % vi(i, j)
D(k + n, k + m) ←˜d % vi(i, j + 1)
D(k + n, k - m) ←˜d % vi(i, j - 1)
D(k + n, k + 1 + m) ←˜c % vi(i + 1, j + 1)
D(k + n, k - 1 + m) ←−˜c % vi(i - 1, j + 1)
D(k + n, k + 1 - m) ←−˜c % vi(i + 1, j - 1)
D(k + n, k - 1 - m) ←˜c % vi(i - 1, j - 1)
8.5. fill in the values for the fixed nodes:
for k = 1 ... n
if NOT(free node(k))
D(k, k) ←1
D(k + n, k + n) ←1
for j = 1 ... 2n
D(k, j) ←0
D(j, k) ←0
D(k + n, j) ←0
D(j, k + n) ←0
8.6. perform an LU decomposition of D to get a lower triangular matrix L and an upper
triangular matrix U.
以下是初始化步骤的流程图：
graph TD;
A[确定参考形状索引 ref] --> B[初始化重网格化计数器为零];
B --> C[复制原始形状图像];
C --> D[连接网格坐标成 2n 维向量];
D --> E[初始化位移场向量为零];
E --> F[创建 3n 维形状向量];
F --> G[创建布尔向量 free node];
G --> H[预计算矩阵 D];
H --> I[LU 分解矩阵 D];
3. 优化步骤
-
进行
niterations次迭代：-
从训练集中随机选择一个形状编号
i（i ≠ ref）。 -
计算目标函数关于第
i个形状样本点位移的梯度gradienti：
-
从训练集中随机选择一个形状编号
gradienti ←
[
∂L / ∂p(i)_1,
...
∂L / ∂p(i)_A,
...
∂L / ∂p(i)_2n
]
where p(i)_A = grid points(A) - displacementi(A)
3. 计算其他各项的分量：
1. 计算平均形状向量。
2. 为每个形状创建形状差异向量。
3. 计算每个样本点连接的所有三角形的面积之和。
4. 计算协方差矩阵 `Covariance`。
5. 获取协方差矩阵的特征向量和特征值。
6. 计算 `∂L / ∂λa`。
7
3. 优化步骤（续）
- 归一化所有特征向量使其具有单位长度。
for a = 1 ... nS - 1
eigenvectora ← eigenvectora / ||eigenvectora||
-
计算
∂λa / ∂Djk。
∂λa / ∂Djk = eigenvectora(j) · eigenvectora(k)
-
计算
δDjk / δSi(xm)的分量：- x 坐标分量：
1 / (nS * sum(int area(m))) * [(nS * δ(i, j) - 1) * centred shape vectork(m)
+ (nS * δ(i, k) - 1) * centred shape vectorj(m)]
- y 和 z 坐标分量通过分别用 `m + n` 和 `m + 2n` 替换 `m` 获得。
-
使用有限差分方案数值近似每个自由网格点的
δSi(x) / δp(i)_A的 x 坐标分量：
perturb the Ath point of the ith shape by a small amount Δ = 10−5 parallel to the x - axis:
re - parameterise by perturbing the x - coordinate of the sample point:
sample point ← [grid point(A) - Δ, grid point(A + n)]T
sample the the perturbed point on the shape:
v ← bilinear interpolation(Grid Points, Shape Imagei, sample point)
create a copy of the current shape vector for the ith shape and replace the Ath point with the perturbed point:
perturbed shape vector ← shape vectori
perturbed shape vector(A) ← v(1)
perturbed shape vector(A + n) ← v(2)
perturbed shape vector(A + 2n) ← v(3)
estimate the derivative as a finite difference:
δSi(x) / δp(i)_A = (perturbed shape vector - shape vectori) / Δ
-
对每个自由网格点的
δSi(x) / δp(i)_A的 y 坐标分量进行同样的操作：
perturb the Ath point of the ith shape by a small amount Δ = 10−5 parallel to the y - axis:
re - parameterise by perturbing the y - coordinate of the sample point:
sample point ← [grid point(A), grid point(A + n) - Δ]T
sample the the perturbed point on the shape:
v ← bilinear interpolation(Grid Points, Shape Imagei, sample point)
create a copy of the current shape vector for the ith shape and replace the Ath point with the perturbed point:
perturbed shape vector ← shape vectori
perturbed shape vector(A) ← v(1)
perturbed shape vector(A + n) ← v(2)
perturbed shape vector(A + 2n) ← v(3)
estimate the derivative as a finite difference:
δSi(x) / δp(i)_(A + n) = (perturbed shape vector - shape vectori) / Δ
-
对于每个形状（
i = 1 ... nS，i ≠ ref）：-
求解线性偏微分方程
-gradienti = D * velocityi以得到第i个形状的速度场velocityi：
-
求解线性偏微分方程
Solve La = -gradienti for a by forward substitution, then solve U * velocityi = a for velocityi by back substitution (where L and U are the pre - computed matrices produced by the LU decomposition computed in the Initialization).
2. 计算第 `i` 个形状的位移场沿 x 和 y 轴的梯度对应的 n 维向量：
for j = 1 ... n
if free node(j)
duxdx(j) ← (displacement(j + m) - displacement(j - m)) / (2h)
duxdy(j) ← (displacement(j + 1) - displacement(j - 1)) / (2h)
duydx(j) ← (displacement(j + m + n) - displacement(j - m + n)) / (2h)
duydy(j) ← (displacement(j + 1 + n) - displacement(j - 1 + n)) / (2h)
else
duxdx(j) ← 0
duxdy(j) ← 0
duydx(j) ← 0
duydy(j) ← 0
3. 计算 n × 2 正则化矩阵 `R`：
for j = 1 ... n
R(j, 1) ← velocity(j) - (velocity(j) * duxdx(j)+ velocity(j + n) * duxdy(j))
R(j, 2) ← velocity(j + n) - (velocity(j) * duydx(j)+ velocity(j + n) * duydy(j))
4. 计算时间步长：
t ← δ / max(||R||)
where the || · || operator returns a vector representing the Euclidian norm of each row of R; a value of δ = 10−3 was used.
5. 计算提议的位移：
for j = 1 : n
proposed displacement(j, 1) ← grid points(j) - displacementi(j) - t * R(j, 1)
proposed displacement(j, 2) ← grid points(j + n) - displacementi(j + n) - t * R(j, 2)
6. 如果 `m * proposed displacement` 的雅可比矩阵在任何点小于 0.5：
- 执行重网格化：
1. 增加第 `i` 个形状的重网格化计数器。
regrid counteri = regrid counteri + 1
2. 将当前变换存储到矩阵中。
Saved Reparameterisationi(regrid counteri, ·) ← grid points - displacementi
3. 根据累积的重新参数化对形状图像进行重采样：
let Cumulative Reparameterisation ← Grid Points;
for k = 1 ... regrid counteri
reshape the saved re - parameterisation into a n × 2 matrix:
for j = 1 ... n
Reshaped Reparameterisation(j, 1) ← Saved Reparameterisationi(k, j)
Reshaped Reparameterisation(j, 2) ← Saved Reparameterisationi(k, j + n)
Reshaped Reparameterisation(j, 3) ← Saved Reparameterisationi(k, j + 2n)
accumulate the re - parameterisations:
for j = 1 ... n
Cumulative Reparameterisation(j, ·) ← bilinear interpolation(Grid Points, Reshaped Reparameterisation, Cumulative Reparameterisation(j, ·));
resample the original shape image to create a new shape image, sampled according to Resized Reparameterisation:
for j = 1 ... n
Shape Imagei(j, ·) ← bilinear interpolation(Grid Points, Original Shape Imagei, Resized Reparameterisation(j, ·));
4. 将第 `i` 个训练示例的位移向量的所有元素设置为零。
displacementi ← (0, 0, ... 0)T
5. 跳转到步骤 2.2。
- 否则，接受位移：
for j = 1 ... n
displacementi(j) ← displacementi(j) + t * Ri(j, 1)
displacementi(j + n) ← displacementi(j + n) + t * Ri(j, 2)
-
更新第
i个形状向量：
for j = 1 ... n
sample point ← [grid points(j) - displacementi(j), grid points(j + n) - displacementi(j + n)]
v ← bilinear interpolate (Grid Points, Shape Imagei, sample point)
shape vectori(j) ← v(1)
shape vectori(j + n) ← v(2)
shape vectori(j + 2n) ← v(3)
以下是优化步骤的流程图：
graph TD;
A[开始迭代] --> B[随机选形状 i];
B --> C[计算梯度 gradienti];
C --> D[计算其他分量];
D --> E[求解线性 PDE 得 velocityi];
E --> F[计算位移场梯度];
F --> G[计算正则化矩阵 R];
G --> H[计算时间步长 t];
H --> I[计算提议位移];
I --> J{雅可比矩阵 < 0.5?};
J -- 是 --> K[执行重网格化];
K --> F;
J -- 否 --> L[接受位移];
L --> M[更新形状向量];
M --> N{迭代结束?};
N -- 否 --> B;
N -- 是 --> O[结束];
4. 后处理步骤
对每个形状（
i = 1 ... nS
）执行以下操作：
1. 保存最终的重新参数化：
regrid counteri ← regrid counteri + 1
Saved Reparameterisationi(regrid counteri, ·) ← grid points - displacementi
- 计算累积的重新参数化：
let Cumulative Reparameterisation ← Grid Points;
for k = 1 ... regrid counteri
reshape the saved re - parameterisation into a n × 2 matrix:
for j = 1 ... n
Reshaped Reparameterisation(j, 1) ← Saved Reparameterisationi(k, j)
Reshaped Reparameterisation(j, 2) ← Saved Reparameterisationi(k, j + n)
Reshaped Reparameterisation(j, 3) ← Saved Reparameterisationi(k, j + 2n)
accumulate the re - parameterisations:
for j = 1 ... n
Cumulative Reparameterisation(j, ·) ← bilinear interpolation(Grid Points, Reshaped Reparameterisation, Cumulative Reparameterisation(j, ·));
-
设置
Reparameterisationi为累积的重新参数化。
Reparameterisationi ← Cumulative Reparameterisationi
最后返回
{Reparameterisationi}
。
整个算法流程涵盖了初始化、优化和后处理三个主要阶段，通过这些步骤可以有效地优化形状图像训练集之间的对应关系。在实际应用中，可以根据具体需求调整算法中的参数，如迭代次数、正则化常数等，以达到更好的效果。
以下是后处理步骤的总结表格：
|步骤|操作内容|
|----|----|
|1|保存最终重新参数化，增加重网格化计数器并存储当前变换|
|2|计算累积重新参数化，包括重塑和累积操作|
|3|设置
Reparameterisationi
为累积的重新参数化|
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
