1// Headless smoke test for the Bloch simulator. Reads a phantom (.phantom via
2// h5wasm) and a golden .seq, runs simulate(), and checks basic physics:
3// - FID: right after a 90° pulse all spins are in phase, so |S| at the first
4// ADC sample ≈ Σρ (≈ number of spins), decaying by T2 thereafter.
5// - GRE: the DC point of k-space (gradient-refocused echo centre) should be
6// the brightest sample in its readout.
7// Run: npm run sim-test
8import * as h5 from 'h5wasm/node'
9import fs from 'node:fs'
10import path from 'node:path'
11import { fileURLToPath } from 'node:url'
12import { parseSeq } from '../src/seq/parseSeq.ts'
13import { simulate } from '../src/sim/simulate.ts'
15const __dirname = path.dirname(fileURLToPath(import.meta.url))
16const root = path.resolve(__dirname, '..')
18function loadPhantomNode(file) {
19 const bytes = fs.readFileSync(file)
20 // h5wasm/node writes to the real cwd; use a temp name and unlink after.
21 const tmp = `.sim-test-${path.basename(file)}.tmp`
22 h5.FS.writeFile(tmp, new Uint8Array(bytes))
23 const f = new h5.File(tmp, 'r')
24 const g = (p) => Float32Array.from(f.get(p).value)
25 const ph = {
26 name: path.basename(file),
27 x: g('position/x'),
28 y: g('position/y'),
29 z: g('position/z'),
30 rho: g('contrast/ρ'),
31 t1: g('contrast/T1'),
32 t2: g('contrast/T2'),
33 t2s: g('contrast/T2s'),
34 dw: g('contrast/Δw'),
35 }
36 ph.ns = ph.x.length
37 f.close()
38 h5.FS.unlink(tmp)
39 return ph
40}
42function mag(sig, i) {
43 return Math.hypot(sig.re[i], sig.im[i])
44}
46let failures = 0
47function check(name, cond, detail) {
48 if (cond) {
49 console.log(` ✓ ${name}${detail ? ' — ' + detail : ''}`)
50 } else {
51 console.log(` ✗ ${name}${detail ? ' — ' + detail : ''}`)
52 failures++
53 }
54}
56async function main() {
57 await h5.ready
58 const cube = loadPhantomNode(path.join(root, 'src/phantom/data/cube.phantom'))
59 console.log(`Phantom: cube, ${cube.ns} spins`)
61 // --- FID ---
62 console.log('\nFID (test-data/fid.seq):')
63 const fidSeq = parseSeq(fs.readFileSync(path.join(root, 'test-data/fid.seq'), 'utf8'))
64 const t0 = Date.now()
65 const fid = simulate(fidSeq, cube)
66 console.log(` simulated in ${Date.now() - t0} ms; ${fid.numReadouts} readouts, ${fid.re.length} samples`)
67 const m0 = mag(fid, 0)
68 const sumRho = cube.rho.reduce((a, b) => a + b, 0)
69 check('all samples finite', fid.re.every((v) => Number.isFinite(v)) && fid.im.every((v) => Number.isFinite(v)))
70 // Right after a 90° pulse all spins are in phase; |S(0)| = Σρ, minus T2 decay
71 // over the ~20 ms of dead time before the first sample. So it should be a
72 // large fraction of Σρ but not exceed it.
73 check('|S(0)| ~ Σρ (in phase)', m0 > 0.7 * sumRho && m0 <= sumRho * 1.001, `|S(0)|=${m0.toFixed(1)}, Σρ=${sumRho.toFixed(0)}`)
74 // Direct T2 validation: within a readout, |S| must decay as exp(-t/T2).
75 const dwell = fidSeq.adcs.get(1).dwell
76 const T2 = cube.t2[0]
77 const K = 400
78 const ratio = mag(fid, K) / m0
79 const expected = Math.exp(-(K * dwell) / T2)
80 check('FID decays at rate 1/T2', Math.abs(ratio - expected) / expected < 0.02, `|S(${K})|/|S(0)|=${ratio.toFixed(4)}, exp(-t/T2)=${expected.toFixed(4)}`)
82 // --- GRE ---
83 console.log('\nGRE (test-data/gre.seq):')
84 const greSeq = parseSeq(fs.readFileSync(path.join(root, 'test-data/gre.seq'), 'utf8'))
85 const t1 = Date.now()
86 const gre = simulate(greSeq, cube, {
87 onProgress: (p) => {
88 if (p.segment === 1 || p.fraction === 1 || p.segment % 20000 === 0)
89 process.stdout.write(`\r progress ${(p.fraction * 100).toFixed(1)}% (seg ${p.segment}/${p.numSegments}, ${p.samplesDone}/${p.numSamples} samples) `)
90 },
91 })
92 process.stdout.write('\n')
93 console.log(` simulated in ${((Date.now() - t1) / 1000).toFixed(1)} s; ${gre.numReadouts} readouts, ${gre.maxSamplesPerReadout} samples/readout`)
94 check('all samples finite', gre.re.every((v) => Number.isFinite(v)) && gre.im.every((v) => Number.isFinite(v)))
95 // For the middle readout, the brightest sample should be near the readout centre (k=0 echo).
96 const midR = Math.floor(gre.numReadouts / 2)
97 const o = gre.offsets[midR]
98 const n = gre.samplesPerReadout[midR]
99 let peak = -1
100 let peakIdx = -1
101 for (let i = 0; i < n; i++) {
102 const m = mag(gre, o + i)
103 if (m > peak) {
104 peak = m
105 peakIdx = i
106 }
107 }
108 check('GRE echo peaks near readout centre', Math.abs(peakIdx - n / 2) < n * 0.2, `peak at col ${peakIdx}/${n}`)
110 console.log(`\n${failures === 0 ? 'All checks passed.' : failures + ' check(s) FAILED.'}`)
111 process.exit(failures === 0 ? 0 : 1)
112}
114main()