|
7 | 7 | import { extractBrain, applyMask, largestComponent, fillHoles3D, smoothMask3D, dilate3D } from './bet.js'; |
8 | 8 | import { segmentTissue } from './tissue.js'; |
9 | 9 | import { registerAffine, invert4x4, resampleLabels } from '../chop/chop.js'; |
10 | | -import { conformVol, probToNative, BET_SHAPE } from './conform.js'; |
| 10 | +import { conformVol, probToNative, labelToNative, BET_SHAPE } from './conform.js'; |
11 | 11 |
|
12 | 12 | const P = (pct, msg) => self.postMessage({ type: 'progress', pct, msg }); |
13 | 13 | const ORT_VER = '1.20.1', ORT_CDN = `https://cdn.jsdelivr.net/npm/onnxruntime-web@${ORT_VER}/dist/`; |
14 | | -let _ort = null, _sess = null; |
15 | | -async function loadModel() { |
16 | | - if (_sess) return _sess; |
17 | | - if (!_ort) { _ort = await import(ORT_CDN + 'ort.wasm.min.mjs'); _ort.env.wasm.wasmPaths = ORT_CDN; } |
18 | | - const url = new URL('model/bet_unet.onnx', import.meta.url); |
19 | | - const buf = await (await fetch(url)).arrayBuffer(); |
20 | | - _sess = await _ort.InferenceSession.create(buf, { executionProviders: ['wasm'] }); |
| 14 | +let _ort = null, _sess = null, _sessBackend = null, _buf = null, _backend = 'wasm', _last = null; |
| 15 | +// Backend is user-selectable (Compute dropdown) so it can be A/B'd on a real GPU: |
| 16 | +// • WASM (CPU): verified identical to the PyTorch reference (agreement 99.98%, Dice 0.9998). Default. |
| 17 | +// • WebGPU (GPU): faster on real hardware, but kernels/precision can differ by device. |
| 18 | +// We create the session with the chosen EP ONLY (no silent fallback) so the reported backend is honest. |
| 19 | +async function loadModel(want) { |
| 20 | + want = (want === 'webgpu') ? 'webgpu' : 'wasm'; |
| 21 | + if (_sess && _sessBackend === want) return _sess; |
| 22 | + // the webgpu bundle ships BOTH backends, so one module serves either choice |
| 23 | + if (!_ort) { |
| 24 | + _ort = await import(ORT_CDN + 'ort.webgpu.min.mjs'); _ort.env.wasm.wasmPaths = ORT_CDN; |
| 25 | + _ort.env.wasm.numThreads = (self.crossOriginIsolated && navigator.hardwareConcurrency) |
| 26 | + ? Math.min(8, navigator.hardwareConcurrency) : 1; // multi-threaded WASM if cross-origin isolated |
| 27 | + } |
| 28 | + if (_sess) { try { await _sess.release?.(); } catch (e) {} _sess = null; } // backend switched → recreate |
| 29 | + if (!_buf) { const url = new URL('model/bet_unet_hires.onnx', import.meta.url); _buf = await (await fetch(url)).arrayBuffer(); } |
| 30 | + if (want === 'webgpu') { |
| 31 | + try { |
| 32 | + _sess = await _ort.InferenceSession.create(_buf, { executionProviders: ['webgpu'] }); |
| 33 | + _sessBackend = 'webgpu'; _backend = 'webgpu'; return _sess; |
| 34 | + } catch (e) { _backend = 'wasm (WebGPU unavailable)'; } // no adapter → honest note, fall to wasm |
| 35 | + } |
| 36 | + _sess = await _ort.InferenceSession.create(_buf, { executionProviders: ['wasm'] }); |
| 37 | + _sessBackend = 'wasm'; |
| 38 | + if (want === 'wasm') _backend = _ort.env.wasm.numThreads > 1 ? `wasm×${_ort.env.wasm.numThreads}` : 'wasm'; |
21 | 39 | return _sess; |
22 | 40 | } |
23 | 41 | const affToMT = a => ({ M: [a[0][0], a[0][1], a[0][2], a[1][0], a[1][1], a[1][2], a[2][0], a[2][1], a[2][2]], T: [a[0][3], a[1][3], a[2][3]] }); |
24 | 42 |
|
25 | | -async function learnedMask(vol) { // vol: {data,dims,affine} -> Uint8 native mask |
26 | | - P(20, 'Loading AI model…'); const sess = await loadModel(); |
| 43 | +async function learnedTissue(vol, want) { // vol: {data,dims,affine} -> Uint8 native label 0=bg,1=CSF,2=GM,3=WM |
| 44 | + P(20, 'Loading AI model…'); const sess = await loadModel(want); |
27 | 45 | const { M, T } = affToMT(vol.affine); |
28 | | - P(45, 'Conforming to model grid…'); |
| 46 | + P(40, 'Conforming to model grid…'); |
29 | 47 | const { x, Tat } = conformVol({ data: vol.data, dims: vol.dims, M, T }); |
30 | | - P(60, 'Running AI brain extraction…'); |
| 48 | + P(55, 'Running AI tissue segmentation…'); // ~40s single-threaded WASM at 192³ |
31 | 49 | const out = await sess.run({ t1: new _ort.Tensor('float32', x, [1, 1, BET_SHAPE[0], BET_SHAPE[1], BET_SHAPE[2]]) }); |
32 | | - const logit = (out.logit || out[Object.keys(out)[0]]).data; |
33 | | - const prob = new Float32Array(logit.length); |
34 | | - for (let i = 0; i < logit.length; i++) prob[i] = 1 / (1 + Math.exp(-logit[i])); |
35 | | - P(80, 'Warping mask to native space…'); |
36 | | - const np = probToNative(prob, Tat, { data: vol.data, dims: vol.dims, M, T }); |
37 | | - const mask = new Uint8Array(np.length); for (let i = 0; i < np.length; i++) mask[i] = np[i] >= 0.5 ? 1 : 0; |
38 | | - return mask; |
| 50 | + const o = out.logit || out[Object.keys(out)[0]]; |
| 51 | + const d = o.data, C = o.dims[1], V = BET_SHAPE[0] * BET_SHAPE[1] * BET_SHAPE[2]; |
| 52 | + P(78, 'Warping to native space…'); |
| 53 | + const label = new Uint8Array(V); |
| 54 | + if (C === 1) { for (let v = 0; v < V; v++) label[v] = d[v] > 0 ? 1 : 0; } // binary fallback (old model) |
| 55 | + else { for (let v = 0; v < V; v++) { let best = d[v], bi = 0; for (let c = 1; c < C; c++) { const val = d[c*V + v]; if (val > best) { best = val; bi = c; } } label[v] = bi; } } |
| 56 | + return labelToNative(label, Tat, { data: vol.data, dims: vol.dims, M, T }); // nearest -> native classes |
| 57 | +} |
| 58 | + |
| 59 | +// ---- full-head 6-tissue model (bg/CSF/GM/WM/skull/scalp) — its own grid + session ---- |
| 60 | +const HEAD_SHAPE = [160, 192, 160], HEAD_VOX = 1.4; |
| 61 | +let _headBuf = null, _headSess = null, _headBackend = null; |
| 62 | +async function loadHeadModel(want) { |
| 63 | + want = (want === 'webgpu') ? 'webgpu' : 'wasm'; |
| 64 | + if (_headSess && _headBackend === want) return _headSess; |
| 65 | + if (!_ort) { |
| 66 | + _ort = await import(ORT_CDN + 'ort.webgpu.min.mjs'); _ort.env.wasm.wasmPaths = ORT_CDN; |
| 67 | + _ort.env.wasm.numThreads = (self.crossOriginIsolated && navigator.hardwareConcurrency) ? Math.min(8, navigator.hardwareConcurrency) : 1; |
| 68 | + } |
| 69 | + if (_headSess) { try { await _headSess.release?.(); } catch (e) {} _headSess = null; } |
| 70 | + if (!_headBuf) { const url = new URL('model/bet_head6.onnx', import.meta.url); _headBuf = await (await fetch(url)).arrayBuffer(); } |
| 71 | + if (want === 'webgpu') { |
| 72 | + try { _headSess = await _ort.InferenceSession.create(_headBuf, { executionProviders: ['webgpu'] }); _headBackend = 'webgpu'; _backend = 'webgpu'; return _headSess; } |
| 73 | + catch (e) { _backend = 'wasm (WebGPU unavailable)'; } |
| 74 | + } |
| 75 | + _headSess = await _ort.InferenceSession.create(_headBuf, { executionProviders: ['wasm'] }); _headBackend = 'wasm'; |
| 76 | + if (want === 'wasm') _backend = _ort.env.wasm.numThreads > 1 ? `wasm×${_ort.env.wasm.numThreads}` : 'wasm'; |
| 77 | + return _headSess; |
| 78 | +} |
| 79 | +// -> Uint8 native 6-class label: 0=bg,1=CSF,2=GM,3=WM,4=skull,5=scalp |
| 80 | +async function learnedHead(vol, want) { |
| 81 | + P(18, 'Loading full-head model…'); const sess = await loadHeadModel(want); |
| 82 | + const { M, T } = affToMT(vol.affine); |
| 83 | + P(38, 'Conforming to head grid…'); |
| 84 | + const { x, Tat } = conformVol({ data: vol.data, dims: vol.dims, M, T }, HEAD_SHAPE, HEAD_VOX); |
| 85 | + P(52, 'Running 6-tissue segmentation…'); |
| 86 | + const out = await sess.run({ t1: new _ort.Tensor('float32', x, [1, 1, HEAD_SHAPE[0], HEAD_SHAPE[1], HEAD_SHAPE[2]]) }); |
| 87 | + const o = out.logit || out[Object.keys(out)[0]]; |
| 88 | + const d = o.data, C = o.dims[1], V = HEAD_SHAPE[0] * HEAD_SHAPE[1] * HEAD_SHAPE[2]; |
| 89 | + P(80, 'Warping to native space…'); |
| 90 | + const label = new Uint8Array(V); |
| 91 | + for (let v = 0; v < V; v++) { let best = d[v], bi = 0; for (let c = 1; c < C; c++) { const val = d[c*V + v]; if (val > best) { best = val; bi = c; } } label[v] = bi; } |
| 92 | + return labelToNative(label, Tat, { data: vol.data, dims: vol.dims, M, T }, HEAD_SHAPE, HEAD_VOX); |
| 93 | +} |
| 94 | + |
| 95 | +// grow GM(2)/WM(3) outward into bright partial-volume voxels (> the CSF<->GM intensity valley), up to |
| 96 | +// `iters` voxels — recovers the pial rim the >50%-GM teacher trims. Stops at dark CSF/skull. |
| 97 | +function pialPush(nat, data, dims, iters) { |
| 98 | + const [nx, ny, nz] = dims, N = nx * ny * nz; let gs = 0, gn = 0, cs = 0, cn = 0; |
| 99 | + for (let i = 0; i < N; i++) { const c = nat[i]; if (c === 2) { gs += data[i]; gn++; } else if (c === 1) { cs += data[i]; cn++; } } |
| 100 | + if (!gn || !cn) return nat; const thr = 0.5 * (gs / gn + cs / cn); const out = Uint8Array.from(nat), nxy = nx * ny; |
| 101 | + for (let it = 0; it < iters; it++) { const add = []; |
| 102 | + for (let z = 1; z < nz - 1; z++) for (let y = 1; y < ny - 1; y++) { const b = nx * (y + ny * z); |
| 103 | + for (let x = 1; x < nx - 1; x++) { const i = b + x; if (out[i] >= 2 || data[i] <= thr) continue; |
| 104 | + if (out[i-1] >= 2 || out[i+1] >= 2 || out[i-nx] >= 2 || out[i+nx] >= 2 || out[i-nxy] >= 2 || out[i+nxy] >= 2) add.push(i); } } |
| 105 | + if (!add.length) break; for (let k = 0; k < add.length; k++) out[add[k]] = 2; } |
| 106 | + return out; |
| 107 | +} |
| 108 | +// derive brain mask + CSF/GM/WM from a native label honoring the BET controls (pial reach, mask type), then post |
| 109 | +function deriveAndPost(vol, natRaw, opts) { |
| 110 | + const [nx, ny, nz] = vol.dims, N = nx * ny * nz, data = vol.data; |
| 111 | + P(90, 'Applying controls…'); |
| 112 | + const pial = Math.max(0, Math.min(3, (opts.pial | 0))); |
| 113 | + const nat = pial > 0 ? pialPush(natRaw, data, vol.dims, pial) : natRaw; |
| 114 | + const whole = opts.maskType === 'whole'; |
| 115 | + let mask = new Uint8Array(N); for (let i = 0; i < N; i++) mask[i] = (whole ? nat[i] > 0 : nat[i] >= 2) ? 1 : 0; |
| 116 | + mask = largestComponent(mask, vol.dims); mask = fillHoles3D(mask, vol.dims); |
| 117 | + const cortex = new Float32Array(N), wm = new Float32Array(N), gm = new Float32Array(N); |
| 118 | + let cs = 0, cn = 0, gs = 0, gn = 0, ws = 0, wn = 0, vox = 0; |
| 119 | + for (let i = 0; i < N; i++) { if (!mask[i]) continue; vox++; const c = nat[i], val = data[i]; |
| 120 | + if (c === 1) { cs += val; cn++; } else if (c === 2) { gm[i] = val; cortex[i] = val; gs += val; gn++; } else if (c === 3) { wm[i] = val; cortex[i] = val; ws += val; wn++; } } |
| 121 | + const centroids = [cn ? cs / cn : 0, gn ? gs / gn : 0, wn ? ws / wn : 0]; |
| 122 | + const note = `AI tissue U-Net (1mm 3-class, ${(100 * vox / N).toFixed(1)}%${pial ? ', pial+' + pial : ''}, ${whole ? 'whole' : 'cortical'}) · ${_backend}`; |
| 123 | + const m8 = (mask instanceof Uint8Array) ? mask : Uint8Array.from(mask); |
| 124 | + P(100, 'Done'); |
| 125 | + self.postMessage({ type: 'done', mask: m8, cortex, wm, gm, voxels: vox, total: N, centroids, note }, [m8.buffer, cortex.buffer, wm.buffer, gm.buffer]); |
39 | 126 | } |
40 | 127 |
|
41 | 128 | self.onmessage = async (e) => { |
42 | | - const m = e.data; if (m.cmd !== 'extract') return; |
| 129 | + const m = e.data; |
| 130 | + if (m.cmd === 'reprocess') { // re-apply BET controls to the cached label — no 40s re-inference |
| 131 | + if (_last) { try { deriveAndPost(_last.vol, _last.nat, m.opts || {}); } catch (err) { self.postMessage({ type: 'error', msg: (err && err.message) || String(err) }); } } |
| 132 | + return; |
| 133 | + } |
| 134 | + if (m.cmd !== 'extract') return; |
43 | 135 | try { |
44 | 136 | const { vol, mni, prior, opts = {} } = m; |
45 | 137 | const [nx, ny, nz] = vol.dims, N = nx * ny * nz; |
46 | 138 | let mask, note; |
| 139 | + // ---- full-head 6-tissue layers: post the native label volume for the peelable head render ---- |
| 140 | + if (opts.method === 'headlayers') { |
| 141 | + const nat = await learnedHead(vol, opts.backend); // 0bg 1CSF 2GM 3WM 4skull 5scalp |
| 142 | + const lab = (nat instanceof Uint8Array) ? nat : Uint8Array.from(nat); |
| 143 | + P(100, 'Done'); |
| 144 | + self.postMessage({ type: 'headdone', label: lab, dims: vol.dims, affine: vol.affine, note: `Full-head 6-tissue U-Net · ${_backend}` }, [lab.buffer]); |
| 145 | + return; |
| 146 | + } |
| 147 | + // ---- learned engine: the 3-class U-Net gives brain + CSF/GM/WM directly (no k-means heuristic) ---- |
47 | 148 | if (opts.method === 'learned') { |
48 | | - mask = await learnedMask(vol); |
49 | | - let v0 = 0; for (let i = 0; i < N; i++) v0 += mask[i]; |
50 | | - note = `AI U-Net (distilled, ${(100 * v0 / N).toFixed(1)}%)`; |
51 | | - } else { |
52 | | - P(6, 'Thresholding & deep-core…'); |
53 | | - const M0 = extractBrain(vol, { method: 'deepcore' }); |
54 | | - mask = M0.mask; note = `deep-core ${(100 * M0.voxels / N).toFixed(1)}%`; |
55 | | - if (opts.refined && mni && prior) { |
56 | | - P(42, 'Registering brain → MNI…'); |
57 | | - const brainVol = { data: applyMask(vol.data, mask), dims: vol.dims, affine: vol.affine }; |
58 | | - const reg = registerAffine(brainVol, mni, {}); |
59 | | - P(70, `Warping OASIS prior back (NMI ${reg.nmi.toFixed(3)})…`); |
60 | | - const S2M = invert4x4(reg.matrix), PT = 255 * (prior.pt ?? 0.20); |
61 | | - const mniMask = new Int32Array(prior.data.length); for (let i = 0; i < mniMask.length; i++) mniMask[i] = prior.data[i] >= PT ? 1 : 0; |
62 | | - let np = resampleLabels({ data: mniMask, dims: prior.dims, affine: prior.affine }, { dims: vol.dims, affine: vol.affine }, S2M); |
63 | | - let pr = new Uint8Array(N); for (let i = 0; i < N; i++) pr[i] = np[i] ? 1 : 0; |
64 | | - const dil = Math.max(1, Math.round((opts.dilMM ?? 6) / Math.min(vol.spacing[0], vol.spacing[1], vol.spacing[2]))); |
65 | | - pr = dilate3D(pr, vol.dims, dil); |
66 | | - const g = new Uint8Array(N); for (let i = 0; i < N; i++) g[i] = (mask[i] && pr[i]) ? 1 : 0; |
67 | | - mask = g; note += ` → OASIS-gated (n=${prior.n}, P≥${prior.pt ?? 0.20})`; |
68 | | - } |
| 149 | + const nat = await learnedTissue(vol, opts.backend); // native label 0=bg,1=CSF,2=GM,3=WM |
| 150 | + _last = { vol: { data: vol.data, dims: vol.dims }, nat }; // cache so control tweaks re-derive instantly |
| 151 | + deriveAndPost(vol, nat, opts); // honors opts.pial / opts.maskType |
| 152 | + return; |
| 153 | + } |
| 154 | + // ---- geometric engines (deep-core / OASIS-gated refined) ---- |
| 155 | + P(6, 'Thresholding & deep-core…'); |
| 156 | + const M0 = extractBrain(vol, { method: 'deepcore' }); |
| 157 | + mask = M0.mask; note = `deep-core ${(100 * M0.voxels / N).toFixed(1)}%`; |
| 158 | + if (opts.refined && mni && prior) { |
| 159 | + P(42, 'Registering brain → MNI…'); |
| 160 | + const brainVol = { data: applyMask(vol.data, mask), dims: vol.dims, affine: vol.affine }; |
| 161 | + const reg = registerAffine(brainVol, mni, {}); |
| 162 | + P(70, `Warping OASIS prior back (NMI ${reg.nmi.toFixed(3)})…`); |
| 163 | + const S2M = invert4x4(reg.matrix), PT = 255 * (prior.pt ?? 0.20); |
| 164 | + const mniMask = new Int32Array(prior.data.length); for (let i = 0; i < mniMask.length; i++) mniMask[i] = prior.data[i] >= PT ? 1 : 0; |
| 165 | + let np = resampleLabels({ data: mniMask, dims: prior.dims, affine: prior.affine }, { dims: vol.dims, affine: vol.affine }, S2M); |
| 166 | + let pr = new Uint8Array(N); for (let i = 0; i < N; i++) pr[i] = np[i] ? 1 : 0; |
| 167 | + const dil = Math.max(1, Math.round((opts.dilMM ?? 6) / Math.min(vol.spacing[0], vol.spacing[1], vol.spacing[2]))); |
| 168 | + pr = dilate3D(pr, vol.dims, dil); |
| 169 | + const g = new Uint8Array(N); for (let i = 0; i < N; i++) g[i] = (mask[i] && pr[i]) ? 1 : 0; |
| 170 | + mask = g; note += ` → OASIS-gated (n=${prior.n}, P≥${prior.pt ?? 0.20})`; |
69 | 171 | } |
70 | 172 | P(86, 'Cleaning mask…'); |
71 | | - mask = largestComponent(mask, vol.dims); mask = fillHoles3D(mask, vol.dims); mask = smoothMask3D(mask, vol.dims, opts.method === 'learned' ? 1 : (opts.refined ? 2 : 1)); |
| 173 | + mask = largestComponent(mask, vol.dims); mask = fillHoles3D(mask, vol.dims); mask = smoothMask3D(mask, vol.dims, opts.refined ? 2 : 1); |
72 | 174 | P(90, 'Segmenting tissue (CSF/GM/WM)…'); |
73 | 175 | const seg = segmentTissue(vol.data, mask, { invert: !!opts.t2 }); |
74 | 176 | const cortex = new Float32Array(N), wm = new Float32Array(N), gm = new Float32Array(N); |
|
0 commit comments