Incorrect gradient interpretation when waveforms do not end in zero

[beorostica]

This pull request is meant to address #320

[cncastillo]

Ok, the problem was that we needed to implement the heuristic that Pulseq uses for the first and last points. As these are not stored in the Pulseq file they need to be "invented" while reading.

We need to implement this part:

gradChannels={'gx','gy','gz'};
gradPrevLast=zeros(1,length(gradChannels));
for iB = 1:length(obj.blockEvents)
    b=obj.getBlock(iB);
    block_duration=obj.blockDurations(iB);
    % we also need to keep track of the event IDs because some Pulseq files written by external software may contain repeated entries so searching by content will fail 
    eventIDs=obj.blockEvents{iB};
    % update the objects by filling in the fields not contained in the
    % pulseq file
    for j=1:length(gradChannels)
        grad=b.(gradChannels{j});
        if isempty(grad)
            gradPrevLast(j)=0;
            continue;
        end
        if strcmp(grad.type,'grad')
            if grad.delay>0 
                gradPrevLast(j)=0;
            end
            if isfield(grad,'first')
                continue;
            end
            grad.first = gradPrevLast(j);
            % is this an extended trapezoid?
            if grad.time_id~=0
                grad.last=grad.waveform(end);
                grad_duration=grad.delay+grad.tt(end);
            else
                % restore samples on the edges of the gradient raster intervals
                % for that we need the first sample 
                odd_step1=[grad.first 2*grad.waveform'];
                odd_step2=odd_step1.*(mod(1:length(odd_step1),2)*2-1);
                waveform_odd_rest=(cumsum(odd_step2).*(mod(1:length(odd_step2),2)*2-1))';
                grad.last = waveform_odd_rest(end);
                grad_duration=grad.delay+length(grad.waveform)*obj.gradRasterTime;
            end
            % bookkeeping
            gradPrevLast(j) = grad.last;
            if grad_duration+eps<block_duration
                gradPrevLast(j)=0;
            end
            % need to recover the amplidute from the library data directly...
            id=eventIDs(j+2);
            amplitude=obj.gradLibrary.data(id).array(1);
            %
            if version_combined>=1004000
                old_data = [amplitude grad.shape_id grad.time_id grad.delay];
            else
                old_data = [amplitude grad.shape_id grad.delay];
            end
            new_data = [amplitude grad.shape_id grad.time_id grad.delay grad.first grad.last];
            update_data(obj.gradLibrary, id, old_data, new_data,'g');
        else
            gradPrevLast(j)=0;
        end
    end
end

As this is arbitrary, we could preserve the current method and create a new one that works because Pulseq's implementation does not fully restore the continuity of the samples.

[cncastillo]

Please fix.

[cncastillo]

Hi @beorostica, please fix what I pointed out in the comments before keeping the work with the PR.

Specially

• Incorrect gradient interpretation when waveforms do not end in zero #321 (comment)
• Incorrect gradient interpretation when waveforms do not end in zero #321 (comment)
• Incorrect gradient interpretation when waveforms do not end in zero #321 (comment)

Another checklist below #321 (comment)

[cncastillo]

Just to help, to dispatch based on the Grad type we need to redefine it like so:

using Parameters
abstract type MRISequenceEvent end
NumOrVect = Union{T, Vector{T}} where {T<:Number}

@with_kw_noshow mutable struct Grad{GradAmp<:NumOrVect, TimeDur<:NumOrVect} <: MRISequenceEvent
    # Grad properties
    A    ::GradAmp = zero(eltype(GradAmp))
    T    ::TimeDur = zero(eltype(TimeDur))
    rise ::eltype(TimeDur) = zero(eltype(T))
    fall ::eltype(TimeDur) = rise
    delay::eltype(TimeDur) = zero(eltype(T))
    first::eltype(GradAmp) = zero(eltype(A))
    last ::eltype(GradAmp) = first
    # Checking validity of Grad
    @assert length(A) >= length(T)
    @assert all(T .>= zero(eltype(T)))
    @assert rise   >= zero(eltype(T))
    @assert fall   >= zero(eltype(T))
    @assert delay  >= zero(eltype(T))
end

# Gradient type aliases for dispatching
TrapezoidalGrad      = Grad{G, T} where {G, T}
UniformlySampledGrad = Grad{Vector{G}, T} where {G, T}
TimeShapedGrad       = Grad{Vector{G}, Vector{T}} where {G, T}
ValidGrad = Union{TrapezoidalGrad{G, T}, UniformlySampledGrad{G, T}, TimeShapedGrad{G, T}} where {G, T}

# MRISequenceEvent functions
# Event duration
dur(g::Grad) = g.delay + g.rise + sum(g.T) + g.fall
dur(G::Vector{Grad{G, T}}) where {G, T} = max(dur.(G)...)
# Event amplitude
ampl(g::TrapezoidalGrad{G,T})      where {G, T} = vcat(zero(G), g.first, g.A, g.A, g.last)
ampl(g::TimeShapedGrad{G,T})       where {G, T} = vcat(zero(G), g.first, g.A, g.last)
ampl(g::UniformlySampledGrad{G,T}) where {G, T} = vcat(zero(G), g.first, g.A, g.last)
# Event time
time(g::TrapezoidalGrad{G,T})      where {G, T} = cumsum(vcat(zero(T), g.delay, g.rise, g.T, g.fall))
time(g::TimeShapedGrad{G,T})       where {G, T} = cumsum(vcat(zero(T), g.delay, g.rise, g.T, g.fall))
time(g::UniformlySampledGrad{G,T}) where {G, T} = begin
    NA = length(g.A) - 1
    ΔT = repeat([g.T / NA], NA) 
    return cumsum(vcat(zero(T), g.delay, g.rise, ΔT, g.fall))
end

This also enables the use Unitful.jl, like g1 = Grad(A=[1,1]*1mT/m, T=1ms, rise=0.5ms, delay=1ms) for example (related to #305 and #139, specifically #4 (comment)). Then to plot the gradient we can use:

using Plots
plot(time(g1), ampl(g1))