codeare 1 - getting started

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

Measurement data (data)
K-Space mask (mask)

Provided in matlab file share/compressedsensing/brain512.mat

Parametric input

FFT class (fft)
Overall NLCG iterations (csiter)
Weight for total variation penalty (tvw)
Weight for data consistency penalty (xfmw)
Maximum CG internal iterations (cgiter)
CG convergence criterium
L₁ weight (l1)
Maximum # of line search iteration (lsiter)
P-norm (pnorm)
Line search brackets (lsa and lsb)
Wavelet family (wl_family)
Wavelet family member (wl_member)

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> 
    <data dname="data" fname="brain512.mat" dtype="float"
        ftype="MATLAB"/> 
    <mask dname="mask" fname="brain512.mat" dtype="float"
        ftype="MATLAB"/>
  </data>

</config>

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 ();
#else 
    /* 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");
#endif	

    if (runcs (con))
        return 0;
    else 
        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,

Initialise:
Intialise global variables, allocate memory, etc. Things that are initiallly setup possible time and memory consuming. For example, allocate the FT operator if enough information is avalable at this point.
Prepare:
Prepare the earlier initialised functionality and globals. For example, reset the Tikhonov regularisation weight to a new value before next processing.
Process:
Actually process the algorithm on the data.
Finalise:
Clean up before leaving.

We start by deriving a class CompressedSensing from ReconStrategy

Declaration

namespace RRStrategy {

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


      private:
    
        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.

codeare1.2

Getting started

Binary input

Parametric input

Writing the package XML file

Binary file

Coding the client

Coding the algorithm

Declaration

Implementation: Init

Implementation: Prepare

Implementation: Process

visit also