2018, Vol. 35 
2. 上海卡勒幅磁共振技术有限公司, 上海 201614;
3. 上海纽迈电子科技有限公司, 上海 200333
2. Shanghai Colorful Magnetic Resonance Technology Corporation Limited, Shanghai 201614, China;
3. Shanghai Niumag Electronic Technology Corporation, Shanghai 200333, China
磁共振成像（Magnetic Resonance Imaging, MRI）系统中的数据处理流程包含很多步骤，如图像重建、后处理、可视化、分析、存储、网络传输等.MRI实验方法丰富多样，相应的图像重建与后处理工作内容也随之变化.一些特殊的磁共振扫描方法需要相应的特殊图像重建、后处理和可视化方法.例如，为了加速MRI扫描速度所采用的部分傅里叶[1]、匙孔成像[2]、螺旋式采集[3]、径向采集[4]、并行采集[5-8]以及近年来新出现的压缩感知（Compressed Sensing, CS）技术[9-12]，都要求在图像重建过程中引入特殊的步骤.而功能成像[13]、波谱成像[14]、扩散张量成像[15]等许多高级MRI方法，也对后处理及可视化方法提出了特殊的要求.随着MRI设备的不断普及以及系统软硬件技术的不断提升，对磁共振数据处理提出了新的挑战.例如随着MRI系统通道数的不断增加，重建数据的流量亦显著增加，这对并行及分布式的数据处理提出了要求；另一方面，随着影像数据的增加，在几何形状分析、纹理分析与直方图分析的基础上发展起来的放射组学（Radiomics）[16, 17]和基于深度学习的计算机辅助诊断（Computer Aided Diagnosis, CAD）系统[18]近年来取得了飞速发展，而这些都是以大量的数据处理算法研究为基础的.
在这种情况下，为上述磁共振图像数据处理提供一个通用的处理框架颇显实用价值：（1）通用的基础框架，以及常用的图像处理、存储、网络传输、可视化等工具，为算法研究者实现并测试算法提供便利；（2）开放的接口，使算法研究者可以方便地分享自己的研究成果；（3）通过提供与商用系统连接的基础设施，可以方便算法研究者在实际的磁共振系统中检验自己的算法，也方便算法研究成果的应用推广[19].
在这种背景下，有研究者提出了多种针对不同用途的磁共振图像处理框架，这些框架将磁共振图像处理过程中的不同步骤都做成独立的处理模块，根据处理任务的不同，选择相应的处理模块并将它们串接成一条流水线，完成特定的数据加工任务.比较知名的开源磁共振图像处理流水线包括：（1）2006年，Cook等[20]提出的DMRI（Diffusion Magnetic Resonance Imaging）数据处理流水线——Camino；（2）2012年，Daducci等[21]提出的用Python编写的DMRI数据处理流水线——The Connectome Mapper；（3）2015年，Uecker等[22]提出的医学图像重建工具箱——BART（Berkeley Advanced Reconstruction Toolbox），其中包含了一些比较前沿的迭代优化算法，它可以作为自己算法的开发平台，也可以借助该平台将自己的算法整合到临床环境中去；（4）2013年，Hansen等[19]提出的Gadgetron是目前最具代表性的开源医学图像重建框架，它不仅可以跨平台使用，还可以与磁共振成像系统连接，在线进行数据处理.
本文开发了一个新的医学图像处理流水线框架，其设计目标和结构与上述流水线有很多共同之处，因此将它命名为“又一条流水线”（Yet Another Pipeline，YAP，https://github.com/MR-algorithms/YAP）.之所以要在现存这么多开源流水线的情况下开发自己的流水线框架，是因为现有的流水线无法满足我们在商用的MRI系统中使用流水线的需要，如：（1）不支持分支结构；（2）处理单元的参数配置方法不灵活，一个流水线中同一个处理模块具有多个实例时，配置较为困难；（3）支持的数据形式不够灵活.而且，我们希望将流水线的功能扩展到MRI系统中的数据采集步骤，这需要更灵活的参数管理机制，以便流水线能够自然地与扫描参数设置步骤衔接.
YAP框架采用了面向接口的设计，具有良好的可扩展性，用户可以独立开发流水线中的处理模块.由于不同版本的C++编译器在处理复合数据结构时采用的内存映射可能会有所不同，因此如果开发不同模块的开发人员使用不同版本的C++编译器，则在这些模块之间就无法传递复合数据.很多商用MRI系统的开发平台，都要求进行序列与重建算法开发所使用的C++编译器的版本必须与MRI系统开发时所用的C++编译器相同.由于系统开发的时间通常比较早，往往使用商用的C++编译器，这就给喜欢使用最新的开源C++编译器的算法研究人员带来了很大的不便.而YAP采用基于接口的设计，由于C++的接口在二进制上兼容，因此算法研究人员可以用任意的C++编译器进行插件开发.
本文首先介绍了YAP基本框架，然后介绍利用该框架完成的几个实际应用，验证了YAP框架的可用性.
1 YAP流水线的结构与实现 1.1 YAP流水线框架结构与处理器YAP框架结构中包含了几个功能不同的模块：
(1) Interfaces中定义了YAP中使用的基本接口，所有插件的功能、流水线中的数据和相关信息的访问都需要通过这些接口进行；
(2) Implement模块为这些接口提供了基本的实现类，使得插件的开发人员可以专注于算法实现的细节，而不是接口实现的细节；
(3) BasicRecon项目中实现了一些基本的磁共振图像数据处理器；
(4) YAP模块则提供了使用YAP流水线的应用程序编程接口(Application Programming Interface，API)；
(5) PipelineTest项目演示了YAP API的基本用法.
MRI数据处理流程中往往包含很多步骤，图 1给出了一条典型的磁共振图像处理流水线，包括原始数据读取、k空间处理（如直流分量去除、填零）、傅里叶变换、取模和保存最后的图像文件.数据处理过程中每个步骤的功能相对独立，可以单独实现，本文中将它们称为数据处理器（Processor）.
|
图 1 简单磁共振图像处理流水线 Figure 1 A basic pipeline for MR image processing |
每个处理器都可以接收特定格式的数据进行处理，并将结果输出到下一个处理器.在YAP框架中，每个处理器可以包含多个输入和输出端口（Port），如图 2中的ModulePhase就包括两个输出端口，Module输出模图，Phase输出相位图；如果处理器只有一个输入或者输出端口，也可以省略名称，缺省的输入端口名称为Input，缺省的输出端口名称为Output.
|
图 2 处理器的基本结构 Figure 2 The structure of processor |
此外，每个处理器还可以具有一些可以配置的属性.如图 2中的FFT处理器通过属性Inverse来确定傅里叶变换的方向，如果Inverse为false，则处理器对输入数据进行反傅里叶变换，反之，进行傅里叶变换；FFT另一个属性InPlace则指定了傅里叶变换的结果是否直接保存在原数组中.
从图 1中还可以看出，YAP支持分支结构，ModulePhase处理器的Module输出，被同时送入JpegExporter和2DDisplay两个模块，即保存图像文件的同时，在2D显示窗口中显示图像.
YAP中的处理器，都需要实现IProcessor接口，通常在一个独立的动态链接文件（DLL文件）中实现若干个处理器，在构建流水线时，可以动态加载这些动态链接文件并使用其中的处理器.
1.2 YAP流水线的构建与运行为了能够简单而快速地构建流水线，本文设计了一种定义流水线的脚本文件（.pipeline文件），可以在脚本文件中编辑流水线及其参数（处理器的属性），而不需要重新编译YAP项目.与图 1中的流水线对应的脚本文件内容如下：
|
基于上面的脚本文件，用户很容易就可以构建图像处理流水线.为了使用户更熟悉流水线的构建，下面介绍几个要点：
|
import语句用于加载包含处理器的DLL文件，在后续的脚本中就可以使用DLL中包含的处理器.
|
上述语句定义了一个类型为CmrRawDataReader的处理器的实例reader，该处理器位于BasicRecon.dll中，同时设置了处理器的DataPath、ChannelCount和ChannelSwitch等属性.属性设置采用有名参数的方式，因此对属性出现的次序没有要求，也不需要指定所有的属性.
|
上述语句中操作符“-＞”用于将端口相连，箭头指向表示数据流动的方向，在使用缺省端口的情况下，可以省略端口.如上述语句，将reader处理器的Output端口的输出送入slice_iterator的Input端口，将module_phase的Module输出端口的输出送入jpeg_exporter的Input输入端口.
|
上面语句中的self表示当前定义的流水线，该语句表示将当前流水线的Input端口设置为reader处理器的Input端口.
有了脚本文件，就可以利用脚本解析器对脚本进行解析并构建流水线，构建完成后，同样可以调用流水线的Input()函数来输入数据，启动流水线运行.完成上述任务的C++代码如下：
|
其中PipelineCompiler根据流水线脚本文件的内容构建流水线.如果流水线脚本文件语法定义正确，PipelineCompiler的CompileFile()函数就可以返回构建好的YAP流水线.如果流水线脚本文件中存在语法错误，处理过程中会抛出带有错误信息的异常，使用PipelineCompiler类的应用程序可以捕获抛出的异常并定位错误在脚本文件中的位置.pipeline的Input()函数通过向流水线中馈送数据驱动流水线工作，在上述脚本中，因为已经通过处理器属性指定了要处理的原始数据文件路径，所以Input()函数并不需要指定数据.
1.3 处理器的创建如前所述，为了最大限度的满足扩展性的要求，YAP采用了基于接口的设计.算法研究人员可以实现IProcessor接口，将处理算法做成YAP中可用的处理器；而构建流水线时，只需要包含处理器的DLL文件，而不需要处理器的源代码.这更便于不愿意开放算法源代码的算法研究人员分享并测试自己的算法.由于不同的编译器处理数据结构时采用的内存映像可能不同，因此YAP流水线中跨模块传递的对象必须通过接口实现，并且接口中只能包括基本数据类型；此外，某个模块分配的内存，只能由该模块自己释放，这也意味着这些对象只能在自由存储（堆）上进行分配.为了保证这一点，YAP中跨模块传递的所有对象均要实现ISharedObject接口：
|
在使用对象之前，首先使用Lock()函数锁定对象，而在结束使用之后，调用Release()函数释放对象，对象可以根据自己的使用情况选择合适的方式释放自己.由于YAP中跨模块传递的对象都需要创建自己的拷贝，所以也要求这些对象使用Clone()函数，返回自己的一个拷贝.
IProcessor接口派生于ISharedObject接口.为了方便用户实现IProcessor接口，YAP提供了ProcessorImpl类，它实现了IProcessor接口中的共有功能.创建处理器通常要创建一个DLL项目，然后按以下步骤进行：
(1) 在项目中添加一个ProcessorImpl的派生类，如Fft2D；
(2) 在类的构造函数中添加处理器的属性，并给属性设置缺省值：
|
(3) 在类构造函数中添加输入、输出端口：
|
其中AddInput()函数添加了名为“Input”的输入端口，第二个参数表示该端口上可以接收二维数据，第三个参数表示端口可接收的数据类型是基于float的复数.而AddOutput()函数则添加了名为“Output”的输出端口.
(4) 实现处理器的Input()方法.Input()函数用于向处理器输入数据，由于流水线是数据驱动的，所以在Input()函数中应对输入数据进行处理并将结果输出到输出端口（这会间接调用与输出端口相连的下一个处理器的Input()函数）.
(5) 在DLL项目的dllmain.cpp文件里面添加如下语句：
|
其中包含的3个宏，是YAP定义的，帮助用户在DLL中创建GetProcessorManager()工厂函数，并作为DLL的导出函数.GetProcessorManager()返回IProcessorContainer接口，利用该接口，可以遍历DLL中的处理器.
1.4 YAP流水线中的数据YAP流水线是数据驱动的，启动流水线或者任意一个处理器工作的方式就是向流水线或处理器输入数据.MRI系统中要处理的数据，通常是一个高维数组，包括频率编码、相位编码、层、通道、回波、累加等维度；根据实验设计，也可以包括扩散方向、时间等更多维度；数据本身也可能是多种类型，例如整数、浮点数或者复数；此外，磁共振图像的描述还需要相关的实验参数，包括空间位置信息.为此YAP定义了用于描述数据的IData接口：
|
IData接口从ISharedObject派生，其中GetDataType()返回代表数据元素类型的整型常量，这些类型包括了DataTypeChar、DataTypeShort、DataTypeDouble等.GetDimensions()函数返回用于描述数据维度信息的IDimensions接口的指针，GetGeometry()函数返回用于描述图像空间定位信息的IGeometry接口的指针，GetVariables()返回用于访问与数据相关联的参数的IVariableContainer接口的指针.限于篇幅，本文就不对这些接口进行描述了.
在YAP流水线中实际传输的数据对象需要实现IData接口派生的模板接口IDataArray＜＞.它的结构如下：
|
其中GetData()函数用于返回实际包含数据的C语言数组.这个数组不能使用delete操作符释放，只能使用数据对象的Release()函数释放.YAP允许构建流水线对不同类型、不同维度的数据进行处理，而组成流水线的处理器能够处理的数据类型和维度也不尽相同，利用处理器构建流水线时，将处理器的输出端口与后续处理器的输入端口相连，前提是要连接的两个端口必须支持兼容类型和维度的数据.
2 YAP流水线的应用 2.1 CS在流水线中的实现CS理论打破奈奎斯特采样定理的限制，只要随机采集k空间的一个子集（2D成像中沿着相位编码方向）就能重建出质量较好的磁共振图像，从而加快MRI扫描速度.CS重建过程实际上是一个优化过程，获得最优解的方法有很多，非线性共轭梯度法[23]就是其中一种.由于非线性共轭梯度法是一个耗时迭代过程，我们基于YAP实现了一个图形处理器（Graphics Processing Unit, GPU）加速版的CS重建算法研究的流水线.在保证重建图像质量的同时，大大减少了重建时间.我们以此为例说明应用YAP的基本过程.
2.1.1 CS处理器的创建我们实现的用于CS算法研究的流水线如图 3(a)所示，图 3(b)为重建结果图，图像大小为512×512，采样率为35%.其中MaskCreator、Subsampling和CompressedSensing等几个处理器是该流水线特有的，所以首先需要创建这几个处理器.
|
图 3 (a)CS重建流水线；(b)使用CS算法重建的磁共振图像，采样率为35% Figure 3 (a) CS reconstruction pipeline; (b) The reconstructed MR image using CS method, with the sampling ratio of 35% |
为此，新建一个DLL项目CompressedSensing_CUDA，并添加处理器类CompressedSensing，这个类使用ProcessorImpl作为基类.除了需要编写CS重建算法的代码之外，最主要的工作集中在Input()函数中：
|
其中CsRecon类实现了基于统一计算设备架构（Compute Unified Device Architecture, CUDA）的CS算法，在处理器内部使用.Input()函数在流水线向处理器输入数据时调用，而CompressedSensing处理器有两个输入端口，Mask用于接收采样模板，Input用于接收实际扫描采集的k空间数据.Input()函数接收到采样模板时，将其保存下来；而接收到k空间数据时，则调用CsRecon::GetReconData()函数，利用GPU完成CS计算，最后调用Feed()函数包含重建结果的数组打包成数据对象output（实现了IData接口），并将output对象馈送给下一个处理器.
MaskCreator和Subsampling处理器的实现方式类似.最后在DLL项目中，用ADD_PROCESSOR宏将3个处理器注册到工厂对象中.
2.1.2 CS流水线的构建CS重建流水线同样可以使用脚本来快速构建：
|
与前一条常规流水线相比，有几点需要注意的地方：
(1) 这条流水线只处理某一层数据；
(2) 导入了额外的处理器模块CompressedSensing_CUDA.dll；
(3) mask_creator的输出，被送入sub_sampling和compressedsensing两个处理器，这是分支结构带来的便利；
(4) compressedsensing具有两个输入端口，所以mask_creator和sub_sampling与之连接时需要显式地指定端口名称；
(5) SamplingMaskCreator与CompressedSensing构造时，使用有名参数使得脚本几乎是自说明的，也减少了进行算法研究时出错的可能性.
2.2 补充说明 2.2.1 应用程序内的处理器为了使用方便，YAP流水线的处理器也可以不放在DLL中，而直接放在应用程序中.此时，需要在进行流水线构建之前，调用以下函数来注册处理器：
|
ImportedProcessorManager类用于管理导入的所有处理器，它也会在YAP流水线需要使用处理器时查找从应用程序内注册的处理器.这个功能对于那些与应用程序交互密切的处理器比较方便，如图像显示或者将流水线中的数据流与应用程序中其他流程连接的处理器.在我们使用YAP流水线构建MRI系统软件时，就使用了这种方法.
2.2.2 分支结构对调试的支持YAP支持流水线中使用同一处理器的多个实例，也支持具有分支结构的流水线，这不仅方便了流水线的构建，也方便了流水线的调试.例如，在上述CS流水线中，如果希望观察欠采的数据，可以对流水线结构进行局部的修改，如图 4(a)所示，而图 4(b)是对原始数据欠采填0，再进行傅里叶变换后的重建图，采样率为35%：
|
图 4 (a) 在CS流水线中增加调试处理器观察中间结果；(b)调试处理器输出的欠采填0重建磁共振图像，采样率为35% Figure 4 (a) Observer the intermediate output of compressed sensing pipeline by branching; (b) MR image output by debugging processor, which is reconstructed from an under-sampled k-space data set with sample ratio of 35% |
为此，只需要在前述脚本中增加以下语句：
|
目前，我们的开源框架已经应用到某公司的MRI分析软件中，它的用户界面如图 5所示，该软件通过控制谱仪进行扫描，然后利用YAP流水线重建出图像，并显示在用户界面上.
|
图 5 应用YAP流水线构建的MRI分析软件界面 Figure 5 User interface of an MRI software using YAP |
本文设计并实现了一个新的、开源的医学图像重建与后处理的流水线YAP，主要特点包括：（1）采用基于接口的设计，可使用插件进行功能扩展，插件无需开放源代码，对编译插件的编译器也没有版本限制；（2）允许用户编写脚本来构建图像处理流水线，编辑与修改都很方便；（3）支持带分支结构的流水线，方便了流水线的构建与调试.
然而，本文设计的框架还处于比较初级的阶段，以后会逐渐完善该框架.接下来的工作内容包括：（1）为框架提供一个图形化的用户界面；（2）开发模拟数据采集的模块，对真实扫描过程进行模拟；（3）对医学影像标准DICOM进行支持；（4）实现流水线与Gadgetron流水线及西门子MRI系统互联；（5）完善对数据库的支持.
| [1] | MCGIBNEY G, SMITH M R, NICHOLS S T, et al. Quantitative evaluation of several partial fourier reconstruction algorithms used in MRI[J]. Magn Reson Med, 1993, 30(1): 51-59. |
| [2] | DUERK J L, LEWIN J S, WU D H. Application of keyhole imaging to interventional MRI:A simulation study to predict sequence requirements[J]. J Magn Reson Imaging, 1996, 6(6): 918-924. |
| [3] | LIAO J R. Real-time image reconstruction for spiral MRI using fixed-point calculation[J]. IEEE T Med Imaging, 2000, 19(7): 690-698. |
| [4] | BLOCK K T, UECKER M, FRAHM J. Undersampled radial MRI with multiple coils. Iterative image reconstruction using a total variation constraint[J]. Magn Reson Med, 2007, 57(6): 1086-1098. |
| [5] | PRUESSMANN K P, WEIGER M, SCHEIDEGGER M B, et al. SENSE:sensitivity encoding for fast MRI[J]. Magn Reson Med, 1999, 42(5): 952-962. |
| [6] | GRISWOLD M A, JAKOB P M, HEIDEMANN R M, et al. Generalized autocalibrating partially parallel acquisitions(GRAPPA)[J]. Magn Reson Med, 2002, 47(6): 1202-1210. |
| [7] | DESHMANE A, GULANI V, GRISWOLD M A, et al. Parallel MR imaging[J]. J Magn Reson Imaging, 2012, 36(1): 55-72. DOI: 10.1002/jmri.23639. |
| [8] | BLAIMER M, BREUER F, MUELLER M, et al. SMASH, SENSE, PILS, GRAPPA:how to choose the optimal method[J]. Top Magn Reson Imaging, 2004, 15(4): 223-236. DOI: 10.1097/01.rmr.0000136558.09801.dd. |
| [9] | CANDÈS E J, ROMBERG J, TAO T. Robust uncertainty principles:Exact signal reconstruction from highly incomplete frequency information[J]. IEEE T Inform Theory, 2006, 52(2): 489-509. DOI: 10.1109/TIT.2005.862083. |
| [10] | LUSTIG M, DONOHO D, PAULY J M. Sparse MRI:The application of compressed sensing for Rapid MR imaging[J]. Magn Reson Med, 2007, 58(6): 1182-1195. |
| [11] |
GAO M, XIE H B, YANG G, et al. A synchronized compressed sensing scan-reconstruction scheme in magnetic resonance imaging[J].
Chinese J Magn Reson, 2016, 33(2): 257-268.
高芒, 谢海滨, 杨光, 等. 压缩感知同步扫描重建及其采样方案的研究[J]. 波谱学杂志, 2016, 33(2): 257-268. |
| [12] |
SONG Y, XIE H B, YANG G, et al. Dictionary learning with segmentation for compressed-sensing magnetic resonance imaging[J].
Chinese J Magn Reson, 2016, 33(4): 559-569.
宋阳, 谢海滨, 杨光, 等. 用于压缩感知磁共振成像的分割字典学习算法[J]. 波谱学杂志, 2016, 33(4): 559-569. |
| [13] | WOOLRICH M W, RIPLEY B D, BRADY M, et al. Temporal autocorrelation in univariate linear modeling of FMRI data[J]. NeuroImage, 2001, 14(6): 1370-1386. DOI: 10.1006/nimg.2001.0931. |
| [14] | KUO L W, CHIANG W Y, YEH F C, et al. Diffusion spectrum MRI using body-centered-cubic and half-sphere sampling schemes[J]. J Neurosci Meth, 2013, 212(1): 143-155. DOI: 10.1016/j.jneumeth.2012.09.028. |
| [15] | LE BIHAN D, MANGIN J F, POUPON C, et al. Diffusion tensor imaging:concepts and applications[J]. J Magn Reson Imaging, 2001, 13(4): 534-546. DOI: 10.1002/(ISSN)1522-2586. |
| [16] | LAMBIN P, RIOS-VELAZQUEZ E, LEIJENAAR R, et al. Radiomics:extracting more information from medical images using advanced feature analysis[J]. Eur J Cancer, 2012, 48(4): 441-446. DOI: 10.1016/j.ejca.2011.11.036. |
| [17] | KUMAR V, GU Y, BASU S, et al. Radiomics:the process and the challenges[J]. Magn Reson Imaging, 2012, 30(9): 1234-1248. DOI: 10.1016/j.mri.2012.06.010. |
| [18] | GREENSPAN H, VAN GINNEKEN B, SUMMERS R M. Guest editorial deep learning in medical imaging:Overview and future promise of an exciting new technique[J]. IEEE T Med Imaging, 2016, 35(5): 1153-1159. DOI: 10.1109/TMI.2016.2553401. |
| [19] | HANSEN M S, SORENSEN T S. Gadgetron:an open source framework for medical image reconstruction[J]. Magn Reson Med, 2013, 69(6): 1768-1776. DOI: 10.1002/mrm.24389. |
| [20] | COOK P A, BAI Y, NEDJATI-GILANI S, et al. Proceedings of the 14th Scientific Meeting of the ISMRM[C]. Seattle: Proc Intl Soc, 2006. |
| [21] | DADUCCI A, GERHARD S, GRIFFA A, et al. The connectome mapper:an open-source processing pipeline to map connectomes with MRI[J]. PloS One, 2012, 7(12): e48121. DOI: 10.1371/journal.pone.0048121. |
| [22] | UECKER M, ONG F, TAMIR J I, et al. The ISMRM 23rd Annual Meeting[C]. Toronto: Proc Intl Soc, 2015. |
| [23] | LUSTIG M, DONOHO D L, SANTORS J M, et al. Compressed sensing MRI[J]. IEEE Signal Proc Mag, 2008, 24(2): 72-82. |