《一种大数据量荧光分子断层成像重建方法.pdf》由会员分享,可在线阅读,更多相关《一种大数据量荧光分子断层成像重建方法.pdf(10页珍藏版)》请在专利查询网上搜索。
1、(10)申请公布号 CN 102871646 A (43)申请公布日 2013.01.16 CN 102871646 A *CN102871646A* (21)申请号 201210293156.4 (22)申请日 2012.08.16 A61B 5/00(2006.01) G06T 9/00(2006.01) (71)申请人 清华大学 地址 100084 北京市海淀区 100084 信箱 82 分箱清华大学专利办公室 (72)发明人 白净 曹旭 张宾 刘飞 王鑫 (74)专利代理机构 北京纪凯知识产权代理有限 公司 11245 代理人 徐宁 关畅 (54) 发明名称 一种大数据量荧光分子断层成。
2、像重建方法 (57) 摘要 本发明涉及一种大数据量荧光分子断层成像 重建方法, 包括以下步骤 : 采用的全角度自由空 间 FMT 成像系统采集待成像物体的荧光图像和 白光图像 ; 采用边缘检测方法提取每一幅白光图 像中待成像物体的边界轮廓线, 得到投影轮廓线 图像 ; 对每一幅投影轮廓线图像均采用滤波反投 影方法依次进行反投影得到待成像物体的三维轮 廓图像 ; 采用有限元方法求解扩散方程的格林函 数 ; 将每一幅荧光图像划分为若干相同大小的子 荧光图像 ; 采用所求得的格林函数建立每一个子 荧光图像所对应的子系统方程 ; 对每一幅荧光图 像划分的每一子荧光图像所对应的子系统方程压 缩后依次按行。
3、排列得到压缩系统方程 ; 求解压缩 系统方程得到待成像体内的荧光标记物分布。本 发明可以广泛应用于大数据量的荧光分子断层成 像重建中。 (51)Int.Cl. 权利要求书 1 页 说明书 6 页 附图 2 页 (19)中华人民共和国国家知识产权局 (12)发明专利申请 权利要求书 1 页 说明书 6 页 附图 2 页 1/1 页 2 1. 一种大数据量荧光分子断层成像重建方法, 包括以下步骤 : 1) 设置一包括有旋转台、 激发光发射装置、 CCD 相机和计算机的全角度自由空间 FMT 成 像系统 ; 2) 采用所述全角度自由空间 FMT 成像系统采集待成像物体的荧光图像和白光图像 ; 3) 。
4、采用边缘检测方法提取每一幅白光图像中待成像物体的边界轮廓线, 得到投影轮廓 线图像 ; 4) 对步骤 3) 中的每一幅投影轮廓线图像均采用滤波反投影方法依次进行反投影得到 待成像物体的三维轮廓图像 ; 5) 采用有限元方法求解扩散方程的格林函数 ; 6) 将每一幅荧光图像划分为相同大小的子荧光图像 ; 7) 根据一阶波恩近似理论, 采用所求得的格林函数建立每一个子荧光图像所对应的子 系统方程, 并将子系统方程的子系统矩阵和子荧光向量组成子增广系统矩阵 ; 8) 对每一子增广系统矩阵采用主成分分析方法进行降维压缩处理得到压缩子系统方 程 ; 9) 对每一幅荧光图像划分的每一子荧光图像所对应的子系。
5、统方程采用步骤 8) 降维压 缩后, 将所有压缩子系统方程依次按行排列得到压缩系统方程 ; 10) 求解压缩系统方程得到待成像体内的荧光标记物分布。 2. 如权利要求 1 所述的一种大数据量荧光分子断层成像重建方法, 包括以下步骤 : 所 述步骤 10) 采用 tikhonov 正则化方法求解压缩系统方程。 权 利 要 求 书 CN 102871646 A 2 1/6 页 3 一种大数据量荧光分子断层成像重建方法 技术领域 0001 本发明涉及一种荧光分子断层成像重建方法, 特别是关于一种基于数据降维技术 的大数据量荧光分子断层成像重建方法。 背景技术 0002 荧光分子断层成像 (Fluor。
6、escence Molecular Tomography,FMT)是一种新兴的光 学断层成像技术, 它是一种高灵敏度、 无电离辐射、 低成本的在体小动物成像方法, 在肿瘤 研究、 药物研发和疾病诊断等领域有着广阔的应用前景。荧光分子断层成像技术是利用荧 光标记物标记小动物体内的特定分子或细胞, 采用合适波段和强度的激发光照射被标记的 小动物时, 小动物体内的荧光标记物发出荧光, 采用一定的装置检测此荧光信号, 并采用相 应的重建方法就可以获得小动物体内的荧光标记物的分布情况, 从而可以在分子和细胞水 平上对正常或异常的生物过程进行观察。 0003 不同于传统的医学断层影像 (CT 或 PET)。
7、 , 由于光在生物组织中传播时会受到生物 组织的强烈散射和吸收作用, 使得荧光分子断层成像的重建成为一个严重的病态问题, 这 种病态问题主要体现在重建图像的分辨率低, 稳定性差, 且易受噪声的干扰, 通过对采集的 大量投影数据进行重建能够有效地克服这种病态问题, 有效提高重建图像的分辨率和稳定 性。随着荧光分子断层成像系统的发展, 电荷耦合器件 (CCD) 被用于投影数据的采集, 使得 每一个投影可以采集到大量的投影数据 (例如, 512*512像素的CCD可以采集262144个投影 数据) , 为了方便地获取更多有效信息, 全角度自由空间的成像模式被广泛应用于荧光分子 断层成像系统, 这种全。
8、角度 (360采集, 通常每 10一个投影, 总共 36 个投影) 采集使得采 集的投影数据急剧增加, 采用大量的投影数据 (大约 107) 进行荧光分子断层重建将会有效 地提高重建图像的分辨率和稳定性, 但是这种大数据量荧光分子断层成像重建非常消耗计 算机内存 (大约 102GB) 和计算时间, 如此巨大地内存消耗量远远超出目前普通 PC 机 (个人 电脑) 的内存容量, 难以在普通 PC 机上实现, 即使随着个人电脑的发展, 内存不够不再是瓶 颈, 巨大的计算时间将严重的制约快速、 实时荧光分子断层成像过程, 因此大数据量荧光分 子断层成像重建采用现有的方法几乎是无法完成的。 发明内容 0。
9、004 针对上述问题, 本发明的目的是提供一种能够有效降低投影数据的数据量, 减少 计算时间, 且能够提高重建图像的分辨率低和稳定性的大数据量荧光分子断层成像重建方 法。 0005 为实现上述目的, 本发明采取以下技术方案 : 一种大数据量荧光分子断层成像重 建方法, 包括以下步骤 : 1) 设置一包括有旋转台、 激发光发射装置、 CCD 相机和计算机的全 角度自由空间 FMT 成像系统 ; 2) 采用所述全角度自由空间 FMT 成像系统采集待成像物体的 荧光图像和白光图像 ; 3) 采用边缘检测方法提取每一幅白光图像中待成像物体的边界轮廓 线, 得到投影轮廓线图像 ; 4) 对步骤 3) 中。
10、的每一幅投影轮廓线图像均采用滤波反投影方法 说 明 书 CN 102871646 A 3 2/6 页 4 依次进行反投影得到待成像物体的三维轮廓图像 ; 5) 采用有限元方法求解扩散方程的格林 函数 ; 6) 将每一幅荧光图像划分为相同大小的子荧光图像 ; 7) 根据一阶波恩近似理论, 采 用所求得的格林函数建立每一个子荧光图像所对应的子系统方程, 并将子系统方程的子系 统矩阵和子荧光向量组成子增广系统矩阵 ; 8) 对每一子增广系统矩阵采用主成分分析方法 进行降维压缩处理得到压缩子系统方程 ; 9) 对每一幅荧光图像划分的每一子荧光图像所对 应的子系统方程采用步骤 8) 降维压缩后, 将所有。
11、压缩子系统方程依次按行排列得到压缩系 统方程 ; 10) 求解压缩系统方程得到待成像体内的荧光标记物分布。 0006 所述步骤 10) 采用 tikhonov 正则化方法求解压缩系统方程。 0007 本发明由于采取以上技术方案, 其具有以下优点 : 1、 本发明将采集的每一幅荧光 图像进行分块, 且将分块的每一子荧光图像所对应的子系统方程均进行数据降维压缩处理 后组建压缩系统矩阵, 例如, 如果采用 512*512 像素的 CCD 相机采集 36 幅荧光图像, 系统矩 阵 W 将达到 36*512*512=9437184 行 N(大约 104) 列, 如此大规模的矩阵已经超出一般计算 机的内存。
12、容量, 但是采用数据降维压缩处理后重新组建系统矩阵, 系统矩阵 W 将被极大地 压缩, 如果压缩率设置为 512, 上述系统矩阵 W 规模将减小为 36*256*2=18432 行 N 列, 因此 极大地降低了大数据量荧光分子断层成像重建所消耗的内存, 使得大数据量荧光分子断层 成像重建可以在普通 PC 机完成。2、 本发明由于采用数据降维压缩处理, 可以使系统矩阵得 到极大压缩, 因此进行荧光分子断层图像重建时可以有效地减少计算时间, 从而实现快速 的大数据量荧光分子断层成像重建。 3、 本发明在对系统矩阵进行降维压缩处理时采用主成 分分析方法, 因此可以最大限度地保留原始大量投影数据的有效。
13、信息, 从而实现以小数据 量规模的计算达到大数据量的重建效果, 提高重建图像的分辨率和稳定性。本发明可以广 泛应用于大数据量的荧光分子断层成像重建中。 附图说明 0008 图 1 是本发明所采用现有的全角度自由空间 FMT 成像系统结构示意图 ; 0009 图 2 是本发明的荧光分子断层成像重建方法流程示意图。 具体实施方式 0010 下面结合附图和实施例对本发明进行详细的描述。 0011 如图 1 所示, 本发明采用现有的全角度自由空间 FMT 成像系统对待成像物体进行 图像采集和处理, 它包括有旋转台 1、 CCD 相机 2、 激发光发射装置 3 和计算机 4, 其中激发 光发射装置 3 。
14、包括卤灯 31、 光纤 32 和激发光滤光片 33 ; 采集图像时, 待成像物体 5 放置在 旋转台 1, 卤灯 31 发射的激发光通过光纤 32 和激发光滤光片 33 传播并照射待成像物体 5, CCD 相机 1 完成待成像物体 1 图像的采集并将其发送到计算机 4 进行处理。 0012 如图 2 所示, 本发明的大数据量荧光分子断层成像重建方法, 包括以下步骤 : 0013 1、 将荧光标记物注射到待成像物体 5 内, 根据所注射的荧光标记物特征, 启动激 发光发射装置 3 发射某一波长和强度的激发光照射待成像物体 5, 此时荧光标记物发出荧 光, CCD 相机 2 配合此荧光波段的滤光片。
15、 6 以全角度、 等间隔模式采集待成像物体 5 的荧光 投影图像, 并将荧光投影图像定义为荧光图像, 在 360的采集范围内, 可以采用间隔角度 为 10进行图像采集, 得到 S=36 幅荧光图像, 也可以根据实际实验需要设置其它的间隔角 说 明 书 CN 102871646 A 4 3/6 页 5 度。 0014 2、 在自然光照条件下 (自然光是指通常的白光) , 关闭激发光发射装置 3, CCD 相机 2 直接以全角度、 等间隔模式采集待成像物体 5 的投影图像, 由于是在白光条件下对待成 像物体进行图像采集, 因此将采集的反映待成像物体几何模型的图像定义为白光图像, 在 360的采集范。
16、围内, 可以采用间隔角度为 10进行图像采集, 得到 36 幅白光图像, 也可以 根据实际实验需要设置其它的间隔角度。 0015 3、 采用边缘检测方法提取每一幅白光图像中待成像物体的边界轮廓线, 得到 36 幅投影轮廓线图像, 它们代表待成像物体 5 不同角度的投影轮廓线。边缘检测方法是根据 图像边界点亮度值变化明显的特征提取图像的边界信息, 边缘检测方法可以根据实际需要 进行选择, 在此不做限定, 只要能够实现上述功能即可。 0016 4、 对步骤 3 得到的每一幅投影轮廓线图像均采用滤波反投影方法依次进行反投 影得到待成像物体的三维轮廓图像, 滤波反投影方法广泛应用于计算机断层成像 (C。
17、T) 领 域, 用于从多角度的投影二维图像中计算出几何体的三维结构图像, 由于滤波反投影方法 为现有方法, 故此过程不再赘述。 0017 5、 采用有限元方法求解扩散方程的格林函数。 0018 FMT 中最广泛应用的数学模型是扩散方程, 它可以描述光在生物组织里的吸收和 散射作用, 其公式为 : 0019 0020 式中,为梯度算子, r 是空间点向量, 是成像空间, G(r) 是格林函数, 表示点光 源激发下生物组织内的光强分布 ; D1/(3(a+s)是生物组织的扩散参数, a和s是 生物组织吸收参数, Ks是激发点 rs处 (s=1S, S 是荧光图像个数) 的激发光强度,是待成 像物体。
18、的外法线向量, Cnd是一个常数, 反映光从待成像物体内通过体表的反射情况, 是 狄拉克函数。 0021 由于待成像物体的几何模型是不规则的, 扩散方程在不规则的待测成像物体几何 模型下是没有解析解的, 因此可以采用现有的有限元方法求解扩散方程 (1) , 有限元方法求 解扩散方程一般采用现有的有限元商业软件完成, 具体步骤为 : 0022 1) 将待成像物体的三维轮廓信息进行处理得到待成像物体的三维几何实体, 具体 过程是首先利用每一层的轮廓线创建其对应的二维几何面, 然后采用所有层的二维几何面 创建三维几何实体。 0023 2) 在步骤 1) 的待成像物体的三维实体内加载扩散方程 (1) 。
19、, 在待成像物体的三维 实体表面加载边界条件。 0024 3) 划分网格, 并组集刚度矩阵, 得到以下的线性方程 : 0025 KG bs(2) 0026 式中, K 是刚度矩阵, bs是激发向量, G 是离散格林函数。 0027 在求得离散格林函数 G 后, 根据一阶波恩近似, 可以建立荧光图像和待成像物体 内的荧光标记物分布的关系为 : 0028 m(rs,rd) G(rs,r)x(r)G(r,rd)d 3r (3) 说 明 书 CN 102871646 A 5 4/6 页 6 0029 式中, m(rs,rd)是激发光在点rs处激发时, 在点rd处所检测到的荧光值 ; x(r)是 待成像。
20、物体内的荧光标记物分布, 也就是 FMT 所要重建的未知数 ; G(rs,r) 是 rs点激发的格 林函数在 r 点的值, G(r,rd) 是 rd点激发的格林函数在 r 点的值。 0030 将离散格林函数 G 代入积分方程公式 (3) 中可以得到 : 0031 0032 式中, ri是待成像物体划分的网格节点, N 是几何体网格节点数, Vi是网格体积 大小, 为了描述方便定义 : 0033 W(rs,rd,ri) G(rs,ri)G(ri,rd) Vi (5) 0034 公式 (4) 可以写成 : 0035 0036 式中, W(rs,rd)是一个包含N个元素的行向量, 表示激发光在rs点。
21、激发时待成像物 体内各个节点对检测点 rd的贡献, 定义为系统向量, 可以写为 : 0037 W(rs,rd) W(rs,rd,rl)W(rs,rd,rN) (7) 0038 每一个检测点对应于一个系统向量, 当激发光在rs点激发时, CCD相机2可以采集 到一幅包含 L(L 是每一幅荧光图像的像素总数, 如果采用 512*512 像素的 CCD 相机采集, L=512*512) 个像素的荧光图像, 每一个像素对应于一个检测点 m(rs,rd)(d=1L) , 这样 就对应L个系统向量, 将这L个系统向量依次按行排列就可以组建成一个矩阵, 按照同样的 方法处理36个角度的投影图像, 得到36个。
22、矩阵, 将这些矩阵依次按行排列就可以组建成最 终的系统矩阵 W, 用于描述荧光投影数据和待求解的荧光标记物分布的数学关系, 相应的荧 光投影数据 m(rs,rd) 依次排列成一个列向量, 定义为荧光向量, 这样就得到 FMT 的系统方 程 : 0039 0040 式中, M 是 S 个荧光图像的像素总个数, 即所有监测点总个数 , 将公式 (8) 简写为 : 0041 m Wx (9) 0042 式中, W 是 FMT 的系统矩阵。随着投影角度和检测点的增多, W 会不断增大, 例如 : 36 个投影用 512*512 像素的 CCD 相机采集荧光图, 那么 W 将达到 36*512*512=。
23、9437184 行 N 列 (N 一般数量级在 104) , 如此大规模的矩阵已经超出一般的计算机的内存容量无法求解, 所以针对这种大数据量的FMT问题不方便直接组建系统矩阵W, 因此将W分块进行PCA压缩 后再组建一个降维压缩的系统矩阵。 0043 6、 对荧光图像进行分块, 将每一幅荧光图像按网格划分为 U 个大小相同的子荧光 图像, 例如 : 对 512*512 的像素荧光图像, 以网格形式可以被划分为 U=256 个大小相同的 32*32 像素子荧光图像 m(rs,Rsub),Rsub(sub=1U) 代表子荧光图像的 1024 个像素。 0044 7、 根据一阶波恩近似理论, 采用 。
24、FMT 的数学模型的解建立每一个子荧光图像所对 应的子系统方程, 并将子系统方程的子系统矩阵和子荧光向量组成一个子增广系统矩阵, 说 明 书 CN 102871646 A 6 5/6 页 7 具体过程为 : 0045 针对每一幅荧光图像的这一系列子荧光图像建立子系统方程为 : 0046 m(rs,Rsub) W(rs,Rsub)x (10) 0047 式中, m(rs,Rsub) 为子荧光向量, W(rs,Rsub) 为子系统矩阵。 0048 将子系统方程的子系统矩阵和子荧光向量组成一个子增广系统矩阵为 : 0049 A(rs,Rsub) W(rs,Rsub),m(rs,Rsub) (11) 。
25、0050 8、 对每一子增广矩阵采用主成分分析方法 (Principal Component Analysis, PCA) 进行降维压缩处理得到压缩子系统方程。 0051 PCA 是一种广泛应用的数据降维方法, 它可以将高维数据投影到较低维空间, 用低 维数据描述高维数据所包含的信息, 从而达到数据降维的目的, 公式 (11) 中的 A(rs,Rsub) 是 一个待降维的数据矩阵, 行是变量个数, 列是观测次数。计算 A(rs,Rsub) 的协方差矩阵的特 征值和特征向量, A(rs,Rsub) 的主成分是以前几个特征值所对应的特征向量组合成的特征 矩阵对数据矩阵的投影。为了有效的压缩 A(r。
26、s,Rsub), 这里只保留前几个主成分分量 : 0052 Wpca(rs,Tsub),pcam(rs,Tsub) PCAA(rs,Rsub) (12) 0053 式中, Wpca(rs,Tsub) 是 PCA 压缩后的子系统矩阵, pcam(rs,Tsub) 是压缩后的子荧光 向量, Tsub代表前几个主成分分量。 0054 通过上述压缩方法实现了 Rsub到 Tsub的压缩, 从而得到压缩子系统方程 : 0055 pcam(rs,Tsub) Wpca(rs,Tsub)x (13) 0056 9、 将每一幅荧光图像划分的 U 个子荧光图像所对应的子系统方程进行上述步骤 8 的降维压缩处理, 。
27、即得到S幅荧光图像所对应的S*U个压缩子系统方程, 将所有压缩子系统 方程依次按行排列得到压缩系统方程 : 0057 0058 式中, V 是子荧光图像总个数, 对于 S=36, U=256, V=S*U=9216。 0059 将公式 (14) 简写为 : 0060 pcam Wpcax (15) 0061 对比公式 (14) 和公式 (8) , 可以看到原始的系统矩阵 W 将被压缩为 Wpca, 如果, Tsub 仅保留前 2 个主成分分量, 对于 S=36,U=256, 则通过上述压缩方法实现了 Rsub(1024) 到 Tsb(2) 的压缩, 压缩率将达到 512。相应的, 行数为 36。
28、*512*512=9437184 的原始系统矩阵 W 将被压缩为行数为 2V=18432 的压缩系统矩阵 Wpca。 0062 10、 求解压缩系统方程得到待成像物体 5 内的荧光标记物分布。 0063 本发明采用 tikhonov 正则化方法求解压缩系统方程以获得荧光标记物的三维重 建图像, tikhonov 正则化方法求解压缩系统方程为 : 0064 0065 式中, 是正则化参数, 一般选取 10-5tr(Wpca)TWpca), tr 表示矩阵的迹, T 表示矩 说 明 书 CN 102871646 A 7 6/6 页 8 阵转置。 0066 上述公式 (16) 优化问题的显式解为 : 0067 x (Wpca)TWpca+2I-1(Wpca)Tpcam (17) 0068 式中, x为求解得到的荧光标记物分布, I 是单位矩阵。 0069 上述各实施例仅用于说明本发明, 其中方法的实施步骤等都是可以有所变化的, 凡是在本发明技术方案的基础上进行的等同变换和改进, 均不应排除在本发明的保护范围 之外。 说 明 书 CN 102871646 A 8 1/2 页 9 图 1 说 明 书 附 图 CN 102871646 A 9 2/2 页 10 图 2 说 明 书 附 图 CN 102871646 A 10 。