基于统一计算设备架构技术的图像处理方法和装置 【技术领域】
本发明涉及图像处理技术领域，特别涉及基于统一计算设备架构(CUDA，Compute Unified Device Architecture)技术的图像处理方法和装置。
背景技术
磁共振成像(MRI)是一种生物磁自旋成像技术，在外加磁场内原子核的自旋运动经射频脉冲冲激后产生信号。通过MRI数据采样技术将这些信号用探测器检测并输入计算机；MRI重建技术对信号进行后处理转换成在屏幕上显示的图像。
MRI数据采集是在空间频率域即K空间中进行的。为了将这些K空间采样数据转换到二维或三维图像空间坐标系如笛卡尔坐标系下，必须对数据进行反傅立叶变换。但是标准的反傅立叶变换算法要求数据是在等间隔坐标系统中均匀采样的，这就要求MRI重建算法在运用反傅立叶变换前将采样数据从K空间变换到等间隔坐标系统如笛卡尔坐标系下以实现均匀采样。这个变换的过程通常是通过网格化算法实现的。
网格化算法通过数据重采样技术将非等间隔采样的K空间数据转换到笛卡尔坐标系下，它把K空间划分为等间隔的直线网格，给每一个网格单元的中心赋一个值，该值等于落在该网格单元内的全部数据之“和”。这里的“和”并不是网格单元内全部数据的简单相加，而是一种卷积插值运算。
网格化算法实际上是一种卷积插值运算，每一个采集到的样本都与一个扩散核进行卷积，该样本的能量在卷积后也就分布到了相邻的网格点上，这样相当于将该样本值在附近的网格点上进行了重采样。这种卷积插值算法的优点是插值过程中用到了全部采样数据。目前，Kaiser-Bessel窗函数已成为MRI网格化算法中普遍采用的卷积函数。
计算统一设备架构(CUDA，Compute Unified Device Architecture)是新一代的图像处理芯片(GPU)架构，可以将计算进行多线程任务分配和管理，可以极大程度地发挥GPU的优势，十分方便地对大规模并行计算进行硬件加速。
现有的，网格化算法用公式表示为：
MRS={[M(u,v)·(S(u,v)S(u,v)*C(u,v))]*C(u,v)}·R(u,v)---(2-1)]]>
式中M(u，v)为磁共振采样数据，MRS为重采样后的数据，S(u，v)为不等间隔的采样函数，R(u，v)表示等间隔直线网格采样函数，C(u，v)为卷积函数(Kaiser-Bessel函数)，·表示矩阵乘法，*表示卷积运算，实际上是非均匀采样密度补偿函数。不等间隔采样函数S(u，v)为复数，它的实部用x表示，虚部用y表示，计算公式如(2-2)-(2-3)所示：
S(u,v).x=-corX2+corXSampSizw*u---(2-2)]]>
S(u,v).y=-corY2+corYScanMatrix*(ScanMatrix-ETL2+v)---(2-3)]]>
其中corX、corY分别代表笛卡尔坐标X、Y轴向归一化长度，SampSize代表采样点数，ScanMatrix代表采样扫描矩阵大小，ETL代表回波链长度，0≤u＜SampSize，0≤v＜ETL。
K空间卷积窗大小[kwidth.x，kwidth.y]计算公式如(2-4)-(2-5)所示：
kwidth.x＝k max.x*convwidth.x/gridsize.x (2-4)
kwidth.y＝k max.y*convwidth.y/gridsize.y (2-5)
其中，(gidsize.x，gridsize.y)表示网格化结果数据规模即网格化结果数据矩阵的大小，(convwidth.x，convwidth.y)表示图像空间卷积窗口大小，(k max.x，k max.y)代表K空间归一化大小。
K空间卷积窗口大小为[kwidth.x，kwidth.y]时，采样轨迹点S(u，v)对应网格化结果矩阵卷积窗口坐标范围[(ix min，iy min)，(ix max，iy max)]的计算公式如(2-6)-(2-9)所示：
ixmin=(int)((S(u,v).x-kwidth.x)*gridsize.xkmax.x+gridcenter.x)---(2-6)]]>
if(ix min＜0)ix min＝0；
ixmax=(int)((S(u,v).x+kwidth.x)*gridsize.xkmax.x+gridcenter.x)+1---(2-7)]]>
if(ix max≥gridsize.x)ix max＝gridsize.x-1；
iymin=(int)((S(u,v).y-kwidth.y)*gridsize.ykmax.y+gridcenter.y)---(2-8)]]>
if(iy min＜0)iy min＝0；
iymax=(int)((S(u,v).y+kwidth.y)*gridsize.ykmax.y+gridcenter.y)+1---(2-9)]]>
if(iy max≥gridsize.y)iy max＝gridsize.y-1；
上述公式中，min、max为求最小、最大值的函数，(gidcenter.x，gridcenter.y)表示网格化结果中心点位置，(convwidth.x，convwidth.y)表示图像空间卷积窗口大小，(kwidth.x，kwidth.y)表示K空间卷积窗口大小，其余变量与前述相同，不再重复说明。
计算网格化结果矩阵卷积坐标窗口[(ix min，iy min)，(ix max，iymax)]内的点与窗口中心点的K空间欧氏距离，计算公式如(2-10)-(2-12)所示：：
dkx=kmax.x*grid(p,q).x-gridcenter.xgridsize.x-S(u,v).x---(2-10)]]>
dky=kmax.y*grid(p,q).y-gridcenter.ygridsize.y-S(u,v).y---(2-11)]]>
dk=dkx2+dky2---(2-12)]]>
其中，(p，q)∈[(x min，y min)，(x max，y max)]，(grid(p，q).x，grid(p，q).y)代表(p，q)点的网格化数据，dk表示窗内元素距离卷积窗中心点的K空间欧氏距离；对dk<(kwidth.x)2+(kwidth.y)2]]>的网格化卷积窗内点，对采样轨迹S(u，v)对应的采样数据点进行卷积运算。
用dataw表示采样轨迹列数，datah表示采样轨迹行数，grid(p，q)表示(p，q)点网格化结果，网格化算法的主要过程为：
(1)预处理：S(u，v)采样轨迹生成、M(u，v)数据采样
(2)令i＝0至datah，每次加1，即，for i＝0→datah
(3)令j＝0至dataw，每次加1，即，for j＝0→dataw
(4)利用公式(2-6)到(2-9)计算采样轨迹点(i，j)的K空间卷积窗大小[kwidth.x，kwidth.y]转换到图像空间对应网格化矩阵坐标窗口[(ix min，iy min)，(ix max，iy max)]；
(5)令p＝iy min→iy max，例如，for p＝iy min→iy max
(6)令q＝ix min→ix max，例如，for q＝ix min→ix max
(7)利用公式(2-10)到(2-12)计算dk，如果dk≥(kwidth.x)2+(kwidth.y)2]]>转(5)；
(8)利用公式(2-1)对结果数据grid(p，q)进行卷积运算，获得网格化结果矩阵。
图1是现有技术中网络化算法卷积窗口计算原理示意图。从图1可知，卷积窗口是通过轨迹计算得到的，即通过轨迹1计算得到卷积窗口I，通过轨迹2计算得到卷积窗口II。CUDA技术可以进行多线程任务分配和管理，通过CUDA实现现有技术的网络化算法，即对网格化算法进行并行处理，上述卷积窗口I和卷积窗口II的重叠部分存在被两组计算(对轨迹1的计算和对轨迹2的计算)同时写入的情况，即在网格化算法的主要过程(4)中，对卷积窗口的运算存在交叉现象，由此导致网格化算法存在写-写冲突。
也就是说，由于CUDA技术是对网格化算法并行处理，因而以上述步骤(2)和(3)为初始数据的计算过程由多线程并行完成，这样，会导致同一组输入数据，应用相同的算法后每次执行结果都存在差异(由窗口I和窗口II的重叠部分导致)，且这种差异没有任何规律性。例如，表1是基于CUDA的网格化算法执行环境，表2是网格化算法在表1所示执行环境下存在写-写冲突的代表性数据：第一列“位置信息”表示网格化点的坐标位置；第二列是CPU网格化结果；第三、四列是基于CUDA的网格化算法运行结果；从表中椭圆标注数据可见，基于CUDA的网格化算法结果每次执行结果都不一致，且没有任何规律。
表1
表2
可见，在执行环境相同、输入参数相同的情况下基于CUDA技术的网格化算法执行结果是不稳定的，即每次执行结果都有一定差异，且该差异是没有任何规律的，因而对导致该算法结果进行反傅立叶变换后得到的图像也是每次都不一致的，且是没有任何规律的，从而导致后续依赖于该图像的处理结果都是不准确的。
【发明内容】
本发明实施例在于提供一种基于统一计算设备架构技术的图像处理方法和装置，以避免并行处理中存在的写-写冲突，使得基于CUDA技术的网格化算法执行结果是稳定地，从而保证最终每次得到完全一致的图像。
本申请提供了一种基于统一计算设备架构技术的图像处理方法，包括：
获取基础数据，根据成像分辨率确定网格化结果数据规模；
根据所述网格化结果数据规模获得K空间卷积窗口及窗内所有元素的坐标；
依次根据所述K空间卷积窗内每一个元素的坐标获得该元素的轨迹值，应用所述每个元素的轨迹值依次和K空间中心点的轨迹值计算欧氏距离；
若所获得的欧氏距离小于第一阈值，则结合所述基础数据对该元素值所对应的基础数据中的采样数据进行卷积计算，获得网格化计算结果数据；
对所述网格化计算结果数据进行反傅里叶变换，获得图像数据。
其中，所述基础数据包括：采样数据、密度补偿数据、采样轨迹数据和卷积核数据。
其中，根据所述网格化结果数据规模获得K空间卷积窗口及窗内所有元素的坐标的步骤包括：根据网格化结果数据规模，创建执行核函数，由所述执行核函数进行以下步骤：
根据网格化结果数据规模，获得网格化结果数据的中点坐标；
根据所述中点坐标计算K空间卷积窗口四个顶点轨迹值；
根据所述K空间卷积窗口顶点轨迹值计算窗内四个顶点的坐标，
根据四个顶点的坐标确定窗内所有元素的坐标。
其中，根据所述中点坐标计算K空间卷积窗口四个顶点轨迹值的步骤包括：
xmin=grid.x-gridcenter.x-convwidth.xgridsize.x*kmax.x]]>
xmax=grid.x-gridcenter.x+convwidth.xgridsize.x*kmax.x]]>
ymin=grid.y-gridcenter.y-convwidth.ygridsize.y*kmax.y]]>
ymax=grid.y-gridcenter.y+convwidth.ygridsize.y*kmax.y]]>
其中，grid.x表示当前网格化点的X方向坐标值；grid.y表示当前网格化点的Y方向坐标值；gridsize.x表示网格化结果X方向大小；gridsize.y表示网格化结果Y方向大小；gridcenter.x表示网格化结果中心点X方向坐标值；gridcenter.y表示网格化结果中心点Y方向坐标值；convwidth.x表示图像空间卷积窗口X方向大小；convwidth.y表示图像空间卷积窗口Y方向大小；kmax.x表示K空间X方向最大值；kmax.y表示K空间Y方向最大值。
其中，根据所述K空间卷积窗口顶点轨迹值计算窗内四个顶点的坐标：
umin=(xmin+corX2)*SampSizecorX]]>
umax=(xmax+corX2)*SampSizecorX]]>
vmin=(ymin+corY2)*ScanMatrixcorY-ScanMatrix-ETL2]]>
vmax=(ymax+corY2)*ScanMatrixcorY-ScanMatrix-ETL2]]>
其中，corX、corY分别代表笛卡尔坐标X、Y轴向归一化长度，SampSize代表采样点数，ScanMatrix代表采样扫描矩阵大小，ETL代表回波链长度，0≤u＜SampSize，0≤v＜ETL。
其中，应用所述每个元素的轨迹值依次和K空间中心点的轨迹值计算欧氏距离的步骤包括：
i)计算一个元素的轨迹值和K空间中心点的轨迹值计算欧氏距离dk，具体为：
dkx=kmax.x*grid.x-gridcenter.xgridsize.x-S(p,q).x]]>
dky=kmax.y*grid.y-gridcenter.ygridsize.y-S(p,q).y]]>
dk=dkx2+dky2]]>
其中，(p，q)∈[(x min，y min)，(x max，y max)]，S(p，q)表示采样函数S(u，v)在(p，q)点的函数值，dk表示窗内元素距离卷积窗中心点的K空间欧氏距离；grid.x表示当前网格化点的X方向坐标值；grid.y表示当前网格化点的Y方向坐标值；gridsize.x表示网格化结果X方向大小；gridsize.y表示网格化结果Y方向大小；gridcenter.x表示网格化结果中心点X方向坐标值；gridcenter.y表示网格化结果中心点Y方向坐标值；k max.x表示K空间X方向最大值；k max.y表示K空间Y方向最大值；
ii)重复步骤i)，依次计算每个元素的轨迹值和K空间中心点的轨迹值计算欧氏距离dk。
其中，所述第一阈值d为：
d=(convwidth.xgridsize.x*kmax.x)2+(convwidth.ygridsize,y*kmax.y)2]]>
其中，convwidth.x表示图像空间卷积窗口X方向大小；convwidth.y表示图像空间卷积窗口Y方向大小，gridsize.x表示网格化结果X方向大小；gridsize.y表示网格化结果Y方向大小，k max.x表示K空间X方向最大值；k max.y表示K空间Y方向最大值。
本申请还提供了一种基于统一计算设备架构技术的图像处理装置，包括：
获取模块，用于获取基础数据；
数据规模确定模块，用于根据成像分辨率确定网格化结果数据规模；
坐标计算模块，用于根据所述网格化结果数据规模获得K空间卷积窗口及窗内所有元素的坐标；
轨迹计算模块，用于依次根据所述K空间卷积窗内每一个元素的坐标获得该元素的轨迹值；
欧氏距离计算模块，用于应用所述每个元素的轨迹值依次和K空间中心点的轨迹值计算欧氏距离；
网格化计算模块，用于在所获得的欧氏距离小于第一阈值时，结合所述基础数据对该元素值所对应的基础数据中的采样数据进行卷积计算，获得网格化计算结果数据；
图像数据计算模块，用于对所述网格化计算结果数据进行反傅里叶变换，获得图像数据。
其中，所述基础数据包括：采样数据、密度补偿数据、采样轨迹数据和卷积核数据。
其中，所述坐标计算模块包括：
中心点坐标计算模块，用于根据网格化结果数据规模，获得网格化结果数据的中点坐标；
顶点轨迹计算模块，用于根据所述中点坐标计算K空间卷积窗口四个顶点轨迹值；
顶点坐标计算模块，用于根据所述K空间卷积窗口顶点轨迹值计算窗内四个顶点的坐标；
窗内元素坐标计算模块，用于根据四个顶点的坐标确定窗内所有元素的坐标。
本发明实施例很好的解决了基于CUDA的网格化算法由并行处理所导致的写-写冲突。CUDA实现结果与CPU实现结果仅在小数据点后2-6位存在误差；输入条件一致时，执行结果也是唯一的，从而使得基于CUDA技术的网格化算法执行结果是稳定的，在保证输入条件一致的情况下最终得到完全一致的图像。
【附图说明】
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是现有技术中网络化算法卷积窗口计算原理示意图；
图2是根据本发明实施例的网格化算法卷积窗计算原理示意图；
图3是根据本发明实施例的基于CUDA技术的图像处理方法流程图；
图4是根据本发明实施例的基于CUDA技术的图像处理装置结构示意图；
图5是CPU网格化算法成像结果图；
图6是应用本发明实施例的基于CUDA的网格化算法成像结果图。
【具体实施方式】
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明实施例重点解决了基于CUDA的网格化算法的写-写冲突，其主要思想是：以网格化结果为起点的任务分配策略。即算法以网格化结果点计算卷积窗口，对落在窗内的采样轨迹S(u，v)和磁共振采样数据M(u，v)卷积计算该点的网格化值。由于网格化点的位置唯一性，其对应卷积窗口也是唯一的，所以在并行处理中不存在写-写冲突，使得基于CUDA技术的网格化算法执行结果是稳定的，从而保证了每次得到完全一致的最终图像。
图2是根据本发明实施例的网格化算法卷积窗计算原理示意图。本实施例中，是以现有的网格化结果点为起始点，例如，利用落在窗口I中的采样轨迹S(u，v)和磁共振采样数据M(u，v)卷积计算该点的网格化值，获得结果数据1，利用落在窗口II中的采样轨迹S(u，v)和磁共振采样数据M(u，v)卷积计算该点的网格化值，获得结果数据2，即在输入数据确定的情况下，输出数据是确定值，因而避免了写-写冲突。也就是说，由于网格化点的位置是唯一的，即便其对应卷积窗口存在交叉现象，并行处理时网格化数据点也不存在写-写冲突。
参见图3，其是根据本发明实施例的基于CUDA技术的图像处理方法流程图，具体包括：
步骤301，获取基础数据，根据成像分辨率确定网格化结果数据规模。
具体的，上述基础数据包括算法输入所需要的所有数据：采样数据M(u，v)、密度补偿数据采样轨迹数据S(u，v)、卷积核数据C(u，v)。这里，这些基础数据需载入GPU显存。
上述网格化结果数据规模是根据载入的基础数据和成像结果分辨率[imagesize.x，imagesize.y]决定的。本发明实施例中，网格化结果数据规模和图像分辨率的大小相等。网格化结果数据规模可以表示为[gridsize.x，gridsize.y]。简单的说，网格化结果数据规模就是指网格化结果矩阵的大小。例如，如果成像分辨率是512*512的，那么网格化结果数据的规模也是512*512。
步骤302，根据所述网格化结果数据规模获得K空间卷积窗口及窗内所有元素的坐标。
具体为：根据网格化结果数据规模，创建执行核函数，由所述执行核函数进行以下步骤：
i)根据网格化结果数据规模，获得网格化结果数据的中点坐标；
ii)根据所述中点坐标计算K空间卷积窗口四个顶点轨迹值；
iii)根据所述K空间卷积窗口顶点轨迹值计算窗内四个顶点的坐标；
iv)根据四个顶点的坐标确定窗内所有元素的坐标。
其中，上述创建执行核函数的具体步骤包括：通过CUDAAPI调用，返回GPU硬件参数信息，根据这些硬件参数信息选择线程创建策略，从而创建执行核函数。该创建执行核函数本身的过程是现有技术，这里不再详述。
当执行核函数创建完毕后，令每个网格化数据点对应一个线程。将网格化结果数据规模按blocksize.x×blocksize.y分块，[blocksize.x，blocksize.y]代表块内线程规模，这样一共产生个线程块，每个网格化数据点(即网格中的一个点)对应一个线程。块内线程不应该多于GPU多处理器组(SM)最多可同时执行的线程个数，但是线程块个数可以多于流处理器(SP)个数。
上述步骤ii)中根据所述中点坐标计算K空间卷积窗口四个顶点轨迹值的具体实现过程是：
xmin=grid.x-gridcenter.x-convwidth.xgridsize.x*kmax.x---(2-13)]]>
xmax=grid.x-gridcenter.x+convwidth.xgridsize.x*kmax.x---(2-14)]]>
ymin=grid.y-gridcenter.y-convwidth.ygridsize.y*kmax.y---(2-15)]]>
ymax=grid.y-gridcenter.y+convwidth.ygridsize.y*kmax.y---(2-16)]]>
上述步骤iii)中根据所述K空间卷积窗口顶点轨迹值计算窗内四个顶点的坐标的具体实现过程是：
umin=(xmin+corX2)*SampSizecorX---(2-17)]]>
umax=(xmax+corX2)*SampSizecorX---(2-18)]]>
vmin=(ymin+corY2)*ScanMatrixcorY-ScanMatrix-ETL2---(2-19)]]>
vmax=(ymax+corY2)*ScanMatrixcorY-ScanMatrix-ETL2---(2-20)]]>
其中，corX、corY分别代表笛卡尔坐标X、Y轴向归一化长度，SampSize代表采样点数，ScanMatrix代表采样扫描矩阵大小，ETL代表回波链长度，0≤u＜SampSize，0≤v＜ETL，未说明变量与前述相同。
步骤303，依次根据所述K空间卷积窗内每一个元素的坐标获得该元素的轨迹值，应用所述每个元素的轨迹值依次和K空间中心点的轨迹值计算欧氏距离。
其中，计算一个元素的轨迹值和K空间中心点的轨迹值计算欧氏距离dk，具体为：
dkx=kmax.x*grid.x-gridcenter.xgridsize.x-S(p,q).x---(2-21)]]>
dky=kmax.y*grid.y-gridcenter.ygridsize.y-S(p,q).y---(2-22)]]>
dk=dkx2+dky2]]>
其中，(p，q)∈[(x min，y min)，(x max，y max)]，S(p，q)表示采样函数S(u，v)在(p，q)点的函数值，dk表示窗内元素与卷积窗口中心点的K空间欧氏距离；，未说明变量与前述相同。
重复前述步骤，即可依次计算每个元素的轨迹值和K空间中心点的轨迹值计算欧氏距离dk。
步骤304，若所获得的欧氏距离小于第一阈值，则结合所述基础数据对该元素值所对应的基础数据中的采样数据进行卷积计算，获得网格化计算结果数据。
其中，第一阈值d的计算公式为：
d=(convwidth.xgridsize.x*kmax.x)2+(convwidth.ygridsize.y*kmax.y)2---(2-24)]]>
结合前述基础数据，对坐标范围[(u min，v min)，(u max，v max)]内满足条件dk＜d的采用轨迹进行卷积运算，得到网格化计算结果数据，即网格化中数据点(grid.x，grid.y)处的值。
步骤305，对所述网格化结果数据进行反傅里叶变换，获得图像数据。至此，得到了基于CUDA的磁共振图像。
应用本发明实施例提供的基于CUDA技术的图像处理方法，很好的解决了基于CUDA的网格化算法写-写冲突。CUDA实现结果与CPU实现结果仅在小数据点后2-6位存在误差；输入条件一致时，执行结果也是唯一的，从而保证在输入条件一致的情况下最终得到图像每次都是完全一致的。
本发明实施例还提供了一种基于计算统一设备架构技术的图像处理装置，参见图4，具体包括：
获取模块401，用于获取基础数据；
数据规模确定模块402，用于根据成像分辨率确定网格化结果数据规模；
坐标计算模块403，用于根据所述网格化结果数据规模获得K空间卷积窗口及窗内所有元素的坐标；
轨迹计算模块404，用于依次根据所述K空间卷积窗内每一个元素的坐标获得该元素的轨迹值；
欧氏距离计算模块405，用于应用所述每个元素的轨迹值依次和K空间中心点的轨迹值计算欧氏距离；
网格化计算模块406，用于在所获得的欧氏距离小于第一阈值时，结合所述基础数据对该元素值所对应的基础数据中的采样数据进行卷积计算，获得网格化计算结果数据；
图像数据计算模块407，用于对所述网格化计算结果数据进行反傅里叶变换，获得图像数据。
上述基础数据包括：采样数据、密度补偿数据、采样轨迹数据和卷积核数据。
上述坐标计算模块403可以具体包括：
中心点坐标计算模块，用于根据网格化结果数据规模，获得网格化结果数据的中点坐标；
顶点轨迹计算模块，用于根据所述中点坐标计算K空间卷积窗口四个顶点轨迹值；
顶点坐标计算模块，用于根据所述K空间卷积窗口顶点轨迹值计算窗内四个顶点的坐标；
窗内元素坐标计算模块，用于根据四个顶点的坐标确定窗内所有元素的坐标。
应用本发明实施例提供的基于CUDA技术的图像处理装置，很好的解决了基于CUDA的网格化算法写-写冲突。CUDA实现结果与CPU实现结果仅在小数据点后2-6位存在误差；输入条件一致时，执行结果也是唯一的，从而保证在输入条件一致的情况下最终得到图像每次都是完全一致的。
在表1所述执行环境下，应用本发明实施例所提供的方法，得到如表3所示基于CUDA的逆向网格化算法数据：
表3
在采用表1所示执行环境下，应用CPU与本发明实施例所述方法得到的图像分别如图5和图6所示。图5是CPU网格化算法成像结果，图5(a)-图5(c)依次采用的是运行的头部数据、静止的头部数据和静止水模数据；图6是基于本发明实施例的CUDA逆向网格化算法成像结果，图6(a)-图6(c)依次采用的是运行的头部数据、静止的头部数据和静止水模数据。从以上三组图像的结果可以看出，CUDA逆向网格化算法重建结果和CPU网格化算法重建结果完全一致。
需要说明的是，本文所涉及公式较多，为简洁凡是公式中所涉及的变量未做说明的，可参考前面公式的说明。
需要说明的是，对于系统实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。
需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个......”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。
本领域普通技术人员可以理解实现上述方法实施方式中的全部或部分步骤是可以通过程序来指定相关的硬件来完成，所述的程序可以存储于计算机可读取存储介质中，这里所称得的存储介质，如：ROM/RAM、磁碟、光盘等。
以上所述仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等，均包含在本发明的保护范围内。