GPU光线投射算法
Ray Tracing(光线跟踪)算法是标量场的一种直接可视化方法。它从将屏幕作为光线接收的一个平面,根据光路可逆原理,从屏幕各个像素上发出若干光线并模拟在体数据中的投射、反射、折射,得到光线的积分,作为屏幕像素的值。由于其模拟光线的物理原理,可以得到很高的图像质量,广泛应用于电影、高质量动画的制作。虽然其原理简单,但由于需要计算大量的光线,其运算量很大,要在现有的硬件条件下达到中等解析度下实时的速度都较为困难。
Ray Casting(光线投射)算法作为光线追踪算法的一个简化,仅考虑光线在体中的直线运动,不考虑折射和反射。其运算较为简单,成像质量能够达到数控仿真的要求。
基本原理
|
|
|
图 7-5 光线投射基本原理 |
图 7-5是光线投射算法的基本原理。所有的光线从光线出发点出发,经过成像平面,再经过体数据。所有的光线都在视锥之内。考虑将光路逆反,从体数据穿过,穿过时根据传递函数(Transfer Function)得到体数据在离散点上的颜色和透明度,并对其进行积分得到最终的颜色值。所有的光线经过成像平面后,最后得到体数据的光线投射算法的图像。
在光线追踪算法中,有两种重要的辐射物理模型:发散和吸收。发散模型从一个点向外发射一束光线或者吸收能量后再进行发射一束光线。发射模型增强总光线强度。吸收模型从一个点接收光线,然后经过衰减后再发射一束光线。吸收模型降低总光线强度。这两种物理模型均可以包含散射的物理模型。
为了简化处理,光线投射模型仅考虑无散射和无折射的物理模型。
光线积分
|
|
|
|
(a) |
(b) |
|
图 7- 6光线积分 |
|
设初始光强为
,如果不考虑光线的衰减,则可以得到眼睛看到的光强
:
|
|
(7- 3) |
考虑在
到
的过程中光线按照按照
被吸收,其中:
|
|
(7- 4) |
称为传递函数。在
处发射光线其强度为
,则在眼睛看到的光强为:
|
|
(7- 5) |
考虑在
到
过程中,每一点都发射光线,则有:
|
|
(7- 6) |
为了能够在计算机上应用,将公式(7- 6)转化为数值计算方法:
设吸收方程为:
|
|
(7- 7) |
用Riemann和进行近似积分,可得:
|
|
(7- 8) |
带入
中,并进行变换可得:
|
|
(7- 9) |
定义不透明度
,可得:
|
|
(7- 10) |
将
离散后看成面积的乘积进行近似,其中c代表颜色值:
|
|
(7- 11) |
将公式(7- 10)带入(7- 11)可得:
|
|
(7- 12) |
公式(7- 12)可以如下进行从后向前积分的递归计算:
|
|
(7- 13) |
或从前向后递归计算:
|
|
(7- 14) |
此时应用公式(7- 13)和公式(7- 14),可进行GPU的光线投射计算。
传递函数
在体数据场中,体数据不直接提供颜色数据。如果要对体数据进行可视化处理,则必须要有一个体数据值和颜色的一个对应关系,这个对应关系称为传递函数。
传递函数的维度可以根据不同的需求定义。一般而言,1D传递函数提供了体数据值和颜色的一一对应关系,适合简单的表达。更复杂的2D传递函数则结合了体数据值、梯度和颜色的关系进行建立,可表达更为丰富的信息。
在多轴数控仿真的应用中,应用1D传递函数就足够表达。
GPU实现
GPU上实现光线投射算法如图 7-7所示,成像平面决定了最终成像的位置。光线从体数据中经过离散的采样(小圆圈代表一次采样),根据公式(7- 14)得到的最终积分结果作为成像平面对应像素的值。
当在体数据中有传统多边形遮挡体时候,体数据的渲染和遮挡体的渲染需要进行特殊处理,增加对遮挡体的位置的计算,以便对体数据实现正确的采样和最终的合成。
|
|
|
图 7-7 GPU实现原理 |
设置成像平面
在实际应用中,OpenGL中摄像机的位置即是眼睛的概念。该摄像机位置位于成像平面之后,也即在屏幕之后,无法被观察到。成像平面由OpenGL创建的窗口得到。此窗口的分辨率等于从成像平面发出的光线数量的多少。在OpenSceneGraph中,封装了各个系统下的OpenGL窗口的创建,仅需要以下代码即可设置好窗口,且代码可以跨平台:
|
osg::Vec2 screenSize(800, 600); //窗口分辨率
osg::ref_ptr<osg::GraphicsContext::Traits> traits = new osg::GraphicsContext::Traits; traits->width = screenSize.x(); traits->height = screenSize.y(); traits->windowDecoration = true; //窗口有边框 traits->doubleBuffer = true; //启动双重缓冲,防止闪烁
osg::ref_ptr<osg::GraphicsContext> gc = osg::GraphicsContext::createGraphicsContext(traits.get()); viewer->getCamera()->setGraphicsContext(gc.get()); viewer->getCamera()->setViewport(0, 0, screenSize.x(), screenSize.y()); //设置视口 |
得到包围盒
发出光线的动作在GPU上通过渲染整个成像平面得到。GPU会对成像平面上每个像素执行相关的着色器代码。光线投射算法中最关键的算法集中在片段着色器中。由于显示流水线的限制,而要使片段着色器能够进行积分工作,必须要在GPU管线中存在一个几何物体。通过计算体数据包围盒的方式可以得到一个长方体形状的几何图形。此图形将作为后续渲染环节中重要的几何物体,提供体数据的位置、前表面、后表面的相关数据。
体数据中的采样采用的是和贴图一样的访问形式,使用的坐标是归一化的坐标。对长方体的八个顶点进行贴图坐标赋值(如图 7-8),对应于体数据的归一化坐标。经过GPU固定管线进行插值之后,各个平面上相关点的贴图坐标就是对应的贴图坐标。
|
|
|
图 7-8 包围盒颜色赋值 |
前表面坐标、后表面坐标、场景物体坐标
为了进行积分操作,首先需要得到光线的路径。在光线投射算法中仅考虑光线沿直线进行的情况,所以得到光线进入包围盒的前表面坐标和后表面坐标即可得到光线的路径。
前表面和后表面坐标存储在GPU显存中,不需要进行显示。应用中设置摄像机为RTT相机:
|
osg::ref_ptr<osg::Camera> RTTCamera = new osg::Camera(); RTTCamera->setClearColor(osg::Vec4(0, 0, 0, 0)); RTTCamera->setClearMask(GL_COLOR_BUFFER_BIT| GL_DEPTH_BUFFER_BIT); RTTCamera->setReferenceFrame(osg::Transform::RELATIVE_RF); //设置相对于主相机的变换 RTTCamera->setRenderOrder(osg::Camera::PRE_RENDER); //设置为预渲染 RTTCamera->setRenderTargetImplementation(osg::Camera::FRAME_BUFFER_OBJECT); //设置渲染到贴图 |
在预渲染过程中,RTT相机将首先渲染包围盒,将包围盒位置和贴图坐标分别存储在贴图中。为了得到位置,必须要对顶点变换做处理。负责进行顶点变换的顶点着色器采用以下代码:
|
varying vec4 pos; varying vec4 tex;
void main() { gl_Position=ftransform(); //顶点进行固定管线转换 pos=gl_ModelViewMatrix*gl_Vertex; //模型视图坐标 tex=gl_MultiTexCoord0; //得到贴图坐标 } |
GPU固定渲染管道中的光栅化单元会对像素进行插值,每个像素会赋值为插值后的数据。负责像素着色的片段着色器采用以下代码:
|
varying vec4 pos; varying vec4 tex;
void main() { gl_FragData[0]=pos; gl_FragData[1]=tex; } |
这样,位置数据就保存在RTT的第一张贴图中,贴图坐标数据保存在RTT的第二张贴图中。
为了得到前后两个表面的数据,需要对包围盒进行两次不同的渲染。
第一遍渲染依照常规渲染,将包围盒的前表面坐标保存在内存中。
第二遍需要对包围盒进行前表面剔除,仅保留后表面,并且深度缓冲区要设置为GREATER使得Z值大的像素进行保留。
|
osg::ref_ptr<osg::StateSet> stateSet = renderBackCamera->getOrCreateStateSet(); stateSet->setMode(GL_CULL_FACE, osg::StateAttribute::ON); //设置OpenGL状态为表面剔除模式 stateSet->setAttributeAndModes(new osg::CullFace(osg::CullFace::FRONT), osg::StateAttribute::OVERRIDE | osg::StateAttribute::ON); //剔除前表面
stateSet->setMode(GL_DEPTH_TEST, osg::StateAttribute::ON); stateSet->setAttributeAndModes(new osg::Depth(osg::Depth::GREATER), osg::StateAttribute::ON | osg::StateAttribute::OVERRIDE); //设置深度缓冲区,保留Z值较大像素。 |
对场景的渲染进行一次,关闭包围盒的渲染,渲染场景物体的前表面位置到贴图中。
|
|
|
|
|
|
(a) 前表面位置贴图 |
(b) 后表面位置贴图 |
(c) 前表面贴图坐标 |
(d) 后表面贴图坐标 |
|
图 7- 9 渲染结果 |
|||
计算光线积分
在7.3.1一节中,已经创建了成像平面,将摄像机设置为2D模式,将一个和成像平面相同大小的长方形渲染到成像平面上,就可以模拟光线的发射。对于每一个像素,GPU会依据着色器的程序进行计算。
首先GPU要对顶点进行变换,由于不需要将长方形进行任何特殊变换,所以按照OpenGL的固定管线进行转换:
|
void main(void) { gl_Position=ftransform(); //调用固定管线计算方法计算 } |
然后GPU对平面进行片段着色。在光线投射算法中,核心的算法和主要的工作量都在这个着色器中进行。下面对片段着色器进行分段描述。
步骤1:设置贴图
|
uniform sampler1D transferFunction; //颜色传递函数 uniform sampler2D texFrontPos; //前位置贴图 uniform sampler2D texFrontTex; //前贴图坐标 uniform sampler2D texBackPos; //后位置贴图 uniform sampler2D texBackTex; //后贴图坐标 uniform sampler2D texScenePos; //场景物体位置 uniform sampler3D volumeData; //体数据 uniform float transparency; //透明度 uniform float sampleDensity; //采样密度 uniform vec4 volScale; //体伸缩系数 uniform vec2 renderSize; //渲染窗口大小 void main() { |
步骤2:得到光线长度
|
vec2 screenTex = gl_FragCoord.xy / renderSize.xy; //得到屏幕的坐标的归一化坐标
//计算光线长度 float rayLength; vec4 sceneFrontPos = texture2D(texScenePos, screenTex); //从场景贴图中获取物体前表面位置 vec4 frontPos = texture2D(texFrontPos, screenTex); //从包围盒前表面贴图中得到前表面位置 vec4 backPos = texture2D(texBackPos, screenTex); //从包围盒后表面贴图中得到后表面位置
//计算前表面到后表面的距离和前表面到场景的距离 float frontToBack = distance(frontPos, backPos); float frontToScene = distance(frontPos, sceneFrontPos); //针对有遮挡的情况,取距离最小的值作为光线的长度。 rayLength = (sceneFrontPos.a == 0) ? frontToBack : min(frontToBack, frontToScene); |
步骤3:计算光线方向
|
vec3 rayDirection; vec4 frontTex = texture2D(texFrontTex, screenTex); //前表面贴图位置 vec4 backTex = texture2D(texBackTex, screenTex); //后表面贴图位置 //由于都是归一化的数据,前表面减去后表面再进行归一化可得到向量的方向。 rayDirection = normalize(frontTex - backTex).xyz; |
步骤4:对光线长度进行伸缩
由于光线是在体数据的归一化坐标系中进行穿行,要将光线的长度转化为体数据坐标系下的长度。
|
rayLength = rayLength / length(rayDirection * volScale.xyz); |
步骤5:计算循环次数和采样递增量
|
vec3 rayPosition = frontTex.xyz; //光线起点 vec3 deltaRayPosition = rayDirection / sampleDensity; //光线采样递增量 rayPosition -= deltaRayPosition * fract(sin(dot(screenTex.xy, vec2(12.9898, 78.233))) * 43758.5453); //对起点进行微小抖动处理,减少最终成像呈等高面的问题 int stepCount = sampleDensity * rayLength; //采样个数 |
根据Shannon 采样定理,采样率要保持在体数据的两倍以上。所以sampleDensity一般设置为1/(2*体数据最大维度)
步骤6:循环积分
|
for (int i = 0; i < stepCount; i++) { vec4 vol = texture3D(volumeData, rayPosition); //体数据 vec4 src = texture1D(transferFunction,vol.a); //应用传递函数,取得该点的颜色值 src.a = src.a * transparency; //透明度调整 //根据公式(7- 14)进行颜色和不透明度的计算,应用GLSL内置函数提速 dst.rgb = mix(src.rgb, dst.rgb, dst.a); dst.a = mix(src.a, 1.0, dst.a); if (dst.a >= 1) break; //提前退出 rayPosition -= deltaRayPosition; //向前一个采样点 } |
步骤7 :屏幕显示
|
gl_FragColor = dst; //将当前片段颜色设置为最终颜色 |
体数据渲染得到图形必须要和传统光栅化的图形进行融合显示,将摄像机设置为融合状态:
|
cameraStateSet->setMode(GL_ALPHA_TEST, osg::StateAttribute::ON | osg::StateAttribute::OVERRIDE); cameraStateSet->setMode(GL_BLEND, osg::StateAttribute::ON | osg::StateAttribute::OVERRIDE); |
融合后体渲染的数据和光栅化数据就同时显示在屏幕上了。
图形仿真结果及演示
根据前面对包络面生成的算法并结合本章的光线投射算法,通过生成例子验证本方案的可行性。
|
|
|
图 7- 10 球铣刀加工包络面 |
图 7- 10展示了球形铣刀在空间进行五轴联动时生成的包络面的情形。其运动过程中包含了A,B轴的转动以及XYZ的移动,为典型的五轴联动运动。
图 7- 11展示了该包络面与工件进行同时进行显示的情形。
图 7- 12展示了包络面与工件进行布尔运算后得到的加工工件的情形。可见工件表面被切削,生成了自由曲面。
|
|
|
|
图 7- 11 球铣刀加工包络面 |
图 7- 12 球铣刀任意路径切割显示 |
通常,加工中心支持多把刀具的加工。图 7- 13展示了两把刀分别进行钻孔和铣削的过程,证明该仿真系统能够适应多种刀具加工。
|
|
|
|
(a)第一把刀进行钻孔 |
(b)钻孔成型结果 |
|
|
|
|
(c)球铣刀加工 |
(d)最后成型 |
|
图 7- 13 两把刀加工演示 |
|
最后,图 7- 14展示了叶轮数据体素化后经过光线追踪算法并进行场景融合显示的实例。
|
|
|
图 7- 14 叶轮体数据和简化机床模型 |








































方向分别布置1个视线方向朝向物体的摄像机。各个摄像机正交布置,获取到6个正交方向的光栅化的图像。这些图像记录了从不同方向获取的位置信息。要判断离散空间一个点
是否是位于物体表面或者物体内一点,仅需要满足:

所在的数据设置为1表示该粒子占据空间,否则设置为0。


,的值。再根据公式(5-1)可得到物体在三维空间离散的结果。

在GPU上对应的线程如图 5-3所示。在CUDA编程模型中,可以用如下代码获得当前线程所对应的点
:


。将t标准化为
后,将各时段
相加,可表示为:

表示实体扫描体G的边界,
表示第k个面,
和
分别表示在扫描起始位置和最终位置的轮廓,则可分解为另一种表达:



,
足够小,并且在
和
之间存在交线
,则上式可以用下式代替:

,方程(4- 5)可写为如下的微分形式:



,则点P必须满足下面条件:

为运动面F在点P的法线,该等式称之为包络点的运动条件。
