## 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)

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
• L1 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"/>
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> pdf;
Matrix<cxfl>  pc;

/* Need at least measured data  */
return false;

/* Optionals */

if (!(MXRead (pdf,  df, "pdf")))  pdf  = 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  ("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
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                 */

};

}```

### 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<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.

gagdetron
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.