CN113570724A - 基于逆球极投影的网格球面保角参数化方法及其应用 - Google Patents

Abstract

基于逆球极投影的网格球面保角参数化方法及其应用，属于计算机图像分析领域。该方法主要是针对零亏格网格球面保角参数化计算量大的现状，对网格进行切割，切成两片亏格为一的网格；然后对两片网格分别进行平面参数化，该参数化是固定边界为单位圆的调和映射，把切缝均匀映射到单位圆上，从而把两个亏格为一的网格曲面分别保角参数化到两个单位圆盘；再对一个圆盘进行复反演，使两个圆盘覆盖整个扩充复平面；然后进行球极投影的逆变换，使扩充复平面上的点映射到黎曼球面；最后从莫比乌斯群里找到使调和能量最小且满足0质心约束的莫比乌斯变换，以使参数化后的网格尽可能保持的原始网格特征。将该方法应用于脑部几何向量分析，实现在不采用外科手术的前提下，了解疾病进展的状态，减轻病人的痛苦。

Description

基于逆球极投影的网格球面保角参数化方法及其应用

技术领域

本发明属于计算机图像分析领域，具体涉及一种基于逆球极投影的网格球面保角参数化方法，该方法可广泛应用于医学影像大脑病态分析领域。

背景技术

上个世纪60年代，人类对灵长类动物猴子的解剖发现，从视网膜到第一级视觉中枢的大脑皮层曲面的映射是保角映射，保角变换的最大特点是局部保持形状，但是忽略面积大小。保角映射本质上就是共形映射。近20年来，随着人们对大脑的深入研究，发现众多神经系统疾病都与大脑某个区域形态异常有着直接关系，例如阿兹海默症一般都伴随着海马体与内侧题叶的萎缩。现代医学的传统方法是医生通过CT或MRI等成像技术获取大脑图像，在不给大脑造成伤害的前提下，对大脑进行测量，获取大脑信息。基于这类成像技术获得的医学图像，使得现代医学在脑科学领域不断发展。但是在海量的医学影像中，大部分需要医生凭借个人专业知识与临床经验去分析，难免会使效率降低，甚至出现误诊情况。在脑科学领域，将健康标准的大脑皮层网格模型与待诊断的大脑皮层网格模型运用本发明的参数化方法，映射到同一参数域，再做对比。这将会大大提升诊断效率与准确率。还有癫痫，儿童自闭症等脑神经性疾病，都可以利用微分几何的知识将相关器官共形映射到球面上，进而加以精确分析比对，判断疾病的类型与发展。

所有这些都要用到曲面的球面参数化技术。随着计算机图形学和微分几何的发展，人们对球面参数化也进行了深入研究，总结来说所有方法都是在微分同胚群中进行优化，问题具有本质的难度，第一计算量大，第二计算复杂度大。

本发明就是变形过程中的配准研究，也就是找到源曲面与目标曲面的一一对应关系，涉及在曲面变形过程中的一对一映射问题，以更好的保持曲面的局部特性。基于逆球极投影的球面参数化方法，该方法主要是针对零亏格网格球面保角参数化计算量大的现状，对网格进行切割，切成两片亏格为一的网格；然后对两片网格分别进行平面参数化，该参数化是固定边界为单位圆的调和映射，把切缝均匀映射到单位圆上，从而把两个亏格为一的网格曲面分别保角参数化到两个单位圆盘；再对一个圆盘进行复反演，使两个圆盘覆盖整个扩充复平面；然后进行球极投影的逆变换，使扩充复平面上的点映射到黎曼球面；最后从莫比乌斯群里找到使调和能量最小且满足0质心约束的莫比乌斯变换，以使参数化后的网格尽可能保持的原始网格特征。本发明通过上述方法，提供了一种可以实现网格球面保角参数化的方法。

发明内容

为了解决上述存在的问题，本发明提供一种基于逆球极投影的网格球面保角参数化方法，该方法利用逆球极投影和调和能量相结合，采用简单的方法避免复杂的运算，实现对任意零亏格曲面进行球面参数化。将该方法可广泛应用于脑部几何向量分析。

本发明目的是通过下述技术方案实现的：基于逆球极投影的网格球面保角参数化方法，包括如下步骤：

步骤一：将大脑医学影像数据图形化，获取数据的网格坐标位置和厚度、所属解剖区域等属性信息，进而获得封闭的大脑网格曲面；

步骤二：对步骤一获取的网格曲面进行切割，切割成两片亏格为一的网格曲面C₁、C₂

步骤三：对步骤二形成的两片网格分别进行平面参数化：以网格的固定边界为单位圆的调和映射，把切缝均匀映射到单位圆上，从而把两个亏格为一的网格分别保角参数化到两个单位圆盘c₁，c₂上；

步骤四，对步骤三获得的c₂进行复反演，使c₁，c₂覆盖整个扩充复平面C；

步骤五，针对步骤四获得的扩充复平面C进行球极投影的逆变换，使扩充复平面上的点映射到黎曼球面；

步骤六，针对步骤五获得的黎曼球面进行莫比乌斯变换，完成网格球面参数化；从莫比乌斯群里找到使调和能量最小且满足0质心约束的莫比乌斯变换，以使参数化后的网格尽可能保持的原始网格特征。

所述的步骤二实现过程如下：

(2.1)最远点采样：对于网格中的点，采用最远点采样算法抽取3个相互距离最远的点N₁，N₂，N₃；

(2.2)初始切缝形成：由最短路径算法获取N₁到N₂、N₂到N₃、N₃到N₁之间的最短路径，并将它们联接在一起形成一个闭合环路作为初始切缝；

(2.3)对切缝的处理：确保在切缝上的任意连续的三个点所构成的三角形不属于网格；具体为：当连续的三个点所构成的三角形属于网格，且出现在切缝的联接点N₁，N₂，N₃附近，则在切缝中删除三点中处于中间位置的点，然后重新遍历直至不存在连续的三个点所构成的三角形属于网格这种情况；

(2.4)网格切割：沿切缝对网格进行切割，将网格切割成切割成两片亏格为1的曲面C₁、C₂。

所述的步骤三实现过程如下：

(3.1)固定平面参数化的边界：统计切缝点的数量n，将单位圆盘c₁，c₂的边界即单位圆均分成n段，将切缝点一一映射到边界上；

(3.2)进行Tutte映射：将曲面C₁、C₂；分别进行固定边界的Tutte映射，边界为(3.1)所述的单位圆；

(3.3)计算调和能量

使调和能量最小，实现参数化的过程中保角，其中α、β、f(u)、f(v)参数关系满足图2所示位置关系；

其中，u、v是原网格上一条边的两个顶点，f(u)、f(v)是u、v经过映射后的点，α、β是原网格上由u、v所构成的边的两个对角；

(3.4)保存平面参数化结果{(X_n，Y_n)},其中X_n为由平面参数化结果的横坐标组成的向量，Y_n为由平面参数化结果的纵坐标组成的向量，n＝{1，2}表示两组平面参数化结果。

所述的步骤四实现过程如下：

(4.1)坐标转换：将参数化结果{(X_n，Y_n)}转换到扩充复平面上，记为复数{X_n+iY_n}；其中

表示虚数；

(4.2)复反演：将{X₂+iY₂}做关于单位圆的复反演

表示z的共轭复数；经过复反衍{(X_n，Y_n)}将无重叠地连续地映射到整个扩充复平面；

其中

表示虚数，x+iy表示一个复数，x为实部，y为虚部；X₂，Y₂分别由圆盘c2经参数化后实部和虚部；{X₂+iY₂}为参数化后有序集合。

所述的步骤五实现过程如下：

(5.1)球极投影：在球极投影下，黎曼球面等价于扩充复平面；将两份平面参数化结果无重叠地连续地一一映射到球面上；

(5.2)求取球面坐标：由扩充重复平面与球面的对应关系

求得球面参数化的坐标。

所述的步骤六实现过程如下：

(6.1)莫比乌斯变换:

其中z，a，b，c，d为满足ad-bc≠0的扩充复数；

(6.2)0质心约束：

其中σ_M网格模型M的面元，

是质心；

由于该机算是非线性的，在实际实施中用以下步骤代替：

(6.2.1)计算M的质心

(6.2.2)对于所有点v∈M，

(6.2.3)对于所有点v∈M，

(6.3)计算调和能量，以(6.2)的网格为输入重复步骤二到步骤六，使调和能量最小。

本发明的有益效果：本方法针对脑科学领域中大脑皮层的病态分析以及计算机图形学领域的纹理映射的球面参数化计算量不断增加的现状，对0亏格网格进行切割，切成两片亏格为一的网格；对两片网格分别进行平面参数化，该参数化是固定边界为单位圆的调和映射，把切缝均匀映射到单位圆上，从而把两个亏格为一的网格分别保角参数化到两个单位圆盘c₁，c₂上；对c₂进行复反演，使c₁，c₂覆盖整个扩充复平面C；进行球极投影的逆变换，使扩充复平面上的点映射到黎曼球面；从莫比乌斯群里找到使调和能量最小且满足0质心约束的莫比乌斯变换，以使参数化后的网格尽可能保持的原始网格特征，进而完成曲面的球面参数化。将该方法应用于脑部皮层厚度分析，首先将病变的大脑按本发明参数化方法映射到球面，再利用莫比乌斯变换将大脑球面的解剖区域对准，根据统计理论，分析其皮层厚度或者是其他几何特征量(形态学张量分析，多元形态学张量)，辅助病变状态分析，实现在不采用外科手术的前提下，了解疾病进展的状态，减轻了病人的痛苦。

附图说明

图1是本发明基于逆球极投影的网格球面保角参数化方法的流程图。

图2是本发明能量调和过程中，保角参数关系图。

具体实施方式

步骤一：将大脑的MRI核磁共振影像通过Freesurfer软件预处理后再图形化，获取数据的网格坐标位置和厚度、所属解剖区域等属性信息，从而得到封闭的大脑皮层网格模型。

步骤二：对步骤一获取的网格曲面进行切割，切割成两片亏格为一的网格曲面C₁、C₂；

具体为：

2.1)最远点采样：对于网格中的点，采用最远点采样算法(Farthest PointSampling简称FPS)抽取3个相互距离最远的点N₁，N₂，N₃；

2.2)初始切缝形成：由最短路径算法(Dijkstra)获取N₁到N₂、N₂到N₃、N₃到N₁之间的最短路径，并将它们联接在一起形成一个闭合环路作为初始切缝；

2.3)对切缝的处理：为使参数化的失真最小，要确保在切缝上的任意连续的三个点所构成的三角形不属于网格。存在连续的三个点所构成的三角形属于网格这种情形主要出现在切缝的联接点N₁，N₂，N₃附近。如果出现这种情形，则在切缝中删除三点中处于中间位置的点，然后重新遍历直至不存在这种情形；

2.4)网格切割：沿切缝对网格进行切割，将网格切割成切割成两片亏格为1的曲面C₁、C₂。

步骤三：对步骤二形成的两片网格分别进行平面参数化：以网格的固定边界为单位圆的调和映射，把切缝均匀映射到单位圆上，从而把两个亏格为一的网格分别保角参数化到两个单位圆盘c₁，c₂上；

具体为：

3.1)固定平面参数化的边界：统计切缝点的数量n，将单位圆盘的边界即单位圆均分成n段，将切缝点一一映射到边界上；

3.2)进行Tutte映射：将C₁、C₂；分别进行固定边界的Tutte映射，边界为3.1)所述的单位圆。

3.3)计算调和能量

使调和能量最小，实现参数化的过程中保角，其中α、β、f(u)、f(v)参数关系满足图2所示位置关系；其中，u、v是原网格上一条边的两个顶点，f(u)、f(v)是u、v经过映射后的点，α、β是原网格上由u、v所构成的边的两个对角。

3.4)保存平面参数化结果{(X_n，Y_n)}，其中X_n为由平面参数化结果的横坐标组成的向量，Y_n为由平面参数化结果的纵坐标组成的向量，n＝{1，2}表示两组平面参数化结果。

步骤四，对步骤三获得的c₂进行复反演，使c₁，c₂覆盖整个扩充复平面C；

具体为：

4.1)坐标转换：将参数化结果{(X_n，Y_n)}转换到扩充复平面上，记为复数{X_n+iY_n}；其中

表示虚数；

4.2)复反演：将{X₂+iY₂}做关于单位圆的复反演

表示z的共轭复数；经过复反衍{(X_n，Y_n)}将无重叠地连续地映射到整个扩充复平面；

其中

表示虚数，x+iy表示一个复数，x为实部，y为虚部；X₂，Y₂分别由圆盘c2经参数化后实部和虚部；{X₂+iY₂}为参数化后有序集合。

步骤五，针对步骤四获得的扩充复平面C进行球极投影的逆变换，使扩充复平面上的点映射到黎曼球面；

具体为：

5.1)球极投影：在球极投影下，黎曼球面等价于扩充复平面。由此，将两份平面参数化结果无重叠地连续地一一映射到球面上；

5.2)求取球面坐标：由扩充重复平面与球面的对应关系

求得球面参数化的坐标。

步骤六，针对步骤五获得的黎曼球面进行莫比乌斯变换，完成网格球面参数化；从莫比乌斯群里找到使调和能量最小且满足0质心约束的莫比乌斯变换，以使参数化后的网格尽可能保持的原始网格特征。

具体为：

6.1)莫比乌斯变换:

其中z，a，b，c，d为满足ad-bc≠0的(扩充)复数。

6.2)0质心约束：

其中σ_M网格模型M的面元，

是质心。

由于该机算是非线性的，在实际实施中用以下步骤代替：

6.2.1)计算M的质心

6.2.2)对于所有点v∈M，

6.2.3)对于所有点v∈M，

6.3)计算调和能量，以(6.2)的网格为输入重复步骤二到步骤六，使调和能量最小。

最后计算参数化后曲面的雅克比矩阵和行列式，进而计算厚度差异，和形态学张量差异等。

本发明的应用原理：

在神经系统疾病的诊断中，将健康标准大脑的网格模型与待诊断的大脑网格模型运用本方法参数化到单位球面。分别将两个大脑网格模型进行球面参数化：首先，将网格进行最远点采样采取三个点N₁，N₂，N₃，由短路径算法获取N₁到N₂、N₂到N₃、N₃到N₁之间的最短路径，并将它们联接在一起形成一个闭合环路作为初始切缝，处理切缝确保切缝上的任意连续的三个点所构成的三角形不属于网格；其次，对两片网格分别进行平面参数化，该参数化是固定边界为单位圆的调和映射，把切缝均匀映射到单位圆上，从而把两个亏格为一的网格分别保角参数化到两个单位圆盘c₁，c₂上；接着，将参数化结果{(X_n，Y_n)}转换到扩充复平面上，记为{X_n+iY_n}，对c₂进行复反演，使c₁，c₂覆盖整个扩充复平面c；然后，进行球极投影的逆变换，使扩充复平面上的点映射到黎曼球面由扩充重复平面与球面的对应关系

求得球面参数化的坐标；最后，从莫比乌斯群里找到使调和能量最小且满足0质心约束的莫比乌斯变换，以使参数化后的网格尽可能保持的原始网格特征。至此，网格模型的保角参数化完成。将两组参数化结果做对准，识别待诊断大脑与健康标准大脑的差异，进而进行诊断。

具体应用例：

利用该方法对正常人与盲人脑部几何向量分析(包括皮层厚度分析，形态学张量分析TBM，多元形态学张量分析mTBM)的过程如下：

现有正常人样本大脑核磁共振影像数据18个，先天盲人样本14个，后天盲人样本11个。

将所有人的大脑皮层数据映射到球面上，通过前述参数化方法分别生成正常人标准解剖大脑曲面、先天盲病人解剖大脑曲面、后天盲病人解剖大脑曲面；将先天盲病人解剖大脑曲面、后天盲病人解剖大脑曲面分别与正常人大脑标准曲面做对比，得出他们视觉区域的厚度差异以及几何张量的差异。先天盲人，左半球视觉区域的皮层平均厚度与正常人和后天盲人差异比较大，后天盲人与正常人在厚度上差异不太明显，而且他们的形态学张量TBM以及多元形态学张量mTBM都是有差异的。大脑左半球面视觉区域比较多，右半球面比较少。具体对比数据如下：

本发明还可以应用于阿尔茨海默病症对海马体结构的检测，从而诊断、判别疾病的发展状态和确定治疗方案。以及对抑郁症等脑部疾病的诊断。

以上对本发明实施例所提供的基于逆球极投影的三角网格参数化方法进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只用于帮助理解本发明的方法及其核心思想；同时对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (6)

1.基于逆球极投影的网格球面保角参数化方法，其特征在于，包括如下步骤：

步骤一：将大脑医学影像数据图形化，获取数据的网格坐标位置和厚度、所属解剖区域等属性信息，进而获得封闭的大脑网格曲面；

步骤二：对步骤一获取的网格曲面进行切割，切割成两片亏格为一的网格曲面C₁、C₂

步骤三：对步骤二形成的两片网格分别进行平面参数化：以网格的固定边界为单位圆的调和映射，把切缝均匀映射到单位圆上，从而把两个亏格为一的网格分别保角参数化到两个单位圆盘c₁，c₂上；

步骤四，对步骤三获得的c₂进行复反演，使c₁，c₂覆盖整个扩充复平面C；

步骤五，针对步骤四获得的扩充复平面C进行球极投影的逆变换，使扩充复平面上的点映射到黎曼球面；

步骤六，针对步骤五获得的黎曼球面进行莫比乌斯变换，完成网格球面参数化；从莫比乌斯群里找到使调和能量最小且满足0质心约束的莫比乌斯变换，以使参数化后的网格尽可能保持的原始网格特征。

2.根据权利要求1所述基于逆球极投影的网格球面保角参数化方法，其特征在于，所述的步骤二实现过程如下：

(2.1)最远点采样：对于网格中的点，采用最远点采样算法抽取3个相互距离最远的点N₁，N₂，N₃；

(2.2)初始切缝形成：由最短路径算法获取N₁到N₂、N₂到N₃、N₃到N₁之间的最短路径，并将它们联接在一起形成一个闭合环路作为初始切缝；

(2.3)对切缝的处理：确保在切缝上的任意连续的三个点所构成的三角形不属于网格；具体为：当连续的三个点所构成的三角形属于网格，且出现在切缝的联接点N₁，N₂，N₃附近，则在切缝中删除三点中处于中间位置的点，然后重新遍历直至不存在连续的三个点所构成的三角形属于网格这种情况；

(2.4)网格切割：沿切缝对网格进行切割，将网格切割成切割成两片亏格为1的曲面C₁、C₂。

3.根据权利要求1所述基于逆球极投影的网格球面保角参数化方法，其特征在于，所述的步骤三实现过程如下：

(3.1)固定平面参数化的边界：统计切缝点的数量n，将单位圆盘c₁，c₂的边界即单位圆均分成n段，将切缝点一一映射到边界上；

(3.2)进行Tutte映射：将曲面C₁、C₂；分别进行固定边界的Tutte映射，边界为(3.1)所述的单位圆；

(3.3)计算调和能量

使调和能量最小，实现参数化的过程中保角，其中α、β、f(u)、f(v)参数关系满足图2所示位置关系；

其中，u、v是原网格上一条边的两个顶点，f(u)、f(v)是u、v经过映射后的点，α、β是原网格上由u、v所构成的边的两个对角；

(3.4)保存平面参数化结果{(X_n，Y_n)}，其中X_n为由平面参数化结果的横坐标组成的向量，Y_n为由平面参数化结果的纵坐标组成的向量，n＝{1，2}表示两组平面参数化结果。

4.根据权利要求1所述基于逆球极投影的网格球面保角参数化方法，其特征在于，所述的步骤四实现过程如下：

(4.1)坐标转换：将参数化结果{(X_n，Y_n)}转换到扩充复平面上，记为复数{X_n+iY_n}；其中

表示虚数；

(4.2)复反演：将{X₂+iY₂}做关于单位圆的复反演

表示z的共轭复数；经过复反衍{(X_n，Y_n)}将无重叠地连续地映射到整个扩充复平面；

其中

表示虚数，x+iy表示一个复数，x为实部，y为虚部；X₂，Y₂分别由圆盘c2经参数化后实部和虚部；{X₂+iY₂}为参数化后有序集合。

5.根据权利要求1所述基于逆球极投影的网格球面保角参数化方法，其特征在于，所述的步骤五实现过程如下：

(5.1)球极投影：在球极投影下，黎曼球面等价于扩充复平面；将两份平面参数化结果无重叠地连续地一一映射到球面上；

(5.2)求取球面坐标：由扩充重复平面与球面的对应关系

求得球面参数化的坐标。

6.根据权利要求1所述基于逆球极投影的网格球面保角参数化方法，其特征在于，所述的步骤六实现过程如下：

(6.1)莫比乌斯变换：

其中z，a，b，c，d为满足ad-bc≠0的扩充复数；

(6.2)0质心约束：

其中σ_M网格模型M的面元，

是质心；

由于该机算是非线性的，在实际实施中用以下步骤代替：

(6.2.1)计算M的质心

(6.2.2)对于所有点v∈M，

(6.2.3)对于所有点v∈M，

(6.3)计算调和能量，以(6.2)的网格为输入重复步骤二到步骤六，使调和能量最小。

Images