Getting started

In this tutorial we walk through the implementation of Compressed Sensing for Cartesian data acquisition and its extension to arbitraty accelerations along arbitrary trajectories.

Binary input

Provided in matlab file share/compressedsensing/brain512.mat

Parametric input

Writing the package XML file

We concentrate the above setup in a package file in XML brain512.xml:

<?xml version="1.0" ?>

  dim:    Dimension
  N[x-z]: Side lengths
  maxint: Maximum optimisation runs
  tvw:    Total variation penalty weight
  xfmw:   Data consistency weight
  ftoper: Fourier transform class
          0: Cartesian FFT
          1: Cartesian SENSE
          2: Cartesian GRAPPA
          3: Non-Cartesian FFT
          4: Non-Cartesian SENSE
  cgiter: Maximum # CG iterations
  cgconv: CG vonvergence criterium
  l1:     L1 weight
  pnorm:  P-Norm 
  lslim:  Line search lim
  lsa:    Line search alpha
  lsb:    Line search beta
  lsto:   Line search TO

< config; csiter="2" tvw="0.002" xfmw="0.005" fft="0" cgiter="8"
    cgconv="5.0e-3" l1="1.0e-15" lsiter="10" pnorm="1" lslim="10"
    lsa="0.01" lsb="0.6" wl_family="0" wl_member="4"/>

    <data dname="data" fname="brain512.mat" dtype="float"
    <mask dname="mask" fname="brain512.mat" dtype="float"


Binary file

The package file above explicitly includes the file locations and the variable names for the matrices data and mask. In this case data and mask reside in the MATLAB file brain512.mat in the same directory as the package file.

This binary input could be extracted from SyngoMR, HDF5 and NIFTI files

Coding the client

The package file must be read, binary data digested and announced to the matrix database.

#include "RemoteConnector.hpp"
#include "LocalConnector.hpp"
#include "IO.hpp"

template <class T> bool runcs (Connector<T>& con) {
    Matrix<cxfl>  indata;
    Matrix<cxfl>  im_dc;
    Matrix<float> mask;
    Matrix<float> pdf;
    Matrix<cxfl>  pc;
    /* Read package configuration */
	con.ReadConfig ("brain512.xml");

    /* Need at least measured data  */
    if (!Read (pdf, con.GetElement("/config/data/data"),  "") 
        return false; 

    /* Optionals */
	Read (mask, con.GetElement("/config/data/mask"), "");
	Read (pdf,  con.GetElement("/config/data/pdf"),  "");
	Read (pc,   con.GetElement("/config/data/pc"),   "");

    if (!(MXRead (pdf,  df, "pdf")))  pdf  = Matrix<float>(1);
    if (!(MXRead (mask, df, "mask"))) mask = Matrix<float>(1);
    if (!(MXRead (pc,   df, "ph")))   pc   = Matrix<cxfl> (1);

    /* Load shared module, and initialise */
    if (con.Init ("CompressedSensing") != OK) {
        printf ("Intialising failed ... bailing out!\n"); 
        return false;

    /* Announce binary data matrices to database of con   */
    con.SetMatrix  ("data", indata); /* Measurement data */
    con.SetMatrix  ("pdf",  pdf);    /* Sensitivities    */
    con.SetMatrix  ("mask", mask);   /* Weights          */
    con.SetMatrix  ("pc",   pc);     /* Phase correction */

    /* Prepare and process */
    con.Prepare  ("CompressedSensing");     
    con.Process  ("CompressedSensing");

    /* Retrieve result from database */
    con.GetMatrix ("im_dc", im_dc);

    /* Clean up module */
    con.Finalise ("CompressedSensing");

    /* Write output */
    if (!(MXDump (im_dc, "ph"))) {
        printf ("Writing output failed ... bailing out!\n")
        return false;

    /* :) */
    return true;


int main (int argc, char** argv) {

#ifdef LOCAL
    /* Locally load the CS module */
    Connector<LocalConnector>  con ();
    /* Remote call to codeare server loading the CS module.
       EXAMPLE: codeserv is resolved by CORBA name server debug level 5*/
    Connector<RemoteConnector> con ("codserv", "5");

    if (runcs (con))
        return 0;
        return 1;


Coding the algorithm

codeare was programmed to keep the data handling and the algorithms on matrices (or ND arrays for that matter) away from the implementational parts as much as possible.

Thus, we shall start by creating a new ReconStrategy by deriving from the base class. We then only need to implement the algorithm in 4 contextually different functions,

We start by deriving a class CompressedSensing from ReconStrategy


namespace RRStrategy {

    class CompressedSensing : public ReconStrategy {
        CompressedSensing  () {};
        virtual ~CompressedSensing () {};
        virtual RRSModule::error_code Process ();
        virtual RRSModule::error_code Init ();
        virtual RRSModule::error_code Process ();
        virtual RRSModule::error_code Finalise ();

        int            m_dim;     /* Image recon dim               */
        int            m_N[3];    /* Image side lengths            */
        int            m_csiter;  /* # global iterations           */
        CGParam        m_cgparam; /* Structure handed over to NLCG */
        int            m_wf;      /* Wavelet family                */
        int            m_wm;      /* Family member                 */



Download CompressedSensing.hpp

Implementation: Init

The implementation of the Init method consists of getting initialising the appropriate local variables with the meta data provided by the package file.

RRSModule::error_code CompressedSensing::Init () {

    for (size_t i = 0; i < 3; i++)
        m_N[i] = 1;

    int wli   = 0;
    int m_fft = 0;

    Attribute ("tvw",       &m_cgparam.tvw);
    Attribute ("xfmw",      &m_cgparam.xfmw);
    Attribute ("l1",        &m_cgparam.l1);
    Attribute ("pnorm",     &m_cgparam.pnorm);
    Attribute ("fft",       &m_cgparam.fft);
    Attribute ("csiter",    &m_csiter);
    Attribute ("wl_family", &m_wf);
    Attribute ("wl_member", &m_wm);
    if (m_wf < -1 || m_wf > 5)
        m_wf = -1;

    Attribute ("cgconv", &m_cgparam.cgconv);
    Attribute ("cgiter", &m_cgparam.cgiter);
    Attribute ("lsiter", &m_cgparam.lsiter);
    Attribute ("lsa",    &m_cgparam.lsa);
    Attribute ("lsb",    &m_cgparam.lsb);

    m_initialised = true;

    return RRSModule::OK;


Implementation: Prepare

DWT, TVOP and FT operators need only be declared once but probably only at a later stage than the initialisation?

RRSModule::error_code CompressedSensing::Prepare () {

    /* Declare TVOP, DWT and FT operators */
    m_cgparam.tvt = new TVOP ();
    m_cgparam.dwt = new DWT (data.Height(), wlfamily(m_wf), m_wm);
    m_cgparam.ft  = (FT<cxfl>*) 
        new DFT<cxfl> (ndims (data)+1, data.Height(), mask, pc);


Implementation: Process

To keep things simple, only the non-accelerated single channel Cartesian code is discussed. The more general code is found here.

Data matrices are retrieved from the database - data, pdf, mask and pc. For convenience references to the prepared transform operators (in this case on std::complex<float> are obtained (dft and dwt).

NOTE: The transform operators can be used like mathematical operators (i.e. mhat = ft * m is the forward and m = ft ->* mhat is the inverse transform).

RRSModule::error_code CompressedSensing::Process () {

    float ma; /* max(abs(data)) */

    /* Get scan data, density compensation, k-space mask and phase
       correction matrices */
    Matrix<cxfl>&  data = GetCXFL ("data");
    Matrix<float>& pdf  = GetRLFL ("pdf" );
    Matrix<float>& mask = GetRLFL ("mask");
    Matrix<cxfl>&  pc   = GetCXFL ("pc");  

    /* Outgoing images are declared to the database by the name im_dc */
    Matrix<cxfl>&  im_dc = 
        AddMatrix ("im_dc", (Ptr<Matrix<cxfl> >) NEW (Matrix<cxfl>  (data.Dim())));

    /* For convenience */
    FT<cxfl>& dft = *m_cgparam.ft;
    DWT& dwt      = *m_cgparam.dwt;
    /* Compensate for k-space coverage density */
    im_dc    = data;
    im_dc   /= pdf;
    /* Reconstruct image with FT operator*/
    im_dc    = dft ->* im_dc;
    /* Normalise data magnitude */
    ma       = max(abs(im_dc));
    im_dc   /= ma;
    data    /= ma;

    /* Wavelet transform with DWT operator */    
    im_dc    = dwt * im_dc;
    /* NLCG runs */
    for (size_t i = 0; i < m_csiter; i++)
        NLCG (im_dc, data, m_cgparam);

    /* Assign outgoing images */
    im_dc    = dwt ->* im_dc * ma;

    /* Return control to client */
    return OK;


The remaining implementation resides in the NLCG optimisation which is declared as a global static function in CompressedSensing.hpp. As an example how coding the algorithmic part is made easy and intuitive to write and read. Note, that this is C++ and that we do not have a shell as provided by matlab including lexer and parser.

Matrix<cxfl> GradTV    (Matrix<cxfl>& x, CGParam& cgp) {

    /* References to the DWT and TVOP operators for convenience */
    DWT&  dwt = *cgp.dwt;
    TVOP& tvt = *cgp.tvt;

    float p   = ((float)cgp.pnorm)/2.0-1.0;

    Matrix&lt;cxfl> dx, g;

    /* dx = tvt(idwt(x)) */
    dx = tvt * (dwt->*x);

    /* g  = p * dx .* (dx .* dx.' + cgp.l1) .^ (p/2-1); */
    g  = dx * conj(dx);
    g += cxfl(cgp.l1);
    g ^= p;
    g *= dx;
    g *= cxfl(cgp.pnorm);

    /* g = dwt(itvt(g)) */
    g  = dwt * (tvt->*g);

    return (cgp.tvw * g);


Yippie! We're done. Obviously we need not read / write data from / to disk, when we are running the module inside the online image reconstruction software of an MRI scanner.

visit also

jemris mri simulator
open-source full-featured parallel tx/rx mri sequence and hardware simulator by tony stoecker, kaveh vahedipour and daniel pflugfelder.

bloch simulator for education in mri and nmr
free educational mri and nmr sofware. multi- os/platform without installation of software. remarkable and beautiful effort by lars hanson from drcmr, hvidovre, dk.

a comparable project to codeare by michael hansen and thomas sorensen.

mri unbound
a collaborative forum for mri data acquisition and image reconstruction maintained by the unbreakable jim pipe.