Writing a Gadget

A Gadget is a Node in the Gadgetron chain, which processes data coming in through an GenericInputChannel and sends the processed data to the next Node in the chain using an

OutputChannel.

The simplest Gadgets to write are PureGadget and ChannelGadget.

PureGadget

A PureGadget is a Gadget which processes Messages one at a time, and holds no state. Examples could be a Gadget which removes oversampling on Acquisitions, or one which takes an Image and performs autoscaling.

A PureGadget inheritss from PureGadget<OUTPUT,INPUT>, where OUTPUT and INPUT are the output type and input type of the Gadget.

AutoScaleGadget.h

/**
    \brief  Autoscales real-type images based on a given max value and the 99th percentile of the data
    \test   Tested by: epi_2d.cfg and others
*/

#pragma once

#include "PureGadget.h"
#include "Types.h"
#include "hoNDArray_math.h"
#include <algorithm>

namespace Gadgetron{
    class AutoScaleGadget : public Core::PureGadget<Core::AnyImage, Core::AnyImage> {
    public:
        using Core::PureGadget<Core::AnyImage,Core::AnyImage>::PureGadget;
        Core::AnyImage process_function(Core::AnyImage image) const override;
    protected:
        NODE_PROPERTY(max_value, float, "Maximum value (after scaling)", 2048);
        NODE_PROPERTY(histogram_bins, unsigned int, "Number of Histogram Bins", 100);
        float current_scale_;
        std::vector<size_t> histogram_;
    };
}

The NODE_PROPERTY macro defines a variable on the AutoScaleGadget which can be set from the XML file defining the chain.

AutoScaleGadget.cpp

#include "AutoScaleGadget.h"

using namespace Gadgetron;
using namespace Gadgetron::Core;

namespace {
	template<class T>
	Image<T> autoscale(Image<T> &image, float max_value, unsigned int histogram_bins) {
		auto &header = std::get<ISMRMRD::ImageHeader>(image);
		auto &data = std::get<hoNDArray<T>>(image);
		auto &meta = std::get<optional<ISMRMRD::MetaContainer>>(image);	
		if (header.image_type == ISMRMRD::ISMRMRD_IMTYPE_MAGNITUDE) { //Only scale magnitude images for now
			auto max = *std::max_element(data.begin(), data.end());
			auto current_scale_ = 1.0;
			auto histogram_ = std::vector<size_t>(histogram_bins);    
			std::fill(histogram_.begin(), histogram_.end(), 0);
			for (auto& d : data) {
				size_t bin = static_cast<size_t>(std::floor((d/max)*histogram_bins));
				if (bin >= histogram_bins) {
					bin = histogram_bins-1;
				}
				histogram_[bin]++;
			}
			//Find 99th percentile
			long long cumsum = 0;
			size_t counter = 0;
			while (cumsum < (0.99*data.get_number_of_elements())) {
				cumsum += (long long)(histogram_[counter++]);
			}
			max = (counter+1)*(max/histogram_bins);
			current_scale_ = max_value/max;
			for (auto& d : data){
				d *= current_scale_;
			}
		}
		return Image<T>(header,data,meta);
	}

    template<class T>
    Image<std::complex<T>> autoscale(Image<std::complex<T>> &image, float max_value, unsigned int histogram_bins) {
        GERROR("Autoscaling image is not well defined for complex images.");
		return image;
    }
}

namespace Gadgetron {
	AnyImage AutoScaleGadget::process_function(AnyImage image) const {
		return visit([&](auto &image) -> AnyImage { return autoscale(image, max_value, histogram_bins); }, image);
	}
	GADGETRON_GADGET_EXPORT(AutoScaleGadget);
}

Note the GADGETRON_GADGET_EXPORT declaration, which produces the code causing the AutoScaleGadget to be loadable by Gadgetron.

ChannelGadget

PureGadget can’t hold any state between different messages and must send one message per input. This makes it easier to reason about and implement, but is also limiting in cases where we want to accumulate multiple messages for processing. In this case, ChannelGadget should be used. If we want to create a Gadget which takes several Acquisitions and reconstruct them, we inherit from ChannelGadget<Acquisition>.

using namespace Gadgetron;
using namespace Gadgetron::Core;
class SimpleRecon : public ChannelGadget<Acquisition>{
    public:
        void process(InputChannel<Acquisition>& in, OutputChannel& out) override {
            ...
        }
}

We can take the messages from the channel either by calling InputChannel::pop() directly, or by using it in

a for loop.

Channels are ranges, meaning they can be used directly with for loops and with algorithm from standard library, such as std::transform() and std::accumulate().

void process(InputChannel<Acquisition>& in, OutputChannel& out) override {
    for (auto acquisition : in ) {

        auto& header = std::get<ISMRMRD::AcquisitionHeader>(acquisition);
        auto& data = std::get<hoNDArray<std::complex<float>>>(acquisition);
        auto& trajectory = std::get<optional<hoNDArray<float>>>(acquisition);
        //Gather acquisitions here
    }
}

Or if you’re using C++17, this would be

void process(InputChannel<Acquisition>& in, OutputChannel& out) override {
    for (auto [header, data, trajectory] : in ) {
        //Gather acquisitions here
    }
}

We want to gather acquisitions until we have enough for a (possibly undersampled) image. The AcquisitionHeader has the

ISMRMRD::_ACQ_LAST_IN_ENCODE_STEP1 flag which we can use as a trigger. By importing channel_algorithms.h, we can write

void process(InputChannel<Acquisition>& in, OutputChannel& out) override {

    auto split_condition = [](auto& message){
      return std::get<ISMRMRD::AcquisitionHeader>(message).isFlagSet(ISMRMRD::ISMRMRD_ACQ_LAST_IN_ENCODE_STEP1);
    };

    for (auto acquisitions : buffer(in,split_condition)) {
        for (auto [header, data, trajectory] : acquisitions ) {
        //Gather acquisitions here
        }
}

#include <gadgetron/Gadget.h>
#include <gadgetron/hoNDFFT.h>
#include <gadgetron/mri_core_utility.h>
#include <gadgetron/ChannelAlgorithms.h>
#include <gadgetron/log.h>
#include <gadgetron/mri_core_coil_map_estimation.h>
using namespace Gadgetron;
using namespace Gadgetron::Core;
using namespace Gadgetron::Core::Algorithm;

class SimpleRecon : public ChannelGadget<Acquisition> {

    public:
        SimpleRecon(const Context& context, const GadgetProperties& params) : ChannelGadget<Acquisition>(params), header{context.header} {

        }

        void process(InputChannel<Acquisition>& in, OutputChannel& out){

            auto recon_size = header.encoding[0].encodedSpace.matrixSize;

            ISMRMRD::AcquisitionHeader saved_header;

            auto split_condition = [](auto& message){
            return std::get<ISMRMRD::AcquisitionHeader>(message).isFlagSet(ISMRMRD::ISMRMRD_ACQ_LAST_IN_ENCODE_STEP1);
            };

            for (auto acquisitions : buffer(in,split_condition)) {

               auto data = hoNDArray<std::complex<float>>(recon_size.x,recon_size.y,recon_size.z,header.acquisitionSystemInformation->receiverChannels.get());
               for ( auto [acq_header, acq_data, trajectories] : acquisitions){
                    saved_header = acq_header;
                    data(slice,acq_header.idx.kspace_encode_step_1,0,slice) = acq_data;
                }

                hoNDFFT<float>::instance()->fft2c(data);

                auto coil_map = coil_map_Inati(data);
                data = coil_combine(data,coil_map,3);

                auto image_header = image_header_from_acquisition(saved_header,header,data);

                out.push(image_header,data);
            }
        }
    private:
        const ISMRMRD::IsmrmrdHeader header;
};

GADGETRON_GADGET_EXPORT(SimpleRecon)