专利名称：一种基于sirt算法且使用cpu+mic协同计算模式的电镜三维重构方法
技术领域：
本发明涉及计算机应用技术领域，具体地说是一种基于SIRT算法且使用CPU+MIC协同计算模式的电镜三维图像重构方法，尤其涉及在MIC众核架构上使用CPU+MIC协同计算模式进行电镜三维重构的方法及装置。
背景技术：
MIC (Many Integrated Core)是Intel公司推出的众核处理器,跟通用的多核至强处理器相比，MIC众核架构具有更小的内核和硬件线程，众核处理器计算资源密度更高，片上通信开销显著降低，更多的晶体管和能量，能够胜任更为复杂的并行应用。Intel MIC产品基于X86架构,基于重核的众核处理器,包含50个以上的核心，以及512bit的向量位宽，双精性能超过lTFlops。MIC拥有极其灵活的编程方式，MIC卡可以作为一个协处理器存在，也可以被看作是一个独立的节点。基本的MIC编程模型是将MIC看作一个协处理器，CPU根据程序的指令，将一部分代码运行在MIC端。此时存在两类设备，即CPU端和MIC众核协处理器端。电子断层三维重构技术(Electron Tomography, ET)是从一个物体的二维投影图经过重构获取物体内部结构的技术，通过获取多个不同角度的二维投影图进行反向重构出所研究对象的三维结构。图像重构无论是在医学诊断，还是在生物样品的成像，以及在地表层析成像等诸多领域中都有着广泛的应用。迭代法可以表述成Ax=p，其中A为非奇异矩阵，X为未知数向量，P为已知向量，迭代法求解目标在于确定X的数值。在迭代重构图像方法中，A为加权因子，P为投影得到的图像像素值，X则为需要重构的图像的像素值。图I说明了采用迭代重构图像方法实现图像重构的过程。该方法通过赋O值或随机值或背投影方法获取用于迭代过程中的初值χ(°)，然后经过多次的重投影和背投影的迭代过程，来求得最后满足精度要求的重构图像的像素值。 SIRT (联合迭代重构)算法在所有图像重构算法中重构效果较好，但其运算量巨大，因此需要提高运算效率，来满足快速图像重构的需求。可见，目前需要针对基于MIC平台设计一种迭代重构图像的方法，通过MIC高效简单快速的计算特性，能够极大地加速迭代重构方法的执行，从而快速地完成图像重构任务。
发明内容
本发明的目的是提供一种使用CPU+MIC协同计算模式的电镜三维图像重构方法。能够基于MIC众核架构并且使用CPU+MIC协同计算的快速图像重构。本发明的目的是按以下方式实现的，包括CPU端，MIC众核协处理器端以及CPU+MIC协同计算模式，其中
(PU端负责将要进行迭代重构图像中的物体分成切片单元，向MIC卡传递重构所需要的值，CPU+MIC协同计算模式的框架搭建以及任务调度和参数初始化工作，而且在整个三维图像重构的计算任务CPU也会以Openmp多线程模式，依次通过重投影算法和背投影算法计算获取切片的重构图像值来动态参与；
MIC众核协处理器负责多线程并行地根据切片的重构图像值和测量获取的该切片投影图像值，依次通过重投影算法和背投影算法计算获取切片的重构图像值，在MIC卡上也采用openmp多线程的方式来运算；
重构步骤如下
(PU端将要进行迭代重构图像中的物体分成切片，向MIC卡上传递重构计算所需要的值，以openmp多线程模式，并行执行迭代重构图像，具体包括
根据物体的大小将该物体分成切片，切片的数目为该物体宽度的像素值SY，切片在投影图像上的元素数为N = SX*ANG_NUM，SX为该物体长度的像素值，ANG_NUM为该物体旋转角度的总数；切片在重构图像上的元素数为M = SX*SZ，SZ为该物体高度的像素值；· CPU端将为图像数据的角度文件，thita等文件在MIC卡上开辟空间，确定其输入输出类型并且将其值offload到MIC卡上，CPU负责CPU+MIC协同计算模式的框架搭建以及任务调度，具体步骤包括
单节点服务器采用由双路6核CPU和2块KNF MIC卡组成的桌面服务器，在CPU+MICs协同计算中，把双路CPU和MIC卡都作为计算设备，这样每个单节点就相当于有3个计算设备，每个设备通过一个OpenMP线程进行控制；
假设单节点上，MIC卡的总数为M个，则节点上共M+1个设备，采用数据动态划分的方式，每次各个设备自动获取切片的输入数据，片数目>=1，数据处理完之后自动获取下一个切片数据，直到所有设备计算完成所有样本切片数据；
在MIC众核协处理器端通过赋O值或随机值或背投影方法获取切片的初始的重构图像值X (O)，然后再通过重投影算法和背投影算法计算获取切片的重构图像值，具体步骤包括
根据赋初值得到的X (O)或上次迭代得到的X (k-Ι)进行重新投影得到一个投影图像值P计；根据赋初值得到的X(O)或上次迭代得到的x(k-l)，测量得到的投影图像值P测和重投影计算得到的投影图像值P计进行反向投影得到重构图像X (k)，背投影主要是进行X值的更新操作，通过多次的重投影和背投影操作可以得到满足精度要求的重构图像；
在图像重构中，测量物体的投影图像像素值P测的过程中，记物体的长为SX，宽为SY，高为SZ，旋转角度的总数为ANG_NUM，垂直于Y方向的面即为切片，共SY个切片，重构的过程以切片为单位，即每次重构一个切片，重构SY个切片就可以重构出物体模型；
分析串行迭代法中的性能瓶颈，找到耗时的重投影和背投影的代码，测试其在整个程序中所占的时间比例；
在每个切片的重构过程中，对于重投影过程，就是根据像素值X求得一个计算得到的投影值P测，P测共有SX*ANG_NUM个元素，这SX*ANG_NUM个元素的求解中没有任何的数据依赖性，以SX*ANG_NUM的并行度进行重投影算法的并行，在内核中，每个线程负责计算P测中的一个元素；
在每个切片的重构过程中，对于背投影过程，也就是求重构图像像素值X的过程，该算法就是根据P计和P测通过背投影的方式求重构图像像素值X，对于每一个切片，X共有SX*SY个元素,这SX*SY个元素的求解中没有任何的数据依赖性，以SX*SY的并行度进行重投影算法的并行，在内核中，每个线程负责计算X中的一个元素。
图I是采用迭代重构图像方法实现图像重构的流程图。
具体实施例方式
·
以下结合附图和优选实施例对本发明的技术方案进行详细地阐述。以下例举的实施例仅用于说明和解释本发明，而不构成对本发明技术方案的限制。本发明首先分析用SIRT串行方法进行电镜三维图像重构的性能瓶颈以及将串行方法移植到其他高性能平台上的难度，找到耗时的热点代码，测试其在电镜三维图像重构的整个过程中所占的时间比例以及分析提高整个方法性能的难度以及开发轴端。测试结果表明在采用SIRT联合迭代法进行图像重构时，大部分的时间均消耗在重投影和背投影两步骤中，这两个步骤不仅执行次数多，而且要计算的数据量大，同时，数据的计算要满足单指令多数据(SIMD, Single Instruction Multiple Data)的特点。因此，本发明将这两步骤作为加入MIC众核协处理器进行并行改进的重点，同时由于CPU端除了进行任务调度等功工作之余，也参与赋初值，再通过重投影算法和背投影算法计算获取切片的重构图像值的工作，即CPU+MIC协同计算的工作模式。。在图像重构中，测量获取物体的投影图像像素值Pm的过程，记物体的长为SX像素，宽为SY像素,高为SZ像素,旋转角度的总数(即旋转次数)为ANG_NUM,垂直于Y方向的每一个面为切片，共有SY个切片。图像重构的过程是以切片为单位进行的，其过程如图I所示，即对每一切片通过迭代多次计算和测量获取切片的投影图像像素值来重构；重构了 SY个切片便重构出该物体模型。
实施例本发明针对上述分析，提供了一种基于SIRT算法且使用CPU+MIC协同计算模式的电镜三维重构并行方法的实施例，包括如下步骤
步骤一 CPU+MIC协同计算框架搭建
在单节点上，共有M+1个设备(一个CPU+M个MIC卡)，采用OpenMP的fork-join模式搭建单节点上的框架，当程序开始执行的时候只有一个主线程存在，需要进行并行计算的时候，主线程派生出附加线程，即启用M+1个OpenMP线程，通过while(l)循环各个主动进行输入数据的分发读写，(ΓΜ-1号进程控制所有MIC设备，M号线程控制CPU设备。如在本专利中，CPU+2MICS平台上，主线程控制输入数据的动态分发，O号线程控制MICO设备，I号线程控制MICl设备，2号线程控制CPU设备。在SIRT算法中，MIC卡上的内存空间足够存放切片数据的大小，所以不需要采取分块处理的方式，每个设备重构完一个切片数据之后，输出重构图像，在每次读取当前切片的时候，需进行上锁操作，保证每片切片处理的独立性和对切片的读写冲突，各个设备间没有相互依赖。单节点上CPU+MICs协同计算的外围框架的伪代码如下所示
//SY为切片数目
int DEVICE_NUM=3 // 为 CPU 和 MIC 卡的总个数int Slice_private=0; //为每个线程定义一个私有的切片变量int Slice=O;//定义一个共享变量，来控制切片的读写
omp_lock_t lck;
#pragma omp parallel for private (Slice_private) , num_threads (DEVICE_NUM) for (thread=。； thread<DEVICE_NUM; thread++)
{
while (I) //死循环控制三个线程不停地进行工作
{
if ( Slice>SY )
{
break;
}
omp_set_lock (&lck);
{
Slice_private=slice;
Slice++; //共享变量++
}
omp_destory_lock(&lck); omp_set_lock (&lck);
Read(Slice_private) ;//读取第slice_private个切片的投影图像值omp_destory_lock(&lck);
if ( thread == 0 ) //mic 0 compute {
#pragma offload target(mic:0) \ in (SX, SZ, ANGLE_NUM, N, ···) \ in (P测length (ANGLE_NUM * SX)) \ in (thita: length(ANGLE_NUM)) \ in (P计length (ANGLE_NUM * SX)) \ inout (X: length (SZ*SX))
BackPro ject_Init_kernel () ; // 背投影赋初值 X(0) for: k=l: N//N代表要迭代的次数
ReProject_kernel () ; //重投影BackPro ject_kernel () ; //背投影写回第sliCe_private个切片重构后的图像值X(N)；
}
else if(thread==DEVICE_NUM-l) //cpu compute {
BackPro ject_Init_kernel () ; // 背投影赋初值 X(0) for: k=l: N//N代表要迭代的次数ReProject—kernel () ; // 重投影 BackPro ject—kernel () ; // 背投影 写回第sliCe_private个切片重构后的图像值X(N)；
}
else //other mic compute
{
#pragma offload target (mic: I)\ in (SX, SZj ANGLE—NUM，Nj ···) \ in (P 测length (ANGLE—NUM * SX)) \in (thita: length (ANGLE—NUM)) \ in (P 计length (ANGLE—NUM * SX)) \ inout (X: length (SZ*SX))
BackProject—Init—kernel () ; // 背投影赋初值 Xic0 for: k=l: N//N代表要迭代的次数
ReProject—kernel () ; // 重投影 BackPro ject—kernel () ; // 背投影 写回第sliCe_private个切片重构后的图像值X(N)；
}
} H end while 步骤二 CPU/MIC内核实现
在CPU+MICs协同计算平台上，每个节点上有双路6核CPU，总共12个核，如果超线程打开的话可以达到24个线程，在CPU内核设计中可以采用OpenMP开启多线程加速计算。在MIC内核设计中，根据MIC卡的核数目设计OpenMP的并行线程数。CPU和MIC采用相同的内核计算，伪代码如下
算法I :SIRT_MIC整体实现 begin
1for: y=0:SY-l //y代表每一个切片
2读取第y个切片的投影图像值；若运行在MIC卡上，则offload
3#pragma offload target \
4in (SX, SZ, ANGLE_NUM, N, ...) \
5in (P测length (ANGLE_NUM * SX)) \
6in (thita: length(ANGLE_NUM)) \
7in (P计length (ANGLE_NUM * SX)) \
8inout (X: length(SZ*SX))
否则
9
{
10背投影赋初值x(c0;
11for: k=l: N //N代表要迭代的次数12重投影 SIRT_ReProject;
13背投影 SIRT_BackProject;
14}
15写回第y个切片重构后的图像值X(N)；end
算法2 BackProject_Init背投影赋初值Begin
—declspec(target(mic)) BackProject_Init_kerneI ()
{
1#pragma omp parallel for private (n, x, . . .) num_threads()
2for z =O :SZ-1 //z为像素点在重构图像中的行号
3for n=0: ANGLE_NUM_1 //n代表旋转角度序号
4for X = O :SX-1 //x为像素点在重构图像中的列号
5根据旋转角度和几何参数确定经过像素点(z，X)的光线(n，i)以及A(z，x);
6//(n, i)表示第η个旋转角度的第i条光线，即总第n*SX+i条光线，以下类同
7// (z, X)表示重构图像中第z*SX+x个像素点，以下类同
8Χ(0) ( z, X) += A(z, x)*P测(n, i);
}
end
算法3 ReProject重投影Begin
—declspec (target(mic)) ReProject_kerneI ()
{
1#pragma omp parallel for private (i，zcur，···)num—threads ()
2for n=0: ANGLE—NUM-I //n表示旋转角度序号 3for zcur=0: SZ //循环步长为cos ([η])，[η]表示第η个选择角度值
4for i = 0:SX-1 //i为像素点在第n个角度下投影图像中的序号
5根据旋转角度和参数确定光线(n，i)经过的像素点(zcur，x)和A(zcur，x);
6 P 计(n，i) += A (zcur, x)*X(k—(zcur, x);
}
end
算法4 : BackProject背投影Begin
—declspec (target(mic)) BackProject—kernel ()
{
1#pragma omp parallel for private (n，x, ···) num—threads ()
2for z =0 SZ-1 //z为像素点在重构图像中的行号
3for n=0: ANGLE—NUM-1 //n代表旋转角度序号
4for X = O SX-1 //x为像素点在重构图像中的列号5根据旋转角度和参数确定经过像素点(z，x)的光线(n，i)以及A(z，x);
6 X(k) (z, x) += A(z, x)*(P#J (n, i)—P计(n， i));
}
步骤三负载均衡设计
(1)SIRT算法的任务级并行
三维投影数据是按照每片进行处理的。而每片数据的大小和格式都是相同的，这样在对他们进行处理所需要的时间几乎是一致的。所以无论是在CPU上的每个核上进行计算还是在多MIC卡上计算，按片来分别处理是非常合理的。CPU上的每个核或者每块MIC卡的处理能力一致，这样将每片切片数据分别划分给每个核或者每块MIC卡上，可以很好的利用硬件资源，并行计算，提高处理速度；· (2)SIRT算法的算法级并行
将每片切片数据进行分割，然后分配到不同的CPU核心上或者MIC卡上来同时计算一片数据，根据SIRT算法的分析，数据间有交互，即存在I/O传输，每片数据都处理的话，积累的频繁I/O传输将会大大影响性能。综上所述，两种级别并行的比较，可见采用任务级并行性能更佳。根据上述并行比较的结果，任务级并行中，MIC卡内存空间足够存放每片数据的大小，对数据的动态划分是受到OpenMP线程的控制，同时使用锁操作，避免数据读写的冲突，使用上锁操作，保证每片切片处理的独立性和对切片的读写冲突，各个设备间没有相互依赖。具体优化I :SIRT_MIC优化，实施过程如下
1)减少MICoffload次数
在MIC编程模式中,在offload引语中直接用in O，out O , inout O等引语进行内存和MIC卡上内存之间的数据拷贝，在MIC中频繁的10会严重影响最后的性能，为了减少10，我们把offload语句放到赋初值语句之前，并把所有的迭代过程放到MIC卡上执行，控制流采用一个线程执行，openmp并行之后的语句在MIC上采用多线程执行的方式；
2)自动向量化
在SIRT算法中，赋初值和背投影两个函数的最内层循环语句都是按行访问的方式，通过编译器可以完全自动化，而重投影函数最内层循环并没有按行访问，并且内层for有数据依赖，无法自动向量化，我们对内层的两个for语句做了顺序的调整，调整之后的语句可以完全自动向量化。本发明在提供了上述方法实施例和系统实施例后，根据电镜三维重构的数据规模，选取合适的CPU和KNF MIC卡分别进行SIRT方法的电镜三维重构的运算测试，其中CPU采用了 Intel公司的双路4核Intel Xeon(E) CPU E5630,每个CPU核可以启动两个CPU线程MIC采用Intel公司的KNF MIC，其详细参数见表I。表I测试环境
权利要求
1.一种基于SIRT算法且使用CPU+MIC协同计算模式的电镜三维重构方法，其特征在于包括CPU端，MIC众核协处理器端以及CPU+MIC协同计算模式，其中 (PU端负责将要进行迭代重构图像中的物体分成切片单元，向MIC卡传递重构所需要的值，CPU+MIC协同计算模式的框架搭建以及任务调度和参数初始化工作，而且在整个三维图像重构的计算任务CPU也会以Openmp多线程模式，依次通过重投影算法和背投影算法计算获取切片的重构图像值来动态参与； MIC众核协处理器负责多线程并行地根据切片的重构图像值和测量获取的该切片投影图像值，依次通过重投影算法和背投影算法计算获取切片的重构图像值，在MIC卡上也采用openmp多线程的方式来运算； 重构步骤如下 (PU端将要进行迭代重构图像中的物体分成切片，向MIC卡上传递重构计算所需要的值，以openmp多线程模式，并行执行迭代重构图像，具体包括 根据物体的大小将该物体分成切片，切片的数目为该物体宽度的像素值SY，切片在投影图像上的元素数为N = SX*ANG_NUM，SX为该物体长度的像素值，ANG_NUM为该物体旋转角度的总数；切片在重构图像上的元素数为M = SX*SZ，SZ为该物体高度的像素值； CPU端将为图像数据的角度文件，thita文件在MIC卡上开辟空间，确定其输入输出类型并且将其值offload到MIC卡上，CPU负责CPU+MIC协同计算模式的框架搭建以及任务调度，具体步骤包括 单节点服务器采用由双路6核CPU和2块KNF MIC卡组成的桌面服务器，在CPU+MICs协同计算中，把双路CPU和MIC卡都作为计算设备，这样每个单节点就相当于有3个计算设备，每个设备通过一个OpenMP线程进行控制； 假设单节点上，MIC卡的总数为M个，则节点上共M+1个设备，采用数据动态划分的方式，每次各个设备自动获取切片的输入数据，片数目>=1，数据处理完之后自动获取下一个切片数据，直到所有设备计算完成所有样本切片数据； 在MIC众核协处理器端通过赋O值或随机值或背投影方法获取切片的初始的重构图像值X(O)，然后再通过重投影算法和背投影算法计算获取切片的重构图像值，具体步骤包括 根据赋初值得到的X(O)或上次迭代得到的X (k-Ι)进行重新投影得到一个投影图像值P计；根据赋初值得到的X(O)或上次迭代得到的x(k-l)，测量得到的投影图像值P测和重投影计算得到的投影图像值P计进行反向投影得到重构图像X (k)，背投影主要是进行X值的更新操作，通过多次的重投影和背投影操作可以得到满足精度要求的重构图像； 在图像重构中，测量物体的投影图像像素值P测的过程中，记物体的长为SX，宽为SY，高为SZ，旋转角度的总数为ANG_NUM，垂直于Y方向的面即为切片，共SY个切片，重构的过程以切片为单位，即每次重构一个切片，重构SY个切片就可以重构出物体模型； 分析串行迭代法中的性能瓶颈，找到耗时的重投影和背投影的代码，测试其在整个程序中所占的时间比例； 在每个切片的重构过程中，对于重投影过程，就是根据像素值X求得一个计算得到的投影值P测，P测共有SX*ANG_NUM个元素，这SX*ANG_NUM个元素的求解中没有任何的数据依赖性，以SX*ANG_NUM的并行度进行重投影算法的并行，在内核中，每个线程负责计算P测中的一个元素； 在每个切片的重构过程中，对于背投影过程，也就是求重构图像像素值X的过程，该算法就是根据P计和P测通过背投影的方式求重构图像像素值X，对于每一个切片，X共有SX*SY个元素,这SX*SY个元素的求解中没有任何的数据依赖性，以SX*SY的并行度进行重投影算法的并行，在内核中，每个线程负责计算X中的一个元素。·
全文摘要
本发明提供一种基于SIRT算法且使用CPU+MIC协同计算模式的电镜三维重构方法，CPU端负责将要进行迭代重构图像中的物体分成切片单元，向MIC卡传递重构所需要的值，CPU+MIC协同计算模式的框架搭建以及任务调度和参数初始化工作，而且在整个三维图像重构的计算任务CPU也会以Openmp多线程模式，依次通过重投影算法和背投影算法计算获取切片的重构图像值来动态参与；MIC众核协处理器负责多线程并行地根据切片的重构图像值和测量获取的该切片投影图像值，依次通过重投影算法和背投影算法计算获取切片的重构图像值，在MIC卡上也采用openmp多线程的方式来运算。
文档编号G06T1/20GK102945566SQ201210398979
公开日2013年2月27日 申请日期2012年10月19日 优先权日2012年10月19日
发明者卢晓伟, 张清 申请人:浪潮电子信息产业股份有限公司