 
 
 // compact voxel array
-__global__ void 
+__global__ void
 // __launch_bounds__ ( KERNTHREADS, 1 )
-compactVoxels(unsigned int * RESTRICT compactedVoxelArray, 
-              const uint2 * RESTRICT voxelOccupied, 
-              unsigned int lastVoxel, unsigned int numVoxels, 
-              unsigned int numVoxelsm1) {
-  unsigned int blockId = (blockIdx.y * gridDim.x) + blockIdx.x;
-  unsigned int i = (blockId * blockDim.x) + threadIdx.x;
-
-  if ((i < numVoxels) && ((i < numVoxelsm1) ? voxelOccupied[i].y < voxelOccupied[i+1].y : lastVoxel)) {
-    compactedVoxelArray[ voxelOccupied[i].y ] = i;
-  }
+compactVoxels(unsigned int* RESTRICT compactedVoxelArray, const uint2* RESTRICT voxelOccupied, unsigned int lastVoxel,
+    unsigned int numVoxels, unsigned int numVoxelsm1) {
+    unsigned int blockId = (blockIdx.y * gridDim.x) + blockIdx.x;
+    unsigned int i = (blockId * blockDim.x) + threadIdx.x;
+
+    if ((i < numVoxels) && ((i < numVoxelsm1) ? voxelOccupied[i].y < voxelOccupied[i + 1].y : lastVoxel)) {
+        compactedVoxelArray[voxelOccupied[i].y] = i;
+    }
 }
 
 
 // version that calculates no surface normal or color,  only triangle vertices
-__global__ void 
+__global__ void
 // __launch_bounds__ ( NTHREADS, 1 )
-generateTriangleVerticesSMEM(float3 * RESTRICT pos, 
-                             const unsigned int * RESTRICT compactedVoxelArray, 
-                             const uint2 * RESTRICT numVertsScanned, 
-                             const float * RESTRICT volume,
-                             uint3 gridSize, float3 voxelSize, 
-                             float isoValue, unsigned int activeVoxels, 
-                             unsigned int maxVertsM3) {
+generateTriangleVerticesSMEM(float3* RESTRICT pos, const unsigned int* RESTRICT compactedVoxelArray,
+    const uint2* RESTRICT numVertsScanned, const float* RESTRICT volume, uint3 gridSize, float3 voxelSize,
+    float isoValue, unsigned int activeVoxels, unsigned int maxVertsM3) {
     unsigned int zOffset = blockIdx.z * gridDim.x * gridDim.y;
     unsigned int blockId = (blockIdx.y * gridDim.x) + blockIdx.x;
     unsigned int i = zOffset * (blockDim.x * blockDim.y) + (blockId * blockDim.x) + threadIdx.x;

//     because we don't strictly enforce coalescing requirements in the
//     layout of the input volume presently.  If we forced X/Y dims to be
//     warp-multiple it would become possible to use wider fetches and
//     a few other tricks to improve global memory bandwidth 
__device__ float sampleVolume(const float * RESTRICT data, 
                              uint3 p, uint3 gridSize) {
    return data[(p.z*gridSize.x*gridSize.y) + (p.y*gridSize.x) + p.x];
}


// sample volume texture at a point p, p CAN NEVER BE OUT OF BOUNDS
__device__ float3 sampleColors(const float3 * RESTRICT data,
                               uint3 p, uint3 gridSize) {
    return data[(p.z*gridSize.x*gridSize.y) + (p.y*gridSize.x) + p.x];
}

// sample volume texture at a point p, p CAN NEVER BE OUT OF BOUNDS
__device__ float3 sampleColors(const uchar4 * RESTRICT data,
                               uint3 p, uint3 gridSize) {
    uchar4 cu = data[(p.z*gridSize.x*gridSize.y) + (p.y*gridSize.x) + p.x];
    float3 cf;
    convert_color(cf, cu);
    return cf;

    // Use shared memory to keep register pressure under control.
    // No need to call __syncthreads() since each thread uses its own
    // private shared memory buffer.
    __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]);

    // output triangle vertices
    unsigned int numVerts = tex1Dfetch(numVertsTex, cubeindex);
    for (int i = 0; i < numVerts; i += 3) {
        unsigned int index = numVertsScanned[voxel].x + i;

        float3* vert[3];
        int edge;
        edge = tex1Dfetch(triTex, (cubeindex * 16) + i);
        vert[0] = &vertlist[(edge * NTHREADS) + threadIdx.x];

        edge = tex1Dfetch(triTex, (cubeindex * 16) + i + 1);
        vert[1] = &vertlist[(edge * NTHREADS) + threadIdx.x];

        edge = tex1Dfetch(triTex, (cubeindex * 16) + i + 2);
        vert[2] = &vertlist[(edge * NTHREADS) + threadIdx.x];

        if (index < maxVertsM3) {
            pos[index] = *vert[0];
            pos[index + 1] = *vert[1];
            pos[index + 2] = *vert[2];
        }
    }
}


__global__ void
// __launch_bounds__ ( KERNTHREADS, 1 )
offsetTriangleVertices(float3* RESTRICT pos, float3 origin, unsigned int numVertsM1) {
    unsigned int zOffset = blockIdx.z * gridDim.x * gridDim.y;
    unsigned int blockId = (blockIdx.y * gridDim.x) + blockIdx.x;
    unsigned int i = zOffset + (blockId * blockDim.x) + threadIdx.x;

    if (i > numVertsM1)
        return;

    float3 p = pos[i];
    p.x += origin.x;
    p.y += origin.y;
    p.z += origin.z;
    pos[i] = p;
}


// version that calculates the surface normal for each triangle vertex
template<class NORMAL>
__global__ void
// __launch_bounds__ ( KERNTHREADS, 1 )
generateTriangleNormals(
    const float3* RESTRICT pos, NORMAL* norm, float3 gridSizeInv, float3 bBoxInv, unsigned int numVerts) {
    unsigned int zOffset = blockIdx.z * gridDim.x * gridDim.y;
    unsigned int blockId = (blockIdx.y * gridDim.x) + blockIdx.x;
    unsigned int i = zOffset + (blockId * blockDim.x) + threadIdx.x;

    if (i > numVerts - 1)
        return;

    float3 n;
    float3 p, p1, p2;
    // normal calculation using central differences
    // TODO
    //p = (pos[i] + make_float3(1.0f)) * 0.5f;
    p = pos[i];
    p.x *= bBoxInv.x;
    p.y *= bBoxInv.y;
    p.z *= bBoxInv.z;
    p1 = p + make_float3(gridSizeInv.x, 0.0f, 0.0f);
    p2 = p - make_float3(gridSizeInv.x, 0.0f, 0.0f);
    n.x = tex3D(volumeTex, p2.x, p2.y, p2.z) - tex3D(volumeTex, p1.x, p1.y, p1.z);
    p1 = p + make_float3(0.0f, gridSizeInv.y, 0.0f);
    p2 = p - make_float3(0.0f, gridSizeInv.y, 0.0f);
    n.y = tex3D(volumeTex, p2.x, p2.y, p2.z) - tex3D(volumeTex, p1.x, p1.y, p1.z);
    p1 = p + make_float3(0.0f, 0.0f, gridSizeInv.z);
    p2 = p - make_float3(0.0f, 0.0f, gridSizeInv.z);
    n.z = tex3D(volumeTex, p2.x, p2.y, p2.z) - tex3D(volumeTex, p1.x, p1.y, p1.z);

    NORMAL no;
    convert_normal(no, n);
    norm[i] = no;
}


// version that calculates the surface normal and color for each triangle vertex
template<class VERTEXCOL, class VOLTEX, class NORMAL>
__global__ void
// __launch_bounds__ ( KERNTHREADS, 1 )
generateTriangleColorNormal(const float3* RESTRICT pos, VERTEXCOL* RESTRICT col, NORMAL* RESTRICT norm,
    const VOLTEX* RESTRICT colors, uint3 gridSize, float3 gridSizeInv, float3 bBoxInv, unsigned int numVerts) {
    unsigned int zOffset = blockIdx.z * gridDim.x * gridDim.y;
    unsigned int blockId = (blockIdx.y * gridDim.x) + blockIdx.x;
    unsigned int i = zOffset + (blockId * blockDim.x) + threadIdx.x;

    if (i > numVerts - 1)
        return;

    float3 p = pos[i];
    p.x *= bBoxInv.x;
    p.y *= bBoxInv.y;
    p.z *= bBoxInv.z;
    // color computation
    float3 gridPosF = p;
    gridPosF.x *= float(gridSize.x);
    gridPosF.y *= float(gridSize.y);
    gridPosF.z *= float(gridSize.z);
    float3 gridPosFloor;
    // Without the offset, rounding errors can occur
    // TODO why do we need the offset??
    gridPosFloor.x = floorf(gridPosF.x + 0.0001f);
    gridPosFloor.y = floorf(gridPosF.y + 0.0001f);
    gridPosFloor.z = floorf(gridPosF.z + 0.0001f);
    float3 gridPosCeil;
    // Without the offset, rounding errors can occur
    // TODO why do we need the offset??
    gridPosCeil.x = ceilf(gridPosF.x - 0.0001f);
    gridPosCeil.y = ceilf(gridPosF.y - 0.0001f);
    gridPosCeil.z = ceilf(gridPosF.z - 0.0001f);
    uint3 gridPos0;
    gridPos0.x = gridPosFloor.x;
    gridPos0.y = gridPosFloor.y;
    gridPos0.z = gridPosFloor.z;
    uint3 gridPos1;
    gridPos1.x = gridPosCeil.x;
    gridPos1.y = gridPosCeil.y;
    gridPos1.z = gridPosCeil.z;

    // compute interpolated color
    float3 field[2];
    field[0] = sampleColors(colors, gridPos0, gridSize);
    field[1] = sampleColors(colors, gridPos1, gridSize);
    float3 tmp = gridPosF - gridPosFloor;
    float a = fmaxf(fmaxf(tmp.x, tmp.y), tmp.z);
    VERTEXCOL c;
    convert_color(c, lerp(field[0], field[1], a));
    col[i] = c;

    // normal calculation using central differences
    float3 p1, p2, n;
    p1 = p + make_float3(gridSizeInv.x, 0.0f, 0.0f);
    p2 = p - make_float3(gridSizeInv.x, 0.0f, 0.0f);
    n.x = tex3D(volumeTex, p2.x, p2.y, p2.z) - tex3D(volumeTex, p1.x, p1.y, p1.z);
    p1 = p + make_float3(0.0f, gridSizeInv.y, 0.0f);
    p2 = p - make_float3(0.0f, gridSizeInv.y, 0.0f);
    n.y = tex3D(volumeTex, p2.x, p2.y, p2.z) - tex3D(volumeTex, p1.x, p1.y, p1.z);
    p1 = p + make_float3(0.0f, 0.0f, gridSizeInv.z);
    p2 = p - make_float3(0.0f, 0.0f, gridSizeInv.z);
    n.z = tex3D(volumeTex, p2.x, p2.y, p2.z) - tex3D(volumeTex, p1.x, p1.y, p1.z);

    NORMAL no;
    convert_normal(no, n);
    norm[i] = no;
}


void allocateTextures(int** d_triTable, unsigned int** d_numVertsTable) {
    cudaChannelFormatDesc channelDescSigned = cudaCreateChannelDesc(32, 0, 0, 0, cudaChannelFormatKindSigned);
    cudaMalloc((void**)d_triTable, 256 * 16 * sizeof(int));
    cudaMemcpy((void*)*d_triTable, (void*)triTable, 256 * 16 * sizeof(int), cudaMemcpyHostToDevice);
    cudaBindTexture(0, triTex, *d_triTable, channelDescSigned);

    cudaChannelFormatDesc channelDescUnsigned = cudaCreateChannelDesc(32, 0, 0, 0, cudaChannelFormatKindUnsigned);
    cudaMalloc((void**)d_numVertsTable, 256 * sizeof(unsigned int));
    cudaMemcpy((void*)*d_numVertsTable, (void*)numVertsTable, 256 * sizeof(unsigned int), cudaMemcpyHostToDevice);
    cudaBindTexture(0, numVertsTex, *d_numVertsTable, channelDescUnsigned);
}


void bindVolumeTexture(cudaArray* d_vol, cudaChannelFormatDesc desc) {
    // set texture parameters
    volumeTex.normalized = 1;
    volumeTex.filterMode = cudaFilterModeLinear;
