CN102945566A - 一种基于sirt算法且使用cpu+mic协同计算模式的电镜三维重构方法 - Google Patents

Abstract

本发明提供一种基于SIRT算法且使用CPU+MIC协同计算模式的电镜三维重构方法，CPU端：负责将要进行迭代重构图像中的物体分成切片单元，向MIC卡传递重构所需要的值，CPU+MIC协同计算模式的框架搭建以及任务调度和参数初始化工作，而且在整个三维图像重构的计算任务CPU也会以Openmp多线程模式，依次通过重投影算法和背投影算法计算获取切片的重构图像值来动态参与；MIC众核协处理器：负责多线程并行地根据切片的重构图像值和测量获取的该切片投影图像值，依次通过重投影算法和背投影算法计算获取切片的重构图像值，在MIC卡上也采用openmp多线程的方式来运算。

Description

一种基于SIRT算法且使用CPU+MIC协同计算模式的电镜三维重构方法

技术领域

本发明涉及计算机应用技术领域, 具体地说是一种基于SIRT算法且使用CPU+MIC协同计算模式的电镜三维图像重构方法，尤其涉及在MIC众核架构上使用CPU+MIC协同计算模式进行电镜三维重构的方法及装置。 

背景技术

MIC(Many Integrated Core)是Intel公司推出的众核处理器，跟通用的多核至强处理器相比，MIC众核架构具有更小的内核和硬件线程，众核处理器计算资源密度更高，片上通信开销显著降低，更多的晶体管和能量，能够胜任更为复杂的并行应用。Intel MIC产品基于X86架构，基于重核的众核处理器，包含50个以上的核心，以及512bit的向量位宽，双精性能超过1TFlops。 

MIC拥有极其灵活的编程方式，MIC卡可以作为一个协处理器存在，也可以被看作是一个独立的节点。基本的MIC编程模型是将MIC看作一个协处理器，CPU根据程序的指令，将一部分代码运行在MIC端。此时存在两类设备，即CPU端和MIC众核协处理器端。 

电子断层三维重构技术(Electron Tomography，ET)是从一个物体的二维投影图经过重构获取物体内部结构的技术，通过获取多个不同角度的二维投影图进行反向重构出所研究对象的三维结构。图像重构无论是在医学诊断，还是在生物样品的成像，以及在地表层析成像等诸多领域中都有着广泛的应用。迭代法可以表述成：Ax=p，其中A为非奇异矩阵，x为未知数向量，p为已知向量，迭代法求解目标在于确定x的数值。在迭代重构图像方法中，A为加权因子，p为投影得到的图像像素值，x则为需要重构的图像的像素值。 

图1说明了采用迭代重构图像方法实现图像重构的过程。该方法通过赋0值或随机值或背投影方法获取用于迭代过程中的初值X⁽⁰⁾，然后经过多次的重投影和背投影的迭代过程，来求得最后满足精度要求的重构图像的像素值。 

SIRT（联合迭代重构）算法在所有图像重构算法中重构效果较好，但其运算量巨大，因此需要提高运算效率，来满足快速图像重构的需求。可见，目前需要针对基于MIC平台设计一种迭代重构图像的方法，通过MIC高效简单快速的计算特性，能够极大地加速迭代重构方法的执行，从而快速地完成图像重构任务。 

发明内容

本发明的目的是提供一种使用CPU+MIC协同计算模式的电镜三维图像重构方法。能够基于MIC众核架构并且使用CPU+MIC协同计算的快速图像重构。 

本发明的目的是按以下方式实现的，包括CPU端， MIC众核协处理器端以及CPU+MIC协同计算模式，其中： 

CPU端：负责将要进行迭代重构图像中的物体分成切片单元，向MIC卡传递重构所需要的值，CPU+MIC协同计算模式的框架搭建以及任务调度和参数初始化工作，而且在整个三维图像重构的计算任务CPU也会以Openmp多线程模式，依次通过重投影算法和背投影算法计算获取切片的重构图像值来动态参与；

MIC众核协处理器：负责多线程并行地根据切片的重构图像值和测量获取的该切片投影图像值，依次通过重投影算法和背投影算法计算获取切片的重构图像值，在MIC卡上也采用openmp多线程的方式来运算； 

重构步骤如下：

CPU端：将要进行迭代重构图像中的物体分成切片，向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 众核协处理器端：通过赋0值或随机值或背投影方法获取切片的初始的重构图像值X(0)，然后再通过重投影算法和背投影算法计算获取切片的重构图像值，具体步骤包括：

根据赋初值得到的X(0)或上次迭代得到的X(k-1)进行重新投影得到一个投影图像值P计；根据赋初值得到的X(0)或上次迭代得到的X(k-1)，测量得到的投影图像值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中的一个元素。

附图说明

图1 是采用迭代重构图像方法实现图像重构的流程图。 

具体实施方式

以下结合附图和优选实施例对本发明的技术方案进行详细地阐述。以下例举的实施例仅用于说明和解释本发明，而不构成对本发明技术方案的限制。 

本发明首先分析用SIRT串行方法进行电镜三维图像重构的性能瓶颈以及将串行方法移植到其他高性能平台上的难度，找到耗时的热点代码，测试其在电镜三维图像重构的整个过程中所占的时间比例以及分析提高整个方法性能的难度以及开发轴端。 

测试结果表明在采用SIRT联合迭代法进行图像重构时，大部分的时间均消耗在重投影和背投影两步骤中，这两个步骤不仅执行次数多，而且要计算的数据量大，同时，数据的计算要满足单指令多数据（SIMD，Single Instruction Multiple Data）的特点。因此，本发明将这两步骤作为加入MIC众核协处理器进行并行改进的重点，同时由于CPU端除了进行任务调度等功工作之余，也参与赋初值，再通过重投影算法和背投影算法计算获取切片的重构图像值的工作，即CPU+MIC协同计算的工作模式。。 

在图像重构中，测量获取物体的投影图像像素值P_测的过程，记物体的长为SX像素，宽为SY像素，高为SZ像素，旋转角度的总数（即旋转次数）为ANG_NUM，垂直于Y方向的每一个面为切片，共有SY个切片。图像重构的过程是以切片为单位进行的，其过程如图1所示，即对每一切片通过迭代多次计算和测量获取切片的投影图像像素值来重构；重构了SY个切片便重构出该物体模型。 

实施例

本发明针对上述分析，提供了一种基于SIRT算法且使用CPU+MIC协同计算模式的电镜三维重构并行方法的实施例，包括如下步骤： 

步骤一：CPU+MIC协同计算框架搭建

在单节点上，共有M+1个设备（一个CPU+M个MIC卡），采用OpenMP的fork-join模式搭建单节点上的框架，当程序开始执行的时候只有一个主线程存在，需要进行并行计算的时候，主线程派生出附加线程，即启用M+1个OpenMP线程， 通过while(1) 循环各个主动进行输入数据的分发读写，0~M-1号进程控制所有MIC设备，M号线程控制CPU设备。

如在本专利中，CPU+2MICs平台上， 主线程控制输入数据的动态分发，0号线程控制MIC0设备，1号线程控制MIC1设备，2号线程控制CPU设备。 

在SIRT算法中，MIC卡上的内存空间足够存放切片数据的大小，所以不需要采取分块处理的方式，每个设备重构完一个切片数据之后，输出重构图像，在每次读取当前切片的时候，需进行上锁操作，保证每片切片处理的独立性和对切片的读写冲突，各个设备间没有相互依赖。 

单节点上CPU+MICs协同计算的外围框架的伪代码如下所示： 

//SY为切片数目

int DEVICE_NUM=3  //为CPU和MIC卡的总个数

int Slice_private=0;    //为每个线程定义一个私有的切片变量

int Slice=0;          //定义一个共享变量，来控制切片的读写

omp_lock_t  lck;

#pragma omp parallel for private (Slice_private) ,num_threads (DEVICE_NUM)

for (thread=0;  thread<DEVICE_NUM;  thread++)

{

while(1) //死循环控制三个线程不停地进行工作

   { 

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))

          BackProject_Init_kernel();  //背投影赋初值X⁽⁰⁾

          for: k=1: N                 //N代表要迭代的次数

            ReProject_kernel();    //重投影

            BackProject_kernel();  //背投影

         写回第slice_private个切片重构后的图像值X^(N)；

  }

else if(thread==DEVICE_NUM-1) //cpu compute

{

            BackProject_Init_kernel();  //背投影赋初值X⁽⁰⁾

          for: k=1: N                 //N代表要迭代的次数

            ReProject_kernel();    //重投影

            BackProject_kernel();  //背投影

         写回第slice_private个切片重构后的图像值X^(N)；

}

else  //other mic compute

{

          #pragma offload target (mic:1)\

          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))

          BackProject_Init_kernel();  //背投影赋初值X⁽⁰⁾

          for: k=1: N                 //N代表要迭代的次数

            ReProject_kernel();    //重投影

            BackProject_kernel();  //背投影

         写回第slice_private个切片重构后的图像值X^(N)；

}

   }  // end while

步骤二：CPU/MIC内核实现

在CPU+MICs协同计算平台上，每个节点上有双路6核CPU，总共12个核，如果超线程打开的话可以达到24个线程，在CPU内核设计中可以采用OpenMP开启多线程加速计算。在MIC内核设计中，根据MIC卡的核数目设计OpenMP的并行线程数。CPU和MIC采用相同的内核计算，伪代码如下：

算法1：SIRT_MIC整体实现

begin

1  for: y=0:SY-1 //y代表每一个切片

2      读取第y个切片的投影图像值;若运行在MIC卡上，则offload

3  #pragma offload target \

4  in (SX, SZ, ANGLE_NUM, N, …) \

5  in (P_测: length(ANGLE_NUM * SX)) \

6  in (thita: length(ANGLE_NUM)) \

7  in (P_计: length(ANGLE_NUM * SX)) \

8  inout (X: length(SZ*SX))

 否则

9

{

10     背投影赋初值X⁽⁰⁾;

11     for: k=1: N //N代表要迭代的次数

12        重投影SIRT_ReProject; 

13        背投影SIRT_BackProject; 

14  }

15     写回第y个切片重构后的图像值X^(N)；

end

算法2：BackProject_Init 背投影赋初值

Begin

__declspec(target(mic)) BackProject_Init_kernel()

{

1  #pragma omp parallel for private (n, x, ...) num_threads()

2  for  z =0：SZ-1 //z为像素点在重构图像中的行号

3      for  n=0: ANGLE_NUM-1 //n代表旋转角度序号

4  for  x = 0：SX-1 //x为像素点在重构图像中的列号

5  根据旋转角度和几何参数确定经过像素点(z, x)的光线 (n, i)以及A(z, x);

6  //(n, i)表示第n个旋转角度的第i条光线，即总第n*SX+i条光线，以下类同

7  // (z, x)表示重构图像中第z*SX+x个像素点，以下类同

8  X⁽⁰⁾( z, x)  +=  A(z, x)*P_测(n, i);

} 

end

算法3： ReProject 重投影

Begin

__declspec(target(mic))  ReProject_kernel()

{

1  #pragma omp parallel for private (i, zcur, …) num_threads()

2  for  n=0:ANGLE_NUM-1 //n表示旋转角度序号

3  for zcur=0: SZ  //循环步长为cos([n]) ,[n]表示第n个选择角度值

4  for  i = 0：SX-1 //i为像素点在第n个角度下投影图像中的序号

5  根据旋转角度和参数确定光线(n,i)经过的像素点(zcur, x)和A(zcur,x); 

6  P_计(n, i) += A(zcur, x)*X^(k-1) (zcur, x);

}

end

算法4： BackProject 背投影

Begin

__declspec(target(mic))  BackProject_kernel()

{

1  #pragma omp parallel for private (n, x, …) num_threads()

2  for  z =0：SZ-1 //z为像素点在重构图像中的行号

3      for  n=0: ANGLE_NUM-1 //n代表旋转角度序号

4         for  x = 0：SX-1 //x为像素点在重构图像中的列号

5  根据旋转角度和参数确定经过像素点(z, x)的光线(n, i)以及A(z, x);

6  X^(k) (z, x) += A(z, x)*(P_测(n, i)-P_计(n, i));

} 

步骤三：负载均衡设计

（1）SIRT算法的任务级并行

三维投影数据是按照每片进行处理的。而每片数据的大小和格式都是相同的，这样在对他们进行处理所需要的时间几乎是一致的。所以无论是在CPU上的每个核上进行计算还是在多MIC卡上计算，按片来分别处理是非常合理的。CPU上的每个核或者每块MIC卡的处理能力一致，这样将每片切片数据分别划分给每个核或者每块MIC卡上，可以很好的利用硬件资源，并行计算，提高处理速度；

（2）SIRT算法的算法级并行

将每片切片数据进行分割，然后分配到不同的CPU核心上或者MIC卡上来同时计算一片数据，根据SIRT算法的分析，数据间有交互，即存在I/O传输，每片数据都处理的话，积累的频繁I/O传输将会大大影响性能。

综上所述，两种级别并行的比较，可见采用任务级并行性能更佳。 

根据上述并行比较的结果，任务级并行中，MIC卡内存空间足够存放每片数据的大小，对数据的动态划分是受到OpenMP线程的控制，同时使用锁操作，避免数据读写的冲突，使用上锁操作，保证每片切片处理的独立性和对切片的读写冲突，各个设备间没有相互依赖。 

具体优化1：SIRT_MIC优化，实施过程如下： 

1）减少MIC offload次数

在MIC编程模式中，在offload引语中直接用in(), out(), inout()等引语进行内存和MIC卡上内存之间的数据拷贝，在MIC中频繁的IO会严重影响最后的性能，为了减少IO，我们把offload语句放到赋初值语句之前，并把所有的迭代过程放到MIC卡上执行，控制流采用一个线程执行，openmp并行之后的语句在MIC上采用多线程执行的方式；

2）自动向量化

在SIRT算法中，赋初值和背投影两个函数的最内层循环语句都是按行访问的方式，通过编译器可以完全自动化，而重投影函数最内层循环并没有按行访问，并且内层for有数据依赖，无法自动向量化，我们对内层的两个for语句做了顺序的调整，调整之后的语句可以完全自动向量化。

本发明在提供了上述方法实施例和系统实施例后，根据电镜三维重构的数据规模，选取合适的CPU和KNF MIC卡分别进行SIRT方法的电镜三维重构的运算测试，其中CPU采用了Intel 公司的双路4核Intel^(R) Xeon^(R) CPU E5630,每个CPU核可以启动两个CPU线程MIC采用Intel 公司的KNF MIC，其详细参数见表1。 

表1 测试环境 

CPU	Intel(R) Xeon(R) CPU E5630 2.53GHz，8核16线程
		Memory	DDR3 1066Mhz 24GB
MIC	Total No of Cores：30，No of Active Cores：30，1.05Ghz，MIC Board Stepping：C0，GDDR5 1500 MHz 2GB
		OS	Red Hat Enterprise Linux Server release 6.0，Linux version 2.6.32-71.el6.x86_64
编译器	Icc

表2 测试用例

ST图像	大小(像素)	旋转角度	样品厚度(像素)	角度文件
					test4a	2048*2048	112	400	从-60度开始每隔一度取一个角度

表3 实验结果

我们通过在MIC的实验机器上，测得test4a的性能如表3.3所示。对于串行程序，MIC程序的整体性能达到25.78倍，其中背投影函数的性能最好，性能超过32倍，而重投影只有22倍的性能。

表3.3 实验结果 

注：bp_init表示背投影赋初值，Rp表示重投影，Bp表示背投影，all表示

由本发明给出的上述测试结果可以看出，在MIC众核架构平台上，利用SIRT方式和MIC众核协处理器强大的并行计算能力，可以在较短周期内，较容易地大幅度的加速三维电镜重构的运算。

以上说明仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权力要求书的保护范围为准。 

除说明书所述的技术特征外，均为本专业技术人员的已知技术。 

Claims (1)

1.一种基于SIRT算法且使用CPU+MIC协同计算模式的电镜三维重构方法，其特征在于包括CPU端， MIC众核协处理器端以及CPU+MIC协同计算模式，其中：

CPU端：负责将要进行迭代重构图像中的物体分成切片单元，向MIC卡传递重构所需要的值，CPU+MIC协同计算模式的框架搭建以及任务调度和参数初始化工作，而且在整个三维图像重构的计算任务CPU也会以Openmp多线程模式，依次通过重投影算法和背投影算法计算获取切片的重构图像值来动态参与；

MIC众核协处理器：负责多线程并行地根据切片的重构图像值和测量获取的该切片投影图像值，依次通过重投影算法和背投影算法计算获取切片的重构图像值，在MIC卡上也采用openmp多线程的方式来运算；    

重构步骤如下：

CPU端：将要进行迭代重构图像中的物体分成切片，向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 众核协处理器端：通过赋0值或随机值或背投影方法获取切片的初始的重构图像值X(0)，然后再通过重投影算法和背投影算法计算获取切片的重构图像值，具体步骤包括：

根据赋初值得到的X(0)或上次迭代得到的X(k-1)进行重新投影得到一个投影图像值P计；根据赋初值得到的X(0)或上次迭代得到的X(k-1)，测量得到的投影图像值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中的一个元素。

Images