   CHECKED_CUDA(cudaMalloc(&deviceMem, scalarField_host.size()));
   CHECKED_CUDA(cudaMemcpy(deviceMem, scalarField_host.data(), scalarField_host.size(), cudaMemcpyHostToDevice));
 
+  std::vector<uint8_t> counts(256, 0);
+  std::vector<uint8_t> mctable(256 * 16, 0xffu);
+
+  for (unsigned c = 0; c < 256; c++) {
+    // 0 -> 000 = 3
+    // 1 -> 001 = 0
+    // 2 -> 010 = 7
+    // 3 -> 011 = 4
+    // 4 -> 100 = 2
+    // 5 -> 101 = 1
+    // 6 -> 110 = 6
+    // 7 -> 111 = 5
+    //    3--------7
+    //   /|       /|
+    //  / |      / |
+    // 2--------6  |
+    // |  1-----|--5
+    // | /      | /   y x
+    // |/       |/    |/
+    // 0--------4     ->z
+    //    4--------5
+    //   /|       /|
+    //  / |      / |
+    // 7--------6  |
+    // |  0-----|--1
+    // | /      | /
+    // |/       |/
+    // 3--------2
+    unsigned n =
+      (((1 << 3) & c) ? (1 << 0) : 0) |
+      (((1 << 0) & c) ? (1 << 1) : 0) |
+      (((1 << 7) & c) ? (1 << 2) : 0) |
+      (((1 << 4) & c) ? (1 << 3) : 0) |
+      (((1 << 2) & c) ? (1 << 4) : 0) |
+      (((1 << 1) & c) ? (1 << 5) : 0) |
+      (((1 << 6) & c) ? (1 << 6) : 0) |
+      (((1 << 5) & c) ? (1 << 7) : 0);
+    unsigned i = 0;
+    for (; i < 16; i++) {
+      const auto t = triangles[c][i];
+      if (t < 0) break;
+
+      unsigned vv[8] = {
+        1,  // 0
+        5,  // 1
+        4,  // 2
+        0,  // 3
+        3,  // 4
+        7,  // 5
+        6,  // 6
+        2,  // 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
+      };
+      mctable[16 * n + i] = (vv[e[t][0]] << 3) | vv[e[t][1]];
+      if (i % 3 == 2) {
+        std::swap(mctable[16 * n + i - 1], mctable[16 * n + i]);
+      }
+    }
+    counts[n] = i;
+  }
 
   // reference code
   fprintf(stderr, "%d %zu\n", nx * ny * nz * 2, scalarField_host.size());
   assert(nx * ny * nz * 4 == scalarField_host.size());
-  const float threshold = 100.f;
+  const float threshold = 0.f;
   std::vector<float> vertices;
   const float* f = reinterpret_cast<const float*>(scalarField_host.data());
   for (unsigned k = 0; k + 1 < nz; k++) {

          {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
        };
        const unsigned e[12][2] = {
          {0,1}, // 0

          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
      }
    }
  }
