一种血管内光声二维图像的重建方法
一种血管内光声二维图像的重建方法
【技术领域】
[0001] 本发明涉及一种利用血管内光声（intravascularphotoacoustic,IVPA)信号反 演重建出血管的轴向截面光吸收分布图像的方法，属于医学成像技术领域。
【背景技术】
[0002] 光声成像是一种新型的无损医学成像技术，它根据不同生物组织对特定波长的激 光具有差别较大的光学吸收系数，进而辐射出不同强度的超声波这一原理，通过对采集到 的组织光声信号进行图像重建，得到组织内部结构信息的一种成像方法。将光声成像技术 应用于血管内部成像，可检测血管壁的组织结构。血管内光声成像技术的原理是将脉冲激 光或幅度调制的连续激光均匀照射到血管管壁上，激发管壁产生光声信号，安置在探头顶 端的超声探测器进行360°圆周旋转采集各个方向的超声回波并传送至计算机，最后采用 适当的算法重建出血管的横截面光吸收分布或者初始光声压强分布图像，获得血管的组织 结构图像和几何形态信息，即图像重建是以这些原始信息为基础重新构筑血管内光声图 像，与其它扫描设备围绕成像物体旋转的光声断层成像系统不同，IVPA是在封闭的血管腔 内成像，其扫描孔径在成像物体内是封闭的，使数据的采集方式受到高度限制，因而无法直 接采用现有的光声图像重建算法。因此，如何进行血管内光声二维图像的重建就成为有关 技术人员面临的课题。
【发明内容】
[0003] 本发明的目的在于针对现有技术之弊端，提供一种血管内光声二维图像的重建方 法，以显示血管内部的组织结构变化，为血管疾病的临床诊断提供有益的参考。
[0004] 本发明所述问题是以下述技术方案实现的：
[0005] -种血管内光声二维图像的重建方法，所述方法首先根据IVPA成像设备采集光 声信号的方式建立初始重建图像；然后对成像导管在各个测量位置采集到的光声数据进行 一次求导、逆卷积以及滤波预处理，去除超声探测器脉冲响应对图像重建的影响；最后根据 超声探测器的有效接收角度和预处理后的光声数据计算每一个网格点对应的灰度值，得到 最终的灰度图像，所述方法包括以下步骤：
[0006] a.建立初始重建图像：
[0007] 接收光声信号的超声探测器位于成像导管顶端，成像平面通过超声探测器并垂直 于成像导管，超声探测器的扫描轨迹为平行于成像平面、半径等于导管半径的圆形轨迹； [0008] 初始重建图像A的宽度和高度均为1，A由MXM个正方形网格组成，第m行、第n 列的网格点编号记为（m，n)，其中m= 1，2，...，M，n= 1，2，...，M;M为奇数，用MXM的矩 阵B表示A，其中矩阵B的第m行、第n列元素的值B(m，n)是图像A中网格点（m，n)的灰度 值，将B初始化为零矩阵，成像导管半径为屯，成像平面所在的坐标系为X-Y平面直角坐标 系，其中坐标原点0是成像导管中心，水平向右的方向为X轴正方向，垂直于X轴向上的方 向为Y轴正方向；
[0009] b.光声信号的预处理
[0010] 成像导管在成像平面内旋转360°共采集K个位置的光声信号，步进角度为（i)= 2 /K，成像设备的采样频率为f，超声探测器在每个测量位置采集到的离散信号长度为T， 将探测器在第i个测量位置实际接收到的信号用1XT的向量^表示，对同时进行求 导、逆卷积和滤波处理，得到预处理后的结果i?:
[0012] 其中，的是P的傅里叶变换；H(?)是探测器的脉冲响应函数h(t)的傅里叶 变换；W(?)为滤波器的窗函数；IFFT( ?)表示快速傅里叶反变换；
[0013] c.重建血管横截面的初始光声压强分布图像
[0014] 对于图像A中到坐标原点的距离大于士且满足外0/2的网格点，计算探测器 在第i个测量位置探测到的网格点（m，n)产生的光声压强值DiOn，!!）：
[0015] DiOn，n) =wi(m，n) ?g, (N)
[0016] 其中，士为成像导管半径；0为超声探测器的有效接收角度；gi"（N)为中第
|个元素的值，N= 1，2,. ..，T;round( ?)表示对输入的数值进行四舍五 入取整运算，
[0018]dx= 1/M是离散网格点之间的距离，
为权重因子，
[0022] 对于第i个测量位置，对图像A的MXM个网格点都进行以上计算，得到MXM的矩 阵Di的所有元素值，然后，对第i+1个测量位置的采样数据重复上述操作，直至完成所有测 量位置的计算，再计算网格点（m，n)产生的初始压强值D(m，n):
[0024] 得到MXM的矩阵D的第m(m= 1，2，...，M)行n(n= 1，2，...，M)列元素的值，最 后，对矩阵D按照下式进行归一化处理：
[0026] 将矩阵D转换为 256 级灰度矩阵B，其中，m=l，2，...，M，n=l，2，...，M;min(D) 和max(D)分别为矩阵D所有元素中的最小值和最大值；B(m，n)为MXM的矩阵B的第m行、 第n列元素的值，显示即可得到重建图像A。
[0027] 上述血管内光声二维图像的重建方法，所述滤波器的窗函数W(?)选取Hanning 窗函数，其定义为
[0029] 其中，w。是截止频率。
[0030] 本发明可根据实际光声信号采集设备的不同，灵活地设置成像导管半径和采样频 率等设备参数，调整重建图像的分辨率，选取适当的滤波函数对原始光声信号进行预处理， 最终精确重构出目标血管横截面的初始光声压强分布图像，为血管疾病的临床诊断提供有 益的参考。
【附图说明】
[0031] 下面结合附图对本发明作进一步详述。
[0032] 图1本发明方法的流程图。
[0033] 图2IVPA系统成像及图像重建示意图。其中图（a)是IVPA系统的成像示意图，图 (b)是IVPA图像重建示意图。
[0034] 图3超声探测器接收血管壁光声信号示意图。
[0035] 图4本发明方法的重建示意图。
[0036] 文中各符号为：A、IVPA重建图像；1、重建图像的宽度/高度；M、重建图像的宽度 /高度对应的离散网格个数；(m，n)、A中第m行第n列的网格点，其中m= 1，2,. . .，M，n= 1，2,. ..，M;B、表示图像A的矩阵；B(m，n)、矩阵B的第m行第n列元素的值，即第（m，n)个 网格点（m，n)的灰度值；屯、成像导管半径；X、Y、0、X-Y平面直角坐标系的横轴、纵轴和坐 标原点；K、成像导管在成像平面内旋转360°共采集光声信号的位置数；巾、步进角度；f、 光声成像设备的采样频率；T、超声探测器在每个测量位置采集到的离散信号长度；^、第i个测量位置采集到的去除探测器脉冲响应影响后的离散光声信号，其中i= 1，2，...，K; 、探测器在第i个测量位置实际接收到的离散光声信号；h(t)、探测器的脉冲响应函数； 、H(?)、$和h(t)的傅里叶变换；IFFT( ?)、快速傅里叶反变换；0、超声探测器的 有效接收角度；c、光声信号在血管内的传播速度；t、探测到光声信号的时刻；r、光声信号 所在弧的半径；^'、预处理后的结果；W(?)、滤波器的窗函数；《。、截止频率；$、探 测器在第i个采样点时接收平面的法向量；iVn)、第（m，n)个网格点到坐标原点（成像导管 中心）的距离；dx、离散网格点之间的距离；^、探测器的第i个测量位置点到第（m，n) 个网格点的向量；g、探测器在第i个测量位置接收平面的法向量；外、向量g与^ 之间的夹角；?、是向量的点积运算
的模；cos、?）、反余弦函 数；riGll，n)、第（m，n)个网格点所在弧线的半径；gi"（N)、P中的第N个元素的值，其中N= 1，2,. . .，T;round( ?)、对数值进行四舍五入的取整运算；wi0n，n)、权重因子；DpMXM的矩阵；Di(m，n)、矩阵Di的第m行n列元素的值，是探测器在第i个测量位置探测到的网格点（m，n) 产生的光声压强值；D、MXM的矩阵；D(m，n)、矩阵D的第m行第n列元素的值，是最终计算 得到的网格点（m，n)产生的初始压强值；min(D)、max(D)、矩阵D所有元素中的最小值和最 大值。
【具体实施方式】
[0037] 下面结合附图1详细说明本发明的步骤：
[0038] (1)初始化：
[0039] 如附图2(a)所示，接收光声信号的超声探测器位于成像导管顶端，成像平面通过 超声探测器并垂直于成像导管。本发明方法将超声探测器看作理想的点探测器，其扫描轨 迹为平行于成像平面、半径等于导管半径的圆形轨迹。
[0040] 首先，如附图2(b)所示，初始重建图像A的宽度和高度均为1，A由MXM个正方形 网格组成，第m行、第n列的网格点编号记为（m,n)，其中m= 1，2,. . .，M，n= 1，2,. ..，M; M为奇数，并且M越大，重建图像的分辨率越高。用MXM的矩阵B表示A，其中矩阵B的第m 行、第n列元素的值B(m，n)是图像A中网格点（m，n)的灰度值，将B初 始化为零矩阵。成 像导管半径为屯，成像平面所在的坐标系为X-Y平面直角坐标系，其中坐标原点0是成像导 管中心，水平向右的方向为X轴正方向，垂直于X轴向上的方向为Y轴正方向。
[0041] (2)、光声信号的预处理
[0042] 成像导管在成像平面内旋转360°共采集K个位置的光声信号，步进角度为小= 2 /K，成像设备的采样频率为f，超声探测器在每个测量位置采集到的离散信号长度为T。 将超声探测器在第i(i= 1，2，...，K)个测量位置采集到的光声信号用1XT的向量i表 不〇
[0043] 实际上，由于探测器的带宽有限，脉冲响应函数不是理想的S函数，而是有一定 的频率响应范围，所以探测器在第i个测量位置实际接收到的信号为
[0045] 其中，h(t)为探测器的脉冲响应函数；表示卷积运算。通过逆卷积运算，可由 艮'求得g,:
[0047] 其中，斤(0)和H(?)分别是f和h(t)的傅里叶变换；IFFT( ?)表示快速傅里叶 反变换。
[0048] 附图3为附图2(a)中的IVPA成像平面，超声探测器的有效接收角度为0，光声信 号在血管内的传播速度为c。非聚焦超声探测器在第i(i= 1，2, ...，K)个测量位置、时刻 t=r/c探测到的光声信号是以当前探测器位置为圆心、半径为r、圆心角为0的弧线上所 有吸收点发出的光声压强信号的积分。由于积分运算具有误差叠加效应，因此本发明方法 将每个测量位置上的对时间取一阶导数，然后再进行加权累加，得到近似的初始光声压 强分布。
[0049] 为了保证算法的稳定性，本发明方法对同时进行求导、逆卷积和滤波处理，得 到f'：
[0051] 其中，W(?)为滤波器的窗函数，本发明方法选取Hanning窗函数，其定义为
[0053] 其中，是截止频率。
[0054] (3)重建血管横截面的初始光声压强分布图像
[0055] 图像A中的网格点（m，n)到坐标原点（即成像导管中心）的距离为
[0057] 其中，dx= 1/M是离散网格点之间的距离。当rL时，说明网格点（m，n)位 于成像导管内，对该类网格点不做计算。反之，当网格点在超声成像导管外部，即iVn)>cU 时，对这些网格点继续进行以下步骤。
[0058] 如附图4所示，当探测器位于第i个测量位置时，其位置坐标为 (d0cos(i<i))，d0sin(i<i)))，该点到网格点（m，n)的向量为：
[0060] 接收平面的法向量为
[0062] $与^之间的夹角为
[0063]
[0064] 其中，"是向量的点积运算；7^和办m,n)分别是％和g/(mn|的模；cos1 ( ?)为反 余弦函数。若外说明网格点（m，n)位于以第i个测量位置为中心、圆心角为0 的扇形有效测量范围内，并且该网格点位于半径为：
[0066] 的弧线上。该弧线上所有网格点产生的光声压强的累加等于中第
个元素的值，即g/'（N)。其中，N= 1，2,. ..，T;round( ?)表示对输入 的数值进行四舍五入取整运算。由gi"（N)根据下式得到
[0067] DiOn，n) =wi(m，n) ?g/'（N) (10)
[0068] 其中，wi0ll，n)为权重因子；D,(m，n)为MXM的矩阵Di的第m行n列元素的值，是探 测器在第i个测量位置探测到的网格点（m，n)被脉冲激光照射后产生的光声压强值。当测 量位置比较密集时，相邻测量位置之间的扇形有效测量范围会发生重叠，为了避免重复计 算重叠区域的网格点在不同测量位置产生的光声压强，本发明方法引入权重因子wi0n，n)
[0070] 由式（11)可知，外越小，权重wi0ll，n)越大，即g/'(iV)对计算网格点（m,n)处的 初始光声压强的影响越大。
[0071] 对于第i个测量位置，对图像A的MXM个网格点都进行以上计算，得到矩阵 所有元素值。此后，对第i+1个测量位置的采样数据重复上述操作，直至完成所有测量位置 的计算。那么
[0073]其中，D(m,n)为MXM的矩阵D的第m(m=1，2, ? ? ?，M)行n(n=1，2, ? ? ?，M)列 元素的值，矩阵中的各元素值对应为最终计算得到的各网格点的初始压强值。
[0074] 最后，将D按照式（13)进行归一化处理转换为256级灰度矩阵B，显示即可得到重 建图像A:
[0076]其中，m=l，2，...，M，n=l，2，...，M;min〇))和max(D)分别为矩阵D所有元素 中的最小值和最大值；B(m，n)为MXM的矩阵B的第m行、第n列元素的值。
【主权项】
1. 一种血管内光声二维图像的重建方法，其特征是，所述方法首先根据IVPA成像设备 采集光声信号的方式建立初始重建图像；然后对成像导管在各个测量位置采集到的光声数 据进行一次求导、逆卷积以及滤波预处理，去除超声探测器脉冲响应对图像重建的影响；最 后根据超声探测器的有效接收角度和预处理后的光声数据计算每一个网格点对应的灰度 值，得到最终的灰度图像，所述方法包括以下步骤： a. 建立初始重建图像： 接收光声信号的超声探测器位于成像导管顶端，成像平面通过超声探测器并垂直于成 像导管，超声探测器的扫描轨迹为平行于成像平面、半径等于导管半径的圆形轨迹； 初始重建图像A的宽度和高度均为1，A由MXM个正方形网格组成，第m行、第η列的 网格点编号记为（m，η)，其中m = 1，2,. . .，Μ，η = 1，2,. . .，M ;Μ为奇数，用MXM的矩阵B 表示Α，其中矩阵B的第m行、第η列元素的值B(m，η)是图像A中网格点（m，η)的灰度值， 将B初始化为零矩阵，成像导管半径为Cltl，成像平面所在的坐标系为X-Y平面直角坐标系， 其中坐标原点O是成像导管中心，水平向右的方向为X轴正方向，垂直于X轴向上的方向为 Y轴正方向； b. 光声信号的预处理 成像导管在成像平面内旋转360°共采集K个位置的光声信号，步进角度为Φ =2π/ Κ，成像设备的采样频率为f，超声探测器在每个测量位置采集到的离散信号长度为Τ，将探 测器在第i个测量位置实际接收到的信号用IXT的向量^表示，对&同时进行求导、逆 卷积和滤波处理，得到预处理后的结果^ :其中，€(的是的傅里叶变换；Η(ω)是探测器的脉冲响应函数h(t)的傅里叶变 换；W(Q)为滤波器的窗函数；IFFT( ·)表示快速傅里叶反变换； c. 重建血管横截面的初始光声压强分布图像 对于图像A中到坐标原点的距离大于Cltl且满足f 0/2的网格点，计算探测器在第 i个测量位置探测到的网格点（m，η)产生的光声压强值Di (m，n): Di (m, n) = wi(m；n) · (N) 其中，Cltl为成像导管半径；Θ为超声探测器的有效接收角度；gi"（N)为中第个元素的值，N = 1，2,. . .，T ;round( ·)表示对输入的数值进行四舍五 入取整运算，dx= 1/M是离散网格点之间的距离，为权重因子，对于第i个测量位置，对图像A的MXM个网格点都进行以上计算，得到MXM的矩阵Di的所有元素值，然后，对第i+Ι个测量位置的采样数据重复上述操作，直至完成所有测量位 置的计算，再计算网格点（m，η)产生的初始压强值D (m，η):得到MXM的矩阵D的第m(m = 1，2,. . .，Μ)行η(η = 1，2,. . .，Μ)列元素的值，最后， 对矩阵D按照下式进行归一化处理：将矩阵D转换为256级灰度矩阵Β，其中，m= 1，2，...，M，n = l，2，...，M;min(D)和 max (D)分别为矩阵D所有元素中的最小值和最大值；B (m, η)为MXM的矩阵B的第m行、第 η列元素的值，显示即可得到重建图像A。2.根据权利要求1所述的血管内光声二维图像的重建方法，其特征是，所述滤波器的 窗函数W(Q)选取Hanning窗函数，其定义为其中，ω。是截止频率。
【专利摘要】一种血管内光声二维图像的重建方法，所述方法首先根据IVPA成像设备采集光声信号的方式建立初始重建图像；然后对成像导管在各个测量位置采集到的光声数据进行一次求导、逆卷积以及滤波预处理，去除超声探测器脉冲响应对图像重建的影响；最后根据超声探测器的有效接收角度和预处理后的光声数据计算每一个网格点对应的灰度值，得到最终的灰度图像。本发明可根据实际光声信号采集设备的不同，灵活地设置成像导管半径和采样频率等设备参数，调整重建图像的分辨率，选取适当的滤波函数对原始光声信号进行预处理，最终精确重构出目标血管横截面的初始光声压强分布图像，为血管疾病的临床诊断提供有益的参考。
【IPC分类】G06T11/00
【公开号】CN104899902
【申请号】CN201510177981
【发明人】孙正, 韩朵朵
【申请人】华北电力大学（保定）
【公开日】2015年9月9日
【申请日】2015年4月14日
【技术领域】
[0001] 本发明涉及一种利用血管内光声（intravascularphotoacoustic,IVPA)信号反 演重建出血管的轴向截面光吸收分布图像的方法，属于医学成像技术领域。
【背景技术】
[0002] 光声成像是一种新型的无损医学成像技术，它根据不同生物组织对特定波长的激 光具有差别较大的光学吸收系数，进而辐射出不同强度的超声波这一原理，通过对采集到 的组织光声信号进行图像重建，得到组织内部结构信息的一种成像方法。将光声成像技术 应用于血管内部成像，可检测血管壁的组织结构。血管内光声成像技术的原理是将脉冲激 光或幅度调制的连续激光均匀照射到血管管壁上，激发管壁产生光声信号，安置在探头顶 端的超声探测器进行360°圆周旋转采集各个方向的超声回波并传送至计算机，最后采用 适当的算法重建出血管的横截面光吸收分布或者初始光声压强分布图像，获得血管的组织 结构图像和几何形态信息，即图像重建是以这些原始信息为基础重新构筑血管内光声图 像，与其它扫描设备围绕成像物体旋转的光声断层成像系统不同，IVPA是在封闭的血管腔 内成像，其扫描孔径在成像物体内是封闭的，使数据的采集方式受到高度限制，因而无法直 接采用现有的光声图像重建算法。因此，如何进行血管内光声二维图像的重建就成为有关 技术人员面临的课题。
【发明内容】
[0003] 本发明的目的在于针对现有技术之弊端，提供一种血管内光声二维图像的重建方 法，以显示血管内部的组织结构变化，为血管疾病的临床诊断提供有益的参考。
[0004] 本发明所述问题是以下述技术方案实现的：
[0005] -种血管内光声二维图像的重建方法，所述方法首先根据IVPA成像设备采集光 声信号的方式建立初始重建图像；然后对成像导管在各个测量位置采集到的光声数据进行 一次求导、逆卷积以及滤波预处理，去除超声探测器脉冲响应对图像重建的影响；最后根据 超声探测器的有效接收角度和预处理后的光声数据计算每一个网格点对应的灰度值，得到 最终的灰度图像，所述方法包括以下步骤：
[0006] a.建立初始重建图像：
[0007] 接收光声信号的超声探测器位于成像导管顶端，成像平面通过超声探测器并垂直 于成像导管，超声探测器的扫描轨迹为平行于成像平面、半径等于导管半径的圆形轨迹； [0008] 初始重建图像A的宽度和高度均为1，A由MXM个正方形网格组成，第m行、第n 列的网格点编号记为（m，n)，其中m= 1，2，...，M，n= 1，2，...，M;M为奇数，用MXM的矩 阵B表示A，其中矩阵B的第m行、第n列元素的值B(m，n)是图像A中网格点（m，n)的灰度 值，将B初始化为零矩阵，成像导管半径为屯，成像平面所在的坐标系为X-Y平面直角坐标 系，其中坐标原点0是成像导管中心，水平向右的方向为X轴正方向，垂直于X轴向上的方 向为Y轴正方向；
[0009] b.光声信号的预处理
[0010] 成像导管在成像平面内旋转360°共采集K个位置的光声信号，步进角度为（i)= 2 /K，成像设备的采样频率为f，超声探测器在每个测量位置采集到的离散信号长度为T， 将探测器在第i个测量位置实际接收到的信号用1XT的向量^表示，对同时进行求 导、逆卷积和滤波处理，得到预处理后的结果i?:
[0012] 其中，的是P的傅里叶变换；H(?)是探测器的脉冲响应函数h(t)的傅里叶 变换；W(?)为滤波器的窗函数；IFFT( ?)表示快速傅里叶反变换；
[0013] c.重建血管横截面的初始光声压强分布图像
[0014] 对于图像A中到坐标原点的距离大于士且满足外0/2的网格点，计算探测器 在第i个测量位置探测到的网格点（m，n)产生的光声压强值DiOn，!!）：
[0015] DiOn，n) =wi(m，n) ?g, (N)
[0016] 其中，士为成像导管半径；0为超声探测器的有效接收角度；gi"（N)为中第
|个元素的值，N= 1，2,. ..，T;round( ?)表示对输入的数值进行四舍五 入取整运算，
[0018]dx= 1/M是离散网格点之间的距离，
为权重因子，
[0022] 对于第i个测量位置，对图像A的MXM个网格点都进行以上计算，得到MXM的矩 阵Di的所有元素值，然后，对第i+1个测量位置的采样数据重复上述操作，直至完成所有测 量位置的计算，再计算网格点（m，n)产生的初始压强值D(m，n):
[0024] 得到MXM的矩阵D的第m(m= 1，2，...，M)行n(n= 1，2，...，M)列元素的值，最 后，对矩阵D按照下式进行归一化处理：
[0026] 将矩阵D转换为 256 级灰度矩阵B，其中，m=l，2，...，M，n=l，2，...，M;min(D) 和max(D)分别为矩阵D所有元素中的最小值和最大值；B(m，n)为MXM的矩阵B的第m行、 第n列元素的值，显示即可得到重建图像A。
[0027] 上述血管内光声二维图像的重建方法，所述滤波器的窗函数W(?)选取Hanning 窗函数，其定义为
[0029] 其中，w。是截止频率。
[0030] 本发明可根据实际光声信号采集设备的不同，灵活地设置成像导管半径和采样频 率等设备参数，调整重建图像的分辨率，选取适当的滤波函数对原始光声信号进行预处理， 最终精确重构出目标血管横截面的初始光声压强分布图像，为血管疾病的临床诊断提供有 益的参考。
【附图说明】
[0031] 下面结合附图对本发明作进一步详述。
[0032] 图1本发明方法的流程图。
[0033] 图2IVPA系统成像及图像重建示意图。其中图（a)是IVPA系统的成像示意图，图 (b)是IVPA图像重建示意图。
[0034] 图3超声探测器接收血管壁光声信号示意图。
[0035] 图4本发明方法的重建示意图。
[0036] 文中各符号为：A、IVPA重建图像；1、重建图像的宽度/高度；M、重建图像的宽度 /高度对应的离散网格个数；(m，n)、A中第m行第n列的网格点，其中m= 1，2,. . .，M，n= 1，2,. ..，M;B、表示图像A的矩阵；B(m，n)、矩阵B的第m行第n列元素的值，即第（m，n)个 网格点（m，n)的灰度值；屯、成像导管半径；X、Y、0、X-Y平面直角坐标系的横轴、纵轴和坐 标原点；K、成像导管在成像平面内旋转360°共采集光声信号的位置数；巾、步进角度；f、 光声成像设备的采样频率；T、超声探测器在每个测量位置采集到的离散信号长度；^、第i个测量位置采集到的去除探测器脉冲响应影响后的离散光声信号，其中i= 1，2，...，K; 、探测器在第i个测量位置实际接收到的离散光声信号；h(t)、探测器的脉冲响应函数； 、H(?)、$和h(t)的傅里叶变换；IFFT( ?)、快速傅里叶反变换；0、超声探测器的 有效接收角度；c、光声信号在血管内的传播速度；t、探测到光声信号的时刻；r、光声信号 所在弧的半径；^'、预处理后的结果；W(?)、滤波器的窗函数；《。、截止频率；$、探 测器在第i个采样点时接收平面的法向量；iVn)、第（m，n)个网格点到坐标原点（成像导管 中心）的距离；dx、离散网格点之间的距离；^、探测器的第i个测量位置点到第（m，n) 个网格点的向量；g、探测器在第i个测量位置接收平面的法向量；外、向量g与^ 之间的夹角；?、是向量的点积运算
的模；cos、?）、反余弦函 数；riGll，n)、第（m，n)个网格点所在弧线的半径；gi"（N)、P中的第N个元素的值，其中N= 1，2,. . .，T;round( ?)、对数值进行四舍五入的取整运算；wi0n，n)、权重因子；DpMXM的矩阵；Di(m，n)、矩阵Di的第m行n列元素的值，是探测器在第i个测量位置探测到的网格点（m，n) 产生的光声压强值；D、MXM的矩阵；D(m，n)、矩阵D的第m行第n列元素的值，是最终计算 得到的网格点（m，n)产生的初始压强值；min(D)、max(D)、矩阵D所有元素中的最小值和最 大值。
【具体实施方式】
[0037] 下面结合附图1详细说明本发明的步骤：
[0038] (1)初始化：
[0039] 如附图2(a)所示，接收光声信号的超声探测器位于成像导管顶端，成像平面通过 超声探测器并垂直于成像导管。本发明方法将超声探测器看作理想的点探测器，其扫描轨 迹为平行于成像平面、半径等于导管半径的圆形轨迹。
[0040] 首先，如附图2(b)所示，初始重建图像A的宽度和高度均为1，A由MXM个正方形 网格组成，第m行、第n列的网格点编号记为（m,n)，其中m= 1，2,. . .，M，n= 1，2,. ..，M; M为奇数，并且M越大，重建图像的分辨率越高。用MXM的矩阵B表示A，其中矩阵B的第m 行、第n列元素的值B(m，n)是图像A中网格点（m，n)的灰度值，将B初 始化为零矩阵。成 像导管半径为屯，成像平面所在的坐标系为X-Y平面直角坐标系，其中坐标原点0是成像导 管中心，水平向右的方向为X轴正方向，垂直于X轴向上的方向为Y轴正方向。
[0041] (2)、光声信号的预处理
[0042] 成像导管在成像平面内旋转360°共采集K个位置的光声信号，步进角度为小= 2 /K，成像设备的采样频率为f，超声探测器在每个测量位置采集到的离散信号长度为T。 将超声探测器在第i(i= 1，2，...，K)个测量位置采集到的光声信号用1XT的向量i表 不〇
[0043] 实际上，由于探测器的带宽有限，脉冲响应函数不是理想的S函数，而是有一定 的频率响应范围，所以探测器在第i个测量位置实际接收到的信号为
[0045] 其中，h(t)为探测器的脉冲响应函数；表示卷积运算。通过逆卷积运算，可由 艮'求得g,:
[0047] 其中，斤(0)和H(?)分别是f和h(t)的傅里叶变换；IFFT( ?)表示快速傅里叶 反变换。
[0048] 附图3为附图2(a)中的IVPA成像平面，超声探测器的有效接收角度为0，光声信 号在血管内的传播速度为c。非聚焦超声探测器在第i(i= 1，2, ...，K)个测量位置、时刻 t=r/c探测到的光声信号是以当前探测器位置为圆心、半径为r、圆心角为0的弧线上所 有吸收点发出的光声压强信号的积分。由于积分运算具有误差叠加效应，因此本发明方法 将每个测量位置上的对时间取一阶导数，然后再进行加权累加，得到近似的初始光声压 强分布。
[0049] 为了保证算法的稳定性，本发明方法对同时进行求导、逆卷积和滤波处理，得 到f'：
[0051] 其中，W(?)为滤波器的窗函数，本发明方法选取Hanning窗函数，其定义为
[0053] 其中，是截止频率。
[0054] (3)重建血管横截面的初始光声压强分布图像
[0055] 图像A中的网格点（m，n)到坐标原点（即成像导管中心）的距离为
[0057] 其中，dx= 1/M是离散网格点之间的距离。当rL时，说明网格点（m，n)位 于成像导管内，对该类网格点不做计算。反之，当网格点在超声成像导管外部，即iVn)>cU 时，对这些网格点继续进行以下步骤。
[0058] 如附图4所示，当探测器位于第i个测量位置时，其位置坐标为 (d0cos(i<i))，d0sin(i<i)))，该点到网格点（m，n)的向量为：
[0060] 接收平面的法向量为
[0062] $与^之间的夹角为
[0063]
[0064] 其中，"是向量的点积运算；7^和办m,n)分别是％和g/(mn|的模；cos1 ( ?)为反 余弦函数。若外说明网格点（m，n)位于以第i个测量位置为中心、圆心角为0 的扇形有效测量范围内，并且该网格点位于半径为：
[0066] 的弧线上。该弧线上所有网格点产生的光声压强的累加等于中第
个元素的值，即g/'（N)。其中，N= 1，2,. ..，T;round( ?)表示对输入 的数值进行四舍五入取整运算。由gi"（N)根据下式得到
[0067] DiOn，n) =wi(m，n) ?g/'（N) (10)
[0068] 其中，wi0ll，n)为权重因子；D,(m，n)为MXM的矩阵Di的第m行n列元素的值，是探 测器在第i个测量位置探测到的网格点（m，n)被脉冲激光照射后产生的光声压强值。当测 量位置比较密集时，相邻测量位置之间的扇形有效测量范围会发生重叠，为了避免重复计 算重叠区域的网格点在不同测量位置产生的光声压强，本发明方法引入权重因子wi0n，n)
[0070] 由式（11)可知，外越小，权重wi0ll，n)越大，即g/'(iV)对计算网格点（m,n)处的 初始光声压强的影响越大。
[0071] 对于第i个测量位置，对图像A的MXM个网格点都进行以上计算，得到矩阵 所有元素值。此后，对第i+1个测量位置的采样数据重复上述操作，直至完成所有测量位置 的计算。那么
[0073]其中，D(m,n)为MXM的矩阵D的第m(m=1，2, ? ? ?，M)行n(n=1，2, ? ? ?，M)列 元素的值，矩阵中的各元素值对应为最终计算得到的各网格点的初始压强值。
[0074] 最后，将D按照式（13)进行归一化处理转换为256级灰度矩阵B，显示即可得到重 建图像A:
[0076]其中，m=l，2，...，M，n=l，2，...，M;min〇))和max(D)分别为矩阵D所有元素 中的最小值和最大值；B(m，n)为MXM的矩阵B的第m行、第n列元素的值。
【主权项】
1. 一种血管内光声二维图像的重建方法，其特征是，所述方法首先根据IVPA成像设备 采集光声信号的方式建立初始重建图像；然后对成像导管在各个测量位置采集到的光声数 据进行一次求导、逆卷积以及滤波预处理，去除超声探测器脉冲响应对图像重建的影响；最 后根据超声探测器的有效接收角度和预处理后的光声数据计算每一个网格点对应的灰度 值，得到最终的灰度图像，所述方法包括以下步骤： a. 建立初始重建图像： 接收光声信号的超声探测器位于成像导管顶端，成像平面通过超声探测器并垂直于成 像导管，超声探测器的扫描轨迹为平行于成像平面、半径等于导管半径的圆形轨迹； 初始重建图像A的宽度和高度均为1，A由MXM个正方形网格组成，第m行、第η列的 网格点编号记为（m，η)，其中m = 1，2,. . .，Μ，η = 1，2,. . .，M ;Μ为奇数，用MXM的矩阵B 表示Α，其中矩阵B的第m行、第η列元素的值B(m，η)是图像A中网格点（m，η)的灰度值， 将B初始化为零矩阵，成像导管半径为Cltl，成像平面所在的坐标系为X-Y平面直角坐标系， 其中坐标原点O是成像导管中心，水平向右的方向为X轴正方向，垂直于X轴向上的方向为 Y轴正方向； b. 光声信号的预处理 成像导管在成像平面内旋转360°共采集K个位置的光声信号，步进角度为Φ =2π/ Κ，成像设备的采样频率为f，超声探测器在每个测量位置采集到的离散信号长度为Τ，将探 测器在第i个测量位置实际接收到的信号用IXT的向量^表示，对&同时进行求导、逆 卷积和滤波处理，得到预处理后的结果^ :其中，€(的是的傅里叶变换；Η(ω)是探测器的脉冲响应函数h(t)的傅里叶变 换；W(Q)为滤波器的窗函数；IFFT( ·)表示快速傅里叶反变换； c. 重建血管横截面的初始光声压强分布图像 对于图像A中到坐标原点的距离大于Cltl且满足f 0/2的网格点，计算探测器在第 i个测量位置探测到的网格点（m，η)产生的光声压强值Di (m，n): Di (m, n) = wi(m；n) · (N) 其中，Cltl为成像导管半径；Θ为超声探测器的有效接收角度；gi"（N)为中第个元素的值，N = 1，2,. . .，T ;round( ·)表示对输入的数值进行四舍五 入取整运算，dx= 1/M是离散网格点之间的距离，为权重因子，对于第i个测量位置，对图像A的MXM个网格点都进行以上计算，得到MXM的矩阵Di的所有元素值，然后，对第i+Ι个测量位置的采样数据重复上述操作，直至完成所有测量位 置的计算，再计算网格点（m，η)产生的初始压强值D (m，η):得到MXM的矩阵D的第m(m = 1，2,. . .，Μ)行η(η = 1，2,. . .，Μ)列元素的值，最后， 对矩阵D按照下式进行归一化处理：将矩阵D转换为256级灰度矩阵Β，其中，m= 1，2，...，M，n = l，2，...，M;min(D)和 max (D)分别为矩阵D所有元素中的最小值和最大值；B (m, η)为MXM的矩阵B的第m行、第 η列元素的值，显示即可得到重建图像A。2.根据权利要求1所述的血管内光声二维图像的重建方法，其特征是，所述滤波器的 窗函数W(Q)选取Hanning窗函数，其定义为其中，ω。是截止频率。
【专利摘要】一种血管内光声二维图像的重建方法，所述方法首先根据IVPA成像设备采集光声信号的方式建立初始重建图像；然后对成像导管在各个测量位置采集到的光声数据进行一次求导、逆卷积以及滤波预处理，去除超声探测器脉冲响应对图像重建的影响；最后根据超声探测器的有效接收角度和预处理后的光声数据计算每一个网格点对应的灰度值，得到最终的灰度图像。本发明可根据实际光声信号采集设备的不同，灵活地设置成像导管半径和采样频率等设备参数，调整重建图像的分辨率，选取适当的滤波函数对原始光声信号进行预处理，最终精确重构出目标血管横截面的初始光声压强分布图像，为血管疾病的临床诊断提供有益的参考。
【IPC分类】G06T11/00
【公开号】CN104899902
【申请号】CN201510177981
【发明人】孙正, 韩朵朵
【申请人】华北电力大学（保定）
【公开日】2015年9月9日
【申请日】2015年4月14日
最新回复(0)