一种探测器晶体位置表的建立方法及系统与流程

本发明涉及正电子发射断层(PET)成像领域，尤其是涉及一种正电子发射断层成像系统的探测器晶体位置表的建立方法及系统。

背景技术：

在PET系统中，最直接影响图像分辨率的部分是探测器部分。PET探测器多采用闪烁晶体耦合光电倍增管的结构组成。正电子湮灭事例经PET探测器像素阵列解码之后产生事例的位置分布，统计后形成反映晶体阵列分布的二维直方图(2D flood histogram image，以下简称FHI)，由于探测器实现过程中的编解码的非线性，使得生成的FHI不规则。高分辨率PET系统相比于传统的PET系统，探测器中的闪烁晶体数目增多，尺寸减小，密度增大，由此导致晶体阵列不规则的问题也更多更严重。以一个探测器单元为例，晶体阵列为16×16的等距规则排列像素阵列，解码之后会出现蝶形、旋转、压缩、扩张等整体阵列轮廓的形变、能量不均、阵列行列扭曲、晶体黏连、晶体缺失以及晶体伪影等问题。若不做任何处理，数据获取时，在线处理无法判断晶体像素的行列坐标，直接导致图像分辨率下降，严重的将会导致图像信息错误。

目前业界广泛使用的做法是使用晶体像素查找表(Crystal Lookup Table，CLT)，对获得的图像进行矫正和纠偏。由于漂移随时间推移会不停恶化，使CLT每隔一段时间就要被重新绘制。现阶段CLT靠手工绘制而得，对图像矫正效果尚可。但PET设备上往往有几万个晶体，手工绘制过程耗时耗力，一台机 器的CLT需要3个有经验的工程师工作一个上午，成为PET设备维护中的主要时间瓶颈之一。

设备维护所花费的时间，不仅影响到企业的人工成本，还影响到客户对设备的使用率。知名的设备供应商为提升售后服务的效率和客户满意度，积极提高CLT绘制的自动化程度。为了既保证矫正效果，又缩短维护时间，目前最常用的方法是半自动绘制CLT，即计算机自动生成一个中间结果，技师在此基础上检查和修正结果，最后再自动生成CLT。中间结果质量越好，人为检查和修正花费的时间越少。

已有的自动绘制CLT算法可以分为晶体中心检测法，晶体行列检测法，和晶体边界检测法。

晶体行列检测法指基于一维投影的方法，先把二维图像投影到两个方向(最常见的是将图像分别投影在X方向和Y方向分别投影)。如此一来，二维检测问题就被简化为一维寻峰问题，求解可以使用一维的极值检测，一维求导定位法，或者更复杂的Gaussian fitting。投影可以使用加法，也可以使用方差。UC Davis的Abhijit Chaudhari提出的方法本质是基于Fourier变换对X，Y方向做投影。基于投影的算法充分利用了晶体阵列(探测器)是方形这个先验条件，优点是基本不受噪声或图像灰度不均匀的影响，能简单地得到全局解，而且反投影之后的晶体阵列的排序坐标很容易得知，不需要后续复杂的晶体映射算法。缺点是过于依赖方形这个全局先验条件，而对局部信息和行列不规则形变考虑不足。首先，由于缺乏局部信息，其所得到的波峰位置缺乏精准性。此外，可利用初步定位波峰位置的方法，后续在每个初始位置的领域内用局部最大值法进行定位细化。尽管如此，当晶体行列倾斜(不完全垂直于X或Y轴)，或者行列弯曲变形较厉害时，投影图像上将无法形成易于区分的 高斯波形，则一维寻峰将出错甚至失败。

晶体中心检测法大致包括基于晶体亮斑特征的检测方法和聚类算法。由于亮斑在图像上有许多明显的视觉特征，可以直接用于晶体检测。一种方式是利用亮斑中心往往是局部最大值(local maxima)的性质在图像中进行探测。另一种方式是探测亮区域(像素值在10’000到100’000之间的区域)的质心点(centroids)作为亮斑中心点。更为严谨的思路是基于分水岭降水模拟算法的二维峰值检测。在图像上直接对每个亮斑进行检测的优点是简单直观，运算快捷，但其受噪声，伪影，图像灰度不均匀性等影响，容易出现漏检和错检。其实际性能往往取决于图像质量和预处理，稳定性是主要问题。而聚类算法则是基于每个亮斑实质是落到一个晶体中光子的分布，把每个电子的落点位置作为一个数据，那么波峰就是数据最密集的位置，即位置数据的聚类中心。这类算法中最常用的是二维高斯混叠模型(Gaussian Mixture Model，GMM)法。由于随机统计量理论上应该以晶体中心为圆心呈二维高斯分布，且不同晶体接收到的电子有重叠部分，因此使用GMM来定位高斯波峰也就比较合适。GMM不依赖晶体阵列是方形这个先验条件，从而对晶体阵列的形变，噪声，不均匀性等等问题均不敏感。但其依赖二维高斯分布这个假设，而现实中的晶体亮斑因形变和融合，形状不一定符合高斯分布。比GMM更泛化的聚类算法被提了出来，如FCM(Fuzzy C-means)和SOFM(Self-Organized Feature Map)等。这类算法使用了更宽松的空间距离约束来代替GMM使用的高斯密度约束，从而不再假设亮斑形状。其中SOFM以神经网络的形式完成的，可通过设置神经元之间的拓扑结构来隐性植入晶体数量及晶体间的空间关系，因此被广泛用于CLT绘制中。聚类算法本质上是一个优化(optimization)问题，其优点是良好的精确性和鲁棒性，而缺点在于其运算复杂度和对伪影的敏感度。一方面，其运算 时间与需要估计晶体数量呈正比，而现代PET海量晶体的事实使这类算法变得不太实用。另一方面，聚类算法往往更容易受伪影的误导。事实上，优化算法迭代的停止条件往往是实际晶体的数量。实际使用中由于伪影十分常见，晶体数量难以准确判定。而伪影形成的鬼点与真实波峰无明显差异，有些甚至比真实点亮度更高。基于投影的算法只有在伪影成行成列时才会出错，而聚类算法缺少空间约束，即使独立的伪影也能使算法结果出错。伪影造成的错探将极大地影响探测结果，并导致后续的晶体映射错位，严重的会导致整个晶体阵列映射错位。

也有算法尝试直接寻找亮斑间的波谷作为分界线，即晶体边界检测法，最常规的就是，反转图像灰度，使波谷成为山脊，然后直接使用分水岭算法。但反转图像中的山脊几乎是受噪声污染最严重的部分，极易导致分水岭算法的过分割。相比之下，邻域标准差NSD(neighborhood standard deviation)方法是更稳定的边界检测法。它使用局部统计矩在局部投影图上做检测，行列形变的影响在局部表现并不显著，而对行列之间的波谷响应很敏感。即便如此，在行列扭曲变形很严重或者晶体发生黏连的时候，NSD的分割效果依然不理想。

尽管已有数量众多的晶体检测算法，但每种算法都有其适用的局限性。在不规则问题众多的高分辨率PET探测器形成的FHI上，很难使用同一种算法解决所有的问题，获得理想的晶体分割准确率。

技术实现要素：

本发明所要解决的技术问题是提供一种探测器晶体位置表的建立方法及系统。

本发明为解决上述技术问题而采用的技术方案是：一种探测器晶体位置表的建立方法，包括以下步骤：根据正电子湮灭事例分布，统计后形成反映晶体阵 列分布的二维直方图；采用多种探测方法分别检测二维直方图中的晶体，并获得相应的数种判定结果；将相应的数种判定结果进行融合，并根据融合后的判定结果的置信度确定在所述融合结果中无需矫正的部分及需要矫正的部分；对融合结果中需要矫正的部分进行矫正，获得矫正后的判定结果；使用矫正后的判定结果对二维直方图中的晶体进行分割，获得晶体位置表。

进一步的，所述多种探测方法包括探测方法一及探测方法二，采用探测方法一检测二维直方图中的晶体，获得第一判定结果，包括晶体的位置信息、行信息及列信息；采用探测方法二，检测二维直方图中的晶体，获得第二判定结果，包括晶体的位置信息、行信息及列信息。

进一步的，所述探测方法一与探测方法二为相同类型或不同类型的探测方法。

进一步的，所述探测方法一为晶体中心探测法，所述探测方法二为晶体边界探测法。

进一步的，所述融合结果包括针对每个探测到的晶体位置计算一个置信度。

进一步的，所述置信度从高到低分级为一级融合、二级融合、三级融合；一级融合是指多个判定结果在相同晶体位置的行列信息都相同，此位置的融合结果无需矫正；二级融合是指多个判定结果在相同位置都探测到了晶体，但其行列信息不相同或部分相同，此位置的融合结果需要矫正不一致的行列信息；三级融合是指多个判定结果中仅少数在此位置探测到了晶体，在融合结果中此位置及其行列信息将被暂时弃用。

进一步的，所述矫正不一致的行列信息是指根据高置信度的晶体位置及其行列信息对低置信度的晶体位置进行行列拟合，从而得到低置信度晶体位置的新的行列信息。

进一步的，所述行列拟合是指基于同行或列的多个高置信度的晶体位置，进行样条拟合，并使其穿过低置信度的晶体位置，从而矫正这些位置的行列信息。

本发明为解决上述技术问题而采用的技术方案是：一种探测器晶体位置表的建立系统，包括：探测器单元，所述探测器单元包括多个闪烁晶体及光电倍增管，所述闪烁晶体用于获得光子信号，光电倍增管用于将光信号转换为电信号；采样单元，用于对所述电信号进行数字采样，获得反映光子事例的位置分布的数据；控制及数据处理单元，用于对光子事例的位置分布的数据统计，并形成反映晶体阵列分布的二维直方图；所述控制及数据处理单元对二维直方图中的晶体采用多种晶体探测方法检测，获得多种判定结果，并将判定结果进行融合，并根据融合结果的置信度对判定结果进行自动矫正；使用矫正后的判定结果对二维直方图中的晶体进行分割，获得晶体位置表。

进一步的，所述探测器晶体位置表的建立系统还包括显示单元，用于显示二维直方图以及探测到的晶体。

进一步的，所述探测器晶体位置表的建立系统还包括存储单元，所述存储单元用于保存融合结果，所述融合结果包括针对每个探测到的晶体位置计算一个置信度。

本发明对比现有技术有如下的有益效果：本发明通过将多种不同的探测方法有效地融合，取长补短，则可以自动矫正大部分的探测错误，提升晶体自动分割的准确率。特别是在设备老化，FHI降质严重时，对比单一探测方法的结果，基于多种探测方法融合的结果能大幅提升自动分割的准确率，从而极大地降低人工检查和修正花费的时间。

【附图说明】

图1为本发明实施例的晶体位置表建立方法示意图；

图2为本发明实施例的对探测结果融合处理示意图；

图3为本发明实施例的二级融合示例图；

图4为多种探测方法融合矫正前后示例图；

图5为本发明实施例探测器晶体位置表建立系统示意图。

【具体实施方式】

下面结合附图和实施例对本发明作进一步的描述。

如图1-5所示，本发明实施例中的PET系统的探测器由探测器单元1、采样单元2、存储单元3、控制及数据处理单元4及显示单元5组成。

探测器单元1包括400到500个晶体阵列，每个阵列中有16*16共256个LYSO晶体，晶体的物理尺寸为2.3*2.3mm²。LYSO晶体阵列耦合光导和位置敏感型光电倍增管。入射光子被探测器单元转换为闪烁脉冲信号，采样单元2中对位置脉冲信号进行数字采样，并将采样结果置于存储单元3，控制及数据处理单元4可以为CPU、DSP、单片机，其用于对存储单元3中的采样结果按照一定的算法或程序处理后产生光子事例的位置分布，进一步统计后形成反映晶体阵列分布的二维直方图(2D flood histogram image，以下简称FHI)。所述二维直方图可在显示单元5上被显示出来。FHI是每个位置接收到的光子数的统计图，图像上每个亮斑(晶体)对应探测器中的一个晶体。

以上是通过全自动绘制方法获得的二维直方图(FHI)或者初始晶体像素查找表(Crystal Lookup Table，CLT)，因为噪声、晶体黏连、伪影以及阵列变形的影响，CLT的自动绘制算法获取的结果，需要人工确认或校正。

为减少人工校正所花费的时间，有必要尽可能提高自动绘制结果的准确性 和鲁棒性，本发明的控制及数据处理单元4对二维直方图中的晶体采用多种探测方法检测，获得多种判定结果，并将判定结果进行融合，并根据融合结果的置信度确定需要对所述二维直方图中进行矫正的晶体；对矫正后的二维直方图中的晶体进行分割，获得晶体位置表。

已有的自动绘制CLT的算法由于其基于不同的原理，在实践中各有其优势和问题。例如通过检测晶体中心来初步定位晶体位置，并结合其空间信息确定每个位置实际对应的晶体的方法(以下简称“中心探测法”)，其优势是一般的晶体能被稳定地探测到，且中心位置准确，但也会因为局面的噪声或伪影而出现错误探测结果，且后续的晶体排序依然会出问题。而通过将FHI投影到一维找晶体行列位置的方法(以下简称“行列探测法”)，能很好的避免噪声或伪影的干扰，但对于阵列变形很敏感。再例如直接寻找晶体的边界(以下简称“边界探测”)，由于往往结合了拓扑信息和图像能量信息，对噪声或伪影均不敏感，在遇到晶体黏连边界不清就无法分辨。

本发明的一个重点是选择属于两种或两种以上不同类型的检测方法对二维直方图(FHI)中的晶体(亮斑、峰或散点)进行检测，并基于检测结果之间的相关程度，对检测结果进行自动矫正。但本发明完全可以拓展用于同一种类型的两种或两种以上检测方法的融合和自动矫正。

本发明一个实施方式中对二维直方图进行矫正的步骤如下：

先将湮灭光子源置于PET系统的检测腔内，湮灭光子源不断发射光子，探测器单元不断接收光子，经统计后最终得到反映晶体阵列分布的二维直方图或二维位置散点图，二维位置散点图中的每个像素记录的是每个坐标位置上检测到的入射γ光子事件数；

采用多种探测方法分别检测二维直方图中的晶体，并获得相应的数(多) 种判定结果；

将相应的数种判定结果进行融合，并根据融合后的判定结果的置信度确定在所述融合结果中无需矫正的部分及需要矫正的部分；

对融合结果中需要矫正的部分进行矫正，获得矫正后的判定结果；

使用矫正后的判定结果对二维直方图中的晶体进行分割，获得晶体位置表。

具体的，所述融合结果针对每个探测到的晶体位置计算一个根据置信度。所述置信度从高到低分级为一级融合、二级融合、三级融合。一级融合是指多个判定结果在相同晶体位置信息的行列信息都相同，此位置的融合结果无需矫正；二级融合是指多个判定结果在相同位置都探测到了晶体，但其行列信息不相同或部分相同，此位置的融合结果需要矫正不一致的行列信息；三级融合是指多个判定结果中仅少数在此位置探测到了晶体，在融合结果中此位置及其行列信息将被暂时弃用。

所述矫正不一致的行列信息是指根据高置信度(一级融合结果)的晶体位置及其行列信息对低置信度(二、三级融合结果)的晶体位置进行行列拟合，从而得到低置信度晶体位置的新的行列信息。本实施例中所述行列拟合是指基于同行或列的多个高置信度的晶体位置，进行样条拟合，并使其穿过低置信度的晶体位置，从而矫正这些位置的行列信息。但自动矫正也可以使用其它方法进行，比如使用带有先验信息的网格(Grid)模板，根据高置信度的晶体位置对网格进行形变，然后根据网格形变后固有的行列信息，在低置信度的晶体位置上矫正其行列信息。

对矫正后的二维直方图中的晶体进行分割，获得晶体位置表。具体的，根据所述晶体中心位置图划分各相邻晶体的分界线，生成晶体像素查找表。

以下以采用两种探测方法为例，对本发明做进一步的具体说明：

各种探测方法所对应的算法或程序预先配置到控制及数据处理单元中，并将为各探测方法打分，可以为每种探测方法设定一个对应的权重，以便于选择。边界探测法

选择两种(或多种)探测方法对其进行检测。原则上，优先选择两种不同类型的检测方法，以增加方法之间的互补性，例如中心探测法和边界探测法。若选择的是同一种类型的两种检测方法，则尽量不选在检测后可能出现相同缺陷的两种检测方法。所述的缺陷是指受噪声、几何形变、伪影等影响所产生的不准确结果，如漏检、错检等。

具体实施例中，采用晶体中心探测法作为探测方法一，检测二维直方图的晶体，获得第一判定结果，所述第一判定结果包括晶体块中心点的位置，中心点的行列信息；采用晶体边界探测法作为探测方法二，检测二维直方图的晶体，获得第二判定结果，所述第二判定结果包括被探测到的边界所包围的位置，及其行列信息。本实施例的中心检测法，可参考中国专利申请公开第CN104809460号，名称为晶体中心位置图生成方法及晶体像素查找表生成方法。本实施例中的边界检测法，可参相关论文“Wei，Q.，et al.A neighborhood standard deviation based algorithm for generating PET crystal position maps.in IEEE Nuclear Science Symposium and Medical Imaging Conference Record(NSS/MIC).2013。”但是需要特别说明的是，以上两种不同类型的检测方法，并不局限于这两篇现有技术文献。

将第一判定结果与第二判定结果进行融合，并根据融合结果的置信度确定是否需要对所述二维直方图进行矫正。对于中心探测法，参与融合的信息有中心点的位置、及其行列信息，对于边界探测法，参与融合的信息有边界所包围的位置，及其行列信息。

分级融合是本发明的核心步骤，分级融合的关键有：选取探测结果中的参与融合的信息、信息的相似性量化分级、对不同置信度信息的处理方式。按照置信度从高到低分级为一级融合、二级融合和三级融合。

一级融合表示所有的判定结果信息都相似，如融合同类型的方法的结果，对中心探测法与中心探测法的融合，探测方法一中的第N行M列的晶体中心与探测方法二中的第N行M列的晶体中心差不多处于同一位置；对边界探测法与边界探测法的融合，探测方法一中的第N行M列的晶体边界所包围的位置与探测方法二中的第N行M列的晶体边界所包围的位置大致位于同一位置或相互覆盖；如融合中心探测法与边界探测法这种不同类型方法的结果，中心探测法结果中的第N行M列的晶体中心正好处在边界探测法中的第N行M列的晶体边界所包围的位置内；上述情况都被视作一级融合。一级融合中的信息置信度最高，在融合结果中予以全部保留。

二级融合表示结果信息中部分相似。具体来讲，结果信息中部分相似首先是位置信息一致，如二种探测方法在同一位置附近都探测到了晶体中心，或一种探测方法找到的某个晶体边缘包围的区域包括另一种方法找到的晶体中心。但是，两种方法仅对该晶体的位置达成一致，对其行列信息判定结果不完全一样。“不完全一样”可以是行(列)信息不同但列(行)信息相同，行列信息皆不同但与周边晶体同属一行(列)，或行列信息皆不同且与周边晶体不属于同一行(列)。依据“不完全一样”的程度，来判断置信度的高低，自动矫正时由高到低依次对结果进行处理。自动矫正处理的依据是基于一级融合结果的样条拟合。下面结合图3来具体说明：

如图3圈中的三个晶体点，探测方法一与探测方法二中在相同位置都探测到了这些晶体点，对该晶体点的列坐标标识结果相同，虽然行坐标的标识结果 不同，但是分别都属于同一行，由此可以推断这三个晶体点很有可能是同一行的，如果出现这种情况的晶体点越多，那么这些晶体点是同一行的置信度就越高。基于一级融合的信息，利用样条(行列)拟合法(spline fitting)拟合行列。根据拟合行列的位置信息来判断这三个晶体点究竟属于哪行，并在融合结果中矫正这三个晶体点的行信息。

当待融合的方法对某个晶体的位置都无法达成一致，即某个位置只出现在一种方法的判定结果中，则有可能此位置不存在晶体。这种情况列为三级融合，对三级融合中的结果只保留位置信息，不保留行列信息。然后根据一级融合和二级融合矫正的结果，以及根据探测器晶体的数量及布置，判断在此保留位置是否真的存在晶体。如果真的存在晶体，则根据行列拟合结果更新其行列信息。最终得到矫正后的判定结果。所述矫正后的判定结果从实际情况来包含：无需矫正的部分晶体、及对需要矫正的部分的晶体进行矫正后的所得到的那部分。

本发明实施例中优选的探测方法有基于局部直方图均衡的SNMD方法(探测结果参图4(a))、基于高斯差分的SNMD方法、基于傅里叶模型的非刚性配准方法(探测结果参图4(b))以及NSD方法(探测结果参图4(c))，其中前两种方法为中心探测法，第三种为行列探测法，最后一种方法为边界探测法。从图4中可以看出，单种探测算法在圈标出的部分出错，但是经过本发明基于多探测方法融合的晶体探测自动矫正算法的融合处理后，最终的融合结果(图4(d))综合考虑了各探测结果的置信度，扬长避短，将所有的晶体都探测排序准确，符合预期结果。

本发明可用于任意多个探测算法的探测结果融合，不限探测法的原理，只需单个算法能独立给出探测结果。

由于不同类型的探测算法得到的判定结果不同，在融合时，可以将一种探测结果信息转换成另一种探测结果，以便于比较。例如，中心探测法的晶体中心位置信息和边界探测法的边界位置信息是可以互转的。使用边界探测法得到的大致的晶体所在范围后，在其子区域内可以通过寻找局部极大值等方法反推晶体的中心位置，以此和其他中心探测法中的晶体块中心点信息进行融合；或者也可以将中心探测法得到的中心位置信息，通过距离变换等方法转为边界信息，以便和边界探测法得到的边界信息比较融合。

在结果分级中，也可以根据信息量的多少，再进行量化和细化，归成不同于三层的融合结构。

本发明的优势在于：用参与融合的探测结果的相似性作为判断结果置信度权值的方式，最大程度地肯定了不同结果之间的“共识”，同时对差异信息设计了一种取舍准则，不盲目地信任和舍弃，从而得到更鲁棒的融合结果。

本发明通过将多种不同的探测方法有效地融合，取长补短，则可以自动矫正大部分的探测错误，提升晶体自动分割的准确率。特别是在设备老化，FHI降质严重时，对比单一探测方法的结果，基于多种探测方法融合的结果能大幅提升自动分割的准确率，从而极大地降低人工检查和修正花费的时间。

虽然本发明已以较佳实施例揭示如上，然其并非用以限定本发明，任何本领域技术人员，在不脱离本发明的精神和范围内，当可作些许的修改和完善，因此本发明的保护范围当以权利要求书所界定的为准。