本发明属于图像处理领域，尤其涉及一种基于cuda的医学图像的三维重建方法。
背景技术：
三维医学数据的体绘制技术在可视化中具有十分重要的意义，由于医学图像数据大多数属于规则的结构化体数据(structuredregularvolumedata)，分布在正方体或长方体组成的三维网格点上，因此本文研究的内容都是围绕规则结构化体数据展开的。对于规则数据场的体绘制本文采用体绘制算法中比较典型的光线投射算法(raycasting)。由于光线投射算法只与投射光线的数量密切相关，因此特别适合于数据场规模大、数据集比较规则的三维数据场的体绘制。光线投射法的原理相对简单，过程易于实现，并且可以很容易地实现透视投影，绘制的图像质量也相对较高，因此适合作为医学图像三维重建的算法。
虽然光线投射算法在规则数据场中有很多的优点，但是也存在着一些不尽完善的地方，如：因为运算量比较大导致绘制速度比较慢，难以满足实时性要求。由于光线投射算法的实现与光线数量有关，处理速度一直是一个难点，单纯在算法上进行的优化很难保证三维重建的实时性。近年来，随着图形技术的发展，很多研究人员都对医学数据的体绘制进行了探索。随着gpu的发展，出现了一些基于gpu的加速算法。kruger和westermann等人在传统的光线投射算法基础上基于gpu并行计算特点对算法进行重构，它使用顶点着色程序计算投射光线参数、通过利用深度测试与阻塞询问模拟循环实现光线积分。christof和kolb等人提出了不透明度剥离的体绘制算法。尽管上述基于gpu的光线投射算法加快了体绘制速度，大大超过了只在cpu上运行的算法，但是为了达到理想的绘制速度而放弃了编程的灵活性。
2006年11月，nvidia公司推出了cuda(通用并行计算)体系构架，该构架是一种新的并行编程模型和指令集的通用计算构架，它是基于高速图像处理单元gpu的，在gpu上高速并发执行，极大的提高了程序算法的运行速度，该体系构架为gpu编程提供了一种全新的软硬件构架。利用nvidiagpu的并行计算引擎，线程、线程块和网格可以进行比cpu更高效的复杂计算任务。cuda集成了cpu和gpu各自的优势，内核部分在gpu上执行，其余部分则继续在cpu上完成。l.marsalek最早在cuda构架中对体绘制进行了加速尝试，并且证明了cuda能有效提高体绘制的速度。以上的研究表明了对医学图像的三维重建用cuda对光线投射算法进行优化加速，可以极大提高运算的效率和算法的执行速度。
在cuda编程环境中，主要包括cpu和gpu两个部分。cpu作为主机，即host端，gpu作为设备，即device端。host端与device端有专用的通道进行数据通信，host端负责对逻辑性事务进行处理，以及对串行化运算的控制；device端负责执行大规模的并行化处理任务。将运行在图像处理单元gpu上的cuda并行计算函数称为核函数，即kernel函数。
技术实现要素：
本发明提供一种基于cuda的医学图像的三维重建方法，在保证了较满意的重建效果的前提下，提高了算法的效率，使得实时绘制更加流畅。
为实现上述目的，本发明采用如下的技术方案：
一种基于cuda的医学图像的三维重建方法，包括以下步骤：
步骤1、体数据映射成三维纹理
将unsignedchar类型的体数据映射成可以被gpu读入的三维纹理，在体数据到三维纹理的映射过程中，可现实多种传递函数，进行灰度绘制时，直接使用体数据值作为三维纹理的灰度值；而进行光照绘制时，将每一个体数据点的梯度值作为三维纹理的rgb颜色值，使用中心差分公式计算梯度值，而将原始的体数据值作为透明度值；
步骤2、确定光线的起始位置和方向
在世界坐标系中计算从视点发出的光线，确定光线的起始位置和起始方向，将体数据视为放置在视空间中的一个立方体，记boxmin＝{xmin,ymin,zmin}为该立方体在视空间下所有定点的x,y,z坐标的最小值，boxmax＝{xmax,ymax,zmax}为最大值，把视点设为光线的起始点，在opengl中，视点固定在原点，在计算光线方向时，把光线在体数据立方体上的进入坐标值作为光线方向的值；
步骤3、图像合成
采用从前向后的图像合成算法，在内核函数中循环对一条光线上的体数据进行遍历、采样、赋色并用改进的blinn-phong光照模型进行光照计算，得出颜色值和不透明度；首先计算出光线在体数据立方体上的离开点；并遍历循环进入点到离开点的光线步进，计算每次步进后所在体数据点的颜色值；最后混合相同光线上所有采样点的值作为该像素点的颜色值输出。
步骤4、绘制计算结果
将最后得到的颜色值与不透明度写入像素缓冲对象(pbo)中，这样就从cuda的内核传入到cpu端，然后用opengl绘制到屏幕即可。
作为优选，步骤4具体为：在写入像素缓冲对象(pbo)过程中，先建立用来当做输出结果的pbo，然后cuda部分注册pbo这个缓冲对象；在kernel程序中通过就可以取得体位置；使用完成后，先把缓冲对象删除，然后再取消缓冲对象的注册；最后在opengl中通过相关函数来显示到屏幕，显示部分的内容分为三方面：处理摄像机的位置矩阵、调用绘制函数进行绘制和绘制pbo。
作为优选，步骤3具体为：
1)计算离开点
判断一条射线与包围盒体是否相交，并计算出相交点。射线方向为(xdydzd)，射线出发点为(xoyozo)，射线方程为
2)计算采样点的颜色值
将同一条光线上的采样点的颜色进行混合，求得最终的颜色值为了实现最前面的效果，使用从前往后的混合，其混合算法使用本文优化后的算法，算法方程如公式(4)所示：
光线终止条件：当光线穿越体数据的距离大于光线与体数据的最大交点的最小值时，或者不透明度接近阈值时，光线结束，得到最终的颜色值与不透明度；
实现灰度效果绘制时，直接把体数据值作为灰度颜色和不透明度，可设定不同的阈值来获得不同的重建图像；实现光照效果效果时，首先计算出梯度，并以此作为体数据的法向量，然后在根据改进的blinn-phong光照模型计算出光照模型处理后的颜色值，改进的光照模型如公式(9)所示：
附图说明
图1为从前到后的图像合成算法原理示意图。
图2为光线投射算法的核函数处理流程。
图3为视点光线与包围盒相交示意图。
图4为分辨率为128*128*53的头部渲染效果图。
具体实施方式
本发明公开一种基于cuda的医学图像的三维重建方法，在保证了较满意的重建效果的前提下，提高了算法的效率，使得实时绘制更加流畅。
1、基于光线投射算法的改进。
虽然光线投射算法是有效的三维重建方法，但渲染速度慢。本发明针对体绘制的具体处理过程：重采样、图像分类及渲染、图像合成等进行研究，对传统的光线投射算法进行了以下改进。
1)对图像合成算法进行优化。
图像合成就是将所有计算出来的采样点的光线按照某种特定的规则累积到一起，最终求得像素点的颜色值与不透明度。这里采用从前到后的图像合成方法，从前向后合成算法的好处是，从i到n时
由于
由于公式(3)的计算涉及五次的加减乘复合运算，在图像合成过程中计算量太大，因此考虑对公式(3)进行优化。已知不透明度α的分布范围是0到1，在从前向后的合成过程中，α越来越趋近为1。把光线初始点的不透明度设为零，即
2)对blinn-phong光照模型的实现进行改进。
在blinn-phong光照模型的实现过程中，通过设置调整参数来增强基于cuda实现的光线投射算法的真实感效果。与传统的blinn-phong光照模型一样，改进的blinn-phong光照模型在实现过程中同样也结合了lambert的漫反射光和标准的高光；光照模型中体素的光照强度同样由环境光、漫反射光以及镜面反射光组成。
传统的blinn-phong光照模型如公式(5)所示：
其中：ia、id和imr分别代表环境光、漫反射光和镜面反射光；n是光源个数。
环境反射光亮度表示为：ia＝kail，ka为物体表面的环境光反射系数(0≤ka≤1)，il为光源强度。
漫反射光亮度表示为：id＝kdil(n·l)，kd为物体表面的漫反射系数(0≤kd≤1)，il为光源强度，n是物体表面的单位法向量，l是物体表面一点指向点光源的单位向量。
镜面反射光亮度表示为：imr＝kmril(n·h)σ，kmr为物体表面的镜面反射系数(0≤kd≤1)，il为光源强度，σ是高光指数(代表物体表面的光泽程度，较光滑的物体表面取较大的值，粗糙的物体表面取较小的值)，n是物体表面的单位法向量，h是光入射方向l和视点方向v的平分向量，也称半角向量。h的计算公式如公式(6)所示：
在blinn-phong光照模型中通过调节各个系数可以调整光照模型中各部分的比例，通常我们设置ka比较小，这样可以降低环境光在光照模型中的比例。blinn-phong光照模型整体上光照实现已经有较好的效果，但是对于要实现明暗轮廓清晰的绘制，实现效果就不太理想。
在需要突出渲染一些细节时，一般采用可变步长。使用比较普遍的可变步长的计算公式如公式(7)所示：
其中：d是不变的标准采样步长，α是调节系数，
其中：δ是梯度计算时的坐标偏移量。
为了能给局部增强光照强度，引入一个新的参数μ，
当δ设置大一些时，
在实际绘制过程中，为了提计算的速度，本文使用单一光源。采用改进的blinn-phong光照模型对光线穿越体数据时的每个体素进行明暗渲染，局部绘制效果更佳，明暗轮廓更清晰，更加具有真实感。
2、基于cuda的绘制。
本发明从编程模型出发，利用gpu的simd，在cuda并行构架下对改进后的光线投射算法进行实现，实现了医学图像的三维重建，绘制效果和实时性都较好。
光线投射算法中最耗时的计算是每条光线上的重采样、颜色和不透明度的累积计算。在cuda构架中，所有的线程都被划分成线程块，cuda是以线程块为执行单位的，所有的线程块组成网格。由于每个线程块的索引号是唯一的，这样就可以使用线程块的索引号对存储在gpu上的数据进行读写。
由于在光线投射算法中，每条光线都是独立计算的，计算结果不影响其它光线的计算值，这样就可以把每条光线作为一个线程来处理。在光线采样过程中，每条光线都是并行计算的。该过程在cuda上的内核处理流程图如图2所示。
光线的生成以及包围盒判断。用线程在网格中的索引(x,y)来表示屏幕上像素的坐标(u,v),这样在视空间下，就建立了线程索引和视平面像素坐标的对应关系。在初始化阶段，把视点作为起始点，把一个视平面上的像素坐标和视点坐标的差作为方向，生成一条穿越体数据的光线。光线方程如公式(10)所示，o为光线起点，d为方向。如果与包围盒最大坐标值相交，则返回光线的入点和出点。
p＝o+td(10)
从入点到出点进行光线步进，计算每次步进后的颜色值。计算当前采样点的梯度值，并作为法向量，按照改进的blinn-phong光照模型计算出光照后的颜色，然后混合该条光线上所有采样点的值作为该点像素的颜色值。
将最后得到的颜色值与不透明度写入像素缓冲对象(pbo)中，这样就从cuda的内核传入到cpu端，然后用opengl绘制到屏幕即可。
实施例1
光线投射算法在cuda上实现的关键是能够利用gpu完成对体数据的遍历、采样和计算，并将得到的采样值进行混合，以产生最终的重建结果。
用cuda实现光线投射算法，需要声明一个内核函数kernel来进行体数据的遍历、光线方向上的等距离采样以及颜色值的累加计算。实现过程中，每条光线相当于一个线程，多个线程可以并行执行，这样就可以并行计算多条光线。利用loadrawfile函数在cpu主机端读入体数据，存放在系统内存中；然后通过系统总线把体数据从系统内存上拷贝到gpu设备端的全局存储器中，根据体数据的大小分配线程块和线程的数目，在内核函数中，每条光线分配一个线程，进行并行计算。本文中显卡的sm最多有3072个寄存器，每个线程需要16个寄存器，则每个线程块最多允许256个线程。cuda在执行程序的时候，以wrap为基本单位，目前的cuda装置中，一个wrap里面有32个thread，分成两组16thread的half-wrap。因此得出线程块维度为16*16，也就是每个线程块中包含256个并行执行的线程。在本文的绘制用例中，体数据规格为128*128*53。
实现步骤如下：
步骤1、体数据映射成三维纹理。将unsignedchar类型的体数据映射成可以被gpu读入的三维纹理(3dtexture)。在体数据到三维纹理的映射过程中，可现实多种传递函数，进行灰度绘制时，直接使用体数据值作为三维纹理的灰度值；而进行光照绘制时，将每一个体数据点的梯度值作为三维纹理的rgb颜色值，使用中心差分公式计算梯度值，而将原始的体数据值作为透明度值。
步骤2、确定光线的起始位置和方向。在世界坐标系中计算从视点发出的光线，确定光线的起始位置和起始方向。为了进行三维重建，将体数据视为放置在视空间中的一个立方体，记boxmin＝{xmin,ymin,zmin}为该立方体在视空间下所有定点的x,y,z坐标的最小值，boxmax＝{xmax,ymax,zmax}为最大值。把视点设为光线的起始点，在opengl中，视点固定在原点，因此在计算光线方向时，可以把光线在体数据立方体上的进入坐标值作为光线方向的值。
步骤3、图像合成。采用从前向后的图像合成算法，在内核函数中循环对一条光线上的体数据进行遍历，采样，赋色并用改进的blinn-phong光照模型进行光照计算，得出颜色值和不透明度。首先计算出光线在体数据立方体上的离开点；并遍历循环进入点到离开点的光线步进，计算每次步进后所在体数据点的颜色值；最后混合相同光线上所有采样点的值作为该像素点的颜色值输出。
1)计算离开点
判断一条射线与包围盒体是否相交，并计算出相交点。射线方向为(xdydzd)，射线出发点为(xoyozo)，射线方程为
计算离开点的代码片段如下：
计算光线与盒体六个平面的交点
float3invr＝make_float3(1.0f)/r.d；
float3tbot＝invr*(boxmin-r.o)；//计算进入点
float3ttop＝invr*(boxmax-r.o)；//计算离开点
计算每个轴平面的最小和最大交点坐标
float3tmin＝fminf(ttop,tbot)；//进入点坐标
float3tmax＝fmaxf(ttop,tbot)；//离开点坐标
floatlargest_tmin＝fmaxf(fmaxf(tmin.x,tmin.y),fmaxf(tmin.x,tmin.z))；
floatsmallest_tmax＝fminf(fminf(tmax.x,tmax.y),fminf(tmax.x,tmax.z))；
tnear＝largest_tmin；//获得最大的tnear
tfar＝smallest_tmax；//获得最小的tfar
其中，r.o是光线原点，r.d是光线方向，设光线参数方程为p＝r.o+t×r.d，则tnear是进入点的t值，tfar是离开点的t值。
2)计算采样点的颜色值
将同一条光线上的采样点的颜色进行混合，求得最终的颜色值为了实现最前面的效果。本文使用从前往后的混合，其混合算法使用本文优化后的算法。算法方程如公式(4)所示：
光线终止条件：当光线穿越体数据的距离大于光线与体数据的最大交点的最小值时，或者不透明度接近阈值时，光线结束。得到最终的颜色值与不透明度。
每条光线采样、合成的部分代码如下：
实现灰度效果绘制时，直接把体数据值作为灰度颜色和不透明度，可设定不同的阈值来获得不同的重建图像。实现光照效果效果时，首先计算出梯度，并以此作为体数据的法向量，然后在根据改进的blinn-phong光照模型计算出光照模型处理后的颜色值。改进的光照模型如公式(9)所示：
有时候要用彩色来呈现或凸出某些部位的话，可以用一个transferfunction来做色彩的对应。
使用一个float4的阵列transferfunc，这里面有九个控制点，分别代表灰阶值不同所要呈现的阶段色阶，介于中间的值，使用线性插值来计算获得。
step4.绘制计算结果。将最后得到的颜色值与不透明度写入像素缓冲对象(pbo)中，这样就从cuda的内核传入到cpu端，然后用opengl绘制到屏幕即可。
在写入像素缓冲对象(pbo)过程中，先建立用来当做输出结果的pbo，然后cuda部分注册pbo这个缓冲对象；这样在kernel程序中通过就可以取得体位置。使用完成后，先把缓冲对象删除，然后再取消缓冲对象的注册。最后在opengl中通过相关函数来显示到屏幕，显示部分的内容分为三方面：处理摄像机的位置矩阵、调用绘制函数进行绘制和绘制pbo。
位置矩阵的作用是可以使用鼠标拖动控制摄像机的旋转和平移，把资料存在viewrotation和
viewtranslation中；位置矩阵的计算通过opengl的modelviewmatrix来计算。
绘制过程包括以下内容：把矩阵复制到常量设备内存、把pbo映射到设备内存d_output、调用kernel以及解除pbo和d_output的映射关系。最后在opengl中，先指定pbo缓冲对象，然后绘制出来。
本发明在intel/xeone52.5ghzcpu和win764为操作系统下，使用c++语言和开放性图形接口opengl在vs2010和cuda5.5中实现光线投射算法的体绘制，在512*512分辨率下达到了50+fps的效果。图4是一个head体数据的绘制效果图。