CN102512140A - 一种定位光学投影断层成像旋转中心的方法 - Google Patents

Abstract

本发明公开了一种光学投影断层成像旋转中心的定位方法。该方法通过将X射线计算机断层成像领域中的基于质心轨迹的旋转中心计算方法引入光学投影断层成像领域，同时结合基于重建图像均方差的旋转中心精细定位方法，能针对不同的实验情况，实现自动、快速、精确的旋转中心定位。

Description

一种定位光学投影断层成像旋转中心的方法

技术领域

本发明涉及到光学投影断层成像(Optical Projection Tomography，简称OPT)技术，特别涉及到一种定位光学投影断层成像旋转中心的方法。 

背景技术

光学投影断层成像技术是利用光线在小尺寸生物体中沿直线传播的特点，发射可见光线穿透样品，然后用相机采集多个角度的样品投影视图，进行三维成像。具体来说，在进行光学投影断层成像时，需要对样品进行多角度扫描，一般采用电控转台对样品进行步进式旋转，每旋转到一个角度采集一幅或多幅投影图像，当样品旋转一周时停止扫描。光学投影断层成像***最终采集到的数据是一系列不同角度下光线穿过样品的二维投影图像，如果将所有投影图像的某一行都提取出来，按照扫描顺序依次按行叠加为一幅图像，就可以得到一个类似正弦曲线的正弦图，每一幅正弦图对应了样品的一张水平重建截面，所有正弦图就对应了样品的三维断层重建体，从正弦图到样品三维断层结构的过程称为光学投影断层成像三维重建。 

光学投影断层成像技术可以实现1-10mm尺度生物样品的结构和分子特异性功能成像，具有分辨率高、结构功能一体化、无辐射、成本低等诸多优点，它可以在小尺度对活体生物进行细胞水平的定性和定量研究，实现生物体的实时、无创、动态、在体成像。但是光学投影断层成像技术在三维成像之前需要预先定位转台旋转轴在探测器上的位置，即旋转中心位置，才能够进行高质量的三维成像。 

目前主流的旋转中心计算方法有四种：第一种方法是样品实验之前通过人工调整转台，使转台旋转中心轴尽量与探测器中心线对齐，其缺点是需要人工参与，准确性难以保证，准备时间长；第二种方法是样品实验前预先扫描已知仿体，利用仿体的特性计算旋转中心，再进行实验样品的扫 描，其缺点是操作繁琐，准备时间长；第三种方法是利用样品扫描获得的正弦图上的高亮曲线计算旋转中心，其缺点是这些高亮曲线并不是在每个样品上都存在的，此方法不具有通用性；第四种方法是对扫描样品得到的正弦图，人工给定一个初始的旋转中心，然后在其附近进行多次重建，人工挑选或程序识别出重建效果较好的位置，其缺点是需要人工给定初始值，操作不便，同时如果初始值距真实值过远，此方法会出现极大误差。综上所述，现有旋转中心计算方法需要人工参与，操作复杂，精度及速度都需改进。 

发明内容

(一)要解决的技术问题 

为解决上述的一个或多个问题，本发明提供了一种光学投影断层成像旋转中心的定位方法，以实现旋转中心的自动定位，提高定位的速度与精度。 

(二)技术方案 

根据本发明的一个方面，提供了一种光学投影断层成像旋转中心的定位方法。该定位方法包括：对于光学投影断层成像的多幅选定正弦图中的每一幅，利用质心轨迹法获取该幅正弦图的旋转中心；将多幅选定正弦图的旋转中心进行平均，获取光学投影断层成像旋转中心的粗定位值；以粗定位值的旋转中心为初始值，设置旋转中心精细搜索区间，在该精密搜索区间对选定正弦图中的每一幅进行多次重建，在一系列重建图像中，选取重建图像均方差最大的旋转中心位置作为该幅正弦图的再定旋转中心，将所有选定正弦图的再定旋转中心进行平均获取光学投影断层成像精细定位旋转中心。 

(三)有益效果 

本发明光学投影断层成像旋转中心的定位方法，充分保障了定位的快速性、自动型、准确性，具体来讲： 

(1)在正弦图选择上，充分考虑到光学投影断层成像扫描数据的噪声影响，仅选择样品信号强的正弦图上进行旋转中心搜索，可以保证旋转中心搜索所基于的数据具有较高信噪比，从数据选择上提高了本发明的搜 索精度； 

(2)在粗定位过程中，直接利用质心轨迹进行自动旋转中心估计，粗定位过程的质心方法仅需要投影数据，不需要人为干预，充分保证了本发明的自动性； 

(3)在细定位过程中，基于粗定位的结果进行再次搜索，并且采用并行处理技术，精细定位在粗定位值附近进行搜索，降低了本发明的误定位概率，保障了最终定位精度，同时也减少了搜索范围，结合并行技术保障了本发明的快速性。 

附图说明

图1为本发明实施例光学投影断层成像旋转中心的定位方法中高样品信号强度正弦图筛选流程图； 

图2为本发明实施例光学投影断层成像旋转中心的定位方法中基于质心轨迹的旋转中心粗定位流程； 

图3为本发明实施例光学投影断层成像旋转中心的定位方法中基于重建图像均方差的旋转中心精细定位流程； 

图4为本发明实施例光学投影断层成像旋转中心的定位方法中针对果蝇蛹实验的旋转中心计算过程示意图；其中，图4a显示了各个选定正弦图用质心轨迹法计算得到的旋转中心和平均粗定位旋转中心；图4b显示了用粗定位旋转中心进行重建得到了三维重建体中的截面；图4c显示了在粗定位中心附近用重建图像均方差法计算得到的再定旋转中心和精细定位旋转中心；图4d显示了利用精细定位旋转中心进行重建得到的三维重建体截面； 

图5为本发明实施例光学投影断层成像旋转中心的定位方法中针对果蝇蛹实验的重建结果。 

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。虽然本文可提供包含特定值的参数的示范，但应了解，参数无需确切等于相应的值，而是可在可接 受的误差容限或设计约束内近似于所述值。 

本发明通过将X射线计算机断层成像领域中的基于质心轨迹的旋转中心计算方法引入光学投影断层成像领域，同时结合基于重建图像均方差的旋转中心精细定位方法，设计出一种光学投影断层成像旋转中心的定位方法，能针对不同的实验情况，实现自动、快速、精确的旋转中心定位。 

本发明的实现方案分为三个主要步骤：高样品信号强度正弦图筛选、基于质心轨迹的旋转中心粗定位、基于重建图像均方差的旋转中心精细定位，其中粗定位和细定位都是在第一步选定的正弦图上进行，以提高定位精度。下面利用果蝇蛹实验描述本发明的步骤，实验中采用活体蛹进行光学投影断层成像，成像***采集的投影视图像素数为500*500，每个像素大小为16微米，样品旋转一周均匀采集的角度数为330，以活体蛹实验为例，本发明的详细步骤如下： 

步骤S1：为提高旋转中心定位精度，本发明预先对光学投影断层成像获得的正弦图进行筛选，提取噪声小、信号强度好的正弦图进行后续计算。 

为了衡量正弦图是否具有较高的信噪比，本发明利用正弦图样品像素个数来相对衡量其信噪比，样品像素个数越多的正弦图，其信号强度越好，噪声越低，进行旋转中心搜索效果越佳。但如果计算所有正弦图的样品像素个数，计算复杂，本发明通过简化处理，仅统计每个正弦图中相互垂直投影角度的两行样品像素平均数，相对定量正弦图的样品信号强度，在较小误差的情况下，大大减少了处理时间。当然，对所有的正弦图进行后续计算也是可以的，只不过耗费的软硬件资源更大一些。 

图1为本发明实施例光学投影断层成像旋转中心的定位方法中高样品信号强度正弦图筛选流程图。以下结合图1，对筛选流程进行详细说明。 

步骤S1-1：如图1(a列)所示，本发明首先筛选出两幅投影光线方向相垂直的投影视图，本发明选择0°和90°投影视图，这两个投影视图的相同行可映射到对应正弦图中投影角度相垂直的两行，用这两个投影图便可评估所有正弦图的样品信号强度；本领域技术人员应当了解，此处选择0°和90°投影视图为最优的方案，实际上，也可以选择三幅或更多幅的投影图来评估正弦图的样品信号强度，而投影角度也可以任意选择，当然，最好是选择具代表性的多幅投影图。 

步骤S1-2：如图1(b列)所示，对选择出两个投影视图进行阈值分割，区分样品和背景像素，进而评估各正弦图的样品信号强度，为了简化计算，设置分割阈值为对应投影视图的像素灰度均值，这样可以简便有效地区分出样品与背景； 

步骤S1-3：在两幅投影图上统计每一行样品像素的个数，其代表该行的相对信号强度，以样品像素个数为横轴，行号为纵轴，可获得如图1(c列)所示的投影视图各行样品信号强度曲线，曲线的高低代表该行样品信号的强弱； 

步骤S1-4：将两幅投影图的信号强度曲线进行平均，得到如图1(d列)所示的平均信号强度曲线，通过之前的表述可知，该曲线对应了每个正弦图的样品信号强度； 

步骤S1-5：如图1(d列)所示，筛选出信号强度在前25％的正弦图，图中的灰色区域为筛选出的正弦图，在果蝇蛹实验中其对应的是第216至第340个正弦图，然后，将这些正弦图进行后续的旋转中心计算，如图1(e列)所示。 

步骤S2：针对步骤S1选定的每一个正弦图，计算质心投影在探测器上的一系列位置，并计算其均值，即得到正弦图的旋转中心，将所有用质心轨迹计算的旋转中心进行平均即为本发明的粗定位旋转中心。图2为本发明实施例光学投影断层成像旋转中心的定位方法中基于质心轨迹的旋转中心粗定位流程。如图2所示，针对每个正弦图，按照下面公式计算其粗定位旋转中心位置 

$\stackrel{\‾}{s_c}=\frac{1}{M}{\mathrm{\Σ}}_{i=0}^{M-1}{s}_m\left({\mathrm{\γ}}_i\right)$

其中M为光学投影断层成像的投影视图数，即投影角度数；γ_i为光学投影断层成像的投影角度，其计算公式为 

i＝0，...，(M-1)；s_m(γ_i)为正弦图在γ_i角度质心投影在探测器上的位置，其计算公式如下： 

$s_m\left({\mathrm{\γ}}_i\right)=\frac{{\mathrm{\Σ}}_{j=1}^N{s}_{j}g(s_j,{\mathrm{\γ}}_i)}{{\mathrm{\Σ}}_{j=1}^{N}g(s_j,{\mathrm{\γ}}_i)}$

其中，s_j为探测器各列的中心位置，N为探测器的列数，g(s_j，γ_i)为正弦图在γ_i角度s_j列处的像素值，本步骤利用质心法轨迹法计算出每个正弦图的旋转中心位置 

将所有用质心轨迹法计算出的正弦图旋转中心进行平均，得到本实施例的粗定位旋转中心。 

步骤S3：以步骤S2粗定位的旋转中心为初始值，设置旋转中心精细搜索区间，在此区间中，对选定的正弦图进行多次重建，在一系列重建图像中，选取重建图像均方差最大的旋转中心位置，即为该幅正弦图的再定旋转中心，将所有正弦图的再定旋转中心进行平均即为本发明的精细定位旋转中心。 

图3为本发明实施例光学投影断层成像旋转中心的定位方法中基于重建图像均方差的旋转中心精细定位流程。图3的上半部分为基于不同的备选旋转中心位置所对应的重建图像；图3的下半部分为针对不同备选旋转中心位置重建图像的均方差值。从位置上讲，上半部分的重建图，与其下方虚线指向均方差的值存在对应关系。 

如图3所示，针对粗定位的旋转中心位置，设置其左右各20个像素范围为旋转中心搜索区间，设置搜索步进为1/8像素值，筛选出320个备选旋转中心位置。针对每个正弦图，利用320个备选旋转中心分别对其进行重建，获得320个重建图像，计算每一个重建图像的均方差F_variance，其计算公式为： 

$F_{\mathrm{variance}}=\frac{{\mathrm{\Σ}}_{y=0}^{N-1}{\mathrm{\Σ}}_{x=0}^{N-1}{(f(x,y)-\stackrel{\‾}{f})}^2}{N^2}$

其中f(x，y)为重建图像在坐标(x，y)处的像素值，N为重建图像的行数和列数， 

为重建图像的像素均值，其计算公式为 

通过上述公式计算出每个正弦图重建图像均方差最大的旋转中心位置，即为该幅正弦图的再定旋转中心，如图3的下半部分可见基于再定旋转中心的重建 结果要优于步骤S2粗定位旋转中心的重建结果。最后，将所有正弦图的再定位旋转中心位置进行平均，获得本发明的精细定位旋转中心。 

针对本发明的旋转中心计算方法，图4为本发明实施例光学投影断层成像旋转中心的定位方法中针对果蝇蛹实验的旋转中心计算过程示意图。其中，图4a显示了各个选定正弦图用质心轨迹法计算得到的旋转中心和平均粗定位旋转中心；图4b显示了用粗定位旋转中心进行重建得到了三维重建体中的截面，可见图像边缘处仅有微弱的伪影，说明粗定位旋转中心已接近真实旋转中心位置；图4c显示了在粗定位中心附近用重建图像均方差法计算得到的再定旋转中心和精细定位旋转中心，可见精细定位旋转中心与粗定位旋转中心仅差两个像素，从侧面证明了两步定位的有效性；图4d显示了利用精细定位旋转中心进行重建得到的三维重建体截面，其相比图4b清晰度得到了极大提高。 

图5为本发明实施例光学投影断层成像旋转中心的定位方法中针对果蝇蛹实验的重建结果。其中，(a)-(d)部分显示了利用本发明在果蝇蛹实验中的一些结果，其中(a)显示了一幅蛹的投影视图；(b)为利用本发明的旋转中心计算方法校正后的三维重建体截面；(c)为三维重建体进行可视化的结果；(d)为果蝇蛹的结构图谱，从图5可见，通过本发明校准旋转中心后的三维重建结果，其内部结构与蛹的结构图谱完全吻合，从实验上验证了本发明的有效性。 

综上所述，本发明提出了一种由粗到细的光学投影断层成像旋转中心的定位方法，充分保障了定位的快速性、自动型、准确性，具体来讲： 

(1)在正弦图选择上，充分考虑到光学投影断层成像扫描数据的噪声影响，仅选择样品信号强的正弦图上进行旋转中心搜索，可以保证旋转中心搜索所基于的数据具有较高信噪比，从数据选择上提高了本发明的搜索精度； 

(2)在粗定位过程中，直接利用质心轨迹进行自动旋转中心估计，粗定位过程的质心方法仅需要投影数据，不需要人为干预，充分保证了本发明的自动性； 

(3)在细定位过程中，基于粗定位的结果进行再次搜索，并且采用 并行处理技术，精细定位在粗定位值附近进行搜索，降低了本发明的误定位概率，保障了最终定位精度，同时也减少了搜索范围，结合并行技术保障了本发明的快速性。 

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。 

Claims (12)

1.一种光学投影断层成像旋转中心的定位方法，其特征在于，包括：

对于光学投影断层成像的多幅选定正弦图中的每一幅，利用质心轨迹法获取该幅正弦图的旋转中心；将所述多幅选定正弦图的旋转中心进行平均，获取光学投影断层成像旋转中心的粗定位值；

以所述旋转中心的粗定位值为初始值，设置旋转中心精细搜索区间，在该精密搜索区间对选定正弦图中的每一幅进行多次重建，在一系列重建图像中，选取重建图像均方差最大的旋转中心位置作为该幅正弦图的再定旋转中心，将所有选定正弦图的再定旋转中心进行平均获取光学投影断层成像精细定位旋转中心。

2.根据权利要求1所述的定位方法，其特征在于，所述利用质心轨迹法获取该幅正弦图的旋转中心的步骤中，按照以下公式计算正弦图的旋转中心位置

${\stackrel{\‾}{s}}_c=\frac{1}{M}{\mathrm{\Σ}}_{i=0}^{M-1}{s}_m\left({\mathrm{\γ}}_i\right),$

其中M为光学投影断层成像的投影视图数，即投影角度数；γ_i为光学投影断层成像的投影角度，其计算公式为

i＝0，...，(M-1)；s_m(γ_i)为正弦图在γ_i角度质心投影在探测器上的位置。

3.根据权利要求2所述的定位方法，其特征在于，所述计算正弦图的旋转中心位置

的步骤中，按照以下公式计算所述s_m(γ_i)的值：

$s_m\left({\mathrm{\γ}}_i\right)=\frac{{\mathrm{\Σ}}_{j=1}^N{s}_{j}g(s_j,{\mathrm{\γ}}_i)}{{\mathrm{\Σ}}_{j=1}^{N}g(s_j,{\mathrm{\γ}}_i)}$

其中，s_j为探测器各列的中心位置，N为探测器的列数，g(s_i，γ_i)为正弦图在γ_i角度s_j列处的像素值。

4.根据权利要求1所述的定位方法，其特征在于，所述选取重建图像均方差最大的旋转中心位置作为该幅正弦图的再定旋转中心的步骤中，按照以下公式计算重建图像的均方差：

$F_{\mathrm{variance}}=\frac{{\mathrm{\Σ}}_{y=0}^{N-1}{\mathrm{\Σ}}_{x=0}^{N-1}{(f(x,y)-\stackrel{\‾}{f})}^2}{N^2}$

其中f(x，y)为重建图像在坐标(x，y)处的像素值，N为重建图像的行数和列数，

为重建图像的像素均值，其计算公式为

5.根据权利要求4所述的定位方法，其特征在于，所述在该精密搜索区间对选定的正弦图进行多次重建，在一系列重建图像中，选取重建图像均方差最大的旋转中心位置作为该幅正弦图的再定旋转中心的步骤中，采用并行方式进行正弦图重建，并获取每一重建图像的均方差。

6.根据权利要求4所述的定位方法，其特征在于，所述以粗定位值为初始值，设置旋转中心精细搜索区间的步骤包括：

以所述粗定位中心为中点，设置其左右各N个像素范围为旋转中心精细搜索区间；

设置搜索步进为M像素值，筛选出[2N/M]个精细搜索点，所述N大于M，[  ]为取整符号。

7.根据权利要求6所述的定位方法，其特征在于，所述N＝20，所述M＝1/8。

8.根据权利要求1至7中任一项所述的定位方法，其特征在于，所述对于光学投影断层成像的多幅选定正弦图中的每一幅正弦图，利用质心轨迹法获取该幅正弦图的旋转中心的步骤之前还包括：

从光学投影断层成像的全部正弦图中，选取信号强度强的正弦图作为所述选定正弦图。

9.根据权利要求8所述的定位方法，其特征在于，所述选取信号强度强的正弦图作为选定正弦图的步骤包括：

提取所述光学投影断层成像全部正弦图的相对信号强度；

按相对信号强度从强到弱的顺序对所述全部正弦图进行排序；

将排序靠前的正弦图作为所述选定正弦图。

10.根据权利要求9所述的定位方法，其特征在于，所述提取光学投影断层成像全部正弦图的信号强度的步骤包括：

从所述光学投影断层成像中筛选出至少两幅投影视图；

对所述至少两幅投影视图进行阈值分割，去除背景区域的像素；

计算去除背景区域像素后的所述至少两幅投影中每行的样品像素个数，将所述样品像素个数作为该行的相对信号强度；

将所述至少两幅投影中对应行的相对信号强度平均，将所述平均值作为对应该行的正弦图的相对信号强度。

11.根据权利要求10所述的定位方法，其特征在于，所述至少两幅投影视图为两幅投影视图，所述两幅投影视图的投影光线方向相互垂直。

12.根据权利要求10所述的定位方法，其特征在于，所述对至少两幅投影视图进行阈值分割的步骤中，所述阈值为投影视图的像素灰度均值。

Images