VisuTwin Canvas
C++ 3D Engine — Metal Backend
Loading...
Searching...
No Matches
metalMarchingCubesPass.cpp
Go to the documentation of this file.
1// SPDX-License-Identifier: Apache-2.0
2// Copyright 2025-2026 Arnis Lektauers
3//
4// Metal compute pass for GPU Marching Cubes -- implementation.
5//
6// Contains the embedded MSL compute kernels (classifyCells + generateVertices)
7// and CPU-side dispatch logic for the two-pass atomic pipeline.
8//
9// Custom shader -- no upstream GLSL equivalent exists.
10//
11#include "metalMarchingCubesPass.h"
12
13#include "metalGraphicsDevice.h"
14#include "metalTexture.h"
15#include "spdlog/spdlog.h"
16
17#include <cstring>
18
19namespace visutwin::canvas
20{
21 namespace
22 {
23 // ── Embedded Metal Shading Language ─────────────────────────────
24 //
25 // Two compute kernels for GPU Marching Cubes isosurface extraction.
26 //
27 // Pass 1 (classifyCells): Classify each cell and count total vertices.
28 // Pass 2 (generateVertices): Generate vertex data with atomic allocation.
29 //
30 // The output vertex format matches the engine's standard VertexData
31 // layout (56 bytes): position(3) + normal(3) + uv0(2) + tangent(4) + uv1(2).
32 //
33 // Lookup tables (edgeTable, triTable) match the CPU implementation
34 // in viz/algorithms/marchingCubes.h (Lorensen & Cline convention).
35 //
36 constexpr const char* MC_COMPUTE_SOURCE = R"(
37#include <metal_stdlib>
38using namespace metal;
39
40// ── Uniform parameters ──────────────────────────────────────────────
41// Must match MCComputeParams in metalMarchingCubesPass.h (80 bytes).
42struct MCParams {
43 uint dimsX, dimsY, dimsZ;
44 float isovalue;
45 float domainMinX, domainMinY, domainMinZ;
46 float _pad0;
47 float domainMaxX, domainMaxY, domainMaxZ;
48 float _pad1;
49 float texelSizeX, texelSizeY, texelSizeZ;
50 uint maxVertices;
51 uint flipNormals;
52 float _pad2[3];
53};
54
55// ── Output vertex layout ────────────────────────────────────────────
56// 56 bytes, matches the engine's common.metal VertexData:
57// attribute(0) float3 position offset 0
58// attribute(1) float3 normal offset 12
59// attribute(2) float2 uv0 offset 24
60// attribute(3) float4 tangent offset 32
61// attribute(4) float2 uv1 offset 48
62struct VertexData {
63 packed_float3 position; // 12 bytes
64 packed_float3 normal; // 12 bytes
65 packed_float2 uv0; // 8 bytes
66 packed_float4 tangent; // 16 bytes
67 packed_float2 uv1; // 8 bytes
68};
69
70// ── Helpers ─────────────────────────────────────────────────────────
71
72// Sample volume at exact voxel coordinate (integer grid position).
73// Uses texel-center addressing: uvw = (gridPos + 0.5) / dims.
74inline float sampleVoxel(texture3d<float> volume, sampler s,
75 int ix, int iy, int iz,
76 float3 invDims)
77{
78 float3 uvw = (float3(ix, iy, iz) + 0.5f) * invDims;
79 return volume.sample(s, uvw).x;
80}
81
82// Convert grid coordinates to world position.
83inline float3 gridToWorld(int gx, int gy, int gz,
84 float3 domainMin, float3 domainMax,
85 float3 texelSize)
86{
87 return domainMin + float3(gx, gy, gz) * texelSize * (domainMax - domainMin);
88}
89
90// Interpolate vertex position along an edge between two corners.
91inline float3 interpolateEdge(float3 p0, float3 p1,
92 float v0, float v1, float iso)
93{
94 float denom = v1 - v0;
95 float t = (abs(denom) < 1.0e-10f) ? 0.5f : clamp((iso - v0) / denom, 0.0f, 1.0f);
96 return mix(p0, p1, t);
97}
98
99// Compute gradient normal at a world position via central differences.
100// Samples the volume at +/- one texel in each axis.
101inline float3 computeGradientNormal(texture3d<float> volume, sampler s,
102 float3 worldPos,
103 float3 domainMin, float3 invDomainSize,
104 float3 invDims, float3 halfTexel,
105 bool flipNormals)
106{
107 // Convert to UVW space, offset by half-texel in each direction
108 float3 uvw = (worldPos - domainMin) * invDomainSize;
109
110 float dx = volume.sample(s, uvw + float3(halfTexel.x, 0, 0)).x
111 - volume.sample(s, uvw - float3(halfTexel.x, 0, 0)).x;
112 float dy = volume.sample(s, uvw + float3(0, halfTexel.y, 0)).x
113 - volume.sample(s, uvw - float3(0, halfTexel.y, 0)).x;
114 float dz = volume.sample(s, uvw + float3(0, 0, halfTexel.z)).x
115 - volume.sample(s, uvw - float3(0, 0, halfTexel.z)).x;
116
117 // Gradient points low→high; normals should point outward (high→low)
118 float3 n = flipNormals ? float3(dx, dy, dz) : float3(-dx, -dy, -dz);
119
120 float len = length(n);
121 return (len > 1.0e-8f) ? (n / len) : float3(0.0f, 1.0f, 0.0f);
122}
123
124// Compute tangent from normal.
125inline float4 computeTangent(float3 normal)
126{
127 // Cross normal with world-up (0,1,0)
128 float3 c = cross(normal, float3(0.0f, 1.0f, 0.0f));
129 float len = length(c);
130
131 if (len < 1.0e-6f) {
132 // Normal parallel to up -- use right axis
133 c = cross(normal, float3(1.0f, 0.0f, 0.0f));
134 len = length(c);
135 }
136
137 float3 t = (len > 1.0e-6f) ? (c / len) : float3(1.0f, 0.0f, 0.0f);
138 return float4(t, 1.0f); // handedness = 1.0
139}
140
141// ── Pass 1: Classify cells and count total vertices ─────────────────
142
143kernel void classifyCells(
144 texture3d<float> volume [[texture(0)]],
145 sampler fieldSampler [[sampler(0)]],
146 constant MCParams& params [[buffer(0)]],
147 constant ushort* edgeTable [[buffer(1)]],
148 constant char* triTable [[buffer(2)]],
149 device atomic_uint* vertexCount [[buffer(3)]],
150 uint3 gid [[thread_position_in_grid]])
151{
152 // Each thread processes one cell. Grid covers (dims-1) cells per axis.
153 if (gid.x >= params.dimsX - 1 ||
154 gid.y >= params.dimsY - 1 ||
155 gid.z >= params.dimsZ - 1) return;
156
157 float3 invDims = float3(1.0f / float(params.dimsX),
158 1.0f / float(params.dimsY),
159 1.0f / float(params.dimsZ));
160
161 // Sample 8 corner values (Lorensen & Cline convention)
162 int x = int(gid.x);
163 int y = int(gid.y);
164 int z = int(gid.z);
165
166 float v0 = sampleVoxel(volume, fieldSampler, x, y, z, invDims);
167 float v1 = sampleVoxel(volume, fieldSampler, x+1, y, z, invDims);
168 float v2 = sampleVoxel(volume, fieldSampler, x+1, y+1, z, invDims);
169 float v3 = sampleVoxel(volume, fieldSampler, x, y+1, z, invDims);
170 float v4 = sampleVoxel(volume, fieldSampler, x, y, z+1, invDims);
171 float v5 = sampleVoxel(volume, fieldSampler, x+1, y, z+1, invDims);
172 float v6 = sampleVoxel(volume, fieldSampler, x+1, y+1, z+1, invDims);
173 float v7 = sampleVoxel(volume, fieldSampler, x, y+1, z+1, invDims);
174
175 float iso = params.isovalue;
176
177 // Build cube index
178 uint cubeIndex = 0;
179 if (v0 >= iso) cubeIndex |= 1u;
180 if (v1 >= iso) cubeIndex |= 2u;
181 if (v2 >= iso) cubeIndex |= 4u;
182 if (v3 >= iso) cubeIndex |= 8u;
183 if (v4 >= iso) cubeIndex |= 16u;
184 if (v5 >= iso) cubeIndex |= 32u;
185 if (v6 >= iso) cubeIndex |= 64u;
186 if (v7 >= iso) cubeIndex |= 128u;
187
188 // Skip if entirely inside or outside
189 ushort edges = edgeTable[cubeIndex];
190 if (edges == 0) return;
191
192 // Count vertices from triangle table
193 uint numVerts = 0;
194 for (int i = 0; i < 16; i += 3) {
195 if (triTable[cubeIndex * 16 + i] == -1) break;
196 numVerts += 3;
197 }
198
199 if (numVerts > 0) {
200 atomic_fetch_add_explicit(vertexCount, numVerts, memory_order_relaxed);
201 }
202}
203
204// ── Pass 2: Generate vertices with atomic allocation ────────────────
205
206kernel void generateVertices(
207 texture3d<float> volume [[texture(0)]],
208 sampler fieldSampler [[sampler(0)]],
209 constant MCParams& params [[buffer(0)]],
210 constant ushort* edgeTable [[buffer(1)]],
211 constant char* triTable [[buffer(2)]],
212 device atomic_uint* vertexCount [[buffer(3)]],
213 device VertexData* vertices [[buffer(4)]],
214 uint3 gid [[thread_position_in_grid]])
215{
216 if (gid.x >= params.dimsX - 1 ||
217 gid.y >= params.dimsY - 1 ||
218 gid.z >= params.dimsZ - 1) return;
219
220 float3 invDims = float3(1.0f / float(params.dimsX),
221 1.0f / float(params.dimsY),
222 1.0f / float(params.dimsZ));
223
224 int x = int(gid.x);
225 int y = int(gid.y);
226 int z = int(gid.z);
227
228 // Sample 8 corners
229 float v0 = sampleVoxel(volume, fieldSampler, x, y, z, invDims);
230 float v1 = sampleVoxel(volume, fieldSampler, x+1, y, z, invDims);
231 float v2 = sampleVoxel(volume, fieldSampler, x+1, y+1, z, invDims);
232 float v3 = sampleVoxel(volume, fieldSampler, x, y+1, z, invDims);
233 float v4 = sampleVoxel(volume, fieldSampler, x, y, z+1, invDims);
234 float v5 = sampleVoxel(volume, fieldSampler, x+1, y, z+1, invDims);
235 float v6 = sampleVoxel(volume, fieldSampler, x+1, y+1, z+1, invDims);
236 float v7 = sampleVoxel(volume, fieldSampler, x, y+1, z+1, invDims);
237
238 float iso = params.isovalue;
239
240 uint cubeIndex = 0;
241 if (v0 >= iso) cubeIndex |= 1u;
242 if (v1 >= iso) cubeIndex |= 2u;
243 if (v2 >= iso) cubeIndex |= 4u;
244 if (v3 >= iso) cubeIndex |= 8u;
245 if (v4 >= iso) cubeIndex |= 16u;
246 if (v5 >= iso) cubeIndex |= 32u;
247 if (v6 >= iso) cubeIndex |= 64u;
248 if (v7 >= iso) cubeIndex |= 128u;
249
250 ushort edges = edgeTable[cubeIndex];
251 if (edges == 0) return;
252
253 // Corner world positions
254 float3 domainMin = float3(params.domainMinX, params.domainMinY, params.domainMinZ);
255 float3 domainMax = float3(params.domainMaxX, params.domainMaxY, params.domainMaxZ);
256 float3 texelSize = float3(params.texelSizeX, params.texelSizeY, params.texelSizeZ);
257 float3 domainSize = domainMax - domainMin;
258
259 float3 p0 = domainMin + float3(x, y, z ) * texelSize * domainSize;
260 float3 p1 = domainMin + float3(x+1, y, z ) * texelSize * domainSize;
261 float3 p2 = domainMin + float3(x+1, y+1, z ) * texelSize * domainSize;
262 float3 p3 = domainMin + float3(x, y+1, z ) * texelSize * domainSize;
263 float3 p4 = domainMin + float3(x, y, z+1) * texelSize * domainSize;
264 float3 p5 = domainMin + float3(x+1, y, z+1) * texelSize * domainSize;
265 float3 p6 = domainMin + float3(x+1, y+1, z+1) * texelSize * domainSize;
266 float3 p7 = domainMin + float3(x, y+1, z+1) * texelSize * domainSize;
267
268 // Compute interpolated edge vertices
269 float3 edgeVerts[12];
270 if (edges & 0x001) edgeVerts[0] = interpolateEdge(p0, p1, v0, v1, iso);
271 if (edges & 0x002) edgeVerts[1] = interpolateEdge(p1, p2, v1, v2, iso);
272 if (edges & 0x004) edgeVerts[2] = interpolateEdge(p2, p3, v2, v3, iso);
273 if (edges & 0x008) edgeVerts[3] = interpolateEdge(p3, p0, v3, v0, iso);
274 if (edges & 0x010) edgeVerts[4] = interpolateEdge(p4, p5, v4, v5, iso);
275 if (edges & 0x020) edgeVerts[5] = interpolateEdge(p5, p6, v5, v6, iso);
276 if (edges & 0x040) edgeVerts[6] = interpolateEdge(p6, p7, v6, v7, iso);
277 if (edges & 0x080) edgeVerts[7] = interpolateEdge(p7, p4, v7, v4, iso);
278 if (edges & 0x100) edgeVerts[8] = interpolateEdge(p0, p4, v0, v4, iso);
279 if (edges & 0x200) edgeVerts[9] = interpolateEdge(p1, p5, v1, v5, iso);
280 if (edges & 0x400) edgeVerts[10] = interpolateEdge(p2, p6, v2, v6, iso);
281 if (edges & 0x800) edgeVerts[11] = interpolateEdge(p3, p7, v3, v7, iso);
282
283 // Gradient normal helpers
284 float3 invDomainSize = float3(1.0f / domainSize.x,
285 1.0f / domainSize.y,
286 1.0f / domainSize.z);
287 float3 halfTexel = 0.5f * invDims;
288 bool doFlip = (params.flipNormals != 0);
289
290 // Count vertices first to do a single atomic allocation
291 uint numVerts = 0;
292 for (int i = 0; i < 16; i += 3) {
293 if (triTable[cubeIndex * 16 + i] == -1) break;
294 numVerts += 3;
295 }
296
297 if (numVerts == 0) return;
298
299 // Atomically allocate space in the output buffer
300 uint baseIdx = atomic_fetch_add_explicit(vertexCount, numVerts, memory_order_relaxed);
301
302 // Safety check: don't write past buffer end
303 if (baseIdx + numVerts > params.maxVertices) return;
304
305 // Emit triangles
306 uint writeIdx = baseIdx;
307 for (int i = 0; i < 16; i += 3) {
308 int e0 = triTable[cubeIndex * 16 + i];
309 if (e0 == -1) break;
310 int e1 = triTable[cubeIndex * 16 + i + 1];
311 int e2 = triTable[cubeIndex * 16 + i + 2];
312
313 float3 positions[3] = { edgeVerts[e0], edgeVerts[e1], edgeVerts[e2] };
314
315 for (int vi = 0; vi < 3; ++vi) {
316 float3 pos = positions[vi];
317 float3 normal = computeGradientNormal(volume, fieldSampler,
318 pos, domainMin, invDomainSize,
319 invDims, halfTexel, doFlip);
320 float4 tangent = computeTangent(normal);
321
322 VertexData vd;
323 vd.position = packed_float3(pos);
324 vd.normal = packed_float3(normal);
325 vd.uv0 = packed_float2(0.0f, 0.0f);
326 vd.tangent = packed_float4(tangent);
327 vd.uv1 = packed_float2(0.0f, 0.0f);
328
329 vertices[writeIdx] = vd;
330 ++writeIdx;
331 }
332 }
333}
334)";
335
336 constexpr uint32_t THREADGROUP_SIZE_X = 4;
337 constexpr uint32_t THREADGROUP_SIZE_Y = 4;
338 constexpr uint32_t THREADGROUP_SIZE_Z = 4;
339
340 // ── Marching Cubes Lookup Tables ────────────────────────────────
341 // Identical to viz/algorithms/marchingCubes.h (Lorensen & Cline convention).
342
343 constexpr uint16_t kEdgeTable[256] = {
344 0x000, 0x109, 0x203, 0x30A, 0x406, 0x50F, 0x605, 0x70C,
345 0x80C, 0x905, 0xA0F, 0xB06, 0xC0A, 0xD03, 0xE09, 0xF00,
346 0x190, 0x099, 0x393, 0x29A, 0x596, 0x49F, 0x795, 0x69C,
347 0x99C, 0x895, 0xB9F, 0xA96, 0xD9A, 0xC93, 0xF99, 0xE90,
348 0x230, 0x339, 0x033, 0x13A, 0x636, 0x73F, 0x435, 0x53C,
349 0xA3C, 0xB35, 0x83F, 0x936, 0xE3A, 0xF33, 0xC39, 0xD30,
350 0x3A0, 0x2A9, 0x1A3, 0x0AA, 0x7A6, 0x6AF, 0x5A5, 0x4AC,
351 0xBAC, 0xAA5, 0x9AF, 0x8A6, 0xFAA, 0xEA3, 0xDA9, 0xCA0,
352 0x460, 0x569, 0x663, 0x76A, 0x066, 0x16F, 0x265, 0x36C,
353 0xC6C, 0xD65, 0xE6F, 0xF66, 0x86A, 0x963, 0xA69, 0xB60,
354 0x5F0, 0x4F9, 0x7F3, 0x6FA, 0x1F6, 0x0FF, 0x3F5, 0x2FC,
355 0xDFC, 0xCF5, 0xFFF, 0xEF6, 0x9FA, 0x8F3, 0xBF9, 0xAF0,
356 0x650, 0x759, 0x453, 0x55A, 0x256, 0x35F, 0x055, 0x15C,
357 0xE5C, 0xF55, 0xC5F, 0xD56, 0xA5A, 0xB53, 0x859, 0x950,
358 0x7C0, 0x6C9, 0x5C3, 0x4CA, 0x3C6, 0x2CF, 0x1C5, 0x0CC,
359 0xFCC, 0xEC5, 0xDCF, 0xCC6, 0xBCA, 0xAC3, 0x9C9, 0x8C0,
360 0x8C0, 0x9C9, 0xAC3, 0xBCA, 0xCC6, 0xDCF, 0xEC5, 0xFCC,
361 0x0CC, 0x1C5, 0x2CF, 0x3C6, 0x4CA, 0x5C3, 0x6C9, 0x7C0,
362 0x950, 0x859, 0xB53, 0xA5A, 0xD56, 0xC5F, 0xF55, 0xE5C,
363 0x15C, 0x055, 0x35F, 0x256, 0x55A, 0x453, 0x759, 0x650,
364 0xAF0, 0xBF9, 0x8F3, 0x9FA, 0xEF6, 0xFFF, 0xCF5, 0xDFC,
365 0x2FC, 0x3F5, 0x0FF, 0x1F6, 0x6FA, 0x7F3, 0x4F9, 0x5F0,
366 0xB60, 0xA69, 0x963, 0x86A, 0xF66, 0xE6F, 0xD65, 0xC6C,
367 0x36C, 0x265, 0x16F, 0x066, 0x76A, 0x663, 0x569, 0x460,
368 0xCA0, 0xDA9, 0xEA3, 0xFAA, 0x8A6, 0x9AF, 0xAA5, 0xBAC,
369 0x4AC, 0x5A5, 0x6AF, 0x7A6, 0x0AA, 0x1A3, 0x2A9, 0x3A0,
370 0xD30, 0xC39, 0xF33, 0xE3A, 0x936, 0x83F, 0xB35, 0xA3C,
371 0x53C, 0x435, 0x73F, 0x636, 0x13A, 0x033, 0x339, 0x230,
372 0xE90, 0xF99, 0xC93, 0xD9A, 0xA96, 0xB9F, 0x895, 0x99C,
373 0x69C, 0x795, 0x49F, 0x596, 0x29A, 0x393, 0x099, 0x190,
374 0xF00, 0xE09, 0xD03, 0xC0A, 0xB06, 0xA0F, 0x905, 0x80C,
375 0x70C, 0x605, 0x50F, 0x406, 0x30A, 0x203, 0x109, 0x000
376 };
377
378 constexpr int8_t kTriTable[256][16] = {
379 {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
380 { 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
381 { 0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
382 { 1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
383 { 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
384 { 0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
385 { 9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
386 { 2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1},
387 { 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
388 { 0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
389 { 1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
390 { 1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1},
391 { 3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
392 { 0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
393 { 3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1},
394 { 9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
395 { 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
396 { 4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
397 { 0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
398 { 4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1},
399 { 1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
400 { 3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1},
401 { 9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
402 { 2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1},
403 { 8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
404 {11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1},
405 { 9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
406 { 4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1},
407 { 3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1},
408 { 1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1},
409 { 4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1},
410 { 4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
411 { 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
412 { 9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
413 { 0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
414 { 8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1},
415 { 1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
416 { 3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
417 { 5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1},
418 { 2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1},
419 { 9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
420 { 0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
421 { 0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
422 { 2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1},
423 {10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1},
424 { 4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1},
425 { 5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1},
426 { 5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1},
427 { 9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
428 { 9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1},
429 { 0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1},
430 { 1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
431 { 9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1},
432 {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1},
433 { 8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1},
434 { 2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1},
435 { 7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1},
436 { 9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1},
437 { 2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1},
438 {11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1},
439 { 9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1},
440 { 5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1},
441 {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1},
442 {11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
443 {10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
444 { 0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
445 { 9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
446 { 1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
447 { 1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
448 { 1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1},
449 { 9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1},
450 { 5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1},
451 { 2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
452 {11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
453 { 0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
454 { 5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1},
455 { 6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1},
456 { 0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1},
457 { 3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1},
458 { 6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1},
459 { 5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
460 { 4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1},
461 { 1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
462 {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1},
463 { 6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1},
464 { 1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1},
465 { 8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1},
466 { 7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1},
467 { 3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
468 { 5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1},
469 { 0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1},
470 { 9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1},
471 { 8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1},
472 { 5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1},
473 { 0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1},
474 { 6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1},
475 {10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
476 { 4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1},
477 {10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1},
478 { 8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1},
479 { 1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1},
480 { 3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1},
481 { 0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
482 { 8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1},
483 {10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1},
484 { 0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1},
485 { 3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1},
486 { 6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1},
487 { 9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1},
488 { 8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1},
489 { 3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1},
490 { 6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
491 { 7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1},
492 { 0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1},
493 {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1},
494 {10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1},
495 { 1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1},
496 { 2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1},
497 { 7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1},
498 { 7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
499 { 2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1},
500 { 2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1},
501 { 1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1},
502 {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1},
503 { 8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1},
504 { 0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
505 { 7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1},
506 { 7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
507 { 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
508 { 3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
509 { 0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
510 { 8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
511 {10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
512 { 1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
513 { 2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
514 { 6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1},
515 { 7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
516 { 7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1},
517 { 2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1},
518 { 1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1},
519 {10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1},
520 {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1},
521 { 0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1},
522 { 7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1},
523 { 6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
524 { 3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1},
525 { 8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1},
526 { 9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1},
527 { 6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1},
528 { 1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1},
529 { 4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1},
530 {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1},
531 { 8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1},
532 { 0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
533 { 1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1},
534 { 1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1},
535 { 8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1},
536 {10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1},
537 { 4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1},
538 {10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
539 { 4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
540 { 0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
541 { 5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
542 {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1},
543 { 9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
544 { 6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1},
545 { 7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1},
546 { 3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1},
547 { 7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1},
548 { 9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1},
549 { 3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1},
550 { 6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1},
551 { 9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1},
552 { 1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1},
553 { 4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1},
554 { 7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1},
555 { 6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1},
556 { 3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1},
557 { 0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1},
558 { 6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1},
559 { 1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1},
560 { 0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1},
561 {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1},
562 { 6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1},
563 { 5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1},
564 { 9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1},
565 { 1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1},
566 { 1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
567 { 1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1},
568 {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1},
569 { 0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
570 {10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
571 {11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
572 {11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1},
573 { 5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1},
574 {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1},
575 {11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1},
576 { 0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1},
577 { 9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1},
578 { 7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1},
579 { 2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1},
580 { 8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1},
581 { 9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1},
582 { 9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1},
583 { 1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
584 { 0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1},
585 { 9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1},
586 { 9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
587 { 5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1},
588 { 5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1},
589 { 0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1},
590 {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1},
591 { 2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1},
592 { 0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1},
593 { 0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1},
594 { 9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
595 { 2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1},
596 { 5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1},
597 { 3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1},
598 { 5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1},
599 { 8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1},
600 { 0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
601 { 8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1},
602 { 9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
603 { 4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
604 { 0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1},
605 { 1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1},
606 { 3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1},
607 { 4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1},
608 { 9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1},
609 {11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1},
610 {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1},
611 { 2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1},
612 { 9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1},
613 { 3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1},
614 { 1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
615 { 4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1},
616 { 4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1},
617 { 4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
618 { 4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
619 { 9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
620 { 3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1},
621 { 0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
622 { 3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
623 { 1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1},
624 { 3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1},
625 { 0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
626 { 3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
627 { 2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1},
628 { 9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
629 { 2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1},
630 { 1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
631 { 1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
632 { 0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
633 { 0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
634 {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}
635 };
636
637 } // anonymous namespace
638
639 // ─── Construction / Destruction ───────────────────────────────────
640
641 MetalMarchingCubesPass::MetalMarchingCubesPass(MetalGraphicsDevice* device)
642 : device_(device)
643 {
644 }
645
646 MetalMarchingCubesPass::~MetalMarchingCubesPass()
647 {
648 if (classifyPipeline_) { classifyPipeline_->release(); classifyPipeline_ = nullptr; }
649 if (generatePipeline_) { generatePipeline_->release(); generatePipeline_ = nullptr; }
650 if (edgeTableBuffer_) { edgeTableBuffer_->release(); edgeTableBuffer_ = nullptr; }
651 if (triTableBuffer_) { triTableBuffer_->release(); triTableBuffer_ = nullptr; }
652 if (counterBuffer_) { counterBuffer_->release(); counterBuffer_ = nullptr; }
653 if (uniformBuffer_) { uniformBuffer_->release(); uniformBuffer_ = nullptr; }
654 if (fieldSampler_) { fieldSampler_->release(); fieldSampler_ = nullptr; }
655 }
656
657 // ─── Lazy Resource Creation ──────────────────────────────────────
658
659 void MetalMarchingCubesPass::ensureResources()
660 {
661 if (resourcesReady_) return;
662
663 auto* mtlDevice = device_->raw();
664 if (!mtlDevice) return;
665
666 // ── Compile compute shaders ────────────────────────────────
667 if (!classifyPipeline_ || !generatePipeline_) {
668 NS::Error* error = nullptr;
669 auto* source = NS::String::string(MC_COMPUTE_SOURCE, NS::UTF8StringEncoding);
670 auto* library = mtlDevice->newLibrary(source, nullptr, &error);
671 if (!library) {
672 spdlog::error("[MetalMarchingCubesPass] Failed to compile MC shaders: {}",
673 error ? error->localizedDescription()->utf8String() : "unknown");
674 return;
675 }
676
677 // Classify pipeline
678 if (!classifyPipeline_) {
679 auto* funcName = NS::String::string("classifyCells", NS::UTF8StringEncoding);
680 auto* function = library->newFunction(funcName);
681 if (!function) {
682 spdlog::error("[MetalMarchingCubesPass] Entry point 'classifyCells' not found");
683 library->release();
684 return;
685 }
686 classifyPipeline_ = mtlDevice->newComputePipelineState(function, &error);
687 function->release();
688 if (!classifyPipeline_) {
689 spdlog::error("[MetalMarchingCubesPass] Failed to create classify pipeline: {}",
690 error ? error->localizedDescription()->utf8String() : "unknown");
691 library->release();
692 return;
693 }
694 }
695
696 // Generate pipeline
697 if (!generatePipeline_) {
698 auto* funcName = NS::String::string("generateVertices", NS::UTF8StringEncoding);
699 auto* function = library->newFunction(funcName);
700 if (!function) {
701 spdlog::error("[MetalMarchingCubesPass] Entry point 'generateVertices' not found");
702 library->release();
703 return;
704 }
705 generatePipeline_ = mtlDevice->newComputePipelineState(function, &error);
706 function->release();
707 if (!generatePipeline_) {
708 spdlog::error("[MetalMarchingCubesPass] Failed to create generate pipeline: {}",
709 error ? error->localizedDescription()->utf8String() : "unknown");
710 library->release();
711 return;
712 }
713 }
714
715 library->release();
716 }
717
718 // ── Lookup table buffers ───────────────────────────────────
719 if (!edgeTableBuffer_) {
720 edgeTableBuffer_ = mtlDevice->newBuffer(
721 kEdgeTable, sizeof(kEdgeTable), MTL::ResourceStorageModeShared);
722 if (!edgeTableBuffer_) {
723 spdlog::error("[MetalMarchingCubesPass] Failed to create edgeTable buffer");
724 return;
725 }
726 }
727
728 if (!triTableBuffer_) {
729 triTableBuffer_ = mtlDevice->newBuffer(
730 kTriTable, sizeof(kTriTable), MTL::ResourceStorageModeShared);
731 if (!triTableBuffer_) {
732 spdlog::error("[MetalMarchingCubesPass] Failed to create triTable buffer");
733 return;
734 }
735 }
736
737 // ── Atomic counter buffer (single uint32) ──────────────────
738 if (!counterBuffer_) {
739 counterBuffer_ = mtlDevice->newBuffer(
740 sizeof(uint32_t), MTL::ResourceStorageModeShared);
741 if (!counterBuffer_) {
742 spdlog::error("[MetalMarchingCubesPass] Failed to create counter buffer");
743 return;
744 }
745 }
746
747 // ── Uniform buffer (MCComputeParams) ───────────────────────
748 if (!uniformBuffer_) {
749 uniformBuffer_ = mtlDevice->newBuffer(
750 sizeof(MCComputeParams), MTL::ResourceStorageModeShared);
751 if (!uniformBuffer_) {
752 spdlog::error("[MetalMarchingCubesPass] Failed to create uniform buffer");
753 return;
754 }
755 }
756
757 // ── Nearest-neighbor sampler for exact voxel reads ─────────
758 if (!fieldSampler_) {
759 auto* desc = MTL::SamplerDescriptor::alloc()->init();
760 desc->setMinFilter(MTL::SamplerMinMagFilterNearest);
761 desc->setMagFilter(MTL::SamplerMinMagFilterNearest);
762 desc->setSAddressMode(MTL::SamplerAddressModeClampToEdge);
763 desc->setTAddressMode(MTL::SamplerAddressModeClampToEdge);
764 desc->setRAddressMode(MTL::SamplerAddressModeClampToEdge);
765 fieldSampler_ = mtlDevice->newSamplerState(desc);
766 desc->release();
767 }
768
769 resourcesReady_ = (classifyPipeline_ && generatePipeline_ &&
770 edgeTableBuffer_ && triTableBuffer_ &&
771 counterBuffer_ && uniformBuffer_ && fieldSampler_);
772
773 if (resourcesReady_) {
774 spdlog::info("[MetalMarchingCubesPass] Resources initialized successfully");
775 }
776 }
777
778 // ─── Per-Extraction Buffer Allocation ───────────────────────────
779
780 MTL::Buffer* MetalMarchingCubesPass::allocateVertexBuffer(uint32_t vertexCount)
781 {
782 auto* mtlDevice = device_->raw();
783 if (!mtlDevice) return nullptr;
784
785 constexpr size_t VERTEX_SIZE = 56; // sizeof(VertexData) in MSL
786 const size_t bufferSize = static_cast<size_t>(vertexCount) * VERTEX_SIZE;
787
788 auto* buffer = mtlDevice->newBuffer(bufferSize, MTL::ResourceStorageModeShared);
789 if (!buffer) {
790 spdlog::error("[MetalMarchingCubesPass] Failed to allocate vertex buffer ({} verts, {:.1f} MB)",
791 vertexCount, static_cast<double>(bufferSize) / (1024.0 * 1024.0));
792 }
793
794 return buffer; // Caller takes ownership (newBuffer returns +1 ref)
795 }
796
797 // ─── GPU Extraction ──────────────────────────────────────────────
798
799 MCExtractResult MetalMarchingCubesPass::extract(Texture* volumeTexture,
800 const MCComputeParams& params)
801 {
802 MCExtractResult result;
803
804 if (!volumeTexture) {
805 spdlog::warn("[MetalMarchingCubesPass] No volume texture provided");
806 return result;
807 }
808
809 ensureResources();
810 if (!resourcesReady_) return result;
811
812 // Get the underlying Metal texture
813 auto* hwTexture = dynamic_cast<gpu::MetalTexture*>(volumeTexture->impl());
814 if (!hwTexture || !hwTexture->raw()) {
815 spdlog::warn("[MetalMarchingCubesPass] Volume texture has no Metal backing");
816 return result;
817 }
818
819 // Upload uniforms
820 std::memcpy(uniformBuffer_->contents(), &params, sizeof(MCComputeParams));
821
822 // Reset atomic counter to 0
823 uint32_t zero = 0;
824 std::memcpy(counterBuffer_->contents(), &zero, sizeof(uint32_t));
825
826 // Grid dimensions for dispatch (one thread per cell)
827 uint32_t cellsX = params.dimsX - 1;
828 uint32_t cellsY = params.dimsY - 1;
829 uint32_t cellsZ = params.dimsZ - 1;
830
831 auto threadgroupsX = (cellsX + THREADGROUP_SIZE_X - 1) / THREADGROUP_SIZE_X;
832 auto threadgroupsY = (cellsY + THREADGROUP_SIZE_Y - 1) / THREADGROUP_SIZE_Y;
833 auto threadgroupsZ = (cellsZ + THREADGROUP_SIZE_Z - 1) / THREADGROUP_SIZE_Z;
834
835 MTL::Size threadgroups(threadgroupsX, threadgroupsY, threadgroupsZ);
836 MTL::Size threadsPerGroup(THREADGROUP_SIZE_X, THREADGROUP_SIZE_Y, THREADGROUP_SIZE_Z);
837
838 // ── Pass 1: Classify cells and count vertices ────────────
839 // No vertex buffer needed — only counts via atomic counter.
840 {
841 auto* commandBuffer = device_->_commandQueue->commandBuffer();
842 if (!commandBuffer) {
843 spdlog::warn("[MetalMarchingCubesPass] Failed to create command buffer for Pass 1");
844 return result;
845 }
846
847 auto* encoder = commandBuffer->computeCommandEncoder();
848 if (!encoder) {
849 spdlog::warn("[MetalMarchingCubesPass] Failed to create compute encoder for Pass 1");
850 return result;
851 }
852
853 encoder->pushDebugGroup(
854 NS::String::string("MC-Classify", NS::UTF8StringEncoding));
855
856 encoder->setComputePipelineState(classifyPipeline_);
857 encoder->setTexture(hwTexture->raw(), 0); // [[texture(0)]]
858 encoder->setSamplerState(fieldSampler_, 0); // [[sampler(0)]]
859 encoder->setBuffer(uniformBuffer_, 0, 0); // [[buffer(0)]]
860 encoder->setBuffer(edgeTableBuffer_, 0, 1); // [[buffer(1)]]
861 encoder->setBuffer(triTableBuffer_, 0, 2); // [[buffer(2)]]
862 encoder->setBuffer(counterBuffer_, 0, 3); // [[buffer(3)]]
863
864 encoder->dispatchThreadgroups(threadgroups, threadsPerGroup);
865
866 encoder->popDebugGroup();
867 encoder->endEncoding();
868 commandBuffer->commit();
869 commandBuffer->waitUntilCompleted();
870 }
871
872 // Read back vertex count from atomic counter
873 uint32_t totalVertices = 0;
874 std::memcpy(&totalVertices, counterBuffer_->contents(), sizeof(uint32_t));
875
876 if (totalVertices == 0) {
877 spdlog::debug("[MetalMarchingCubesPass] No vertices generated (isovalue={:.4f})",
878 params.isovalue);
879 result.vertexCount = 0;
880 result.success = true;
881 return result;
882 }
883
884 spdlog::debug("[MetalMarchingCubesPass] Pass 1: {} vertices ({} triangles)",
885 totalVertices, totalVertices / 3);
886
887 // ── Allocate a fresh vertex buffer for this extraction ────
888 // Each extract() call produces its own buffer — caller takes ownership.
889 auto* vbuffer = allocateVertexBuffer(totalVertices);
890 if (!vbuffer) return result;
891
892 // Update maxVertices in uniforms to actual count (for safety check in shader)
893 MCComputeParams pass2Params = params;
894 pass2Params.maxVertices = totalVertices;
895 std::memcpy(uniformBuffer_->contents(), &pass2Params, sizeof(MCComputeParams));
896
897 // Reset counter for Pass 2 (used for atomic allocation)
898 std::memcpy(counterBuffer_->contents(), &zero, sizeof(uint32_t));
899
900 // ── Pass 2: Generate vertices ────────────────────────────
901 {
902 auto* commandBuffer = device_->_commandQueue->commandBuffer();
903 if (!commandBuffer) {
904 spdlog::warn("[MetalMarchingCubesPass] Failed to create command buffer for Pass 2");
905 vbuffer->release();
906 return result;
907 }
908
909 auto* encoder = commandBuffer->computeCommandEncoder();
910 if (!encoder) {
911 spdlog::warn("[MetalMarchingCubesPass] Failed to create compute encoder for Pass 2");
912 vbuffer->release();
913 return result;
914 }
915
916 encoder->pushDebugGroup(
917 NS::String::string("MC-Generate", NS::UTF8StringEncoding));
918
919 encoder->setComputePipelineState(generatePipeline_);
920 encoder->setTexture(hwTexture->raw(), 0); // [[texture(0)]]
921 encoder->setSamplerState(fieldSampler_, 0); // [[sampler(0)]]
922 encoder->setBuffer(uniformBuffer_, 0, 0); // [[buffer(0)]]
923 encoder->setBuffer(edgeTableBuffer_, 0, 1); // [[buffer(1)]]
924 encoder->setBuffer(triTableBuffer_, 0, 2); // [[buffer(2)]]
925 encoder->setBuffer(counterBuffer_, 0, 3); // [[buffer(3)]]
926 encoder->setBuffer(vbuffer, 0, 4); // [[buffer(4)]]
927
928 encoder->dispatchThreadgroups(threadgroups, threadsPerGroup);
929
930 encoder->popDebugGroup();
931 encoder->endEncoding();
932 commandBuffer->commit();
933 commandBuffer->waitUntilCompleted();
934 }
935
936 // Transfer ownership to caller — they must release or adopt this buffer.
937 result.vertexBuffer = vbuffer;
938 result.vertexCount = totalVertices;
939 result.success = true;
940
941 return result;
942 }
943
944 // ─── Multi-Isovalue Batch Extraction ────────────────────────────
945
946 MCBatchResult MetalMarchingCubesPass::extractBatch(
947 Texture* volumeTexture,
948 const MCComputeParams& baseParams,
949 const std::vector<float>& isovalues,
950 const std::vector<bool>& flipNormals)
951 {
952 MCBatchResult batchResult;
953
954 if (!volumeTexture || isovalues.empty()) {
955 spdlog::warn("[MetalMarchingCubesPass] extractBatch: invalid args");
956 return batchResult;
957 }
958
959 ensureResources();
960 if (!resourcesReady_) return batchResult;
961
962 auto* hwTexture = dynamic_cast<gpu::MetalTexture*>(volumeTexture->impl());
963 if (!hwTexture || !hwTexture->raw()) {
964 spdlog::warn("[MetalMarchingCubesPass] Volume texture has no Metal backing");
965 return batchResult;
966 }
967
968 // Grid dimensions (shared across all layers)
969 const uint32_t cellsX = baseParams.dimsX - 1;
970 const uint32_t cellsY = baseParams.dimsY - 1;
971 const uint32_t cellsZ = baseParams.dimsZ - 1;
972
973 const auto threadgroupsX = (cellsX + THREADGROUP_SIZE_X - 1) / THREADGROUP_SIZE_X;
974 const auto threadgroupsY = (cellsY + THREADGROUP_SIZE_Y - 1) / THREADGROUP_SIZE_Y;
975 const auto threadgroupsZ = (cellsZ + THREADGROUP_SIZE_Z - 1) / THREADGROUP_SIZE_Z;
976
977 MTL::Size threadgroups(threadgroupsX, threadgroupsY, threadgroupsZ);
978 MTL::Size threadsPerGroup(THREADGROUP_SIZE_X, THREADGROUP_SIZE_Y, THREADGROUP_SIZE_Z);
979
980 const size_t numLayers = isovalues.size();
981 batchResult.layers.resize(numLayers);
982
983 // ── Pass 1: Count vertices for all layers ────────────────
984 std::vector<uint32_t> vertexCounts(numLayers, 0);
985
986 for (size_t i = 0; i < numLayers; ++i) {
987 MCComputeParams layerParams = baseParams;
988 layerParams.isovalue = isovalues[i];
989 if (!flipNormals.empty() && i < flipNormals.size()) {
990 layerParams.flipNormals = flipNormals[i] ? 1 : 0;
991 }
992
993 std::memcpy(uniformBuffer_->contents(), &layerParams, sizeof(MCComputeParams));
994
995 uint32_t zero = 0;
996 std::memcpy(counterBuffer_->contents(), &zero, sizeof(uint32_t));
997
998 auto* commandBuffer = device_->_commandQueue->commandBuffer();
999 if (!commandBuffer) continue;
1000 auto* encoder = commandBuffer->computeCommandEncoder();
1001 if (!encoder) continue;
1002
1003 encoder->setComputePipelineState(classifyPipeline_);
1004 encoder->setTexture(hwTexture->raw(), 0);
1005 encoder->setSamplerState(fieldSampler_, 0);
1006 encoder->setBuffer(uniformBuffer_, 0, 0);
1007 encoder->setBuffer(edgeTableBuffer_, 0, 1);
1008 encoder->setBuffer(triTableBuffer_, 0, 2);
1009 encoder->setBuffer(counterBuffer_, 0, 3);
1010 encoder->dispatchThreadgroups(threadgroups, threadsPerGroup);
1011 encoder->endEncoding();
1012 commandBuffer->commit();
1013 commandBuffer->waitUntilCompleted();
1014
1015 std::memcpy(&vertexCounts[i], counterBuffer_->contents(), sizeof(uint32_t));
1016 }
1017
1018 // ── Allocate per-layer vertex buffers ────────────────────
1019 for (size_t i = 0; i < numLayers; ++i) {
1020 if (vertexCounts[i] > 0) {
1021 batchResult.layers[i].vertexBuffer = allocateVertexBuffer(vertexCounts[i]);
1022 }
1023 }
1024
1025 // ── Pass 2: Generate vertices for all layers ─────────────
1026 for (size_t i = 0; i < numLayers; ++i) {
1027 if (vertexCounts[i] == 0 || !batchResult.layers[i].vertexBuffer) {
1028 batchResult.layers[i].vertexCount = 0;
1029 continue;
1030 }
1031
1032 MCComputeParams layerParams = baseParams;
1033 layerParams.isovalue = isovalues[i];
1034 layerParams.maxVertices = vertexCounts[i];
1035 if (!flipNormals.empty() && i < flipNormals.size()) {
1036 layerParams.flipNormals = flipNormals[i] ? 1 : 0;
1037 }
1038
1039 std::memcpy(uniformBuffer_->contents(), &layerParams, sizeof(MCComputeParams));
1040
1041 uint32_t zero = 0;
1042 std::memcpy(counterBuffer_->contents(), &zero, sizeof(uint32_t));
1043
1044 auto* commandBuffer = device_->_commandQueue->commandBuffer();
1045 if (!commandBuffer) continue;
1046 auto* encoder = commandBuffer->computeCommandEncoder();
1047 if (!encoder) continue;
1048
1049 encoder->setComputePipelineState(generatePipeline_);
1050 encoder->setTexture(hwTexture->raw(), 0);
1051 encoder->setSamplerState(fieldSampler_, 0);
1052 encoder->setBuffer(uniformBuffer_, 0, 0);
1053 encoder->setBuffer(edgeTableBuffer_, 0, 1);
1054 encoder->setBuffer(triTableBuffer_, 0, 2);
1055 encoder->setBuffer(counterBuffer_, 0, 3);
1056 encoder->setBuffer(batchResult.layers[i].vertexBuffer, 0, 4);
1057 encoder->dispatchThreadgroups(threadgroups, threadsPerGroup);
1058 encoder->endEncoding();
1059 commandBuffer->commit();
1060 commandBuffer->waitUntilCompleted();
1061
1062 batchResult.layers[i].vertexCount = vertexCounts[i];
1063 }
1064
1065 batchResult.success = true;
1066 return batchResult;
1067 }
1068
1069} // namespace visutwin::canvas
visutwin::canvas::MetalMarchingCubesPass::extractBatch
MCBatchResult extractBatch(Texture *volumeTexture, const MCComputeParams &baseParams, const std::vector< float > &isovalues, const std::vector< bool > &flipNormals={})
Definition metalMarchingCubesPass.cpp:946
visutwin::canvas::MetalMarchingCubesPass::extract
MCExtractResult extract(Texture *volumeTexture, const MCComputeParams &params)
Definition metalMarchingCubesPass.cpp:799
visutwin::canvas::Texture
GPU texture resource supporting 2D, cubemap, volume, and array formats with mipmap management.
Definition texture.h:57
visutwin::canvas::Texture::impl
gpu::HardwareTexture * impl() const
Definition texture.h:101
visutwin::canvas::MCBatchResult
Result of a multi-isovalue batch extraction.
Definition metalMarchingCubesPass.h:63
visutwin::canvas::MCComputeParams::maxVertices
uint32_t maxVertices
Safety cap for output buffer.
Definition metalMarchingCubesPass.h:44
visutwin::canvas::MCComputeParams::dimsZ
uint32_t dimsZ
Volume dimensions in voxels.
Definition metalMarchingCubesPass.h:37
visutwin::canvas::MCComputeParams::flipNormals
uint32_t flipNormals
0 = outward (low-to-high), 1 = inward.
Definition metalMarchingCubesPass.h:45
visutwin::canvas::MCExtractResult::vertexCount
uint32_t vertexCount
Number of vertices (triangles*3).
Definition metalMarchingCubesPass.h:57
visutwin::canvas::MCExtractResult::vertexBuffer
MTL::Buffer * vertexBuffer
56-byte VertexData per vertex. Caller owns.
Definition metalMarchingCubesPass.h:56