DOI: 10.12677/aam.2025.144136, PDF, HTML, XML,  被引量   
作者: 吴霖磊：浙江师范大学数学科学学院，浙江 金华
关键词: 保拓扑结构；非刚性；三维医学图像配准；变分模型；Topology-Preserving； Non-Rigid； 3D Image Registration； Variation Model

文章引用：吴霖磊. 保拓扑结构的非刚性医学图像配准方法[J]. 应用数学进展, 2025, 14(4): 21-32. https://doi.org/10.12677/aam.2025.144136

1. 引言

1.1. 研究背景及意义

图像配准，也称作图像匹配、图像融合，在图像处理领域发挥着重要的作用，已经成为图像处理领域的重要研究课题。图像配准的目的是找到一个最优的几何变换来对齐相应的图像数据，这些图像对可以来自于不同的成像时间，不同的成像机制，或者不同的成像角度。通过这种技术，我们可以得到两幅或多组图像的对应信息，比较分析这些目标图像，得到我们想要的图像信息。目前，图像配准在遥感图像处理、计算机视觉、机器人导航等多个领域都有非常广泛的研究与应用。

随着数字技术的发展，数字图像技术愈发成熟，医学影像在疾病诊断、疾病预防、临床手术和术后恢复等方面发挥着越来越重要的作用。与此同时，随着医学成像设备的快速发展，使得医学影像能保存更加丰富与复杂的患者原始信息，这不仅能帮助医生更好地诊断，同时也为我们分析医学影像数据提供了更多的信息与可能性。但是在医学诊疗过程中，为了得到患者更加全面和准确的医学信息，医生往往会同时使用多种医学设备。这就使得医生需要综合这些不同的医学影像来得到所需的信息，因此，图像配准在医学图像处理领域至关重要。无论是不同患者的医学影像对齐，还是同一患者的不同设备医学影像，亦或是不同时期的医学影像对齐，都需要使用医学图像配准的技术。

对于微小的形变或者简单的几何变换模型(比如线性)的形变，已经有许多已知和成熟的方法[1] [2]，但是在当前时代下，面对现今各个领域更加复杂、更具挑战性的问题，这些简单的配准方法已经无法满足我们的需求。特别是在医学图像领域[3] [4]，如在临床应用中，由于器官形变的复杂性和高度非线性，普通的刚体配准(线性配准)已经很难达到应用要求。因此，非刚性图像配准成为了一个热点研究方向。

非刚性图像配准的目的是在不同的图像之间寻找一个具有局部变化变形的最优映射。它在图像处理的许多应用中起着最基础的作用，特别是现代医学图像分析。对于非刚性医学图像配准，由于人体组织中广泛存在的非刚体性质，我们需要寻找可能来自不同成像模态的两个给定图像的非刚性映射。特别地，由于目标器官的现实结构，我们需要保证配准过程保持目标器官的拓扑结构不发生变化，即在形变过程中，目标器官不会发生空间折叠。基于上述事实，目前国内外学者提出的医学图像配准模型都会考虑到以下几个方面的问题：(1) 配准之后的图像需要与参考图像尽可能相似，这部分需要根据图像的模态选择合适的相似性测度，从而对齐图像。(2) 基于不同的实际需求及先验信息，需要选择合适的正则项，如光滑正则项，弹性正则项等。在进行医学图像配准时，根据不同的应用场景，形变需要满足的性质则不同。(3) 所得的图像配准模型必须保证得到的形变是一一对应的，也即不发生空间折叠。否则会破坏器官的解剖结构，反而加大医生的诊疗难度，更可能导致医生的诊疗失误。(4) 综合上述三点，需要根据不同的实际任务，建立合适的图像配准模型，与此同时，图像配准特别是医学图像配准需要一定的计算效率，因此需要选择适合的优化方法，高速、准确地求解所提出的图像配准模型。

综上所述，医学图像配准在医学影像处理领域占据着重要的地位，且在临床等应用中，需要我们提出更加复杂的保持拓扑结构的非刚性医学图像配准模型，以更好地帮助医生分析医学影像数据，更精准地进行医学诊疗。

1.2. 国内外研究现状分析

对于非刚性医学图像配准，国内外研究学者基于不同的优化方法已经提出了许多解决方案，然而在医学影像临床应用中仍是具有挑战性的研究领域，特别是需要在医学图像配准中加入保拓扑结构的先验知识。由于微分同胚映射能保证映射的光滑映射和一一映射，故可以根据微分同胚建立合适的非刚性医学图像配准模型，即在理论上保拓扑结构的非刚性配准[5]-[7]。在实际处理中，为了在数值上求解这种保拓扑结构映射，就要保证得到的映射是一对一可逆映射，即使产生的映射形变，其雅各比行列式为正。目前，已经成为数值上解决保拓扑结构映射的主要方法。

在医学图像配准中，给定一个固定的模板图片 $T:\Omega \subset R^n\to R$ 和一个运动图片 $M:\Omega \subset R^n\to R$ ，其中 $\Omega$ 是 $R^n$ 上的一些图像域，通过形变场 $y\left(x\right):\Omega \subset R^n\to R^n$ 将运动图像配准到模板图像的变分模型可以概括为：

$\underset{y}{\min }L\left(y\right):=S\left(M\circ y,T\right)+\alpha R\left(y\right)$

其中 $S\left(M\circ y,T\right)$ 是评估形变后的运动图像 $M\circ y$ 和固定模板图像T之间差异的保真度度量，正则项 $R\left(y\right)$ 则强迫优化形变场满足特定的约束条件，比如形变场的光滑性与保持拓扑结构。保真项 $S\left(M\circ y,T\right)$ 的选择通常取决于M和T的模态，同模态的问题广泛采用平方差的和的保真项(SSD)：

$S\left(M\circ y,T\right)=\frac{1}{2}{\Vert M\circ y-T\Vert }^2:=\frac{1}{2}{\displaystyle {\int }_{\Omega }{\left(M\left(y\left(x\right)\right)-T\left(x\right)\right)}^2\text{d}x}$

而当图像为不同成像模态时，一方面，可以定义更为复杂的保真项来衡量两幅图像之间的差异，如归一化互相关(NCC)、交互信息(MI)以及比较 $M\circ y$ 和T的归一化梯度场等[8] [9]；另一方面，也可以将不同模态图像投影到统一图像特征空间然后加以比较。为了使提出的优化模型在数学上有良好的表现，在正则项上给形变增加额外的先验信息，许多研究工作已经提出了在不同变换先验信息上的正则项，如弹性正则项[10]，基于粘度的正则项[11]，光滑正则项[12] [13]，以及更复杂的保持拓扑正则项，该正则项明确保证y的雅各比行列式大于0，即 $\det \nabla y>0$ [7] [14]。

最近的研究表明[15]-[17]，拟共形映射理论为拓扑保持映射提供了另一个适当的数学描述。拟共形映射表明，如果Beltrami系数 $\mu$ 满足 $\mu \left(y\left(x\right)\right)<1$ ，则对应的映射是微分同胚的，即 $\left|\mu \left(y\right)\right|<1\iff \det \nabla y>0$ 。这意味着我们可以将具有挑战性的硬性约束条件 $\det \nabla y>0$ ，通过施加 $\mu \left(y\right)<1$ 约束来隐式保持，这种形式的约束条件在数值实验中更容易实现。但是，Beltrami系数只定义在二维的复空间中，而我们研究的多数医学影像为三维或更高维。最新的研究已经提出多种将经典Beltrami系数从二维扩展到n维的方法，并提出了相应的正则项[15] [18]。

Zhang和Chen在三维空间提出了一种新的Beltrami系测度 $N\left(y\right)$ [18]，并推广到一个n维拟共形映射f上，但是需要考虑n的奇偶性，且在2维情形与传统Beltrami系数并不完全等价。为此，Huang和Chen等人直接从2维Beltrami系数进行拓展，将传统的Beltrami系数从2维推广到n维，并提出了新的保拓非刚性医学图像配准模型：

$\underset{y}{\min }L\left(y\right):=D\left(M\circ y,T\right)+\frac{{\alpha }_1}{2}{\displaystyle {\int }_{\Omega }{\Vert \nabla \left(y-x\right)\Vert }_F^2\text{d}x}+\frac{{\alpha }_2}{2}{\displaystyle {\int }_{\Omega }{\Vert {\nabla }^2\left(y-x\right)\Vert }_F^2\text{d}x}+\beta {\displaystyle {\int }_{\Omega }\Phi \left(N\left(x\right)\right)\text{d}x}$

其中新的测度 $N\left(f\right)$ 定义为：

$N\left(f\right)=\frac{{\Vert \nabla f\Vert }_F^2-n\text{sign}\left(\det \nabla f\right){\left(\det \nabla f\right)}^{2/n}}{{\Vert \nabla f\Vert }_F^2+n\text{sign}\left(\det \nabla f\right){\left(\det \nabla f\right)}^{2/n}}$

且惩罚函数选取 $\Phi \left(v\right)=v^2/\left({\left(v-1\right)}^2+1\right)$ 。该模型能很好地保持目标图像的拓扑结构并能达到较好的配准效果。

但是Huang和Chen等人提出的类Beltrami系数正则项的惩罚函数的选取并不是最优的，在优化过程中只能迫使N尽可能趋近于0。而事实上，我们需要约束 $N<1$ 才能保证拓扑结构，但从图像上来看，只有 $\Phi \left(v\right)<1$ 时， $N<1$ 才会成立。

1.3. 本文的研究内容及创新点

针对Huang和Chen等人提出的基于类Beltrami系数的保拓扑结构非刚性医学图像配准模型，我们改进了其类Beltrami系数正则项的惩罚函数，并提出了新的保拓扑结构的非刚性医学图像配准模型，并分析证明了该数学模型解的存在性。同时，我们使用多尺度金字塔的思想方法，运用高斯牛顿法迭代求解所提出的模型。最后，我们在多个不同的三维医学图像配准任务上进行了大量的实验，实验结果表明我们提出的保拓扑结构的非刚性医学图像配准模型在保证拓扑结构的同时，能得到更好的配准效果。

2. 理论基础

2.1. 经典Beltrami方程

直观上讲，拟共形映射，也称为广义微分同胚，将一个域上的无穷小圆映射到目标域上的无穷小椭圆，它代表一个保持定向的同胚映射，具有有界共形失真。

定义2.1 在二维空间中，存在映射 $f:D\to D^{\prime }$ ，其中D和 $D^{\prime }$ 表示复空间C上的两个域，则f称为拟共性映射，如果f满足经典的Beltrami方程：

$\frac{\partial f}{\partial \bar{z}}=\mu \left(z\right)\frac{\partial f}{\partial z},$

其中复值勒贝格可测函数 $\mu \left(z\right)$ 满足 ${\left|\mu \right|}_{\infty }<1$ ，并称之为Beltrami系数。根据拟共形映射的理论，当Beltrami系数 $\mu \left(z\right)=0$ 时，此时映射f则是局部共形的，反之，映射f则是不共形的。因此，我们可以使用Beltrami系数来评价映射f的局部共形性，通过局部无穷小椭圆的离心率和方向。

可以看出，在上述Beltrami方程中，复值 $\bar{z}$ 在二维复空间C中正好与z垂直。因此，由f产生的共形失真可以用Beltrami系数 $\mu \left(z\right)$ 来表示，其中最大放大角是 $\arg \left(\mu \right)/2$ ，相应的放大系数为 $1+\left|\mu \right|$ ，而最大收缩角度是 $\left(\arg \left(\mu \right)-\pi \right)/2$ ，相应的收缩系数为 $1-\left|\mu \right|$ 。因此，可以通过以下公式计算局部失真：

$K_d\left(f\right):=\frac{1+\left|\mu \right|}{1-\left|\mu \right|},$

这基本上量化了在拟共形映射下，一个无穷小球体形变为一个无穷小椭球体的扭曲程度。

因此，Beltrami系数提供了关于映射性质的所有信息。另一方面，如果给定一个Beltrami系数 $\mu \left(z\right):C\to C$ ，且满足条件 ${\left|\mu \right|}_{\infty }\le 1$ ，那么在复空间上总存在一个拟共形映射，它在分布意义上满足Beltrami方程。

2.2. n维拟共形映射

从其他角度来看，容易发现二维映射f的最大放大系数和最大收缩系数正是由 $\nabla f^{\text{T}}\nabla f$ 的最大特征值和最小特征值分别给出的，且映射f是共形的当且仅当 $\nabla f^{\text{T}}\nabla f$ 的这两个特征值相等，即 $K_d\left(f\right)=1$ 。这启发了Lee等人将这种拟共形度量 $K_d\left(f\right)$ 从二维空间扩展到n维空间。

如果f是一个n维的拟共形映射， ${\lambda }_1,\cdots ,{\lambda }_n$ 是矩阵 $\nabla f^{\text{T}}\nabla f$ 的特征值，那么有

${\left|\nabla f\right|}_F^2=\text{Tr}\left(\nabla f^{\text{T}}\nabla f\right)={\lambda }_1+\cdots +{\lambda }_n,$

${\left(\det \nabla f\right)}^2=\det \left(\nabla f^{\text{T}}\nabla f\right)={\lambda }_1\cdots {\lambda }_n.$

且有均值不等式：

${\left({\lambda }_1\cdots {\lambda }_n\right)}^{\frac{1}{n}}\le \left({\lambda }_1+\cdots +{\lambda }_n\right)/n,$

通过观察这一点，LLL模型通过计算 $\nabla f^{\text{T}}\nabla f$ 特征值的算术平均数与几何平均数的比值来评估一个n维映射f的局部失真，即

$K\left(f\right)=\{\begin{array}{l}\frac{1}{n}\frac{{\Vert \nabla f\Vert }_F^2}{{\left(\det \nabla f\right)}^{2/n}}\text{if}\det \nabla f>0\\ +\infty \text{otherwise}\end{array}$

显然，根据均值不等式，这个拟共形映射f的度量 $K\left(f\right)$ 大于或等于1，即 $K\left(f\right)\ge 1$ ，且当 $K\left(f\right)=1$ 时映射f成为共形映射，当且仅当 $\nabla f^{\text{T}}\nabla f$ 的特征值全部相等。因此， $K\left(f\right)$ 的值恰当地描述了f与共形映射之间的距离。

在基于地标的三维配准模型(LLL)中，他们提出了基于测度的新的正则项：

$\underset{y}{\min }\frac{{\left|\nabla y\right|}_F^2}{3{\left(\det \nabla y\right)}^{2/3}}+\frac{\alpha }{2}{\left|\text{Δ}y\right|}_2^2\text{s}\text{.t}.\det y>0,y\left(p_i\right)=q_i,i=1,\cdots ,m$

其中 $y\left(x\right)$ 是一个三维映射， $p_i$ 和 $q_i$ 是对应的地标点。然而，当我们求解这个配准数学模型时，我们仍然需要在每一步保证具有挑战性的硬性约束条件： $\det y>0$ 。

2.3. n维类Beltrami系数方法

观察Beltrami系数 $\mu$ 和它的失真度量K之间的关系，可以得到

$\left|\mu \right|=\frac{K-1}{K+1}.$

显然，对任意的 $K\ge 1$ ，在 $\left|\mu \right|$ 和K之间有一个一对一的映射，且当 $K\left(f\right)=1$ 时 $\mu =0$ ，此时映射f成为共形映射。这意味着在n维空间中， $\left|\mu \right|$ 和K都可以用来量化拟共形映射的局部形变。受到这个想法的启发，ZC等人提出了类似Beltrami系数的度量 $N_1\left(y\right)$ (ZCO1)，以正则化一个三维映射y，使得：

$N_1\left(y\right)=\frac{{\left|\nabla y\right|}_F^2-3{\left(\det \nabla y\right)}^{2/3}}{{\left|\nabla y\right|}_F^2+3{\left(\det \nabla y\right)}^{2/3}},$

这仅仅是模仿了 $\mu$ 。为了进行比较，ZC等人同时使用了另一个度量 $N_2\left(y\right)=K\left(y\right)$ (ZCO2)。

然而，对于这两个度量，我们不能通过直接控制它们的值以保持变换的拓扑性质，我们仍然需要在每一步都满足硬性约束条件 $\det \nabla \left(y\right)>0$ 。为了避免这个困难的约束条件，文献[15]提出了另一个与三维连续可微映射 $y\left(x\right)$ 相关的类Beltrami系数的量 $N\left(y\right)$ ，即

$N_2\left(y\right)=K\left(y\right)=\frac{1}{3}\frac{{\left|\nabla y\right|}_F^2}{{\left(\det \nabla y\right)}^{2/3}}.$

作为y的局部共形失真度量，并证明了它的关键数学性质：

$N\left(y\right)<1\iff \det \nabla y>0,$

类似于二维Beltrami系数。对于新的度量ZC，我们只需要隐式地保持约束 $N\left(y\right)<1$ ，而不需要满足 $\det \nabla y>0$ 。

三维共形失真度量也可以推广到n维拟共形映射f，如下所示：

$N\left(f\right)=\{\begin{array}{l}\frac{{\left|\nabla f\right|}_F-\sqrt{n}{\left(\det \nabla f\right)}^{\frac{1}{n}}}{{\left|\nabla f\right|}_F+\sqrt{n}{\left(\det \nabla f\right)}^{\frac{1}{n}}}\text{if}n\text{is}\text{odd}\\ \frac{{\left|\nabla f\right|}_F-\sqrt{n}{\left(\det \nabla f\right)}^{\frac{1}{n}}}{{\left|\nabla f\right|}_F+\sqrt{n}{\left(\det \nabla f\right)}^{\frac{1}{n}}}\text{if}n\text{is}\text{even}\text{and}\det \nabla f>0\\ \infty \text{otherwise}\end{array}$

然而，由于需要考虑n的奇偶性，这种n维类Beltrami系数的度量的形式相对复杂，因为对于偶数n要求 $\det \nabla y>0$ 。同时，新的度量与当 $n=2$ 时的经典Beltrami系数 $\left|\mu \right|$ 不一致。

为了解决上述ZC模型的问题，HC等人提出了一种新的更通用的类Beltrami系数的度量方法：

$N\left(f\right)=\frac{{\Vert \nabla f\Vert }_F^2-n\text{sign}\left(\det \nabla f\right){\left(\det \nabla f\right)}^{2/n}}{{\Vert \nabla f\Vert }_F^2+n\text{sign}\left(\det \nabla f\right){\left(\det \nabla f\right)}^{2/n}}.$

与ZC类Beltrami系数度量相比，这个方法有两个巨大的优势：首先，该方法不需要考虑n的奇偶性，具有更为统一的形式；其次，当 $n=2$ 时，该方法与经典Beltrami系数相一致。与此同时，他们证明了这个度量的关键数学性质：

$N\left(y\right)<1\iff \det \nabla y>0,$

这可以看作是二维Beltrami系数在n维空间的自然推广。

基于新的类Beltrami系数度量，他们提出了以下HC图像配准模型：

$\underset{y}{\min }L\left(y\right):D\left(M\circ y,T\right)+\frac{{\alpha }_1}{2}{\displaystyle {\int }_{\Omega }{\Vert \nabla \left(y-x\right)\Vert }_F^2\text{d}x}+\frac{{\alpha }_2}{2}{\displaystyle {\int }_{\Omega }{\Vert {\nabla }^2\left(y-x\right)\Vert }_F^2\text{d}x}+\beta {\displaystyle {\int }_{\Omega }\Phi \left(N\left(x\right)\right)\text{d}x}$

其中，第一项数据保真项评估运动图像 $M\circ y$ 与模板图像T之间的差异，第二项和第三项保证变换的光滑性，第四项中的惩罚函数 $\Phi \left(v\right)=v^2/\left({\left(v-1\right)}^2+1\right)$ 施加微分同胚条件 $N\left(y\right)<1$ ，即 $\det \nabla y>0$ 。

在上述将二维Beltrami系数自然推广到n维空间的研究中，已取得了令人满意的成果。但是在上述最新提出的HC模型中，为了施加微分同胚(即保持拓扑结构)的条件，需要通过选取合适惩罚函数来实现，观察可以发现，他们选取的惩罚函数并不是最合理的，仍有许多改进空间。

3. 提出的保拓扑结构非刚性医学图像配准模型

在这一章中，我们首先提出了一个新的惩罚函数，使得我们的新的变分模型满足保拓扑结构的条件 $N\left(y\right)<1$ 。然后，基于我们提出的新的惩罚函数，我们进一步提出了一个新的保拓扑结构的图像配准模型并证明了这个模型的解的存在性。

3.1. 提出的新惩罚函数

在C. Huang等人提出的模型中，她们选择了惩罚函数：

$\Phi \left(v\right)=\frac{v^2}{{\left(v-1\right)}^2+1}$

通过最小化惩罚函数 $\Phi \left(N\left(y\right)\right)$ 使得 $N\left(y\right)<1$ ，并得到了不错的保拓扑结构的配准效果。

在我们的工作中，我们提出了一个新的惩罚函数：

$\text{Φ}\left(v\right)=\log \left[1+\max \left(0,v-1\right)\right]$

与HC惩罚函数相比，我们使用了对数类型的惩罚函数。为了使提出的惩罚函数具有上界，我们增加了最小值函数限制：

$\Phi \left(v\right)=\min \left\{2,\log \left[1+\max \left(0,v-1\right)\right]\right\}$

通过最小化我们的惩罚函数，可以迫使Beltrami类系数 $N\left(y\right)<1$ ，且取得更好的保拓扑配准效果。

两个惩罚函数的图像如图1所示。可以发现，在HC惩罚函数中，如果v小于1，则惩罚函数 $\Phi \left(v\right)$ 也会小于1，因此通过最小化惩罚函数 $\Phi \left(v\right)$ ，HC模型可以尽可能地使v小于1。然而，HC惩罚函数并不能明确的保证 $v<1$ 。我们可以观察到，随着v的逐渐增大， $\varphi \left(v\right)$ 也会逐渐减小。因此，在最小化HC模型时可能会产生一些问题。因此，我们构造了一个新的惩罚函数，如图所示。容易发现，当我们最小化惩罚函数 $\Phi \left(v\right)$ 时，我们可以确保v逐渐减小，特别地，当v小于1时，惩罚函数 $\Phi \left(v\right)$ 将会急剧减小为0。因此，相较于HC惩罚函数，我们的惩罚函数能更好地施加微分同胚的条件 $N\left(y\right)<1$ 。

Figure 1. HC penalty function (left), our penalty function (right)

图1. HC惩罚函数(左)，提出的惩罚函数(右)

3.2. 新的保拓扑结构的图像配准模型

基于我们提出来的新的保拓扑结构的惩罚函数，我们提出了新的保拓扑结构的非刚性医学图像配准模型：

$\underset{y}{\min }L\left(y\right):=D\left(M\circ y,T\right)+\frac{{\alpha }_1}{2}{\displaystyle {\int }_{\Omega }{\Vert \nabla \left(y-x\right)\Vert }_F^2\text{d}x}+\frac{{\alpha }_2}{2}{\displaystyle {\int }_{\Omega }{\Vert {\nabla }^2\left(y-x\right)\Vert }_F^2\text{d}x}+\beta {\displaystyle {\int }_{\Omega }\Phi \left(N\left(x\right)\right)\text{d}x}$

其中第一项保真项 $D\left(M\circ y,T\right)$ 用来保证模板图像T和运动图像 $M\circ y$ 的配准相似性，第二项和第三项用来保证形变场y的光滑性，最后一项则用来保证微分同胚的条件 $N\left(y\right)<1$ ，也就是 $\det \nabla \left(y\right)>0$ ，从而保证该模型能保持拓扑结构，通过最小化我们的惩罚函数：

$\Phi \left(v\right)=\min \left\{2,\log \left[1+\max \left(0,v-1\right)\right]\right\}$

在我们的工作中，我们主要关注来自不同模态的三维医学图像的非刚性配准问题，在实际情况与处理中，这个问题仍具有很大的挑战性。在这样的背景下，我们通过比较模板图像T和运动图像 $M\circ y$ 的归一化的梯度场来构建保真项，也就是比较梯度的方向：

$D\left(M\circ y,T\right)={\displaystyle {\int }_{\Omega }{\left|{\nabla }^{*}M\circ y-{\nabla }^{*}T\right|}^2\text{d}x}+{\left(\left|{\nabla }^{*}M\circ y\right|+\left|{\nabla }^{*}T\right|-\left|{\nabla }^{*}M\circ y+{\nabla }^{*}T\right|\right)}^2\text{d}x$

其中 ${\nabla }^{\text{*}}$ 表示归一化梯度算子。

在我们的研究中，我们认为所有给定的图像在 $\Omega$ 上是连续可微且有紧支撑的，而且形变y是二阶可导的。现在，我们分析我们提出的模型在函数空间 $W^{2,2}\left(\Omega \right)$ 上解的存在性。

不失一般性，我们有下述两个事实：

定义3.1 如果 $\Omega \subset R^3$ 是一个开集，假设函数 $f:\Omega \times R^3\times R^{3\times 3}\times R^{3\times 3\times 3}\to \left[0,+\infty \right)$ 满足下列条件：

(i) $f\left(x,\cdot ,\cdot ,\cdot \right)$ 在 $x\in \Omega$ 上几乎处处连续。

(ii) $f\left(x,y,\psi ,\Theta \right)$ 对任意的 $\left(y,\psi ,\Theta \right)\in R^3\times R^{3\times 3}\times R^{3\times 3\times 3}$ 是关于x可测的。

则称函数f是Caratheodory函数。

引理 3.2 如果 $\Omega \subset R^3$ 是一个开集，假设函数 $f:\Omega \times R^3\times R^{3\times 3}\times R^{3\times 3\times 3}\to \left[0,+\infty \right)$ 满足下列条件：

(i) f是Caratheodory函数。

(ii) $f\left(x,y,\psi ,\Theta \right)$ 是关于 $\Theta$ 的拟凸函数。

(iii) $0\le f\left(x,y,\psi ,\Theta \right)\le a\left(x\right)+C\left({\left|y\right|}^2+{\left|\psi \right|}^2+{\left|\Theta \right|}^2\right)$ ，其中 $a\left(x\right)\in L^1\left(\Omega \right),C\ge 0$ 。

那么 $L\left(y\right)={\displaystyle {\int }_{\Omega }f}\left(x,y,\psi ,\Theta \right)\text{d}x$ 是 $W^{2,2}\left(\Omega \right)$ 上的弱下半连续函数(表示为wlsc)。

实际上，我们可以将我们提出的模型重新表示为：

$L\left(y\right)={\displaystyle {\int }_{\Omega }f}\left(x,y,\nabla y,{\nabla }^{2}y\right)\text{d}x,$

其中

$f\left(x,y,\psi ,\Theta \right)=D\left(M\circ y,T\right)+\frac{{\alpha }_1}{2}{\left|\psi -I\right|}^2+\frac{{\alpha }_2}{2}{\left|\Theta \right|}^2+\beta \Phi \left(\frac{|\psi {|}_F^2-3{\text{sign}}_{ϵ}\left(\det \psi \right){\left(\det \psi \right)}^{2/3}}{|\psi {|}_F^2+3{\text{sign}}_{ϵ}\left(\det \psi \right){\left(\det \psi \right)}^{2/3}}\right)$

解空间定义为 $W=\left\{y\in W^{2,2}\left(\Omega \right):y\left(x\right)=x\text{on}\partial \Omega \right\}$ ，然后我们可以得到下述结论。

命题 3.3 如果模板图像T和运动图像M在 $\Omega$ 上具有紧支撑且连续可微，则我们提出的模型目标函数是 $W^{2,2}\left(\Omega \right)$ 上的弱下半连续函数(wlsc)。

证明. 为了证明这个结论，我们需要满足引理的三个条件：

(i) 首先，我们假定模板图像T和运动图像M是连续可微的，因此它们的归一化梯度场也是连续的。所以我们容易得到 $f(\cdot )$ 是Caratheodory函数。

(ii) 其次， $f\left(x,y,\psi ,\Theta \right)$ 是关于 $\Theta$ 的凸函数，因此也是拟凸函数。

(iii) 最后，由于 $\left|{\nabla }^{\text{*}}\left(M\circ y\right)\right|\le 1$ 以及 $\left|{\nabla }^{\text{*}}\left(T\right)\le 1\right|$ ，我们有

${\left(\left|{\nabla }^{*}\left(M\circ y\right)\right|+\left|{\nabla }^{*}\left(T\right)\right|-\left|{\nabla }^{*}\left(M\circ y\right)+{\nabla }^{*}\left(T\right)\right|\right)}^2\le {\left(\left|{\nabla }^{*}\left(M\circ y\right)\right|+\left|{\nabla }^{*}\left(T\right)\right|\right)}^2\le 4,$

且对于任意的v，有 $\Phi \left(v\right)=\min \left(2,\log \left(1+\max \left(0,v-1\right)\right)\right)\le 2$ ，因此可以得到

$\begin{array}{c}f\left(x,y,\psi ,\Theta \right)=\frac{{\alpha }_1}{2}{\left|\psi -I\right|}^2+\frac{{\alpha }_2}{2}{\left|\Theta \right|}^2+{\left|{\nabla }^{*}\left(M\circ y\right)-{\nabla }^{*}\left(T\right)\right|}^2\\ +{\left(\left|{\nabla }^{*}\left(M\circ y\right)\right|+\left|{\nabla }^{*}\left(T\right)\right|-\left|{\nabla }^{*}\left(M\circ y\right)+{\nabla }^{*}\left(T\right)\right|\right)}^2\\ +\beta \Phi \left(\frac{|\psi {|}_F^2-2{\text{sign}}_{ϵ}\left(\det \psi \right){\left(\det \psi \right)}^{2/n}}{|\psi {|}_F^2+2{\text{sign}}_{ϵ}\left(\det \psi \right){\left(\det \psi \right)}^{2/n}}\right)\\ \le \frac{{\alpha }_1}{2}{\left|\psi \right|}^2+\frac{{\alpha }_2}{2}{\left|\Theta \right|}^2+2\beta +2+4\\ \le \frac{\alpha }{2}\left({\left|y\right|}^2+{\left|\psi \right|}^2+{\left|\Theta \right|}^2\right)+2\beta +6\end{array}$

其中 $\alpha =\max \left({\alpha }_1,{\alpha }_2\right)$ 。因此，函数 $f(\cdot )$ 满足 $0\le f$ 以及 $f\le \alpha +C\left({\left|y\right|}^2+{\left|\psi \right|}^2+{\left|\Theta \right|}^2\right)$ ，也意味着满足引理的条件(iii)： $a\left(x\right)\equiv a=2\beta +6$ 以及 $C=\alpha /2$ 。

因此，函数 $f(\cdot )$ 是 $W^{2,2}\left(\text{Ω}\right)$ 上的弱下半连续函数。

基于引理3.2和命题3.3，我们可以证明我们提出的模型解的存在性。

定理3.4 如果T和M在 $\Omega$ 上连续可微和紧支撑，则我们提出的模型优化问题在空间W上至少存在一个解，其中 $W=\left\{y\in W^{2,2}\left(\Omega \right):y\left(x\right)=x\text{on}\partial \Omega \right\}$ 。

证明. 因为 $L\left(y\right)$ 有下界0，所以可以找到 $L(\cdot )$ 的一个最小序列 ${\left(y_n\right)}_{n\in N}\subset W$ ，使得

$L\left(y_n\right)\underset{\left[n\to \infty \right]}{\to }{m}_1:=\underset{y\in W}{\inf }L\left(y\right)$

由于 $G{F}_{TM}\left(y\right)$ 的正则性以及广义的庞加莱不等式，存在常数 $C_1,C_2\in R$ 使得

$L\left(y\right)\ge C_1{\Vert y\Vert }_{W^{2,2}}^2+C_2$

又因为 $L\left(Id\right)=G{F}_{TM}\left(Id\right)$ 是有限的，我们可以假设 ${\left(L\left(y_n\right)\right)}_{n\in N}$ 以一个常数 $m_2>0$ 为界，确保序列 ${\left(y_n\right)}_{n\in N}$ 在W中一致有界：

$m_2\ge L\left(y_n\right)\ge C_1{\Vert y_n\Vert }_{W^{2,2}}^2+C_2$

$W^{2,2}\left(\Omega \right)$ 是一个自反空间，因此存在一个序列，记为 ${\left(y_{n_l}\right)}_{n_l\in N}$ 满足在 $W^{2,2}\left(\Omega \right)$ 空间， $y_{n_l}\underset{n\to \infty }{\to }{y}^{*}$ 弱收敛。应用引理3.3中的函数 $L\left(y\right)$ 的弱下半连续性(wlsc)，我们有：

$\underset{y\in W}{\inf }L\left(fy\right)=\underset{n\to \infty }{\lim }L\left(y_n\right)=\underset{n_l\to \infty }{\lim }L\left(y_{n_l}\right)\ge L\left(y^{*}\right)\ge \underset{y\in W}{\inf }L\left(y\right)$

因此， $y^{\text{*}}\in W$ 。

4. 数值实验

在这一章，我们研究了在具体场景中我们的配准方法的实际效果，既包括了简单的人造三维图像，也包括具有挑战性的临床医学图像。由于本文提出的方法是基于HC方法改进的配准方法，且HC方法在配准方面的精度以及效率优于主流的方法，因此实验以HC方法为主要的对比基准。

关于模型中的参数 ${\alpha }_1,{\alpha }_2$ 以及 $\beta$ 的选择，本文手动选取参数，使两种方法尽可能得到最好的配准效果。

首先，在章节4.1，我们通过两张三维人工图像测试了同模态图像的配准效果。然后，在章节4.2和4.3我们在多模态图像上进行了实验。在章节4.2中，我们使用了40个病人三维大脑核磁图像。在章节4.3中，我们配准了27个病人的多期相腹部CT图像。最后两个实验都是极具挑战性的多模态场景，能够很好地展示我们模型的实际配准效果。

4.1. 三维人工图像配准实验

首先，我们人工生成了三维图像，并在三维人工图像上进行了配准实验，并与HC方法进行了比较。显然，我们的方法与HC方法相比刚性配准都取得了更优的配准效果。而且相比于HC方法，我们的方法的配准表现有轻微的领先。在表1中我们展示了详细的结果。可以发现，在DSC，Re_SSD和NCC中，我们的方法均优于HC方法，在时间方面则与HC方法相当。表中的指标： $\min \det \nabla$ 表示计算得到的形变场中的最小的雅各比行列式，可以观察到我们和HC的方法求得的最小雅各比行列式均大于0，也就意味着我们的方法与HC方法一样，都能保持形变后的拓扑结构不发生变化。

Table 1. Experiments result on synthetic image registration

表1. 人工合成图像配准实验结果

	DSC	Re_SSD	NCC	$\min \det \nabla$	$\max \det \nabla$	Times
提出的方法	0.918	1.87%	0.899	0.183	2.868	37
HC方法	0.914	2.03%	0.892	0.206	2.908	36

4.2. 三维大脑核磁配准实验

在这个实验中，我们使用了40个来自不同项目的三维大脑核磁图像进行配准实验。具体来说，我们使用第一张图像作为固定图像，并使用剩下的39张图像作为运动图像。因此，总共进行了39次配准实验。

由于这40个三维大脑核磁图像来自不同的病人，因此不同于人工图像配准实验，这个实验是一个具有挑战性的多模态图像配准实验。正如表2中展示的结果，我们可以发现我们的方法仍然能胜任该实验任务，并在配准性能表现上优于HC方法。我们的方法在DSC和NCC中依然优于HC方法，在时间上则与HC方法相当。与此同时，我们的方法和HC方法的 $\min \det \nabla$ 均大于0，也就意味着仍然保持了拓扑结构。

Table 2. Experiments result on 3D brain MRI registration

表2. 三维大脑核磁配准实验结果

	DSC	NCC	$\min \det \nabla$	$\max \det \nabla$	Times
提出的方法	0.616 ± 0.039	0.958 ± 0.008	0.466 ± 0.061	1.977 ± 0.175	163 ± 44
HC方法	0.614 ± 0.040	0.957 ± 0.012	0.438 ± 0.070	1.981 ± 0.189	162 ± 44

4.3. 多期相腹部CT图像配准实验

在本研究中，我们对从不同扫描阶段获得的三维多期腹部CT图像进行了多模态图像配准实验。在图2中，我们展示了腹部CT的四个不同时期的图像。具体来说，我们使用了27个患者的数据，对于每一个患者，我们将动脉期图像作为固定图像，并将剩余三个时期的图像配准到动脉期图像上，因此，我们共进行了81次配准实验。特别的，在本次实验中，关于DSC，我们选择肝脏病变区域作为我们的感兴趣区域。

由于造影剂的注射，不同时期的CT图像会呈现不同的灰度特征，导致不同时期的图像像素之间缺乏强度一致性，并且成像过程中的生理运动(如呼吸)会导致不同时期的图像发生扭曲，因此多期腹部CT图像配准任务是一项十分具有挑战性的多模态医学图像配准任务。详细结果如表3所示，可以发现，我们的方法在DSC和NCC方面仍然领先于HC方法，但是在时间上有轻微滞后。当然，我们的方法和HC方法依然保证了拓扑结构，因为 $\min \det \nabla$ 都大于0。

Figure 2. Four-phase abdomen CT image

图2. 四期腹部CT图像

Table 3. Experiments result on multi-phase abdomen CT image registration

表3. 多期相CT图像配准的实验结果

	DSC	NCC	$\min \det \nabla$	$\max \det \nabla$	Times
提出的方法	0.870 ± 0.068	0.964 ± 0.017	0.089 ± 0.061	4.158 ± 1.233	263 ± 74
HC方法	0.868 ± 0.069	0.960 ± 0.021	0.065 ± 0.053	3.168 ± 0.905	241 ± 99

5. 总结与展望

在本文研究中，基于HC方法，我们主要改进了可以保持拓扑结构的三维医学图像配准模型。结果显示，改进后的模型在几乎相同的时间内取得了更好的性能。实际上，在保拓扑结构的图像配准实验中，与其他几种常见的主流方法相比，HC方法在准确性和效率方面已经取得了非常好的成果。

我们认为，HC配准模型的惩罚函数并不是最合理的选择。正如我们在3.1节中所讨论，HC惩罚函数背后存在着一些问题。相比之下，我们提出的惩罚函数能更有效地保持拓扑结构。在所有实验中，我们提出的基于新的惩罚函数的配准模型在配准精度上都达到了更好的效果，在收敛速度上都达到了近乎相同的时间。在简单的场景中，如三维人工图像配准实验(章节4.1)，我们的方法和HC方法在精度和效率上都能得到令人满意的结果。在实际更复杂和困难的现实场景中，相比于HC方法我们的方法仍然取得了最好的表现(章节4.2，章节4.3)。特别是在章节4.3中，我们选择肝脏部分的病灶作为我们的感兴趣区域，这个区域的研究在临床诊断中非常有意义。因此在章节4.3中，我们使用病灶区域来评估配准效果。正如在表中展示的结果，病灶区域的DSC结果意味着我们的配准模型在临床诊断中能很好地应用且取得了比HC方法更好的效果。此外，正如所有实验所示，我们的方法可以有效地保持图像的拓扑结构。

最后，我们提出的惩罚函数仍然不是最优的选择。在将来的研究中，我们将探索惩罚函数的更多可能性。另一方面，我们将会继续探索能保持拓扑结构的图像配准模型。

参考文献

[1]	Malik, A.S., Boyko, O., Aktar, N. and Young, W.F. (2001) A Comparative Study of MR Imaging Profile of Titanium Pedicle Screws. Acta Radiologica, 42, 291-293. https://doi.org/10.1034/j.1600-0455.2001.042002291.x
[2]	Modersitzki, J. (2009) FAIR: Flexible Algorithms for Image Registration. Society for Industrial and Applied Mathematics.
[3]	Sotiras, A., Davatzikos, C. and Paragios, N. (2013) Deformable Medical Image Registration: A Survey. IEEE Transactions on Medical Imaging, 32, 1153-1190. https://doi.org/10.1109/tmi.2013.2265603
[4]	Chen, K., Lui, L.M. and Modersitzki, J. (2019) Image and Surface Registration. In: Handbook of Numerical Analysis, Vol. 20, Elsevier, 579-611. https://doi.org/10.1016/bs.hna.2019.07.001
[5]	Beg, M.F., Miller, M.I., Trouvé, A. and Younes, L. (2005) Computing Large Deformation Metric Mappings via Geodesic Flows of Diffeomorphisms. International Journal of Computer Vision, 61, 139-157. https://doi.org/10.1023/b:visi.0000043755.93987.aa
[6]	Droske, M. and Rumpf, M. (2004) A Variational Approach to Nonrigid Morphological Image Registration. SIAM Journal on Applied Mathematics, 64, 668-687. https://doi.org/10.1137/s0036139902419528
[7]	Han, H. and Wang, Z. (2020) A Diffeomorphic Image Registration Model with Fractional-Order Regularization and Cauchy-Riemann Constraint. SIAM Journal on Imaging Sciences, 13, 1240-1271. https://doi.org/10.1137/19m1260621
[8]	Haber, E. and Modersitzki, J. (2006) Intensity Gradient Based Registration and Fusion of Multi-Modal Images. In: Larsen, R., Nielsen, M. and Sporring, J., Eds., Medical Image Computing and Computer-Assisted Intervention—MICCAI 2006. Lecture Notes in Computer Science, Springer, 726-733. https://doi.org/10.1007/11866763_89
[9]	Maes, F., Collignon, A., Vandermeulen, D., Marchal, G. and Suetens, P. (1997) Multimodality Image Registration by Maximization of Mutual Information. IEEE Transactions on Medical Imaging, 16, 187-198. https://doi.org/10.1109/42.563664
[10]	Broit, C. (1981) Optimal Registration of Deformed Images. University of Pennsylvania.
[11]	Christensen, G.E., Rabbitt, R.D. and Miller, M.I. (1996) Deformable Templates Using Large Deformation Kinematics. IEEE Transactions on Image Processing, 5, 1435-1447. https://doi.org/10.1109/83.536892
[12]	Fischer, B. and Modersitzki, J. (2002) Fast Diffusion Registration. Contemporary Mathematics, 313, 117-128.
[13]	Frohn-Schauf, C., Henn, S. and Witsch, K. (2007) Multigrid Based Total Variation Image Registration. Computing and Visualization in Science, 11, 101-113. https://doi.org/10.1007/s00791-007-0060-2
[14]	Zhang, J. and Chen, K. (2015) Variational Image Registration by a Total Fractional-Order Variation Model. Journal of Computational Physics, 293, 442-461. https://doi.org/10.1016/j.jcp.2015.02.021
[15]	Zhang, D. and Chen, K. (2020) 3D Orientation-Preserving Variational Models for Accurate Image Registration. SIAM Journal on Imaging Sciences, 13, 1653-1691. https://doi.org/10.1137/20m1320006
[16]	Lam, K.C. and Lui, L.M. (2014) Landmark-and Intensity-Based Registration with Large Deformations via Quasi-Conformal Maps. SIAM Journal on Imaging Sciences, 7, 2364-2392. https://doi.org/10.1137/130943406
[17]	Lee, Y.T., Lam, K.C. and Lui, L.M. (2015) Landmark-Matching Transformation with Large Deformation via N-Dimensional Quasi-Conformal Maps. Journal of Scientific Computing, 67, 926-954. https://doi.org/10.1007/s10915-015-0113-5
[18]	Huang, C., Chen, K., Huang, M., Kong, D. and Yuan, J. (2024) Topology-Preserving Image Registration with Novel Multi-Dimensional Beltrami Regularization. Applied Mathematical Modelling, 125, 539-556. https://doi.org/10.1016/j.apm.2023.09.033