1// Turn a ParsedSeq into plottable timelines, mirroring pulseq's
2// Sequence.waveforms_and_times conventions:
3// - trapezoids -> 4 vertices (3 when flat == 0), skipped when empty
4// - arbitrary gradients on the centers raster -> edge points from
5// first/last plus the raster-center samples
6// - extended trapezoids (time shape on raster edges) -> vertices as stored
7// - RF samples at (i+0.5)*rfRaster or the time shape; phase includes the
8// phase/frequency offsets
9// Series are piecewise-linear polylines; NaN values break the line.
10import type { ParsedSeq } from './types.ts'
12export interface Series {
13 t: Float64Array
14 v: Float64Array
15}
17export interface BlockSpan {
18 index: number
19 id: number
20 start: number
21 duration: number
22}
24export interface AdcSpan {
25 blockIndex: number
26 start: number
27 end: number
28 num: number
29 dwell: number
30}
32export interface RfPulseSpan {
33 blockIndex: number
34 start: number
35 end: number
36 /** Absolute time of the pulse center (NaN for v1.4 files) */
37 center: number
38 use: string
39}
41export interface Reconstructed {
42 duration: number
43 blockSpans: BlockSpan[]
44 gx: Series
45 gy: Series
46 gz: Series
47 /** |B1| in Hz */
48 rfMag: Series
49 /** rad, NaN between pulses */
50 rfPhase: Series
51 adcSpans: AdcSpan[]
52 rfSpans: RfPulseSpan[]
53 stats: {
54 maxGrad: { gx: number; gy: number; gz: number }
55 maxSlew: { gx: number; gy: number; gz: number }
56 rfCount: number
57 adcCount: number
58 adcSamples: number
59 }
60}
62const EPS = 1e-9
64class SeriesBuilder {
65 t: number[] = []
66 v: number[] = []
67 push(t: number, v: number) {
68 this.t.push(t)
69 this.v.push(v)
70 }
71 breakLine() {
72 if (this.v.length > 0 && !Number.isNaN(this.v[this.v.length - 1])) {
73 this.push(this.t[this.t.length - 1], NaN)
74 }
75 }
76 build(): Series {
77 return { t: Float64Array.from(this.t), v: Float64Array.from(this.v) }
78 }
79}
81function wrapToPi(x: number): number {
82 const w = ((x + Math.PI) % (2 * Math.PI)) - Math.PI
83 return w < -Math.PI ? w + 2 * Math.PI : w
84}
86export function reconstruct(seq: ParsedSeq): Reconstructed {
87 const { rasters } = seq
88 const gradBuilders = {
89 gx: new SeriesBuilder(),
90 gy: new SeriesBuilder(),
91 gz: new SeriesBuilder(),
92 }
93 const rfMag = new SeriesBuilder()
94 const rfPhase = new SeriesBuilder()
95 const blockSpans: BlockSpan[] = []
96 const adcSpans: AdcSpan[] = []
97 const rfSpans: RfPulseSpan[] = []
99 const channels = ['gx', 'gy', 'gz'] as const
101 // Zero anchors so empty channels still draw a baseline
102 for (const ch of channels) gradBuilders[ch].push(0, 0)
103 rfMag.push(0, 0)
105 let t0 = 0
106 for (let bi = 0; bi < seq.blocks.length; bi++) {
107 const block = seq.blocks[bi]
108 blockSpans.push({ index: bi, id: block.id, start: t0, duration: block.duration })
110 for (const ch of channels) {
111 const gid = ch === 'gx' ? block.gxId : ch === 'gy' ? block.gyId : block.gzId
112 if (gid === 0) continue
113 const grad = seq.grads.get(gid)
114 if (!grad) continue
115 const b = gradBuilders[ch]
116 const gt0 = t0 + grad.delay
117 if (grad.kind === 'trap') {
118 if (Math.abs(grad.flat) > EPS) {
119 b.push(gt0, 0)
120 b.push(gt0 + grad.rise, grad.amp)
121 b.push(gt0 + grad.rise + grad.flat, grad.amp)
122 b.push(gt0 + grad.rise + grad.flat + grad.fall, 0)
123 } else if (Math.abs(grad.rise) > EPS && Math.abs(grad.fall) > EPS) {
124 b.push(gt0, 0)
125 b.push(gt0 + grad.rise, grad.amp)
126 b.push(gt0 + grad.rise + grad.fall, 0)
127 }
128 // else: empty gradient, skip
129 } else {
130 const shape = seq.shapes.get(grad.shapeId)
131 if (!shape) continue
132 const n = shape.numSamples
133 const raster = rasters.gradient
134 let tt: Float64Array
135 if (grad.timeShapeId === 0) {
136 tt = new Float64Array(n)
137 for (let i = 0; i < n; i++) tt[i] = (i + 0.5) * raster
138 } else {
139 const timeShape = seq.shapes.get(grad.timeShapeId)
140 if (!timeShape) continue
141 tt = new Float64Array(n)
142 for (let i = 0; i < n; i++) tt[i] = timeShape.samples[i] * raster
143 }
144 const onCenters = Math.abs(tt[0] - 0.5 * raster) < 1e-6 * raster
145 if (onCenters) {
146 // v1.4 files carry no first/last; fall back to the edge samples
147 const first = Number.isNaN(grad.first) ? grad.amp * shape.samples[0] : grad.first
148 const last =
149 Number.isNaN(grad.last) ? grad.amp * shape.samples[n - 1] : grad.last
150 const shapeDur = Math.ceil((tt[n - 1] - EPS) / raster) * raster
151 b.push(gt0, first)
152 for (let i = 0; i < n; i++) b.push(gt0 + tt[i], grad.amp * shape.samples[i])
153 b.push(gt0 + shapeDur, last)
154 } else {
155 // extended trapezoid: vertices as stored
156 for (let i = 0; i < n; i++) b.push(gt0 + tt[i], grad.amp * shape.samples[i])
157 }
158 }
159 }
161 if (block.rfId !== 0) {
162 const ev = seq.rf.get(block.rfId)
163 const magShape = ev ? seq.shapes.get(ev.magShapeId) : undefined
164 if (ev && magShape) {
165 const n = magShape.numSamples
166 const raster = rasters.rf
167 const phaseShape = ev.phaseShapeId !== 0 ? seq.shapes.get(ev.phaseShapeId) : undefined
168 let tt: Float64Array
169 if (ev.timeShapeId === 0) {
170 tt = new Float64Array(n)
171 for (let i = 0; i < n; i++) tt[i] = (i + 0.5) * raster
172 } else {
173 const timeShape = seq.shapes.get(ev.timeShapeId)
174 tt = new Float64Array(n)
175 for (let i = 0; i < n; i++) tt[i] = (timeShape?.samples[i] ?? 0) * raster
176 }
177 const rt0 = t0 + ev.delay
178 const end = rt0 + Math.ceil((tt[n - 1] - EPS) / raster) * raster
179 rfMag.push(rt0, 0)
180 rfPhase.breakLine()
181 for (let i = 0; i < n; i++) {
182 const t = rt0 + tt[i]
183 rfMag.push(t, Math.abs(ev.amp * magShape.samples[i]))
184 const ph = 2 * Math.PI * (phaseShape?.samples[i] ?? 0) + ev.phase + 2 * Math.PI * ev.freq * tt[i]
185 rfPhase.push(t, wrapToPi(ph))
186 }
187 rfPhase.breakLine()
188 rfMag.push(end, 0)
189 rfSpans.push({
190 blockIndex: bi,
191 start: rt0,
192 end,
193 center: Number.isNaN(ev.center) ? NaN : rt0 + ev.center,
194 use: ev.use,
195 })
196 }
197 }
199 if (block.adcId !== 0) {
200 const adc = seq.adcs.get(block.adcId)
201 if (adc) {
202 const start = t0 + adc.delay
203 adcSpans.push({
204 blockIndex: bi,
205 start,
206 end: start + adc.num * adc.dwell,
207 num: adc.num,
208 dwell: adc.dwell,
209 })
210 }
211 }
213 t0 += block.duration
214 }
216 for (const ch of channels) gradBuilders[ch].push(t0, 0)
217 rfMag.push(t0, 0)
219 const stats = {
220 maxGrad: { gx: 0, gy: 0, gz: 0 },
221 maxSlew: { gx: 0, gy: 0, gz: 0 },
222 rfCount: rfSpans.length,
223 adcCount: adcSpans.length,
224 adcSamples: adcSpans.reduce((acc, a) => acc + a.num, 0),
225 }
226 const series = {
227 gx: gradBuilders.gx.build(),
228 gy: gradBuilders.gy.build(),
229 gz: gradBuilders.gz.build(),
230 }
231 for (const ch of channels) {
232 const { t, v } = series[ch]
233 for (let i = 0; i < v.length; i++) {
234 if (Number.isNaN(v[i])) continue
235 stats.maxGrad[ch] = Math.max(stats.maxGrad[ch], Math.abs(v[i]))
236 if (i > 0 && !Number.isNaN(v[i - 1]) && t[i] > t[i - 1] + EPS) {
237 stats.maxSlew[ch] = Math.max(stats.maxSlew[ch], Math.abs((v[i] - v[i - 1]) / (t[i] - t[i - 1])))
238 }
239 }
240 }
242 return {
243 duration: t0,
244 blockSpans,
245 ...series,
246 rfMag: rfMag.build(),
247 rfPhase: rfPhase.build(),
248 adcSpans,
249 rfSpans,
250 stats,
251 }
252}