CN102426696B

CN102426696B - 光学投影断层成像运动伪影校正方法

Info

Publication number: CN102426696B
Application number: CN 201110326234
Authority: CN
Inventors: 朱守平; 陈冬梅; 梁继民; 陈多芳; 屈晓超; 赵恒�; 李军; 陈雪利; 侯彦宾; 田捷
Original assignee: Xidian University
Current assignee: Xidian University
Priority date: 2011-10-24
Filing date: 2011-10-24
Publication date: 2013-04-17
Anticipated expiration: 2031-10-24
Also published as: CN102426696A

Abstract

本发明公开了光学投影断层成像运动伪影校正方法，使用不基于特征点的估计样品运动轨迹的方法，其实现步骤如下，(1)获取投影数据；(2)计算投影数据的零阶矩；(3)计算投影数据的一阶矩；(4)计算投影数据的质心；(5)获取运动参数；(6)计算运动量；(7)运动伪影校正。本发明中利用数据一致性条件建立样品与投影数据之间的对应关系，采用多项式估计被扫描样品的运动，进而直接从投影数据中估计出样品的运动信息。本发明可应用于样品的投影断层重建，可以提高光学断层成像***的空间分辨率，减少图像伪影。

Description

光学投影断层成像运动伪影校正方法

技术领域

本发明属于医学影像处理领域，更进一步涉及一种对光学投影断层成像***的运动伪影进行校正的方法。该方法可用于光学投影断层成像的图像处理。

背景技术

光学投影断层成像(以下简称OPT)是一种新型的光学分子影像成像技术，其成像原理和X射线计算机断层成像的原理类似。OPT可获得扫描样品结构像，又能利用荧光染料或荧光蛋白进行分子特异性标记，实现分子特性成像，并且设备成本低，使用方便。在OPT扫描过程中，被扫描样品的运动以及OPT***的机械不稳定性都会导致重建图像中产生运动伪影。运动伪影会导致***空间分辨率下降、产生图像伪影甚至使重建图像无法正常使用。

为了抑制运动伪影，需要在图像重建过程中进行校正。主要方法通过加权以减少一致性最差的投影数据对图像的贡献，消除或减弱运动带来的影响。

伊利克股份有限公司在其申请的专利“CT扫描中运动伪影的减少”(专利申请号200780053157.9，公开号CN101842807A)中公开了一种丢弃包含运动结构的投影图像以获得投影图像的一致子集的方法。该方法是通过自动检测进而标识出需要被丢弃的有角区域，该方法可以在处理偶发性运动的情况下进行图像伪影校正。该方法仍存在的不足是，这种方法只能抑制小幅度运动引起的伪影，由于丢弃了包含运动结构的投影图像，会降低图像信噪比。

U.J.Birk等人在“Correction for specimen movement and rotation errorsfor in-vivo Optical Projection Tomography，Biomedical Optics Express，vol.1，no.1，pp.87-96，2010.”中提出了通过准确估计出样品的运动轨迹，从而在图像重建过程中进行补充以实现运动伪影的校正。主要方法是使用基于特征点的运动估计方法。该方法需要样品的投影数据中具有明显易辨识的特征点，然后利用手工、自动或半自动的方法提取这些特征点，利用特征点的运动轨迹估计扫描样品的运动信息。该方法的不足是，需要估计特征点的运动轨迹，不适用于不易辨识的特征点。

发明内容

本发明的目的在于克服上述现有技术的不足，提出了一种光学投影断层成像运动伪影校正的方法。采用多项式函数对被扫描样品的运动进行建模，利用数据一致性条件建立样品与投影数据之间的对应关系，进而直接从投影数据中估计出样品的运动信息，应用于样品的投影断层重建。

为实现上述目的，本发明的具体步骤如下：

(1)获取投影数据

1a)照射源对固定在电控旋转台上的样品进行水平投影断层扫描；

1b)使用探测器采集样品的投影数据。

(2)计算机计算下式中投影数据的零阶距：

C_{0} (t_{i}) = {&Integral;}_{- \infty}^{\infty} g (θ (t_{i}), l) dl

其中，C₀(t_i)为投影数据的零阶距，

为正负无穷范围内的一维积分，g(θ(t_i)，l)为步骤(1)中采集的投影数据，θ(t_i)为在t_i时刻的样品的扫描角度，l为探测器各个像素与探测器中心像素的距离。

(3)计算机计算下式中投影数据的一阶距：

C_{1} (t_{i}) = {&Integral;}_{- \infty}^{\infty} g (θ (t_{i}), l) ldl

其中，C₁(t_i)为投影数据的一阶距，

(4)计算机计算下式中投影数据的质心：

Q_{i} = \frac{C_{1} (t_{i})}{C_{0} (t_{i})}

其中，Q_i为投影数据的质心，C₁(t_i)为投影数据的一阶距，C₀(t_i)为投影数据的零阶距。

(5)获取运动参数

5a)计算机按照下式建立投影数据质心与扫描样品质心在探测器上投影的一致性条件方程：

Q_{i} = (p_{1,0} + p_{1,1} t_{i} + . . . + p_{1, N} t_{i}^{N}) \cos θ (t_{i}) + (p_{2,0} + p_{2,1} t_{i} + . . . + p_{2, N} t_{i}^{N}) \sin θ (t_{i})

其中，Q_i为投影数据的质心，p_1，0，p_1，1，...，p_1，N，p_2，0，p_2，1，...，p_2，N是待确定的运动参数，t_i为扫描时刻，N为运动参数多项式的阶数，cosθ(t_i)为t_i时刻扫描角度的余弦值，sinθ(t_i)为t_i时刻扫描角度的正弦值；

5b)当投影数据个数大于待定的运动参数的个数时，使用最小二乘法从步骤5a)中得到的一致性条件方程中获取运动参数。

(6)当样品的运动曲线是光滑连续的并且样品只发生平动不发生转动时，计算机按照下面的公式计算每个扫描时刻样品的运动量：

\{\begin{matrix} d_{x} (t_{i}) = p_{1,1} t_{i} + . . . + p_{1, N} {t_{i}}^{N} \\ d_{y} (t_{i}) = p_{2,1} t_{i} + . . . + p_{2, N} {t_{i}}^{N} \end{matrix}

其中，d_x(t_i)和d_y(t_i)分别为样品的水平运动分量和竖直运动分量的多项式表达式，t_i为扫描时刻，p_1，1，...，p_1，N，p_2，1，...，p_2，N为步骤(5)中获取的运动参数，N为运动参数多项式的阶数。

(7)运动伪影校正：将样品的投影重建图像中减去步骤(6)中获取的样品运动量，实现光学投影断层重建运动伪影校正。

本发明与现有的技术相比具有以下优点：

第一，本发明采用估计样品运动轨迹的方法对投影图像进行处理，由于重建时使用了所有的投影数据，克服了现有技术降低图像信噪比的缺点，使得本发明重建图像的信噪比较高。

第二，本发明估计样品运动轨迹时使用的是不基于特征点的运动估计方法，克服了现有技术需要明显易辨识特征点的缺点，能够直接从投影数据中提取样品的运动参数信息，不要求投影数据中具有明显易于辨识的特征点，因此本发明使用起来比较方便、快捷。

附图说明

图1为本发明的流程图；

图2为本发明获取的样品运动量随时间变化的估计曲线；

图3为未经运动伪影校正的重建结果；

图4为本发明实施例的重建结果示意图。

具体实施方式

下面结合附图对本发明作进一步的详细描述。

结合附图1对本发明的步骤作进一步的描述。

步骤1，获取投影数据：

首先，照射源对固定在电控旋转台上的样品进行水平投影断层扫描，照射源采用激光器，并使用远心透镜将光线扩束为平行光照射样品。

然后，使用探测器采集样品的投影数据。本发明采集的方法是：样品每旋转一个角度采集一个投影数据，将每次采集的投影数据记录在计算机中。

步骤2，计算机计算下式中投影数据的零阶距：

C_{0} (t_{i}) = {&Integral;}_{- \infty}^{\infty} g (θ (t_{i}), l) dl

其中，C₀(t_i)为投影数据的零阶距，

为正负无穷范围内的一维积分，g(θ(t_i)，l)为步骤1中采集的投影数据，θ(t_i)为在t_i时刻的样品的扫描角度，l为探测器各个像素与探测器中心像素的距离。

步骤3，计算机计算下式中投影数据的一阶距：

C_{1} (t_{i}) = {&Integral;}_{- \infty}^{\infty} g (θ (t_{i}), l) ldl

其中，C₁(t_i)为投影数据的一阶距，

步骤4，计算机计算下式中投影数据的质心：

Q_{i} = \frac{C_{1} (t_{i})}{C_{0} (t_{i})}

步骤5，获取运动参数

首先，计算机按照下式建立投影数据质心与扫描样品质心在探测器上投影的一致性条件方程：

Q_{i} = (p_{1,0} + p_{1,1} t_{i} + . . . + p_{1, N} t_{i}^{N}) \cos θ (t_{i}) + (p_{2,0} + p_{2,1} t_{i} + . . . + p_{2, N} t_{i}^{N}) \sin θ (t_{i})

其中，Q_i为投影数据的质心，p_1，0，p_1，1，...，p_1，N，p_2，0，p_2，1，...，p_2，N是待确定的运动参数，t_i为扫描时刻，N为运动参数多项式的阶数，cosθ(t_i)为t_i时刻扫描角度的余弦值，sinθ(t_i)为t_i时刻扫描角度的正弦值。

然后，当投影数据个数大于待定的运动参数的个数时，使用最小二乘法从一致性条件方程中获取运动参数。最小二乘法是通过最小化采集的数据与实际值之间的误差平方和来获取运动参数。

步骤6，当样品的运动曲线是光滑连续的并且样品只发生平动不发生转动时，计算机按照下面的公式计算每个扫描时刻样品的运动量：

\{\begin{matrix} d_{x} (t_{i}) = p_{1,1} t_{i} + . . . + p_{1, N} {t_{i}}^{N} \\ d_{y} (t_{i}) = p_{2,1} t_{i} + . . . + p_{2, N} {t_{i}}^{N} \end{matrix}

其中，d_x(t_i)和d_y(t_i)分别为样品的水平运动分量和竖直运动分量的多项式表达式，t_i为扫描时刻，p_1，1，...，p_1，N，p_2，1，...，p_2，N为步骤5中获取的运动参数，N为运动参数多项式的阶数。

步骤7，运动伪影校正：将样品的投影重建图像中减去步骤6中获取的样品运动量，实现光学投影断层重建运动伪影校正。

下面结合附图2、附图3、附图4对本发明的重建结果做进一步的描述。

附图2为本发明获取的样品运动量随时间变化的估计曲线。其中，横坐标为时间，单位为秒；纵坐标为距离，单位为毫米；红色虚线为样品运动量水平分量随时间的变化曲线；蓝色实线为样品运动量垂直分量随时间的变化曲线。

附图2是本发明的步骤6获取的样品运动量在时间-距离坐标系中做出的随时间变化的曲线。附图2说明本发明使用的不基于特征点的运动估计方法，不要求投影数据中具有明显易于辨识的特征点，能够直接从投影数据中提取样品的运动参数信息获取样品的运动量。

附图3是未经运动伪影校正的重建结果。其中，扫描样品是固定在一个毛细玻璃管中的线虫，旋转角度为360°，采集的投影数据的个数是500个。

附图4是本发明实施例的重建结果示意图。其中，扫描样品是固定在一个毛细玻璃管中的线虫，旋转角度为360°，采集的投影数据的个数是500个。

附图4是用本发明的步骤进行了伪影校正后获取的重建结果。

将本发明的重建效果附图4和未经运动伪影校正的重建结果附图3进行比较，可以看出重建结果中的运动伪影得到了校正，清晰的得到了样品的重建结果，说明本发明克服了现有技术对投影数据有限制和需要明显易辨识特征点的缺点，方便、快捷的完成光学投影断层成像的运动伪影校正。

Claims

1.一种光学投影断层成像运动伪影校正方法，包括如下步骤：

(1)获取投影数据

1b)使用探测器采集样品的投影数据；

(2)计算机计算下式中投影数据的零阶距：

C_{0} (t_{i}) = {&Integral;}_{- \infty}^{\infty} g (θ (t_{i}), l) dl

其中，C₀(t_i)为投影数据的零阶距，

为正负无穷范围内的一维积分，g(θ(t_i)，l)为步骤(1)中采集的投影数据，θ(t_i)为在t_i时刻的样品的扫描角度，l为探测器各个像素与探测器中心像素的距离；

(3)计算机计算下式中投影数据的一阶距：

C_{1} (t_{i}) = {&Integral;}_{- \infty}^{\infty} g (θ (t_{i}), l) ldl

其中，C₁(t_i)为投影数据的一阶距，

(4)计算机计算下式中投影数据的质心：

Q_{i} = \frac{C_{1} (t_{i})}{C_{0} (t_{i})}

其中，Q_i为投影数据的质心，C₁(t_i)为投影数据的一阶距，C₀(t_i)为投影数据的零阶距；

(5)获取运动参数

Q_{i} = (p_{1,0} + p_{1,1} t_{i} + . . . + p_{1, N} t_{i}^{N}) \cos θ (t_{i}) + (p_{2,0} + p_{2,1} t_{i} + . . . + p_{2, N} t_{i}^{N}) \sin θ (t_{i})

5b)当投影数据个数大于待定的运动参数的个数时，使用最小二乘法从步骤5a)中得到的一致性条件方程中获取运动参数；

\{\begin{matrix} d_{x} (t_{i}) = p_{1,1} t_{i} + . . . + p_{1, N} {t_{i}}^{N} \\ d_{y} (t_{i}) = p_{2,1} t_{i} + . . . + p_{2, N} {t_{i}}^{N} \end{matrix}

其中，d_x(t_i)和d_y(t_i)分别为样品的水平运动分量和竖直运动分量的多项式表达式，t_i为扫描时刻，p_1，1，...，p_1，N，p_2，1，...，p_2，N为步骤(5)中获取的运动参数，N为运动参数多项式的阶数；

2.根据权利要求1所述的光学投影断层成像运动伪影校正方法，其特征在于：照射源采用激光器。

3.根据权利要求1所述的光学投影断层成像运动伪影校正方法，其特征在于：步骤1b)所述的采集方法是：样品每旋转一个角度采集一个投影数据，将每次采集的投影数据记录在计算机中。

4.根据权利要求1所述的光学投影断层成像运动伪影校正方法，其特征在于：步骤5b)中所述的最小二乘法是通过最小化采集的数据与实际值之间的误差平方和来获取运动参数。