   for (unsigned k = 0; k + 1 < nz; k++) {
     for (unsigned j = 0; j + 1 < ny; j++) {
       for (unsigned i = 0; i + 1 < nx; i++) {
-
-        const unsigned v[8][3] = {
-          {0, 0, 1}, // 0
-          {1, 0, 1}, // 1
-          {1, 0, 0}, // 2
-          {0, 0, 0}, // 3 
-          {0, 1, 1}, // 4 
-          {1, 1, 1}, // 5
-          {1, 1, 0}, // 6
-          {0, 1, 0}, // 7
-        };
-
-       
-
-        unsigned vv[8] = {
-          0b001,  // 0
-          0b101,  // 1
-          0b100,  // 2
-          0b000,  // 3
-          0b011,  // 4
-          0b111,  // 5
-          0b110,  // 6
-          0b010,  // 7
-        };
-        const unsigned e[12][2] = {
-          {0,1}, // 0
-          {2,1}, // 1
-          {3,2}, // 2
-          {3,0}, // 3
-          {4,5}, // 4
-          {5,6}, // 5
-          {6,7}, // 6
-          {7,4}, // 7
-          {0,4}, // 8
-          {1,5}, // 9
-          {2,6}, // 10
-          {3,7}  // 11
-        };
-
-
-
-
-#if 1
         unsigned mccase =
           (f[((k + 0) * ny + (j + 0)) * nx + (i + 0)] < threshold ? (1 << 0) : 0) |
           (f[((k + 0) * ny + (j + 0)) * nx + (i + 1)] < threshold ? (1 << 1) : 0) |

          vertices.push_back((2 * j + 2*((r >> 4) & 1) + ((r >> 1) & 1)) / (2.f * ny));
          vertices.push_back((2 * k + 2*((r >> 5) & 1) + ((r >> 2) & 1)) / (2.f * nz));
        }
#else
        unsigned mccase =
          (f[((k + v[0][2]) * ny + (j + v[0][1])) * nx + (i + v[0][0])] < threshold ? (1 << 0) : 0) |
          (f[((k + v[1][2]) * ny + (j + v[1][1])) * nx + (i + v[1][0])] < threshold ? (1 << 1) : 0) |
          (f[((k + v[2][2]) * ny + (j + v[2][1])) * nx + (i + v[2][0])] < threshold ? (1 << 2) : 0) |
          (f[((k + v[3][2]) * ny + (j + v[3][1])) * nx + (i + v[3][0])] < threshold ? (1 << 3) : 0) |
          (f[((k + v[4][2]) * ny + (j + v[4][1])) * nx + (i + v[4][0])] < threshold ? (1 << 4) : 0) |
          (f[((k + v[5][2]) * ny + (j + v[5][1])) * nx + (i + v[5][0])] < threshold ? (1 << 5) : 0) |
          (f[((k + v[6][2]) * ny + (j + v[6][1])) * nx + (i + v[6][0])] < threshold ? (1 << 6) : 0) |
          (f[((k + v[7][2]) * ny + (j + v[7][1])) * nx + (i + v[7][0])] < threshold ? (1 << 7) : 0);
        for (unsigned l = 0; l < 15; l++) {
          const auto r = triangles[mccase][l];
          if (r < 0) break;

          vertices.push_back((2 * i + v[e[r][0]][0] + v[e[r][1]][0]) / (2.f * nx));
          vertices.push_back((2 * j + v[e[r][0]][1] + v[e[r][1]][1]) / (2.f * ny));
          vertices.push_back((2 * k + v[e[r][0]][2] + v[e[r][1]][2]) / (2.f * nz));
        }
#endif
      }
    }
  }

  CHECKED_CUDA(cudaMalloc(&deviceMem, scalarField_host.size()));
  CHECKED_CUDA(cudaMemcpy(deviceMem, scalarField_host.data(), scalarField_host.size(), cudaMemcpyHostToDevice));

  assert(nx * ny * nz * 4 == scalarField_host.size());
  const float threshold = 0.f;
  std::vector<float> vertices;
