KomaMRIFiles: Include .seq produced by pypulseq or pulseq

[beorostica]

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.

• Test files generated with MATLAB Pulseq (1.4)
• Test files generated with Python PyPulseq (1.3) #367

[curtcorum]

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]

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
end

and 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);

Modified_plot_seq_Pulseq.zip