医学图像融合的几何方法与算法评价：立方体相交、数字几何、医学成像的研究

285理论计算机科学电子笔记46（2001）网址：http://www.elsevier.nl/locate/entcs/volume46.html24页体素空间相交的几何Jean-PierreReveill`es1，2LLAIC法国克莱蒙费朗大学摘要当前众多数字医疗设施和3D技术（正电子发射断层扫描、磁共振成像、同步辐射、雷达、立体成像等）3D物体的形状通常可以通过相同的立方体（或体素）被同化为R3空间的将通过两个不同的过程获得的单个对象的两个这样的视图相关联，以便将它们的信息融合在新的图像上，需要对这两个拼接进行处理，特别是对于其尺寸接近所采用的技术的分辨率的对象关键词：立方体相交，连分数，数字几何，医学成像。1介绍医学图像融合问题是将三维数字图像看作由相同的立方体构成的R3我们提出了一个几何方法来解决这个问题，通过计算两个家庭的重叠体素的交集。两个3D图像的融合问题仍然是医学成像中的一个问题，主要是当这些图像的小关键部分必须匹配时，例如，尾状核和壳核等大脑区域就是病理学家感兴趣的几个类似区域包含在相当小的框中，有时每个方向不超过10个像素，并且由于技术原因，医疗设施的精度是固定的，没有办法放大以获得更好的视图。 只有详细的数学研究 可能对外科医生有帮助1ContratACIno23duMinist'eredelaRechercheetdelaTechnologie2电子邮件地址：reveil@llaic.u-clermont1.fr2001年由ElsevierScienceB出版。 诉 在CCBY-NC-ND许可下开放访问。REVEILLE`S286到目前为止，这个问题通常使用简单的图像处理技术来解决，其不准确性是众所周知的。本文提出的公式和算法可以作为评价现有融合方法质量的参考工具，特别是对设计新的融合方法有一定的帮助。不同3D图像的融合取决于一个更普遍的问题，即在固定精度的离散空间中数字坐标的变化，这种情况在我们现在的数字时代比比皆是给定放置在同一nD对象上的不同分辨率和大小的两个规则网格，数字坐标改变问题在于找到与两个参数化的整数坐标和附加到重叠体素的值相关的公式。如果有关粗坐标和精坐标的公式可以很容易地建立起来，那么只有当两个图像分辨率相当不同时，它们的精度才是可以接受的;当分辨率彼此接近时，它们就不足够了。在后一种情况下，必须确定相交多面体的实际体积并将其混合到插值公式中。让我们描述我们的方法和本文的内容数字坐标变换问题是多面体相交与格子周期性混合的问题。这个问题在一维中是微不足道的，在二维中丰富得多，并且随着维数的增加变得越来越复杂。这主要是由于相交问题的复杂性不断增加，同时也是由于没有明确的方法来利用格的周期性和立方体的简单性当n≥3时，结构似乎可以减少总工作量2D的情况是有指导意义的，因为基本的交叉问题很简单，并且晶格周期性的好处可以很容易地在连分数的帮助下解释。不幸的是，这种晶格性质并不像在更高的维度中那样普遍化。然而，这种2D的情况下给出了宝贵的指示，其中更高的维度的概念是值得介绍的，像平面截面的体素空间研究的第5，6部分，公式给出平面截面的立方体的第8部分和二次立方体的第13部分。这种立方体具有二次坐标顶点，这使得代数计算机系统能够精确地计算边和面的交点。这反过来又给出了形式上表示两个二次立方体相交的顶点的过程;这在第14部分中公开但是，到了这一点，在一般位置相交的两个立方体的看似简单的问题，令人困惑的问题出现了：应该使用为任何凸多面体设计的铝出租，还是应该尝试，相反，利用立方体特别是简单的结构和众多的对称性？对这个相交问题的直接考虑表明，沿着这两种方式中的任何一种的天真方法都会导致大量的计算，特别是不等式检查。看来，立方体的特殊结构，主要是它们的对称性（在第9、10、11部分中使用），导致需要检验的不等式的数量大大减少。本文给出了立方体平面截面的不含不等式的更精确的精确公式，REVEILLE`S287−ǁǁǁǁ(-b（a）（a、b）3D凸小平面相交算法（参见部分15），用于两个体素的相交的有效过程。一个很好的结果，正式的方法是删除处理的所有特殊情况下，从编程这些退化的情况下被包含在分析公式。2二维数字坐标变换问题当同一物体的不同方向和分辨率的两个图像的信息（例如灰度级）必须匹配时，就会出现这个问题这两种类型的像素可以被描述为小和大。我们假设Fig. 1.数字坐标改变。小像素由整数向量（m，0）和（0，m）定义（即，它们是单位正方形的m×m块），大像素由整数向量定义（a，b）和（−b，a）。很容易看到坐标为（x，y）的最大g像素大致对应于小像素（斧由Mbx+ayM）. 类似的公式给出作为小像素的坐标的函数的大像素的坐标可以写着当欧几里得范数（a，b）/m不太小（即大于10）时，这些简单的数字坐标变换公式的精度是足够的，但是当（a，b）/m很小时，这意味着两个网格具有相似的尺寸时，这些简单的数字坐标变换公式的精度是粗略的。在这种情况下，必须计算小像素和大像素的交叉区域，并将其考虑到插值公式中。但是通过蛮力计算像素的许多交叉点是必要的;试图优化这项任务是明智的，连续分数允许这样做。向量（a，b）定义大像素为整数，让我们假设a和b是互质的。设φ表示由下式引导的欧几里得直线：、REVEILLE`S288LLL −−{|联系我们{|联系我们（a、b）。所有的整数点k（a，b），其中k∈Z都属于n，同样对于由（−b，a）指向的欧几里得直线，也是如此（见图2）。①的人。事实上，顶点为（0，0），（a，b），（ab，b+a），（b，a）的正方形（这是大像素之一）平铺数字平面Z2，并且这种平铺是双周期的。因此，小像素和大像素的相交方案已经包含在其重叠的所有小像素的相交中。这减少了对覆盖边缘的小像素的考虑。显然，这些覆盖是4连接的数字线，其参数化可以容易地获得;这样的示例覆盖由fig说明。 2.让我们简单地说一下，如果向量（a，b）被任何实向量−→v所代替，会发生什么。首先，这并不限制应用，因为实向量−→vmay可以任意近似的整数向量，（见[7]）;众所周知，这样的近似可以获得与连续分数相关联，到−→v，（参见[8]）。第二，如果−→v的peα是无理数，则前一个有界平方图块变成无限的，并且相交方案通过由小像素覆盖整条线及其垂直线来因此，在数字坐标变换中出现的相交方案总是由覆盖欧几里得线段的4连通数字线给出;这些是有界的，如果−→v是整数，而在opposite情况下是无界的。图二.覆盖像素。如果α是无理数或大有理数，则该方案可能难以描述，因为在第一种情况下不存在覆盖像素的明显参数化，而在第二种情况下4-连通参数化可能涉及大整数。当α是无理数时，Klein这个图只是两个凸壳的一对，（x，y）x0的 和yαx类似地，（x，y）x0的和0yαx。 克莱因图中给出的两条折线y = αx通常被称为克莱因漏斗，见图2。3.第三章。众所周知（见[4]，[8]），Klein克莱因漏斗的优良命题2.1直线y = αx与单位像素的相交方案与克莱因漏斗的相交方案相同。这一结果表明，为了得到直线y=αx的相交方案与单位像素，它的表面构造的克莱因REVEILLE`S289§漏斗表面的两条多边形线）。由于克莱因这显示了连分数如何允许容易地确定同一对象的两个不同2D视图864202 4 6 8 10 12 14图3.第三章。克莱因任何小像素和大像素的交点面积公式都可以得到;第11节和第12节研究了一个密切相关的问题。3三维数字坐标变换问题3D图像的每个体素被同化为单位立方体，这样的图像自然是空间R3的平铺，其可以被称为体素空间或体素平铺。这个问题可以解决应用一个双循环的过程给出相交多面体的两个立方体，其中每一个运行到一个体素的集合但这种方法可能会错过存在于立方体的相交多面体之间的最终周期性，类似于为2D数字坐标变化问题所锯的周期性。虽然本文不会研究晶格周期性的所有方面（与多维连分式相关的方面只在16中概述），但考虑到的那些方面，称为体素空间平面截面，需要以某种方式偏离明显的方式。研究体素平铺的平面截面的几何形状及其与两个立方体的相交计算的关系将占据本文的其余部分。REVEILLE`S290∩∈F<$∈F ∈FFF4体素相交如果两个（并且只有两个）立方体的相交问题从理论的角度来看可能被认为是毫无价值的，那么如果它们属于两个体素平铺，则情况并非如此，因为在这种情况下，它们具有在它们的相交处再现的强规则性属性。如前所述，两个立方体的相交可以被认为是两个凸多面体相交的特殊情况，也可以被认为是两个具有特殊结构的物体的相交。如果第一种观点在计算几何中是众所周知的（[5]，[9]，[11]），第二种观点，据我们所知，从未发表过，即使计算机图像的一些技术，如3D裁剪，（[6]），似乎接近它。一般的计算几何算法，它使用相当复杂的数据结构来描述凸多面体，并不适合于目前的情况有两个主要原因。首先，它们避免了退化的情况：第一个立方体的顶点属于一个面或第二个立方体的边，共线边，与一个面共面的边，共面的面等等。第二个是，在精确算术或计算机代数系统中实现它们仍然是领先的（[11]），目前的实现与浮动点系统破坏了格的性质，这禁止他们的研究。这解释了我们的方法的主要特点，其目的是一个正式的描述两个立方体相交的多面体，即得到其顶点的精确代数值，因为它是一些晶格周期隐藏的地方但是，如果不使用立方体的特殊性，如它们的少量面和丰富的对称群，工作的计算机科学家知道理论算法和编码之间的漫长道路本文不仅讨论了体素求交算法的几何问题，还讨论了其编码问题。这段代码的效率，主要是删除了几乎所有发生在主循环中的控制测试，是我们这样做的第一个原因;第二个原因是保持在一个方便的数字系统中：二次代数系统。该算法的主要思想是确定两个立方体相交多面体的面。如果C和CJ是两个在一般位置的立方体，J表示它们的面的集合，则多面体K=CCJ的任何面要么由面f支撑，要么由一个fjj支撑。如果K的面由面f支撑，则它是两个多边形的交点，即supp（f）CJ和正方形f，反之亦然。f的支撑平面与立方CJ的交点）。事实上，两个族中的一个族的体素的支撑平面切割另一个族的无限数量的立方体。但是这些支撑平面是周期性间隔的，因此体素平铺的周期性应该出现在其平面截面中，该平面截面是2D普通平铺。合理的剖面假设大大简化了这些二维诱导瓦片，我们能够给出明确的公式，其顶点，反过来，大大减少了REVEILLE`S291我们的相交算法的复杂性，这最后一个被减少到多边形相交计算，一个相当简单的任务。见图4。体素求交算法原理。观察到两个立方体C和CJ的相交多面体K具有至多12个面，这些面是其中一个立方体的正方形面与另一个立方体的由它们的支撑平面引起的平面截面的公共部分，相交算法为：在两个立方体的所有12个面上循环由该面的支撑平面确定另一个立方体的平面截面计算两个多边形的相交：平面截面和正方形面。图4示出了左立方体的一个面以及其支撑平面与右侧支撑平面的相交。正方形面和相交多边形的公共部分是多面体K=C<$CJ的面。5体素空间体素拼接的平面截面的几何学研究在医学成像中也称为信息提取，它是我们求交算法中的一个关键部分从图5可以看出，问题是描述所有的交叉点一个给定的平面与平铺的所有体素填充空间R3。 即我们想要描述该平面的所有相交多边形与它遇到的每个体素。但是，众所周知，几何对象（如多边形和多面体）的近似知识破坏了它们的大部分理论性质，避免了对退化情况的完整研究和实现。 这就是为什么我们不仅要寻找这些多边形顶点的数值近似，而且要寻找它们坐标的解析表达式REVEILLE`S292图五、体素平铺的平面部分图5示出了由平面（P）切割的一些体素，它们在平面（P）上引起普通的2D平铺（平面的这种平铺不能与由体素引起的R3如下图6所示，这种二维平铺是通过将平面与三个等距平面族相交而构建的，这些等距平面的方程为x=k，y=l，z=m，并且其中k，l，m是整数（这里我们认为小的体素是单位立方体）。图六、由截面上的体素引起的2D平铺让我们假设我们的相交平面（P）通过0并且有一个整数REVEILLE`S293≤ ≤≤||||法向量（a，b，c），其中0 a bc，且gcd（a，b，c）= 1。因此方程（P）为ax+by+cz= 0。约束条件规定（a，b，c）为整数似乎令人惊讶，但我们记得，R3中的任何实方向u都具有任意接近的丢番图近似，而且，很容易找到整数a，b，c，使得角度（u，（a，b，c））任意小。我们还将看到，如果使用八面体群Oh，约束0≤a≤b≤c可以放松，而不会产生大量的控制检验。 最后一点将在下面治疗。在平面（P）上的基上进行的一个简单计算表明，上述平铺是由下列四个参数族引起的行，x像往常一样是可变的。d1（x，a，b，c，k）=k（a2+b2+c2）/sqrt（a2+b2）√d2（x，a，b，c）=（ax−kd3（x，a，b，c）=−（bx+k（a2+b2））sqrt（a2+b2+c2）/bc√（a2+b2））sqrt（a2+b2+c2）/ac图6示出了由方程为3x+7y+13z= 0的平面（P）上的空间体素引起的平铺，图6是使用这三个线族绘制的。6平面平铺的算术化表示周期性。仔细观察图6，可以发现周期的存在，即相同的瓦片。 这是我们关于（P）法向量（a，b，c）的算术假设的（好）结果之一。问题是要表达这些时期，是两个平移向量和一个周期的所有瓦片;这就是算术的用武之地每个体素都连接到其左下底顶点（x，y，z），我们用f（x，y，z）表示线性形式ax+by+cz。显然，整数点M=（x，y，z）和MJ=（xJ，yJ，zJ）到平面（P）的距离相等，当且仅当f（x，y，z）= f（xJ，yJ，zJ），并且M和MJ在由（P）界定的相同半空间中，并且当且仅当f（x，y，z）= f（xJ，yJ，zJ）时，M和M j到（P）的距离相同。我们由此推出分别附着在M和MJ上的立方体C和CJ与平面（P）有相同的交多边形当且仅当f（x，y，z）= f（xJ，yJ，zJ）.因此，平面（P）平铺中的周期性等价于解决丢番图问题(1)ax+by+cz= 0也就是找到这个方程的整数解（x，y，z）。虽然有几篇论文给出了完整的证明和算法来解决一般的线性丢番图方程组，（例如见[7]），我们将简要描述一个算法，由于Blanklin（见[3]），它解决了限制REVEILLE`S294等式（1）的情况。理解和使用这个算法应该是清楚的，一旦我们提醒，解决方案（1）使一个二维格包含在（P）一个格有许多基，但2个向量的自由族不一定是基。给定丢番图方程ax+bx+cz= 0，我们建立矩阵a B C1 0 00 1 00 0 1并应用由第一行条目a、b、c驱动的列操作，以便取消它们之间的两个值。位于所得矩阵的第1行的两个零点之下的第2、 3、4行中的值给出了Z3的两个向量，这两个向量构成了（1）的解的基础。让我们解决丢番图方程3x+7y+13z= 0来说明这个过程，这是一个简单的整数克隆高斯算法。以下矩阵显示了必须完成的连续步骤。1−41 000 1可以验证。整数参数方程（1）求解方法。如果M=（x，y，z）s平面（P），则所有多边形（在奥利瓦。它仍然需要描述D.所有这些瓷砖都是我们用 F1表示的时刻，3713311001 0 01 −2 −47 20 1 00 −3 10 0 1 0 0 1 −1因此，一个可能的基是（7，-3， 0），（2， 1，-1）作为rea。使用相同的思想，他（她）将能够推导丢番图方程ax+by+cz=k，其中k是另一个（提示：开始求解ax+by+cz= 1）。已知与下面的平面晶格的eq周期性问题相关联的晶格的基{v1，v2}是整数点，使得在M切割处附着的体素整数点，其中具有sa平移的附加体素切割（P）是点(2)M+k1v1+k2v2其中k1和k2是任意整数。7相交p的显式描述一旦获得（P）的平铺的周期向量，精确地，每个平铺出现在给定单位立方体的一个周期平面部分内。设立方（C）为[0，1]3，[0，1]为R的单位区间.REVEILLE`S295∈≤R3的方向子集满足0≤a≤b≤c（三）（四）F1类似地，F2表示子集F2a+b≤c0≤a≤b≤ca+b≥c定理7.1当k从−∞到∞变化时，单位立方体与平面ax + by + cz = k相交的顶点序列与（a，b，c）无关，只要它在F1中。当然，当（a，b，c）F2时，这个结果也成立，但交叉点的顺序改变了.在我们给出顶点交叉的精确顺序之前，先对这个定理说几句话，以解释它的目的。无论切割立方体的平面是什么，所得到的相交多边形的顶点都取自该平面与立方体边的12条支撑线的12个交点。但是立方体的平面（非退化）截面的边数大于3且小于6，这意味着必须消除一些先前的相交点以获得正确的多边形。对于每个多边形，这个定理首先将帮助确定它在12个可能的点中的顶点，（注意，选择随多边形类型而变化），其次，但这是程序员例如，如果在所有多边形中只考虑三角形，很明显只有8种情况（接近立方体顶点）;但是如果没有前一个定理，其他多边形类型的情况数量将更难找到，更糟糕的是，给出顶点坐标的公式的整个编程将非常无聊和难以调试。如果平面（P）平行于自身移动，在约束a+b c下，可以立即看到（C）顶点交叉的顺序由下式给出：八个值f（v1），f（v2），.，f（v8）在（C）个顶点上由形式f所取。立方（C）=[0，1]3的定义给出了这8个值的评估，它们是0，a，b，c，a+b，a+c，b+c，a+b+c，并且假设（a，b，c）∈F1立即表明它们的递增排序是：(5)0< a b a +b c a+c b+c a+b+c由于f与到平面（P）的距离成比例，因此这组不等式表示（C）顶点与（P）位移交叉的顺序是：（0，0，0），（1，0，0），（0，1，0），（1，1，0），（0，0，1），（1，0，1），（0，1，1），（1，1，1）。由于（C）在（0，0，0）处相交的边是[（0，0，0），（1，0，0）]，[（0，0，0），（0，1，0）]和[（0，0，0），（0，0，1）]，我们推出，如果（P）使得f值满足REVEILLE`S296∩≤∩≤0f（x，y，z）=0）则eps：= 1：u：=k其他eps：= 1：u：=kend if;如果是（c a+b），则如果是（0=u）且是（u l （a+b c）），则eps hex（u）elif是（l（a+bc）=u）并且是（u l（a b+c）），则eps笔（u）elif是（l （a b+c）=u）并且是（ul（c a+b）），则[医]前列腺素其他eps三（u）结束，如果否则如果是（0=u）且是（u l（c ab）），则每股面值elif是（l （c a b）=u）并且是（u l（a b+c）），则eps笔（u）elif是（l（a b+c）=u）并且是（ul（c a+b）），则[医]前列腺素其他eps三（u）结束if;结束if;endif;endproc;本文给出了当参数a，b，c满足F1或F2时，即（a，b，c）F，集合F是（C）对称群的基本域时，当然，对于其他情况也存在类似的公式，但使用控制测试来找到它们将非常困难。尴尬：代码编写和调试将是漫长而乏味的任务。使用八面体群，这是一REVEILLE`S301个立方体的对称群，将简化这一推广。REVEILLE`S302S≤ ≤≤∈9正方体的八面体对称群这个群，记为Oh，可以用8阶群（Z/2Z）3与6阶三字母排列群S3的乘积来标识因此，Oh的阶数等于48。虽然Oh通常借助于立方体的对称性（如（C））来解释，但我们将脱离这一过程，采用一个等价的方法，其中Oh将任何空间方向约化为属于其基本域F的规范方向，该域由不等式0a b c定义。接下来，O h的第一个因子，即群（Z/2 Z）3，将对应于Z 3向量的三个坐标a，b，c的明显的8符号排列，而第二个因子3将映射到这三个值的6个排列。图8.第八条。八面体群的几何视图所以空间的任何方向都等价于F通过一个元素的方向 的H。这组在48个三角形中归纳出（C）边界的三角剖分。 三角形是基本域F与（C）边界的交点：{（x，y，z）|最大值（|X|、|y|、|z|）= 1}也将由F表示。换句话说，对于Z3的任何向量v，这样g.v属于F. 当然，g−1.g.v=v，也就是说v可以很容易地恢复。10执行组Oh我们详细介绍了这一点，这可能是令人费解的读者不熟悉正式的计算系统，也因为存在几个解决方案，其中一些导致效率低下的编码。我们保留的解决方案使用组Oh的线性表示，其由3× 3矩阵的集合给出，其中每行和每列都有一个，REVEILLE`S303−−−×××→≤只有一个非零项等于1或1。这是一个简单的证明，证明O（3）（3维正交群）的这个子群同构于Oh。 重要的一点是找到元素g，它将减少给定的向量v，如上所述。给出与v=（a，b，c）相关联的矩阵g的算法工作如下：• 建立43表边界的3第一行等于的3个单位矩阵到（a，b，c），（与前一个相同的结构开始Blanklenship al-出租m）.• 将列乘以其第一个元素的符号（如果a=6，则将列1乘以①的人。 因此，第一行的所有元素都变为正数（或零）下面的1号将保留这个珍贵的标志注意给空元素指定符号1（而不是像某些编译器那样指定符号0）。• 然后交换列，以便第一行按递增顺序排序。结果第一行以下的3 × 3矩阵h是所寻找的g元素的逆，即转置（h），v = w属于F，h，w = v.由于O h包含在O（3）中，我们有g=h−1=转置（h）。保持矩阵h或转置（h）=g在数学上是等价的，但从编码的角度来看，w作为先前算法的第一行给出，不需要g;相反，h可以将w（以及当（a，b，c）∈F时找到的所有数据）转换回v的域。该算法的伪代码，给出与向量v=[v[1]，v[2]，v[3]]可以写成如下。LinearRepr：=proc（v）localM;M：=[[v[1]，sgn（v[1]），0， 0]，[v[2]， 0，sgn（v[2]），0]，[v[3]， 0， 0，sgn（v[3]）]];M：=sort（M，（x，y）is（abs（x[1]）abs（y[1]）;[abs（M[1][1]），abs（M[2][1]），abs（M[3][1]）]，array（[[M[1][2]，M[2][2]，M[3][2]]，[M[1][3]，M[2][3]，M[3][3]]，[M[1][4]，M[2][4]，M[3][4]）;结束进程：11体素空间截面与数字平面的关系。有几种方法可以从3D图像中提取平面截面。第一个是上面指出的，并开始切割给定平面的所有体素，导致已经研究的平面平铺。该方法的第二步包括将这些图块映射到与该图块重叠的像素这种对应关系的精确处理假设确定瓦片和像素的所有交点，这可以使用两个凸多边形的交点算法来完成REVEILLE`S304§∈××- 不超过-≤图9.第九条。立方体的平面截面是一个周期的瓦片（见15）。我们不会深入讨论这种第一种提取方法的所有细节，而是使用数字平面提出一种更简单的方法必须指定体素和整数点之间的对应关系。连接到点（x，y，z）的体素，记为V（x，y，z），是区间[x，x+1][y，y+1][z，z+1]的卡累利阿乘积，因此（x，y，z）是它的左下底顶点。因此，体素V（x，y，z）切割平面（P），其方程为ax+by+cz= 0，（条件（a，b，c）F仍然满足）当且仅当：ax+by+cz≤0且0≤a（x+1）+b（y+1）+c（z+1）因此，相交（P）的体素属于由下式定义的厚数字平面：−（a + b + c）≤ ax + by + cz ≤ 0。这组体素是6-连接的，并且比由c ax + by + cz 0定义的朴素18-连接数字平面更厚（其名称）（注释c = max（a，b，c））。这表明最多存在a+b+ c个不同的瓦片，并且使用朴素数字平面不足以提取3D图像中的切片;使用6-连接的厚数字平面获得更好的结果。然而，使用以c/2ax +by+cz c/2定义的欧几里得平面为中心的明智的18连接数字平面可以得到很好的近似，因为这些平面将给出最大的瓦片。