CN106108942B

CN106108942B - 基于OpenCL的并行ARFI成像方法

Info

Publication number: CN106108942B
Application number: CN201610450849.8A
Authority: CN
Inventors: 尹皓; 石丹; 刘东权
Original assignee: Sichuan University
Current assignee: Sichuan University
Priority date: 2016-06-21
Filing date: 2016-06-21
Publication date: 2018-11-30
Anticipated expiration: 2036-06-21
Also published as: CN106108942A

Abstract

本发明涉及超声成像领域，特别涉及一种基于OpenCL的并行ARFI成像方法。本发明提供的基于OpenCL语言，可基于各种显卡GPU进行计算的ARFI成像方法，通过将传统的ARFI成像步骤中的位移估算、中值滤波、扫描转换等步骤都应用OpenCL理念将其分解为指定项数的工作组及工作项，通过每个工作项对应一个或多个待处理数据的方式，并行处理ARFI成像过程中的I、Q数据，从而大幅度的提高了数据的处理速度，提高了ARFI成像中数据处理效率。

Description

基于OpenCL的并行ARFI成像方法

技术领域

本发明涉及超声成像领域，特别涉及一种基于OpenCL的并行ARFI成像方法。

背景技术

声辐射力脉冲(Acoustic Radiation Force Impulse,ARFI)弹性成像是过去十几年间提出并快速发展起来的一项新的弹性成像技术。ARFI弹性成像通过向生物组织发射短时程脉冲(几十～几百微秒)聚焦于组织感兴趣区域(Region Of Interest,ROI)来产生辐射力，使其产生瞬时、微米级的形变位移，位移的大小取决于组织弹性模量的大小；弹性模量大，组织的相对较硬，激励所引起的位移小；反之，弹性模量小，组织相对较软，激励所引起的位移大。然后发射超声检测脉冲序列同时在时间和空间两个维度上对这些位移进行追踪和监测，通过信号相关性算法估算出组织位移，并对其进行分析和成像。声辐射力弹性成像方法直接将辐射力聚焦在感兴趣组织位置，可以穿透浅表组织直接抵达并作用在较深部组织上，拓展了传统超声弹性成像的范围。

例如，申请号为2014106604147，名称为《基于声辐射力回波位移检测方法及成像方法》的专利提供了一种基于声辐射力的回拨位移检测方法及成像方法，其通过交替发射两种脉冲分别用于激励组织和检测组织响应。该脉冲序列由至少一个参考脉冲、一个或多个激励脉冲以及多个连续的检测脉冲组成，图1给出了典型ARFI脉冲序列的示例，其包含24个脉冲信号，其中包括4个激励脉冲，且激励脉冲与检测脉冲交替排列，激励结束后有16个连续检测脉冲，第一个出现的检测脉冲即为参考脉冲。图2给出了ARFI成像的基本流程，当获取检测脉冲的回波数据后，该专利通过自相关Loupas算法计算所述回波在时间取样窗口内的平均位移速度，进而获得回波在时间方向上的平均位移，进一步的通过中值滤波、扫描转换等步骤将其处理成像，但是在该专利中，该检测方法是基于CPU依次对各个维度上的各个回波数据来进行计算，其计算效率极其低下。

发明内容

本发明的目的在于克服现有技术中只能依靠CPU依次对各个维度上的回波数据进行计算，导致计算效率极其低下的问题，提供一种基于OpenCL语言的，可基于各种显卡GPU进行计算，且可大幅度提高计算效率的并行ARFI成像方法。

为了实现上述发明目的，本发明提供了以下技术方案：

一种基于OpenCL的并行ARFI成像方法，位移估算的步骤采用如下方法完成：

(1-1)根据扫描线条数gx及每条扫描线中的检测脉冲个数gy设置工作组为2维的(gx,gy)，每个工作组对应某一时刻的检测脉冲；每个工作组中设置(ix,iy)＝(64×k,1)个工作项，k为2以上整数。

(1-2)将获取的声辐射力回波信号I、Q数据按照深度方向顺序存储；将I、Q数据从主机端读取到设备端，如果每个检测脉冲在深度方向上的采样点个数不满足128字节的整数倍，则通过补零使其满足。

(1-3)将每个检测脉冲的采样点分配入对应的工作组中的各工作项中，由该工作组中各工作项基于自相关Loupas算法计算所述回波在每个采样点处的四个自相关参数并存储；四个自相关参数的计算公式分别是：

其中，M是深度方向的取样窗口大小，N是检测脉冲时间方向的取样窗口大小。

(1-4)同一工作组中的各个采样点的自相关参数均计算完成后，工作项基于自相关Loupas算法计算所述回波在深度方向取样窗口内的平均位移

其中，c是声速，f_c是发射信号中心频率，f_s是采样频率；

(1-5)将计算得出的回波平均位移数据按照深度方向顺序存储。

进一步的，步骤(1-3)中每个检测脉冲中的各个采样点的四个自相参数计算完成后存入共享存储器，其大小为nSample*4*sizeof(float)，其中nSample表示采样点的个数，4表示公式中包含四项自相关参数。

进一步的，每个采样点数据包括I、Q两组数据；采用数据类型float2进行存储,float2表示2个连续存放的32-bit的浮点数，分别存放I和Q数据。

进一步的，所述主机端为CPU，所述设备端为GPU；所述I、Q数据读入GPU后存入全局存储器中。

进一步的，中值滤波的步骤采用如下方法完成：

(2-1)设定工作组大小为(i，j)其中i、j为相等或不等的8的整数倍；设置个工作组，符号表示求不小于给定实数的最小整数；nPoint表示每个检测脉冲在深度方向上的平均位移个数,i表示在深度方向上的工作项个数，j表示在扫描线方向上的工作项个数；每一个工作项对应一个特定时刻的回波位移数据。

(2-2)截取任意时刻的平均位移数据，从主机端将数据读取到设备端，每一个工作项处理其对应的一个位移数据；以当前点为中心以xsize*ysize为取样窗口大小，计算当前点在邻域范围内的中位数，将每个取样窗口的中心点值用中位数替换；xsize,ysize为3以上奇数；所有取样窗口数据存储在共享存储器中，大小为i*j*xsize*ysize*sizeof(float)。

(2-3)将当前点值计算结果按照深度方向顺序存储。

进一步的，在从主机端向设备端读取数据时，如果每个扫描线该时刻深度方向上的平均位移数据不满足128字节的整数倍，则通过补零使其达到对齐访问的要求。

进一步的，采用中位数算法(MedianOfMeidans)计算每个采样组中的中位数数值；而不必对每一个点的所有邻域数据进行排序，计算时，采用宏定义并结合循环展开方式进行，采用宏定义的方式省去了函数调用过程，通过循环展开方式减少了控制指令，提高了算法效率。

进一步的，扫描转换的步骤中，采用如下方法实现：

(3-1)设定工作组大小为(x,y)，x，y为相等或不等的16的整数倍；设置个工作组，每个工作项对应一个目标像素点,ImageW表示目标图像水平像素，其对应扫描线方向，ImageH表示目标图像垂直像素，其对应深度方向。

(3-2)目标像素点坐标转换至原始数据坐标，所述原始数据坐标是指将步骤(2-3)的计算结果按照深度方向为横轴、扫描线方向为纵轴排列的坐标；根据每个目标像素点对应的一个以上的原始数据坐标对应的数据计算该目标像素点的像素值。

(3-3)将步骤(3-2)的计算结果按照目标像素点坐标的分布存储。

进一步的，步骤(3-2)中是通过双线性插值或者三次样条插值法计算得到当前点的像素值。

进一步的，步骤(3-2)中，还包括设置最大位移阈值T，将获取的像素值归一化到[0,255]的像素范围内的步骤。

进一步的，步骤(3-2)输出的数据类型为uchar或者伪彩int。

与现有技术相比，本发明的有益效果：本发明提供的基于OpenCL语言，可基于各种显卡GPU进行计算的ARFI成像方法，通过将传统的ARFI成像步骤中的位移估算、中值滤波、扫描转换等步骤都应用OpenCL理念将其分解为指定项数个工作组及工作项，通过每个工作项对应一待处理数据的方式，并行处理ARFI成像过程中的I、Q数据，从而大幅度的提高了数据的处理速度，提高了ARFI成像中数据处理效率。

同时，基于OpenCL的处理算法可以适用于多核心处理器(CPU)、图形处理器(GPU)、数字信号处理器(DSP)、以及现场可编程门阵列(FPGA)等其他并行处理器中，从而使得本方法可以在基于上述处理器的硬件中实行。

附图说明：

图1现有技术中典型的脉冲序列。

图2是现有技术中典型的ARFI成像流程。

图3是为本发明提供的基于OpenCL的并行ARFI成像方法流程图。

图4是本发明中，位移估算步骤中每个工作项需要计算的采样点示意图。

具体实施方式

下面结合附图及具体实施例对本发明作进一步的详细描述。但不应将此理解为本发明上述主题的范围仅限于以下的实施例，凡基于本发明内容所实现的技术均属于本发明的范围。

实施例1：本实施例提供一种基于OpenCL的并行ARFI成像方法，位移估算的步骤采用如下方法完成：

S110：根据扫描线条数gx及每条扫描线中的检测脉冲个数gy设置工作组为2维的(gx,gy)，每个工作组对应某一时刻的检测脉冲；每个工作组中设置(ix,iy)＝(64×k,1)个工作项，k为2以上整数；本实施例中，每个工作组中设置256个工作项，另外一些实施例中，该值也可取128或者512,需要根据特定设备进行配置；应注意的是，检测脉冲方向又称慢时间方向；另外，本步骤只需在执行S130前设置完成即可，其与S120并无特定执行顺序限定。

S120：将获取的声辐射力回波信号I、Q数据按照深度方向顺序存储；进一步的，每个采样点数据包括I、Q两组数据；采用数据类型float2进行存储；将I、Q数据从主机端CPU读取到设备端GPU，所述I、Q数据读入GPU后存入全局存储器中；应注意的是，读取时，如果每个检测脉冲在深度方向上的采样点个数不满足128字节的整数倍，则通过补零使其满足；这是因为，并行计算核心kernel运行时，待处理的IQ数据需要从Host主机端(CPU)写入Device设备端(GPU)的全局存储器(global memory)。整个NDRange中的所有work-item(线程，又称工作项)都可以读写全局存储器的任意位置；为了能够高效的访问全局存储器，读取和写入数据必须对齐，并且一个work-group中的work-iem以合并访问(coalesced access)的方式读取连续的存储空间，否则读写会被编译器拆分为多次操作，影响算法执行效率；而不同工作项是沿着深度方向对应分配的，这保证了合并访问，同一个工作组(work-group)中的work-item访问的数据是连续的；因此在同一检测脉冲中，如果深度方向上的采样点数不是128字节的整数倍，为了使数据对齐，需要在深度方向补零以使其满足数据对齐的要求。

S130：将每个检测脉冲的采样点分配入对应的工作组中的各工作项中，计算是，由该工作组中各工作项基于自相关Loupas算法计算所述回波在每个采样点处的四个自相关参数并存储；四个自相关参数的计算公式分别是：

其中，M是深度方向的取样窗口大小，N是是检测脉冲时间方向的取样窗口大小；该四个自相关参数实际是基于自相关Loupas算法计算回波平均位移公式中的分母分子中的四个参数项，这是因为在计算深度方向平均位移时，是通过在深度采样频率快时间方向选取一个大小为M的取样窗口；在采样包方向(即慢时间方向或者说检测脉冲方向)选取一个大小为N，本实施例中，N＝2的取样窗口，其中一个对应脉冲序列中的第一个参考脉冲，另一个对应采样包方向PRF某个时刻的检测脉冲；上式中，分子分母对于每个采样点都是在深度方向作求和运算，对于相邻的采样点来说，其大部分计算是重复的；因此，在采用OpenCL算法时，我们为了方便并行处理，在每个工作项(work-item)计算相应采样点的位移时，只计算公式分子、分母求和项中的一项，即，只计算上述四个自相关参数；待所有采样点计算完成后，每个工作项(work-item)再对取样窗口中的点进行求和运算，具体如图4所示IQ(m,0)和IQ(m+1,0)对应参考脉冲在深度方向的两个点，IQ(m,1)和IQ(m+1,1)对应检测脉冲在采样包方向某个时刻的两个点。上式中分子对应IQ(m,0)与IQ(m,1)之间作互相关，分母对应IQ(m,0)与IQ(m+1,0)之间、以及IQ(m,1)与IQ(m+1,1)之间作自相关，这样每个工作项(work-item)在处理当前采样点时只涉及4个点之间的相关性计算。即，IQ(m,0)、IQ(m+1,0)、IQ(m,1)、IQ(m+1,1)。

应注意的是，大多数情况下，每个检测脉冲在深度方向上的采样点数均大于每个工作组中的工作项数，如，每个检测脉冲在深度方向上可能存在1000个采样点，而本实施例中，与每个检测脉冲对应的各个工作组中，每个工作组包含256个工作项，此时，通过每次并行处理256个采样点相关计算，计算完成后，重新顺序进行下一个256个采样点的计算，这样，大概经过四轮的计算，该检测脉冲在深度方向上的1000个采样点即可全部计算完毕，而采用传统的计算方式，则需要大概1000轮计算才能处理完所有采样点。

S140：同一工作组中的各个采样点的自相关参数均计算完成后，工作项基于自相关Loupas算法计算所述回波在深度方向取样窗口内的平均位移

其中，c是声速，f_c是发射信号中心频率，f_s是采样频率；本步骤中，每个工作项都需要在深度方向求和以计算位移，数据重复访问率高。为了提高数据访问效率，同时考虑到设备端的资源有限，为了让尽可能多的工作组(work-group)处于活动状态，本实施例中，我们使用共享存储器(shared local memory)来存储每个工作项(work-item)的计算结果，其字节大小为nSample*4*sizeof(float)，其中nSample表示采样点个数，4表示4个自相关参数；因为工作组work-group中的work-item处理的是沿深度方向连续的存储空间，不存在共享存储器的访问冲突bank conflict。

S150：将计算得出的回波平均位移数据按照深度方向顺序存储(本步骤图3中未显示)。

进一步的，中值滤波是一种非线性的平滑滤波器，可以有效的抑制脉冲噪声，并能够保留信号的边缘信息，防止边缘模糊的发生。中值滤波的基本原理是使用邻域内各点的像素值的中值作为中心像素的输出值。本实施例中，中值滤波的步骤采用如下方法完成：

S210：设定工作组大小为(16,8)；设置个工作组,即个工作组，nPoint表示每个检测脉冲在深度方向上的平均位移个数,通常，nPoint会略少于采样点数Sample,16表示在深度方向上的工作项个数，8表示在扫描线方向上的工作项个数；每一个工作项对应一个特定时刻的回波位移数据。

S220：截取任意时刻的所有平均位移数据(实际应用中，由用户选择任意一个时间点，则系统截取该时间点的所有平均位移数据进行成像处理)，同样的，待处理的位移数据需要从Host主机端(CPU)写入Device设备端(GPU)的全局存储器(global memory)。整个NDRange中的所有工作项(work-item，线程)都可以读写全局存储器的任意位置；为了能够高效的访问全局存储器，读取和写入数据必须对齐；将各个平均位移数据置入对应工作项；本实施例中，以5*5为取样窗口大小，沿深度方向采用中位数算法(MedianOfMedians)计算各取样窗口的中位数，而不必对每一个点的所有邻域数据进行排序，将每个取样窗口的中心点值用中位数替换；本实施例中，具体计算时，采用宏定义并结合循环展开方式进行，采用宏定义的方式省去了函数调用过程，通过循环展开方式减少了控制指令，提高了算法效率。

应注意的是，本实施例中，工作组work-group中的每一个work-item读取其对应的点的邻域点信息，然后计算其中值；邻域的大小即取样窗口大小，常见的可取值为(3,3),(5,3)和(5,5)等；对于每个取样窗口来说，其待处理的点数量规模有限。由于每个点都需要计算邻域中的数值中值，因此全局存储器的重复访问率很高，为了提高数据访问效率，本实施例中，我们采用共享存储器(shared local memory)来存储一个work-group中每个work-item的邻域点的信息，其字节大小为i*j*xsize*ysize*sizeof(float)，i表示在深度方向上的工作项个数，j表示在扫描线方向上的工作项个数；另外二维子空间的边缘点需要读取其邻域点的信息，超出边界的邻域点取值为零，这样就包含了二维子空间所有点的邻域信息。因为工作组work-group中的work-item处理的是沿深度方向连续的存储空间，不存在共享存储器的访问冲突bank conflict。

S230：将当前点值计算结果按照深度方向顺序存储(本步骤图3中未显示)。

进一步的，将每个取样窗口的中位数数值作为当前点的值。

进一步的，扫描转换的步骤中，采用如下方法实现：

S310：本实施例中，设立工作项的二维维度为(16,16)，设置个工作组，每个工作项对应一个目标像素点，ImageW表示目标图像水平像素，其对应扫描线方向，ImageH表示目标图像垂直像素，其对应深度方向。

S320：目标像素点坐标转换至原始数据坐标，所述原始数据坐标是指将步骤S230的计算结果按照深度方向为横轴、扫描线方向为纵轴排列的坐标；根据每个目标像素点对应的一个以上的原始数据坐标对应的数据通过双线性插值或者三次样条插值法计算得到当前点的像素值。应注意的是，本实施例中，本步骤输入数据是float类型的，输出数据是uchar类型的，由于目标图像数据(输出数据)的存储方式与原始数据(输入数据)不一致，会导致读取全局存储器时不能保证一个work-group中的work-item都能够进行合并访问，但是在ARFI成像模式下系统扫描的线密度很稀疏，多数情况下坐标变换后映射到原始数据坐标系统的数据分布是相邻的，能够满足有限次合并访问的要求，而且现有GPU设备支持全局存储器的数据缓存，实验证明算法的时间效率并没有明显下降，非常接近，因此可以在一个kernel中同时实现扫描转换和二维图像转置。

S330：将S320的计算结果按照目标像素点坐标的分布也就是行(图像宽度)优先方式存储(本步骤图3中未显示)。

进一步的，S320中，还包括设置最大位移阈值T，将获取的像素值归一化到[0,255]的像素范围内的步骤。

另外，本申请中，S110、S210及S310均不代表限定任何执行顺序，其完全可以在执行其他步骤(其他步骤指S120、S130、S140、S150、S220、S230、S320及S330)前提前先行设定完成，也可以在执行后续步骤前才完成，如在执行S120前完成S110的设定，在执行S220前完成S210的设定，在执行S320千完成S310的设定。

Claims

1.一种基于OpenCL的并行ARFI成像方法，其特征在于，位移估算的步骤采用如下方法完成：

(1-1)根据扫描线条数gx及每条扫描线中的检测脉冲个数gy设置工作组为2维的(gx,gy)，每个工作组对应某一时刻的检测脉冲；每个工作组中设置(ix,iy)＝(64×k,1)个工作项，k为2以上整数；

(1-2)将获取的声辐射力回波信号I、Q数据按照深度方向顺序存储；将I、Q数据从主机端读取到设备端，如果每个检测脉冲在深度方向上的采样点个数不满足128字节的整数倍，则通过补零使其满足；

其中，m为M中的取值，M是深度方向的取样窗口大小，n为N中的取值，N是检测脉冲时间方向的取样窗口大小；I、Q均为声辐射力回波信号数据分量；

其中，c是声速，f_c是发射信号中心频率，f_s是采样频率；

2.如权利要求1所述的ARFI成像方法，其特征在于，步骤(1-3)中每个检测脉冲中的各个采样点的四个自相参数计算完成后存入共享存储器，其大小为nSample*4*sizeof(float)，其中nSample表示采样点的个数，4表示公式中包含四项自相关参数。

3.如权利要求1所述的ARFI成像方法，其特征在于，每个采样点数据包括I、Q两组数据；采用数据类型float2进行存储。

4.如权利要求1所述的ARFI成像方法，其特征在于，所述主机端为CPU，所述设备端为GPU；所述I、Q数据读入GPU后存入全局存储器中。

5.如权利要求1所述的ARFI成像方法，其特征在于，进一步包括中值滤波的步骤，并采用如下方法完成：

(2-1)设定工作组大小为(i,j)其中i、j为相等或不等的8的整数倍；设置个工作组；nPoint表示每个检测脉冲在深度方向上的平均位移个数,i表示在深度方向上的工作项个数，j表示在扫描线方向上的工作项个数；每一个工作项对应一个特定时刻的回波位移数据；

(2-2)截取任意时刻的平均位移数据，从主机端将数据读取到设备端，每一个工作项处理其对应的一个位移数据；以当前点为中心以xsize*ysize为取样窗口大小，计算当前点在邻域范围内的中位数，将每个取样窗口的中心点值用中位数替换；xsize,ysize为3以上奇数；所有取样窗口数据存储在共享存储器中，大小为i*j*xsize*ysize*sizeof(float)；

(2-3)将当前点值计算结果按照深度方向顺序存储。

6.如权利要求5所述的ARFI成像方法，其特征在于，在从主机端向设备端读取数据时，如果每个扫描线该时刻深度方向上的平均位移数据不满足128字节的整数倍，则通过补零使其达到对齐访问的要求。

7.如权利要求5所述的ARFI成像方法，其特征在于，采用中位数算法计算每个点的中位数数值，计算时，采用宏定义并结合循环展开方式进行。

8.如权利要求5所述的ARFI成像方法，其特征在于，进一步包括扫描转换的步骤，并采用如下方法实现：

(3-1)设定工作组大小为(x,y)，x，y为相等或不等的16的整数倍；设个工作组，每个工作项对应一个目标像素点,ImageW表示目标图像水平像素，其对应扫描线方向，ImageH表示目标图像垂直像素，其对应深度方向；

(3-2)目标像素点坐标转换至原始数据坐标，所述原始数据坐标是指将步骤(2-3)的计算结果按照深度方向为横轴、扫描线方向为纵轴排列的坐标；根据每个目标像素点对应的一个以上的原始数据坐标对应的数据计算该目标像素点的像素值；

(3-3)将步骤(3-2)的计算结果按照目标像素点坐标的分布存储。

9.如权利要求8所述的ARFI成像方法，其特征在于，步骤(3-2)中是通过双线性插值或者三次样条插值法计算得到当前点的像素值。

10.如权利要求8所述的ARFI成像方法，其特征在于，步骤(3-2)中，还包括设置最大位移阈值T，将像素值归一化到[0,255]的像素范围内的步骤。