GPU光线投射算法

Ray Tracing（光线跟踪）算法是标量场的一种直接可视化方法。它从将屏幕作为光线接收的一个平面，根据光路可逆原理，从屏幕各个像素上发出若干光线并模拟在体数据中的投射、反射、折射，得到光线的积分，作为屏幕像素的值。由于其模拟光线的物理原理，可以得到很高的图像质量，广泛应用于电影、高质量动画的制作。虽然其原理简单，但由于需要计算大量的光线，其运算量很大，要在现有的硬件条件下达到中等解析度下实时的速度都较为困难。

Ray Casting（光线投射）算法作为光线追踪算法的一个简化，仅考虑光线在体中的直线运动，不考虑折射和反射。其运算较为简单，成像质量能够达到数控仿真的要求。

基本原理

图 7-5是光线投射算法的基本原理。所有的光线从光线出发点出发，经过成像平面，再经过体数据。所有的光线都在视锥之内。考虑将光路逆反，从体数据穿过，穿过时根据传递函数（Transfer Function）得到体数据在离散点上的颜色和透明度，并对其进行积分得到最终的颜色值。所有的光线经过成像平面后，最后得到体数据的光线投射算法的图像。

在光线追踪算法中，有两种重要的辐射物理模型：发散和吸收。发散模型从一个点向外发射一束光线或者吸收能量后再进行发射一束光线。发射模型增强总光线强度。吸收模型从一个点接收光线，然后经过衰减后再发射一束光线。吸收模型降低总光线强度。这两种物理模型均可以包含散射的物理模型。

为了简化处理，光线投射模型仅考虑无散射和无折射的物理模型。

光线积分

设初始光强为，如果不考虑光线的衰减，则可以得到眼睛看到的光强：

考虑在到的过程中光线按照按照被吸收，其中:

称为传递函数。在处发射光线其强度为，则在眼睛看到的光强为：

考虑在到过程中，每一点都发射光线，则有：

为了能够在计算机上应用，将公式（7- 6）转化为数值计算方法：

设吸收方程为：

用Riemann和进行近似积分，可得：

带入中，并进行变换可得：

定义不透明度，可得：

将离散后看成面积的乘积进行近似，其中c代表颜色值：

将公式（7- 10）带入（7- 11）可得：

公式（7- 12）可以如下进行从后向前积分的递归计算：

或从前向后递归计算：

此时应用公式（7- 13）和公式（7- 14），可进行GPU的光线投射计算。

传递函数

在体数据场中，体数据不直接提供颜色数据。如果要对体数据进行可视化处理，则必须要有一个体数据值和颜色的一个对应关系，这个对应关系称为传递函数。

传递函数的维度可以根据不同的需求定义。一般而言，1D传递函数提供了体数据值和颜色的一一对应关系，适合简单的表达。更复杂的2D传递函数则结合了体数据值、梯度和颜色的关系进行建立，可表达更为丰富的信息。

在多轴数控仿真的应用中，应用1D传递函数就足够表达。

GPU实现

GPU上实现光线投射算法如图 7-7所示，成像平面决定了最终成像的位置。光线从体数据中经过离散的采样（小圆圈代表一次采样），根据公式（7- 14）得到的最终积分结果作为成像平面对应像素的值。

当在体数据中有传统多边形遮挡体时候，体数据的渲染和遮挡体的渲染需要进行特殊处理，增加对遮挡体的位置的计算，以便对体数据实现正确的采样和最终的合成。

设置成像平面

在实际应用中，OpenGL中摄像机的位置即是眼睛的概念。该摄像机位置位于成像平面之后，也即在屏幕之后，无法被观察到。成像平面由OpenGL创建的窗口得到。此窗口的分辨率等于从成像平面发出的光线数量的多少。在OpenSceneGraph中，封装了各个系统下的OpenGL窗口的创建，仅需要以下代码即可设置好窗口，且代码可以跨平台：

得到包围盒

发出光线的动作在GPU上通过渲染整个成像平面得到。GPU会对成像平面上每个像素执行相关的着色器代码。光线投射算法中最关键的算法集中在片段着色器中。由于显示流水线的限制，而要使片段着色器能够进行积分工作，必须要在GPU管线中存在一个几何物体。通过计算体数据包围盒的方式可以得到一个长方体形状的几何图形。此图形将作为后续渲染环节中重要的几何物体，提供体数据的位置、前表面、后表面的相关数据。

体数据中的采样采用的是和贴图一样的访问形式，使用的坐标是归一化的坐标。对长方体的八个顶点进行贴图坐标赋值（如图 7-8），对应于体数据的归一化坐标。经过GPU固定管线进行插值之后，各个平面上相关点的贴图坐标就是对应的贴图坐标。

前表面坐标、后表面坐标、场景物体坐标

为了进行积分操作，首先需要得到光线的路径。在光线投射算法中仅考虑光线沿直线进行的情况，所以得到光线进入包围盒的前表面坐标和后表面坐标即可得到光线的路径。

前表面和后表面坐标存储在GPU显存中，不需要进行显示。应用中设置摄像机为RTT相机：

在预渲染过程中，RTT相机将首先渲染包围盒，将包围盒位置和贴图坐标分别存储在贴图中。为了得到位置，必须要对顶点变换做处理。负责进行顶点变换的顶点着色器采用以下代码：

GPU固定渲染管道中的光栅化单元会对像素进行插值，每个像素会赋值为插值后的数据。负责像素着色的片段着色器采用以下代码：

这样，位置数据就保存在RTT的第一张贴图中，贴图坐标数据保存在RTT的第二张贴图中。

为了得到前后两个表面的数据，需要对包围盒进行两次不同的渲染。

第一遍渲染依照常规渲染，将包围盒的前表面坐标保存在内存中。

第二遍需要对包围盒进行前表面剔除，仅保留后表面，并且深度缓冲区要设置为GREATER使得Z值大的像素进行保留。

对场景的渲染进行一次，关闭包围盒的渲染，渲染场景物体的前表面位置到贴图中。

计算光线积分

在7.3.1一节中，已经创建了成像平面，将摄像机设置为2D模式，将一个和成像平面相同大小的长方形渲染到成像平面上，就可以模拟光线的发射。对于每一个像素，GPU会依据着色器的程序进行计算。

首先GPU要对顶点进行变换，由于不需要将长方形进行任何特殊变换，所以按照OpenGL的固定管线进行转换：

然后GPU对平面进行片段着色。在光线投射算法中，核心的算法和主要的工作量都在这个着色器中进行。下面对片段着色器进行分段描述。

步骤1：设置贴图

步骤2：得到光线长度

步骤3：计算光线方向

步骤4：对光线长度进行伸缩

由于光线是在体数据的归一化坐标系中进行穿行，要将光线的长度转化为体数据坐标系下的长度。

步骤5：计算循环次数和采样递增量

根据Shannon 采样定理，采样率要保持在体数据的两倍以上。所以sampleDensity一般设置为1/(2*体数据最大维度)

步骤6：循环积分

步骤7 ：屏幕显示

体数据渲染得到图形必须要和传统光栅化的图形进行融合显示，将摄像机设置为融合状态：

融合后体渲染的数据和光栅化数据就同时显示在屏幕上了。

图形仿真结果及演示

根据前面对包络面生成的算法并结合本章的光线投射算法，通过生成例子验证本方案的可行性。

图 7- 10展示了球形铣刀在空间进行五轴联动时生成的包络面的情形。其运动过程中包含了A,B轴的转动以及XYZ的移动，为典型的五轴联动运动。

图 7- 11展示了该包络面与工件进行同时进行显示的情形。

图 7- 12展示了包络面与工件进行布尔运算后得到的加工工件的情形。可见工件表面被切削，生成了自由曲面。

通常，加工中心支持多把刀具的加工。图 7- 13展示了两把刀分别进行钻孔和铣削的过程，证明该仿真系统能够适应多种刀具加工。

最后，图 7- 14展示了叶轮数据体素化后经过光线追踪算法并进行场景融合显示的实例。

3D空间离散点的值可以视为一个3D标量场。标量场可视化的经典算法包括Marching Cubes算法（简称MC算法）和光线投射（Ray-casting）算法。前者将标量场转换为多边形片进行显示，后者通过光线投射算法直接进行显示。

一直以来，研究人员都追求如何快速的进行MC算法的计算，并将视线由传统的CPU转变到了GPU上。早期由于图形硬件的限制，导致Marching Cubes算法很难应用到GPU领域。经过研究人员的不断努力，提出了简化的Marching Tetrahedron方法，但其成像速度及质量较差。在GPU越来越通用化以后，通过应用Render To Texture和GPU上的流式缩减，在GPU上实现MC算法[35]，并且达到实时运算的速度。7.1节中将重点阐述此方法的应用。

Ray Casting（光线投射）算法是Ray Tracing（光线追逐）算法的一个简化，它以屏幕像素作为光线的出发点，向标量场投射光线，并对光线经过的标量值进行积分，通过积分的结果决定像素的颜色、透明度并显示在屏幕。此算法如果利用CPU计算由于涉及到大量像素的操作，不具有实时性。所以大量的算法都基于GPU的计算。

下面将分别介绍两种算法并展示其在五轴数控仿真方面的应用。

GPU Marching Cubes算法

目前基于GPU的算法中，由于可以将所得到的结果直接进行渲染，不必经过CPU和GPU之间的交互，所以利用Histogram Pyramid（简称HistoPyramid）在GPU上进行流式缩减和扩展的算法能够在得到最高的性能[35]。

Marching Cubes算法

在标量场中提取等值面的算法中，文献[36]提到的Marching Cubes算法是目前应用最多的算法。它将一个大小为H*W*D的3D标量场表示为(H-1)*(W-1)*(D-1)个MC单元（图 7-1）组成的立体网格。每个单元置于8个标量值之中，其值可以通过算法应用需要对8个标量值进行处理。这样一个MC单元就有一个标量值，以及8个包围它的8个标量值。如果依次遍历这个立体网格，按照设定的阈值可以区分MC单元属于物体内还是外，再根据MC单元值和8个标量值的不同组合可以得到组成这个MC单元的多个三角形的形状。所有的三角形组合起来就可以得到标量场等值提取所表示的面。

图 7-1是一个典型的MC单元，根据阈值将8个标量值分为了里（bd）和外（acefgh）。MC单元根据8个标量值的处于内和外的不同，具有256种组合。这些组合可以通过基本的15种组合经过变换得到。由于每个单元的计算可以单独计算，所以MC算法并行性较高，能够应用GPU进行加速。

在图 7-1所示的这种MC类型中，平面和MC单元相交形成了一个四边形，阈值的取值和标量值的值决定了交点所在位置。例如采用线性插值，当阈值iso=0.4，c=0.8，d=0.32时，点P(c,d)在直线cd上，直线CD:CP=0.4。如果采用取中间点算法，则得到P(c,d)在直线cd的中点上。但这种取中点的算法没有用线性插值所得到的结果平滑，而且由于目前的GPU运算速度很快，线性插值的性能损失可忽略不计。

HistoPyramid

HistoPyramid是一种类似于三角金字塔的结构，其最底层结构为输入流，输入流中的一个元素代表了某种特定类型的数量，在金子塔的顶端，代表了所有输入流元素所有类型的数量的和，这个和决定了输出流元素的数量。

如果输入流中的元素代表了0个特定类型，则这个输入流元素就被丢弃，这样的生成输出流就被压缩了。如果输入流元素代表了1个以上的特定类型，则在输入流元素中会有1个以上的元素输出，这样就实现了数据的扩张。

HistoPyramid有两个不同的阶段。第一个阶段创建一个HistoPyramid，它的数据结构类似于MipMap，为金字塔型，其计算方式也类似于MipMap的计算，只是将MipMap的平均值算法改为了求和算法。图 7-2是真实显存中的数据。每一级的大小都是前一级的1/4，直到大小变为1*1为止。最底层的一级称为"基础层"，最上层的一级称之为"顶层"。

第二个阶段为生成输出流阶段。输出流的生成逆向的从"顶层"向下查找，直到最后在"基础层"中找到输入流中对应元素的索引后，在根据特定规则生成相应元素。

构造基础层

由于目前的显卡不能渲染到3D贴图，或者通过软件实现渲染速度很慢，所以在构造基础层时将提供的标量场3D贴图进行平坦化，使之成为一个标准的2D贴图。由于标量场的一个元素中仅存储有一个float数据，而显卡在处理2D贴图时按照传统图形应用的RGBA像素进行提取，一次性可以取得四个数据。所以将标量场中每四个相邻数据存储到一个RGBA像素中，可以大大节约存储空间，并且得到更高的效率。

每一个RGBA像素代表了2*2*1个MC单元，所以获得一个像素所代表的4个MC单元的分类，需要提取3*3*2个MC单元的边角的值。阈值的大小和MC单元边角的标量场的数据决定了MC单元的分类和所需三角形的数量。为了减少存储空间，由于最多需要15个顶点就可以表示一个MC单元，一个RGBA像素可以利用整数部分存储顶点数量，利用小数部分存储MC单元的类型。这样GPU的流式缩减操作可以将整数部分保留，进行加法，得到最终输出内的顶点的数量。当在基础层决定顶点的类型时，将小数部分提取出来即可。

流式缩减

执行流式缩减操作将得到顶点的数量。在基础层中已经计算好每个MC单元所拥有的顶点数量，将其相加即可得到所有的顶点。

在GPU上实现此操作需要进行多步。此操作类似于生成MipMap的计算：第n张贴图的边长是第n-1张贴图边长的1/4大小。第n张贴图的某一个像素点(x,y)的值由下式决定：

当最后一张贴图的边长为1*1时，流式缩减操作结束。由于将顶点数量存储在RGBA中，此时仍不能得到顶点总数量。此时，存储的RGBA值将会传回CPU，并将其RGBA值相加，得到总数量t。

下面以一个简化的HistoPyramid为例（图 7- 3）介绍流式缩减操作。 BaseLevel含有16个元素，将红框中的4个数据相加得到2，存入level=1的(0,0)位置，其他12个元素按照此规则进行。此后level=1的所有元素相加，得到顶层的数值10。如果以RGBA进行存储，则进行到level=1时就要将数据传回CPU，然后CPU将RGBA的值相加，得到最终元素个数10。

流式扩展

此例中计算得知最终的output中有10个元素。要得到这10个元素，首先根据得到的总元素个数t生成一个k值列表，范围从0到t-1以1为步进递增。

从k值得到最后的坐标和序号的过程叫做"遍历"。设当前k值为m，当前层数为l，顶层为t，当前像素位置为p。

初始时l=t，k=m ，p指向最顶层元素，此时只有一个元素。

l=l-1，p指向l-1层的中心位置，其附近有4个值，记这四个值为a,b,c,d。四个数值可以得到如下数值范围：

然后得到k所在的范围。例如如果k在C范围内，则我们将p指向c元素，并且将k减去范围C的起始范围，这里k=k-(a+b)。

重复过程2，直到l=0

此时p指向k所对应的输入元素的索引，k表示元素的第k+1个值。例如k指向的输入中包含两个元素R，S。当k=0时代表输出元素中的元素是R，k=1时代表S。

以图 7- 3为例，l=t，k=4时，第一次从顶层出发，k=4不变，p指向level=1的中心。此时a=2,b=3,c=1,d=4。得到4个范围：

A=[0,2],B=[2,5),C=[5,6),D=[6,10)

k=4落到B区间中。此时k=4-2=2，l=t-1=1，p指向元素3。

p在level=0的图中得到a=0,b=1,c=0,d=2，此时得到4个范围：

A=[0,0),B=[0,1),C=[1,1),D=[1,3)。

k=2落到区间D中。此时k=2-1=1，l=l-1=0，p指向元素2

此时l=0，结束遍历。剩余的k=1表示输出元素k=4时输出输入元素坐标（3,1）中的第二个元素。

通过计算k=3和k=4的值可以发现，其都指向输入元素坐标（3,1），但其输入元素序号不同，这样就实现了流式扩展。

注意输入元素中值为0的部分在输出元素中并没有被输出，这样就实现了数据的压缩。

HistoPyramid应用于MC算法

HistoPyramid应用重点在于基础层的结构。在MC算法中，涉及到MC单元的顶点数量和单元类型。将浮点型数据的整数部分存储顶点数量，小数部分存储单元类型即可。建立好基础结构以后，OpenGL执行多遍渲染以便执行流式缩减，将基础层的整数部分相加，得到三角形的数量总和t。此后，CPU根据t创建一个GPU输入的缓冲区存储k值，大小即为t。输入缓冲可以用Vertex Buffer Object实现。此外，还需要创建一个大小为t的输出缓冲区。

GPU将对缓冲区的每个k值进行遍历操作。每一次遍历完成后会得到MC单元位置和顶点索引。这时再查找基础层里面对应MC单元存储的单元类型，进行查表。MC单元的类型及顶点数量存储在一个256*16的2D贴图中，其中横坐标的256代表256种类型，16代表15种MC基础类型。由于GPU在存储时按照4的整数倍存储效率最高，所以用16来代表15种基础类型，剩余一个可忽略。这样就可生成MC单元的顶点，根据阈值通过插值可以得到顶点坐标。

得到的顶点坐标可直接应用于GPU的三角形渲染，只需要将顶点连接起来，并赋予相应的颜色值和法线坐标（应用于光照）即可。图 7-4是使用MC算法生成的复杂曲面的包络面，体数据分辨率256*256*256，8800GTS上可达40fps以上。

作为扫描体组成的一部分，包络面的生成可根据其数学表达式或几何意义进行计算。文献[24]结合离散造型法和实体造型法这两种造型方法的优点，对刀具进行离散采样的预处理，然后将刀具上的离散采样点沿着扫描路径运动，生成扫描元曲线。再在扫描路径上建立起来的若干截平面与扫描元曲线相交后得到的交点集，构造平面网格，求解出平面网格的外轮廓顶点后将其分组，分别进行曲线拟合后再组合成刀具扫描体的包络曲线。通过这些包络曲线进行放样，构建刀具扫描体的包络曲面，最终完成刀具扫描体的完整构造。其算法新颖，所得包络面质量较高，但计算颇为复杂，特别是涉及自交情况时更为复杂。文献[25]中脱离了传统利用曲面族包络面方程的计算，提到了一种针对五轴数控刀具的简化的运算方法。该算法将刀具分为入点集、出点集和切点集，并对切点集的运动进行计算。文献[26-30]都使用相似的思路对路径进行离散，采用小直线段进行逼近，然后生成包络面的曲面，对其进行体素化再和工件进行布尔计算。

根据包络面形成的几何意义，是物体在三维空间中移动所占用的空间。早期对包络面生成的算法主要集中在基于体积的扫描的算法。但由于随着精度的提高，早期计算机没有相应的存储空间并且运算速度较低，导致算法无法大量应用。文献[31-33]均使用移动立方体的方法进行运算，其运算量大，但由于不受视点的限制，算法并行性较高，多在高性能计算机上有应用。文献[34]将算法简化为在单视点上进行操作，减少了计算和存储量，但移动视点会导致重新计算。

本文将根据工程需要，介绍两种不同的算法，分别是直接移动立方体的算法和光栅化包络面的方法。前者算法较为简单，工作量主要集中在对显存的存储，需要对算法进行大量细致的优化才能达到较高的效率并且由于离散化后精度不高，进行旋转和平移操作后会丢失掉部分包络面实体，精度较差。后者将工作量集中在3D显示大量刀具移动体拼合的包络面上，对显卡的几何运算能力要求较高，但生成的离散体更为精确。

移动刀具离散体生成包络面离散体

由GPU计算得到的刀具离散体存储在GPU的显存中。如果要计算包络面离散体，则对所有的体素都进行离散路径上的位置变换并对相应空间进行赋值即可。

考虑将刀具的一个横截面向x正方向进行移动，刀具的横截面如图 6- 1右边浅色圆形所示。要得到包络体的截面，仅需要将离散后的刀具进行位置变换即可。

但移动一个完整的刀具截面所需要的内存写入量较大，并且移动前后需要重复操作大量的数据。考虑到包络面的成型仅需要物体表面离散体经过移动形成，所以在得到刀具离散体时可以将内部体素去掉，仅留下对包络面成型有影响的部分。如图 6- 1左边的浅色圆圈表示刀具外围的体素，亮色圆弧代表其在移动过程中对包络面离散体产生的贡献。

在CUDA实际计算中，将经过处理后的刀具外表面离散体视为一个一维数组，数组元素由体素坐标（x,y,z）和体素值决定。传递给CUDA的参数包含刀具中心位置坐标（cx,cy,cz），以及刀具起始坐标及姿态（tx0,ty0,tz0,a0,b0）和刀具终点坐标及姿态（tx1,ty1,tz1,a1,b1）。CUDA核函数进行移动的插值计算。由于该算法实用性不高，在此不列出相关算法。

光栅化刀具包络面生成离散体

文献[27-30]将包络面的生成按照特定的刀具分类，对每种刀具都应用特定的算法来计算包络面。但由于算法没有充分的考虑到刀具的多样性，所以每增加一种刀具就要增加一种算法，对于特殊的刀具路径还需要做进一步的处理。

由于GPU仅根据所得到的6个面的离散位置信息来确定最终体素的值，而离散贴图的精度有限，如果初始位置和末位置距离足够近，则最终产生的表面在特定分辨率下为一个连续表面。所以可以离散的移动刀具的几何实体来产生在贴图上连续的表面，避免包络面生成的复杂数学运算，并且可以适应任何路径。根据此思路，图 6- 2展示了离散步数和生成的表面之间的关系，图 6- 3展示了复杂表面的生成，此生成的表面将用于GPU生成扫描体并应用于下一步的计算中。

 

概述

体素化是一种用一系列的体素（voxel）逼近连续体的方法[20]。这些体素将被三维光栅化赋值，并应用到后期的处理过程中[21]。

文献[22]提到使用一对深度缓冲来进行体素化的方法。但这种算法仅能适用于凸多面体，适用面不广。此算法后来被文献[23]中得到改进。它使用GPU从6个方位渲染，通过得到Z-buffer来进行计算。多个Z-buffer能够适用于凸多面体和大多数非多面体，并且充分利用了硬件加速。但由于进一步的体生成算法在CPU上进行，需要将得到的图形从GPU上传回CPU，而且CPU负责大量的赋值运算，故整体速度并没有得到提升。

本文延续以上的思路并进行了改进。首先利用OpenGL Shader渲染几何体得到6个方位的表面位置信息和表面的法线信息。然后利用GPU的高运算速度和存储带宽，使用CUDA在GPU上直接进行离散体的计算，生成的离散体直接存储在GPU的内存中。生成的离散体经过转化为3D纹理可直接应用于基于GPU的Marching Cubes算法或者光线投射算法。整个算法都集中在GPU上进行处理，省去了GPU到CPU的数据转移，大大增强了算法的性能。

获得正交光栅化贴图

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

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

在OpenGL中摄像机有透视摄像机和平行摄像机两种。透视摄像机的视域体是锥形，通过其成像会得到近大远小的图形，符合人眼的视觉；平行摄像机是平行正六面体视域体，得到的图像是与摄像机和物体距离无关的图形，通常适合于CAD领域。

由于需要得到与物体距离无关的图形，并且要保证像平面的像素和空间像素的一一映射关系，摄像机必须要首先设置成平行摄像机模式。在OSG中，我们应用osg::Camera类的setProjectionMatrixAsOrtho2D来来初始化摄像机。以X+方向的摄像机为例：

这样就设置好了一个左下角坐标为(0,0)，右上角为(1,1)的摄像机。其他方向以此类推。

为了让GPU能够存储图像，6个摄像机的图像不能直接输出到屏幕，而是利用OpenGL的RTT（Render To Texture）技术将图像渲染到缓冲区中。此缓冲区连接到一个2D贴图，使得离散化时能够从这些2D贴图中直接获取数据。在OpenGL 2.0规范中，能够使一个摄像机连接多个2D贴图，使得可以通过一次渲染获得多种数据。

在OSG中，可以如下设置摄像机将摄像机连接到存储离散点位置信息的maskTexture和存储法向的normalTexture。

摄像机的windowSize属性设置了像平面像素大小，该像素大小直接决定了最后离散体的精度。maskTexture和normalTexture的大小必须要和摄像机的windowSize相符合。为了能够得到高精度的maskTexture和normalTexture，必须要将贴图的属性设置为32位浮点类型，并用RGB来存储x,y,z值。

在文献[23]中应用了Z-buffer进行光栅化，但由于Z-buffer仅仅记录了深度信息，无法获取关于被离散物体更多的信息。为了增加灵活性以及脱离传统的固定管线，所以本文引入Vertex Shader（顶点着色器）和Fragment Shader（片段着色器）进行空间位置的转换和最终位置的写入。

为了能计算空间位置，必须要将顶点进行适当的转换进而得到顶点在世界坐标系下的空间位置。世界坐标系是一个应用程序的概念，出现在OpenGL物体坐标系和眼坐标系之间。如果将模型的视点矩阵分为模型变换和视点变换两次，则对物体空间坐标进行模型变换将得到世界坐标，对世界坐标进行视点变换将得到眼坐标。由于OpenGL不在世界坐标系中执行任何操作，所以正规的OpenGL变换流水线中没有这种坐标系。不过，我们依然可以人为的将转换分为两步得到顶点的世界坐标系。

根据世界坐标系的概念，我们首先需要得到摄像机的视点矩阵的逆，并与模型视点矩阵进行相乘，得到模型矩阵，即可获知世界坐标系。

Vertex Shader中，通过camera->getInverseViewMatrix()函数得到ViewMatrix的逆并进行赋值。

顶点转换处理完成后，通过硬件光栅化执行插值处理后，Fragment Vertex将逐点进行存储。

经过一次渲染过程以后，6个摄像机将会得到12张贴图。如果不需要法向信息，可将normalTexture相关的代码省略。

GPU渲染后得到12张离散信息图存储在显存中(图 5-2)，将用于下一步计算3D离散体。

CUDA计算3D离散体

在不考虑法线的情况下，通过6个摄像机的离散图得到，的值。再根据公式（5-1）可得到物体在三维空间离散的结果。

在CPU上实现的伪代码如下

以上经过优化的CPU伪代码使用了嵌套循环，在CPU上实现该函数效率较低，填充一个典型的256*256*256的离散体需要数百毫秒的时间。

对于GPU来说，GPU可以利用其超强的并行性模拟嵌套for循环，对离散体的每一个粒子进行赋值。

下面以256*256*256的离散体作为例子，说明该算法在CUDA上的实现。

CUDA编程中，两个重要概念是gridDim和blockDim，通过它们的不同组合，可以实现各种类似于for循环的算法。如果我们将blockDim设置为一维，即blockDim.y=1，blockDim.x=256，gridDim设置为2维，即gridDim.x=256，gridDim.y=256，则可以模拟三层嵌套循环为立方体离散体赋值，如图 5-3。

GPU在运行kernel时，会出发大量线程。每一个block内可以有多达720个线程并发，每一个SM上可以同时运行多个block。在上述立方体离散赋值模型中，空间一个点在GPU上对应的线程如图 5-3所示。在CUDA编程模型中，可以用如下代码获得当前线程所对应的点：

下面，为了简化叙述，将6个方位的贴图中的一个点视为如下结构：

整幅贴图视为结构ImageData的一个集合，在C语言中可定义为一个一维数组：ImageData *

离散体的一个粒子由如下结构定义：

整个离散体可视为VoxData的一个集合：VoxData *

则可由6个贴图和需要赋值的离散体得到GPU kernel的函数声明：

这里为kernel函数传入了6个方向的贴图的指针，和最终写入的离散体的指针。kernel函数并行的在GPU内执行，根据threadId和blockId可以最终确定当前线程所需要写入离散体的位置。要得到当前点对应的6个方向贴图的数据，仅需要访问各个对应ImgDevice数组内的值即可。例如需要访问x+方向的数据，则执行以下数组访问即可：

所以，根据上述描述和对应CPU的伪代码，可得下列GPU kernel代码：

注意，由于GPU乱序执行会大大降低运算性能，故这里没有利用和CPU伪代码一样的预判断zn和zp是否等于0的代码，而是直接取得所有的值。

在主程序中，CPU需要调用GPU的代码运行。CUDA大大简化了其调用过程。一个256*256*256的体可以调用如下代码进行计算：

性能

我们在VC++2005环境下对多个模型进行GPU上的的离散化，并测试其性能。测试所用硬件平台为：

CPU：Intel Core 2 Duo Q6600，频率2.4G

GPU：NVIDIA 8800GTS 320M

操作系统：Windows 7

测试一

叶轮由共计31934个顶点，32697个三角片。体素化分辨率为256*256*256，体素化后有效体素为16%。

测试完全由GPU生成6个方向的离散位置贴图和法向贴图，分辨率均为256*256，平均帧率保持在95fps左右。如去掉法向贴图，平均帧数可保持在160fps左右。每一个贴图平均耗时0.98ms，生成6个方向贴图耗时5.88ms。

如果将叶轮替换为1*1*1的正方体盒子，则可得知由GPU光栅化三角片带来的性能损失。此时生成12个贴图的平均帧率将保持在90fps左右。由叶轮多边形带来的性能损失在5%以下，可以忽略不计。

在CPU性能测试中，将6个贴图传回主内存后，然后测试算法性能。由6个贴图生成最终的离散体的CPU算法中，使用和GPU类似的循环算法，需要0.277s，使用优化后的代码跳过空白区域，需要0.124s。

在GPU性能测试中，6个贴图始终保持在GPU内存。生成离散体需要26.98ms。

测试二

直径0.5的球型，由3DSMAX创建，细分值64，离散化后有效空间占51.9%。

从6个贴图生成最终的离散体，非优化算法CPU耗时320ms，优化算法消耗325ms。优化算法由于粒子有效空间太大，加上if判断语句的性能损失，反而比非优化算法性能差。

GPU算法消耗47.48ms。

测试三

1*1*1的立方体，离散化后有效空间占100%

非优化算法CPU耗时373ms，优化算法657ms。也是由于if语句性能损失造成算法性能差。

GPU算法消耗73.59ms。对比测试二，其时间几乎消耗在对体数据的赋值上。

从以上测试可以看出：

简单的GPU算法的效率要比CPU效率高出数倍。如果考虑从GPU传送贴图到CPU上进行处理的时间，并且优化内存读取的方法，其效率还会更高。

GPU算法的时间消耗和离散化后有效空间比率成正比关系，时间消耗在GPU对显存的操作。

在铣床加工过程中，工件装在工作台上或分度头等附件上，铣刀旋转为主运动，辅以工作台或铣头的进给运动，工件即可获得所需的加工表面。如果将工件看成固定的物体，其所在空间位置固定不变，将工作台面的运动转换为刀具的移动，移动过程中所占用的位置就是刀具的扫描体。

在加工之前，需要根据工件的形状进行建模，由专门的数控软件根据工件形状和最终生成的零件生成相应的刀具路径，然后根据刀具路径生成相应的G代码。五轴联动系统具有五个自由度，根据不同的型号的机床，可以进行各种形式的多轴联动，加工极为复杂的表面。

但正是由于高自由度和各种类型机床的不同形式，软件生成的加工路径可能造成机床部件和加工件的干涉，或者刀具与机床部件的干涉，导致加工无法进行。而现有的路径规划软件都未考虑刀具夹具等因素。所以，仿真系统模拟机床的加工过程就显得尤为重要。

一个典型的数控系统由代码翻译为机床各轴的运动进行驱动。仿真机床也由同样的G代码进行驱动。由于各种机床的G代码有所不同，需要进行定制开发，并且G代码翻译成刀位点的算法已经很成熟，所以本文不再进行G代码翻译的讨论。

G代码翻译后经过离散插值可得到的刀位信息可控制仿真机床的各轴进行运动，可得到各个时态刀具的位置、速度、方向、加速度信息。在不考虑机械力的相互作用的前提下，机床、刀具和工件均可视为刚形体，不产生变形，这样仅需知道刀位点的位置、速度、方向就可以得到刀具的姿态，并进行刀具扫描体的计算。

另外，有别于三轴数控加工，五轴数控加工由于其加工自由度增加，传统的应用于三轴的刀具扫描体算法不再适合于五轴数控加工。如果能够找到适当的方法，将刀具在空间连续运动所生成的扫描体求出，再用三维布尔运算将工件减去包络面，就可以得到最终成品形状。

扫描体生成算法

刀具在空间中的运动可以视为扫描操作，根据扫描操作路径形态，可以将一个复杂的扫描路径分解为平移扫描，旋转扫描及一般三维扫描^[17]。其中平移扫描和旋转扫描在工业造型应用较多。

平移扫描

记扫描母体为对象A，沿着一条自由扫描路径P做平移操作，生成扫描体的过程称之为平移扫描。其数学上可计算A与B的Minkowski和来得到。计算Minkowski和的实质就是沿着多边形的边对多面体进行扫描操作^[18]。

旋转扫描

旋绕扫描定义为以多面体为扫描母体单轴或者放射状的绕着某点进行旋转，以此进行扫描体构造。Korein提出求扫描体母体的极限特征，然后用分段式平面逼近求得扫描体的外轮廓曲面^[19]。

一般三维扫描

由于一般三维扫描的扫描母体、扫描路径、以及坐标的标架的不同，算法的复杂程度往往有很大的差别。在三轴数控当中，仅需应用平移扫描和旋转扫描即可完成扫描体的构造。在五轴数控中，由于涉及到AB轴的自由旋转，导致扫描体会随着刀具姿态发生变化，平移扫描和旋转扫描不能完成这类扫描体的构造。对于此类扫描体的构造，大量的研究应用包络面理论作为扫描成型的基础。

包络理论

扫描体的具体定义是一个表面或者一个实体在一定时间内沿某一轨迹运动时所占用的几何空间。参与运动的物体称之为扫描机基体G。G的位置可以由位置函数L(t)确定，则G在时刻t的几何轮廓可表示为。将t标准化为后，将各时段相加，可表示为：

此公式对体素表示的实体可进行离散计算。对于B-rep表达的扫描基体G，用表示实体扫描体G的边界，表示第k个面，和分别表示在扫描起始位置和最终位置的轮廓，则可分解为另一种表达：

通过此公式，可得到扫描体边界扫描构造的扫描体。

在五轴数控加工中，刀具旋转体的形状可用参数化曲面表示，而参数化曲面可转换为隐式表达。将刀具旋转体离散为数个单独的面f，其具有如下内在形式：

当该面按照路径运动时，将生成一簇连续运动面F。对于某一时刻t，其F可用下式表达：

设扫描基体移动，足够小，并且在和之间存在交线，则上式可以用下式代替：

考虑，方程（4- 5）可写为如下的微分形式：

满足下面两个方程的点C的集合称为面f与族F有关的特征点。

当时间在一定时间域内变化时，如果族F的特征点存在，且构成一个面，则该面称之为族F的包络面。消掉时间参数t，可得到包络面E的隐式方程如下：

从几何学上讲，包络线包括了族F沿相应特征点的每个运动面。这意味着运动面F上每一点的切面都与包络面E重合。设在时刻t，包络面E和面F上有一点P，其瞬时速度为，则点P必须满足下面条件：

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

所以，整个曲面族F 的包络面可定义为：

该式几何意义明显，可用于扫描实体包络的求解。当以实体作为扫描母体时，相应的扫描实体包络里包含了母体的外轮廓面以及各条边扫描生成的扫描面。但当遇到包络面退化或者自相交情况时，需要进一步进行检查和处理，算法较为复杂，要得到较高的鲁棒性比较困难[11]。

基于CPU的扫描算法

基于边界的扫描

基于B-rep描述的边界的扫描体的构造依赖于扫描母体的形状及扫描路径，此类造型的优势在于数学理论充分，几何意义明显，精度较高，但算法复杂，对于一些退化情况和自交问题较难处理。在边界扫描的实践中，Minkowski和于包络理论适用于一般扫描体的边界计算。由于算法的复杂性，此类算法的并行度不高，很难应用于多核和GPU的编程环境。

基于体积的扫描

基于体积的扫描操作将扫描母体视为有限的离散体近似构成，类似于现实物体的构成方式。将扫描母体的所有离散体经过刚性变换后得到的位置进行合并，可得到最终扫描体的近似解。

由于此类扫描方法基于离散空间的方式，适应面更广，对于路径和扫描体的形状没有任何特定的要求。简单的提高离散的精度可相应提高扫描的精度，可根据应用场所，应用光线投射（Ray-casting）算法或Marching Cubes算法进行显示。相对于基于边界扫描，算法避免了自交处理的繁琐算法，算法简单易行。由于其针对单个粒子进行计算，不涉及粒子之间的计算，所以算法的并行度高，可应用于多核编程环境。但通常计算量很庞大，精度与计算量成立方比，当今CPU的运算能力和内存带宽不足以支持更高精度的离散。