Readout Gradient

In this final part, we will add gradient pulses and get a functional FLASH sequence. The first pulse we are going to add is provided by the msl::Sampling object. As for RF pulses, the msl::Sampling contains a set of real-time events, and can be played in the run function of the kernel. As for the slices and the phase cycling, we will start by passing the sampling object, and its accompanying trajectory counter, to the kernel:

// ...
#include <msl/samplings.h>
// ...
class Kernel: public msl::graph::AbstractNode
{
public:
    DECLARE_POINTERS(Kernel);
    static Pointer New(
        std::string const & sampling, std::string const & trajectories,
        std::string const & slices, std::string const & phaseCycling);
    // ...
private:
    msl::DictionaryItem<msl::Sampling::Pointer> _sampling;
    msl::DictionaryItem<msl::Counter> _trajectories;
    // ...
    
    Kernel(
        std::string const & sampling, std::string const & trajectories,
        std::string const & slices, std::string const & phaseCycling);
};

// ...
Kernel::Pointer
Kernel
::New(
    std::string const & sampling, std::string const & trajectories,
    std::string const & slices, std::string const & phaseCycling)
{
    return Pointer(new Kernel(sampling, trajectories, slices, phaseCycling));
}
// ...
Kernel
::Kernel(
    std::string const & sampling, std::string const & trajectories,
    std::string const & slices, std::string const & phaseCycling)
: AbstractNode(), _sampling(this, sampling), _trajectories(this, trajectories),
    _slices(this, slices), _phaseCycling(this, phaseCycling)
{
    // Nothing else
}

NLSStatus
FLASH
::initialize(SeqLim & limits)
{
    // ...
    auto trajectories = msl::graph::Loop::New(
        "trajectories", {
            Kernel::New("sampling", "trajectories", "slices", "phaseCycling"),
            msl::graph::Action::New([](msl::Dictionary::Pointer d) {
                ++d->get<msl::PhaseCycling>("phaseCycling"); })
        }
    );
    // ...
}

We then need to prepare the sampling: this requires setting the phase equal to that of the RF pulse (to correctly demodulate the signal during the readout), and the slice. In addition, we will specify the maximal amplitude and slew rate of the gradient using a msl::GradientSpecs object, as well as the time at which the excitation pulse occurs. This last field is required to correctly determine the time at which the center of the k-space is reached. As for the RF pulse, once the sampling is configured, we call its prepare function, wrapped in the error-handling ON_ERROR_RETURN_STATUS macro.

// ...
#include <msl/GradientSpecs.h>
// ...
 
NLSStatus
Kernel
::prepare(MrProt & protocol, SeqLim & limits, MrProtocolData::SeqExpo & exports)
{
    // ...
    auto const gradientSpecs = msl::GradientSpecs::current(protocol);
    dynamic_cast<msl::samplings::Cartesian &>(*this->_sampling())
        .setPhase(rfPhase)
        .setTimeOffset(this->_excitation.peakTime())
        .setGradientSpecs(gradientSpecs)
        .setSlice(slice);
    ON_ERROR_RETURN_STATUS(this->_sampling()->prepare(protocol, limits, exports));
    // ...
}

Finally, we run the sampling object in the run function of the kernel. As for the RF pulse, it will require the phase of the RF pulse and the slice. In addition, it also requires the index of the current trajectory, and an update of the readout based on the current index.

NLSStatus
Kernel
::run(MrProt & protocol, SeqLim & limits, MrProtocolData::SeqExpo & exports)
{
    auto const & slice = this->_slices().item();
    auto const index = this->_trajectories().index();
    
    // ...
    (*this->_sampling()).setPhase(rfPhase).setIndex(index).setSlice(slice);
    ON_ERROR_RETURN_STATUS(
        this->_sampling()->updateReadout(protocol, limits, exports));
    
    // ...
    ON_ERROR_RETURN_STATUS(this->_excitation.run(protocol, limits, exports));
    ON_ERROR_RETURN_STATUS(this->_sampling()->run(protocol, limits, exports));
    // ...
}

After simulation, you will see the readout gradient being played, and the synchronized ADC object, also provided by the sampling object. Note that for Cartesian sampling, all readouts are the same: the displacement in the k-space happening during the readout is constant, only the displacement before the readout, to reach the start of the current line will vary. Calling msl::Sampling::setIndex and msl::Sampling::updateReadout does nothing for this specific samping, but it is a good practice to call these functions regardless of the actual sampling type.

Moving through the k-Space

Right after the excitation pulse, the k-space coordinates are 0 for \(k_x\) and \(k_y\), but non-zero for \(k_z\) due to the part of the slice selection gradient which occurs after the peak of the excitation pulse. For our Cartesian FLASH sequence, we must be at a specific point in the k-space before starting the readout-out:

on the \(k_x\) axis, the center of the readout must be at the center of the k-space: the start of the readout must be at the "left edge" of the k-space (if the readout amplitude is positive, otherwise at the "right edge")
the \(k_y\) and \(k_z\) coordinates depend on the current line, and, for a 3D sequence, the current partition.

Moving to the start of the trajectory has then a fixed part on the \(k_x\) axis and on the \(k_z\) axis, and a part depending on the current trajectory on the \(k_y\) and on the \(k_z\) axes.

In msl, each sampling object generate msl::Readout objects, which contain an ADC object (using the raw sREADOUT object from IDEA) for signal recording, and a msl::GradientPulse for the accompanying gradient. The gradient area incurred before the echo is returned by the msl::Readout::areaBeforeEcho function: it is a 3D vector, i.e. a msl::Vector3d. If the echo must be at the center of the readout, then the motion on the \(k_x\) axis must be equal to -readout.areaBeforeEcho().

The msl::GradientPulse objects can also report area between specific time points using the msl::GradientPulse::area, msl::GradientPulse::areaFrom, and msl::GradientPulse::areaTo functions. Real-time objects, as RF pulses and gradient pulses provide timing informations throught the msl::RealTimeEvents::startTime and msl::RealTimeEvents::startTime functions. Using these functions, we can compute the motion needed to compensate the slice-selection gradient occuring after the pulse peak:

- this->_excitation.gradient().areaFrom(this->_excitation.peakTime() - this->_excitation.startTime());

The addition of the fixed part of the move-to-begin gradient can then be added to the kernel:

// ...
#include <msl/Vector.h>
 
class Kernel: public msl::graph::AbstractNode
{
    // ...
private:
    // ...
    msl::Vector3d _toBeginFixed;
};

NLSStatus
Kernel
::prepare(MrProt & protocol, SeqLim & limits, MrProtocolData::SeqExpo & exports)
{
    // ...
    ON_ERROR_RETURN_STATUS(this->_sampling()->prepare(protocol, limits, exports));
    
    // Constant part of move-to-begin
    auto const & readout = this->_sampling()->readout();
    this->_toBeginFixed =
        // Move to begin on the readout axis
        - readout.areaBeforeEcho()
        // Compensate the slice-selection gradient offset
        - this->_excitation.gradient().areaFrom(
            this->_excitation.peakTime() - this->_excitation.startTime());
    // ...
}

The varying part of the to-begin motion depends on the current line and partition (for 3D sequences): the gradient area will be proportional to the line index, and, for a given maximal amplitude, will require a longer time for lines and partitions closer to the edge of the k-space. This means that, if we can fit the gradient lobe corresponding to extremal lines in the available time, then any other line will also fit in the same constraints. The current k-space information is available in the msl::KSpace class: we can add the following the the prepare function of the kernel.

#include <msl/KSpace.h>
 
NLSStatus
Kernel
::prepare(MrProt & protocol, SeqLim & limits, MrProtocolData::SeqExpo & exports)
{
    // ...
    
    // Only prepare the most extreme k-space line: if this works, then any other
    // will work. Due to the way k-space is defined, this is always based on
    // nMin
    msl::KSpace const kSpace(protocol);
    auto const extremalKYZ = kSpace.k({0, kSpace.nMin()[1], kSpace.nMin()[2]});
    
    // Prepare with a gradient area of this->_toBeginFixed + extremalKYZ
    // ...
}

To actually run the gradient pulse, we will need to add a msl::GradientPulse object to the kernel class:

// ...
#include <msl/GradientPulse.h>
// ...
 
class Kernel: public msl::graph::AbstractNode
{
    // ...
private:
    // ...
    msl::Vector3d _toBeginFixed;
    msl::GradientPulse _toBegin;
    // ...
};

We will then prepare this gradient pulse to be the quickest while accomodating the larges k-space displacement possible, and set its start time so that it ends when the readout starts:

NLSStatus
Kernel
::prepare(MrProt & protocol, SeqLim & limits, MrProtocolData::SeqExpo & exports)
{
    // ...
    auto const extremalKYZ = kSpace.k({0, kSpace.nMin()[1], kSpace.nMin()[2]});
    this->_toBegin = msl::gradient_pulses::quickest(
        this->_toBeginFixed + extremalKYZ, gradientSpecs);
    this->_toBegin.setEndTime(readout.startTime());
}

In the run method, we will then get the index of the current point and convert it to a k-space vector:

NLSStatus
Kernel
::run(MrProt & protocol, SeqLim & limits, MrProtocolData::SeqExpo & exports)
{
    // ...
    auto const point = 
        dynamic_cast<msl::samplings::Cartesian &>(*this->_sampling())
        .enabledPoints()[index];
    
    msl::KSpace const kSpace(protocol);
    auto const kPoint = kSpace.k(
        msl::Vector3l{kSpace.nCenter()[0], point[0], point[1]}
        - kSpace.nCenter());
    // ...
}

The value of kPoint needs to be added to _toBeginFixed, and we can update the _toBegin gradient pulse, while keeping its duration.

NLSStatus
Kernel
::run(MrProt & protocol, SeqLim & limits, MrProtocolData::SeqExpo & exports)
{
    // ...
    auto const gradientSpecs = msl::GradientSpecs::current(protocol);
    this->_toBegin.update(msl::gradient_pulses::fixedDuration(
        this->_toBeginFixed + kPoint, this->_toBegin, gradientSpecs));
    // ...
}

The last step is to run this gradient lobe in the real-time events block:

NLSStatus
Kernel
::run(MrProt & protocol, SeqLim & limits, MrProtocolData::SeqExpo & exports)
{
    // ...
    ON_ERROR_RETURN_STATUS(this->_excitation.run(protocol, limits, exports));
    ON_ERROR_RETURN_STATUS(this->_toBegin.run());
    ON_ERROR_RETURN_STATUS(this->_sampling()->run(protocol, limits, exports));
    // ...
}

Gradient Spoiling

Similarly to the to-begin gradient, we need another gradient lobe after the echo: it must go back to the center of the k-space on the \(k_y\) and \(k_z\) axes in order to respect the steady-state conditions, and it must apply spoiling on the \(k_x\) axis.

The spoiling part is fixed and the back-to-center part varies according to lines and partitions, so we will use a similar design as before: add a msl::Vector3d for the spoiler and a msl::GradientPulse for both fixed and varying parts.

class Kernel: public msl::graph::AbstractNode
{
    // ...
private:
    // ...
    // Return to the center of the k-space (y and z) and spoil (x)
    msl::Vector3d _spoiler;
    msl::GradientPulse _toCenterAndSpoil;
    // ...
};

In the prepare function of the kernel, we define the spoiler based on the readout, and the gradient pulse to be the shortest based on extremal k-space lines.

NLSStatus
Kernel
::prepare(MrProt & protocol, SeqLim & limits, MrProtocolData::SeqExpo & exports)
{
    // ...
    
    // Move back to k-space center (y and z) and spoil (x)
    this->_spoiler = readout.gradient().area()/2.;
    this->_toCenterAndSpoil = msl::gradient_pulses::quickest(
        this->_spoiler - extremalKYZ, gradientSpecs);
    // Merge the ramp down of the readout and the spoiler
    this->_toCenterAndSpoil.setStartTime(
        readout.endTime() - readout.gradient().rampDuration());
    
    // ...
}

In the run function, we update the gradient with respect to the current line, and run it in the real-time events block.

NLSStatus
Kernel
::run(MrProt & protocol, SeqLim & limits, MrProtocolData::SeqExpo & exports)
{
    // ...
    
    this->_toCenterAndSpoil.update(msl::gradient_pulses::fixedDuration(
        this->_spoiler - kPoint, this->_toCenterAndSpoil, gradientSpecs));
    
    // ...
    ON_ERROR_RETURN_STATUS(this->_toCenterAndSpoil.run());
    // ...
}

Checking the Kernel Timing

The RF and gradient pulses are now complete for the FLASH sequence. The kernel is however not yet robust: if the user a very low echo time or a very high resolution, the kernel will pass the preparation step, but the real-time events may overlap. There are two constraints to check that this will not happen:

the start of the readout must happen after the end of the excitation
the end of the readout must happend before the repetition time

Using the msl::RealTimeEvents::startTime and msl::RealTimeEvents::endTime functions, these constraints can be checked and an error status returned if needed:

// ...
#include <MrImaging/libSBB/libSBBmsg.h>
// ...
 
NLSStatus
Kernel
::prepare(MrProt & protocol, SeqLim & limits, MrProtocolData::SeqExpo & exports)
{
    // ...
    
    // Check timing: readout must start after excitation has ended and must end
    // before TR
    if(readout.startTime() < this->_excitation.endTime())
    {
        return MRI_SBB_SBB_NEGATIV_TEFILL;
    }
    if(readout.endTime() > protocol.tr()[0])
    {
        return MRI_SBB_SBB_NEGATIV_TRFILL;
    }
    
    return MRI_SEQ_SEQU_NORMAL;
}

Linking the Kernel and the Sequence

The final step is to fill the reconstruction-related information of the acquired data in the measurement data header or MDH. The MDH should be independent of the kernel as the kernel does not have the full information regarding the sequence graph and its loop structure: the code which fills the MDH will then be in the sequence. However, the MDH is contained in the sREADOUT object, itself in the msl::Readout object, which is prepared in the kernel. We can solve this by storing the sREADOUT in the registry, and adding a new parameter to the kernel constructor.

class Kernel: public msl::graph::AbstractNode
{
public:
    DECLARE_POINTERS(Kernel);
    static Pointer New(
        std::string const & sampling, std::string const & trajectories,
        std::string const & slices, std::string const & phaseCycling,
        std::string const & adc);
    // ...
private:
    // ...
    msl::DictionaryItem<sREADOUT *> _adc;
    // ...
    Kernel(
        std::string const & sampling, std::string const & trajectories,
        std::string const & slices, std::string const & phaseCycling,
        std::string const & adc);
};

Kernel::Pointer
Kernel
::New(
    std::string const & sampling, std::string const & trajectories,
    std::string const & slices, std::string const & phaseCycling,
    std::string const & adc)
{
    return Pointer(new Kernel(sampling, trajectories, slices, phaseCycling, adc));
}
// ...
Kernel
::Kernel(
    std::string const & sampling, std::string const & trajectories,
    std::string const & slices, std::string const & phaseCycling,
    std::string const & adc)
: AbstractNode(), _sampling(this, sampling), _trajectories(this, trajectories),
    _slices(this, slices), _phaseCycling(this, phaseCycling), _adc(this, adc)
{
    // Nothing else
}

We then store the correct value in the prepare function of the kernel: as opposed to the preceding uses of the registry, the informations flows out of the kernel.

NLSStatus
Kernel
::prepare(MrProt & protocol, SeqLim & limits, MrProtocolData::SeqExpo & exports)
{
    // ...
    // Update the ADC in the dictionary
    this->_adc() = &readout.adc();
    
    // Check timing: readout must start after excitation has ended and must end
    // before TR
    // ...
}

We need to adapt the sequence class to reflect this:

NLSStatus
FLASH
::initialize(SeqLim & limits)
{
    // ...
    auto trajectories = msl::graph::Loop::New(
        "trajectories", {
            Kernel::New(
                "sampling", "trajectories", "slices", "phaseCycling", "adc"),
            msl::graph::Action::New([](msl::Dictionary::Pointer d) {
                ++d->get<msl::PhaseCycling>("phaseCycling"); })
        }
    );
    // ...
    this->_root->registry()->insert({
        // ...
        {"adc", (sREADOUT*)(nullptr)}
    });
    // ...
}

Finally, we can update the MDH from the sequence, by adding an msl::graph::Action node in the sequence graph, before calling the kernel.

class FLASH: public msl::Sequence
{
    // ...
private:
    static void updateMDH(MrProt & protocol, msl::Dictionary::Pointer registry);
};

NLSStatus
FLASH
::initialize(SeqLim & limits)
{
    auto trajectories = msl::graph::Loop::New(
        "trajectories", {
            msl::graph::Action::New(FLASH::updateMDH),
            Kernel::New(
                "sampling", "trajectories", "slices", "phaseCycling", "adc"),
            msl::graph::Action::New([](msl::Dictionary::Pointer d) {
                ++d->get<msl::PhaseCycling>("phaseCycling"); })
        }
    );
    // ...
}

The updateMDH function will

Set the slice and average counters (the line and partition counters are set in msl::samplings::Cartesian)
Tell the reconstruction to perform the Fourier transform once all the data has been received, using the setPhaseFT and setPartitionFT functions
Notify the system when the current line is the first or the last

Its full code is:

void
FLASH
::updateMDH(MrProt & protocol, msl::Dictionary::Pointer registry)
{
    auto & adc = *registry->get<sREADOUT*>("adc");
    auto & mdh = adc.getMDH();
    
    auto const & trajectories = registry->get<msl::Counter>("trajectories");
    auto const & slices = registry->get<msl::SliceConstIterator>("slices");
    auto const & averages = registry->get<msl::Counter>("averages");
    
    // Set non sampling-related counters (sampling-related counters are handled
    // by the sampling)
    mdh.setCslc(uint16_t(slices.index()));
    mdh.setCacq(uint16_t(averages.index()));
    
    // NOTE: FT flags require both sampling-specific and sequence-specific
    // information (here last trajectory and last average).
    mdh.setPhaseFT(trajectories.last() && averages.last());
    mdh.setPartitionFT(trajectories.last() && averages.last());
    
    // Set the first/last flags
    mdh.setLastScanInMeas(trajectories.last() && slices.last() && averages.last());
    // NOTE: this assumes that there is only one concatenation
    mdh.setLastScanInConcat(mdh.isLastScanInMeas());
    
    mdh.setLastScanInSlice(trajectories.last());
    mdh.setFirstScanInSlice(trajectories.first());
    
    mdh.addToEvalInfoMask(MDH_ONLINE);
}

Full Code