Loading Raw Data - sickkids-mri/mrrecon GitHub Wiki

Basic Usage

The mrrecon.data.reader.DataLoader class can be used to load raw data.

from mrrecon.data.reader import DataLoader

filename = '/home/user/mrirawdata.dat'

data = DataLoader(filename).run()

This returns a dictionary data containing k-space data, calibration data, and other parameters relevant for image reconstruction. To print a short summary of the items in data, the mrrecon.data.reader.print_dict_summary function can be used.

from mrrecon.data.reader import print_dict_summary

print_dict_summary(data)

Out:

NAME               TYPE                                 SHAPE OR VALUE
kspace             <class 'numpy.ndarray'> complex64    (26, 2, 2100, 512)
nx                 <class 'int'>                        256
ny                 <class 'int'>                        256
nz                 <class 'int'>                        1
fovx               <class 'float'>                      256.0
fovy               <class 'float'>                      256.0
fovz               <class 'float'>                      4.0
dx                 <class 'float'>                      1.0
dy                 <class 'float'>                      1.0
dz                 <class 'float'>                      4.0
tr                 <class 'float'>                      12.16
te                 <class 'float'>                      3.67
ti                 <class 'float'>                      300.0
flipangle          <class 'float'>                      20.0
venc               <class 'float'>                      150.0
veldir             <class 'int'>                        4
weight             <class 'float'>                      80.0
times              <class 'numpy.ndarray'> float32      (4200,)
user_float         <class 'numpy.ndarray'> float32      (24, 4200)
noise              <class 'numpy.ndarray'> complex64    (26, 256, 512)

The names in the left column are the dictionary keys. The second column shows the variable type, i.e. the type of scalar or type of array, and the third column shows the value of the scalar or shape of the array.

Custom Loading

It is often necessary to adjust DataLoader for new sequences with different data shapes and different parameters necessary for image reconstruction. The DataLoader class generally contains 5 steps:

  1. Reading the Siemens raw data file with twixtools. DataLoader uses twixtools for interfacing with the Siemens raw data file and parsing headers. Though twixtools is general enough for loading all types of acquisitions, the output of twixtools is too 'raw' and not as simple as it should be for image reconstruction. The next 4 steps will process the output from twixtools into a convenient dictionary called self.data. Step 1 should not require any modification.

  2. Storing scan data in NumPy arrays. A Siemens raw data file may contain multiple 'scans' or 'measurements', e.g. noise measurements or k-space measurements. This step will loop through each scan and store the acquired data in NumPy arrays. self.data['kspace'] is a list that holds k-space arrays, and self.data['calib'] is a list that holds calibration data arrays. In most cases there is only 1 array in self.data['kspace'] and 0 or 1 array in self.data['calib']. This step is also very general and should not require any modification.

  3. Picking out relevant scan parameters from the header (e.g. image resolution, TR, VENC). Only a small handful of parameters are required for image reconstruction, whereas the Siemens raw data file contains hundreds of parameters in a complicated structure. Picking out relevant scan parameters and storing them in self.data will make referencing them later a lot easier. All parameters can still be accessed from the original output of twixtools (self.scan_list).

  4. Reading the mini data header for per line data values (e.g. time stamps, custom user-defined data). Similarly, the mini data header (MDH) may contain data necessary for image reconstruction (e.g. cardiac gating). If steps 3 and 4 require modification, it may make the most sense to modify DataLoader directly.

  5. Reformatting the data. This final step transposes and reshapes arrays, since their format from previous steps may still be suboptimal. This step is likely to vary the most between different sequences. The example below will demonstrate how this method can be modified.

Example

The original DataLoader class was written to reformat 2D radial phase contrast data. To modify the data reformatting step for radial 4D flow data, inherit DataLoader and override the _reformat method.

import mrrecon as mr

class Flow4DLoader(mr.data.reader.DataLoader):  # Inherits methods from original DataLoader
    def _reformat(self):  # Overrides old method
        # New code goes here
        return

As with 2D phase contrast, there is only one k-space scan and one noise scan in 4D flow, so the single array in the lists created in step 2 should be taken out.

        self.data['kspace'] = self.data['kspace'][0]
        self.data['calib'] = self.data['calib'][0]

Since the calibration scan is actually a measurement of noise, rename the key from 'calib' to 'noise'.

        self.data['noise'] = self.data.pop('calib')

The native shape for the arrays produced in step 2 is (ncoils, nlines, nro), where nro is the number of readouts or number of samples along a line of k-space. However, the different velocity encoded lines are interleaved in the nlines dimension. The following code will store the different velocity encodes separately for easier referencing:

        nv = 4  # Number of velocity encodes

        # Reshape
        (ncoils, nlines, nro) = self.data['kspace'].shape

        tmp = np.empty((ncoils, nv, int(nlines/nv), nro), dtype=np.complex64)

        for v in range(nv):
            tmp[:, v, :, :] = self.data['kspace'][:, v::nv, :]

        self.data['kspace'] = tmp
        tmp = None  # Clear memory

Lastly, the time stamps are natively stored with poor precision. If necessary, they can be recalculated at higher precision, using the first time stamp as reference.

        # Recalculate times at higher precision
        time0 = self.data['times'][0]
        times = np.linspace(time0,
                            time0 + (nlines - 1) * (self.data['tr'] / nv),
                            num=nlines, dtype=np.float32)
        self.data['times'] = times

The new loader class with the modified _reformat method looks like this:

import mrrecon as mr

class Flow4DLoader(mr.data.reader.DataLoader):  # Inherits methods from original DataLoader
    def _reformat(self):  # Overrides old method
        self.data['kspace'] = self.data['kspace'][0]
        self.data['calib'] = self.data['calib'][0]

        self.data['noise'] = self.data.pop('calib')

        nv = 4  # Number of velocity encodes

        # Reshape
        (ncoils, nlines, nro) = self.data['kspace'].shape

        tmp = np.empty((ncoils, nv, int(nlines/nv), nro), dtype=np.complex64)

        for v in range(nv):
            tmp[:, v, :, :] = self.data['kspace'][:, v::nv, :]

        self.data['kspace'] = tmp
        tmp = None  # Clear memory

        # Recalculate times at higher precision
        time0 = self.data['times'][0]
        times = np.linspace(time0,
                            time0 + (nlines - 1) * (self.data['tr'] / nv),
                            num=nlines, dtype=np.float32)
        self.data['times'] = times
        return

Note: For extremely large k-space arrays, reshaping may be expensive and it may make sense to modify _fill_array so that the desired shape is created from step 2.

The new class can be used like this:

filename = '/home/user/4dflowdata.dat'

data = Flow4DLoader(filename).run()