改进的数字重建影像快速生成算法

罗博博，贠照强，卢振泰，林仁回，张 娟，冯前进

（南方医科大学 生物医学工程学院，广东 广州510515）

0 引 言

数 字 重 建 影 像 （digitally reconstructed radiographs，DRRs）是对CT图像进行模拟投影方式产生的虚拟X线图像，广泛应用于计算机辅助外科手术、图像引导介入治疗和图像引导放射治疗等领域［1］。特别地，在骨科临床，骨折手术治疗涉及的解剖结构相对复杂，常规手术风险较大，应用导航技术可提高手术的安全性和准确性。DRRs的重建已经成为临床骨科手术导航中的一项十分关键的关键技术。

目前生成DRRs的算法主要包括：光线投影 （ray casting）［2］、摇晃抛雪球 （wobbled splatting）［3］、剪切变形分解（shear－warp factorization）［4］、跨 图 形 学 （the transgraph）［5］、自适应蒙特卡洛体绘制 （adaptive monte carlo volume rendering）［6］、衰 减 场 （attenuation fields）［7］、逐 步 衰 减 场 （progressive attenuation fields）［8］等算法。DRRs的生成是一个耗时相当长的过程，对于N×N×N的图像数据，传统的算法，如投影 （ray casting）算法、抛雪球 （SPLATTING）算法和剪切变形 （SHEAR－warp）等算法的计算复杂度为O（N3）［9］，为满足实时性要求，Daniel B.Russakoff等提出了基于光场的数字重建影像方法，将光场理论引入DRRs的形成过程，大幅度提高了DRRs的生成速度，但该方法内存消耗量大、建立的光场数据有大量冗余、预处理时间长；D.Ruijters等提出了基于图形处理器 GPU （graphic processing unit）的DRRs生成方法［10］，大幅提升了DRRs图像的生成效率。

而在临床骨科实时图像引导手术中，基于灰度的2D／3D医学图像配准给临床医生提供一个到人体内部虚拟的、非侵入式的窗口，使医生能够看到一个肉眼无法直接看到的解剖与手术器械的三维空间相对位置关系，使得手术精确、安全进行。为了快速、准确获取病灶和手术器械的三维空间相对位置关系，常要生成数百张DRRs，并做到实时配准。因此快速DRRs生成方法的研究成为骨科手术导航系统的关键问题。

降采样是对过采样或正常采样的数据进行二次抽取或者降低频率二次采样数据的过程，这样可以减少计算量，节约计算时间，是实时处理中的一种常用的方法。通用并行计算架构 （compute unified device architecture）是一种将GPU作为数据并行计算设备的软硬件体系，可以对数据进行高度的并行计算，在处理图形数据和复杂算法方面拥有比CPU更高的效率。为保证临床骨科手术安全、顺利进行，采用只含骨性结构的DRRs实现术中实时图像引导，本文结合多种加速技术，提出一种改进的算法——基于光线投影算法的DRRs生成原理，采用统一并行计算架构CUDA加速，再引入大步长与只采样骨性结构相结合的降采样，快速生成只含骨性结构的DRRs，并采用B.SrinivasaReddy和 B.N.Chatterji提 出 的 基 于 Fourier－Meilin变换和相位相关相似性测度的医学图像配准方法将通过改进算法得到的DRRs和临床X线图像进行医学图像配准验证。实验结果表明，只采样骨性结构同时不超过3倍 （甚至5倍）原步长采样时，改进的算法能快速生成只含骨性结构的DRRs，也能使得获取的DRRs和临床X线精确配准。但采样量需根据Volume数据大小而控制，Volume数据量越大，采样范围更大，改进算法的提速优势更为明显；如果超出理想采样范围，该改进算法获取的DRRs与临床X线图像配准效果不理想，不宜应用于临床。

1 改进的算法原理及方法

DRRs的生成是体绘制算法的一种典型应用，它是模拟X线穿透CT体元，经过衰减和吸收后投影到成像平面的过程，本文使用透视投影，如图1所示。

DRRs生成过程中，如果入射光源的强度为I0，任何一光线穿过 “人体”后到达投影平面的强度为I，连续型积分表达式［11］为

式中：μ（x）——对应组织的衰减系数，0——光线经过体数据的穿入点，D——光线经过体数据的穿出点。因为人体组织是非均匀的，不同的组织对应不同的衰减系数，离散型表达式［11］为

下标i为沿着射线路径穿过的第i种组织的索引，xi是穿过第i种组织的有效长度，μi是对应第i种组织的衰减系数，衰减系数μi与CT值的转换公式参考如下［12］

图1 光线投影算法

式中：μw——水的衰减系数，HUi——第i种组织的CT值，范围为 ［－1000，3096］，单位为Hounsfield Unit。

光线投影算法生成DRRs的具体过程描述如下：光线穿过包含CT体数据的长方体，沿着光线的行进方向以一定大小的步长采样CT值，把得到的CT值转换成相应的衰减系数μi，累加此方向上的衰减系数，由此得到所有射线到达投影板上的强度I，最后形成DRRs。

上述DRRs生成过程中，每条穿透CT体元的模拟X线互不交叉，且得到投影板上像素值的计算形式相同，此特点正好与统一计算设备架构CUDA的并行计算体系相吻合，因此，借助CUDA技术进行高度并行计算，实现多条光线的同时绘制，大幅度加快了DRRs的形成。同时，穿透CT体元的任何模拟X线路径都可表示为u→＋Δt v→，u→是模拟X线的起点，v→是模拟X线的方向，采样开始于模拟X线进入CT体元的起点tin，结束于光线穿出CT体元的终点tout。而降采样是实时处理中一种常用方法。在体绘制中，它主要包括增大采样步长Δt、增大起始点tin或减小终点tout、根据特征抽取特定信息和减少穿透CT体元的模拟X线的数量等方式。DRRs生成时，如果采用增大采样步长Δt和只采样骨头组织相结合的降采样，将大大减少DRRs过程的计算量，快速生成只含骨性结构的DRRs。通过对颈部、胸部、腰部、腿部等多幅骨科临床CT数据的骨窗和灰度直方图分析，得知人体软组织的灰度值小于860，骨组织的灰度值大于860。如果将整幅图像的灰度值归一化为0－1，采用灰度阈值采样的方法，即采样过程中只采样灰度值大于0.2670的像素，同时增大采样步长 （3倍甚至5倍原步长），结合CUDA硬件加速，快速获取只含骨性结构的DRRs。本文针对临床骨科手术导航系统，该临床应用需获取只含骨性结构的DRRs来实现术中实时图像引导。多次实验证明，采用该改进的算法——基于光线投影算法的DRRs生成原理，采用硬件CUDA加速，同时引入大步长与只采样骨性结构相结合的降采样生成DRRs，可以大大减少重建DRRs的计算量，大幅度节约了生成DRRs的时间［13］，满足了临床骨科手术导航系统的实时性要求。

2 基于Fourier－Meilin变换和相位相关的图像配准

B.SrinivasaReddy和 B.N.Chatterji提出的 基于 Fourier－Meilin变换和相位相关相似性测度的医学图像配准方法［14］具有运算量小、抗干扰能力强、配准精确度高和鲁棒性强的等优点，为验证DRRs图像质量，本文采用此方法对获取的DRRs和临床的X线图像进行配准。它将时域上的医学图像灰度经过傅里叶变换到频域上再求频谱相位相关性——互功率谱特性，然后对目标互功率谱进行反傅里叶变换找到相应的峰值位置，从而确定配准参数。具体过程如下：

假设只存在平移变换时，f1（x，y）和f2（x，y）为两个图像信号，它们满足关系 （即f2（x，y）是由f1（x，y）经过平移得到）

根据傅里叶变换的性质可得

F1（u，v）和F2（u，v）分别为f1（x，y）和f2（x，y）经过傅里叶变换的结果，它们的互功率谱为

式中：F＊1（u，v）——F1（u，v）的复共轭，的傅里叶反变换为一个二维脉冲函数δ（x－x0，y－y0），相位相关法就是求取式（6）的反傅里叶变换，然后找到峰值位置最终确定配准参数。

当平移、旋转、缩放参数同时存在时，设第一幅图像为f1（x，y），第二幅图像f2（x，y）是由第一幅图像f1（x，y）平移、旋转、缩放得到，可表示为

式中：σ——缩放尺度，θ0——旋转角度，（x0，y0）——平移参量；配准就是在参考图像f1（x，y）中心处取一个图像直角坐标尺寸为N×N的小区域I1（x，y），配准图像f2（x，y）任何可能的位置截取一个与I1（x，y）大小相同的小区域I2（x，y），其次在对数极坐标表示下，对I1（x，y）和I2（x，y）分别进行傅里叶变换后求两者的互功率谱，最后对互功率谱进行傅里叶反变换，如果产生一个二维脉冲信号，则找出峰值位置估计出配准参数；否则配准失败。

3 结果与分析

3.1 数字重建影像与采样时间结果

本文按上述方法，在一台CPU主频为3.20GHZ、内存为2G、显卡为Nvidia Quadro 600、显存为1G的PC机上处理了A、B、C三组16位CT数据，分辨率分别为250×512×512、831×512×512、1440×512×512，以下分别为A、B、C数据的DRRs效果图和采样时间的结果（λ＝0.01为原步长，B为只采样骨性结构）。

图2（a）为只借助CUDA加速进行原步长采样获取的A数据DRR结果，图2（b）至图2（f）为基于改进的算法获取的A数据DRRs结果。

图2 不同采样方式获取的A数据DRRs

表1～表3为DRRs生成过程中不同采样方式分别获取一张A、B、C数据DRR的时间对照表。

表1 不同采样方式获取一张A数据DRR的时间对照

图3（a）为只借助CUDA加速进行原步长采样获取的B数据DRR结果，图3（b）至图3（f）为基于改进的算法获取的B数据DRRs结果。

图3 不同采样方式获取的B数据DRRs

表2为B数据DRRs生成过程中不同采样方式获取一张DRR的时间对照表。

表2 不同采样方式获取一张B数据DRR的时间对照

图4（a）为只借助CUDA加速进行原步长采样获取的C数据DRR结果，图4（b）至图4（j）为基于改进的算法获取的C数据DRRs结果。

表3为C数据DRRs生成过程中不同采样方式获取一张DRR的时间对照表。

图4 不同采样方式获取的C数据DRRs

表1结果显示，3λ＋B采样方式获取A数据DRRs相比λ采样方式获取A数据DRRs的时间减少了53.60ms；表2结果显示，3λ＋B采样方式获取B数据DRRs相比λ采样方式获取B数据DRRs的时间减少了67.90ms；表3结果显示，5λ＋B采样方式获取C数据DRRs相比λ采样方式获取C数据DRRs的时间减少了100.60ms。由此得知，通过增大采样步长与只采样骨性结构相结合的降采样可以更大幅度提高DRRs的生成速度，且数据量越大，提速优势更明显。

表3 不同采样方式获取一张C数据DRR的时间对照

3.2 不同采样方式得到的DRRs与临床X线图像的配准

改进的算法可以大幅度提高DRRs的生成速度，同时保证了DRRs与临床X线图像的精确配准。为验证获取的DRRs图 像 质 量，本 文 采 用 B.SrinivasaReddy和B.N.Chatterji提出的基于Fourier－Meilin变换和相位相关相似性测度的医学图像配准方法对DRRs进行配准验证，此处将A数据DRRs和A数据对应的临床X线、B数据DRRs和B数据对应的临床X线、C数据DRRs和C数据对应的临床X线配准 （用原步长采样方式得到的模拟X线图像代替临床X线图像，并设为参考图像，降采样得到的图像为浮动图像）。表4（a）显示了A数据的DRRs和对应临床X线图像配准的结果，表4（b）显示了B数据DRRs和对应临床X线图像配准的结果，表4（c）显示了C数据的DRRs和对应临床X线图像配准的结果；图5（a）、5图（b）、图5（c）分别显示了采样方式与配准旋转角度平均偏差、X方向平均偏差、Y方向平均偏差的关系。

表4 A、B、C数据DRRs分别与对应临床X线图像的配准结果

（c）C数据的DRRs和临床X线图像配准结果

图5 不同采样方式与配准偏差关系

表4（a）和图5（a）～图5（c）结果显示，λ＋B、2λ＋B、3λ＋B方式获取的A数据DRRs分别与A数据对应的临床X线图像配准的误差和λ方式获取的A数据DRRs与A数据对应的临床X线图像配准的误差很相近（如果计算出的变换参数和真实值相差不超过1个像素或1度，即认为配准成功），以上误差均在临床可接受的范围内，分析可知增大采样步长和只采样骨头相结合的降采样不影响获取的A数据DRRs和A数据对应的临床X线图像的精确配准，但实验结果证明，如果继续减少采样量 （小于3λ＋B的采样量），获取的A数据DRRs和A数据对应的临床X线图像的配准将产生很大的偏差，很难保证实时图像引导手术治疗的准确性，不宜采用；按上所述，表4（b）和图5（a）～图5（c）结果表明，增大采样步长和只采样骨头相结合的降采样不影响获取的B数据DRRs和B数据对应的临床X线图像的精确配准，但实验结果证明，如果继续减少采样量 （小于3λ＋B的采样量），获取的B数据DRRs和B数据对应的临床X线图像的配准将产生很大的偏差，很难保证实时图像引导手术治疗的准确性，不宜采用；表4（c）和图5（a）～图5（c）结果表明，增大采样步长和只采样骨头相结合的降采样不影响获取的C数据DRRs和C数据对应的临床X线图像的精确配准，但实验结果证明，如果继续减少采样量 （小于5λ＋B的采样量），获取的C数据DRRs和C数据对应的临床X线图像的配准将产生很大的偏差，很难保证实时图像引导手术治疗的准确性，不宜采用。

由上可知，当采样量控制在理想范围内 （一般不低于3或5λ＋B的采样量）时，改进的算法可以更大幅度加快DRRs的形成，也能使得获取的DRRs和临床X线精确配准；如果超出理想采样范围，获取的DRRs不能保证和临床X线的精确配准。于此，采样量的选择需根据Volume数据量的大小而控制，需做到既能加快DRRs的形成，又不影响DRRs和临床X线的精确配准；并且Volume数据量越大，采样范围更大，改进的算法的优势会更为明显。

4 结束语

本文提出了一种改进的数字重建影像生成算法——基于光线投影算法的DRRs生成原理，采用统一并行计算架构CUDA加速，再引入大步长和只采样骨性结构相结合的降采样，快速生成只含骨性结构的DRRs，并将改进算法获取的DRRs和临床X线图像进行了医学图像配准验证。实验结果表明，当采样量控制在理想范围内时，本文采用改进的算法获取DRRs不但在提速上优势很明显而且能保证获取的DRRs与临床X线图像的精确配准，从而实现实时图像引导临床骨科等手术，满足了临床骨科手术导航系统等临床应用的实时性和精确性的要求。

［1］Markelj P，TomaeviD，Likar B，et al.A review of 2D／3D registration methods for image－guided interventions ［J］.Medical Image Analysis，2012，16 （3）：642－661.

［2］Deutschmann H，Steininger P，Nairz O，et al. “Augmented reality”in conventional simulation by projection of 3－D structures into 2－D images ［J］.Strahlentherapie und Onkologie，2008，184 （2）：93－99.

［3］Jakob Spoerk，Helmar Bergmann，Felix Wanschitz，et al.Fast DRR splat rendering using common consumer graphics hardware ［J］.Med Phys，2007，34 （11）：4302－4308.

［4］Hadjira Bentoumi，Pascal Gautron，Kadi Bouatouch.GPU－based volume rendering for medical imagery ［J］.World Academy of Science，Engineering and Technology，2010，4 （1）：419－125.

［5］Gu Y，Wang C.Transgraph：Hierarchical exploration of transition relationships in time－varying volumetric data ［J］.IEEE Trans Vis Comput Graph，2011，17 （12）：2015－2024.

［6］Gao Y，Haapasalo M，Shen Y，et al.Development of virtual simulation platform for investigation of the radiographic features of periapical bone lesion ［J］.Journal of Endodontics，2010，36 （8）：1404－1409.

［7］Baumhauer M，Feuerstein M，Meinzer H P，et al.Navigation in endoscopic soft tissue surgery：Perspectives and limitations［J］.Journal of Endourology，2008，22 （4）：751－766.

［8］Mori S，Kobayashi M，Kumagai M，et al.Development of a GPU－based multithreaded software application to calculate digitally reconstructed radiographs for radiotherapy ［J］.Radiological Physics and Technology，2009，2 （1）：40－45.

［9］Li X，Zhou L，Zhen X，et al.The generation of digitally reconstructed radiographs with six parameters ［C］／／4th International Conference on Bioinformatics and Biomedical Engineering.IEEE，2010：1－4.

［10］Fluck O，Vetter C，Wein W，et al.A survey of medical image registration on graphics hardware ［J］.Computer Methods and Programs in Biomedicine，2011，104 （3）：45－57.

［11］Ruijters D，ter Haar Romeny B M，Suetens P.GPU－accelerated digitally reconstructed radiographs ［C］／／Proceedings of the Sixth IASTED International Conference on Biomedical Engineering，2008：431－435.

［12］Jaffray D，Kupelian P，Djemil T，et al.Review of imageguided radiation therapy ［J］.Expert Review of Anticancer Therapy，2007，7 （1）：89－103.

［13］Osama M Dorgham，Stephen D Laycock，Mark H Fisher.GPU accelerated generation of digitally reconstructed radiographs for 2D／3Dimage registration ［J］.IEEE transactions on Biomedical Engineering，2012，59 （9）：2594－2603.

［14］Lin Y H，Chen C H.Template matching using the parametric template vector with translation，rotation and scale invariance［J］.Pattern Recognition，2008，41 （7）：2413－2421.