janniklasrose / rwcdti

Random Walk for cardiac Diffusion Tensor Imaging

Random Walk for cardiac Diffusion Tensor Imaging

Compiling

The code is using MATLAB, but performance-critical functions have been written in C in the form of MEX functions. They need to be compiled to be available at runtime. To do this, execute compile.m inside MATLAB.

Running a simulation

To run a simulation, it is easiest to specify all parameters using a configuration file (using the YAML format). Example configuration files are provided in the config directory.

For small scale simulations with a low number of walkers, executing on a local machine (laptop/desktop) will be sufficient. If large scale jobs are desired, parallelisation is almost certainly required. The pbs.sh file should give an indication how to run these on the HPC cluster. For storing simulation data of significant size, it is best to use the more portable HDF5 format which is generally faster than using compressed MAT files.

Publication

Citation

If this software is useful to you, please consider citing the original paper:

Rose, JN, Nielles-Vallespin, S, Ferreira, P, Firmin, DN, Scott, AD, Doorly, DJ. Novel insights into in-vivo diffusion tensor cardiovascular magnetic resonance using computational modelling and a histology‐based virtual microstructure. Magn Reson Med. 2019; 81: 2759–2773. DOI:10.1002/mrm.27561

The tagged releases are archived on Zenodo.

Dataset

Simulation data (both input geometries and results) from the MRM paper is published here.

rwcdti's Issues

Correct the alpha timing parameters in the example sequences

The values in sequences.yml are currently set to 1 as placeholders. They needs to be updated to the correct values used for simulations in the paper.

We should probably provide a small (ASCII) diagram or documentation to explain how all sequence parameters are defined/interpreted, especially because alpha is not a standard parameter.

Restrict step according to the dimension of the simulation

If the user wants to run 1D or 2D simulations the step is added to all three components and is not restricted to either one/two components. Add substrate.dimension in the substrate.yml config file and pass in this information into function one_dt when the step is created. (Currently the dimension of the step is hardcoded to dim=3)

Handle input units

The user should be able to run simulations in their units of choice. While the units in the example config files are scaled SI units (it is encouraged to keep track of this somewhere, in this case it is done here), an external geometry file may be in SI units or some other unit depending on how it was generated. Pulse sequence parameters are often expressed in scanner units.

In general, we assume the user knows what they're doing. As long as the units are consistent, there is no need to convert to SI. For the geometry, a scale option in combination with the file option inside config.yml should do the trick to specify a different scale.

• Add support for scale option when reading a geometry from file
• Clarify in documentation that no assumption regarding units is made by the software

Set up testing framework

Continuous integration is now possible with MATLAB and it supports GitHub Actions.

To start:

• Provide config .yaml files that fully define the test cases
• Execute them as scripts
• GitHub-hosted runners have multiple cores, so both serial and parallel functionality can and should be tested

Advanced

More advanced testing is possible with the Testing Frameworks available in MATLAB.

Parallel performance

If I recall correctly, a significant time is spent remapping the output of the parallel computation. This should be profiled to see if it's worth refactoring.

There is also a warning inside the parfor loop regarding func. This can be mitigated by introducing an anonymous function fn = @(i) func(i, varargin{:}); before the loop, but I'm not sure if this is the right way to go...

Implement calculation of compartment residence time

Improve documentation

• Document UI
  • Add proper documentation to each function (including minimal doc for private functions)
  • Consistent naming scheme: dromedaryCase for functions (even private ones), CamelCase for class properties
• Tidy & review some functions
  • Substrate.Transform, Substrate.Transform.global2local, Substrate.Transform.local2global
  • ParticleWalker.seedParticlesInBox
  • ParticleWalker.one_dt

Implement boundary conditions

Allow the user to choose boundary conditions at the end of [Lx,Ly,Lz]. Currently , the user can choose two types of geometries: full or block. In both cases, the walker uses function one_dt() and checks if the walker is stepping out of the bounding box [Lx,Ly,Lz] . Periodic boundary would mean that the walker is free to step into the following block whereas reflect would not allow walkers crossing from one box to the other.

Implement a config parameter inside substrate.yml in order to be able to switch between these two options. The parameter should be called boundary: reflect or boundary:periodic.

Property validation and access permissions

Consider adding validators and updating access permissions