Volume Visualization with Raycasting(GPU+OPENGL)

<script src="http://widgets.amung.us/classic.js" type="text/javascript"></script> <script type="text/javascript"> </script>

Volume Visualization with Raycasting

 

前面已经介绍了scalar visualization和vector visualization，现在介绍一下volume visualization，中文里，我们称为体渲染或体可视化。虽然体可视化也是一种scalar可视化，但是它和一般scalar可视化有很大的区别。把它称为体可视化，是因为我们处理的数据对象是一个容积，用的最多的就是医学上，比如CT，MRI，PET，SPECT等等。

比如有一个对头部扫描的数据，拥有256×256×256个体素(voxel)，可以把这个数据想象成一个长宽高为256的一个正方体，我们称为bounding box。在这个bounding box里包含了一个人头的CT扫描来的数据。我们可以想象这个正方体由256×256×256个体素组成，每一个体素就表示一个数据，如下图所示。

 

现在的问题是，我们如何在一个2D平面上来显示这个正方体内所包含的数据。它的基本思想很简单，我们可以想象从相机位置出发并通过一个平面的射线。这个平面，就是我们所要显示画面的平面。这条射线并不是随意穿过这个平面，而是依次穿过这个平面画面的每一个像素(pixel)，如下图所示。

 

 

在上图中可以看到，从眼睛位置发射出射线r，然后通过一平面的所有像素和体数据相交于f和l。如果射线和平面的交点为p(也就是像素)，那么可以得到射线r的参数方程:。有了射线上点，再通过一函数，这个函数称为ray function，它的作用主要是把我们得到的射线点的位置转换成2D图像像素对应的RGB值,。这个函数说明了把体数据转换为2D图像上的一像素的RGB值的过程。函数F的定义域就是射线和体数据的焦点。这中体渲染的方法称为Raycasting。

   

下面来看看如何具体实现这个过程。首先，要得到这个射线的方程，我们需要射线的方向。这里有两种获得方向的方法，一是用2D画面上像素的位置减去相机的位置；一是用bounding box的后表面像素位置减去前面表像素位置。如上图中绿色表面上的像素减去蓝色表面上的像素。得到射线上取样点的scalar值后，利用这个scalar值进行混合(blending)，也就是对每一个样点的值进行混合，就可以得到在2D画面对应的像素的值。如果最终像素用C表示，那么

其中n表示射线上总的样点数，i表示第i个样点，α表示透明度。

 

下面是主要代码

首先要从文件中读取出体数据，然后生成一个3D贴图。

void createTex3D(const char* path) { data=new unsigned char[voxcelNum]; ifstream file(path, ios::in | ios::binary); file.read((char*)data, voxcelNum*sizeof(unsigned char)); file.close(); glPixelStorei(GL_UNPACK_ALIGNMENT, 1); glGenTextures(1, &VolumeTex3D); glBindTexture(GL_TEXTURE_3D, VolumeTex3D); glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE); glTexImage3D(GL_TEXTURE_3D, 0, GL_LUMINANCE, width,height,depth,0, GL_LUMINANCE, GL_UNSIGNED_BYTE, data); useSamplerParameter(CgEffect, "VolumeTex3D", VolumeTex3D); }

 

然后生成bounding box的前后表面，由于使用GPU来计算，所以要把它所谓参数传入cgfx文件中。

 

//生成bounding box的前后表面，作为GPU的参数 glGenTextures(1, &BackTex2D); glBindTexture(GL_TEXTURE_2D, BackTex2D); glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE); glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, 512, 512, 0, GL_RGBA, GL_UNSIGNED_BYTE, NULL); useSamplerParameter(CgEffect, "BackTex2D", BackTex2D); glGenTextures(1, &FrontTex2D); glBindTexture(GL_TEXTURE_2D, FrontTex2D); glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE); glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, 512, 512, 0, GL_RGBA, GL_UNSIGNED_BYTE, NULL); useSamplerParameter(CgEffect, "FrontTex2D", FrontTex2D);

 

最后使用GPU来完成我们raycasting的工作

float4 RayCastingFS( vfConn input ) : COLOR0 { //把bounding box的坐标转换为0-1之间 float2 texCoord = ((input.texCoord.xy / input.texCoord.w) + 1.0) * 0.5; float3 front = tex2D( FrontTex2D, texCoord ).xyz; //前面表 float3 back = tex2D( BackTex2D, texCoord ).xyz; //后表面 float3 dir = normalize( back - front); //获得方向 float4 pos = float4( front, 0 ); float4 src = float4(0.0, 0.0, 0.0, 0.0); float4 dst = float4(0.0, 0.0, 0.0, 0.0); float val = 0.0; float3 step = dir * 0.02; //raycasting 的主要部分，对每个体素进行混合并计算出最后的像素的值 //这里alpha的值设置为0.5，取样步长为0.02 for(int i = 0; i < 300; ++i) { val = tex3D(VolumeTex3D, pos); src = (float4)val; src.a *= 0.5f; src.xyz *= src.a; dst = (1.0f - dst.a) * src + dst; if( 0.95f < dst.a) break; pos.xyz += step; if( 1.0f < pos.x || 1.0f < pos.y || 1.0f < pos.z ) break; } return dst; }

 

---

 

 

下面是可视化的效果

 

 

 

 

需要原始数据的朋友，可以到这里免费下载

http://www.gris.uni-tuebingen.de/edu/areas/scivis/volren/datasets/datasets.html

*原创文章，转载请注明出处*