diff --git a/src/editor/seqEditorProvider.ts b/src/editor/seqEditorProvider.ts index 72eae3a..701d7c7 100644 --- a/src/editor/seqEditorProvider.ts +++ b/src/editor/seqEditorProvider.ts @@ -1,5 +1,5 @@ /** - * Custom Text Editor Provider for Pulseq .seq files. + * Custom readonly editor provider for Pulseq .seq files. * * Registered as `seqeyes.sequenceViewer` — opens automatically when the user * opens a `.seq` file. The provider: @@ -84,9 +84,17 @@ export async function exportKspaceToDirectoryForTest( return await writeKspaceArtifacts(sourceUri, saveUri, packageVersion, defaultStem); } +class SeqDocument implements vscode.CustomDocument { + constructor(public readonly uri: vscode.Uri) { } + + dispose(): void { + // This viewer is read-only and does not hold native resources. + } +} + // ─── Provider class ─────────────────────────────────────────────────────── -export class SeqEditorProvider implements vscode.CustomTextEditorProvider { +export class SeqEditorProvider implements vscode.CustomReadonlyEditorProvider { /** Register the provider with VS Code. */ static register(ctx: vscode.ExtensionContext): vscode.Disposable { @@ -98,10 +106,18 @@ export class SeqEditorProvider implements vscode.CustomTextEditorProvider { constructor(private readonly _ctx: vscode.ExtensionContext) { } - // ── resolveCustomTextEditor ────────────────────────────────────── + openCustomDocument( + uri: vscode.Uri, + _openContext: vscode.CustomDocumentOpenContext, + _token: vscode.CancellationToken, + ): SeqDocument { + return new SeqDocument(uri); + } - async resolveCustomTextEditor( - doc: vscode.TextDocument, + // ── resolveCustomEditor ────────────────────────────────────────── + + async resolveCustomEditor( + doc: SeqDocument, panel: vscode.WebviewPanel, _token: vscode.CancellationToken, ): Promise { @@ -123,10 +139,9 @@ export class SeqEditorProvider implements vscode.CustomTextEditorProvider { }; postProgress('start', 0, 'Reading file\u2026'); - const text = Buffer.from(await vscode.workspace.fs.readFile(uri)).toString('utf8'); - - postProgress('parse', 5, 'Parsing Pulseq sequence\u2026'); - const seq = parseSequenceText(text); + const seq = await readAndParseSequence(uri, () => { + postProgress('parse', 5, 'Parsing Pulseq sequence\u2026'); + }); postProgress('timing', 10, 'Detecting TR/TE timing\u2026'); const timing = detectSequenceTiming(seq); @@ -218,9 +233,8 @@ export class SeqEditorProvider implements vscode.CustomTextEditorProvider { } }; - // ── Initial load: validate, set full UI, show progress, then send data ── + // ── Initial load: set full UI, show progress, then send data ── try { - parseSequenceText(Buffer.from(await vscode.workspace.fs.readFile(doc.uri)).toString('utf8')); panel.webview.html = getWebviewContent(0); // Give the webview a moment to parse its new HTML, then start progress panel.webview.postMessage({ type: 'progress', phase: 'start', percent: 0, text: 'Preparing\u2026' }); @@ -491,6 +505,13 @@ function serializeBlocks(blocks: DecodedBlock[]): object[] { }); } +async function readAndParseSequence(uri: vscode.Uri, didRead: () => void) { + const fileBytes = await vscode.workspace.fs.readFile(uri); + didRead(); + const text = Buffer.from(fileBytes.buffer, fileBytes.byteOffset, fileBytes.byteLength).toString('utf8'); + return parseSequenceText(text); +} + function serializeGrad(g: DecodedGradWaveform): Record { return { s: g.startTime, d: g.duration, diff --git a/src/editor/webview/assets/derived-series.js b/src/editor/webview/assets/derived-series.js index 492529b..6b39331 100644 --- a/src/editor/webview/assets/derived-series.js +++ b/src/editor/webview/assets/derived-series.js @@ -151,3 +151,111 @@ function forEachDerivedPoint(series,viewStart,viewEnd,maxPoints,visit){ emitRange(Math.max(i0,lastBucket*size),i1); return emitted; } + +/* Block-indexed min/max summaries for RF, gradients, and ADC occupancy. */ +function createWaveformOverview(blocks){ + if(!blocks||!blocks.length)return null; + var level=buildWaveformOverviewLevel(blocks,8),levels=[level]; + while(level.count>1){level=mergeWaveformOverviewLevel(level);levels.push(level);} + return{levels:levels,blockCount:blocks.length,pointPrefix:createWaveformPointPrefixes(blocks)}; +} + +function createEmptyWaveformOverviewLevel(count,bucketSize){ + var level={count:count,bucketSize:bucketSize,t0:new Float64Array(count),t1:new Float64Array(count), + rfMin:new Float64Array(count),rfMax:new Float64Array(count),gxMin:new Float64Array(count),gxMax:new Float64Array(count), + gyMin:new Float64Array(count),gyMax:new Float64Array(count),gzMin:new Float64Array(count),gzMax:new Float64Array(count), + adcStart:new Float64Array(count),adcEnd:new Float64Array(count)}; + var mins=[level.rfMin,level.gxMin,level.gyMin,level.gzMin,level.adcStart]; + var maxs=[level.rfMax,level.gxMax,level.gyMax,level.gzMax,level.adcEnd]; + for(var m=0;m1){var step=Math.max(1,Math.ceil(adc.n/200));adcPhase=Math.floor(adc.n/step)+1;} + prefix.phase[i+1]=prefix.phase[i]+rfPhase+adcPhase; + prefix.gx[i+1]=prefix.gx[i]+overviewGradientPointCount(block.gx); + prefix.gy[i+1]=prefix.gy[i]+overviewGradientPointCount(block.gy); + prefix.gz[i+1]=prefix.gz[i]+overviewGradientPointCount(block.gz); + prefix.adc[i+1]=prefix.adc[i]+(adc?1:0); + } + return prefix; +} + +function overviewGradientPointCount(gradient){ + return gradient&&gradient.ty!=='none'&&gradient.t&&gradient.w?Math.min(gradient.t.length,gradient.w.length):0; +} + +function includeOverviewValues(values,minArray,maxArray,bucket){ + if(!values)return; + for(var i=0;imaxArray[bucket])maxArray[bucket]=value; + } +} + +function buildWaveformOverviewLevel(blocks,bucketSize){ + var count=Math.ceil(blocks.length/bucketSize),level=createEmptyWaveformOverviewLevel(count,bucketSize); + for(var bucket=0;bucketlevel.adcEnd[bucket])level.adcEnd[bucket]=adcEnd; + } + } + } + return level; +} + +function mergeWaveformOverviewLevel(child){ + var count=Math.ceil(child.count/4),level=createEmptyWaveformOverviewLevel(count,child.bucketSize*4); + var channels=[['rfMin','rfMax'],['gxMin','gxMax'],['gyMin','gyMax'],['gzMin','gzMax'],['adcStart','adcEnd']]; + for(var bucket=0;bucketlevel[maxKey][bucket])level[maxKey][bucket]=child[maxKey][i]; + } + } + } + return level; +} + +function waveformVisiblePointCount(overview,key,startBlock,endBlock){ + if(!overview||!overview.pointPrefix||!overview.pointPrefix[key])return 0; + var prefix=overview.pointPrefix[key],start=Math.max(0,Math.min(startBlock,prefix.length-1)),end=Math.max(start,Math.min(endBlock,prefix.length-1)); + return prefix[end]-prefix[start]; +} + +function selectWaveformOverview(overview,startBlock,endBlock,maxBuckets){ + if(!overview||endBlock<=startBlock)return null; + var selected=overview.levels[overview.levels.length-1]; + for(var i=0;i=2&&i<=4)||(i>=8&&i<=10)){ctx.beginPath();ctx.moveTo(M.l,cy(vi));ctx.lineTo(w-M.r,cy(vi));ctx.stroke();}} ctx.setLineDash([]); } -function drawBlockBounds(w,h,s){ +function drawBlockBounds(w,h,vs,ve,s){ ctx.strokeStyle=s.getPropertyValue('--ax').trim();ctx.lineWidth=0.6;ctx.setLineDash([3,6]); - for(var i=0;iw-M.r)continue; ctx.beginPath();ctx.moveTo(x,M.t);ctx.lineTo(x,h-M.b);ctx.stroke(); ctx.fillStyle=s.getPropertyValue('--ax').trim();ctx.font='9px monospace';ctx.textAlign='center'; @@ -232,13 +233,80 @@ function drawBlocks(vs,ve,s){ adc:s.getPropertyValue('--adc').trim(),adf:s.getPropertyValue('--adf').trim(), tr:s.getPropertyValue('--tr').trim(),fg:s.getPropertyValue('--fg').trim() }; - if(rows[0]>=0)drawRfBlocks(range.start,range.end,rows[0],ch,colors,vs,ve); - if(rows[1]>=0)drawPhaseBlocks(range.start,range.end,rows[1],ch,colors,vs,ve); - if(rows[2]>=0)drawGradientBlocks(range.start,range.end,'gx',rows[2],2,ch,colors.gx,vs,ve); - if(rows[3]>=0)drawGradientBlocks(range.start,range.end,'gy',rows[3],3,ch,colors.gy,vs,ve); - if(rows[4]>=0)drawGradientBlocks(range.start,range.end,'gz',rows[4],4,ch,colors.gz,vs,ve); - if(rows[5]>=0)drawAdcBlocks(range.start,range.end,rows[5],ch,colors,vs,ve); - if(rows[6]>=0)drawTriggerBlocks(range.start,range.end,rows[6],ch,colors,vs,ve); + var pixelBudget=Math.max(1,Math.floor(plotWidth())); + var overview=selectWaveformOverview(waveformOverview,range.start,range.end,pixelBudget); + function useOverview(key){return !!overview&&waveformVisiblePointCount(waveformOverview,key,range.start,range.end)>pixelBudget;} + if(rows[0]>=0){if(useOverview('rf'))drawRfOverview(overview,rows[0],ch,colors,vs,ve);else drawRfBlocks(range.start,range.end,rows[0],ch,colors,vs,ve);} + if(rows[1]>=0){if(useOverview('phase'))drawPhaseSampled(range.start,range.end,rows[1],ch,colors,vs,ve,pixelBudget);else drawPhaseBlocks(range.start,range.end,rows[1],ch,colors,vs,ve);} + if(rows[2]>=0){if(useOverview('gx'))drawGradientOverview(overview,'gx',rows[2],2,ch,colors.gx,vs,ve);else drawGradientBlocks(range.start,range.end,'gx',rows[2],2,ch,colors.gx,vs,ve);} + if(rows[3]>=0){if(useOverview('gy'))drawGradientOverview(overview,'gy',rows[3],3,ch,colors.gy,vs,ve);else drawGradientBlocks(range.start,range.end,'gy',rows[3],3,ch,colors.gy,vs,ve);} + if(rows[4]>=0){if(useOverview('gz'))drawGradientOverview(overview,'gz',rows[4],4,ch,colors.gz,vs,ve);else drawGradientBlocks(range.start,range.end,'gz',rows[4],4,ch,colors.gz,vs,ve);} + if(rows[5]>=0){if(useOverview('adc'))drawAdcOverview(overview,rows[5],ch,colors,vs,ve);else drawAdcBlocks(range.start,range.end,rows[5],ch,colors,vs,ve);} + if(rows[6]>=0&&range.end-range.start<=pixelBudget)drawTriggerBlocks(range.start,range.end,rows[6],ch,colors,vs,ve); +} + +function drawRfOverview(summary,vi,ch,colors,vs,ve){ + var level=summary.level,y=cy(vi),scale=ch*.9/channelRange(0); + rowClip(vi,ch,function(){ + ctx.strokeStyle=colors.rf;ctx.lineWidth=1;ctx.beginPath(); + for(var i=summary.first;ive||level.rfMin[i]===Infinity)continue; + var x=t2x(.5*(level.t0[i]+level.t1[i])),y0=y+ch*.45-level.rfMin[i]*scale,y1=y+ch*.45-level.rfMax[i]*scale; + if(Math.abs(y1-y0)<1){y0-=.5;y1+=.5;}ctx.moveTo(x,y0);ctx.lineTo(x,y1); + } + ctx.stroke(); + }); +} + +function drawPhaseSampled(start,end,vi,ch,colors,vs,ve,maxPoints){ + var prefix=waveformOverview.pointPrefix.phase,first=prefix[start],last=prefix[end],count=last-first; + if(count<=0)return; + var y=cy(vi),scale=ch*.9/channelRange(1),sampleCount=Math.min(maxPoints,count),sampleStep=count/sampleCount; + rowClip(vi,ch,function(){ + for(var sampleIndex=0;sampleIndex>1;if(prefix[mid+1]<=ordinal)lo=mid+1;else hi=mid;} + if(lo>=end)continue; + var block=BL[lo],rf=block.rf,rfCount=rf&&rf.t&&rf.p?Math.min(rf.t.length,rf.p.length):0; + var local=ordinal-prefix[lo],time=NaN,value=NaN,isAdc=false; + if(local1){ + var adc=block.adc,adcStep=Math.max(1,Math.ceil(adc.n/200)),adcSample=(local-rfCount)*adcStep; + var adcStart=adc.s+adc.d,adcEnd=adcStart+adc.n*adc.dw; + time=adcSampleve)continue; + ctx.fillStyle=isAdc?colors.adc:colors.rf; + ctx.fillRect(t2x(time)-.5,y+ch*.45-value*scale-.5,1,1); + } + }); +} + +function drawGradientOverview(summary,key,vi,ci,ch,color,vs,ve){ + var level=summary.level,minValues=level[key+'Min'],maxValues=level[key+'Max'],y=cy(vi),scale=ch*.4/channelRange(ci); + rowClip(vi,ch,function(){ + ctx.strokeStyle=color;ctx.lineWidth=1;ctx.beginPath(); + for(var i=summary.first;ive||minValues[i]===Infinity)continue; + var x=t2x(.5*(level.t0[i]+level.t1[i])),y0=y-minValues[i]*scale,y1=y-maxValues[i]*scale; + if(Math.abs(y1-y0)<1){y0-=.5;y1+=.5;}ctx.moveTo(x,y0);ctx.lineTo(x,y1); + } + ctx.stroke(); + }); +} + +function drawAdcOverview(summary,vi,ch,colors,vs,ve){ + var level=summary.level,y=cy(vi); + rowClip(vi,ch,function(){ + ctx.fillStyle=colors.adf;ctx.strokeStyle=colors.adc;ctx.lineWidth=1;ctx.beginPath();var hasRect=false; + for(var i=summary.first;ive)continue; + var x0=t2x(Math.max(vs,level.adcStart[i])),x1=t2x(Math.min(ve,level.adcEnd[i])); + ctx.rect(x0,y-ch*.28,Math.max(1,x1-x0),ch*.56);hasRect=true; + } + if(hasRect){ctx.fill();ctx.stroke();} + }); } function visibleBlockRange(vs,ve){ diff --git a/src/editor/webview/assets/interaction.js b/src/editor/webview/assets/interaction.js index f4c13a0..92c548f 100644 --- a/src/editor/webview/assets/interaction.js +++ b/src/editor/webview/assets/interaction.js @@ -2,7 +2,7 @@ Mouse interaction ═══════════════════════════════════════════════════════════════════════ */ -var _touchTooltipTimer=null,_pointerFrame=0,_pendingPointer=null; +var _touchTooltipTimer=null,_pointerFrame=0,_pendingPointer=null,mmDrag=false; mc.addEventListener('wheel',function(e){e.preventDefault(); var changed=false; @@ -137,7 +137,7 @@ function placeTooltip(cx,cy){ function showTooltipAt(cx,cy,ct){ var ch=cH(),vc=visChannels(),vi2=Math.floor((cy-mc.getBoundingClientRect().top-M.t)/ch); if(vi2>=0&&vi2=BL[i].s&&ct<=BL[i].s+BL[i].d){found=BL[i];break}} + var found=findBlockAtTime(ct); var lines=[]; if(found){ var blockDt=Math.max(0,Math.min(found.d,ct-found.s)); @@ -164,6 +164,15 @@ function showTooltipAt(cx,cy,ct){ tt.style.display='none'; } +function findBlockAtTime(t){ + var lo=0,hi=BL.length; + while(lo>1;if(BL[mid].s+BL[mid].d=BL[i].s&&t<=BL[i].s+BL[i].d)return BL[i]; + } + return null; +} + /* ── Toolbar buttons ──────────────────────────────────────────────────── */ document.getElementById('openBtn').onclick=function(){ if(vscApi){vscApi.postMessage({command:'openFile'});} diff --git a/src/editor/webview/assets/state.js b/src/editor/webview/assets/state.js index 809b097..35c7968 100644 --- a/src/editor/webview/assets/state.js +++ b/src/editor/webview/assets/state.js @@ -7,7 +7,7 @@ var cc=document.getElementById('cc'),mc=document.getElementById('mc'),ctx=mc.get tuSel=document.getElementById('tu'),guSel=document.getElementById('gu'); var exportBtn=document.getElementById('exportKspaceBtn'); var pnsBtn=document.getElementById('pnsBtn'); -var BL=[],TD=0,GR=1e-5,RR=1e-6,AR=1e-7,BR=1e-5; // blocks, duration, rasters [s] +var BL=[],waveformOverview=null,TD=0,GR=1e-5,RR=1e-6,AR=1e-7,BR=1e-5; // blocks, duration, rasters [s] var M={t:8,r:30,b:22,l:92}; // margins var CH=['RF','\u03c6','Gx','Gy','Gz','ADC','Trig','PNS','M1x','M1y','M1z']; var chColors=['var(--rf)','var(--rf)','var(--gx)','var(--gy)','var(--gz)','var(--adc)','var(--tr)','var(--fg)','var(--gx)','var(--gy)','var(--gz)']; @@ -206,6 +206,7 @@ window.addEventListener('message',function(e){ } if(m.type==='sequenceData'){ BL=m.blocks||[];TD=m.totalDuration||0;GR=m.gradRaster||1e-5; + waveformOverview=createWaveformOverview(BL); RR=m.rfRaster||RR;AR=m.adcRaster||AR;BR=m.blockRaster||BR; blockPos=m.blockPositions||[]; mmCache=null; // invalidate minimap cache on new data @@ -488,7 +489,8 @@ function drawMinimap(){ mmCtx.strokeStyle=s.getPropertyValue('--fg').trim(); mmCtx.globalAlpha=0.10;mmCtx.lineWidth=0.4; mmCtx.beginPath(); - for(var tr=0;tr<=seqTiming.trCount;tr++){ + var trStep=Math.max(1,Math.ceil(seqTiming.trCount/Math.max(1,W))); + for(var tr=0;tr<=seqTiming.trCount;tr+=trStep){ var tx=tr*seqTiming.trTimeSec*scM; if(tx<=W){mmCtx.moveTo(tx,0);mmCtx.lineTo(tx,H);} } diff --git a/src/extension.ts b/src/extension.ts index e2d0209..db8b48e 100644 --- a/src/extension.ts +++ b/src/extension.ts @@ -28,8 +28,7 @@ export function activate(context: vscode.ExtensionContext): void { }), vscode.commands.registerCommand('seqeyes.openSequenceViewer', async (uri?: vscode.Uri) => { if (!uri) { - const editor = vscode.window.activeTextEditor; - if (editor?.document.fileName.endsWith('.seq')) uri = editor.document.uri; + uri = getActiveSequenceUri(); } if (uri) { await vscode.commands.executeCommand('vscode.openWith', uri, 'seqeyes.sequenceViewer'); @@ -69,3 +68,20 @@ export function activate(context: vscode.ExtensionContext): void { export function deactivate(): void { console.log('SeqEyes Plugin deactivated'); } + +function getActiveSequenceUri(): vscode.Uri | undefined { + const editor = vscode.window.activeTextEditor; + if (editor?.document.fileName.endsWith('.seq')) { + return editor.document.uri; + } + + const input = vscode.window.tabGroups.activeTabGroup.activeTab?.input; + if (input instanceof vscode.TabInputText && input.uri.path.endsWith('.seq')) { + return input.uri; + } + if (input instanceof vscode.TabInputCustom && input.uri.path.endsWith('.seq')) { + return input.uri; + } + + return undefined; +} diff --git a/src/pulseq/m1.ts b/src/pulseq/m1.ts index 27d9878..482bb47 100644 --- a/src/pulseq/m1.ts +++ b/src/pulseq/m1.ts @@ -132,42 +132,27 @@ function normalizeReferenceMode(mode: M1ReferenceMode | undefined): M1ReferenceM } function collectGradientSeries(blocks: DecodedBlock[], channel: 'gx' | 'gy' | 'gz'): GradientSeries { - const time: number[] = []; - const value: number[] = []; + const series: GradientSeries = { time: [], value: [] }; for (const block of blocks) { const grad = block[channel] as DecodedGradWaveform | undefined; if (!grad?.timePoints || !grad.waveform) continue; const n = Math.min(grad.timePoints.length, grad.waveform.length); for (let i = 0; i < n; i++) { - time.push(grad.timePoints[i]); - value.push(grad.waveform[i]); + appendGradientPoint(series, grad.timePoints[i], grad.waveform[i]); } } - return sanitizeGradientSeries(time, value); + return series; } -function sanitizeGradientSeries(time: number[], value: number[]): GradientSeries { - const pairs: Array<[number, number]> = []; - const n = Math.min(time.length, value.length); - for (let i = 0; i < n; i++) { - const t = time[i]; - const v = value[i]; - if (Number.isFinite(t) && Number.isFinite(v)) pairs.push([t, v]); +function appendGradientPoint(series: GradientSeries, t: number, value: number): void { + if (!Number.isFinite(t) || !Number.isFinite(value)) return; + const last = series.time.length - 1; + if (last >= 0 && Math.abs(t - series.time[last]) <= TIME_EPS) { + series.value[last] = 0.5 * (series.value[last] + value); + } else if (last < 0 || t > series.time[last]) { + series.time.push(t); + series.value.push(value); } - pairs.sort((a, b) => a[0] - b[0]); - - const outT: number[] = []; - const outV: number[] = []; - for (const [t, v] of pairs) { - const last = outT.length - 1; - if (last >= 0 && Math.abs(t - outT[last]) <= TIME_EPS) { - outV[last] = 0.5 * (outV[last] + v); - continue; - } - outT.push(t); - outV.push(v); - } - return { time: outT, value: outV }; } function collectRfEvents(blocks: DecodedBlock[], warnings: string[]): RfEvent[] { @@ -237,6 +222,30 @@ function walkM1( let currentT = tReset; let unsignedM0 = 0; let unsignedM1 = 0; + let gradientIndex = -1; + + const seekGradient = (t: number): void => { + while (gradientIndex + 1 < gradient.time.length + && gradient.time[gradientIndex + 1] <= t + TIME_EPS) { + gradientIndex++; + } + }; + + const sampleGradient = (t: number): number => { + const n = gradient.time.length; + if (n === 0 || t < gradient.time[0] - TIME_EPS || t > gradient.time[n - 1] + TIME_EPS) return 0; + seekGradient(t); + if (gradientIndex < 0) return 0; + if (gradientIndex >= n - 1 || Math.abs(t - gradient.time[gradientIndex]) <= TIME_EPS) { + return gradient.value[gradientIndex]; + } + const t0 = gradient.time[gradientIndex]; + const t1 = gradient.time[gradientIndex + 1]; + if (!(t1 > t0)) return gradient.value[gradientIndex]; + const alpha = (t - t0) / (t1 - t0); + return gradient.value[gradientIndex] + + alpha * (gradient.value[gradientIndex + 1] - gradient.value[gradientIndex]); + }; const reportedM1At = (t: number): number => { if (referenceMode === 'observationTime') return sign * (unsignedM1 - (t - tReset) * unsignedM0); @@ -246,10 +255,13 @@ function walkM1( const advanceTo = (targetT: number): void => { if (!(targetT > currentT + TIME_EPS)) return; while (currentT < targetT - TIME_EPS) { - let nextT = nextGradientBreakpoint(gradient.time, currentT, targetT); + seekGradient(currentT); + let nextT = gradientIndex + 1 < gradient.time.length + ? Math.min(targetT, gradient.time[gradientIndex + 1]) + : targetT; if (!(nextT > currentT)) nextT = targetT; - const ga = sampleGradientAt(gradient, currentT); - const gb = sampleGradientAt(gradient, nextT); + const ga = sampleGradient(currentT); + const gb = sampleGradient(nextT); const [m0Seg, m1Seg] = integrateLinearSegment(currentT, nextT, tReset, ga, gb); unsignedM0 += m0Seg; unsignedM1 += m1Seg; @@ -294,38 +306,6 @@ function walkM1( return { t: outT, m1: outM1 }; } -function sampleGradientAt(gradient: GradientSeries, t: number): number { - const n = gradient.time.length; - if (n <= 0 || t < gradient.time[0] || t > gradient.time[n - 1]) return 0; - if (n === 1 || t <= gradient.time[0]) return gradient.value[0]; - if (t >= gradient.time[n - 1]) return gradient.value[n - 1]; - let lo = 0; - let hi = n - 1; - while (hi - lo > 1) { - const mid = (lo + hi) >> 1; - if (gradient.time[mid] <= t) lo = mid; - else hi = mid; - } - const t0 = gradient.time[lo]; - const t1 = gradient.time[hi]; - if (!(t1 > t0)) return gradient.value[lo]; - const alpha = (t - t0) / (t1 - t0); - return gradient.value[lo] + alpha * (gradient.value[hi] - gradient.value[lo]); -} - -function nextGradientBreakpoint(times: number[], t: number, target: number): number { - if (times.length <= 1 || t >= times[times.length - 1]) return target; - let lo = 0; - let hi = times.length; - const threshold = t + TIME_EPS; - while (lo < hi) { - const mid = (lo + hi) >> 1; - if (times[mid] <= threshold) lo = mid + 1; - else hi = mid; - } - return lo < times.length ? Math.min(target, times[lo]) : target; -} - function integrateLinearSegment(a: number, b: number, tRef: number, ga: number, gb: number): [number, number] { const h = b - a; if (!(h > 0)) return [0, 0]; diff --git a/src/pulseq/pns.ts b/src/pulseq/pns.ts index 16b5609..a76d600 100644 --- a/src/pulseq/pns.ts +++ b/src/pulseq/pns.ts @@ -130,18 +130,6 @@ export function calculatePns( const nSamples = Math.floor(ntMax - ntMin + 1.0); if (nSamples < 2) return invalidPns('Too few samples for PNS computation.'); - const tAxis = new Float64Array(nSamples); - const gxTpm = new Float64Array(nSamples); - const gyTpm = new Float64Array(nSamples); - const gzTpm = new Float64Array(nSamples); - for (let i = 0; i < nSamples; i++) { - const tSec = (ntMin + i) * dtSec; - tAxis[i] = tSec; - gxTpm[i] = interpLinearZero(waves[0], tSec) / gammaHzPerT; - gyTpm[i] = interpLinearZero(waves[1], tSec) / gammaHzPerT; - gzTpm[i] = interpLinearZero(waves[2], tSec) / gammaHzPerT; - } - const longestTauMs = Math.max( hardware.x.tau1Ms, hardware.x.tau2Ms, hardware.x.tau3Ms, hardware.y.tau1Ms, hardware.y.tau2Ms, hardware.y.tau3Ms, @@ -151,13 +139,9 @@ export function calculatePns( const preCount = Math.max(0, Math.round(zptSec / (4.0 * dtSec))); const postCount = Math.max(0, Math.round(zptSec / dtSec)); - const gxPadded = padSamples(gxTpm, preCount, postCount); - const gyPadded = padSamples(gyTpm, preCount, postCount); - const gzPadded = padSamples(gzTpm, preCount, postCount); - - const stimX = safePnsModel(diff(gxPadded, dtSec), dtSec, hardware.x); - const stimY = safePnsModel(diff(gyPadded, dtSec), dtSec, hardware.y); - const stimZ = safePnsModel(diff(gzPadded, dtSec), dtSec, hardware.z); + const stimX = calculatePnsAxis(waves[0], ntMin, nSamples, preCount, postCount, dtSec, gammaHzPerT, hardware.x); + const stimY = calculatePnsAxis(waves[1], ntMin, nSamples, preCount, postCount, dtSec, gammaHzPerT, hardware.y); + const stimZ = calculatePnsAxis(waves[2], ntMin, nSamples, preCount, postCount, dtSec, gammaHzPerT, hardware.z); const hasAnyNonTrap = blocks.some(block => ( block.gx?.type === 'arb' || block.gy?.type === 'arb' || block.gz?.type === 'arb' @@ -165,63 +149,83 @@ export function calculatePns( const hasAnyLabelExt = blocks.some(block => !!(block.labelSets?.length || block.labelIncs?.length)); const shift = hasAnyNonTrap || hasAnyLabelExt ? 1 : 0; - const selectedX: number[] = []; - const selectedY: number[] = []; - const selectedZ: number[] = []; - const selectedT: number[] = []; + let selectedCount = 0; for (let origIdx = 0; origIdx < nSamples; origIdx++) { const paddedIdx = preCount + origIdx; let stimIdx = paddedIdx - shift; - if (shift > 0 && hasAnyLabelExt && origIdx === tAxis.length - 1) { + if (shift > 0 && hasAnyLabelExt && origIdx === nSamples - 1) { stimIdx = Math.min(paddedIdx, stimX.length - 1); } if (stimIdx < 0 || stimIdx >= stimX.length || stimIdx >= stimY.length || stimIdx >= stimZ.length) continue; - selectedX.push(stimX[stimIdx]); - selectedY.push(stimY[stimIdx]); - selectedZ.push(stimZ[stimIdx]); - selectedT.push(tAxis[origIdx]); + selectedCount++; } - const timeSec = new Float64Array(selectedX.length); - const pnsX = new Float64Array(selectedX.length); - const pnsY = new Float64Array(selectedX.length); - const pnsZ = new Float64Array(selectedX.length); - const pnsNorm = new Float64Array(selectedX.length); + const timeSec = new Float64Array(selectedCount); + const pnsX = new Float64Array(selectedCount); + const pnsY = new Float64Array(selectedCount); + const pnsZ = new Float64Array(selectedCount); + const pnsNorm = new Float64Array(selectedCount); let ok = true; - for (let i = 0; i < selectedX.length; i++) { - const xNorm = 0.01 * selectedX[i]; - const yNorm = 0.01 * selectedY[i]; - const zNorm = 0.01 * selectedZ[i]; + let selectedIndex = 0; + for (let origIdx = 0; origIdx < nSamples; origIdx++) { + const paddedIdx = preCount + origIdx; + let stimIdx = paddedIdx - shift; + if (shift > 0 && hasAnyLabelExt && origIdx === nSamples - 1) { + stimIdx = Math.min(paddedIdx, stimX.length - 1); + } + if (stimIdx < 0 || stimIdx >= stimX.length || stimIdx >= stimY.length || stimIdx >= stimZ.length) continue; + const xNorm = 0.01 * stimX[stimIdx]; + const yNorm = 0.01 * stimY[stimIdx]; + const zNorm = 0.01 * stimZ[stimIdx]; const norm = Math.sqrt(xNorm * xNorm + yNorm * yNorm + zNorm * zNorm); - timeSec[i] = selectedT[i]; - pnsX[i] = xNorm; - pnsY[i] = yNorm; - pnsZ[i] = zNorm; - pnsNorm[i] = norm; + timeSec[selectedIndex] = (ntMin + origIdx) * dtSec; + pnsX[selectedIndex] = xNorm; + pnsY[selectedIndex] = yNorm; + pnsZ[selectedIndex] = zNorm; + pnsNorm[selectedIndex] = norm; if (norm >= 1.0) ok = false; + selectedIndex++; } return { valid: true, ok, timeSec, pnsX, pnsY, pnsZ, pnsNorm }; } export function safePnsModel(dgdt: Float64Array, dtSec: number, hw: PnsAxisHardware): Float64Array { - const absDgdt = new Float64Array(dgdt.length); - for (let i = 0; i < dgdt.length; i++) absDgdt[i] = Math.abs(dgdt[i]); + return runPnsModel(dgdt.length, index => dgdt[index], dtSec, hw); +} + +function runPnsModel( + length: number, + derivativeAt: (index: number) => number, + dtSec: number, + hw: PnsAxisHardware, +): Float64Array { const dtMs = dtSec * 1000.0; - const lp1 = lowpassTau(dgdt, hw.tau1Ms, dtMs); - const lp2 = lowpassTau(absDgdt, hw.tau2Ms, dtMs); - const lp3 = lowpassTau(dgdt, hw.tau3Ms, dtMs); - const stim = new Float64Array(dgdt.length); + const alpha1 = lowpassAlpha(hw.tau1Ms, dtMs); + const alpha2 = lowpassAlpha(hw.tau2Ms, dtMs); + const alpha3 = lowpassAlpha(hw.tau3Ms, dtMs); + const stim = new Float64Array(length); const denom = hw.stimLimit > 0 ? hw.stimLimit : 1; - for (let i = 0; i < dgdt.length; i++) { - const s1 = hw.a1 * Math.abs(lp1[i]); - const s2 = hw.a2 * lp2[i]; - const s3 = hw.a3 * Math.abs(lp3[i]); + let lp1 = 0; + let lp2 = 0; + let lp3 = 0; + for (let i = 0; i < length; i++) { + const derivative = derivativeAt(i); + lp1 = alpha1 * derivative + (1.0 - alpha1) * lp1; + lp2 = alpha2 * Math.abs(derivative) + (1.0 - alpha2) * lp2; + lp3 = alpha3 * derivative + (1.0 - alpha3) * lp3; + const s1 = hw.a1 * Math.abs(lp1); + const s2 = hw.a2 * lp2; + const s3 = hw.a3 * Math.abs(lp3); stim[i] = ((s1 + s2 + s3) / denom) * hw.gScale * 100.0; } return stim; } +function lowpassAlpha(tauMs: number, dtMs: number): number { + return tauMs <= 0 || dtMs <= 0 ? 1 : dtMs / (tauMs + dtMs); +} + function invalidPns(error: string): PnsResult { return { valid: false, @@ -341,82 +345,69 @@ function hasValidWeights(hw: PnsAxisHardware): boolean { return Math.abs(hw.a1 + hw.a2 + hw.a3 - 1.0) <= 1e-2 && hw.stimLimit > 0; } -function lowpassTau(input: Float64Array, tauMs: number, dtMs: number): Float64Array { - const out = new Float64Array(input.length); - if (!input.length) return out; - if (tauMs <= 0 || dtMs <= 0) { - out.set(input); - return out; - } - const alpha = dtMs / (tauMs + dtMs); - out[0] = alpha * input[0]; - for (let i = 1; i < input.length; i++) out[i] = alpha * input[i] + (1.0 - alpha) * out[i - 1]; - return out; -} - function collectGradientSeries(blocks: DecodedBlock[], channel: 'gx' | 'gy' | 'gz'): GradientSeries { - const time: number[] = []; - const value: number[] = []; + const series: GradientSeries = { time: [], value: [] }; for (const block of blocks) { const grad = block[channel] as DecodedGradWaveform | undefined; if (!grad?.timePoints || !grad.waveform) continue; const n = Math.min(grad.timePoints.length, grad.waveform.length); for (let i = 0; i < n; i++) { - time.push(grad.timePoints[i]); - value.push(grad.waveform[i]); + appendGradientPoint(series, grad.timePoints[i], grad.waveform[i]); } } - return sanitizeGradientSeries(time, value); + return series; } -function sanitizeGradientSeries(time: number[], value: number[]): GradientSeries { - const pairs: Array<[number, number]> = []; - const n = Math.min(time.length, value.length); - for (let i = 0; i < n; i++) { - if (Number.isFinite(time[i]) && Number.isFinite(value[i])) pairs.push([time[i], value[i]]); +function appendGradientPoint(series: GradientSeries, t: number, value: number): void { + if (!Number.isFinite(t) || !Number.isFinite(value)) return; + const last = series.time.length - 1; + if (last >= 0 && Math.abs(t - series.time[last]) <= TIME_EPS) { + series.value[last] = 0.5 * (series.value[last] + value); + } else if (last < 0 || t > series.time[last]) { + series.time.push(t); + series.value.push(value); } - pairs.sort((a, b) => a[0] - b[0]); - const outT: number[] = []; - const outV: number[] = []; - for (const [t, v] of pairs) { - const last = outT.length - 1; - if (last >= 0 && Math.abs(t - outT[last]) <= TIME_EPS) { - outV[last] = 0.5 * (outV[last] + v); - continue; - } - outT.push(t); - outV.push(v); - } - return { time: outT, value: outV }; } -function interpLinearZero(series: GradientSeries, t: number): number { - const n = series.time.length; - if (!n || t < series.time[0] || t > series.time[n - 1]) return 0; - if (n === 1 || t <= series.time[0]) return series.value[0]; - if (t >= series.time[n - 1]) return series.value[n - 1]; - let lo = 0; - let hi = n - 1; - while (hi - lo > 1) { - const mid = (lo + hi) >> 1; - if (series.time[mid] <= t) lo = mid; - else hi = mid; - } - const t0 = series.time[lo]; - const t1 = series.time[hi]; - if (!(t1 > t0)) return series.value[lo]; - const alpha = (t - t0) / (t1 - t0); - return series.value[lo] + alpha * (series.value[hi] - series.value[lo]); -} +function calculatePnsAxis( + series: GradientSeries, + ntMin: number, + nSamples: number, + preCount: number, + postCount: number, + dtSec: number, + gammaHzPerT: number, + hardware: PnsAxisHardware, +): Float64Array { + const sampleGradient = createGradientSampler(series); + const totalSamples = preCount + nSamples + postCount; + const paddedValue = (index: number): number => { + if (index < preCount || index >= preCount + nSamples) return 0; + const rasterIndex = index - preCount; + return sampleGradient((ntMin + rasterIndex) * dtSec) / gammaHzPerT; + }; -function padSamples(input: Float64Array, preCount: number, postCount: number): Float64Array { - const out = new Float64Array(preCount + input.length + postCount); - out.set(input, preCount); - return out; + let previous = paddedValue(0); + return runPnsModel(Math.max(0, totalSamples - 1), index => { + const current = paddedValue(index + 1); + const derivative = (current - previous) / dtSec; + previous = current; + return derivative; + }, dtSec, hardware); } -function diff(input: Float64Array, dtSec: number): Float64Array { - const out = new Float64Array(Math.max(0, input.length - 1)); - for (let i = 0; i < out.length; i++) out[i] = (input[i + 1] - input[i]) / dtSec; - return out; +function createGradientSampler(series: GradientSeries): (t: number) => number { + let index = -1; + return (t: number): number => { + const n = series.time.length; + if (n === 0 || t < series.time[0] || t > series.time[n - 1]) return 0; + while (index + 1 < n && series.time[index + 1] <= t + TIME_EPS) index++; + if (index < 0) return 0; + if (index >= n - 1 || t <= series.time[index] + TIME_EPS) return series.value[index]; + const t0 = series.time[index]; + const t1 = series.time[index + 1]; + if (!(t1 > t0)) return series.value[index]; + const alpha = (t - t0) / (t1 - t0); + return series.value[index] + alpha * (series.value[index + 1] - series.value[index]); + }; } diff --git a/src/pulseq/reader.ts b/src/pulseq/reader.ts index 23e4b4b..08b672d 100644 --- a/src/pulseq/reader.ts +++ b/src/pulseq/reader.ts @@ -28,25 +28,29 @@ import { /** Parse a .seq file from its text content. */ export function parseSequenceText(text: string): PulseqSequence { - const lines = text.split(/\r?\n/); const seq = createEmptySequence(); const seenSections = new Set(); + const shapeParser = new ShapeSectionParser(seq); let sectionName: string | null = null; let sectionLines: string[] = []; - for (const line of lines) { + forEachLine(text, line => { const m = line.match(/^\[(\w+)\]$/); if (m) { - if (sectionName) dispatchSection(seq, sectionName, sectionLines); + if (sectionName === 'SHAPES') shapeParser.finish(); + else if (sectionName) dispatchSection(seq, sectionName, sectionLines); sectionName = m[1]; seenSections.add(sectionName); sectionLines = []; + } else if (sectionName === 'SHAPES') { + shapeParser.consume(line); } else { sectionLines.push(line); } - } - if (sectionName) dispatchSection(seq, sectionName, sectionLines); + }); + if (sectionName === 'SHAPES') shapeParser.finish(); + else if (sectionName) dispatchSection(seq, sectionName, sectionLines); // Compute versionCombined AFTER parsing [VERSION] seq.versionCombined = makeVersionCombined( @@ -58,9 +62,25 @@ export function parseSequenceText(text: string): PulseqSequence { return seq; } +function forEachLine(text: string, visit: (line: string) => void): void { + let start = 0; + while (start <= text.length) { + let end = text.indexOf('\n', start); + if (end < 0) end = text.length; + const contentEnd = end > start && text.charCodeAt(end - 1) === 13 ? end - 1 : end; + visit(text.slice(start, contentEnd)); + if (end === text.length) break; + start = end + 1; + } +} + // ─── Section dispatcher ─────────────────────────────────────────────────── function dispatchSection(seq: PulseqSequence, name: string, lines: string[]): void { + if (name === 'SHAPES') { + parseShapes(seq, lines); + return; + } const valid = lines.filter(l => { const t = l.trim(); return t && !t.startsWith('#'); }); switch (name) { case 'VERSION': parseVersion(seq, valid); break; @@ -71,7 +91,6 @@ function dispatchSection(seq: PulseqSequence, name: string, lines: string[]): vo case 'TRAP': parseTrapGrads(seq, valid); break; case 'ADC': parseADC(seq, valid); break; case 'EXTENSIONS': parseExtensions(seq, valid); break; - case 'SHAPES': parseShapes(seq, lines); break; // [SIGNATURE] — intentionally ignored } } @@ -588,53 +607,82 @@ function parseRFShimSpecs(seq: PulseqSequence, lines: string[]): void { // ─── Shapes ─────────────────────────────────────────────────────────────── function parseShapes(seq: PulseqSequence, lines: string[]): void { - let i = 0; - while (i < lines.length) { - const t = lines[i].trim(); - if (!t || t.startsWith('#') || t.startsWith('[')) { i++; continue; } + const parser = new ShapeSectionParser(seq); + for (const line of lines) parser.consume(line); + parser.finish(); +} - const m = t.match(/^shape_id\s+(\d+)/); - if (!m) { i++; continue; } - const shapeId = +m[1]; - i++; +class ShapeSectionParser { + private shapeId = 0; + private numSamples = 0; + private raw = new Float64Array(); + private rawCount = 0; - // Find num_samples - let numSamples = 0; - while (i < lines.length) { - const l = lines[i].trim(); - if (!l || l.startsWith('#')) { i++; continue; } - const nm = l.match(/^num_samples\s+(\d+)/); - if (nm) { numSamples = +nm[1]; i++; break; } - if (l.match(/^shape_id\s+\d+/) || l.startsWith('[')) break; - i++; + constructor(private readonly seq: PulseqSequence) { } + + consume(line: string): void { + const trimmed = line.trim(); + if (!trimmed || trimmed.startsWith('#')) return; + + const shapeMatch = /^shape_id\s+(\d+)/.exec(trimmed); + if (shapeMatch) { + this.storeCurrent(); + this.shapeId = Number(shapeMatch[1]); + return; } - if (numSamples <= 0) continue; - - // Read sample values - const vals: number[] = []; - while (i < lines.length && vals.length < numSamples) { - const l = lines[i].trim(); - if (l.match(/^shape_id\s+\d+/) || l.startsWith('[')) break; - if (!l || l.startsWith('#')) { i++; continue; } - for (const n of l.split(/\s+/).map(Number).filter(x => !isNaN(x))) { - if (vals.length < numSamples) vals.push(n); - } - i++; + + const countMatch = /^num_samples\s+(\d+)/.exec(trimmed); + if (countMatch) { + this.numSamples = Number(countMatch[1]); + this.raw = new Float64Array(Math.min(this.numSamples, 1024)); + this.rawCount = 0; + return; } - if (vals.length === 0) continue; + if (this.shapeId <= 0 || this.numSamples <= 0 || this.rawCount >= this.numSamples) return; - storeShape(seq, shapeId, numSamples, vals); + if (!/\s/.test(trimmed)) { + this.appendRawValue(trimmed); + return; + } + for (const field of trimmed.split(/\s+/)) { + this.appendRawValue(field); + if (this.rawCount >= this.numSamples) break; + } + } + + finish(): void { + this.storeCurrent(); } -} -function storeShape(seq: PulseqSequence, id: number, num: number, raw: number[]): void { - // SeqEyes: decompress if run‑length encoded, otherwise use raw values as‑is. - // NO normalisation, NO clamping — amplitude shapes are already [0,1], - // time shapes are in grad‑raster units (can be large integers). - const decompressed = raw.length === num - ? new Float64Array(raw) - : decompressShape(raw, num); - seq.shapes.set(id, { numSamples: num, samples: decompressed }); + private ensureRawCapacity(): void { + if (this.rawCount < this.raw.length) return; + const nextLength = Math.min(this.numSamples, Math.max(1, this.raw.length * 2)); + const expanded = new Float64Array(nextLength); + expanded.set(this.raw); + this.raw = expanded; + } + + private appendRawValue(field: string): void { + const value = Number(field); + if (!Number.isFinite(value)) return; + this.ensureRawCapacity(); + this.raw[this.rawCount++] = value; + } + + private storeCurrent(): void { + if (this.shapeId > 0 && this.numSamples > 0 && this.rawCount > 0) { + // Uncompressed shapes reuse the preallocated buffer. Compressed shapes + // only expose their populated prefix to the existing decompressor. + const samples = this.rawCount === this.numSamples + ? this.raw + : decompressShape(this.raw.subarray(0, this.rawCount), this.numSamples); + this.seq.shapes.set(this.shapeId, { numSamples: this.numSamples, samples }); + } + this.shapeId = 0; + this.numSamples = 0; + this.raw = new Float64Array(); + this.rawCount = 0; + } } // ─── Raster times from definitions ────────────────────────────────────────