Merge develop branch into main (#16)

* Ignore blocks containing only labels
On branch develop
	modified:   isdelayblock.m
	modified:   seq2ceq.m

* Wrapper for seq2ceq
On branch develop
	new file:   seq2seg.m

* Squashed commit of the following:

    Support for multiple ADC events

    Add number of ADC samples to .mod header
    Store dwell time (in us) in .mod header int array
    	modified:   ceq2ge.m

    Allow multiple ADC events
    	modified:   compareblocks.m

On branch develop
	modified:   matlab/ceq2ge.m
	modified:   matlab/compareblocks.m

* Rename 'group' to 'segment'
On branch develop
	modified:   pulCeq.h

* Working on alternate figure explaining sequence model
On branch develop
	new file:   model2.svg

* On branch develop
	modified:   model2.svg

* On branch develop
	modified:   model2.svg

* Add physio trigger (#15)

Add cardiac/physio trigger to ceq.loop array, and write it to scanloop.txt.

On branch develop_trig
	modified:   ceq2ge.m
	modified:   getdynamics.m
	modified:   seq2ceq.m
	new file:   getblocktype.m

Co-authored-by: Jon-Fredrik Nielsen <[email protected]>

* Bug fix
On branch develop
	modified:   ceq2ge.m

---------

Co-authored-by: Jon-Fredrik Nielsen <[email protected]>

matlab/ceq2ge.m

@@ -111,6 +111,11 @@ function ceq2ge(ceq, sysGE, ofname, varargin)
         npre = round(b.adc.delay/raster);
         rfres = round(b.adc.numSamples*b.adc.dwell/raster);
         readoutFile = modFiles{p};
+        dwell = round(b.adc.dwell*1e6);
+        numSamples = b.adc.numSamples;
+    else
+        dwell = [];
+        numSamples = [];
     end
 
     % Set nChop, which is the number of samples to trim from beginning and end of RF/ADC window
@@ -125,7 +130,7 @@ function ceq2ge(ceq, sysGE, ofname, varargin)
     if ~isDelayBlock 
         toppe.writemod(sysGE, 'ofname', modFiles{p}, ...
             'rf', rf(:), 'gx', grad.x(:), 'gy', grad.y(:), 'gz', grad.z(:), ...
-            'nChop', nChop);
+            'nChop', nChop, 'hdrints', [dwell numSamples]);
     end
 end
 
@@ -144,15 +149,16 @@ function ceq2ge(ceq, sysGE, ofname, varargin)
     end
 
     % trigger out
+    trigpos = -1;  % default: no trigger
     if isfield(b, 'trig') & ~arg.ignoreTrigger
-        if b.trig.delay + eps < 100e-6
-            warning('Requested trigger time too short. Setting to 100us');
-            trigpos = 100;  % us
-        else
-            trigpos = round(b.trig.delay*1e6);   % us
+        if ~strcmp(b.trig.channel, 'physio1')
+            if b.trig.delay + eps < 100e-6
+                warning('Requested trigger time too short. Setting to 100us');
+                trigpos = 100;  % us
+            else
+                trigpos = round(b.trig.delay*1e6);   % us
+            end
         end
-    else
-        trigpos = -1;    % no trigger
     end
 
     rf = toppe.readmod(modFiles{p});
@@ -259,6 +265,7 @@ function ceq2ge(ceq, sysGE, ofname, varargin)
         'textra',      0, ...  
         'waveform',    1, ...
         'trigout',     trigout, ...
+        'trig',        ceq.loop(n, 11), ...
         'core', i);
 
 end

matlab/compareblocks.m

@@ -113,12 +113,13 @@
         issame = false; return;
     end
 
-    if (adc1.numSamples == adc2.numSamples | ...
-        adc1.dwell == adc2.dwell | ...
+    if (adc1.numSamples == adc2.numSamples & ...
+        adc1.dwell == adc2.dwell & ...
         adc1.delay == adc2.delay)
         issame = true;
     else
-        error(sprintf('blocks %d, %d: ADC numsamples/dwell/delay must be same for all ADC events', n1, n2));
+        issame = false;
+        %error(sprintf('blocks %d, %d: ADC numsamples/dwell/delay must be same for all ADC events', n1, n2));
     end
 
 return

matlab/getblocktype.m

@@ -0,0 +1,41 @@
+function T = getblocktype(block)
+% function T = getblocktype(block)
+%
+% Returns a vector with flags indicating block type/properties.
+%
+%  T = [<has rf/gradient/adc> ...
+%       <has TRID label> ...
+%       <has physio1 trigger> ...
+%       <is pure delay block>]
+
+minDelayBlockDuration = 20e-6 + 100*eps; 
+
+% Are waveforms and/or ADC events defined?
+T(1) = isempty(block.rf) & isempty(block.adc) & ...
+       isempty(block.gx) & isempty(block.gy) & isempty(block.gz);
+T(1) = ~T(1);
+
+% Does the block define a TRID (start of segment) label?
+T(2) = false;
+if isfield(block, 'label')
+    for ii = 1:length(block.label)
+        if strcmp(block.label(ii).label, 'TRID')
+            T(2) = true;
+            break;
+        end
+    end
+end
+
+% Does the block contain a cardiac trigger?
+T(3) = false;
+if isfield(block, 'trig')
+    if strcmp(block.trig.channel, 'physio1')
+        T(3) = true;
+    end
+end
+
+% Is this a pure delay block?
+T(4) = ~T(1) & block.blockDuration >= minDelayBlockDuration;
+
+%assert( ~ (T(1) | block.blockDuration < minDelayBlockDuration), ...
+%       sprintf('Duration of block containing rf/gradient/adc events must exceed %es', minDelayBlockDuration) );

matlab/getdynamics.m

@@ -1,7 +1,11 @@
-function loop = getdynamics(block, segmentID, parentBlockID, parentBlock)
+function loop = getdynamics(block, segmentID, parentBlockID, parentBlock, physioTrigger)
 % Return vector containing waveform amplitudes, RF/ADC phase, etc,
 % for a Pulseq block, in physical (Pulseq) units.
 
+if nargin < 5
+    physioTrigger = false;
+end
+
 % defaults
 rfamp = 0;
 rfphs = 0;
@@ -42,5 +46,5 @@
     recphs = block.adc.phaseOffset;
 end
 
-loop = [segmentID parentBlockID rfamp rfphs rffreq amp.gx amp.gy amp.gz recphs block.blockDuration];
+loop = [segmentID parentBlockID rfamp rfphs rffreq amp.gx amp.gy amp.gz recphs block.blockDuration physioTrigger];
 

matlab/isdelayblock.m

@@ -3,5 +3,5 @@
 % 
 % Returns true if block contains no waveforms
 
-idb = isempty(block.rf) & isempty(block.adc) & ...
+idb = block.blockDuration > 0 & isempty(block.rf) & isempty(block.adc) & ...
         isempty(block.gx) & isempty(block.gy) & isempty(block.gz);

matlab/seq2ceq.m

@@ -9,7 +9,7 @@
 %
 % Input options with defaults
 %   nMax                  [1]            Only parse the first nMax blocks in the .seq file  [all blocks]
-%   ignoreSegmentLables   true/false     Treat each block as a segment. Use with caution! [false]
+%   ignoreSegmentLabels   true/false     Treat each block as a segment. Use with caution! [false]
 %   verbose               true/false     Print some info to the terminal [false]
 %
 % Output
@@ -54,11 +54,13 @@
 % parent blocks = unique up to a scaling factor, or phase/frequency offsets.
 % Contains waveforms with maximum amplitude across blocks.
 % First find unique blocks, then determine and set max amplitudes.
-% parentBlockIDs = [1 nMax], vector of parent block IDs for all blocks
+% parentBlockIDs = [nMax], vector of parent block IDs for all blocks
 
-parentBlockIndex(1) = 1;  % first block is unique by definition
+parentBlockIndex = []; 
 
-fprintf('seq2ceq: Getting block %d/%d', 1, ceq.nMax); prev_n = 1; % Progress update trackers
+parentBlockIDs = [];
+
+fprintf('seq2ceq: Getting block %d/%d', 1, ceq.nMax); prev_n = 1;
 for n = 1:ceq.nMax
     if ~mod(n, 500) || n == ceq.nMax
         for ib = 1:strlength(sprintf('seq2ceq: Getting block %d/%d', prev_n, ceq.nMax))
@@ -69,20 +71,31 @@
         if n == ceq.nMax, fprintf('\n'), end
     end
 
-    % Pure delay blocks are handled separately
     b = seq.getBlock(n);
+
+    % Skip blocks with zero duration (only contains label(s))
+    if b.blockDuration < 2*eps
+        parentBlockIDs(n) = -1;
+        continue;
+    end
+
+    % Pure delay blocks are handled separately
     if isdelayblock(b)
         parentBlockIDs(n) = 0;
         continue;
     end
 
+    if isempty(parentBlockIndex)
+        parentBlockIndex(1) = n;
+    end
+
     for p = 1:length(parentBlockIndex)
         n2 = parentBlockIndex(p);
         IsSame(p) = compareblocks(seq, blockEvents(n,:), blockEvents(n2,:), n, n2);
     end
     if sum(IsSame) == 0
         if arg.verbose
-            fprintf('\nFound new block on line %d\n', n);
+            fprintf('\nFound new parent block on line %d\n', n);
         end
         parentBlockIndex(p+1) = n;  % add new block to list
         parentBlockIDs(n) = p+1;
@@ -105,7 +118,7 @@
     ceq.parentBlocks{p}.amp.gz = 0;
 
     for n = 1:ceq.nMax
-        if parentBlockIDs(n) ~= p
+        if parentBlockIDs(n) ~= p | parentBlockIDs(n) == -1
             continue; 
         end
         block = seq.getBlock(n);
@@ -148,16 +161,28 @@
 
 
 %% Get segment (block group) definitions
-% Segments defined by their first occurrence in the .seq file.
+% Segments are defined by their first occurrence in the .seq file
 previouslyDefinedSegmentIDs = [];
 if ~arg.ignoreSegmentLabels
     segmentIDs = zeros(1,ceq.nMax);  % keep track of which segment each block belongs to
     for n = 1:ceq.nMax
         b = seq.getBlock(n);
 
+        if parentBlockIDs(n) == -1
+            continue;    % skip
+        end
+
+        % Get segment ID label (TRID) if present
         if isfield(b, 'label') 
-            if strcmp(b.label.label, 'TRID')   % marks start of segment
-                activeSegmentID = b.label.value;
+            hasTRIDlabel = false;
+            for ii = 1:length(b.label)
+                if strcmp(b.label(ii).label, 'TRID')
+                    hasTRIDlabel = true;
+                    break;
+                end
+            end
+            if hasTRIDlabel  % marks start of segment
+                activeSegmentID = b.label(ii).value;
 
                 if ~any(activeSegmentID == previouslyDefinedSegmentIDs)
                     % start new segment
@@ -171,7 +196,7 @@
         end
 
         if ~exist('firstOccurrence', 'var')
-            error('First block must contain a segment ID');
+            %error('First block must contain a segment ID');
         end
 
         % add block to segment
@@ -220,37 +245,47 @@
 ceq.nGroups = length(ceq.groups);
 
 %% Get dynamic scan information
-ceq.loop = zeros(ceq.nMax, 10);
+ceq.loop = zeros(ceq.nMax, 11);
+physioTrigger = false;
 for n = 1:ceq.nMax
     b = seq.getBlock(n);
+    % set cardiac trigger
+    T = getblocktype(b);
+    physioTrigger = T(3);
     p = parentBlockIDs(n); 
     if p == 0  % delay block
         ceq.loop(n,:) = getdynamics(b, segmentID2Ind(segmentIDs(n)), p);
-    else
-        ceq.loop(n,:) = getdynamics(b, segmentID2Ind(segmentIDs(n)), p, ceq.parentBlocks{p});
+    elseif p > 0
+        ceq.loop(n,:) = getdynamics(b, segmentID2Ind(segmentIDs(n)), p, ceq.parentBlocks{p}, physioTrigger);
     end
 end
 
+
+%% Remove zero-duration (label-only) blocks from ceq.loop
+ceq.loop(parentBlockIDs == -1, :) = [];
+ceq.nMax = size(ceq.loop,1);
+
+
 %% Check that the execution of blocks throughout the sequence
 %% is consistent with the segment definitions
 n = 1;
 while n < ceq.nMax
     i = ceq.loop(n, 1);  % segment id
+
+    if (n + ceq.groups(i).nBlocksInGroup) > ceq.nMax
+        break;
+    end
+
+    % loop through blocks in segment
     for j = 1:ceq.groups(i).nBlocksInGroup
+
+        % compare parent block id in ceq.loop against block id in ceq.groups(i)
         p = ceq.loop(n, 2);  % parent block id
         p_ij = ceq.groups(i).blockIDs(j);
         if p ~= p_ij
-            %b = ceq.parentBlocks{p};
-            %b_ij = ceq.parentBlocks{p_ij};
             warning(sprintf('Sequence contains inconsistent segment definitions. This may occur due to programming error (possibly fatal), or if an arbitrary gradient resembles that from another block except with opposite sign or scaled by zero (which is probably ok). Expected parent block ID %d, found %d (block %d)', p_ij, p, n));
         end
+
         n = n + 1;
     end
 end
-
-%fprintf('\n');
-
-% function returns true if blocks are same except for zero-scaling of RF and/or gradients 
-function sub_compareblocks(b1, b2)
-
-return

matlab/seq2seg.m

@@ -0,0 +1,4 @@
+% Call 'seq2ceq' by the name 'seq2seg' 
+function seg = seq2seg(seqarg, varargin)
+
+seg = seq2ceq(seqarg, varargin{:});

src/pulCeq.h

@@ -57,34 +57,34 @@ typedef struct {
     PulseqADC  adc;
     PulseqTrig trig;
 
-    /* vectors available for use as needed by the client program */
-    int nVal1;                /* number of int values. Must be defined. */
-    int* val1;                /* to be allocated dynamically by the client program */
-    int nVal2;                /* number of float values. Must be defined. */
-    float* nVal2;             /* to be allocated dynamically by the client program */
+    /* arrays available for use as needed by the client program */
+    int nVal1;            /* number of int values. Must be defined. */
+    int* val1;            /* to be allocated dynamically by the client program */
+    int nVal2;            /* number of float values. Must be defined. */
+    float* nVal2;         /* to be allocated dynamically by the client program */
 } PulseqBlock; 
 
-/* Struct containing block IDs for all block groups */
+/* Struct containing block IDs for all segments */
 typedef struct {
-    int groupID;
-    int  nBlocksInGroup;    /* number of blocks in group */
-    int* blockIDs;          /* block id's in this group */
-} BlockGroup;
+    int  segmentID;
+    int  nBlocks;        /* number of blocks in segment */
+    int* blockIDs;       /* block id's in this segment */
+} Segment;
 
 /* struct containing entire sequence definition */
 typedef struct {
     int nParentBlocks; 
     PulseqBlock* parentBlocks;
 
-    int nGroups;    
-    BlockGroup* groups;      /* optional */
+    int nSegments;    
+    Segment* segments;      /* optional */
 
     float** loop    /* Dynamic scan settings: waveform amplitudes, phase offsets, etc.
-                       loop[n] = [groupID blockID rfamp gxamp gyamp gzamp rfFreqOffset rfPhaseOffset ...]
-                       units:    [int     int     T     mT/m  mT/m  mT/m  Hz           rad           ...]
+                       loop[n] = [segmentID blockID rfamp gxamp gyamp gzamp rfFreqOffset rfPhaseOffset ...]
+                       units:    [int       int     T     mT/m  mT/m  mT/m  Hz           rad           ...]
                     */
 
-   int nMax;         /* number of blocks (rows in BLOCKS section) in .seq file */
+    int nMax;         /* number of blocks (rows in BLOCKS section) in .seq file */
 } Ceq;
 
 /* function prototypes that this specification implements */