现代医学统计学(第2版)
上QQ阅读APP看书,第一时间看更新

第六章 统计方法在医学图像分析中的应用

金声 1 赵耐青 2
1澳大利亚悉尼大学计算机系
2复旦大学公共卫生学院卫生统计学教研室
医学图像是二维随机信号。随机信号研究中有许多共同的问题:噪音的过滤、信号的还原以及信号的采样等。此外还有许多仅与高维信号有关的特殊问题,比如图像分割、像素聚类等。这一章将讨论医学图像的问题,特别将讨论统计方法在这个领域中的应用。
第一节 概述
医学影像学由于能得到丰富的信息及采用更多的方法而快速发展,这些方法包括:核磁共振(MRI)、X线传输影像(X-ray)、计算机断层(CT)、超声(2维和3维)、正电子发射断层(PET)、单光子计算机断层(SPECT)、磁源成像(MSI)、电子源成像(ESI)、X线钼靶(MG)、曲面体层全颌摄影(OPG)等。
MRI在临床诊断和生物医学研究中是一种非常有效的非侵袭性技术,它采用了核磁共振(NMR)这一非常著名的化学、物理、分子生物学分析方式。MRI的主要技术是获得解剖影像,并且还能提供组织、流体弥散、运动的物理-化学状态信息。磁共振分光术(MRS)可提供化学/组成信息。MRI引发脑组织、脊髓、肌肉骨骼系统影像检查的革命性进展。超强的软组织对照和间隙分辨使MRI用于许多中枢和整形疾病的检查。
X线是通过加速电子与一种靶物质(通常为钨)相撞击而产生。X线的偏转和吸收程度在人体不同的组织和骨骼中是不同的,吸收量依赖于组织构成。例如,骨密度物质吸收的X线量比软组织(如肌肉,脂肪,血液)明显得多。偏转的量依赖于组织中电子密度。因此,由于照射骨或金属的X线到达X线相片的光子量比软组织少,前者的X线相片看起来较亮。
CT自20世纪70年代起开始广泛运用,并被认为是医学科学的一次重大进步。X线、CT检查可获得气体、软组织和骨骼位置的解剖学信息。三维影像是通过围绕病人旋转的X线发射器和对不同角度传送强度计算来获得的。
当前医学所用的超声是一种即时断层影像方式。它不仅能获得反射界面(器官和结构的内部)的即时断层影像,而且能获得组织和血液运动的即时断层影像。
PET的历史可追溯到20世纪50年代早期,当时波士顿的科学家第一次认识到放射活性同位素的特殊分级可产生医学影像。在大多数的放射活性同位素通过释放γ-射线和电子衰竭的同时,一小部分通过释放质子衰竭。质子可看作正电子。由于X线CT辅助重建运算法的发展和核子探测技术的进步,PET技术引起了广泛的兴趣并加速了发展。到20世纪80年代中叶,PET成为一种医学诊断的工具,用于人体代谢的动态研究和脑部活动的研究。在脑部和其他人体组织的局部代谢和中枢受体活性研究中,PET的敏感性是其他技术的上百万倍。相反,核磁共振对于解剖研究和流体或血管研究更准确。另外,虽然核磁共振光谱对于组织化学组成的评价有独到之处,但是纳米水平的测量不如毫米水平的测量,而人体中大多数受体蛋白却集中于纳米水平。正电子发射断层是理想的影像学技术。PET的主要临床应用在于脑部、乳房、心脏、肺、结肠和直肠肿瘤的发现。另一种运用是通过心肌代谢成像来评估冠状动脉疾病。
SPECT和PET一样,通过注入人体的放射活性同位素的浓集来获得图像。SPECT的资料开始于20世纪60年代早期,D. E. Kuhl和R. Q. Edwards首先介绍了发散横断层理念,早于PET,X线CT或MRI。
心脏和脑部的神经元中铁元素的电流产生体外磁场。这些磁场可通过SQUID(超导量子干扰装置)列来测量,探测器置于心脏或胸部之上或旁边。头部的磁场记录被称为磁脑成像(MEG),同样心脏磁场记录被称为磁心成像(MCG)。磁源成像(MSI)是通过体外磁场重建心脏或脑部电流的总称。
电子源成像(ESI)是通过测量头皮或躯干电势进行脑部或心脏的电子活性重建成像的一种影像方式。标准的脑电图(EEG),心电图(ECG)和心脏向量图(VCG)技术受到它们提供脑部和心脏局部电子活动或局限生物电子事件信息能力的限制。非侵袭性的ESI要求同时用20个或更多的电极记录脑部电势,用100至250个躯干电极来勾画心脏的体表电势。
X线钼靶(MG)是一种有效的诊断乳房肿瘤的方法。低剂量的X线乳房摄片可用于发现无症状的早期乳腺癌。它还可以用于获取诊断性的乳房摄片。乳房穿刺定位可为外科手术提供定位信息和细致的组织信息。
曲面体层全颌摄影(OPG)和测颅X线片是最新用于牙齿或牙齿矫正评价的技术。
本章第二节讨论压缩,第三节滤波,第四节分割和第五节配准,最后是结论。
第二节 应用图像的统计特性进行采样和压缩
以计算机为依托的高级医学图像技术(如PET)已经在现代医学研究和诊断中扮演重要角色并起到重要作用。然而伴随这些重要技术而来的问题是图像数据量的增长。例如,应用CTI951扫描仪作20至30时间点的常规动态PET研究一般需要31幅128×128像素的截面图像。这个图像将包含11×10 6个4位数的数据,近似占用22兆存储器空间。随着目前PET图形分辨率的改善,图像数据所需的存储器容量将进一步增加。因此发展有效的图像压缩技术是非常有意义的,它有助于当前医学数字化、图像库管理和远程医学的迅速发展。
运用特定环境下动态PET图像和生理追踪子运动建模的生理运动知识,本节将介绍对动态PET图像数据进行压缩而基本不丢失数据的一种新奇算法。这个压缩算法由3个步骤组成:①利用优化图像采样点,在时域中压缩数据;②通过聚类分析,在空间域中压缩数据;③利用标准图像压缩技术建立图像压缩的索引。该算法被用于压缩二氟脱氧葡萄([ 18 F]2-fluoro-deoxy-glucose,FDG)追踪子的临床人类大脑PET图像。这项技术同样可以用到不同追踪的其他PET图像中。这项技术 [1]已在临床动态PET图像上实现。实验结果将用来展示压缩技术的性能和重建图像的质量。
一、追踪子运动建模和功能性图像
追踪子动态技术被广泛地应用到PET图像中提取描述人体动态过程的有用信息。这些信息通常用一个数学模型 ut| p)描述,其中 t=1,2,…, T以及 p是这个模型的参数。这些参数描述追踪的起始位置、传送和生物化学转变。模型的驱动函数是血浆输入函数,通常由血浆采样值得到 [1]。从PET所定义的组织活动时曲线(time activity curve,TAC)或输出函数得到测量值,并记为 z it)( t=1,2,…, Ti=1,2,…, I),其中 t表示离散采样时间, i表示对应在这个图像区域中的第 i个像素。动态PET图像分析的目的是为在这个图像区域的每个像素获得追踪子的TACs和参数估计。然后应用这些参数来诠释生理参数。比如,大脑局部的葡萄糖新陈代谢的比率(local cerebral metabolic rate of glucose,LCMRGlc)。
传统的建模方法需要完整的PET投影数据。以像素与像素之间的联系为基础,应用某个快速估计算法 [1,2,3],通过参数估计产生功能图像。在本节,利用Patlak方法 [2,4]来产生LCMRGlc功能图像,以比较原始数据和压缩数据的准确度。
二、在时域和空间域的采样和压缩
根据追踪动力学的统计特性,采样和压缩方案由3步骤组成 [5]
步骤1:使用最优化的图像采样方案在时域进行数据压缩
在动态PET研究中,时域帧的可靠性直接受采样方案和获得数据所需持续时间的影响。获取数据的持续时间越长和帧扫描速度越快,时域帧图像就越可靠。然而,为了从动态进程中获得大量的信息,需要一定数量的时域帧。最近有报道,最小时域帧数的要求为模型中要估计的参数个数 (式6)。基于上述原因,一种为获得PET数据的最大信息量而自动产生最优化采样方案(optimal image sampling schedule,OISS)的算法应运而生 (式6)。这种算法利用累计PET测量值或对PET测量值积分。
在OISS的设计中,我们提出一个新的基于Fisher信息矩阵的目标函数以限制动态信息的损失 [7]。这个目标函数用于区别各个实验协议与采样方案之间的差别。由于OISS能直接获取PET投影数据,时域帧数可以减少,以至相应减少数据存储空间。进而言之,时域帧数的减少可进一步减少图像重建的计算量。详细算法可查阅文献 (式6)
步骤2:通过聚类分析,在时域压缩数据
先验知识在示踪动态模型中以PET示踪子累积测量相对时间序列的形式表现。利用聚类分析,从模型中能提取全图的TAC(image-wide TACs),并根据它们的运动相似性进一步分成不同的TAC组,对应于不同的组织区域。
聚类分析的目的在于根据相似的特征度量,把全图的TAC曲线 z it)( i=1,2,…, I)分到 C j个聚类组( j=1,2,…, JJI)。高度不相似的TAC曲线将属于不同的组 [8]。应该注意的是每个TAC必须唯一属于一个聚类组。在本文中,临床动态PET图像数据分类采用的是基于欧几里得距离度量为指标的系统-聚合的聚类算法。
利用聚类分析的结果产生一个索引表,它包含每个聚类的平均TAC和一个索引图像。索引图像代表从一个聚类映射到一个像素TAC位置的映射。这个索引图像与索引表共同形成了压缩时域和空间数据的基础。可区分的PET聚类组的总数通常不会超过64。因此一个8位索引图像足够表达这些聚类的映射。
步骤3:索引图像的压缩
本文中,用一种无损失压缩的方法来进一步压缩索引图像。用可移植的网络图形(portable network graphics,PNG) [9]格式来压缩并且存储从聚类分析所获得的索引图像。目前定义和实施的PNG编码技术是以一个32KB活动窗口中的缩小/放大为基础的。选用PNG格式是由于它的可移植性、灵活性和合法性而优于其他无损失图像压缩文件格式。PNG数据格式的细节可在文献 [9]中找到。
人的动态FDG-PET大脑研究是利用一个8环、15片扫描仪(GE/Scanditronix PC4096-15WB)进行的。这种扫描仪在可见区的中心包含4096检测器,输出轴和转换轴的最大图形分辨率的一半(FWHM)为6.5mm宽。在静脉注射200mBq与400mBq之间(近似0.5mg)的FDG之后立即开始动脉的血采样。在开始的第一个2分钟时间内,以0.25分钟的时间间隔进行8次血采样(每次2~3ml),然后分别在2.5、3、3.5、7、10、15、20、30、60、90、120分钟进行。为了测定血浆FDG和“冷的”葡萄糖浓度,这些采样样品是立即放置在冰上随后进行血浆分离。图6-1(a)显示的图像是从一个患者研究的第15个平面的一套时域帧图像。由于在开始的几幅图像所对应的示踪浓度较低,这些图像均经过加工以便看得清楚些。图6-1(b)显示的图像是在时域内,对图6-1(a)用上述压缩方法得到的5幅时域帧图像(等级度量的)。
图6-1 压缩结果
第三节 用各向异性统计扩散减少噪声
在量子物理、材料科学、流体动力学、核科学、医学和物理化学领域中广泛使用扩散过程。Perona和Malik [10,11]把该方法引入图像处理领域并提出多尺度平滑和边缘检测方案。它能在消除噪声成分的同时保存有效高频成分,也就是边缘检测 [12]
扩散是一个迭代过程。扩散的程度取决于扩散的阈值,即对比度的截断值。当对比度高于阈值时,在扩散过程中对比度将加强;对比度低于阈值,扩散中对比度将平滑消除。阈值的选定对于过滤过程是重要的。然而,阈值因图像而变化。这个问题与不同区域有不同的对比度以及同一图像同一区域的亮度失真混在一起。因此需要有一个选择阈值的适宜标准。
扩散进程中的阈值与这个图像边缘的对比度紧密相关。通过分析局部对比度选择阈值。在低对比度的图像中,特别含有噪声并且信噪比(signal-noise ratio,SNR)低的情况下,区域之间的对比度差别不明显,要确定阈值是很困难的。当信号中含有噪声时,主要困难为随机信号的分布未知和多个未知干扰的组合。在大多数情形,区域的直方图显示一个单峰。许多自动的阈值选择机制要求双峰(Bi-peak)的直方图如Tsai、Chen [13]和Bhandari等 [14]。在许多情形下,双峰直方图或多峰(multi-peak)的直方图不一定存在。Luijendijk [15]建议利用4连结区域的计数建立两个直方图,进行自动阈值选择。Tseng和Huang [16]建议利用边缘信息(即沿边缘区间的亮度)选择阈值。Nagawa和Rosenfeld [17]用直方图数据拟合2个高斯函数,Cho等 [18]应用偏移校正因子。Glaseby [19]利用迭代改善了它们的结果。然而正态分布的假设是很牵强的而且校正也不能弥补这个重要缺陷。另外,迭代使计算很费时。
确定阈值的另一个困难是图像亮度失真。直方图分析假定在图像中感兴趣的目标或区域的所有像素具有相似的灰度。然而,对大多数图像这假定是不成立的。Rodriguez和Mitchell [20]使用2阶段提取背景的一种自适应阈值方法。第一步使用全局的阈值提取各区域的结构,第二步局部优化。Parker [21]在一个目标里找到一个种子像素后利用局部阈值产生一个区域。Spann和Horne [22]在四叉树结构中从低分辨率到高分辨率产生几个区域。自适应方案是与亮度失真抗衡的一个合适方法。然而,上述试凑方法没有一个坚强的理论基础。
本节应用中心极限定理描述了一个自适应扩散方案。利用回归分析从含许多对象的单峰直方图中,分离出一个局部窗口中主要对象的分布。这种分离有助于自动确定阈值。我们应用这个算法从X线血管造影图(X-ray angiogram,XRA)提取大脑动脉的图像。这个算法对于单峰且没有谷底的分布直方图很有效。肾的显微镜图像中有多个可视物体(visual objects)且各个可视物体之间的对比度差别很小。该算法也被用于这些显微镜图像的滤波。这个方案显示完全自动过滤的过程是能够实现的。对于有纹理图案并且噪声分布未知的图像,该方法很有效,而传统滤波方案(如小波滤噪声)则显得力所不及 [23]
一、非线性的各向异性扩散
低通滤波常被用来滤除噪声。大多数滤波器是各向同性的(各方向函数变化对应相同)。小心地审视这类问题,我们注意到沿着一边的梯度不是各向同性的。垂直于边缘的梯度值是最大的并且沿着边缘扩散。因此在平行于边缘的方向时,适当地增加平滑功能并且在垂直于边缘的方向停止平滑。非线性的各向异性扩散的函数形式为:
(6-1)
其中 Ixyt)表示信号, gI)是扩展函数的梯度。
经常用到的2个扩展函数是
(6-2)
(6-3)
扩散滤波器的计算可以由差分运算进行
扩展扩散使得区域内平滑优先于越过边界的平滑。这方法的基础是通过选择合适的局部扩展强度抑制跨边界的平滑。函数中的 κ值起到很大作用。如果 κ值设得太高,则滤波器将起平滑滤波器的作用在横越边界处扩展。当 κ值设得太低,则小幅度推进导致许多迭代。在某些 κ值下,滤波产生额外边界。因此,在设计中的一个重要问题是选定 κ
二、对比度截断值(cut-off)的选定
需要处理的图像经常在许多亮度层的对比度非常低。对于这样的图像,决定适当的阈值是困难的。图6-2显示大脑动脉的一幅XRA图像(a)以及对应的单峰直方图(b)。从直方图上选定阈值是不明确的。试凑法是不可行的。我们提出了根据区域利用回归分析进行动态选择阈值的方法。
图6-2 在背景图像上的直方图分析
1.由回归分析和似然比分类选择阈值
我们方案的理论基础是中心极限定理。因为低对比度以及背景比例过高,所以从该图像背景中分割出大脑动脉很困难。虽然图像背景像素直方图分布是未知的,然而根据中心极限定理可知,若x 1,x 2,…,x n是独立同分布的随机变量且数学期望为μ以及方差为有限值σ 2,则当 n充分大时, y= 渐近服从N( μσ 2/ n[24]。利用正态分布的回归模型可以使前景直方图和背景直方图分开,图6-3中,阴影区显示背景直方图,深黑的区域显示前景直方图。在分开直方图以后,很容易选择阈值进行图像分割以及分析前景对象。
图6-3 用正态分布的回归分析分割直方图
从部分的直方图获得回归分析所需的样本数据。我们计算直方图的平均值,采用方差较小的一半,然后我们找到那一半直方图的众数。采样数据 h iiS,是上述直方图中与此众数同侧的数据(相对于均数而言),回归分析获得
然而,当背景像素的个数不是足够大时,在回归中应用正态分布是不合适的。图6-4显示另一幅XRA图像(a)和它的直方图(b)。图6-4(c)是在整个背景区域上用正态分布回归后所得的直方图。在2个峰之间没有出现我们所期望的谷底,意味着没有清楚地分离这些均数。对于这种情况,我们改用Rayleigh分布作回归分析。概率论表明,当 n不是足够大时, 满足Rayleigh分布。
图6-4 用回归分析进行分割
应用Rayleigh分布的回归分析得到
2.提取对比度的截断值
扩散过程紧紧依赖于(6-2)和(6-3)函数中的参数 κ的临界值。这参数 κ与对比度有关。下列关于提取 κ的讨论是以(6-2)的扩展函数为基础的,但通过 κ除以 很容易转化为(6-3)。虽然我们能对图像上的一个可视对象的分布获得一个合适估计,例如,XRA图像中的背景,但是我们没有其他对象的信息,例如,XRA图像中的血管。估计2个对象之间的平均对比度是很困难的。在我们从前景直方图中分离背景直方图后,我们应用似然分类法 [25]来分离两个对象的像素。把这两个直方图作为似然分类中两个类别的概率分布。我们从下列等式计算 κ值。
其中, N l是灰色级别为 l的像素个数, P是以邻近像素对为元素的集合,这些像素对中的两个元素分别属于两个不同聚类。这计算可以限制在一个局部区域内。如果两个邻近像素属于不同聚类,则它们的差值将集中反映在一个差值的直方图中。可以从在一个局部区域内差值的众数提取对比度。
第四节 医学图像的分割
分割是把图像分成各个成分或各个部分的过程。它经常是图像分析的首要步骤。有效的分割通常能使图像分析最终成功。正是这个原因,世界各国的研究人员开发了许多分割技术 [26]。亮度图像的分割通常涉及4个主要途径:阈值法、边界检测法、区域法和混合法。
阈值技术 [27]基于一个假定,即在某个范围以内的所有的像素属于同一等级。这种方法忽视了所有的图像空间信息,因而不能很好处理含有噪声或边界模糊的图像。
基于边界的方法有时称为边缘探测法 [28],因为它们假定像素值在两个区域之间的边界线出现快速变化。基本的方法是对图像进行梯度滤波。滤波器的许多高值提供了候选区域边界线,这些候选区域边界线必须修改后产生闭曲线,以代表这些区域的边界线。
区域分割算法假定在同一个区域内的邻近像素有相似的亮度值。最著名的分离-合并技术 [29]就是根据齐性准则进行的。它包括种子区域增长 [30]和非种子区域增长。
混合法综合了上述一种准则或多种的准则。这方法包括形态上的分水岭 [31]分割,变量顺序值曲面拟合 [32]和动态轮廓 [33]方法。
本节所介绍是应用统计特性进行图像分割的两个方法。
一、应用期望-最大EM算法进行概率分割
现已证实基于亮度作MR图像分类是有问题的,即使采用现代技术也是如此。扫描内和扫描间亮度不均匀是困难的一个共同根源。虽然有报道称一些方法成功地校正了扫描内不均匀的问题,但这些方法需要对单个扫描进行指导。这节将描述一个被称为自适应分割的新方法,该方法利用组织亮度特性和亮度不均匀的知识校正和分割MR图像。使用EM算法导出的这种方法能够较正确分割组织类型,也能够较好地显示MRI造影数据。一个包括1000个大脑扫描图像的研究已证实该方法是有效的。我们将描述在下列类型的图像中分割大脑的实施和结果:横断面(双回波自旋回波)和冠状面(3DFT梯度回波T1-加权)的图像全部用传统的头部线圈得到;矢状面的图像用表面线圈得到。当分割灰质和白质时,自适应分割的准确性能与人工分割相比并比监控多变量分割更接近人工分割。
MRI图像在形态上的现代应用经常要求把体积图像分割成各种组织类型。人们常将统计分类方法应用于亮度信号,并结合形态图像处理来实现这种组织分割 [34,35,36,37]
现已证实,当即使采用先进技术如非参数,多信道方法处理MR图像时,传统的以MR图像的亮度为基础的分类是有问题的。由于RF线圈或获得序列(例如,在梯度回波图像中人工磁化系数)的原因,层内扫描亮度不均匀性是困难的公共来源。虽然MRI图像可以在视觉上显得均匀,层内扫描的不齐性经常扰乱以亮度为基础的分割法。在理想的情形中,大脑中的白质和灰质很容易区分,因为这些组织类型显示出互不相同的信号亮度。在实际中,空间亮度不齐性的幅度经常足以引起与这些组织分类有关联的信号亮度分布显著地重叠。另外,MR设备操作的条件和状态影响观测的亮度也可引起内部扫描亮度显著的不均匀性。而这一点经常需要在每个扫描基础上进行人工调整和改善。
层内和层间扫描MRI亮度不均匀性用一个空间变化因子建模,并称该因子为增益变量(gain field)。它放大亮度数据。对亮度作对数转变使它成为一个人为可加的有偏变量。如果增益变量是已知的,则通过应用以亮度为基础的传统分割方法校正数据从而相对较容易地估计组织分类。同样,如果组织分类是已知的,则通过直接比较亮度预测值与亮度观察值也很容易估计增益变量。然而,如果没有增益信息就不容易估计组织类型,同样没有组织类型信息就不容易估计增益变量。用迭代算法进行估计增益变量和组织类型(典型的情况下,迭代5次至10次就收敛了)的可能情况。
用Bayesian方法近似估计偏差变量。因为偏差变量描述了在MR亮度数据的对数转化后的人工增益,因此先对亮度数据取对数。
Y i= gX i)=(ln([ X i1),ln([ X i2),…,ln([ X im)) T
(6-4)
其中 X i是在MRI观察图像中第 i体素(体积的像素,voxel)信号亮度, m是MRI图像信号的维数。其他的统计逼近类似于MRI图像的亮度分割方法 [36,37],观察值的分布在建模时作为正态分布(包含一个显式偏差信息变量):
(6-5)
其中
m维正态分布方差为 Y i是在第 i个体素的亮度观察值的对数值。
Γ i是第 i个体素所属的组织分类;
μx)是组织分类为 x上的亮度均数;
ψ x是组织分类 x的协方差矩阵;
β i是在第 i个体素的偏倚变量。
Y iμx)和 β i表示 m维列向量, ψ x表示 m× m矩阵。每个体素的每个亮度的对数值都分别对应偏倚信息变量值。通俗地说,在已知给定组织分类和由正态分布的中心位置在偏移平均亮度所给出的这个分类的偏倚信息变量的情况下,方程2表示观察到一个特殊图像亮度的概率。
P(Γ i)表示在组织分类上的一个固定的先验概率分布。若这个概率在这些组织分类上是均匀分布的,就可以用最大似然值逼近这个组织分类的分量。Kamber等 [38]描述了在大脑组织分类的空间变化的先验概率密度。这样结构对使用这样的模型也许很有利。
整个偏倚信息变量记为 β=( β 0β 1,…, β n -1T,其中 n是体素的个数。这个偏倚信息变量为 n维均数为0的正态分布先验概率密度模型。这个模型允许我们捕获出现在那些亮度不均匀性的光滑表面。
Pβ)= G ψββ
(6-6)
其中
n维正态分布密度函数。 ψ β是整个偏倚信息变量 βn× n协方差矩阵。虽然 ψ β太多的元素以致在实际中无法直接利用。但当 ψ β选定后,结合 ψ β可以得到较容易处理的估计量。
假定偏倚变量和组织分类是独立的,则要考虑亮度的不均匀性的原因是否来自设备。利用亮度和组织分类的联合概率和条件概率关系可以得到
(6-7)
以及在计算出整个组织分类的边缘概率后得到亮度的条件概率
(6-8)
这个表达式可以简洁地改写为
(6-9)
其中加权数W ij的定义如下:
(6-10)
式中的下标 ij分别对应体素的顺序编号和组织的分类编号,以及第 j个组织分类的亮度均数为
μ jμ(第 j个组织分类)
而残差均数定义为
(6-11)
协方差倒数的平均数为
(6-12)
这一节中统计建模的结果可以具体归结为一个非线性最优化的问题并可用下列公式表示偏倚变量的估计问题:
(6-13)
这个结果表明:最优化的解取决于亮度观察值与每个组织分类的亮度均数之间的残差均数、各组织分类的亮度协方差均数和偏倚变量的协方差。
对非线性估计量方程(6-9)用期望-最大(expectation-maximization,EM)算法得到偏倚变量估计。交替使用公式(6-10)至(6-13)构成了EM算法的迭代序列。
(6-14)
(6-15)
换而言之,给出偏倚变量的估计值的情况下,用(6-14)估计权重 w ij;反过来,在给出权重 w ij的情况下,用(6-15)估计偏倚信息变量 ,两个方程交替使用构成迭代序列直至收敛。
自适应分割方法可以应用到旋转回波图像和梯度回波图像。举例所示图像为冠状面(3DFT梯度回波T1-加权)图像。在这一节所显示的所有MR图像都是使用General Electric Signa 1.5 Tesla临床MR图像器(General Electric Medical Systems,Milwaukee,WI)获得。在第3节所描述的各向异性扩张滤波器被用于信号预处理阶段以达到降低噪声的目的。
图6-5(a)显示输入的图像是冠状面(3DFT梯度回波T1-加权)中的一个薄切片。脑组织的ROI是由人工产生的。图6-5(b)显示最终偏倚信息变量的收敛估计值。当偏倚校正的最大值与最小值的差为10左右,输入数据的最大值是85。图6-5(c)显示由自适应分割方法所得的分割结果。
图6-5 利用期望-最大分割
注意到在右临时区域的显著地改善。在初始的分割时,白质在二值化的图像中完全消失。
二、无种子区域的增长
无种子区域增长(unseeded region growing,URG)类似于种子区域增长,只不过无需选择种子:种子由分割过程中自动产生。因此这个方法实现了在区域基础上的自动分割且有很强的抗干扰性。
在分割初始过程中预置一些仅包含一个图像素的区域A 1,以后的连续分割过程状况由 n个等同性质的区域 A 1A 2,…, A n构成的一个集合组成。设T是由所有未分配的图像素构成的一个集合,而且 T中的任一元素至少在其中一个区域的边界上。
其中 Nx)是点 x周围最邻近的像素构成的集合。差分测度定义为:
其中,记 gx)为在点 x的图像值, iNx)与 A i.相交区域的下标索引。
增长过程包括了一个点 zT和区域 A i,其中 i∈[1, n]并且满足
δzA i)小于预定的阈值 t,则这个像素加入到 A i,反之,我们必须选择在大体上最相似的区域A,且满足
δzA)﹤ t,则我们指定这个像素到 A。如果上述两个条件都没有满足,则显然这个像素显著地不同于当前已找到的所有区域。因此认定一个新的区域 A n +1,且归入已找到的区域之列,并将点 z作为其初始值。上述三种情况,一旦像素加入这个区域,被指定区域的统计量必须更新。
URG分割过程是迭代过程,重复上述过程直到所有像素都分配到区域中为止。为了确保符合均匀性准则的正确操作,当每次区域统计量变化时,区域增长操作需要确定最佳像素。详细的运作可以参照文献 [39]。分割结果在图6-6显示。
图6-6 无种子区域增长分割
第五节 利用3维Monte Carlo模拟改进图像配准的置信区间
临床诊断和治疗通常需要多种配准图像。大多数的医学图像配准方法 [40,41,42]以某个损失函数的最大化或最小化作为全局的优化匹配。这些函数通常是二个图像集合上一定同质特征的距离平方和。例如,两个图像集合的同质点对之间的距离平方和 [43],头帽方法中头部CT,MR和PET图像的皮肤表面之间的距离 [44],PET图像的像素值与MR图像的模拟像素值的绝对差值 [45],像素值与它们在同类组织中的均数之比 [46]等。然而,大多数的损失函数不直接反映目标的实际位置和估计位置之间的距离,即目标配准误差(target registration error,TRE)。大多数医学应用需要评估准确性和精密性的方法来判断他们的结果。Woods等 [46]提出用内部一致性度量为MRI数据定位准确性设置标准。几乎所有的配准准确度评估方法可以分为二大类:由视觉检查进行定性评估和参照金标准进行定量评估。前者需要特殊的专家和广泛的经验,而后者要求非常准确的金标准,但不易实现。在同一个准则下,不同的方法常常是不可比较的。
借用非线性回归分析的思想 [47],我们可以把图像配准的问题归结为用非线性最小二乘法通过一个图像集合(视为函数)到另一个图像集合(视为数据)的最优拟合进行变换参数估计。对于最小二乘估计方法,可以假定损失函数在当前参数值附近的邻域内是线性的。因此我们可以使用下列方程式 [47]计算置信区间或区域:
θ- θ 0)∑( f′)≤( σ 2n-1)) Fpn- p,1- α
(6-16)
其中 F是在 F检验中置信水平 α所对应的临界值, σ 2n-1)是参数估计值所在的位置的残差平方和(配准损失函数),Σ( f′)表示参照模型图像关于变换参数的导数之和。 θ是对应于置信水平的参数, θ 0是由配准过程中所产生的最优化参数值。
因为(6-16)所用到的所有数据点应该是相互独立的,然而图像上的数据点是相关的,图像上的数据点总数不能作为样本量 n,所以独立数据点的有效个数需要进行估计。
为了确定在置信区间估计中所用的独立数据点的有效个数,我们首先用一个正态条件下的Monte Carlo模拟研究。根据模拟结果选用相同的 n作为95%和90%的置信区间估计中独立数据点的有效个数。在这个研究的其他部分中的各种模拟条件下,我们进一步考查了选择 n的有效性。
用Monte Carlo方法模拟2维PET图像,接着进行模拟图像的配准处理。用变换参数估计的分布来评估90%,95%和99%置信区间和参数空间分布的一致性。在重建图像前先查看两个图像灰白质比不相符合是否影响置信区间,用2∶1,3∶1和4∶1的灰质与白质的比对2维Hoffman脑模型的2维灰质和白质窦腔照片分割图 [48]进行组合。然后使用含有各种滤波器(如,Hanning,Ramp,Butter-worth,Ham,Parzen和Shepp-Logan滤波器)的向后投影重建程序,重建128×128像素的图像。采用了各种空间移位量(即旋转0.3、0.8、1.2和3.3度,平移0.16、0.8、1.6和2.4毫米)。模拟了各种水平的泊松分布噪声(即总计数为5×10 5,1×10 6和2×10 6)。在配准前,对两个图像集合均应用了5mm的FWHM和高斯光滑滤波器。选用Powell的算法 [49]作为最优化方法在极端噪声条件和大对比度差异的情况中,残差平方和(RSS)由两部分组成:系统的误差和由于统计噪声的误差:
RSS=RSS system+RSS noise
(6-17)
系统误差的形成是由于两图像之间的原始差异、不恰当的配准方法和程序的精度错误等。这种误差独立于最初的移位和噪声。残差平方和的第二部分(RSS noise)是由于统计噪声。假如系统误差相对大于噪声项,即噪声信号非常低以及不相符的灰白比相当高的情况下,则残差平方和的估计中的系统误差部分需要进行校正。
因为在RSS中的系统成分对空间平滑的敏感性远低于(6-17)中其他成分对空间平滑的敏感性,当参数找到时,通过把平滑过滤器应用于两组相对大FWHMs的图像来估计RSS的系统成分。除去系统成分后,(6-17)中的RSS提供噪声成分的估计。
在统计的回归分析中计算的置信区间与图像定位的变换参数样本分布的模拟结果是一致的。移位、重建处理、噪声程度或示踪物分布的变化对计算置信区间的有效性没有什么影响。在校正了估计的残差平方和中的系统误差后,即使噪声非常大和两个图像集合有大的分布差异,置信区间也能正确地计算。因为多度量配准可以视为一个图像集合与另一个从其他图像形态所得的模拟图像的单度量配准,这个方法也可以适用于多度量配准。因此,评估配准结果的精度不需要专家的目测和核实。结果显示,使用统计置信区间有可能提供单个图像配准自动和客观的评估。
第六节 结论
我们尝试简要地介绍了统计方法在一般图像处理中的应用,特别是在医学图像中的应用。我们讨论的内容包括了图像采样、压缩、滤波、分割和配准。在理论上讨论了方法以及用结果验证了方法。在许多信号处理应用中,统计方法是有力的工具。我们希望这些概况能为在图像处理中进一步使用统计方法提供一些借鉴。
致谢
本章第一节得到龚高全医生的帮助,特表谢意。
参考文献
1.Huang SC,Phelps ME,Hoffman EJ,et al. Non-invasive determination of local cerebral metabolic rate of glucose in man,American Journal of Physiology,1980,238,E69-E82.
2.Patlak CS,Blasberg RG,Fenstermacher J. Graphical evaluation of blood to brain transfer constants from multiple-time uptake data. Journal of Cerebral Blood Flow and Metabolism,1983,3,1-7.
3.Feng D,Ho D,Chen K,et al. An evaluation of the algorithms for constructing local cerebral metabolic rates of glucose tomographical maps using positron emission tomography dynamic date,IEEE Trans. on Medical Imaging,1995,14(4),697-710.
4.Patlak CS,Blasberg RG. Graphical evaluation of blood to brain transfer constaints from multiple-time uptake data generalizations. Journal of Cerebral Blood Flow and Metabolism,1985,5:584-590.
5.Cai W,Feng D,Fulton R. Clinical investigation of a knowledge-based data compression algorithm for dynamic neurologic FDG-PET images. Proceedings of the 20th Annual International Conference of the IEEE Engineering in Medicine and Biology Society(EMBS’98),1998,vol.20,part3,pp1270-1273,Hong Kong,Oct 29-Nov 1.
6.Li X,Feng D,Chen K. Optimal image sampling schedule:A new effective way to reduce dynamic image storage space and functional image processing time,IEEE Trans. on Medical Imaging,1996,15,710-718.
7.Cobelli C,Ruggeri A,DiStefano Ⅲ J J,et al. Optimal design of multioutput sampling schedules-software and applications to endocrine-metabolic and pharmacokinetic models,IEEE Trans. on Biomedical Engineering,1985,32(4),249-256.
8.Ciaccio EJ,Dunn SM,Akay M. Biosignal pattern recognition and interpretation systems:Methods o classification. IEEE Engineering in Medicine and Biology,1994,13,129-135.
9.Crocker LD. PGN:The portable network graphic format. Dr. Dobb’s Journal,1995,36-49,July.
10.Perona P,Malik J. Scale-space and edge detection using anisotropic diffusion. Proc. IEEE Workshop Computer Vision,Miami,1987,16-22.
11.Perona P,Malik J. Scale-space and edge detection using anisotropic diffusion. IEEE Trans,1990,PAMI 12:629-639.
12.Alvarez L,Mazorra L. Signal and image restoration using shock filters and anisotropic diffusion. SIAM Journal on Numerical Analysis,1994,31(2):590-594.
13.Tsai D M,Chen Y. A fast histogram-clustering approach for multi-level thresholding. Pattern Recognition Letter,1992,13:245-252.
14.Bhandari D,Pal NR,Majumder D D. Fuzzy divergence,probability measure of fuzzy events and image thresholding. Pattern Recognition Letter,1992,13:857-867.
15.Luijendijk H. Automatic threshold selection using histograms based on the count of 4-connected regions. Pattern Recognition Letter,1991,12:219-228.
16.Tseng D C,Huang M Y. Automatic thresholding based on human visual perception. Image and Vision Computing,1993,11:539-548.
17.Nagawa Y,Rosenfeld A. Some experiments on variable thresholding. Pattern Recognition,1979,11:191-204.
18.Cho S,Haralick R M,Yi S. Improvement of Kittle and Illingworth’s minimum error thresholding. Pattern Recognition,1989,22:609-617.
19.Glaseby C A. An analysis of histogram based thresholding algorithms. CVGIP:Graphical Models and Image Processing,1993,55:532-533.
20.Rodriguez A A,Mitchell O R. Image segmentation by successive background extraction. Pattern Recognition,1991,24:409-420.
21.Parker J R. Gray level thresholding in badly illuminated images,IEEE Trans. PAMI,1991,13:813-819.
22.Spann M,Horne C. Image segmentation using a dynamic thresholding pyramid. Pattern Recognition,1989,22:719-732.
23.Donoho D L. De-noising by soft-thresholding. IEEE Trans. Info,1995,Theory IT-41,613-627.
24.Ross S M. Introduction to Probability and Statistics for Engineers and Scientists. New York:John Wiley and Sons,1987.
25.Tou J T,Gonzalez R C. Pattern Recognition Principle. Addison-Wesley,1972.
26.Haralick R M,Shapiro L G. Image segmentation techniques. Computer Graphics Image Processing,1985,29:100-132.
27.Sahoo P K,Soltani S,Wong A K C. A survey of threhsolding techniques. Computer Graphics Image Processing,1988,41:230-260.
28.Davis L S. A survey of edge detection techniques. Computer Graphics Image Processing,1975,4:248-270.
29.Horowitz S L,Pavlidis T. Picture segmentation by a directed split-and-merge procedure,Proc. 2 nd Int. Joint Conf. On Pattern Recognition,1974,424-433.
30.Adams R,Bischof L. Seeded region growing,IEEE Trans Pattern Anal Mach Intel,1994,16(式6):641-647.
31.Meyer F,Beucher S. Morphological segmentation. Journal of Visual Communication And Image Representation,1990,1:21-46.
32.Besl P J,Jain R C. Segmentation through variable-order surface fitting. IEEE Trans Pattern Anal Mach Intel,1988,10(2):167-192.
33.Kass M,Witkin A,Terzonpoulos D. Snakes:active contour models. London:Proc Int Conf On Computer Vision,1987.
34.Vannier M,Butterfield R,Jordan D,et al. Multi-spectral analysis of magnetic resonance images. Radiology,1985,154:221-224.
35.Kohn M,Tanna N,Herman G.,et al. Analysis of brain and cerebrospinal fluid volumes with MR imaging. Radiology,1991,178:115-122.
36.Gerig G.,Kuoni W,Kikinis R,et al. Medical imaging and computer vision:an integrated approach for diagnosis and planning. Proc 11’th DAGM Symposium,1989,425-443.
37.Cline HE,Lorensen WE,Kikinis R,et al Three-dimensional segmentation of MR images of the head using probability and connectivity. JCAT,1990,14(式6):1037-1045.
38.Kamber M,Collins D,Shinghal R,et al. Model-based 3D segmentation of multiple sclerosis lesions in dual-echo MRI data. SPIE Vol. 1808,Visualization in Biomedical Computing,1992.
39.Lin J Z,Jin J S,Hugo T. Unseeded region growing,Proc. Sydney:Workshop on Visual Information Processing 2000,2001.
40.Van den Elsen P A,Pol E. J. D,Viergever M A. Medical image matching-a review with classification. IEEE Eng Med Biol,1993,12:26-39.
41.Maurer CR,Fitzpatrick JM. A review of medical image registration,In Maciunas,R. J.(ed.),Interactive Imageguided Neurosurgery. Parkridge:American Association of Neurological Surgeons,1993,17-44.
42.Maintz J. B. A,Viergever M A. A survey of medical image registration. Medical Image Analysis,1998,2(1):1-36.
43.Evans A C,Marrett S,Collins L,et al. Anatomical-functional correlative analysis of the human brain using three dimensional imaging systems,In Schneider,R. H.,Dwyer Ⅲ,S. J. and Jost,R. G.(eds),Medical Imaging:Image Processing,SPIE Press,Bellingham,WA.,vol.1092,pp. 264-274.
44.Pelizzari C A,Chen GTY. Accurate three-dimensional registration of CT,PET and/or MR images of the brain. J Comput Assis Tomogr,1989,13:20-26.
45.Lin KP,Huang SC. A general technique for interstudy registration of multifunction and multimodality images. IEEE Tran Nucl Sci,1994,41(式6):2850-55.
46.Woods RP,Grafton S T,Holmes C J,et al. Automated image registration,I. General methods and intrasubject,intramodality validation,J. Comp. Assist. Tomogr,1998,22(1):139-152.
47.Draper N R. Applied regression analysis. 2nd Ed. New York:Wiley,1981.
48.Hoffman E J,Cutler P D,Guerrero T M,et al. Assessment of accuracy of PET utilizing a 3-D phantom to simulate the activity distribution of[18F]fluorodeoxyglucose uptake in the human brain. J. Cereb. Blood Flow Metab,1991,11:17-25.
49.Powell M J D. An efficient method for finding the minimum of a function of several variables without calculating derivatives. Comput J,1964,7:155-163.
50.Wells III WM,Grimson W E L,Kikinis R,et al. Adaptive segmentation of MRI data. IEEE Transactions on Medical Imaging,1996,15(4):429-443.
作者简介
金声,博士,1982年毕业于上海交通大学,后获新西兰奥塔戈大学博士。曾任澳大利亚悉尼大学计算机系副教授及研究生委员会主任、新南威尔士大学客座副教授及可视信息处理实验室主任、纽卡斯尔大学教授。他是多媒体技术和可视信息检索及处理方面的国际知名专家。目前已发表百余篇学术论文、出版多部专著。他创办一个高科技公司并获澳大利亚科技园1999年最佳新企业奖。他同时是多家公司(包括莫托罗拉,CA,扫描世界,超软等)的顾问。他曾任麻省理工学院、加州大学洛杉矶分校、香港理工大学和清华大学访问教授、澳华科技协会副会长。
赵耐青,教授,现为复旦大学公共卫生学院生物统计学教研室主任、博士生导师,中国医学数学会副理事长、中国卫生信息学会常务理事,中国卫生信息学会卫生统计教学专业委员会副主任委员;中国卫生信息学会统计理论与方法专业委员会常务委员;现场统计研究会生物医学统计分会副理事长。1983年获复旦大学数学系学士学位,1996年获澳大利亚Newcastle University生物统计硕士学位。主要研究方向为:临床流行病学和流行病学中的卫生统计问题,时间序列分析。主编教育部十五规划教材《医学统计学》,发表论文数十篇。