使用 OpenGL 和 OpenCL 进行拾取选择

引言

在任何3D应用程序中，无论是视频游戏还是CAD工具，最基本的操作之一就是通过单击来选择一个对象。对用户来说，这是一个不费脑子的事情，他们期望拾取选择能够快速准确地完成。但对于开发者来说，编写这个操作的难度一点也不小；它需要对计算几何和高性能编程有扎实的理解。具体来说，拾取选择涉及到找到与给定射线相交的第一个多边形（三角形）。本文的目的是提出一种使用OpenCL（开放计算语言）实现高速拾取选择的方法。与大多数在CPU上执行的拾取选择例程不同，本文的代码旨在在GPU上执行。GPU在通用计算方面不如CPU，但它们能通过并行处理快速计算数学运算。本文的目的是展示如何利用GPU的并行处理来检测射线-三角形相交。

在撰写本文时，我假设您以前从未编写过拾取选择例程，并且对拾取选择与射线-三角形相交的关系一无所知。因此，本文的第一部分描述了拾取选择的底层理论。第二部分解释了如何在C和OpenCL中实现拾取选择。

图1展示了OpenGL中的示例渲染，它包含十个选择目标。右侧的图形说明了结果：橙色球体被染成了白色，因为用户点击了它。

图1：拾取选择目标，选择前和选择后

本文的目的是解释这个拾取选择例程是如何在C和OpenCL中实现的。但在您深入研究代码之前，您应该对底层理论有扎实的理解。

1. 拾取选择背后的理论

尽管对用户来说很简单，但拾取选择对程序员来说是一个复杂的问题。实现这个过程需要对坐标系、矩阵-向量运算和重心坐标有透彻的了解。本节描述的方法与Tomas Moller和Ben Trumbore的论文《Fast, Minimum Storage Ray/Triangle Intersection》中所描述的方法相同。您可以在此处下载该论文。

1.1 将鼠标点击转换为射线

当用户在渲染窗口中点击时，事件处理例程会收到鼠标指针的(x, y)位置。这个位置以窗口坐标给出，其中x范围从0到窗口的像素宽度，y范围从0到窗口的像素高度。要理解拾取选择是如何工作的，您需要了解这些坐标与其他3D渲染中使用的坐标系的关系。

1. 对象坐标 - 最初定义OpenGL图形的3D顶点以对象坐标给出。模型视图变换将模型相对于观察者定位，并将对象坐标转换为眼睛坐标。
2. 眼睛坐标 - 在这个坐标系中，观察者位于原点，并沿着-z轴看向。投影变换定义了观察区域，并将眼睛坐标转换为裁剪坐标。
3. 裁剪坐标 - 在这个坐标系中，x、y和z坐标的范围是从-w到w，其中w是裁剪坐标。透视除法将裁剪坐标除以w，并将裁剪坐标转换为归一化设备坐标。
4. 归一化设备坐标(NDCs) - 在这个系统中，x、y和z坐标的范围是从-1到1。视口变换将NDCs缩放到窗口大小，并将NDCs转换为窗口坐标。
5. 窗口坐标 - 这些坐标在x方向上从0到窗口的像素宽度，在y方向上从0到窗口的像素高度。

理想情况下，我们会知道模型中每个图形的窗口坐标或NDCs。但这些坐标在渲染过程中计算得出，应用程序无法访问。但我们有图形的对象坐标。这就引出了我们的第一个问题：在我们确定用户选择了哪个图形之前，我们需要确保鼠标点击（窗口坐标）和图形的3D顶点（对象坐标）在同一个坐标系中。

有两种选择。我们可以将点击的窗口坐标转换为对象坐标，或者将模型中每个图形的每个顶点转换为窗口坐标。第一种选择需要的计算量少得多，所以我们将使用这种方法。

将窗口坐标转换为NDCs是直接的。NDCs在每个维度上的范围都是-1到1，所以我们可以通过从两个坐标中减去最大维数的一半，然后将结果除以最大维数的一半，来将点击的位置转换为NDC系统。图2展示了这一点是如何工作的。左侧，窗口的尺寸是400x300，用户点击了点(320, 195)。

图2：窗口坐标和归一化设备坐标 (NDCs)

注意：图中窗口坐标在左上角为(0, 0)。这是GLUT约定，对于其他工具集可能不适用。

在图的右侧，用户的二维点(320, 195)被描绘成一个三维箭头，从初始点(0.6, -0.3, -1.0)延伸出来。从几何学上讲，这个箭头被称为射线，其初始点被称为原点，记为O。因为O是以归一化设备坐标给出的，我们将称之为O_ndc，等于(0.6, -0.3, -1.0)。同样，射线的方向记为D_ndc，因为射线指向+z方向，D_ndc等于(0, 0, 1)。构成射线的点集可以用数学表示为O_ndc + tD_ndc，其中t表示从O到射线上某点的距离。请注意，t总是大于零。

1.2 将射线转换为对象坐标

我们在归一化设备坐标中定义了一个射线，但这还不够。为了看到射线与模型中的对象相交的位置，我们需要将O_ndc和D_ndc计算到对象坐标系中。为了理解这一点是如何工作的，了解模型视图和投影变换是如何执行的是很重要的。假设一个顶点v的对象坐标是(v_x, v_y, v_z)。模型视图变换是通过将v乘以一个称为模型视图矩阵的矩阵来执行的。如果我们称这个矩阵为M，模型视图变换的结果等于M * v。

同样，一个点的投影变换是通过将该点乘以透视矩阵P来执行的。因此，要将v从对象坐标转换，应用程序会计算P * (M * v)。矩阵乘法是结合的，所以这等于(P * M) * v。P和M的乘积也是一个矩阵，我们将称之为模型视图投影矩阵，或MVP矩阵。

我们的目标是将O_ndc和D_ndc转换为对象坐标系，而不是从对象坐标系转换。因此，我们需要反转MVP矩阵表示的变换。实际上，我们可以通过将O_ndc和D_ndc乘以MVP矩阵的逆，记为MVP^-1，来实现这一点。用方程表示如下：

O_obj = MVP^-1 * O_ndc
D_obj = MVP^-1 * D_ndc

计算矩阵的逆是复杂的。幸运的是，OpenGL数学库提供了求逆函数。以下代码展示了它的工作原理：

mvp_inverse = glm::inverse(mvp_matrix); 

此时，我们已经推导出了从用户鼠标点击导出射线并将其转换为对象坐标（O_obj + tD_obj）的方法。接下来，我们将看看如何计算该射线与模型中的三角形的相交位置。为了简化起见，我们将删除下标，将射线标识为O + tD。

1.3 重心坐标

OpenGL通过绘制三角形来渲染三维表面，图1中的渲染包含数千个三角形。当射线从其原点向外传播时，它可能与零个、一个或多个三角形相交。为了拾取选择的目的，我们想要距离射线原点最近的三角形，因为它必须属于离用户最近的图形。用伪代码表示，算法如下：

t_min = 1000
for each figure in the model
   for each triangle in the figure
      if the ray (O + tD) intersects the triangle
         if t is less than t_min
            set t_min equal to t
            store the figure containing the triangle
         end if
      end if
   end for
end for
if t_min = 1000
   no intersection
else
   the selected figure = figure corresponding to t_min 

目前，我们关心的是确定“射线是否与三角形相交”。要了解射线-三角形相交是如何计算的，理解重心坐标的概念很重要。这些坐标唯一地标识了形状上的点，并沿着重心轴进行测量。每个三角形有三个重心轴，每个轴通过一个顶点和对面边的中点。图3描绘了一个通过顶点K的重心轴。

图3：重心轴

在这个图中，虚线是顶点K的重心轴。相应的重心坐标记为k。在顶点K处，k等于1，在对面边的中点处，k等于0。在对面边的外侧，k小于零，在K的右侧，k大于1。

图4展示了三角形KLM的所有三个重心轴。给出的相应重心坐标（(k, l, m)三元组）标识了三角形内部及其边上的点。

图4：重心坐标

在处理重心坐标时，有两点需要记住：

1. 如果一个点的所有坐标都介于0和1之间，那么这个点就位于形状内部。
2. 无论点位于何处，其重心坐标之和始终为1。因此，任何一个坐标都可以被其他两个坐标之和减去1来代替。也就是说，(k, l, m)始终等于(k, l, 1 - k - l)。

出现了一个重要问题。如何将重心坐标(k, l, m)转换为矩形坐标(x, y, z)？答案很简单：将每个重心坐标乘以其对应顶点的坐标。例如，如果一个点的重心坐标是(k, l, m)，那么它的矩形坐标可以通过以下方程计算：

P = kK + lL + mM

举个例子就能清楚这一点。假设三角形KLM的顶点是K = (13, 3)，L = (9, 10)，M = (1, 1)。下图计算了点P0、P1、P2和P3的矩形坐标。

图5：将重心坐标转换为矩形坐标

现在我们知道了如何在数学上表示三角形上的点，我们可以研究如何确定一个射线是否与三角形上的点相交以及在哪里相交。这是下一节讨论的主题。

1.4 计算射线-三角形相交

如果我们用O + tD表示射线，用kK + lL + mM表示三角形的点，我们可以用以下方程计算它们的交集：

O + tD = kK + lL + mM

重心坐标(k, l, m)可以表示为(k, l, 1 - k - l)，所以我们可以将这个方程改写如下：

O + tD = kK + lL + (1 - k - l)M

这给出了三个未知数：t，k和l。通过重新排列项，我们得到以下线性方程组：

为了简化这个方程，我们将用E替换K - M，用F替换L - M，用G替换O - M。这就得到了以下结果：

克莱姆法则可以很容易地求解t、k和l。得到的方程如下：

在这个方程中，竖线表示行列式。也就是说，|A|是A的行列式。行列式可以通过公式|A B C| = -(A x C)·B = -(C x B)·A来计算。重新排列项得到以下关系：

通过测试k和l，我们可以找出点是否位于K、L和M构成的三角形上。然后，如果t是所有测试三角形中最小的正值，我们就可以确定该三角形属于用户选择的图形。

2. 在代码中实现射线-三角形相交

现在我们已经研究了拾取选择背后的数学原理，了解如何将该算法在代码中实现非常重要。本节首先介绍执行射线-三角形相交的常规C代码。然后，我将展示如何使用OpenCL加速这一点。但在我展示代码之前，我需要介绍ColGeom数据结构，它存储了模型的顶点数据。

2.1 ColGeom数据结构

本文讨论了COLLADA和ColGeom数据结构，它存储了3D模型中图形的数据。这个数据结构定义如下：

struct ColGeom {
   std::string name;         // The ID attribute of the <geometry>element
   SourceMap map;            // Contains data in the <source />elements
   GLenum primitive;         // Identifies the primitive type, such as GL_LINES
   int index_count;          // The number of indices used to draw elements
   unsigned int* indices;    // The index data from the <p> element
};</geometry>

map字段，类型为SourceMap，包含结构所代表图形的几何数据。此类型定义如下：

typedef std::map<std::string,> SourceMap;</std::string,> 

map将源的语义名称与其数据匹配。此名称可以取POSITION、NORMALS或TEXCOORDS等值。SourceData元素包含对象网格的实际数据，定义如下：

struct SourceData {
   GLenum type;              // The data type of the mesh data, such as GL_FLOAT
   unsigned int size;        // Size of the mesh data in bytes
   unsigned int stride;      // Number of data values per group
   void* data;               // Mesh data
}; 

在本文中，data数组包含表示顶点坐标的浮点值。根据stride字段的规定，这些顶点以三个一组（x、y和z）访问。indices字段告诉我们顶点是如何组合成三角形的。例如，如果前六个索引是{0, 1, 2, 1, 0, 3}，这意味着第一个三角形由第一个、第二个和第三个顶点组成，第二个三角形由第二个、第一个和第四个顶点组成。顺序很重要——顶点的方向（顺时针或逆时针）决定了哪些三角形会被渲染器剔除，哪些不会。

2.2 在C中实现拾取选择

现在我们知道了如何访问图形数据，我们可以创建一个例程来实现通过射线-三角形相交进行拾取选择。以下代码展示了如何在C中实现这一点。请注意，鼠标坐标以x和y给出，网格的索引数组通过geom_vec[i].indices访问，顶点数据通过geom_vec[i].map["POSITION"].data访问。

// Compute O and D in object coordinates
glm::vec4 origin = mvp_inverse * glm::vec4(
        (x-half_width)/half_width, (half_height-y)/half_height, -1, 1);
glm::vec4 dir = mvp_inverse * glm::vec4(0, 0, 1, 0);
O = glm::vec3(origin.x, origin.y, origin.z);
D = glm::normalize(glm::vec3(dir.x, dir.y, dir.z));
// Iterate through the figures in the model
tmin = 1000.0;
for(i=0; i<num_objects; i++) {
   data_array = (float*)geom_vec[i].map["POSITION"].data;
   // Iterate through the triangles in the figure
   for(j=0; j<geom_vec[i].index_count; j+=3) {
      index = geom_vec[i].indices[j]*3;
      // Read the first point of Triangle KLM
      K = glm::vec3(data_array[index],
                    data_array[index+1],
                    data_array[index+2]);
      // Read the second point of Triangle KLM
      index = geom_vec[i].indices[j+1]*3;
      L = glm::vec3(data_array[index],
                    data_array[index+1],
                    data_array[index+2]);
      // Read the third point of Triangle KLM
      index = geom_vec[i].indices[j+2]*3;
      M = glm::vec3(data_array[index],
                    data_array[index+1],
                    data_array[index+2]);
      // Compute vectors E, F, and G
      E = K - M;
      F = L - M;
      G = O - M;
      // Compute the vector containing t and the coordinates k and l
      tkl = 1/glm::dot(glm::cross(D, F), E) *
               glm::vec3(glm::dot(glm::cross(G, E), F),
                         glm::dot(glm::cross(D, F), G),
                         glm::dot(glm::cross(G, E), D));
      // Check if t and the intersection point (k, l) are acceptable
      if(tkl.x < tmin[i] && tkl.y > 0.0f && tkl.z > 0.0f && (tkl.y + tkl.z) < 1) {
         tmin = tkl.x;
         sphere = i;
      }
   }
} 

这段代码包含两个循环：一个循环遍历模型中的每个图形，另一个循环遍历图形中的每个三角形。每次处理一个三角形时，计算出的距离tkl.x都会与最小距离tmin进行比较。如果tkl.x小于tmin并且交点位于三角形内部，则输出值sphere将被设置为对象索引i。一旦模型中的所有三角形都处理完毕，就必须测试tmin。如果tmin小于1000，则sphere的值将标识用户选择的对象。如果tmin等于1000，则意味着射线没有与任何三角形相交，用户的鼠标点击没有落在任何图形上。

2.3 在OpenCL中实现拾取选择

前面展示的C代码能够完成任务，但如果模型中有许多图形，遍历所有三角形可能会花费大量时间。这时OpenCL就能派上用场了。一个称为内核的OpenCL例程可以同时被数千个线程执行。因此，如果我们一次性将一个图形的所有三角形发送到GPU，工作项将并行处理它们，从而实现比常规C应用程序更好的性能。

本文附带的应用程序使用OpenGL渲染图1中所示的模型，并使用OpenCL执行拾取选择。整个应用程序工作流程如下：

1. 应用程序启动时，它会读取一个COLLADA文件（spheres.dae）并为十个球体中的每一个初始化一个ColGeom结构。对于每个球体，应用程序还创建了一个顶点数组对象（VAO）、一个索引缓冲区对象（IBO）和两个顶点缓冲区对象（VBO）。第一个VBO存储球体顶点的坐标，第二个存储每个顶点处的法线向量分量。
2. 当用户在窗口内单击时，GLUT回调函数（mouse）将鼠标事件的窗口坐标转换为对象坐标，并调用名为execute_selection_kernel的函数。
3. execute_selection_kernel函数从O和D创建内核参数，然后遍历模型中的每个图形。在每次迭代中，函数创建两个OpenCL缓冲区对象：一个来自包含图形顶点位置的VBO，另一个来自包含图形索引的IBO。在将这些缓冲区对象设置为内核参数后，函数执行内核。
4. pick_sphere内核实现了射线-三角形相交。更准确地说，每个工作项使用O、D、K、L和M为给定三角形计算E、F和G。然后计算并测试重心坐标k和l。如果这些坐标具有可接受的值，工作项将计算t，它确定了从射线原点到三角形的距离。每个工作组中的工作项完成其计算后，工作组确定最小的t值并将其返回给应用程序。
5. 应用程序找到最小的t值并确定相应的球体。它将该球体的颜色设置为白色，从而指示已选中。

以下代码展示了pick_selection内核中的OpenCL指令。每个工作项读取三角形的三个坐标，其目标是计算一个值t_loc，它标识了从射线原点到三角形的距离。如果射线不与三角形相交，则从原点到三角形的距离保持在10000.0。

t_loc[get_local_id(0)] = 10000.0f;
if(get_global_id(0) < 528) {
   /* Read coordinates of triangle vertices */
   indices = vload3(get_global_id(0), ibo);
   K = vload3(indices.x, vbo);
   L = vload3(indices.y, vbo);
   M = vload3(indices.z, vbo);
   /* Compute vectors */
   E = K - M;
   F = L - M;
   /* Compute and test determinant */
   t_test = dot(cross(D.s012, F), E);
   if(t_test > 0.0001f) {
      /* Compute and test k */
      G = O.s012 - M;
      k = dot(cross(D.s012, F), G);
      if(k > 0.0f && k <= t_test) {
         /* Compute and test l */
         l = dot(cross(G, E), D.s012);
         if(l > 0.0f && k + l <= t_test) {
            /* Compute distance from ray to triangle */
            t_loc[get_local_id(0)] = dot(cross(G, E), F)/t_test;
         }
      }
   }
}
barrier(CLK_LOCAL_MEM_FENCE); 

如果初始行列式(D x F)·E为负，则三角形朝向模型的背面而不是正面。因此，它不是测试的候选。工作项通过测试行列式（t_test）是否为负来检查此条件。

但是，如果行列式为正且重心坐标（k和l）有效，则工作项将t_loc设置为((G x E)·F)/行列式。之后，每个工作组中的第一个工作项查找最小的t_loc值，并将其传递给主机应用程序。

3. 结论

拾取选择使用户能够通过2D鼠标点击与3D模型进行交互。这个过程通过射线-三角形相交在数学上进行了建模。本文介绍了该方法的理论，并展示了如何使用常规C和OpenCL实现该理论。

在GPU上执行OpenCL内核的主要优点是，内核可以同时由数百或数千个工作项执行。本文中的拾取选择代码利用了这一点，当内核在GPU上执行时，工作项会并行处理图形的三角形。然而，通过将每个图形的顶点发送到GPU的单个内核中，可以进一步提高性能。

4. 使用代码

本文的代码存档包含运行应用程序所需的源文件。它还包含球体的COLLADA文件（sphere.dae）和项目的Makefile。

历史

提交编辑审批：2013/7/24