1// A from-scratch isochromat Bloch simulator for pulseq sequences.
2//
3// The magnetisation of every spin is evolved through the sequence in the
4// rotating frame. The timeline is split at every gradient vertex, every RF
5// sample, and every ADC sample (the union is `boundaries`), so within each
6// segment the gradients are linear and the RF is ~one raster long. Two regimes:
7//
8// * Free precession (no RF in the segment): the effective field is purely
9// longitudinal, so each spin's transverse magnetisation just rotates about
10// z by the gradient phase dφ = 2π·(g·r)·dt + Δw·dt and relaxes. Exact for
11// a linear gradient, so these segments can be as long as a whole readout
12// dwell or delay — this is what makes the simulation fast.
13// * Excitation (RF present): the full 3-D effective field
14// Ω = (2π·B1·cosθ, 2π·B1·sinθ, 2π·g·r + Δw) rotates each spin (Rodrigues).
15//
16// Signal at each ADC sample = Σ_j ρ_j·(Mx+iMy), demodulated by the receiver
17// phase (ADC phase + frequency offset). Units: pulseq gradients are Hz/m and RF
18// amplitude Hz, positions m, so g·r and B1 are already in Hz.
19import type { ParsedSeq } from '../seq/types.ts'
20import { reconstruct } from '../seq/reconstruct.ts'
21import type { Series } from '../seq/reconstruct.ts'
22import type { Phantom } from '../phantom/phantomTypes.ts'
24export interface RawSignal {
25 /** Real part of the acquired signal, one entry per ADC sample (all readouts concatenated). */
26 re: Float64Array
27 im: Float64Array
28 /** Number of ADC events (readout lines). */
29 numReadouts: number
30 /** Samples in each readout, length numReadouts. */
31 samplesPerReadout: Int32Array
32 /** Prefix offsets into re/im, length numReadouts+1. */
33 offsets: Int32Array
34 maxSamplesPerReadout: number
35}
37export interface SimProgress {
38 fraction: number
39 segment: number
40 numSegments: number
41 samplesDone: number
42 numSamples: number
43}
45export interface SimOptions {
46 onProgress?: (p: SimProgress) => void
47 /** How often (in segments) to invoke onProgress. */
48 progressEvery?: number
49}
51const TWO_PI = 2 * Math.PI
53/** Linear interpolation of a (t,v) polyline at ascending query times (moving cursor, O(n)). */
54function sampleAtTimes(series: Series, times: Float64Array): Float64Array {
55 const { t, v } = series
56 const n = t.length
57 const out = new Float64Array(times.length)
58 if (n === 0) return out
59 let p = 0
60 for (let i = 0; i < times.length; i++) {
61 const time = times[i]
62 if (time <= t[0]) {
63 out[i] = v[0]
64 continue
65 }
66 if (time >= t[n - 1]) {
67 out[i] = v[n - 1]
68 continue
69 }
70 while (p < n - 1 && t[p + 1] < time) p++
71 const ta = t[p]
72 const tb = t[p + 1]
73 if (tb === ta) {
74 out[i] = v[p + 1]
75 } else {
76 const f = (time - ta) / (tb - ta)
77 out[i] = v[p] + f * (v[p + 1] - v[p])
78 }
79 }
80 return out
81}
83/** One ADC sample's recording target: which readout, which column, and receiver phase. */
84interface AdcSampleRec {
85 segment: number // record after the segment ending at this sample's time
86 readout: number
87 col: number
88 cosPhi: number // demod: multiply signal by exp(-i φ_rx) = (cosPhi - i sinPhi)... stored as cos/sin of φ_rx
89 sinPhi: number
90}
92interface Schedule {
93 boundaries: Float64Array
94 midpoints: Float64Array
95 dt: Float64Array
96 gxMid: Float64Array
97 gyMid: Float64Array
98 gzMid: Float64Array
99 rfMagMid: Float64Array
100 rfActive: Uint8Array
101 records: AdcSampleRec[]
102 numReadouts: number
103 samplesPerReadout: Int32Array
104 offsets: Int32Array
105 numSamples: number
106 rfPhase: Series
107 rfMag: Series
108}
110function buildSchedule(seq: ParsedSeq): Schedule {
111 const rec = reconstruct(seq)
113 // ADC sample schedule, straight from the parsed events (so we keep phase/freq).
114 const adcSampleTimes: number[] = []
115 const adcSampleReadout: number[] = []
116 const adcSampleCol: number[] = []
117 const adcSamplePhi: number[] = []
118 const samplesPerReadout: number[] = []
119 let readout = 0
120 let t0 = 0
121 for (let bi = 0; bi < seq.blocks.length; bi++) {
122 const block = seq.blocks[bi]
123 if (block.adcId !== 0) {
124 const adc = seq.adcs.get(block.adcId)
125 if (adc && adc.num > 0) {
126 const start = t0 + adc.delay
127 samplesPerReadout.push(adc.num)
128 for (let i = 0; i < adc.num; i++) {
129 const tk = start + adc.dwell * (i + 0.5)
130 // Receiver phase: constant ADC phase offset + frequency offset ramp.
131 const phi = adc.phase + TWO_PI * adc.freq * (tk - start)
132 adcSampleTimes.push(tk)
133 adcSampleReadout.push(readout)
134 adcSampleCol.push(i)
135 adcSamplePhi.push(phi)
136 }
137 readout++
138 }
139 }
140 t0 += block.duration
141 }
142 const numReadouts = readout
143 const numSamples = adcSampleTimes.length
145 // Union of all event boundaries (gradient vertices, RF samples, ADC samples, endpoints).
146 const times = new Set<number>()
147 times.add(0)
148 times.add(rec.duration)
149 for (const s of [rec.gx, rec.gy, rec.gz, rec.rfMag]) {
150 for (let i = 0; i < s.t.length; i++) if (!Number.isNaN(s.t[i])) times.add(s.t[i])
151 }
152 for (const s of rec.rfSpans) {
153 times.add(s.start)
154 times.add(s.end)
155 }
156 for (const tk of adcSampleTimes) times.add(tk)
157 const boundaries = Float64Array.from(times)
158 boundaries.sort()
160 const numSegments = boundaries.length - 1
161 const midpoints = new Float64Array(numSegments)
162 const dt = new Float64Array(numSegments)
163 for (let i = 0; i < numSegments; i++) {
164 midpoints[i] = 0.5 * (boundaries[i] + boundaries[i + 1])
165 dt[i] = boundaries[i + 1] - boundaries[i]
166 }
168 const gxMid = sampleAtTimes(rec.gx, midpoints)
169 const gyMid = sampleAtTimes(rec.gy, midpoints)
170 const gzMid = sampleAtTimes(rec.gz, midpoints)
171 const rfMagMid = sampleAtTimes(rec.rfMag, midpoints)
173 // Mark segments that fall inside an RF pulse span.
174 const rfActive = new Uint8Array(numSegments)
175 for (const span of rec.rfSpans) {
176 for (let i = 0; i < numSegments; i++) {
177 if (midpoints[i] > span.start && midpoints[i] < span.end) rfActive[i] = 1
178 }
179 }
181 // Map each ADC sample time to the segment it ends, via a value->index lookup.
182 const indexOf = new Map<number, number>()
183 for (let i = 0; i < boundaries.length; i++) indexOf.set(boundaries[i], i)
184 const records: AdcSampleRec[] = []
185 for (let k = 0; k < numSamples; k++) {
186 const bIndex = indexOf.get(adcSampleTimes[k])
187 if (bIndex === undefined || bIndex === 0) continue
188 const phi = adcSamplePhi[k]
189 records.push({
190 segment: bIndex - 1,
191 readout: adcSampleReadout[k],
192 col: adcSampleCol[k],
193 cosPhi: Math.cos(phi),
194 sinPhi: Math.sin(phi),
195 })
196 }
197 // Group records by the segment they fire after.
198 records.sort((a, b) => a.segment - b.segment)
200 const spr = Int32Array.from(samplesPerReadout)
201 const offsets = new Int32Array(numReadouts + 1)
202 for (let r = 0; r < numReadouts; r++) offsets[r + 1] = offsets[r] + spr[r]
204 return {
205 boundaries,
206 midpoints,
207 dt,
208 gxMid,
209 gyMid,
210 gzMid,
211 rfMagMid,
212 rfActive,
213 records,
214 numReadouts,
215 samplesPerReadout: spr,
216 offsets,
217 numSamples,
218 rfPhase: rec.rfPhase,
219 rfMag: rec.rfMag,
220 }
221}
223export function simulate(seq: ParsedSeq, phantom: Phantom, opts: SimOptions = {}): RawSignal {
224 const sched = buildSchedule(seq)
225 const { boundaries, midpoints, dt, gxMid, gyMid, gzMid, rfMagMid, rfActive, records } = sched
226 const numSegments = boundaries.length - 1
227 const ns = phantom.ns
228 const { x, y, z, rho, t1, t2, dw } = phantom
230 // Magnetisation, initialised at thermal equilibrium (Mz = ρ).
231 const Mx = new Float64Array(ns)
232 const My = new Float64Array(ns)
233 const Mz = new Float64Array(ns)
234 for (let j = 0; j < ns; j++) Mz[j] = rho[j]
236 const re = new Float64Array(sched.numSamples)
237 const im = new Float64Array(sched.numSamples)
239 // Fast path when relaxation is uniform across the phantom (our built-ins are).
240 let uniformRelax = true
241 for (let j = 1; j < ns; j++) {
242 if (t1[j] !== t1[0] || t2[j] !== t2[0]) {
243 uniformRelax = false
244 break
245 }
246 }
247 const r1u = ns > 0 ? 1 / t1[0] : 0
248 const r2u = ns > 0 ? 1 / t2[0] : 0
250 const progressEvery = opts.progressEvery ?? 256
251 let recPtr = 0
253 for (let i = 0; i < numSegments; i++) {
254 const dti = dt[i]
255 if (dti > 0) {
256 if (rfActive[i]) {
257 // --- Excitation: full 3-D rotation about the effective field ---
258 const b1 = rfMagMid[i]
259 const phase = sampleSeriesScalar(sched.rfPhase, midpoints[i])
260 const ph = Number.isNaN(phase) ? 0 : phase
261 const w1 = TWO_PI * b1
262 const wx = w1 * Math.cos(ph)
263 const wy = w1 * Math.sin(ph)
264 const gx = gxMid[i]
265 const gy = gyMid[i]
266 const gz = gzMid[i]
267 const e1 = uniformRelax ? Math.exp(-dti * r1u) : 0
268 const e2 = uniformRelax ? Math.exp(-dti * r2u) : 0
269 for (let j = 0; j < ns; j++) {
270 const wz = TWO_PI * (gx * x[j] + gy * y[j] + gz * z[j]) + dw[j]
271 const wmag = Math.sqrt(wx * wx + wy * wy + wz * wz)
272 let mx = Mx[j]
273 let my = My[j]
274 let mz = Mz[j]
275 if (wmag > 0) {
276 // Rotate by θ = -wmag·dt about n = (wx,wy,wz)/wmag (sign matches free precession).
277 const theta = -wmag * dti
278 const c = Math.cos(theta)
279 const s = Math.sin(theta)
280 const inv = 1 / wmag
281 const nx = wx * inv
282 const ny = wy * inv
283 const nz = wz * inv
284 const dot = nx * mx + ny * my + nz * mz
285 // Rodrigues: m' = m c + (n×m) s + n (n·m)(1-c)
286 const crx = ny * mz - nz * my
287 const cry = nz * mx - nx * mz
288 const crz = nx * my - ny * mx
289 const k = dot * (1 - c)
290 mx = mx * c + crx * s + nx * k
291 my = my * c + cry * s + ny * k
292 mz = mz * c + crz * s + nz * k
293 }
294 const e2j = uniformRelax ? e2 : Math.exp(-dti / t2[j])
295 const e1j = uniformRelax ? e1 : Math.exp(-dti / t1[j])
296 Mx[j] = mx * e2j
297 My[j] = my * e2j
298 Mz[j] = mz * e1j + rho[j] * (1 - e1j)
299 }
300 } else {
301 // --- Free precession: rotate about z by the gradient phase, then relax ---
302 const ax = TWO_PI * gxMid[i] * dti
303 const ay = TWO_PI * gyMid[i] * dti
304 const az = TWO_PI * gzMid[i] * dti
305 const e1 = uniformRelax ? Math.exp(-dti * r1u) : 0
306 const e2 = uniformRelax ? Math.exp(-dti * r2u) : 0
307 for (let j = 0; j < ns; j++) {
308 const dphi = ax * x[j] + ay * y[j] + az * z[j] + dw[j] * dti
309 const c = Math.cos(dphi)
310 const s = Math.sin(dphi)
311 const mx = Mx[j]
312 const my = My[j]
313 const e2j = uniformRelax ? e2 : Math.exp(-dti / t2[j])
314 const e1j = uniformRelax ? e1 : Math.exp(-dti / t1[j])
315 Mx[j] = (mx * c + my * s) * e2j
316 My[j] = (-mx * s + my * c) * e2j
317 Mz[j] = Mz[j] * e1j + rho[j] * (1 - e1j)
318 }
319 }
320 }
322 // Record any ADC samples that fire at the end of this segment.
323 while (recPtr < records.length && records[recPtr].segment === i) {
324 const r = records[recPtr]
325 let sre = 0
326 let sim = 0
327 for (let j = 0; j < ns; j++) {
328 sre += rho[j] * Mx[j]
329 sim += rho[j] * My[j]
330 }
331 // recorded = (sre + i·sim) · exp(-i φ_rx)
332 const idx = sched.offsets[r.readout] + r.col
333 re[idx] = sre * r.cosPhi + sim * r.sinPhi
334 im[idx] = -sre * r.sinPhi + sim * r.cosPhi
335 recPtr++
336 }
338 if (opts.onProgress && (i % progressEvery === 0 || i === numSegments - 1)) {
339 opts.onProgress({
340 fraction: numSegments > 0 ? (i + 1) / numSegments : 1,
341 segment: i + 1,
342 numSegments,
343 samplesDone: recPtr,
344 numSamples: sched.numSamples,
345 })
346 }
347 }
349 let maxSamplesPerReadout = 0
350 for (let r = 0; r < sched.numReadouts; r++)
351 maxSamplesPerReadout = Math.max(maxSamplesPerReadout, sched.samplesPerReadout[r])
353 return {
354 re,
355 im,
356 numReadouts: sched.numReadouts,
357 samplesPerReadout: sched.samplesPerReadout,
358 offsets: sched.offsets,
359 maxSamplesPerReadout,
360 }
361}
363/** Point sample of a (t,v) polyline (binary search); used only at RF midpoints. */
364function sampleSeriesScalar(series: Series, time: number): number {
365 const { t, v } = series
366 const n = t.length
367 if (n === 0) return NaN
368 if (time <= t[0]) return v[0]
369 if (time >= t[n - 1]) return v[n - 1]
370 let lo = 0
371 let hi = n - 1
372 while (hi - lo > 1) {
373 const mid = (lo + hi) >> 1
374 if (t[mid] <= time) lo = mid
375 else hi = mid
376 }
377 const ta = t[lo]
378 const tb = t[hi]
379 const va = v[lo]
380 const vb = v[hi]
381 if (Number.isNaN(va) || Number.isNaN(vb)) return Number.isNaN(va) ? vb : va
382 if (tb === ta) return vb
383 return va + ((time - ta) / (tb - ta)) * (vb - va)
384}