     // 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];
+    __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]);
+    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) {
+    for (int i = 0; i < numVerts; i += 3) {
         unsigned int index = numVertsScanned[voxel].x + i;
 
-        float3 *vert[3];
+        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];
+        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];
+            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) {
+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);
+    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);
+    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) {
+void bindVolumeTexture(cudaArray* d_vol, cudaChannelFormatDesc desc) {
     // set texture parameters
     volumeTex.normalized = 1;
     volumeTex.filterMode = cudaFilterModeLinear;



// classify voxel based on number of vertices it will generate one thread per two voxels
template <int gridis3d, int subgrid>
__global__ void 
// __launch_bounds__ ( KERNTHREADS, 1 )
classifyVoxel(uint2 * RESTRICT voxelVerts, 
              const float * RESTRICT volume,
              uint3 gridSize, unsigned int numVoxels, float3 voxelSize,
              uint3 subGridStart, uint3 subGridEnd,
              float isoValue) {
    uint3 gridPos;
    unsigned int i;

    // Compute voxel indices and address using either 2-D or 3-D 
    // thread indexing depending on the caller-provided gridis3d parameter
    if (gridis3d) {
      // Compute voxel index using 3-D thread indexing
      // compute 3D grid position
      gridPos.x = blockIdx.x * blockDim.x + threadIdx.x;
      gridPos.y = blockIdx.y * blockDim.y + threadIdx.y;
      gridPos.z = blockIdx.z * blockDim.z + threadIdx.z;

      // safety check
      if (gridPos.x >= gridSize.x || 
          gridPos.y >= gridSize.y || 
          gridPos.z >= gridSize.z)
        return;

      // compute 1D grid index
      i = gridPos.z*gridSize.x*gridSize.y + gridPos.y*gridSize.x + gridPos.x;
    } else {
      // Compute voxel index using 2-D thread indexing
      unsigned int blockId = (blockIdx.y * gridDim.x) + blockIdx.x;
      i = (blockId * blockDim.x) + threadIdx.x;

      // safety check
      if (i >= numVoxels)
        return;

      // compute current grid position
      gridPos = calcGridPos(i, gridSize);
    }

    // If we are told to compute the isosurface for only a sub-region

 * \author Michael Krone <michael.krone@visus.uni-stuttgart.de>
 * \author John Stone <johns@ks.uiuc.edu>
 *
 * \copyright
 * (C) Copyright 1995-2019 The Board of Trustees of the
 * University of Illinois, All Rights Reserved \n
 * UIUC Open Source License:
 * http://www.ks.uiuc.edu/Research/vmd/plugins/pluginlicense.html
 */

