CN102282587A - 用于处理图像的系统和方法 - Google Patents

Abstract

一种处理多个时间分离的图像的方法，包括：在各个图像中选择多个成像单元；测量在各个成像单元中的时域差；以及选择在阈限以上的时域差。

Description

用于处理图像的系统和方法

相关申请的交叉引用

本申请要求于2008年10月2日提交的、名为“用于处理图像的系统和方法(system and method for processing images)”的、李(Lee)的美国临时申请No.61/102,206的权益，该美国临时申请的内容通过引用而并入本申请中。

技术领域

本发明总体上涉及图像处理，并且更具体地，涉及处理通过使用动态成像技术而获得的图像数据的系统和方法。

背景技术

医学成像涵盖用于为了临床目的创建人体的图像的技术和过程，包括用于诊断或监视疾病的医疗程序。医学成像技术已经发展为涵盖很多图像记录技术，所述图像记录技术包括电子显微镜法、荧光透视法、磁共振成像(MRI)、核医学、光声成像、正电子发射断层摄影术(PET)、投影射线摄影术、热成像法、计算机断层摄影术(CT)和超声。

医学成像可以包含，使用被称为造影剂或造影材料的复合物，以提高图像中体内结构的可见度。

医学成像原本局限于采集单个静态图像以捕获身体的器官/区域的解剖结构。目前，更复杂精密的成像技术的使用允许进行动态研究，所述动态研究提供可以表征生理或病理生理信息的图像的时域序列。

动态医学成像涉及获取过程，即，随时间的推移而拍摄受关注的器官/区域/身体的许多“快照(snapshot)”，以便捕捉随时间变化的行为，例如，造影剂的分布，且因而捕捉具体的生物状态(疾病、不适、生理现象等)。随着医学成像的速度和数字性的发展，此获取数据会具有极大的瞬时清晰度且会导致大量数据。

医学成像技术已广泛用于改善对诸如癌症、心脏疾病、大脑失调和心血管疾病的病情的诊断和护理。大部分评估表明，数以百万的生命由于这些医学成像技术而得以挽救或者显著得到改善。然而，必须对这类医学成像技术给患者带来的辐射暴露风险加以考虑。

在这点上，由于CT中所递送的辐射剂量比更为常见的常规x射线成像过程大，医学诊断中CT越来越多的使用已经对被暴露的患者增加的罹患癌症风险进行了突出的关注。对于两相CT灌注成像方案而言，则尤其如此。为了说明目的，将剂量-长度乘积用于由可市购的CT扫描仪(GE Healthcare)提供的方案，CT中风系列的有效剂量为10mSv，其中4.9mSv是由两相CT灌注成像方案贡献的，其中所述CT中风系列由以下组成：1)非增强CT扫描，以排除出血；2)CT血管造影术，定位导致中风的梗塞；和3)两相CT灌注成像方案，以用血流量和血容量来定义缺血区，并用血脑屏障渗透表面积乘积(BBB-PS)来预测出血性转变(HT)。预测CT中风系列对于每暴露的100,000个男性和女性急性缺血性中风(AIS)患者分别会引发80和131例额外癌症，并且仅两相CT灌注成像方案就被预测为导致了一半的额外癌症(Health Risks from Exposure to Low Levels of Ionizing Radiation：BEIRVII.The National Academies Press，Washington DC，2006(暴露于低水平电离辐射的健康风险BEIR VII.国家学术出版社华盛顿特区2006))。

除非减少两相CT灌注成像方案的有效剂量，否则医学成像的益处将会由于对癌症引发的担忧而被破坏。此外，对辐射暴露风险的担忧延伸到其他医学成像技术，尤其是对于儿科患者或正经受重复暴露的那些患者。

例如，在2007年3月6日授权的第7,187,794号美国专利中，描述了使用统计滤波技术以增大低剂量辐射成像中的信噪比。然而，目前，统计滤波的应用限于图像构建之前的投影数据。若干不利之处与投影数据的统计滤波相关。例如，因为来自动态序列的各个图像一般是由～1,000个投影重建而成，所以使用投影数据工作时计算负担高。另外，滤波过程必须被“硬布线”至扫描仪的图像重建管线中。如此一来，投影数据的统计滤波器由于其通常专用于扫描仪平台而缺乏灵活性。

不论是否已被觉察或证实，对辐射暴露风险的担忧将必须得到解决，以使患者对医学成像技术的长期安全放心。据此，需要允许辐射剂量的减少的医学成像技术。

因此，本发明的目的是要提供一种新颖的用于处理图像的系统和方法。

发明内容

在一方面，提供有一种处理多个时间分离的图像的方法，包括：在各个图像中选择多个成像单元；为各个成像单元确定时域差；以及选择在阈限以上的时域差。

在另一方面，提供有一种处理多个时间分离的图像的系统，包括：接口，用于接收多个时间分离的图像；和处理器，用于在各个图像中选择多个成像单元，为各个成像单元确定时域差，以及选择在阈限以上的时域差。

在又一方面，提供有一种包含处理多个时间分离的图像的计算机程序的计算机可读介质，所述计算机程序包括：用于在各个图像中选择多个成像单元的计算机程序代码；用于为各个成像单元确定时域差的计算机程序代码；和用于选择在阈限以上的时域差的计算机程序代码。

附图说明

现在，将仅通过示例的方式、参照附图来描述实施例，在附图中：

图1是图像处理方法的流程图；

图2示意性图示了多个时间分离的图像中的成像单元的定位；

图3A至3E示出了用统计滤波技术处理过的、来自CT图像的脑血流量图；

图4图示性示出了使用图3E中勾画出的区域为图3A至3D中所示的图像所计算的血流量品质因数(Figure of Merit)；以及

图5是用于实施图1中所示图像处理方法的系统的图示。

具体实施方式

动态医学成像涉及采集过程，即，随时间的推移而拍摄受关注的器官/区域/身体(即靶区)的许多“快照”，以便捕捉随时间变化的行为，例如造影剂的吸取和/或消除。此采集过程导致多个时间分离的图像的产生。本申请中所描述的方法和系统涉及这样的时间分离的图像的处理。

现在转到图1，图示出了根据该方法在时间分离的图像的处理期间所执行的步骤。多个时间分离的图像可以在使用或者不使用注射的造影剂的情况下从患者的动态医学成像获得。造影剂可以例如根据靶区的结构或生理状态，来选择性地增加图像中靶区的对比度。例如，可以将对特定器官、疾病、状态或生理过程具有生物物理、分子、基因或细胞亲和性的复合物注射到患者体内。选择这样的造影剂，使之具有下述属性：通过改变成像条件，例如，通过改变图像对比度，对给定成像技术提供增强的信息，以反映体内复合物的行为。这可以通过局域性部位处增大的X射线衰减(例如对于CT/X射线)、被改变的顺磁特性(例如对于MRI)或者对于核医疗/PET使用放射性同位素来实现。用于许多成像技术的造影剂是公知的。在一些情况下，对比增强并不依赖于所注入的造影剂，所述一些情况例如：在磁共振成像(MRI)中使用“黑血”或“白血”，其中，使用特定脉冲序列以改变血的磁饱和且因而改变其在图像中的呈现；或者使用带标记的MRI序列，所述带标记的MRI序列改变具体组织或流体的磁性能。所注入的造影剂或者具有改变了的性能的组织或流体可以全部被认为是“成像剂”。

多个时间分离的图像不限于任何具体的成像技术，并且包括例如使用磁共振成像(MRI)、计算机断层摄影术(CT)、核医学(NM)或正电子发射断层摄影术(PET)的动态医学成像。多个时间分离的图像也不限于具体图像类型。例如，可以处理灰度或彩色图像。

在该方法期间，在各个时间分离的图像内，选择多个成像单元(步骤110)。“成像单元”可以为用于将图像分成等同的部分的任何所期望的单元，包括例如像素、多个像素、部分像素、体素、多个体素、部分体素等。

成像单元数据由例如数字值的值来表示，并且因而可以量化和测量。一般而言，为时间分离的图像中的对应的成像单元测量图像强度。对于使用造影剂的动态成像而言，可以关于造影剂浓度确定时域差。造影剂浓度的时域差可以被表示为作为时间的函数的对比增强的图示。

图1中示出了用于确定造影剂浓度的时域差并选择时域差的两种备选方案。在图1所示的一个备选方案中，一旦在步骤110从各个时间分离的图像中选择了多个成像单元，就为各个成像单元确定时域差(步骤115)。然后根据所观察到的时域差的程度来对所确定的时域差或其表示进行排序(步骤120)。然后选择在阈限以上的时域差(步骤125)。然后基于所选的时域差来构建处理后的图像(步骤150)。

在图1所示的另一备选方案中，不对时域差排序。而是，利用统计技术，且更具体地，利用盲信号分离统计技术，来分析成像单元的时域数据，以选择在阈限以上的时域差。例如，关于使用造影剂的动态CT成像，一旦在步骤110已经从各个时间分离的图像中选出了多个成像单元，就建立与平均曲线相关的对比增强图示的所有对的协变的协方差矩阵(步骤130)并使用盲信号分离统计技术进行分析(步骤135)，诸如对协方差矩阵应用主分量分析。基于统计分析来计算特征向量/特征值对(步骤140)。然后所得到的在阈限以上的特征向量/特征值对基于它们的特征值被选择出来(步骤145)。然后利用所选的特征向量来构建处理后的图像，并且使用所选特征向量的不同权重来构建各成像单元(步骤160)。

可选地，可以配准多个时间分离的图像，以使时间分离的图像中的对应位置彼此相关。仅需要多个时间分离的图像的配准来校正显著失配。否则，在时间分离的图像之间的轻微移动可以为本申请中所描述的处理方法所容许甚至可以被消除。

作为另一可选步骤，可以对多个时间分离的图像进行预处理，使得时间分离的图像由可量化的值或代码表示。例如，可以利用已知数字化技术而由数字值表示时间分离的图像。如果多个时间分离的图像的成像数据是以数字形式提供的，则无需数字化步骤。

现在将参照图2来进一步描述造影剂浓度的时域差的确定。图2示出了一系列6个按照时间顺序排列的时间分离的图像t₁至t₆。在各个时间分离的图像上标记有成像单元的4×6阵列。在此示例中，选择了五(5)个成像单元IU1至IU5。如从比较整个系列的时间分离的图像t₁至t₆清楚可知，对于各个时间分离的图像在阵列的相同位置(即第1列第2行)处选择成像单元IU1。相似地，贯穿该系列时间分离的图像在成像单元阵列的对应位置处选择成像单元IU2至IU5中的各个。据此，可以将贯穿时间分离的图像在其对应位置处的各个成像单元的值绘制为时间的函数，以获得时域曲线。绘制步骤由连接成像单元IU1、IU2和IU5的对应位置的箭头图示出。然后可以对各个成像单元的时域曲线进行分析，以确定时域差。例如，可以计算所获得的全部时域曲线的平均曲线。然后从各个时域曲线减去平均曲线，以获得各个所选成像单元的时域差曲线。将予以理解的是，所标记的阵列和成像单元IU1至IU5的选择、以及表示将成像单元IU1、IU2和IU5的值作为时间的函数进行绘图的箭头，仅用于说明目的，并且可以按照需要选择任何数目的成像单元并对其进行分析。此外，通过从成像单元的时域曲线减去平均曲线而进行的时域差的确定仅用于说明目的，并且可以想到分析时域曲线的任何其他方法，例如使用导数和/或统计技术。

现在，将使用来自CT灌注研究的时间分离的图像作为示例，提供用于多个时间分离的图像的盲信号分离统计分析的数学基础。

给定来自CT灌注研究的一系列n个动态(即时间分离的)图像，n个动态图像可以紧凑地表示为矩阵：

$\tilde{X}={\left[\begin{array}{cccccc}{\tilde{x}}_1^T& {\tilde{x}}_2^T& \·\·\·& \·\·\·& {\tilde{x}}_{p-1}^T& {\tilde{x}}_p^T\end{array}\right]}^T$

其中，p为CT图像中的像素数(即p＝512×512、256×256等)，且

为n×1向量，其元素为n个图像中的像素值或来自第i个像素的时间-增强曲线。注意，

的第(i，j)个元素为

的第j个元素x_ij。各个

也可以解释为单个时间-增强曲线

(由n个随机变量χ₁，χ₂，…，χ_n-1，χ_n组成的随机过程)的重复观察。

符号解释如下：

(i)向量由具有排成一列的元素的希腊或罗马字母字符上方的颚化符号(～，tilde)来表示；

(ii)矩阵由大写的希腊或罗马字母字符上方的颚化符号来表示；

(iii)向量或矩阵的转置由符号上的上标T来表示；以及

(iv)随机变量由希腊字母字符表示，而随机变量的观察值(样本)由对应的罗马字母字符表示；希腊和罗马字母字符之间的对应关系将在下列讨论中如下所示使用：

$\mathrm{\χ}\↔x$

$\mathrm{\ψ}\↔y$

$\mathrm{\ζ}\↔z$

(v)由n个随机变量组成的随机过程由n×1向量表示。

从n个动态图像中消除噪声的一个途径是，找到p×n矩阵

的“展开式(expansion)”，并且然后按照某个预定标准截短展开式。常用的展开式为的奇异值分解：

$\tilde{X}=\tilde{U}\·\tilde{L}\·{\tilde{V}}^T---\left(1\right)$

其中，

为p×p正交矩阵，

为n×n正交矩阵，且

为p×n对角矩阵并将

的奇异值作为其对角元素。令l₁，…，l_n为

的奇异值(在不失一般性的情况下假定

具有满列秩n，且p＞n)。则：

其中，

为p×1向量且其元素除第i个元素为l_i以外均为零。通常，奇异值l₁，…，l_n按降序排列，并且小奇异值由于其是由噪声引起的因而可以被丢弃。减少

中的噪声的简单方式是只保留r＜n个奇异值。则：

如果截短后的奇异值序列被用以根据等式(1)来重构

则：

$\tilde{X}\≈\left[\begin{array}{cccc}{\tilde{u}}_1& {\tilde{u}}_2& \·\·\·& {\tilde{u}}_p\end{array}\right]\·\left[\begin{array}{cccccc}{\tilde{l}}_1& \·\·\·& {\tilde{l}}_r& {\tilde{0}}_{r+1}& \·\·\·& {\tilde{0}}_n\end{array}\right]\·\left[\begin{array}{c}{\tilde{v}}_1^T\\ \·\\ \·\\ \·\\ {\tilde{v}}_r^T\\ {\tilde{v}}_{r+1}^T\\ \·\\ \·\\ \·\\ {\tilde{v}}_n^T\end{array}\right]$

其中，

为

的第i列，且

为

的第i列。估算展开式，

得到：

$\tilde{X}\≈\left[\begin{array}{cccc}{\tilde{u}}_1& {\tilde{u}}_2& \·\·\·& {\tilde{u}}_p\end{array}\right]\·\left[\begin{array}{cccccc}{\tilde{l}}_1& \·\·\·& {\tilde{l}}_r& {\tilde{0}}_{r+1}& \·\·\·& {\tilde{0}}_n\end{array}\right]\·\left[\begin{array}{c}{\tilde{v}}_1^T\\ \·\\ \·\\ \·\\ {\tilde{v}}_r^T\\ {\tilde{v}}_{r+1}^T\\ \·\\ \·\\ \·\\ {\tilde{v}}_n^T\end{array}\right]$

$\≈\left[\begin{array}{cccc}{\tilde{u}}_1& {\tilde{u}}_2& \·\·\·& {\tilde{u}}_p\end{array}\right]\·[{\tilde{l}}_1\·{\tilde{v}}_1^T+\·\·\·+{\tilde{1}}_r\·{\tilde{v}}_r^T]\≈\left[\begin{array}{cccc}{\tilde{u}}_1& {\tilde{u}}_2& \·\·\·& {\tilde{u}}_p\end{array}\right]\·(\left[\begin{array}{c}{l}_1\\ \·\\ \·\\ \·\\ 0_r\\ 0_{r+1}\\ \·\\ \·\\ \·\\ 0_p\end{array}\right]\·{\tilde{v}}_1^T+\·\·\·+\left[\begin{array}{c}0\\ \·\\ \·\\ \·\\ l_r\\ 0_{r+1}\\ \·\\ \·\\ \·\\ 0_p\end{array}\right]\·{\tilde{v}}_r^T)$

$\≈\left[\begin{array}{cccc}{\tilde{u}}_1& {\tilde{u}}_2& \·\·\·& {\tilde{u}}_p\end{array}\right]\·\left[\begin{array}{c}{l}_1\·{\tilde{v}}_1^T\\ \·\\ \·\\ \·\\ l_r\·{\tilde{v}}_r^T\\ 0_{r+1}\\ \·\\ \·\\ \·\\ 0_p\end{array}\right]=\underset{l=1}{\overset{r}{\mathrm{\Σ}}}{l}_i\·{\tilde{u}}_i\·{\tilde{v}}_i^T---\left(2\right)$

由此，在

的展开式中，保留r个奇异值意味着只需要

和

的前r列。

然而，如这里用公式表达的

的奇异值分解(SVD)(展开式)在对来自CT灌注的动态图像去噪时存在两个主要缺点，即：

1.矩阵的大小可以大至(256)²×40-60，从而使得在计算上难以找到

的SVD；并且

2.虽然丢弃小奇异值的标准具有一般意义，但难以验证用以保留或丢弃奇异值的阈值。

现在将描述用以获得

的SVD的替代方法。如上所述，来自各个像素(或像素块)的时间-增强曲线

可以视为底层(underlying)随机过程(时间-增强曲线)

的重复样本(观察值)。则随机过程的样本均值

可以计算为：

$\stackrel{\‾}{x}=\frac{1}{p}\·\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{\tilde{x}}_i---\left(3\right)$

按照定义为n×1向量(样本均值的时间-增强曲线)，且根据等式(3)，的第j个分量为：

其中

为

的第j个分量

来自

的零均值随机过程为

其中

为

的总体均值。对于

的各个观察值(样本) 的对应的观察值(样本)由下式给出：

${\tilde{y}}_i={\tilde{x}}_i-\stackrel{\‾}{x}.$

的样本均值

(一个n×1向量)由下式给出：

$\stackrel{\‾}{y}=\frac{1}{p}\·\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{\tilde{y}}_i=\frac{1}{p}\·\underset{i=1}{\overset{p}{\mathrm{\Σ}}}({\tilde{x}}_i-\stackrel{\‾}{x})=(\frac{1}{p}\·\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{\tilde{x}}_i)-\stackrel{\‾}{x}=\tilde{0}$

其中，由于为零均值随机过程，所以如所预期，

为n×1零向量。

令

为p×n矩阵，其各行为 ${\tilde{y}}_i^{T}i=1,\·\·\·,p.$ 则：

${\tilde{Y}}^T\·\tilde{Y}=\left[\begin{array}{cccccc}{\tilde{y}}_1& {\tilde{y}}_2& \·\·\·& \·\·\·& {\tilde{y}}_{p-1}& {\tilde{y}}_p\end{array}\right]\·\left[\begin{array}{c}{\tilde{y}}_1^T\\ {\tilde{y}}_2^T\\ \·\\ \·\\ \·\\ \·\\ \·\\ \·\\ {\tilde{y}}_{p-1}^T\\ {\tilde{y}}_p^T\end{array}\right]=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{\tilde{y}}_i\·{\tilde{y}}_i^T$

其中，

为ψ_i i＝1，…，n的样本方差，且

为ψ_i和ψ_j i，j＝1，…，n的样本协方差。

或者，

为

的样本协方差矩阵并且具有n×n的维数。在CT灌注中n一般为40-60。

ψ_i，i＝1，…，n，的线性组合ζ可以定义为：

$\mathrm{\ζ}={\tilde{a}}^T\·\widetilde{\mathrm{\ψ}}$

其中，

为给定的n×1系数向量。ζ为随机变量(非过程)，且ζ的观察值z_i i＝1，…，p按照ζ的定义而由下式给出：

$z_i={\tilde{a}}^T\·{\tilde{y}}_i$

ζ的样本均值

由下式给出：

$\stackrel{\‾}{z}=\frac{1}{p}\·\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{\tilde{z}}_i=\frac{1}{p}\·\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{\tilde{a}}^T\·{\tilde{y}}_i={\tilde{a}}^T\·\frac{1}{p}\·\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{\tilde{y}}_i=0$

由此，ζ也为零均值随机变量。ζ的样本方差由下式给出：

$\mathrm{sample\; var}\left(\mathrm{\ζ}\right)=\frac{1}{p-1}\·\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{\left(z_i\right)}^2=\frac{1}{p-1}\·\underset{i=1}{\overset{p}{\mathrm{\Σ}}}({\tilde{a}}^T\·{\tilde{y}}_i)\·{({\tilde{a}}^T\·{\tilde{y}}_i)}^T=\frac{1}{p-1}\·\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{\tilde{a}}^T\·{\tilde{y}}_i\·{\tilde{y}}_i^T\·\tilde{a}$

$=\frac{1}{p-1}\·{\tilde{a}}^T\·(\underset{i=1}{\overset{p}{\mathrm{\Σ}}}\·{\tilde{y}}_i\·{\tilde{y}}_i^T)\·\tilde{a}={\tilde{a}}^T\·\tilde{S}\·\tilde{a}$

现在将就确定系数向量

以最大化ζ的样本方差

进行研究。由于可以通过按比例将

放大常数倍而使

与所需要的一样大，所以可以只通过施加诸如

的归一化条件来进行优化。为了使满足

的

最大化，标准途径是使用拉格朗日乘子技术。即寻求使下式最大化：

${\tilde{a}}^T\·\tilde{S}\·\tilde{a}-\mathrm{\λ}\·({\tilde{a}}^T\·\tilde{a}-1)$

其中，λ为拉格朗日乘子。对

求微分得出：

$\tilde{S}\·\tilde{a}-\mathrm{\λ}\·\tilde{a}=\tilde{0}$

由此，λ为的特征值，并且

为对应的特征向量，且

的特征值也就是使满足

的

最大化，

为

的具有最高特征值的特征向量。

令λ₁为对应特征向量为

的最大特征值，则

被称为的第一样本主分量。其为具有由

定义的观察值的随机变量。z_i1也被称为第一主分量上的第i个观察值的得分值。类似地，第二主分量

在满足与ζ₁不相关和归一化条件

的情况下，使得ζ₂的样本方差

最大化。对应于ζ₂的观察值为

ζ₁和ζ₂的样本协方差为：

$\mathrm{sample\; cov}({\mathrm{\ζ}}_1,{\mathrm{\ζ}}_2)=\frac{1}{p-1}\·\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{z}_{i1}\·z_{i2}=\frac{1}{p-1}\·\underset{i=1}{\overset{p}{\mathrm{\Σ}}}({\tilde{a}}_1^T\·{\tilde{y}}_i)\·{({\tilde{a}}_2^T\·{\tilde{y}}_i)}^T=\frac{1}{p-1}\·\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{\tilde{a}}_1^T\·{\tilde{y}}_i\·{\tilde{y}}_i^T\·{\tilde{a}}_2$

$=\frac{1}{p-1}\·{\tilde{a}}_1^T\·(\underset{i=1}{\overset{p}{\mathrm{\Σ}}}\·{\tilde{y}}_i\·{\tilde{y}}_i^T)\·{\tilde{a}}_2={\tilde{a}}_1^T\·\tilde{S}\·{\tilde{a}}_2$

$=\frac{1}{p-1}\·\underset{i=1}{\overset{p}{\mathrm{\Σ}}}({\tilde{a}}_2^T\·{\tilde{y}}_i)\·{({\tilde{a}}_1^T\·{\tilde{y}}_i)}^T={\tilde{a}}_2^T\·\tilde{S}\·{\tilde{a}}_1$

ζ₂与ζ₁不相关的要求等同于要求：

${\tilde{a}}_2^T\·\tilde{S}\·{\tilde{a}}_1={\tilde{a}}_2^T\·{\mathrm{\λ}}_1\·{\tilde{a}}_1={\mathrm{\λ}}_1\·[{\tilde{a}}_2^T\·{\tilde{a}}_1]=0\⇒{\tilde{a}}_2^T\·\tilde{S}\·{\tilde{a}}_1=0$ 且 ${\tilde{a}}_2^T\·{\tilde{a}}_1=0$

${\tilde{a}}_1^T\·\tilde{S}\·{\tilde{a}}_2={(\tilde{S}\·{\tilde{a}}_1)}^T\·{\tilde{a}}_2={\mathrm{\λ}}_1\·[{\tilde{a}}_1^T\·{\tilde{a}}_2]=0\⇒{\tilde{a}}_1^T\·\tilde{S}\·{\tilde{a}}_2=0$ 且 ${\tilde{a}}_1^T\·{\tilde{a}}_2=0$

这样，下列任何条件可以被用以规定ζ₂与ζ₁不相关：

${\tilde{a}}_2^T\·\tilde{S}\·{\tilde{a}}_1=0$ ${\tilde{a}}_2^T\·{\tilde{a}}_1=0$

${\tilde{a}}_1^T\·\tilde{S}\·{\tilde{a}}_2=0$ ${\tilde{a}}_1^T\·{\tilde{a}}_2=0.$

使满足ζ₂与ζ₁不相关且

的条件的ζ₂的样本方差

最大化等同于最大化：

其中，条件

被用以规定ζ₂与ζ₁不相关并且λ₂和为拉格朗日乘子。对

求微分，得到：

通过乘以

得出：

由于前面两项为零且所以其简化为

因此，

$\tilde{S}\·{\tilde{a}}_2-{\mathrm{\λ}}_2\·{\tilde{a}}_2=0$

所以λ₂再一次是

的特征值且

是对应的特征向量。

再者，

或者ζ₂的样本方差为λ₂。假定

没有重复的特征值，则λ₂不能等于λ₁。如果等于，则会有

从而违背了的约束条件。因而，λ₂为第二大特征值，且

为对应的特征向量。

可以看出，对第三个、第四个、……、第n个样本主分量使用相同技术，系数向量

为对应于第三大、第四大、……、以及最小特征值λ₃，λ₄，…，λ_n的

的特征向量。此外，主分量的样本方差也由

的特征值给出。

将

定义为主分量上的观察值的得分值的p×n矩阵。在此情况下，第i行由所有主分量ζ₁，ζ₂，…，ζ_n-1，ζ_n上的

的第i个观察值(即

)的得分值构成，即：

$={(\left[\begin{array}{c}{\tilde{a}}_1^T\\ {\tilde{a}}_j^T\\ {\tilde{a}}_n^T\end{array}\right]\·\left[\begin{array}{ccccc}{\tilde{y}}_1& \·\·\·& {\tilde{y}}_i& \·\·\·& {\tilde{y}}_p\end{array}\right])}^T={\left[\begin{array}{ccccc}{\tilde{y}}_1^T& \·\·\·& {\tilde{y}}_i^T& \·\·\·& {\tilde{y}}_{p_p}^T\end{array}\right]}^T\·\left[\begin{array}{ccccc}{\tilde{a}}_1& \·\·\·& {\tilde{a}}_j& \·\·\·& {\tilde{a}}_n\end{array}\right]$

$=\tilde{Y}\·\tilde{A}---\left(4a\right)$

其中，

为n×n矩阵，其各列为

的特征向量并且正交。

最大化后的第一主分量ζ₁由下式表示：

${\tilde{a}}_1^T\·\tilde{S}\·{\tilde{a}}_1=\frac{1}{p-1}\·{\tilde{a}}_1^T\·(\underset{i=1}{\overset{p}{\mathrm{\Σ}}}\·{\tilde{y}}_i\·{\tilde{y}}_i^T)\·{\tilde{a}}_1=\frac{1}{p-1}\·{\tilde{a}}_1^T\·(\underset{i=1}{\overset{p}{\mathrm{\Σ}}}({\tilde{x}}_i-\stackrel{\‾}{x})\·{({\tilde{x}}_i-\stackrel{\‾}{x})}^T)\·{\tilde{a}}_1$

为来自第i个像素(像素块)的时间-增强曲线(TEC)与来自全部像素(像素块)的TEC的均值的偏差。由此，特征向量：表示来自全部像素的TEC组中的第一、第二、第三、第四、……、以及最小主导地位的时间-增强性能。来自第i个像素(像素块)的TEC中的特征向量

的载荷(loading)为：

${\tilde{a}}_1^T\·{\tilde{y}}_i,{\tilde{a}}_2^T\·{\tilde{y}}_i,{\tilde{a}}_3^T\·{{\tilde{y}}_i\widetilde{,a}}_4^T\·{\tilde{y}}_i,\·\·\·,{\tilde{a}}_n^T\·{\tilde{y}}_i$

且为矩阵

的第i行，而来自全部像素(像素块)的TEC中的特征向量

的载荷为：

${\tilde{a}}_i^T\·{\tilde{y}}_1,{\tilde{a}}_i^T\·{\tilde{y}}_2,{\tilde{a}}_i^T\·{\tilde{y}}_3,{\tilde{a}}_i^T\·{\tilde{y}}_4,\·\·\·,{\tilde{a}}_i^T\·{\tilde{y}}_p$

且为矩阵

的第i列。因此，

的第i列给出了对应于第i大特征值λ_i的特征向量的载荷图(loading map)。

为了对来自CT灌注研究的动态图像去噪，需要对特征向量的载荷图的计算求逆。需要由截短后的一系列特征向量

重构原始动态图像的平滑版本

如等式(1)中所示，p×n矩阵

可以通过奇异值分解而展开为：

$\tilde{Y}=\tilde{U}\·\tilde{L}\·{\tilde{V}}^T$

其中，

为p×p正交矩阵，

为n×n正交矩阵，且

为p×n对角矩阵并将

的奇异值作为其对角线元素。令l₁，…，l_n为

的奇异值(在不失一般性的情况下假定

具有满列秩n，且p＞n)。使用奇异值分解：

$\tilde{S}=\frac{1}{p-1}\·{\tilde{Y}}^T\·\tilde{Y}=\frac{1}{p-1}\·{(\tilde{U}\·\tilde{L}\·{\tilde{V}}^T)}^T\·\tilde{U}\·\tilde{L}\·{\tilde{V}}^T=\frac{1}{p-1}\·\tilde{V}\·{\tilde{L}}^T\·{\tilde{U}}^T\·\tilde{U}\·\tilde{L}\·{\tilde{V}}^T$

$=\frac{1}{p-1}\·\tilde{V}\·{\tilde{L}}^T\·\tilde{L}\·{\tilde{V}}^T$

将等式两边乘以得到：

$\tilde{S}\·\tilde{V}=\frac{1}{p-1}\·\tilde{V}\·{\tilde{L}}^T\·\tilde{L}---\left(5\right)$

其中，

为p×1向量，其元素除等于l_i的第i个元素以外均为零，即为

且

${{\tilde{l}}_i}^T\·{\tilde{l}}_j=\left\{\begin{array}{cc}{l}_i^2& i=j\\ 0& i\≠j\end{array}\right.$

则，

${\tilde{L}}^T\·\tilde{L}=\left[\begin{array}{c}{{\tilde{l}}_1}^T\\ \·\\ \·\\ \·\\ {{\tilde{l}}_i}^T\\ \·\\ \·\\ \·\\ {{\tilde{l}}_n}^T\end{array}\right]\·\left[\begin{array}{ccccc}{\tilde{l}}_1& \·\·\·& {\tilde{l}}_i& \·\·\·& {\tilde{l}}_n\end{array}\right]=\left[\begin{array}{c}{\tilde{m}}_1^T\\ \·\\ \·\\ \·\\ {\tilde{m}}_i^T\\ \·\\ \·\\ \·\\ {\tilde{m}}_n^T\end{array}\right]$

其中 ${\tilde{m}}_i^T=\left[\begin{array}{ccccc}0& \·\·\·& l_i^2& \·\·\·& 0\end{array}\right].$

令

其中为n×1向量。

则，

$\tilde{V}\·{\tilde{L}}^T\·\tilde{L}=\left[\begin{array}{ccccc}{\tilde{v}}_1& \·\·\·& {\tilde{v}}_i& \·\·\·& {\tilde{v}}_n\end{array}\right]\·\left[\begin{array}{c}{\tilde{m}}_1^T\\ \·\\ \·\\ \·\\ {\tilde{m}}_i^T\\ \·\\ \·\\ \·\\ {\tilde{m}}_n^T\end{array}\right]=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\tilde{v}}_i\·{\tilde{m}}_i^T=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\tilde{v}}_i\·\left[\begin{array}{ccccc}0& \·\·\·& l_i^2& \·\·\·& 0\end{array}\right]$

$=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left[\begin{array}{ccccc}0& \·\·\·& l_i^2\·{\tilde{v}}_i& \·\·\·& 0\end{array}\right]=\left[\begin{array}{ccccc}{l}_1^2\·{\tilde{v}}_1& \·\·\·& l_i^2\·{\tilde{v}}_i& \·\·\·& l_n^2\·{\tilde{v}}_n\end{array}\right]$

代入等式(5)中得到：

$\tilde{S}\·\tilde{V}=\frac{1}{p-1}\·\tilde{V}\·{\tilde{L}}^T\·\tilde{L}=\left[\begin{array}{ccccc}\frac{l_1^2\·{\tilde{v}}_1}{p-1}& \·\·\·& \frac{l_i^2\·{\tilde{v}}_i}{p-1}& \·\·\·& \frac{l_n^2\·{\tilde{v}}_n}{p-1}\end{array}\right]$

$\⇒\left[\begin{array}{ccccc}\tilde{S}\·{\tilde{v}}_1& \·\·\·& \tilde{S}\·{\tilde{v}}_i& \·\·\·& \tilde{S}\·{\tilde{v}}_n\end{array}\right]=\left[\begin{array}{ccccc}\frac{l_1^2\·{\tilde{v}}_1}{p-1}& \·\·\·& \frac{l_i^2\·{\tilde{v}}_i}{p-1}& \·\·\·& \frac{l_n^2\·{\tilde{v}}_n}{p-1}\end{array}\right]$

$\⇒\tilde{S}\·{\tilde{v}}_i=\frac{l_i^2}{p-1}\·{\tilde{v}}_i,i=1,\·\·\·,n$

根据对特征向量和特征值的讨论有：

${\tilde{v}}_i={\tilde{a}}_i\⇒\tilde{V}=\tilde{A}$

$\frac{l_i^2}{p-1}={\mathrm{\λ}}_i$

或者

的奇异值等于

的特征值乘以(p-1)的平方根，其中

$\tilde{S}=\frac{1}{p-1}\·{\tilde{Y}}^T\·\tilde{Y},$ 即 $l_i=\sqrt{(p-1)\·{\mathrm{\λ}}_i}$

由此， ${\tilde{Y}}^T\·\tilde{Y}\·{\tilde{v}}_i=l_i^2\·{\tilde{v}}_i$ 或 ${\tilde{Y}}^T\·\tilde{Y}\·{\tilde{a}}_i=l_i^2\·{\tilde{a}}_i,i=1,\·\·\·,n---\left(5a\right)$

接下来，考虑

${\tilde{S}}^{*}=\frac{1}{p-1}\·\tilde{Y}\·{\tilde{Y}}^T=\frac{1}{p-1}\·\tilde{U}\·\tilde{L}\·{\tilde{V}}^T\·{(\tilde{U}\·\tilde{L}\·{\tilde{V}}^T)}^T=\frac{1}{p-1}\·\tilde{U}\·\tilde{L}\·{\tilde{V}}^T\·\tilde{V}\·{\tilde{L}}^T\·{\tilde{U}}^T$

$=\frac{1}{p-1}\·\tilde{U}\·L\·{\tilde{L}}^T\·{\tilde{U}}^T---\left(6\right)$

(注意，

不是

的转置矩阵)

将等式两边乘以

得到：

${\tilde{S}}^{*}\·\tilde{U}=\frac{1}{p-1}\·\tilde{U}\·\tilde{L}\·{\tilde{L}}^T$

$\tilde{L}\·{\tilde{L}}^T=\left[\begin{array}{ccccc}{\tilde{l}}_1& \·\·\·& {\tilde{l}}_i& \·\·\·& {\tilde{l}}_n\end{array}\right]\·\left[\begin{array}{c}{{\tilde{l}}_1}^T\\ \·\\ \·\\ \·\\ {{\tilde{l}}_i}^T\\ \·\\ \·\\ \·\\ {\tilde{l}}_n\end{array}\right]=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\tilde{l}}_i\·{{\tilde{l}}_i}^T$

因而，

令

为如所定义的p×1向量，则

$\tilde{L}\·{\tilde{L}}^T=\left[\begin{array}{c}{\tilde{m}}_1^{*T}\\ \·\\ \·\\ \·\\ {\tilde{m}}_i^{*T}\\ \·\\ \·\\ \·\\ {\tilde{m}}_n^{*T}\\ 0_{n+1}\\ \·\\ \·\\ \·\\ 0_p\end{array}\right]$

令

其中

为p×1向量。

则，

$\tilde{U}\·\tilde{L}\·{\tilde{L}}^T=\left[\begin{array}{cccccccc}{\tilde{u}}_1& \·\·\·& {\tilde{u}}_i& \·\·\·& {\tilde{u}}_n& {\tilde{u}}_{n+1}& \·\·\·& {\tilde{u}}_p\end{array}\right]\·\left[\begin{array}{c}{\tilde{m}}_1^{*T}\\ \·\\ \·\\ \·\\ {\tilde{m}}_i^{*T}\\ \·\\ \·\\ \·\\ {\tilde{m}}_n^{*T}\\ 0_{n+1}\\ \·\\ \·\\ \·\\ 0_p\end{array}\right]=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\tilde{u}}_i\·{\tilde{m}}_i^{T}l$

$=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\tilde{u}}_i\·\left[\begin{array}{ccccc}0& \·\·\·& l_i^2& \·\·\·& 0_P\end{array}\right]+\underset{i=n+1}{\overset{p}{\mathrm{\Σ}}}{\tilde{u}}_i\·\left[\begin{array}{ccccc}0& \·\·\·& 0_i& \·\·\·& 0_p\end{array}\right]$

$=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left[\begin{array}{ccccc}0& \·\·\·& l_i^2\·{\tilde{u}}_i& \·\·\·& 0_p\end{array}\right]=\left[\begin{array}{cccccccc}{l}_1^2\·{\tilde{u}}_1& \·\·\·& l_i^2\·{\tilde{u}}_i& \·\·\·& l_n^2\·{\tilde{u}}_n& 0_{n+1}& \·\·\·& 0_p\end{array}\right]$

代入等式(6)中得到：

${\tilde{S}}^{*}\·\tilde{U}=\frac{1}{p-1}\·\tilde{U}\·\tilde{L}\·{\tilde{L}}^T=\left[\begin{array}{cccccccc}\frac{l_1^2\·{\tilde{u}}_1}{p-1}& \·\·\·& \frac{l_i^2\·{\tilde{u}}_i}{p-1}& \·\·\·& \frac{l_n^2\·{\tilde{u}}_n}{p-1}& 0_{n+1}& \·\·\·& 0_p\end{array}\right]$

$\⇒\begin{array}{c}\left[\begin{array}{cccccccc}{\tilde{S}}^{*}\·{\tilde{u}}_1& \·\·\·& {\tilde{S}}^{*}\·{\tilde{u}}_i& \·\·\·& {\tilde{S}}^{*}\·{\tilde{u}}_n& 0_{n+1}& \·\·\·& 0_p\end{array}\right]=\left[\begin{array}{cccccc}\frac{l_1^2\·{\tilde{u}}_1}{p-1}& \·\·\·& \frac{l_i^2\·{\tilde{u}}_i}{p-1}& \·\·\·& \frac{l_n^2\·{\tilde{u}}_n}{p-1}& 0_{n+1}\end{array}\right.\end{array}$

$\⇒{\tilde{S}}^{*}\·{\tilde{u}}_i=\left\{\begin{array}{cc}\frac{l_i^2}{p-1}\·{\tilde{u}}_i& i\≤n\\ 0\·{\tilde{u}}_i& n+1\≤i\≤p\end{array}\right.$

由此，

为对应于第i个特征值

的的特征向量，

也碰巧为

的第i个特征值，同时全部

为具有零特征值的

的特征向量。并且

$\tilde{Y}\·{\tilde{Y}}^T\·{\tilde{u}}_i=l_i^2\·{\tilde{u}}_i,i=1,\·\·\·,n---\left(6a\right)$

如上所述，等式(5a)为：

${\tilde{Y}}^T\·\tilde{Y}\·{\tilde{v}}_i=l_i^2\·{\tilde{v}}_i$

将两边乘以

得到：

$\tilde{Y}\·{\tilde{Y}}^T\·\tilde{Y}\·{\tilde{v}}_i=l_i^2\·\tilde{Y}\·{\tilde{v}}_i$

与等式(6a)对比：

其中c为常数        (7a)

类似地，将等式(6a)两边乘以得到；

${\tilde{Y}}^T\·\tilde{Y}\·{\tilde{Y}}^T\·{\tilde{u}}_i=l_i^2\·{\tilde{Y}}^T\·{\tilde{u}}_i$

与等式(5a)对比：

其中c为常数        (7b)

由等式(7b)：

${\tilde{v}}_i=\frac{1}{c^{*}}\·{\tilde{Y}}^T\·{\tilde{u}}_i$

代入等式(7a)中得到：

$\tilde{Y}\·{\tilde{Y}}^T\·{\tilde{u}}_i=c^{*}\·c\·{\tilde{u}}_i$

代入等式(6a)中得到：

$l_i^2\·{\tilde{u}}_i=c^{*}\·c\·{\tilde{u}}_i$

因此， $c^{*}\·c=l_i^2$

在不失一般性的情况下，将c^*设为：

$c^{*}=c=l_i=\sqrt{(p-1)\·{\mathrm{\λ}}_i}i=1,\·\·\·,n$ 且n＜p

由此，对于

的奇异值分解：

$\tilde{Y}=\tilde{U}\·\tilde{L}\·{\tilde{V}}^T$

其中

为p×p正交矩阵：

且为p×1向量；

为n×n正交矩阵：

且

为n×1向量；

且

为p×n对角矩阵，且将

的奇异值l₁，…，l_n作为其对角线元素；在不失一般性的情况下还假定

具有满列秩n，且p＞n。

示出有：

其中

为对应于特征值

的的第i个特征向量；并且由等式(7c)有

${\tilde{u}}_i=\frac{1}{c}\tilde{Y}\·{\tilde{v}}_i=\frac{1}{\sqrt{(p-1)\·{\mathrm{\λ}}_i}}\·\tilde{Y}\·{\tilde{a}}_i$

如上所述，为了对除噪，可以在第r个奇异值后面截短SVD展开式：

$\tilde{Y}\≈\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{l}_i\·{\tilde{u}}_i\·{\tilde{v}}_i^T=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{l}_i\·{\tilde{u}}_i\·{\tilde{a}}_i^T=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{l}_i\·\frac{1}{\sqrt{(p-1)\·{\mathrm{\λ}}_i}}\·\tilde{Y}\·{\tilde{a}}_i\·{\tilde{a}}_i^T=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}\tilde{Y}\·{\tilde{a}}_i\·{\tilde{a}}_i^T---\left(8\right)$

$\⇒\tilde{Y}\≈\tilde{Y}\·\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{\tilde{a}}_i\·{\tilde{a}}_i^T$

等式(8)也提供了如下列所解释的非常简单的几何解释。

$\tilde{Y}\≈\tilde{Y}\·\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{\tilde{a}}_i\·{\tilde{a}}_i^T=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}\tilde{Y}\·{\tilde{a}}_i\·{\tilde{a}}_i^T=\left(\begin{array}{ccccc}\tilde{Y}\·{\tilde{a}}_1& \tilde{Y}\·{\tilde{a}}_2& \·\·\·& \tilde{Y}\·{\tilde{a}}_{r-1}& \tilde{Y}\·{\tilde{a}}_r\end{array}\right)\·\left(\begin{array}{c}{\tilde{a}}_1^T\\ {\tilde{a}}_2^T\\ \·\\ \·\\ \·\\ {\tilde{a}}_{r-1}^T\\ {\tilde{a}}_r^T\end{array}\right)=\tilde{Y}\·{\tilde{A}}_r\·{\tilde{A}}_r^T---\left(9\right)$

其中，

为n×r矩阵，其第i列为

且由等式(4a)得出：

${\tilde{Z}}_r=\tilde{Y}\·{\tilde{A}}_r$

并且

的行和列具有与

的行和列相同的解释。

改写等式(9)：

$\tilde{Y}\≈\tilde{Y}\·\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{\tilde{a}}_i\·{\tilde{a}}_i^T=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}\tilde{Y}\·{\tilde{a}}_i\·{\tilde{a}}_i^T=\left(\begin{array}{ccccc}\tilde{Y}\·{\tilde{a}}_1& \tilde{Y}\·{\tilde{a}}_2& \·\·\·& \tilde{Y}\·{\tilde{a}}_{r-1}& \tilde{Y}\·{\tilde{a}}_r\end{array}\right)\·\left(\begin{array}{c}{\tilde{a}}_1^T\\ {\tilde{a}}_2^T\\ \·\\ \·\\ \·\\ {\tilde{a}}_{r-1}^T\\ {\tilde{a}}_r^T\end{array}\right)$

等式(10)表明，第i个像素(像素块)的TEC为r个特征向量的加权和，且该加权的权重由第i个像素TEC在特征向量上的载荷来给定。

特征向量

可以被认为在n维空间中形成正交基。在

上的投影为

且

可以重组为

用以对

去噪的算法如下：

1.计算对应于

的p×n零均值矩阵，

其中

${\tilde{y}}_i={\tilde{x}}_i-\stackrel{\‾}{x}$ 且 $\stackrel{\‾}{x}=\frac{1}{p}\·\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{\tilde{x}}_i$

2.计算n×n样本协方差矩阵

$\tilde{S}=\frac{1}{p-1}\·{\tilde{Y}}^T\·\tilde{Y}$

3.计算n对特征值和特征向量，

$({\mathrm{\λ}}_i,{\tilde{a}}_i)i=1,\·\·\·,n$

4.通过只保留最大的r个特征值来计算n×n平滑矩阵：

$\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{\tilde{a}}_i\·{\tilde{a}}_i^T$

5.通过计算以下乘积来平滑

${\tilde{Y}}^{*}=\left[\begin{array}{ccccc}{\tilde{y}}_1^{*}& \·\·\·& {\tilde{y}}_i^{*}& \·\·\·& {\tilde{y}}_p^{*}\end{array}\right]\≈\tilde{Y}\·\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{\tilde{a}}_i\·{\tilde{a}}_i^T$

6.如下计算

的平滑版本：

${\tilde{X}}^{*}=\left[\begin{array}{ccccc}{\tilde{y}}_1^{*}+\stackrel{\‾}{x}& \·\·\·& {\tilde{y}}_i^{*}+\stackrel{\‾}{x}& \·\·\·& {\tilde{y}}_p^{*}+\stackrel{\‾}{x}\end{array}\right]$

现在将就实验结果来说明本申请中所描述的方法的示例。针对实现减少CT扫描剂量的能力来对本申请中所描述的主分量分析(PCA)进行测试。

现在将基于CT灌注研究来概述PCA方法，所述CT灌注研究产生捕获通过大脑的造影剂的通路的多个时间分离的图像。对应于该图像序列中的各个像素，生成将该像素中的造影剂浓度表示为时间的函数的曲线。还生成等于该曲线与全部像素曲线的均值之差的差曲线(difference curve)。计算所有差曲线对的协变的协方差矩阵并且其表示与平均曲线有关的全部可能偏差。主分量分析法分析协方差矩阵，以通过找出表示关于平均曲线的大部分偏差的差曲线的线性组合，而找出这些曲线之中的共有时域图案(pattern)。可以看出，这样的线性组合或者主分量，是协方差矩阵的特征向量，并且具有最大特征值的特征向量是全部差曲线之中占最大主导地位的时域特征，具有第二大特征值的特征向量为占第二大主导地位的时域特征，依此类推，则具有小特征值的特征向量来自于CT图像中的噪声。推论是，各个差曲线可以表示为全部主分量之和，并且通过丢弃具有小特征值的分量能够有效抑制噪声。为了对图像的时间序列去噪，将平均曲线加至通过使用具有大特征值的前几个主分量而重组的全部差曲线。

尽管目前可用的统计滤波器在图像重建之前消除了投影数据中的噪声，但主分量分析法消除了来自重建后的图像的噪声。下面描述的实验结果显示出执行PCA致使减少CT成像的剂量。

就CT灌注应用的辐射剂量减少，相对于传统的滤波反投影(FBP)法来测试PCA法。在健康的猪身上对PCA对比FBP在使用CT灌注来测量脑灌注上的剂量节约进行了测试，其中，对该健康的猪扫描四次，该四次扫描分别用从190mA改变到100mA、75mA和50mA的x射线管电流而保持其他扫描参数相同。

图3示出了用下述参数和方法获得的CT灌注图像且利用J-W模型而计算的猪的脑血流量图：(A)190mA且FBP重建；(B)100mA，FBP重建且PCA；(C)70mA，FBP重建且PCA；以及(D)50mA，FBP重建且PCA。图3(E)示出了用于计算图4中所示血流量图的品质因数的受关注区域(3个勾勒出的圈)。这些区域重叠在示出了通过猪的大脑的冠状断面的CT图像上。暗区域为颅骨和颌骨。

图4示出了图3(E)中所示三个区域中的血流量的标准差/均值(品质因数)：70mA且PCA比190mA且FBP好，而50mA PCA比190mA且FBP差(p＜0.05)。图3中血流量图的图像质量也反映了此发展(progression)。因此，PCA法能够减少用以在CT灌注研究中获得图像并且相比于仅使用FBP而言维持所得到的血流量图的质量所需要的辐射剂量。更具体来说，与仅使用FBP相比，PCA法使辐射剂量减少至将近三分之一。因为血脑屏障保持完整无缺，所以没有展示出PCA法在生成BBB-PS图时对剂量减少的效果。然而，由于血流量、血容量和BBB-PS是在J-W模型中一起计算的，所以预期利用血流量图示出的剂量减少适用于BBB-PS图。

本申请中所描述的主分量分析法已利用C++而实现于计算机可执行软件应用程序中。在诸如例如个人计算机(PC)的通用计算设备上用以处理九十五(95)个512×512的CT脑图像所需的时间少于两(2)分钟。

本申请中所描述的方法可以在诸如图5中所示的系统400的硬件上加以实现。输入510接收多个时间分离的图像，该多个时间分离的图像可以是之前存储的，从成像设备直接接收的，或者从远程位置经由有线或无线通信链接传送的。通用计算设备505接收来自输入的时间分离的图像并执行软件应用程序，从而如上所述处理时间分离的图像。通用计算设备505通过显示器520、键盘515和/或鼠标或者其他指向设备525而与用户进行接口交互。结果，例如所构建的处理后的图像，可以呈现在显示器510和/或输出至任何适当的输出设备530，例如，打印机、存储设备或者用于将结果传送给远程位置的通信设备。

本申请中所描述的软件应用程序可以作为单机应用程序运行，或者可以包含在其他可用应用程序中，以向那些应用软件提供增强的功能。所述软件应用程序可以包括包含例程、程序、对象构件、数据结构等的程序模块，并且可以实现为存储在计算机可读介质上的计算机可读程序代码。所述计算机可读介质是能够存储数据的任何数据存储设备，而所述数据之后可由计算机系统读取。计算机可读介质的示例包括，例如，只读存储器、随机存取存储器、CD-ROM、磁带和光学数据存储设备。所述计算机可读程序代码也可以分布在包含耦接的计算机系统的网络上，从而计算机可读程序代码可以分布式方式存储和执行。

上述实施例旨在作为示例，并且本领域技术人员可以在不脱离由所附的权利要求书所限定的本发明的范围的情况下对其进行改变和修改。

Claims (28)

1.一种处理多个时间分离的图像的方法，包括：

在各个图像中选择多个成像单元；

为各个成像单元确定时域差；以及

选择在阈限以上的时域差。

2.如权利要求1所述的方法，进一步包括：在所述选择之前对各个测得的时域差排序。

3.如权利要求1或2所述的方法，进一步包括：使用所选择的时域差来构建处理后的图像。

4.如权利要求1至3中任一项所述的方法，其中，各个成像单元为一个像素和一个体素中的一种。

5.如权利要求1至3中任一项所述的方法，其中，各个成像单元为多个像素和多个体素中的一种。

6.如权利要求1所述的方法，其中，所述时域差基于对比增强。

7.如权利要求6所述的方法，其中，所述时域差为造影剂的差别。

8.如权利要求7所述的方法，其中，对比增强基于所述造影剂的浓度。

9.如权利要求7所述的方法，其中，对比增强基于所述造影剂的激活。

10.如权利要求1所述的方法，进一步包括：利用盲信号分离统计技术来表征所测得的时域差。

11.如权利要求1所述的方法，其中，确定所述时域差包括计算所述成像单元的时域差的协方差矩阵。

12.如权利要求11所述的方法，进一步包括：利用盲信号分离统计技术分析所述协方差矩阵，以确定所述时域差之中的共有时域图案。

13.如权利要求12所述的方法，进一步包括：分析所述协方差矩阵的特征向量。

14.如权利要求13所述的方法，其中，所述特征向量的特征值表示时域差的程度。

15.如权利要求14所述的方法，其中，所述选择时域差包括选择在所述阈限以上的特征值。

16.如权利要求1至15中任一项所述的方法，其中，所述阈限是能够调节的。

17.如权利要求12所述的方法，其中，所述盲信号分离统计技术选自由下述技术构成的组：主分量分析、独立分量分析、盲源解卷积、卡亨南-拉维变换和霍特林变换。

18.如权利要求1至17中任一项所述的方法，其中，所述多个时间分离的图像是使用动态医学成像技术而获得的。

19.如权利要求18所述的方法，其中，所述动态医学成像技术包括造影剂的使用。

20.如权利要求18所述的方法，其中，所述时间分离的图像是使用具有减少的辐射剂量的方案而获得的。

21.如权利要求18所述的方法，其中，所述时间分离的图像是使用具有减少的造影剂量的方案而获得的

22.如权利要求18所述的方法，其中，所述动态医学成像技术为磁共振成像(MRI)、功能性MRI、正电子发射断层摄影术、核医学或计算机断层摄影术。

23.如权利要求18所述的方法，其中，所述动态医学成像技术为计算机断层摄影术(CT)。

24.如权利要求23所述的方法，其中，所述时间分离的图像是使用具有减少的辐射剂量的CT灌注方案而获得的。

25.如权利要求23所述的方法，其中，所述时间分离的图像是使用具有减少的造影剂量的CT灌注方案而获得的。

26.一种处理多个时间分离的图像的系统，包括：

接口，用于接收多个时间分离的图像；和

处理结构，用于在各个图像中选择多个成像单元、为各个成像单元确定时域差以及选择在阈限以上的时域差。

27.一种包含处理多个时间分离的图像的计算机程序的计算机可读介质，所述计算机程序包括：

用于在各个图像中选择多个成像单元的计算机程序代码；

用于为各个成像单元确定时域差的计算机程序代码；和

用于选择在阈限以上的时域差的计算机程序代码。

28.一种用于根据权利要求1至25中任一项所述的方法来处理多个时间分离的图像的系统。

Images