Skip to content

Commit

Permalink
ENH: Outlier-robust SRR algorithm, Python 3 compatibility
Browse files Browse the repository at this point in the history
  • Loading branch information
Michael Ebner committed Sep 10, 2018
2 parents 3260823 + 14cf6ca commit a124799
Show file tree
Hide file tree
Showing 57 changed files with 4,392 additions and 846 deletions.
3 changes: 3 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -9,3 +9,6 @@ data/
*.pyc
*iterate.dat
*.egg-info
*.idea
*.wiki
/.project
52 changes: 39 additions & 13 deletions .gitlab-ci.yml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
# -----------------------------------Set Up------------------------------------
variables:
PY_VERSION: 2
CMICLAB: 0
PY_VERSION: 3
PRIVATE: 0
TMPDIR: ./tmp
DATADIR: /home/mebner/data/ci/FetalBrain
VENV: pysitk-test-py${PY_VERSION}
Expand Down Expand Up @@ -42,12 +42,12 @@ before_script:
- cp -v ${ITK_DIR}/Wrapping/Generators/Python/WrapITK.pth ${py_sitepkg}

- cd $cwd_dir
# If cmiclab is used:
# If PRIVATE is used:
# add CI_JOB_TOKEN for cloning dependent repositories in requirements.txt
# (https://docs.gitlab.com/ee/user/project/new_ci_build_permissions_model.html#dependent-repositories)
- >
(if [ ${CMICLAB} == 1 ];
then sed -i -- "s#github.com/gift-surg#gitlab-ci-token:${CI_JOB_TOKEN}@cmiclab.cs.ucl.ac.uk/GIFT-Surg#g" requirements.txt;
(if [ ${PRIVATE} == 1 ];
then sed -i -- "s#github.com/gift-surg#gitlab-ci-token:${CI_JOB_TOKEN}@PRIVATE.cs.ucl.ac.uk/GIFT-Surg#g" requirements.txt;
fi);
# install requirements
- pip install -r requirements.txt
Expand Down Expand Up @@ -85,9 +85,10 @@ reconstruct_volume_tk1l2:
# - master
script:
- >
python niftymic_reconstruct_volume.py
--dir-input ${DATADIR}/input_data
niftymic_reconstruct_volume
--filenames ${DATADIR}/input_data/axial.nii.gz ${DATADIR}/input_data/coronal.nii.gz ${DATADIR}/input_data/sagittal.nii.gz
--dir-output ${TMPDIR}
--suffix-mask _mask
--verbose 0
--isotropic-resolution 2
--reconstruction-type TK1L2
Expand All @@ -97,14 +98,16 @@ reconstruct_volume_tk1l2:
tags:
- gift-adelie


reconstruct_volume_huberl2:
# only:
# - master
script:
- >
python niftymic_reconstruct_volume.py
--dir-input ${DATADIR}/input_data
niftymic_reconstruct_volume
--filenames ${DATADIR}/input_data/axial.nii.gz ${DATADIR}/input_data/coronal.nii.gz ${DATADIR}/input_data/sagittal.nii.gz
--dir-output ${TMPDIR}
--suffix-mask _mask
--verbose 0
--isotropic-resolution 2
--reconstruction-type HuberL2
Expand All @@ -120,8 +123,10 @@ reconstruct_volume_from_slices:
# - master
script:
- >
python niftymic_reconstruct_volume_from_slices.py
--dir-input ${DATADIR}/motion_correction_oriented
niftymic_reconstruct_volume_from_slices
--filenames ${DATADIR}/input_data/axial.nii.gz ${DATADIR}/input_data/coronal.nii.gz ${DATADIR}/input_data/sagittal.nii.gz
--dir-input-mc ${DATADIR}/motion_correction_oriented
--suffix-mask _mask
--reconstruction-space ${DATADIR}/SRR_stacks3_TK1_lsmr_alpha0p03_itermax10_oriented.nii.gz
--dir-output ${TMPDIR}
--verbose 0
Expand All @@ -130,14 +135,35 @@ reconstruct_volume_from_slices:
tags:
- gift-adelie

run_reconstruction_pipeline:
# only:
# - master
script:
- >
python niftymic/application/run_reconstruction_pipeline.py
--filenames ${DATADIR}/input_data/axial.nii.gz
--dir-input-templates ${DATADIR}/templates
--dir-output ${TMPDIR}
--suffix-mask _mask
--gestational-age 28
--bias-field-correction 1
--two-step-cycles 0
--iter-max 1
--run-data-vs-simulated-data 1
--verbose 0
tags:
- gift-adelie

param_study_huberl2:
# only:
# - master
script:
- recon_type=HuberL2
- >
python niftymic_run_reconstruction_parameter_study.py
--dir-input ${DATADIR}/motion_correction_oriented
niftymic_run_reconstruction_parameter_study
--filenames ${DATADIR}/input_data/axial.nii.gz ${DATADIR}/input_data/coronal.nii.gz ${DATADIR}/input_data/sagittal.nii.gz
--dir-input-mc ${DATADIR}/motion_correction_oriented
--suffix-mask _mask
--reference ${DATADIR}/SRR_stacks3_TK1_lsmr_alpha0p03_itermax10_oriented.nii.gz
--dir-output ${TMPDIR}/param_study
--alpha-range 0.01 0.05 2
Expand Down
80 changes: 33 additions & 47 deletions README.md
Original file line number Diff line number Diff line change
@@ -1,28 +1,30 @@
# Motion Correction and Volumetric Image Reconstruction of 2D Ultra-fast MRI

NiftyMIC is a Python-based open-source toolkit for research developed within the [GIFT-Surg][giftsurg] project to reconstruct an isotropic, high-resolution volume from multiple, possibly motion-corrupted, stacks of low-resolution 2D slices. The framework relies on slice-to-volume registration algorithms for motion correction and reconstruction-based Super-Resolution (SR) techniques for the volumetric reconstruction.
Please note that currently **only Python 2** is supported.

The algorithm and software were developed by [Michael Ebner][mebner] at the [Translational Imaging Group][tig] in the [Centre for Medical Image Computing][cmic] at [University College London (UCL)][ucl].

If you have any questions or comments, please drop an email to `[email protected]`.

## Features
Several methods have been implemented to solve the **Super-Resolution Reconstruction** (SRR) problem
Several methods have been implemented to solve the **Robust Super-Resolution Reconstruction (SRR)** problem
<!-- ```math -->
<!-- \vec{x}^* := \text{argmin}_{\vec{x}} \Big[\sum_{k=1}^K \sum_{i=1}^{N_k} \varrho\big( (A_k\vec{x} - \vec{y}_k)_i^2 \big) + \alpha\,\text{Reg}(\vec{x}) \Big] -->
<!-- \vec{x}^* := \text{argmin}_{\vec{x}\ge 0} \Big[\sum_{k\in\mathcal{K}_\sigma} \sum_{i=1}^{N_k} \varrho\big( (A_k\vec{x} - \vec{y}_k)_i^2 \big) + \alpha\,\text{Reg}(\vec{x}) \Big] -->
<!-- ``` -->
<p align="center">
<img src="http://latex.codecogs.com/svg.latex?%5Cvec%7Bx%7D%5E%2A%3A%3D%5Ctext%7Bargmin%7D_%7B%5Cvec%7Bx%7D%7D%5CBig%5B%5Csum_%7Bk%3D1%7D%5EK%5Csum_%7Bi%3D1%7D%5E%7BN_k%7D%5Cvarrho%5Cbig%28%28A_k%5Cvec%7Bx%7D-%5Cvec%7By%7D_k%29_i%5E2%5Cbig%29%2B%5Calpha%5C%2C%5Ctext%7BReg%7D%28%5Cvec%7Bx%7D%29%5CBig%5D">
<img src="https://latex.codecogs.com/svg.latex?\vec{x}^*&space;:=&space;\text{argmin}_{\vec{x}\ge&space;0}&space;\Big[\sum_{k\in\mathcal{K}}&space;\sum_{i=1}^{N_k}&space;\varrho\big(&space;(A_k\vec{x}&space;-&space;\vec{y}_k)_i^2&space;\big)&space;&plus;&space;\alpha\,\text{Reg}(\vec{x})&space;\Big]">
</p>

<!--to obtain the (vectorized) high-resolution 3D MRI volume $`\vec{x}\in\mathbb{R}^N`$ from multiple, possibly motion corrupted, low-resolution stacks of (vectorized) 2D MR slices $`\vec{y}_k \in\mathbb{R}^{N_k}`$ with $`N_k\ll N`$ for $`k=1,...,\,K`$-->
<!--for a variety of regularizers $`\text{Reg}`$ and data loss functions $`\varrho`$.-->
<!--The linear operator $`A_k := D_k B_k W_k`$ represents the comd operator descri the (rigid) motion $`W_k`$, the blurring operator $`B_k`$ and the downsampling operator $`D_k`$.-->
to obtain the (vectorized) high-resolution 3D MRI volume ![img](http://latex.codecogs.com/svg.latex?%5Cvec%7Bx%7D%5Cin%5Cmathbb%7BR%7D%5EN) from multiple, possibly motion corrupted, low-resolution stacks of (vectorized) 2D MR slices ![img](http://latex.codecogs.com/svg.latex?%5Cvec%7By%7D_k%5Cin%5Cmathbb%7BR%7D%5E%7BN_k%7D) with ![img](http://latex.codecogs.com/svg.latex?N_k%5Cll%7BN%7D) for ![img](http://latex.codecogs.com/svg.latex?k%3D1%2C...%2C%5C%2CK)
to obtain the (vectorized) high-resolution 3D MRI volume ![img](http://latex.codecogs.com/svg.latex?%5Cvec%7Bx%7D%5Cin%5Cmathbb%7BR%7D%5EN) from multiple, possibly motion corrupted, low-resolution stacks of (vectorized) 2D MR slices ![img](http://latex.codecogs.com/svg.latex?%5Cvec%7By%7D_k%5Cin%5Cmathbb%7BR%7D%5E%7BN_k%7D) with ![img](http://latex.codecogs.com/svg.latex?N_k%5Cll%7BN%7D) for ![img](https://latex.codecogs.com/svg.latex?k\in\mathcal{K}=\\{1,\dots,K\\})
for a variety of regularizers ![img](http://latex.codecogs.com/svg.latex?%5Ctext%7BReg%7D) and data loss functions ![img](http://latex.codecogs.com/svg.latex?%5Cvarrho).
The linear operator ![img](http://latex.codecogs.com/svg.latex?A_k%3A%3DD_kB_kW_k) represents the combined operator describing the (rigid) motion ![img](http://latex.codecogs.com/svg.latex?W_k), the blurring operator ![img](http://latex.codecogs.com/svg.latex?B_k) and the downsampling operator ![img](http://latex.codecogs.com/svg.latex?D_k).

The toolkit relies on an iterative motion-correction/reconstruction approach whereby **complete outlier rejection** of misregistered slices is achieved by iteratively solving the SRR problem for a slice-index set
![img](https://latex.codecogs.com/svg.latex?\mathcal{K}_\sigma:=\\{1\le&space;k&space;\le&space;K:&space;\text{Sim}(A_k^{\iota}\vec{x}^{\iota},&space;\vec{y}^{\iota}_k)\ge\sigma\\}\subset&space;\mathcal{K}) containing only slices that are in high agreement with their simulated counterparts projected from a previous high-resolution iterate ![img](http://latex.codecogs.com/svg.latex?\vec{x}^{\iota}) according to a similarity measure ![img](http://latex.codecogs.com/svg.latex?\text{Sim}) and parameter ![img](http://latex.codecogs.com/svg.latex?\sigma>0). In the current implementation, the similarity measure ![img](http://latex.codecogs.com/svg.latex?\text{Sim}) corresponds to Normalized Cross Correlation (NCC).

---

The provided **data loss functions** ![img](http://latex.codecogs.com/svg.latex?%5Cvarrho) are motivated by [SciPy](https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.optimize.least_squares.html) and allow for robust outlier rejection. Implemented data loss functions are:
Expand All @@ -48,7 +50,6 @@ The **available regularizers** include
* First-order Tikhonov (TK1): ![img](http://latex.codecogs.com/svg.latex?%5Ctext%7BReg%7D%28%5Cvec%7Bx%7D%29%3D%5Cfrac%7B1%7D%7B2%7D%5CVert%5Cnabla%5Cvec%7Bx%7D%5CVert_%7B%5Cell%5E2%7D%5E2)
* Isotropic Total Variation (TV): ![img](http://latex.codecogs.com/svg.latex?%5Ctext%7BReg%7D%28%5Cvec%7Bx%7D%29%3D%5Ctext%7BTV%7D_%5Ctext%7Biso%7D%28%5Cvec%7Bx%7D%29%3D%5Cbig%5CVert%7C%5Cnabla%5Cvec%7Bx%7D%7C%5Cbig%5CVert_%7B%5Cell%5E1%7D)
* Huber Function: ![img](http://latex.codecogs.com/svg.latex?%5Ctext%7BReg%7D%28%5Cvec%7Bx%7D%29%3D%5Cfrac%7B1%7D%7B2%5Cgamma%7D%5Cbig%7C%7C%5Cnabla%5Cvec%7Bx%7D%7C%5Cbig%7C_%7B%5Cgamma%7D)

---

Additionally, the choice of finding **optimal reconstruction parameters** is facilitated by the [Numerical Solver Library (NSoL)][nsol].
Expand All @@ -59,19 +60,17 @@ NiftyMIC supports medical image registration and volumetric reconstruction for u
**NiftyMIC is not intended for clinical use**.

## How to cite
If you use this software in your work, please cite [Ebner et al., 2018][citation].
If you use this software in your work, please cite

* Ebner, M., Wang, G., Li, W., Aertsen, M., Patel, P. A., Melbourne, A., Doel, T., David, A. L., Deprest, J., Ourselin, S., & Vercauteren, T. (2018). An Automated Localization, Segmentation and Reconstruction Framework for Fetal Brain MRI. In Medical Image Computing and Computer-Assisted Intervention -- MICCAI 2018. Springer.
* Ebner, M., Chung, K. K., Prados, F., Cardoso, M. J., Chard, D. T., Vercauteren, T., & Ourselin, S. (2018). Volumetric reconstruction from printed films: Enabling 30 year longitudinal analysis in MR neuroimaging. NeuroImage, 165, 238–250.


## Installation

NiftyMIC is currently supported for **Python 2 only** and was tested on
NiftyMIC was developed in Ubuntu 16.04 and Mac OS X 10.12 and tested for Python 2.7.12 and 3.5.2.

* Mac OS X 10.10 and 10.12
* Ubuntu 14.04 and 16.04

NiftyMIC builds on a couple of additional libraries developed within the [GIFT-Surg][giftsurg] project including
It builds on a couple of additional libraries developed within the [GIFT-Surg][giftsurg] project including
* [NSoL][nsol]
* [SimpleReg][simplereg]
* [PySiTK][pysitk]
Expand All @@ -91,21 +90,16 @@ Provided the input MR image data in NIfTI format (`nii` or `nii.gz`), NiftyMIC c
### Volumetric MR Reconstruction from Motion Corrupted 2D Slices
Leveraging a two-step registration-reconstruction approach an isotropic, high-resolution 3D volume can be generated from multiple stacks of low-resolution slices.

Examples for basic usage are:
```
niftymic_reconstruct_volume \
--dir-input dir-with-multiple-stacks \
--dir-output output-dir \
--suffix-mask _mask
```
An example for a basic usage reads
```
niftymic_reconstruct_volume \
--filenames path-to-stack1.nii.gz ... path-to-stackN.nii.gz \
--dir-output output-dir \
--suffix-mask _mask
```
whereby complete outlier removal during SRR is activated by default (`--outlier-rejection 1`).

The obtained motion-correction transformations can be stored for further processing, e.g. by using `niftymic_reconstruct_volume_from_slices.py` to solve the SRR problem for a variety of different regularization and data loss function types.
The obtained motion-correction transformations can be used for further processing, e.g. by using `niftymic_reconstruct_volume_from_slices.py` to solve the SRR problem for a variety of different regularization and data loss function types.

### SRR Methods for Motion Corrected (or Static) Data

Expand All @@ -115,47 +109,36 @@ After performed motion correction (or having static data in the first place),
1. parameter studies can be performed to find optimal reconstruction parameters.

#### 1. SRR from Motion Corrected (or Static) Slice Acquisitions
Solve the SRR problem for motion corrected data:
```
niftymic_reconstruct_volume_from_slices \
--dir-input dir-to-motion-correction \
--dir-output output-dir \
--reconstruction-type HuberL2 \
--alpha 0.003
```

Solve the SRR problem for motion corrected data (or static data if `--dir-input-mc` is omitted):
```
niftymic_reconstruct_volume_from_slices \
--dir-input dir-to-motion-correction \
--filenames path-to-stack1.nii.gz ... path-to-stackN.nii.gz \
--dir-input-mc dir-to-motion_correction \
--dir-output output-dir \
--reconstruction-type TK1L2 \
--alpha 0.03
```

Solve the SRR problem for static data:
```
niftymic_reconstruct_volume_from_slices \
--filenames path-to-stack1.nii.gz ... path-to-stackN.nii.gz \
--dir-input-mc dir-to-motion_correction \
--dir-output output-dir \
--reconstruction-type HuberL2 \
--alpha 0.003
--suffix-mask _mask
```
```
niftymic_reconstruct_volume_from_slices \
--filenames path-to-stack1.nii.gz ... path-to-stackN.nii.gz \
--dir-output output-dir \
--reconstruction-type TK1L2 \
--alpha 0.03 \
--suffix-mask _mask
--alpha 0.003
```

Slices that were rejected during the `niftymic_reconstruct_volume` run are recognized as outliers based on the content of `dir-input-mc` and will not be incorporated during the volumetric reconstruction.


#### 2. Parameter Studies to Determine Optimal SRR Parameters
The optimal choice for reconstruction parameters like the regularization parameter or data loss function can be found by running parameter studies. This includes L-curve studies and direct comparison against a reference volume for various cost functions.

Example are:
```
niftymic_run_reconstruction_parameter_study \
--dir-input dir-to-motion-correction \
--filenames path-to-stack1.nii.gz ... path-to-stackN.nii.gz \
--dir-input-mc dir-to-motion_correction \
--dir-output dir-to-param-study-output \
--reconstruction-type HuberL2 \
--reference path-to-reference-volume.nii.gz \
Expand All @@ -165,7 +148,8 @@ niftymic_run_reconstruction_parameter_study \
```
```
niftymic_run_reconstruction_parameter_study \
--dir-input dir-to-motion-correction \
--filenames path-to-stack1.nii.gz ... path-to-stackN.nii.gz \
--dir-input-mc dir-to-motion_correction \
--dir-output dir-to-param-study-output \
--reconstruction-type TVL2 \
--reference path-to-reference-volume.nii.gz \
Expand All @@ -175,7 +159,8 @@ niftymic_run_reconstruction_parameter_study \
```
```
niftymic_run_reconstruction_parameter_study \
--dir-input dir-to-motion-correction \
--filenames path-to-stack1.nii.gz ... path-to-stackN.nii.gz \
--dir-input-mc dir-to-motion_correction \
--dir-output dir-to-param-study-output \
--reconstruction-type TK1L2 \
--reference path-to-reference-volume.nii.gz \
Expand All @@ -193,7 +178,7 @@ niftymic_show_parameter_study \
```

## Licensing and Copyright
Copyright (c) 2017, [University College London][ucl].
Copyright (c) 2018, [University College London][ucl].
This framework is made available as free open-source software under the [BSD-3-Clause License][bsd]. Other licenses may apply for dependencies.


Expand All @@ -202,11 +187,12 @@ This work is partially funded by the UCL [Engineering and Physical Sciences Rese

## References
Associated publications are
* [[Ebner2018]][citation] Ebner, M., Chung, K. K., Prados, F., Cardoso, M. J., Chard, D. T., Vercauteren, T., & Ourselin, S. (2018). Volumetric reconstruction from printed films: Enabling 30 year longitudinal analysis in MR neuroimaging. NeuroImage, 165, 238–250.
* [EbnerWang2018] Ebner, M., Wang, G., Li, W., Aertsen, M., Patel, P. A., Melbourne, A., Doel, T., David, A. L., Deprest, J., Ourselin, S., & Vercauteren, T. (2018). An Automated Localization, Segmentation and Reconstruction Framework for Fetal Brain MRI. In Medical Image Computing and Computer-Assisted Intervention -- MICCAI 2018. Springer
* [[Ebner2018]](https://www.sciencedirect.com/science/article/pii/S1053811917308042) Ebner, M., Chung, K. K., Prados, F., Cardoso, M. J., Chard, D. T., Vercauteren, T., & Ourselin, S. (2018). Volumetric reconstruction from printed films: Enabling 30 year longitudinal analysis in MR neuroimaging. NeuroImage, 165, 238–250.
* [[Ranzini2017]](https://mski2017.files.wordpress.com/2017/09/miccai-mski2017.pdf) Ranzini, M. B., Ebner, M., Cardoso, M. J., Fotiadou, A., Vercauteren, T., Henckel, J., Hart, A., Ourselin, S., and Modat, M. (2017). Joint Multimodal Segmentation of Clinical CT and MR from Hip Arthroplasty Patients. MICCAI Workshop on Computational Methods and Clinical Applications in Musculoskeletal Imaging (MSKI) 2017.
* [[Ebner2017]](https://link.springer.com/chapter/10.1007%2F978-3-319-52280-7_1) Ebner, M., Chouhan, M., Patel, P. A., Atkinson, D., Amin, Z., Read, S., Punwani, S., Taylor, S., Vercauteren, T., and Ourselin, S. (2017). Point-Spread-Function-Aware Slice-to-Volume Registration: Application to Upper Abdominal MRI Super-Resolution. In Zuluaga, M. A., Bhatia, K., Kainz, B., Moghari, M. H., and Pace, D. F., editors, Reconstruction, Segmentation, and Analysis of Medical Images. RAMBO 2016, volume 10129 of Lecture Notes in Computer Science, pages 3–13. Springer International Publishing.

[citation]: https://www.sciencedirect.com/science/article/pii/S1053811917308042

[mebner]: http://cmictig.cs.ucl.ac.uk/people/phd-students/michael-ebner
[tig]: http://cmictig.cs.ucl.ac.uk
[bsd]: https://opensource.org/licenses/BSD-3-Clause
Expand Down
12 changes: 10 additions & 2 deletions install_cli.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,13 +2,21 @@
import sys
import os
import re
import six


DIR_ROOT = os.path.dirname(os.path.abspath(__file__))
DIR_CPP = os.path.join(DIR_ROOT, "niftymic", "cli")
DIR_CPP_BUILD = os.path.join(DIR_ROOT, "build", "cpp")


##
# Compile and install the cpp-code associated with NiftyMIC.
#
# Prior to running `python install_cli.py` set the environment variable
# accordingly. E.g. `export NIFTYMIC_ITK_DIR=path-to-ITK-build`. Moreover, make
# sure Boost is installed, e.g. `sudo apt install libboost-all-dev`
# \date 2018-01-30 10:00:40+0000
#
def main(prefix_environ="NIFTYMIC_"):

# Get current working directory
Expand All @@ -20,7 +28,7 @@ def main(prefix_environ="NIFTYMIC_"):
environment_vars = {p.match(f).group(1): p.match(f).group(0)
for f in os.environ.keys() if p.match(f)}
cmake_args = []
for k, v in environment_vars.iteritems():
for k, v in six.iteritems(environment_vars):
cmake_args.append("-D %s=%s" % (k, os.environ[v]))
cmake_args.append(DIR_CPP)

Expand Down
Loading

0 comments on commit a124799

Please sign in to comment.