Marching Cubes 算法三探

Marching Cubes 算法再探

CUDA Samples
MarchingCubes
- workflow
- Code
- - data structure
  - computeIsosurface
  - launch_classifyVoxel
  - - classifyVoxel
  - ThrustScanWrapper
  - launch_compactVoxels
  - - compactVoxels
  - launch_generateTriangles2
  - - generateTriangles2
- improvements

初探再探三探，现在来到了第三关：CUDA Samples中的Marching Cubes

如果对MC有兴趣，建议先走第一关和第二关（此关不涉及MC原理）

CUDA Samples

CUDA Sampels在安装CUDA时选装，默认安装位置为C:\ProgramData\NVIDIA Corporation\CUDA Samples\下，且在安装时建立了快捷方式Browse CUDA Samples，可以直接在开始菜单或通过工具搜索到
在这里插入图片描述

还贴心准备了各个vs版本的解决方案，选择对应版本打开
在这里插入图片描述

注意，如果在安装时，将CUDA samples 装在了C盘，那么需要以管理员启动VS，否则会有访问权限报错，如

在这里插入图片描述

不知道为是不是个例，第一次以管理员身份启动时，比不以管理员身份启动慢非常多

在这里插入图片描述

MarchingCubes

在解决方案目录树中搜索MarchingCubes可以找到其官方示例，右键设为启动项目，就可以开始使用了

在这里插入图片描述

先跑一下看有无报错

在这里插入图片描述
不出意外的话都能成功运行如下

可以通过鼠标右键交互获得不同效果如
在这里插入图片描述

workflow

marchingCubes.cpp中的注释有说明算法核心步骤

/*
    行进立方体算法
    该示例使用行进立方体算法从体积数据集中提取几何等值面。它利用了 Thrust 库中的扫描（前缀求和）函数来执行流压缩。类似的技术可用于需要每个线程产生可变大小输出的问题。

    算法包含几个阶段：

    1. 执行 "classifyVoxel" 内核
    在每个体素的角上评估体数据，并计算每个体素将产生的顶点数量。
    每个体素使用一个线程执行。
    它写入两个数组 - voxelOccupied 和 voxelVertices 到全局内存。
    voxelOccupied 是一个标志，指示体素是否非空。

    2. 扫描 "voxelOccupied" 数组（使用 Thrust 扫描）
    从 GPU 读回非空体素的总数到 CPU。
    这是排他性扫描最后一个值和最后一个输入值之和。

    3. 执行 "compactVoxels" 内核
    压缩 voxelOccupied 数组以去除空体素。
    这允许我们仅对占用的体素运行复杂的 "generateTriangles" 内核。

    4. 扫描 voxelVertices 数组
    这给出了每个体素顶点数据的起始地址。
    我们从 GPU 读回生成的总顶点数到 CPU。

    注意：通过使用自定义扫描函数，我们可以将上述两个扫描操作合并为一个操作。

    5. 执行 "generateTriangles" 内核
    仅在占用的体素上运行。
    再次查找场值并生成三角形数据，
    使用扫描的结果将输出写入正确的地址。
    行进立方体查找表存储在 1D 纹理中。

    // 6. 渲染几何体
    // 利用得到的顶点进行渲染
*/

通过断点调试，也可以得到算法默认执行路线（省略glut部分）如下

在这里插入图片描述

可以看到核心是computeIsosurface函数，5个步骤与上面算法5个阶段对应，接下来看看核心代码

Code

结合GPT注释来理解代码

众所周知，打Boss仅靠个人拳脚是很优雅取胜的，这时如果有趁手的技能、兵器在身或得力帮手在侧，那就从容许多。“广智救我！”是本吗喽玩黑悟时的呼声。“GPT救我！”是本菜狗看代码时的心声。

data structure

先简单看一下数据结构

//网格相关
uint3 gridSizeLog2 = make_uint3(5, 5, 5);		// 网格大小的对数值 5 表示 2^5
uint3 gridSizeShift;							// 用于快速计算网格索引，通过移位操作来代替除法
uint3 gridSize;									// 存储的是实际的网格大小，即 gridSize.x = 2^5 等
uint3 gridSizeMask;								// 用于掩码操作，通常用来获取网格索引中的特定位，大小为gridSize-1

float3 voxelSize;								//每个体素的物理大小，voxelSize 用于将体素索引转换为实际空间坐标
uint numVoxels    = 0;		// 体素的总数量				
uint maxVerts     = 0;		// 可以生成的最大顶点数量
uint activeVoxels = 0;		// 活跃体素的数量，即那些参与等值面提取的体素
uint totalVerts   = 0;		// 生成的顶点总数

float isoValue      = 0.2f;		// 等值面值与下面的dIsoValue结合动态变化以形成动画
float dIsoValue     = 0.005f;

// 设备数据
GLuint posVbo, normalVbo;		// OpenGL 顶点缓冲对象，用于存储顶点位置和法线
GLint  gl_Shader;				// OpenGL 着色器程序的 ID
struct cudaGraphicsResource *cuda_posvbo_resource, *cuda_normalvbo_resource; // handles OpenGL-CUDA exchange 用于 CUDA 和 OpenGL 之间交换数据的资源句柄

float4 *d_pos = 0, *d_normal = 0;	// 指向设备端存储顶点位置和法线的指针

uchar *d_volume = 0;				// 体数据的指针		
uint *d_voxelVerts = 0;				// 每个体素生成的顶点数
uint *d_voxelVertsScan = 0;			// 用于存储 d_voxelVerts 的前缀和
uint *d_voxelOccupied = 0;			// 标识哪些体素是活跃的
uint *d_voxelOccupiedScan = 0;		// d_voxelOccupied 的前缀和。
uint *d_compVoxelArray;				// 压缩后的体素数组

// 查找表
uint *d_numVertsTable = 0;			// 存储每个体素的三角形顶点数量
uint *d_edgeTable = 0;				// 存储边对应的顶点索引
uint *d_triTable = 0;				// 存储三角形的顶点索引

computeIsosurface

#define DEBUG_BUFFERS 0  // debug 开关，用以输出一些信息
///
//! 运行计算的 CUDA 部分
///
void
computeIsosurface()
{
    int threads = 128;						// 每个线程块（block）中有 128 个线程
    dim3 grid(numVoxels / threads, 1, 1);  	//一维网格，只设置了 x 轴上的线程块数目(y,z维度为1)

	// CUDA 设备对每个维度上的线程块数目有限制，最大通常不超过 65535  
    if (grid.x > 65535)
    {
        grid.y = grid.x / 32768;			// 一维不够那就升维
        grid.x = 32768;
    }

    // 计算每个体素需要的顶点数量
    launch_classifyVoxel(grid, threads,
                         d_voxelVerts, d_voxelOccupied, d_volume,
                         gridSize, gridSizeShift, gridSizeMask,
                         numVoxels, voxelSize, isoValue);
#if DEBUG_BUFFERS
    printf("voxelVerts:\n");
    dumpBuffer(d_voxelVerts, numVoxels, sizeof(uint));
#endif

#if SKIP_EMPTY_VOXELS
    // 扫描体素占用数组
    ThrustScanWrapper(d_voxelOccupiedScan, d_voxelOccupied, numVoxels);

#if DEBUG_BUFFERS
    printf("voxelOccupiedScan:\n");
    dumpBuffer(d_voxelOccupiedScan, numVoxels, sizeof(uint));
#endif

    // 读取回最后一个元素以计算非空体素的总数
    // 由于使用的是前缀扫描（exclusive scan），总数为扫描结果的最后一个值加上输入数组的最后一个值（因为输入数组的最后一个值没有倍计算到前缀和中）
    {
        uint lastElement, lastScanElement;
        checkCudaErrors(cudaMemcpy((void *) &lastElement,
                                   (void *)(d_voxelOccupied + numVoxels-1),
                                   sizeof(uint), cudaMemcpyDeviceToHost));
        checkCudaErrors(cudaMemcpy((void *) &lastScanElement,
                                   (void *)(d_voxelOccupiedScan + numVoxels-1),
                                   sizeof(uint), cudaMemcpyDeviceToHost));
        activeVoxels = lastElement + lastScanElement;	// 非空体素的总数
    }

    if (activeVoxels == 0)
    {
        // 如果没有有效体素，则返回
        totalVerts = 0;
        return;
    }

    // 压缩体素索引数组
    launch_compactVoxels(grid, threads, d_compVoxelArray, d_voxelOccupied, d_voxelOccupiedScan, numVoxels);
    getLastCudaError("compactVoxels failed");

#endif // SKIP_EMPTY_VOXELS

    // 扫描体素顶点计数数组
    ThrustScanWrapper(d_voxelVertsScan, d_voxelVerts, numVoxels);

#if DEBUG_BUFFERS
    printf("voxelVertsScan:\n");
    dumpBuffer(d_voxelVertsScan, numVoxels, sizeof(uint));
#endif

    // 读取返回的顶点总数 
    {
        uint lastElement, lastScanElement;
        checkCudaErrors(cudaMemcpy((void *) &lastElement,
                                   (void *)(d_voxelVerts + numVoxels-1),
                                   sizeof(uint), cudaMemcpyDeviceToHost));
        checkCudaErrors(cudaMemcpy((void *) &lastScanElement,
                                   (void *)(d_voxelVertsScan + numVoxels-1),
                                   sizeof(uint), cudaMemcpyDeviceToHost));
		// 所有体素中所需生成的顶点的总数
        totalVerts = lastElement + lastScanElement;
    }

    // 生成三角形，写入顶点缓冲区
    if (!g_bValidate)
    {
        size_t num_bytes;
        // 将 OpenGL 缓冲区对象映射到 CUDA 地址空间
        // DEPRECATED 旧方法（已弃用）: checkCudaErrors(cudaGLMapBufferObject((void**)&d_pos, posVbo));
        // 使用 CUDA 图形互操作 API 将 OpenGL 顶点缓冲区映射到 CUDA 内存空间
        checkCudaErrors(cudaGraphicsMapResources(1, &cuda_posvbo_resource, 0));
        checkCudaErrors(cudaGraphicsResourceGetMappedPointer((void **)&d_pos, &num_bytes, cuda_posvbo_resource));

		// 将法线缓冲区映射到 CUDA 地址空间
        // DEPRECATED: checkCudaErrors(cudaGLMapBufferObject((void**)&d_normal, normalVbo));
        checkCudaErrors(cudaGraphicsMapResources(1, &cuda_normalvbo_resource, 0));
        checkCudaErrors(cudaGraphicsResourceGetMappedPointer((void **)&d_normal, &num_bytes, cuda_normalvbo_resource));
    }

#if SKIP_EMPTY_VOXELS
	// 计算网格维度，跳过空的体素
    dim3 grid2((int) ceil(activeVoxels / (float) NTHREADS), 1, 1);
#else
    dim3 grid2((int) ceil(numVoxels / (float) NTHREADS), 1, 1);
#endif

    // 解决 grid2.x 超过 65535 的问题
    while (grid2.x > 65535)
    {
        grid2.x /= 2;
        grid2.y *= 2;
    }

#if SAMPLE_VOLUME
    launch_generateTriangles2(grid2, NTHREADS, d_pos, d_normal,
                              d_compVoxelArray,
                              d_voxelVertsScan, d_volume,
                              gridSize, gridSizeShift, gridSizeMask,
                              voxelSize, isoValue, activeVoxels,
                              maxVerts);
#else
    launch_generateTriangles(grid2, NTHREADS, d_pos, d_normal,
                             d_compVoxelArray,
                             d_voxelVertsScan,
                             gridSize, gridSizeShift, gridSizeMask,
                             voxelSize, isoValue, activeVoxels,
                             maxVerts);
#endif

    if (!g_bValidate)
    {
        // DEPRECATED:      checkCudaErrors(cudaGLUnmapBufferObject(normalVbo));
        checkCudaErrors(cudaGraphicsUnmapResources(1, &cuda_normalvbo_resource, 0));
        // DEPRECATED:      checkCudaErrors(cudaGLUnmapBufferObject(posVbo));
        checkCudaErrors(cudaGraphicsUnmapResources(1, &cuda_posvbo_resource, 0));
    }
}

关于65535 为 $2^16-1$ ， $2^16$ 即64K。

launch_classifyVoxel

写入两个数组 - voxelOccupied 和 voxelVertices 到全局内存。
    voxelOccupied 是一个标志，指示体素是否非空

extern "C" void
launch_classifyVoxel(dim3 grid, dim3 threads, uint *voxelVerts, uint *voxelOccupied, uchar *volume,
                     uint3 gridSize, uint3 gridSizeShift, uint3 gridSizeMask, uint numVoxels,
                     float3 voxelSize, float isoValue)
{
    // 每个体素计算所需的顶点数
    classifyVoxel<<<grid, threads>>>(voxelVerts, voxelOccupied, volume,
                                     gridSize, gridSizeShift, gridSizeMask,
                                     numVoxels, voxelSize, isoValue);
    // 检查 classifyVoxel CUDA 核函数是否正确执行
    getLastCudaError("classifyVoxel failed");
}

查找表在 ‪C:\ProgramData\NVIDIA Corporation\CUDA Samples\v10.2\2_Graphics\marchingCubes\tables.h中
这里涉及到一个纹理绑定操作

classifyVoxel

// 根据每个体素将生成的顶点数量对体素进行分类
// 每个线程处理一个体素
__global__ void
classifyVoxel(uint *voxelVerts, uint *voxelOccupied, uchar *volume,
              uint3 gridSize, uint3 gridSizeShift, uint3 gridSizeMask, uint numVoxels,
              float3 voxelSize, float isoValue)
{
    // 计算当前线程处理的体素索引
    uint blockId = __mul24(blockIdx.y, gridDim.x) + blockIdx.x;
    uint i = __mul24(blockId, blockDim.x) + threadIdx.x;

    // 根据体素索引计算体素在网格中的三维位置
    uint3 gridPos = calcGridPos(i, gridSizeShift, gridSizeMask);

    // 读取相邻网格顶点处的场值
    #if SAMPLE_VOLUME
        float field[8];
        field[0] = sampleVolume(volume, gridPos, gridSize);
        field[1] = sampleVolume(volume, gridPos + make_uint3(1, 0, 0), gridSize);
        field[2] = sampleVolume(volume, gridPos + make_uint3(1, 1, 0), gridSize);
        field[3] = sampleVolume(volume, gridPos + make_uint3(0, 1, 0), gridSize);
        field[4] = sampleVolume(volume, gridPos + make_uint3(0, 0, 1), gridSize);
        field[5] = sampleVolume(volume, gridPos + make_uint3(1, 0, 1), gridSize);
        field[6] = sampleVolume(volume, gridPos + make_uint3(1, 1, 1), gridSize);
        field[7] = sampleVolume(volume, gridPos + make_uint3(0, 1, 1), gridSize);
    #else
        // 如果没有 SAMPLE_VOLUME 选项，则使用自定义场函数计算场值
        float3 p;
        p.x = -1.0f + (gridPos.x * voxelSize.x);
        p.y = -1.0f + (gridPos.y * voxelSize.y);
        p.z = -1.0f + (gridPos.z * voxelSize.z);

        float field[8];
        field[0] = fieldFunc(p);
        field[1] = fieldFunc(p + make_float3(voxelSize.x, 0, 0));
        field[2] = fieldFunc(p + make_float3(voxelSize.x, voxelSize.y, 0));
        field[3] = fieldFunc(p + make_float3(0, voxelSize.y, 0));
        field[4] = fieldFunc(p + make_float3(0, 0, voxelSize.z));
        field[5] = fieldFunc(p + make_float3(voxelSize.x, 0, voxelSize.z));
        field[6] = fieldFunc(p + make_float3(voxelSize.x, voxelSize.y, voxelSize.z));
        field[7] = fieldFunc(p + make_float3(0, voxelSize.y, voxelSize.z));
    #endif

    // 计算每个顶点是否在等值面内或外的标志
    uint cubeindex;		// 8bit 案例编码
    cubeindex =  uint(field[0] < isoValue);
    cubeindex += uint(field[1] < isoValue) * 2;
    cubeindex += uint(field[2] < isoValue) * 4;
    cubeindex += uint(field[3] < isoValue) * 8;
    cubeindex += uint(field[4] < isoValue) * 16;
    cubeindex += uint(field[5] < isoValue) * 32;
    cubeindex += uint(field[6] < isoValue) * 64;
    cubeindex += uint(field[7] < isoValue) * 128;

    // 从纹理中读取顶点数量，即查表
    // tex1Dfetch 是一个 CUDA 内置函数，用于从绑定到纹理内存的数组中按索引读取数据。
    uint numVerts = tex1Dfetch(numVertsTex, cubeindex);

    // 如果体素索引在有效范围内，则记录顶点数量和占用状态
    if (i < numVoxels)
    {
        voxelVerts[i] = numVerts;          // 记录体素的顶点数量
        voxelOccupied[i] = (numVerts > 0); // 如果顶点数量大于 0，表示该体素被占用
    }
}

ThrustScanWrapper

调用 Thrust 库中的 exclusive_scan 函数，对输入数据进行前缀和操作。exclusive_scan 的结果不包括输入的第一个元素（即第一个元素的结果为0），而之后的每个元素都是之前所有元素的累加和。
用途
1.确定顶点偏移量:
通过计算前缀和，可以得到每个包含等值面的体素之前有多少个顶点。
这对于后续构建顶点数组非常重要，因为它可以帮助确定每个体素生成的顶点应该放置在顶点数组中的哪个位置。
2.内存管理:
通过知道每个体素之前有多少个顶点，可以更有效地管理内存空间，避免不必要的内存复制和移动。

extern "C" void ThrustScanWrapper(unsigned int *output, unsigned int *input, unsigned int numElements)
{
    thrust::exclusive_scan(thrust::device_ptr<unsigned int>(input),
                           thrust::device_ptr<unsigned int>(input + numElements),
                           thrust::device_ptr<unsigned int>(output));
}

launch_compactVoxels

将包含等值面的体素重新排列到一个新的数组中（压缩），从而去除那些未被占用的体素。这样可以节省内存空间并提高后续处理的效率。

extern "C" void
launch_compactVoxels(dim3 grid, dim3 threads, uint *compactedVoxelArray, uint *voxelOccupied, uint *voxelOccupiedScan, uint numVoxels)
{
    compactVoxels<<<grid, threads>>>(compactedVoxelArray, voxelOccupied,
                                     voxelOccupiedScan, numVoxels);
    getLastCudaError("compactVoxels failed");
}

compactVoxels


// 压缩体素数组
__global__ void
compactVoxels(uint *compactedVoxelArray, uint *voxelOccupied, uint *voxelOccupiedScan, uint numVoxels)
{
    // 计算当前线程在整个网格中的唯一ID
    uint blockId = __mul24(blockIdx.y, gridDim.x) + blockIdx.x;
    uint i = __mul24(blockId, blockDim.x) + threadIdx.x;

    // 如果当前体素是被占用的，并且索引在有效范围内
    if (voxelOccupied[i] && (i < numVoxels))
    {
        // 根据扫描的结果，将当前体素的索引存储到压缩后的体素数组中
        compactedVoxelArray[ voxelOccupiedScan[i] ] = i;
    }
}

launch_generateTriangles2

extern "C" void
launch_generateTriangles2(dim3 grid, dim3 threads,
                          float4 *pos, float4 *norm, uint *compactedVoxelArray, uint *numVertsScanned, uchar *volume,
                          uint3 gridSize, uint3 gridSizeShift, uint3 gridSizeMask,
                          float3 voxelSize, float isoValue, uint activeVoxels, uint maxVerts)
{
    generateTriangles2<<<grid, NTHREADS>>>(pos, norm,
                                           compactedVoxelArray,
                                           numVertsScanned, volume,
                                           gridSize, gridSizeShift, gridSizeMask,
                                           voxelSize, isoValue, activeVoxels,
                                           maxVerts);
    getLastCudaError("generateTriangles2 failed");
}

generateTriangles2

注意到重新计算了体素顶点的8bit编码，且在顶点插值时似乎是在12个边上都进行了计算（那么会不会涉及到后续去重的问题？）
(cubeindex * 16) + i：cubeindex * 16 计算出这个特定 cubeindex 对应的条目在 triTable 中的起始位置（因为每个 cubeindex 对应最多 16 个边编号）。i 是当前正在处理的边的偏移量（以 3 为步长递增）。

// 使用行进立方体为每个体素生成三角形
// 这个版本是计算每个三角形的平面法线，用后缀2区分（另一版本是从场函数插值法线）
__global__ void
generateTriangles2(float4 *pos, float4 *norm, uint *compactedVoxelArray, uint *numVertsScanned, uchar *volume,
                   uint3 gridSize, uint3 gridSizeShift, uint3 gridSizeMask,
                   float3 voxelSize, float isoValue, uint activeVoxels, uint maxVerts)
{
    // 计算当前线程处理的体素索引
    uint blockId = __mul24(blockIdx.y, gridDim.x) + blockIdx.x;
    uint i = __mul24(blockId, blockDim.x) + threadIdx.x;

    // 确保索引不超过活跃体素数
    if (i > activeVoxels - 1)
    {
        i = activeVoxels - 1;
    }

    // 获取当前体素的索引，如果启用了 SKIP_EMPTY_VOXELS，使用压缩后的体素数组
#if SKIP_EMPTY_VOXELS
    uint voxel = compactedVoxelArray[i];
#else
    uint voxel = i;
#endif

    // 计算体素在3D网格中的位置（索引）
    uint3 gridPos = calcGridPos(voxel, gridSizeShift, gridSizeMask);

    // 计算体素在3D空间中的实际位置
    float3 p;
    p.x = -1.0f + (gridPos.x * voxelSize.x);
    p.y = -1.0f + (gridPos.y * voxelSize.y);
    p.z = -1.0f + (gridPos.z * voxelSize.z);

    // 计算体素的8个顶点在3D空间中的位置
    float3 v[8];
    v[0] = p;
    v[1] = p + make_float3(voxelSize.x, 0, 0);
    v[2] = p + make_float3(voxelSize.x, voxelSize.y, 0);
    v[3] = p + make_float3(0, voxelSize.y, 0);
    v[4] = p + make_float3(0, 0, voxelSize.z);
    v[5] = p + make_float3(voxelSize.x, 0, voxelSize.z);
    v[6] = p + make_float3(voxelSize.x, voxelSize.y, voxelSize.z);
    v[7] = p + make_float3(0, voxelSize.y, voxelSize.z);

    // 采样体素的标量场值（或计算场值）
#if SAMPLE_VOLUME
    float field[8];
    field[0] = sampleVolume(volume, gridPos, gridSize);
    field[1] = sampleVolume(volume, gridPos + make_uint3(1, 0, 0), gridSize);
    field[2] = sampleVolume(volume, gridPos + make_uint3(1, 1, 0), gridSize);
    field[3] = sampleVolume(volume, gridPos + make_uint3(0, 1, 0), gridSize);
    field[4] = sampleVolume(volume, gridPos + make_uint3(0, 0, 1), gridSize);
    field[5] = sampleVolume(volume, gridPos + make_uint3(1, 0, 1), gridSize);
    field[6] = sampleVolume(volume, gridPos + make_uint3(1, 1, 1), gridSize);
    field[7] = sampleVolume(volume, gridPos + make_uint3(0, 1, 1), gridSize);
#else
    float field[8];
    field[0] = fieldFunc(v[0]);
    field[1] = fieldFunc(v[1]);
    field[2] = fieldFunc(v[2]);
    field[3] = fieldFunc(v[3]);
    field[4] = fieldFunc(v[4]);
    field[5] = fieldFunc(v[5]);
    field[6] = fieldFunc(v[6]);
    field[7] = fieldFunc(v[7]);
#endif

    // 重新计算体素的二进制标志，判断哪个顶点在等值面之上
    uint cubeindex;
    cubeindex =  uint(field[0] < isoValue);
    cubeindex += uint(field[1] < isoValue)*2;
    cubeindex += uint(field[2] < isoValue)*4;
    cubeindex += uint(field[3] < isoValue)*8;
    cubeindex += uint(field[4] < isoValue)*16;
    cubeindex += uint(field[5] < isoValue)*32;
    cubeindex += uint(field[6] < isoValue)*64;
    cubeindex += uint(field[7] < isoValue)*128;

    // 计算等值面与体素相交的顶点
#if USE_SHARED
    // 使用共享内存以避免使用局部变量
    __shared__ float3 vertlist[12*NTHREADS];

    vertlist[threadIdx.x] = vertexInterp(isoValue, v[0], v[1], field[0], field[1]);
    vertlist[NTHREADS+threadIdx.x] = vertexInterp(isoValue, v[1], v[2], field[1], field[2]);
    vertlist[(NTHREADS*2)+threadIdx.x] = vertexInterp(isoValue, v[2], v[3], field[2], field[3]);
    vertlist[(NTHREADS*3)+threadIdx.x] = vertexInterp(isoValue, v[3], v[0], field[3], field[0]);
    vertlist[(NTHREADS*4)+threadIdx.x] = vertexInterp(isoValue, v[4], v[5], field[4], field[5]);
    vertlist[(NTHREADS*5)+threadIdx.x] = vertexInterp(isoValue, v[5], v[6], field[5], field[6]);
    vertlist[(NTHREADS*6)+threadIdx.x] = vertexInterp(isoValue, v[6], v[7], field[6], field[7]);
    vertlist[(NTHREADS*7)+threadIdx.x] = vertexInterp(isoValue, v[7], v[4], field[7], field[4]);
    vertlist[(NTHREADS*8)+threadIdx.x] = vertexInterp(isoValue, v[0], v[4], field[0], field[4]);
    vertlist[(NTHREADS*9)+threadIdx.x] = vertexInterp(isoValue, v[1], v[5], field[1], field[5]);
    vertlist[(NTHREADS*10)+threadIdx.x] = vertexInterp(isoValue, v[2], v[6], field[2], field[6]);
    vertlist[(NTHREADS*11)+threadIdx.x] = vertexInterp(isoValue, v[3], v[7], field[3], field[7]);
    __syncthreads();	// 等待线程同步
#else
    // 使用局部变量存储顶点插值结果
    float3 vertlist[12];

    vertlist[0] = vertexInterp(isoValue, v[0], v[1], field[0], field[1]);
    vertlist[1] = vertexInterp(isoValue, v[1], v[2], field[1], field[2]);
    vertlist[2] = vertexInterp(isoValue, v[2], v[3], field[2], field[3]);
    vertlist[3] = vertexInterp(isoValue, v[3], v[0], field[3], field[0]);

    vertlist[4] = vertexInterp(isoValue, v[4], v[5], field[4], field[5]);
    vertlist[5] = vertexInterp(isoValue, v[5], v[6], field[5], field[6]);
    vertlist[6] = vertexInterp(isoValue, v[6], v[7], field[6], field[7]);
    vertlist[7] = vertexInterp(isoValue, v[7], v[4], field[7], field[4]);

    vertlist[8] = vertexInterp(isoValue, v[0], v[4], field[0], field[4]);
    vertlist[9] = vertexInterp(isoValue, v[1], v[5], field[1], field[5]);
    vertlist[10] = vertexInterp(isoValue, v[2], v[6], field[2], field[6]);
    vertlist[11] = vertexInterp(isoValue, v[3], v[7], field[3], field[7]);
#endif

    // 查询三角形顶点数量
    uint numVerts = tex1Dfetch(numVertsTex, cubeindex);
	
	// 每3个顶点构建一个三角形
    for (int i=0; i<numVerts; i+=3)	
    {
        uint index = numVertsScanned[voxel] + i;

        float3 *v[3];
        uint edge;
        // 获取边的索引
        edge = tex1Dfetch(triTex, (cubeindex*16) + i);
#if USE_SHARED	// 根据边索引从顶点数组中找到相应顶点的地址
        v[0] = &vertlist[(edge*NTHREADS)+threadIdx.x];
#else
        v[0] = &vertlist[edge];
#endif

        edge = tex1Dfetch(triTex, (cubeindex*16) + i + 1);
#if USE_SHARED
        v[1] = &vertlist[(edge*NTHREADS)+threadIdx.x];
#else
        v[1] = &vertlist[edge];
#endif

        edge = tex1Dfetch(triTex, (cubeindex*16) + i + 2);
#if USE_SHARED
        v[2] = &vertlist[(edge*NTHREADS)+threadIdx.x];
#else
        v[2] = &vertlist[edge];
#endif

        // 计算三角形的平面法线
        float3 n = calcNormal(v[0], v[1], v[2]);

        // 将顶点位置和法线存储到输出数组
        if (index < (maxVerts - 3))
        {
            pos[index] = make_float4(*v[0], 1.0f);
            norm[index] = make_float4(n, 0.0f);

            pos[index+1] = make_float4(*v[1], 1.0f);
            norm[index+1] = make_float4(n, 0.0f);

            pos[index+2] = make_float4(*v[2], 1.0f);
            norm[index+2] = make_float4(n, 0.0f);
        }
    }
}

为什么使用make_float4而非make_float3？——大概率与OpenGL有关。
GPT：
1.内存对齐和访问效率：在 CUDA 和许多 GPU 硬件架构中，使用 float4 可以提高内存访问的对齐性和效率。
2.硬件优化：许多 GPU 处理器对 float4 类型有硬件级的优化，特别是在处理 SIMD操作时，可以同时处理四个浮点数。
3.齐次坐标的使用：在图形学中，使用 float4 是为了支持齐次坐标系。在三维图形渲染中，通常使用 4x4 的矩阵进行仿射变换。
4.与图形 API 的兼容性：大多数图形 API（例如 OpenGL、DirectX 等）都采用 float4 作为顶点数据格式，因为它可以表示位置、颜色、法线等属性，并与着色器程序中的输入类型一致。

它这里的插值也是线性插值

// compute interpolated vertex along an edge
__device__
float3 vertexInterp(float isolevel, float3 p0, float3 p1, float f0, float f1)
{
    float t = (isolevel - f0) / (f1 - f0);
    return lerp(p0, p1, t);
}

// lerp
// - linear interpolation between a and b, based on value t in [0, 1] range
inline __device__ __host__ float3 lerp(float3 a, float3 b, float t)
{
    return a + t*(b-a);
}

improvements

CUDA Samples 中的 Marching Cubes 示例进行了多项改进，以提高算法的性能和可扩展性。以下是一些主要的改进点：

使用 CUDA 的并行计算能力

并行处理每个体素: Marching Cubes 算法天然适合并行处理，因为每个体素可以独立计算。CUDA 实现通过将每个体素的处理分配给一个线程，使得大量体素可以同时进行处理，大幅提高了计算速度。
高效内存管理: 使用 CUDA 的共享内存、纹理内存等优化数据访问，减少了全局内存访问的开销，提高了整体的执行效率。

压缩空体素（Compact Voxel Array）

跳过空体素: 在 CUDA 实现中，使用了紧凑体素数组来存储活跃的体素（即需要生成三角形的体素），避免了对空体素的无效计算。这大幅减少了计算量。
扫描和压缩阶段: 利用前缀和扫描技术，标记并压缩所有有效体素的索引，从而跳过空体素。这种方法通过减少冗余计算进一步提高了效率。

使用查找表优化计算

查找表加速: Marching Cubes 算法中，顶点的插值和三角形生成依赖于查找表。CUDA 实现使用了纹理内存来存储这些查找表，并通过 tex1Dfetch 快速访问数据，减少了复杂的条件判断和重复计算。

Shared Memory 的利用

减少内存访问延迟: 在某些实现中，利用共享内存来存储临时计算结果，例如在计算三角形顶点位置时，将插值结果存储在共享内存中，从而减少对全局内存的访问，提高了线程块内计算的效率。

Surface Extraction 的改进

表面法线计算: CUDA 实现中不仅提取了表面的几何信息，还并行计算了法线。这一步骤对于后续的光照计算至关重要，也是在传统 CPU 实现中较为耗时的部分。
分块处理: 在生成三角形时，算法对网格进行分块处理，以保证处理的数据块不会超过 GPU 的限制。这种方法增强了算法的可扩展性，使其能够处理更大的数据集。

动态分配和映射资源

使用 CUDA 图形互操作: 样例代码展示了如何将生成的顶点和法线数据直接映射到 OpenGL VBO（顶点缓冲区对象）中，允许在 CUDA 计算后无缝渲染。这样避免了 CPU 和 GPU 之间的数据传输开销，提高了整体性能。

优化的线程调度

线程块的优化布局: 通过优化线程块和网格的配置，确保每个 CUDA 核心都能得到充分利用，减少了线程间的同步和调度开销。

通过这些改进，CUDA 实现的 Marching Cubes 比传统 CPU 实现要快得多，特别是在处理大规模数据时。样例代码展示了如何利用 GPU 的并行计算能力来加速常见的计算几何任务，同时保持高效的内存和资源管理。

补充一下关于gridSizeMask的使用，当时看到这个很懵，在CPU版本里没有接触到这个

掩码操作通常用于快速计算数组索引，特别是当索引基于坐标时。在这种情况下，掩码可以用来获取体素坐标中的低有效位（Least Significant Bits, LSBs）。掩码操作通常与按位与（bitwise AND）运算结合使用。
假设我们有一个三维坐标 (x, y, z)，并且我们想要计算这个坐标的体素索引。如果体素网格的大小为 32×32×32，那么我们可以使用掩码操作来获取每个坐标轴上的低 5 位，这会给出该坐标在体素网格中的位置。

例如，对于坐标 (17, 25, 9)，我们可以这样计算体素索引

uint3 gridSizeMask = make_uint3(31, 31, 31); // 31 = 00011111 in binary
uint x = 17;
uint y = 25;
uint z = 9;
uint index = ((x & gridSizeMask) * 32 * 32) + ((y & gridSizeMask) * 32) + (z & gridSizeMask);

这里，& 是按位与运算符，x & gridSizeMask 将保留 x 的低 5 位，并丢弃高 27 位。然后，我们通过>将每个坐标乘以适当的系数并将它们相加以计算出体素索引。
计算过程：

x & gridSizeMask = 17 & 31 = 00010001 & 00011111 = 00010001 // (二进制形式)，转换成十进制为 17。
y & gridSizeMask = 25 & 31 = 00011001 & 00011111 = 00011001 // (二进制形式)，转换成十进制为 25。
z & gridSizeMask = 9 & 31 = 00001001 & 00011111 = 00001001 // (二进制形式)，转换成十进制为 9。

接着，我们将这些值组合起来得到体素索引：

17 * 32 * 32 = 17632
25 * 32 = 800
9 = 9

因此，最终的体素索引为 17632 + 800 + 9 = 18441。
这种方式非常高效，因为它避免了除法运算，而只使用位操作和加法。在图形处理和并行计算中，这种方法非常常见，尤其是在需要快速访问内存的情况下。

整个Marching Cubes过程其实相当明了，但这个示例中涉及到了OpenGL图形绘制和CUDA相关操作，需要多花点时间琢磨。（怪虽已倒，但经验还需慢慢消化。）

先到这儿，你通关！