Title : Protyping at the scanner with MATLAB part 2
Schedule : June 25, 2020 | 16:00-17:00
Speakers : Stanislas Rapacchi & Aurélien Trotier
[TOC]
The Gadgetron responds to two major issues in MRI:
- prototyping: how to develop a new reconstruction and associate it with an existing or developing sequence.
- deployment: how to deploy a sequence and reconstruction on several sites for a clinical study
The Gadgetron also offers software flexibility (choice of language used) and hardware flexibility (choice of reconstruction hardware: from simple PC to Cloud). Specific aera of research impose constraints on reconstruction time or latency. Typically deployment on clinical sites or interventional imaging are two scenarios where the use of C ++ will be preferable. For most of the other thematics, languages such as Python or Matlab are more accessible and particularly adapted to research computing problems (ex: matrix calculation, linear algebra). On one side, Python is quite popular in image processing (itk, vtk) and in machine learing (keras, tensor flow). On the other side, Matlab remains the envirronment of choice for rapid prototyping in signal and image processing. This session will target such needs.
To do the tutorial, you need to install three components:
- MATLAB R2018b and above (Superior to R2017a is also possible but need a modification of gadgetron)
- gadgetron
- gadgetron-matlab
Gadgetron
Detailed installation instructions have been summarized here. But basically, on Ubuntu you need to run the following line:
sudo add-apt-repository ppa:gradient-software/experimental
sudo apt-get update
sudo apt-get install gadgetron-all
Matlab installed in Windows and Gadgetron in WSL
You need to install gadgetron from the source code (just reproduce a Docker file on your computer[Should we provide a complete script?]) and edit the file gadgetron/connection/stream/external/Matlab.cpp
Replace line 15 :
boost::process::search_path("matlab"),'
by
boost::process::search_path("matlab.exe"),
Recompile gadgetron.
And in WSL, before running gadgetron:
export GADGETRON_EXTERNAL_PORT
export GADGETRON_EXTERNAL_MODULE
export WSLENV = GADGETRON_EXTERNAL_PORT:GADGETRON_EXTERNAL_MODULE
then run gadgetron!
Note: a script could do that for you
Matlab version between R2017a and R2018a
You need to install gadgetron from the source code and edit the file gadgetron/connection/stream/external/Matlab.cpp
Replace line 16 :
boost::process::args={"-batch", "gadgetron.external.main"},'
by
boost::process::args={"--nosplash", "--nodesktop", "-r", "\"gadgetron.external.main; exit\""},
Need to be tested i'm not sure how to pass the args
Recompile gadgetron.
Gadgetron-Matlab
Detailed installation and how to use matlab is available here. Basically you can search gadgetron in the Add-On manager.
Optional
The following programm will be used at the end of the tutorial. You can skip this part at the beginning. ismrmrdviewer BART
To verify the Matlab-Gadgetron connection is working, one can type gadgetron --info. Matlab should be supported.
More throroughly, after installing gadgetron it is always a good idea to run the integration test. Here we will run the matlab test.
To do so move to the folder gadgetron/test/integration/
Download all the datasets with :
python get_data.py
and then run the matlab test with :
python run_tests.py cases/external_matlab_tiny_example.cfg
If it is working you will see an output with this indication at the end :
Test status: Passed
1 tests passed. 0 tests failed. 0 tests skipped.
We will use a 3D MP2RAGE sequence with a variable density poisson undersampling mask acquired on a 3T Prisma from Siemens.
Data is available at this link : Need to be change
One dataset (without noise calibration) is available:
- brain, 0.8mm isotropic, acceleration factor = 20.
The data has been converted with siemens_to_ismrmrd, we will not discuss data conversion here. This will be the object of the following readings.
- Review formalism of Matlab scripts and code
- Classes in +folder
- Data types definition
- Retrieve protocol headers (former xml) in Matlab
- Connection next actions
- to understand how data is stored under matlab dependinfg of your xml configuration file :
- Before AcquisitionAccumulateTriggerGadget
- After AcquisitionAccumulateTriggerGadget
- After BucketToBufferGadget
- to create a new matlab gadget using :
- the matlab debugger :
- launching it without a matlab session open :
- to send back images
- to call BART for advanced reconstruction method
- to become familiar with the Cartesian reconstruction pipeline
- to create new python gadget from scratch
- to create a new xml configuration file
- data manipulation (readout, kspace, image)
- to call BART from a python gagdet
- to call SigPy from a python gagdet
import numpy as np
import gadgetron
import ismrmrd
def EmptyPythonGadget(connection):
for acquisition in connection:
# DO SOMETHING
connection.send(acquisition)
The function is responsible for receiving all messages from the previous gadget and for sending a new message to the next gadget using connection.send() It may or may not interact with the information contained in the message.
Create a new directory named GT_Lecture3
and open two terminals at the location.
mkdir GT_Lecture3
cd GT_Lecture3
Create the file my_first_python_gadget.py then copy the previous class.
We can add the following message before connection.send() that we are going through it.
print("so far, so good")
We will now create a new xml file named external_python_tutorial.xml
. Add the following content into external_python_tutorial.xml
<?xml version="1.0" encoding="utf-8"?>
<gadgetronStreamConfiguration xsi:schemaLocation="http://gadgetron.sf.net/gadgetron gadgetron.xsd"
xmlns="http://gadgetron.sf.net/gadgetron"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
<!--
Gadgetron generic recon chain for 2D and 3D cartesian sampling
-->
<!-- reader -->
<reader><slot>1008</slot><dll>gadgetron_mricore</dll><classname>GadgetIsmrmrdAcquisitionMessageReader</classname></reader>
<reader><slot>1026</slot><dll>gadgetron_mricore</dll><classname>GadgetIsmrmrdWaveformMessageReader</classname></reader>
<!-- writer -->
<writer><slot>1022</slot><dll>gadgetron_mricore</dll><classname>MRIImageWriter</classname></writer>
</gadgetronStreamConfiguration>
To call our gadget, we have to add it to the reconstruction chain which is currently empty. For this add the following lines after the MRIImageWriter
<stream>
<external>
<execute name="my_first_python_gadget" target="EmptyPythonGadget" type="python"/>
<configuration/>
</external>
</stream>
Nothing to do.
To run the reconstruction chain, you'll need to run Gadgetron, and the Gadgetron ISMRMRD client.
Start Gadgetron:
$ gadgetron
Run the ISMRMRD client:
$ cd path/to/gadgetron/test/integration
$ gadgetron_ismrmrd_client -f in.h5 -C external_python_tutorial.xml
gadgetron_ismrmrd_client -f -C python_passthrough_tutorial.xml
Les données sont écrites par défaut dans un fichier out.h5.
Utiliser le script python pour les lire et afficher les images
Ajouter les lignes suivantes dans la fonction process.
self.counter=self.counter+1;
print(self.counter)
Puis rejouer l'étape compilation et installation et lancer la reco.
Ajouter les lignes suivantes dans la fonction process et import numpy as np en entête
print(np.shape(data))
Puis rejouer l'étape compilation et installation et lancer la reco.
Ajouter les lignes suivantes
print(header)
Nous allons obtenir ce type de message
2404
(256, 12)
version: 1
flags: 2097152
measurement_uid: 26
scan_counter: 2404
acquisition_time_stamp: 22684240
physiology_time_stamp: 8545183, 0, 0
number_of_samples: 256
available_channels: 12
active_channels: 12
channel_mask: 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
discard_pre: 0
discard_post: 0
center_sample: 128
encoding_space_ref: 0
trajectory_dimensions: 0
sample_time_us: 1.600000023841858
position: 0.0, -4.842615127563477, -75.69005584716797
read_dir: -6.123031769111886e-17, 1.0, 0.0
phase_dir: 1.0, 6.123031769111886e-17, 0.0
slice_dir: 0.0, 0.0, 1.0
patient_table_position: 0.0, 0.0, -1259016.0
idx: kspace_encode_step_1: 23
kspace_encode_step_2: 0
average: 0
slice: 0
contrast: 0
phase: 0
repetition: 2
set: 0
segment: 1
user: 0, 0, 0, 0, 0, 32, 0, 0
user_int: 0, 0, 0, 0, 0, 0, 0, 0
user_float: 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
Récupérer le numero de repetition et de slice avec les lignes suivantes:
slice = header.idx.slice
repetition= header.idx.repetition
e1=header.idx.kspace_encode_step_1
segment=header.idx.segment
Et ajouter la condition suivante pour ne pas transmettre tous les readouts des autres coupes et autres répétitions.
if (repetition>0 and slice>0):
print(self.counter_send, " slice: ",slice , " rep: ", repetition, " e1: ", e1," segment: ", segment)
Ajouter un second compteur nommé self.counter_send qui s'incrémente dans cette boucle.
self.myBuffer = None
self.enc = None
Dans process config ajouter ceci pour récupérer la taille de la amtrice:
self.header = ismrmrd.xsd.CreateFromDocument(conf) self.enc = self.header.encoding[0]
if self.myBuffer is None:
channels = acq.active_channels
eNz = self.enc.encodedSpace.matrixSize.z
eNy = self.enc.encodedSpace.matrixSize.y
eNx = self.enc.encodedSpace.matrixSize.x
self.myBuffer = np.zeros(( int(eNx),eNy,channels),dtype=np.complex64)
self.myBuffer[:,e1,:] = data
from matplotlib import transform
from ismrmrdtools import transform
if (e1==96)
plt.figure(1)
plt.imshow(np.abs(self.myBuffer[:,:,0]))
Pas forcément nécessaire de tout redévelopper
Il existe de nombreux gadgets en python et/ou en C++ qui permettent d'aller plus vite.
TODO: installation sans GPU / avec GPU docker avec/sans Python / Matlab
The data structures in the gadgetron vary during reconstruction. It is important to differenciate, the class or common structures
- used to store a unit of readout that would feed into a buffer
- used to store a unit of data that would feed into a reconstruction
- used to store an array of reconstructed data
Each of them are defined in a C++ and have equivalent in Python. Additionnal structure are also present and you can create new one
Le gadget python recevra deux messages associés qui contiennent le AcquisitionHeader et les données sous forme de matrice hoNDArray< std::complex >. En python le hoNDArray est la matrice multidimensionnel ndarray issue de la librairie numpy
process(self, header, data):
print(type(header))
print(type(data))
En imagerie cartesienne, deux gadgets jouent un role fondamental : AcquisitionAccumulateTriggerGadget et BucketToBufferGadget.
Ces gadgets servent à bufferiser les readouts afin de construire le kspace. En IRM, les dimensionalités sont très nombreuses:
- kx (RO)
- ky (E1)
- kz (E2)
- channels (CHA)
- average
- repetition
- segment
- contrast
- phase
- set
- slice (SLC)
- ...
Par convention, nous aurons en entrée des matrices de dimentions [RO, CHA] et en sortie de BucketToBufferGadget des matrices de dimensions [RO E1 E2 CHA N S SLC]. Les dimensions N et S peuvent être choisie arbitrairement par l'utilisateur.
Il est fort interessant de se positionner après ces gadgets ou sont automatiquement triés les données kspaces, que ce soit les lignes de calibration en imagerie parallèle ou les lignes sous échantionnées.
Les données de calibration si elles sont présentes sont accessibles via la structure suivante:
buffer.ref.data
buffer.ref.header
Les données "standard" sont accessibles via la structure suivante:
buffer.data.data
buffer.data.header
Attention, la taille des headers est associée la taille des données, les headers sont généralement différents, par ex la position des coupes changent suivant la direction SLC. Nous avons donc maintenant une matrice hondarray de acquisitionheader de dimensions [E1 E2 N S SLC]. Les headers étant identique suivant la direction de readout et pour tous les éléments d'antennes.
Le gadget python correspondant recevra un message qui contient une structure nommé IsmrmrdReconBit nommée IsmrmrdReconBit et IsmrmrdDataBuffered. Le gadget python doit donc contenir les includes suivants.
from gadgetron import Gadget,IsmrmrdDataBuffered, IsmrmrdReconBit, SamplingLimit,SamplingDescription, IsmrmrdImageArray
process(self, buffer):
Pour aller plus loin, voici les classes C++ correspondantes.
struct IsmrmrdReconBit
{
public:
IsmrmrdDataBuffered data_;
boost::optional<IsmrmrdDataBuffered> ref_;
}
struct IsmrmrdDataBuffered
{
public:
//7D, fixed order [E0, E1, E2, CHA, N, S, LOC]
hoNDArray< std::complex<float> > data_;
//7D, fixed order [TRAJ, E0, E1, E2, N, S, LOC]
boost::optional<hoNDArray<float>> trajectory_;
// 6D, density weights [E0, E1, E2, N, S, LOC]
boost::optional<hoNDArray<float> > density_;
//5D, fixed order [E1, E2, N, S, LOC]
hoNDArray< ISMRMRD::AcquisitionHeader > headers_;
SamplingDescription sampling_;
}
Nous allons donc maintenant créer un nouveau gadget nommé MyFirstDataBufferedGadget.
import sys
import ismrmrd
import ismrmrd.xsd
from gadgetron import Gadget,IsmrmrdDataBuffered, IsmrmrdReconBit, SamplingLimit,SamplingDescription, IsmrmrdImageArray
class MyFirstDataBufferedGadget(Gadget):
def __init__(self,next_gadget=None):
super(MyFirstDataBufferedGadget,self).__init__(next_gadget=next_gadget)
self.my_value = 0
self.my_list = []
self.my_matrix =[]
def process_config(self, conf):
# do allocation
pass
def process(self, message):
# get information from the message
# modify the message
# send the message to the next gadget
self.put_next(message)
return 0