Open beorostica opened 6 months ago
curtcorum
commented
6 months ago It could possibly be good to leverage the tests for related KomaMRI or julia functions (read or generate seq or other files) from pulseq or pypulseq directly by sending KomaMRI input or output visa versa to those tests.
cncastillo
commented
5 months ago This was done in https://github.com/JuliaHealth/KomaMRI.jl/blob/master/KomaMRIFiles/test/runtests.jl#L53-L82
To include more tests, we just need to include Pulseq files in KomaMRIFiles/test/test_files/pulseq_read_comparison. We still need to add some examples for files generated with PyPulseq (#367).
To generate the corresponding .mat I used this script:
%% Just read a .seq and write the .mat
clc
cd('/home/ccp/Documents/KomaMRI.jl/KomaMRIFiles/test/test_files/pulseq_read_comparison/')
sys = mr.opts(); % Default system
seq = mr.Sequence(sys); % Default sequence
files = dir('*.seq');
for i = 1:length(files)
seq.read(files(i).name);
seq_name = split(files(i).name, '.');
seq.setDefinition('FileName', seq_name{1});
seq.plot() %<---- Modified version that saves the generated samples
endand the modified seq.plot function
diff --git a/matlab/+mr/@Sequence/Sequence.m b/matlab/+mr/@Sequence/Sequence.m
index acf3d18..e7b45b4 100644
--- a/matlab/+mr/@Sequence/Sequence.m
+++ b/matlab/+mr/@Sequence/Sequence.m
@@ -1255,6 +1255,17 @@ classdef Sequence < handle
end
end
%for iB=1:size(obj.blockEvents,1)
+ % Generating struct to save event samples to then compare them with Koma
+ koma_gamma = 42577468.8;
+ sequence = cell(1, length(obj.blockEvents));
+ for i = 1:length(obj.blockEvents)
+ sequence{i} = struct('rf', struct('t', [], 'A', []),...
+ 'df', struct('t', [], 'A', []),...
+ 'gx', struct('t', [], 'A', []),...
+ 'gy', struct('t', [], 'A', []),...
+ 'gz', struct('t', [], 'A', []),...
+ 'adc',struct('t', [], 'A', []));
+ end
for iB=1:length(obj.blockEvents)
block = obj.getBlock(iB);
isValid = t0+obj.blockDurations(iB)>timeRange(1) && t0<=timeRange(2);
@@ -1286,22 +1297,28 @@ classdef Sequence < handle
label_legend_2plot=[];
end
end
+ % Saving the ADC points to compare with Koma
+ sequence{iB}.adc.t = tFactor * (t0 + t);
+ sequence{iB}.adc.A = ones(size(t));
end
if ~isempty(block.rf)
rf=block.rf;
[tc,ic]=mr.calcRfCenter(rf);
sc=rf.signal(ic);
- if max(abs(diff(rf.t)-rf.t(2)+rf.t(1)))<1e-9 && length(rf.t)>100
- % homogeneous sampling and long pulses -- use lower time resolution for better display and performance
- dt=rf.t(2)-rf.t(1);
- st=round(obj.sys.gradRasterTime/dt);
- t=rf.t(1:st:end);
- s=rf.signal(1:st:end);
- else
+ %if max(abs(diff(rf.t)-rf.t(2)+rf.t(1)))<1e-9 && length(rf.t)>100
+ % % homogeneous sampling and long pulses -- use lower time resolution for better display and performance
+ % dt=rf.t(2)-rf.t(1);
+ % st=round(obj.sys.gradRasterTime/dt);
+ % t=rf.t(1:st:end);
+ % s=rf.signal(1:st:end);
+ %else
t=rf.t;
s=rf.signal;
- end
+ %end
sreal=max(abs(imag(s)))/max(abs(real(s)))<1e-6; %all(isreal(s));
+ % Saving raw waveform to compare with Koma
+ s_saved = [0; s; 0];
+ t_saved = [t(1); t; t(end)];
if abs(s(1))~=0 % fix strangely looking phase / amplitude in the beginning
s=[0; s];
t=[t(1); t];
@@ -1318,6 +1335,12 @@ classdef Sequence < handle
plot(tFactor*(t0+t+rf.delay), abs(s),'Parent',ax(2));
plot(tFactor*(t0+t+rf.delay), angle(s*exp(1i*rf.phaseOffset).*exp(1i*2*pi*t *rf.freqOffset)), tFactor*(t0+tc+rf.delay), angle(sc*exp(1i*rf.phaseOffset).*exp(1i*2*pi*tc*rf.freqOffset)),'xb', 'Parent',ax(3));
end
+ % Saving the RF points to compare with Koma
+ sequence{iB}.rf.t = tFactor * (t0 + t_saved + rf.delay);
+ sequence{iB}.rf.A = s_saved * exp(1i*rf.phaseOffset) / koma_gamma;
+ t_saved = [t(1); t(1); t(end); t(end)];
+ sequence{iB}.df.t = tFactor * (t0 + t_saved + rf.delay);
+ sequence{iB}.df.A = [0; rf.freqOffset; rf.freqOffset; 0];
%plot(tFactor*(t0+t+rf.delay), angle( exp(1i*rf.phaseOffset).*exp(1i*2*pi*t *rf.freqOffset)), tFactor*(t0+tc+rf.delay), angle( exp(1i*rf.phaseOffset).*exp(1i*2*pi*t(ic)*rf.freqOffset)),'xb', 'Parent',ax(3));
end
gradChannels={'gx','gy','gz'};
@@ -1331,9 +1354,29 @@ classdef Sequence < handle
%t=grad.delay + [0; grad.t + (grad.t(2)-grad.t(1))/2; grad.t(end) + grad.t(2)-grad.t(1)];
t= grad.delay+[0; grad.tt; grad.shape_dur];
waveform=1e-3* [grad.first; grad.waveform; grad.last];
+ % Saving waveform to compare with Koma
+ t_saved = t;
+ wave_saved = [grad.first; grad.waveform; grad.last];
else
t=cumsum([0 grad.delay grad.riseTime grad.flatTime grad.fallTime]);
waveform=1e-3*grad.amplitude*[0 0 1 1 0];
+ % Saving waveform to compare with Koma
+ t_saved = t(2:end);
+ wave_saved = grad.amplitude*[0 1 1 0];
+ end
+ if sum(abs(waveform)) ~= 0
+ gi_t = tFactor * (t0 + t_saved);
+ gi_A = wave_saved / koma_gamma;
+ if j == 1
+ sequence{iB}.gx.t = gi_t;
+ sequence{iB}.gx.A = gi_A;
+ elseif j == 2
+ sequence{iB}.gy.t = gi_t;
+ sequence{iB}.gy.A = gi_A;
+ elseif j == 3
+ sequence{iB}.gz.t = gi_t;
+ sequence{iB}.gz.A = gi_A;
+ end
end
plot(tFactor*(t0+t),waveform,'Parent',ax(3+j));
end
@@ -1342,6 +1385,9 @@ classdef Sequence < handle
t0=t0+obj.blockDurations(iB);%mr.calcDuration(block);
end
+ filename = strcat('/home/ccp/Documents/KomaMRI.jl/KomaMRIFiles/test/test_files/pulseq_read_comparison/', obj.definitions('FileName'), '.mat');
+ display(filename)
+ save(filename, 'sequence')
% Set axis limits and zoom properties
dispRange = tFactor*[timeRange(1) min(timeRange(2),t0)];
arrayfun(@(x)xlim(x,dispRange),ax);
It is a good idea to go deeper in the creation of code tests. This really will help to accept the contribution of other guys. For instance, tests for KomaMRIFiles for pulseq files .seq.