diff --git a/.gitignore b/.gitignore index 0e6abe09..86729b63 100644 --- a/.gitignore +++ b/.gitignore @@ -9,3 +9,6 @@ data/ *.pyc *iterate.dat *.egg-info +*.idea +*.wiki +/.project diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 41df8997..350f3fe2 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -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} @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/README.md b/README.md index 871af361..27df9572 100644 --- a/README.md +++ b/README.md @@ -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 `michael.ebner.14@ucl.ac.uk`. ## 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 - +

- +

-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: @@ -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]. @@ -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] @@ -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 @@ -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 \ @@ -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 \ @@ -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 \ @@ -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. @@ -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 diff --git a/install_cli.py b/install_cli.py index 8498843c..f1abb823 100644 --- a/install_cli.py +++ b/install_cli.py @@ -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 @@ -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) diff --git a/niftymic/application/convert_nifti_to_dicom.py b/niftymic/application/convert_nifti_to_dicom.py new file mode 100644 index 00000000..02b6b109 --- /dev/null +++ b/niftymic/application/convert_nifti_to_dicom.py @@ -0,0 +1,149 @@ +## +# \file convert_nifti_to_dicom.py +# \brief Script to convert a 3D NIfTI image to DICOM. +# +# \author Michael Ebner (michael.ebner.14@ucl.ac.uk) +# \date Mar 2018 +# + +# NOTES (quite some candidates were tried to get a working solution): +# +# - nifti2dicom (https://packages.ubuntu.com/xenial/science/nifti2dicom; tested version 0.4.11): +# Although nifti2dicom allows the import of a DICOM header from a template +# (-d) not all tags would be set correctly. E.g. if DOB is not given at +# template, it would just be set to 01.01.1990 which would prevent the +# resulting dcm file to be grouped correctly with the original data. +# Moreover, annoying tags like 'InstitutionName' are set to their predefined +# value which cannot be deleted (only overwritten). +# Apart from that, only a relatively small selection of tags can be edited. +# However, it does a good job in creating a series of 2D DICOM slices from a +# NIfTI file (including correct image orientation!). +# +# - medcon (https://packages.ubuntu.com/xenial/medcon; tested version 0.14.1): +# A single 3D dcm file can be created but image orientation is flawed +# when created from a nifti file directly. +# However, if a 3D stack is created from a set of 2D dicoms, the orientation +# stays correct. +# +# - pydicom (https://github.com/pydicom/pydicom; tested version 1.0.2): +# Can only read a single 3D dcm file. In particular, it is not possible +# to read a set of 2D slices unless a DICOMDIR is provided which is not +# always guaranteed to exist (I tried to create it from 2D slices using +# dcmmkdir from dcmtk and dcm4che -- neither seemed to work reliably) +# Once the dicom file is read, pydicom does a really good job of updating +# DICOM tags + there are plenty of tags available to be chosen from! +# Saving a single 3D DICOM file is very easy too then. + +import os +import pydicom + +import pysitk.python_helper as ph + +from niftymic.utilities.input_arparser import InputArgparser +from niftymic.definitions import DIR_TMP + + +COPY_DICOM_TAGS = { + # important for grouping + "PatientID", + "PatientName", + "PatientBirthDate", + "StudyInstanceUID", + + # additional information + "StudyID", + "AcquisitionDate", + "PatientSex", + "MagneticFieldStrength", + "Manufacturer", + "ManufacturerModelName", + "Modality", + "StudyDescription", +} + + +def main(): + + input_parser = InputArgparser( + description="Convert obtained SRR from nifti to dicom format", + ) + input_parser.add_filename(required=True) + input_parser.add_option( + option_string="--template", + type=str, + required=True, + help="Template DICOM to extract relevant DICOM tags.", + ) + input_parser.add_dir_output(required=True) + input_parser.add_label( + help="Label used for series description of DICOM output.", + default="SRR") + + args = input_parser.parse_args() + input_parser.print_arguments(args) + + # Prepare for final DICOM output + ph.create_directory(args.dir_output) + path_to_output = os.path.join(args.dir_output, "%s.dcm" % args.label) + + # Prepare for intermediate output + dir_output_2d_slices = os.path.join(DIR_TMP, "dicom_slices") + ph.create_directory(dir_output_2d_slices, delete_files=True) + + # Create set of 2D DICOM slices from 3D NIfTI image + # (correct image orientation!) + ph.print_title("Create set of 2D DICOM slices from 3D NIfTI image") + cmd_args = ["nifti2dicom"] + cmd_args.append("-i %s" % args.filename) + cmd_args.append("-o %s" % dir_output_2d_slices) + cmd_args.append("-y") + ph.execute_command(" ".join(cmd_args)) + + # Combine set of 2D DICOM slices to form 3D DICOM image + # (image orientation stays correct) + ph.print_title("Combine set of 2D DICOM slices to form 3D DICOM image") + cmd_args = ["medcon"] + cmd_args.append("-f %s/*.dcm" % dir_output_2d_slices) + cmd_args.append("-o %s" % path_to_output) + cmd_args.append("-c dicom") + cmd_args.append("-stack3d") + cmd_args.append("-n") + cmd_args.append("-qc") + cmd_args.append("-w") + ph.execute_command(" ".join(cmd_args)) + + # Update all relevant DICOM tags accordingly + ph.print_title("Update all relevant DICOM tags accordingly") + print("") + dataset = pydicom.dcmread(path_to_output) + dataset_template = pydicom.dcmread(args.template) + + # Copy tags from template (to guarantee grouping with original data) + update_dicom_tags = {} + for tag in COPY_DICOM_TAGS: + try: + update_dicom_tags[tag] = getattr(dataset_template, tag) + except: + update_dicom_tags[tag] = "" + + # Additional tags + update_dicom_tags["InstitutionName"] = "UCL, WEISS" + update_dicom_tags["SeriesDescription"] = args.label + update_dicom_tags["ImageComments"] = "*** NOT APPROVED ***" + update_dicom_tags["AccessionNumber"] = "1" + update_dicom_tags["SeriesNumber"] = "0" + + for tag in sorted(update_dicom_tags.keys()): + value = update_dicom_tags[tag] + setattr(dataset, tag, value) + ph.print_info("%s: %s" % (tag, value)) + + dataset.save_as(path_to_output) + print("") + ph.print_info("3D DICOM image written to %s" % path_to_output) + + return 0 + + +if __name__ == '__main__': + main() diff --git a/niftymic/application/correct_bias_field.py b/niftymic/application/correct_bias_field.py index b9e924ba..63f7da42 100755 --- a/niftymic/application/correct_bias_field.py +++ b/niftymic/application/correct_bias_field.py @@ -56,16 +56,14 @@ def main(): "the width of the Gaussian deconvolution.", default=0.15, ) - input_parser.add_log_script_execution(default=1) + input_parser.add_log_config(default=1) input_parser.add_verbose(default=0) args = input_parser.parse_args() input_parser.print_arguments(args) - # Write script execution call - if args.log_script_execution: - input_parser.write_performed_script_execution( - os.path.abspath(__file__)) + if args.log_config: + input_parser.log_config(os.path.abspath(__file__)) # Read data data_reader = dr.MultipleImagesReader( diff --git a/niftymic/application/correct_intensities.py b/niftymic/application/correct_intensities.py index 1b724904..00efa64a 100755 --- a/niftymic/application/correct_intensities.py +++ b/niftymic/application/correct_intensities.py @@ -37,7 +37,7 @@ def main(): input_parser.add_suffix_mask(default="_mask") input_parser.add_search_angle(default=180) input_parser.add_prefix_output(default="IC_") - input_parser.add_log_script_execution(default=1) + input_parser.add_log_config(default=1) input_parser.add_option( option_string="--registration", type=int, @@ -49,10 +49,8 @@ def main(): args = input_parser.parse_args() input_parser.print_arguments(args) - # Write script execution call - if args.log_script_execution: - input_parser.write_performed_script_execution( - os.path.abspath(__file__)) + if args.log_config: + input_parser.log_config(os.path.abspath(__file__)) if args.reference in args.filenames: args.filenames.remove(args.reference) diff --git a/niftymic/application/multiply_stack_with_mask.py b/niftymic/application/multiply_stack_with_mask.py new file mode 100644 index 00000000..88b80c25 --- /dev/null +++ b/niftymic/application/multiply_stack_with_mask.py @@ -0,0 +1,56 @@ +## +# \file mulitply_stack_with_mask.py +# \brief Script to stack/reconstruction with multiply template mask. +# +# \author Michael Ebner (michael.ebner.14@ucl.ac.uk) +# \date October 2017 +# + +import os + +import pysitk.python_helper as ph +import pysitk.simple_itk_helper as sitkh + +import niftymic.base.stack as st +from niftymic.utilities.input_arparser import InputArgparser +from niftymic.definitions import DIR_TEMPLATES + + +def main(): + + input_parser = InputArgparser( + description="Multiply stack/reconstruction with template mask.", + ) + + input_parser.add_filename(required=True) + input_parser.add_gestational_age(required=True) + input_parser.add_dir_input_templates(default=DIR_TEMPLATES) + input_parser.add_verbose(default=1) + input_parser.add_dir_output() + input_parser.add_prefix_output(default="Masked_") + args = input_parser.parse_args() + input_parser.print_arguments(args) + + template = os.path.join( + args.dir_input_templates, + "STA%d.nii.gz" % args.gestational_age) + template_mask = os.path.join( + args.dir_input_templates, + "STA%d_mask_dil.nii.gz" % args.gestational_age) + + stack = st.Stack.from_filename(args.filename, template_mask, + extract_slices=False) + stack_masked = stack.get_stack_multiplied_with_mask() + stack_masked.set_filename(args.prefix_output + stack.get_filename()) + + if args.dir_output is None: + dir_output = os.path.dirname(args.filename) + else: + dir_output = args.dir_output + stack_masked.write(dir_output) + + if args.verbose: + sitkh.show_stacks([stack, stack_masked], segmentation=stack) + +if __name__ == '__main__': + main() diff --git a/niftymic/application/reconstruct_volume.py b/niftymic/application/reconstruct_volume.py index 12f577b6..df494c95 100755 --- a/niftymic/application/reconstruct_volume.py +++ b/niftymic/application/reconstruct_volume.py @@ -24,6 +24,7 @@ import niftymic.registration.niftyreg as niftyreg import niftymic.registration.simple_itk_registration as regsitk import niftymic.utilities.data_preprocessing as dp +import niftymic.utilities.intensity_correction as ic import niftymic.utilities.segmentation_propagation as segprop import niftymic.utilities.volumetric_reconstruction_pipeline as pipeline import niftymic.utilities.joint_image_mask_builder as imb @@ -37,7 +38,6 @@ def main(): # Set print options for numpy np.set_printoptions(precision=3) - # Read input input_parser = InputArgparser( description="Volumetric MRI reconstruction framework to reconstruct " "an isotropic, high-resolution 3D volume from multiple stacks of 2D " @@ -48,12 +48,12 @@ def main(): "this region will then be reconstructed by the SRR algorithm which " "can substantially reduce the computational time.", ) - input_parser.add_dir_input() - input_parser.add_filenames() + input_parser.add_filenames(required=True) + input_parser.add_filenames_masks() input_parser.add_dir_output(required=True) input_parser.add_suffix_mask(default="_mask") input_parser.add_target_stack_index(default=0) - input_parser.add_search_angle(default=90) + input_parser.add_search_angle(default=45) input_parser.add_multiresolution(default=0) input_parser.add_shrink_factors(default=[2, 1]) input_parser.add_smoothing_sigmas(default=[1, 0]) @@ -67,53 +67,40 @@ def main(): input_parser.add_dilation_radius(default=3) input_parser.add_extra_frame_target(default=10) input_parser.add_bias_field_correction(default=0) - input_parser.add_intensity_correction(default=0) + input_parser.add_intensity_correction(default=1) input_parser.add_isotropic_resolution(default=None) - input_parser.add_log_script_execution(default=1) + input_parser.add_log_config(default=1) input_parser.add_subfolder_motion_correction() input_parser.add_provide_comparison(default=0) input_parser.add_subfolder_comparison() input_parser.add_write_motion_correction(default=1) input_parser.add_verbose(default=0) input_parser.add_two_step_cycles(default=3) - input_parser.add_use_masks_srr(default=1) + input_parser.add_use_masks_srr(default=0) input_parser.add_boundary_stacks(default=[10, 10, 0]) + input_parser.add_metric(default="Correlation") + input_parser.add_metric_radius(default=10) input_parser.add_reference() input_parser.add_reference_mask() - + input_parser.add_outlier_rejection(default=1) + input_parser.add_threshold_first(default=0.6) + input_parser.add_threshold(default=0.7) + args = input_parser.parse_args() input_parser.print_arguments(args) - # Write script execution call - if args.log_script_execution: - input_parser.write_performed_script_execution( - os.path.abspath(__file__)) + if args.log_config: + input_parser.log_config(os.path.abspath(__file__)) # Use FLIRT for volume-to-volume reg. step. Otherwise, RegAladin is used. - use_flirt_for_v2v_registration = True + use_flirt_for_v2v_registration = 1 # --------------------------------Read Data-------------------------------- ph.print_title("Read Data") - - # Neither '--dir-input' nor '--filenames' was specified - if args.filenames is not None and args.dir_input is not None: - raise IOError( - "Provide input by either '--dir-input' or '--filenames' " - "but not both together") - - # '--dir-input' specified - elif args.dir_input is not None: - data_reader = dr.ImageDirectoryReader( - args.dir_input, suffix_mask=args.suffix_mask) - - # '--filenames' specified - elif args.filenames is not None: - data_reader = dr.MultipleImagesReader( - args.filenames, suffix_mask=args.suffix_mask) - - else: - raise IOError( - "Provide input by either '--dir-input' or '--filenames'") + data_reader = dr.MultipleImagesReader( + file_paths=args.filenames, + file_paths_masks=args.filenames_masks, + suffix_mask=args.suffix_mask) if len(args.boundary_stacks) is not 3: raise IOError( @@ -142,7 +129,6 @@ def main(): segmentation_propagator=segmentation_propagator, use_cropping_to_mask=True, use_N4BiasFieldCorrector=args.bias_field_correction, - use_intensity_correction=args.intensity_correction, target_stack_index=args.target_stack_index, boundary_i=args.boundary_stacks[0], boundary_j=args.boundary_stacks[1], @@ -184,8 +170,8 @@ def main(): else: vol_registration = niftyreg.RegAladin( registration_type="Rigid", - use_fixed_mask=True, - use_moving_mask=True, + # use_fixed_mask=True, + # use_moving_mask=True, use_verbose=False, ) v2vreg = pipeline.VolumeToVolumeRegistration( @@ -201,6 +187,33 @@ def main(): else: time_registration = ph.get_zero_time() + # ---------------------------Intensity Correction-------------------------- + if args.intensity_correction: + ph.print_title("Intensity Correction") + intensity_corrector = ic.IntensityCorrection() + intensity_corrector.use_individual_slice_correction(False) + intensity_corrector.use_reference_mask(True) + intensity_corrector.use_stack_mask(True) + intensity_corrector.use_verbose(False) + + for i, stack in enumerate(stacks): + if i == args.target_stack_index: + ph.print_info("Stack %d: Reference image. Skipped." % (i + 1)) + continue + else: + ph.print_info("Stack %d: Intensity Correction ... " % (i + 1), + newline=False) + intensity_corrector.set_stack(stack) + intensity_corrector.set_reference( + stacks[args.target_stack_index].get_resampled_stack( + resampling_grid=stack.sitk, + interpolator="NearestNeighbor", + )) + intensity_corrector.run_linear_intensity_correction() + stacks[i] = intensity_corrector.get_intensity_corrected_stack() + print("done (c1 = %g) " % + intensity_corrector.get_intensity_correction_coefficients()) + # ---------------------------Create first volume--------------------------- time_tmp = ph.start_timing() @@ -254,20 +267,26 @@ def main(): reg_type="TK1", minimizer="lsmr", alpha=args.alpha_first, - iter_max=args.iter_max_first, + iter_max=np.min([args.iter_max_first, args.iter_max]), verbose=True, use_masks=args.use_masks_srr, ) if args.two_step_cycles > 0: + if args.metric == "ANTSNeighborhoodCorrelation": + metric_params = {"radius": args.metric_radius} + else: + metric_params = None + registration = regsitk.SimpleItkRegistration( moving=HR_volume, use_fixed_mask=True, use_moving_mask=True, use_verbose=args.verbose, interpolator="Linear", - metric="Correlation", + metric=args.metric, + metric_params=metric_params, use_multiresolution_framework=args.multiresolution, shrink_factors=args.shrink_factors, smoothing_sigmas=args.smoothing_sigmas, @@ -289,6 +308,8 @@ def main(): cycles=args.two_step_cycles, alpha_range=[args.alpha_first, args.alpha], verbose=args.verbose, + outlier_rejection=args.outlier_rejection, + threshold_range=[args.threshold_first, args.threshold], ) two_step_s2v_reg_recon.run() HR_volume_iterations = \ @@ -297,6 +318,7 @@ def main(): two_step_s2v_reg_recon.get_computational_time_registration() time_reconstruction += \ two_step_s2v_reg_recon.get_computational_time_reconstruction() + stacks = two_step_s2v_reg_recon.get_stacks() else: HR_volume_iterations = [] @@ -306,10 +328,38 @@ def main(): stack.write( os.path.join(args.dir_output, args.subfolder_motion_correction), - write_mask=True, - write_slices=True, + write_stack=False, + write_mask=False, + write_slices=False, write_transforms=True, - suffix_mask=args.suffix_mask, + ) + + if args.outlier_rejection: + deleted_slices_dic = {} + for i, stack in enumerate(stacks): + deleted_slices = stack.get_deleted_slice_numbers() + deleted_slices_dic[stack.get_filename()] = deleted_slices + ph.write_dictionary_to_json( + deleted_slices_dic, + os.path.join( + args.dir_output, + args.subfolder_motion_correction, + "rejected_slices.json" + ) + ) + + if args.outlier_rejection: + deleted_slices_dic = {} + for i, stack in enumerate(stacks): + deleted_slices = stack.get_deleted_slice_numbers() + deleted_slices_dic[stack.get_filename()] = deleted_slices + ph.write_dictionary_to_json( + deleted_slices_dic, + os.path.join( + args.dir_output, + args.subfolder_motion_correction, + "rejected_slices.json" + ) ) # ------------------Final Super-Resolution Reconstruction------------------ @@ -357,8 +407,7 @@ def main(): segmentation=HR_volume, show_comparison_file=args.provide_comparison, dir_output=os.path.join( - args.dir_output, - args.subfolder_comparison), + args.dir_output, args.subfolder_comparison), ) # Summary diff --git a/niftymic/application/reconstruct_volume_from_slices.py b/niftymic/application/reconstruct_volume_from_slices.py index b22b2759..c0a32bb1 100755 --- a/niftymic/application/reconstruct_volume_from_slices.py +++ b/niftymic/application/reconstruct_volume_from_slices.py @@ -12,14 +12,15 @@ import numpy as np import os -# Import modules -import niftymic.base.data_reader as dr +import pysitk.python_helper as ph +import pysitk.simple_itk_helper as sitkh + import niftymic.base.stack as st +import niftymic.base.data_reader as dr import niftymic.reconstruction.admm_solver as admm +import niftymic.utilities.intensity_correction as ic import niftymic.reconstruction.primal_dual_solver as pd import niftymic.reconstruction.tikhonov_solver as tk -import pysitk.python_helper as ph -import pysitk.simple_itk_helper as sitkh from niftymic.utilities.input_arparser import InputArgparser @@ -36,15 +37,16 @@ def main(): "an isotropic, high-resolution 3D volume from multiple " "motion-corrected (or static) stacks of low-resolution slices.", ) - input_parser.add_dir_input() - input_parser.add_filenames() - input_parser.add_image_selection() + input_parser.add_filenames(required=True) + input_parser.add_filenames_masks() + input_parser.add_dir_input_mc() input_parser.add_dir_output(required=True) input_parser.add_prefix_output(default="SRR_") input_parser.add_suffix_mask(default="_mask") input_parser.add_target_stack_index(default=0) input_parser.add_extra_frame_target(default=10) input_parser.add_isotropic_resolution(default=None) + input_parser.add_intensity_correction(default=1) input_parser.add_reconstruction_space(default=None) input_parser.add_minimizer(default="lsmr") input_parser.add_iter_max(default=10) @@ -52,7 +54,7 @@ def main(): input_parser.add_data_loss(default="linear") input_parser.add_data_loss_scale(default=1) input_parser.add_alpha( - default=0.02 # TK1L2 + default=0.01 # TK1L2 # default=0.006 #TVL2, HuberL2 ) input_parser.add_rho(default=0.5) @@ -61,54 +63,72 @@ def main(): input_parser.add_iterations(default=15) input_parser.add_subfolder_comparison() input_parser.add_provide_comparison(default=0) - input_parser.add_log_script_execution(default=1) + input_parser.add_log_config(default=1) + input_parser.add_use_masks_srr(default=0) input_parser.add_verbose(default=0) + args = input_parser.parse_args() input_parser.print_arguments(args) - # Write script execution call - if args.log_script_execution: - input_parser.write_performed_script_execution( - os.path.abspath(__file__)) + if args.log_config: + input_parser.log_config(os.path.abspath(__file__)) + + if args.reconstruction_type not in ["TK1L2", "TVL2", "HuberL2"]: + raise IOError("Reconstruction type unknown") # --------------------------------Read Data-------------------------------- ph.print_title("Read Data") - # Neither '--dir-input' nor '--filenames' was specified - if args.filenames is not None and args.dir_input is not None: - raise IOError( - "Provide input by either '--dir-input' or '--filenames' " - "but not both together") - - # '--dir-input' specified - elif args.dir_input is not None: - data_reader = dr.ImageSlicesDirectoryReader( - path_to_directory=args.dir_input, - suffix_mask=args.suffix_mask, - image_selection=args.image_selection) - - # '--filenames' specified - elif args.filenames is not None: - data_reader = dr.MultipleImagesReader( - args.filenames, suffix_mask=args.suffix_mask) - - else: - raise IOError( - "Provide input by either '--dir-input' or '--filenames'") - - if args.reconstruction_type not in ["TK1L2", "TVL2", "HuberL2"]: - raise IOError("Reconstruction type unknown") + data_reader = dr.MultipleImagesReader( + file_paths=args.filenames, + file_paths_masks=args.filenames_masks, + suffix_mask=args.suffix_mask, + dir_motion_correction=args.dir_input_mc, + ) data_reader.read_data() stacks = data_reader.get_data() ph.print_info("%d input stacks read for further processing" % len(stacks)) + # ---------------------------Intensity Correction-------------------------- + if args.intensity_correction: + ph.print_title("Intensity Correction") + intensity_corrector = ic.IntensityCorrection() + intensity_corrector.use_individual_slice_correction(False) + intensity_corrector.use_stack_mask(True) + intensity_corrector.use_reference_mask(True) + intensity_corrector.use_verbose(False) + + for i, stack in enumerate(stacks): + if i == args.target_stack_index: + ph.print_info("Stack %d: Reference image. Skipped." % (i + 1)) + continue + else: + ph.print_info("Stack %d: Intensity Correction ... " % (i + 1), + newline=False) + intensity_corrector.set_stack(stack) + intensity_corrector.set_reference( + stacks[args.target_stack_index].get_resampled_stack( + resampling_grid=stack.sitk, + interpolator="NearestNeighbor", + )) + intensity_corrector.run_linear_intensity_correction() + stacks[i] = intensity_corrector.get_intensity_corrected_stack() + print("done (c1 = %g) " % + intensity_corrector.get_intensity_correction_coefficients()) + # Reconstruction space is given isotropically resampled target stack if args.reconstruction_space is None: recon0 = \ stacks[args.target_stack_index].get_isotropically_resampled_stack( resolution=args.isotropic_resolution, extra_frame=args.extra_frame_target) + recon0 = recon0.get_cropped_stack_based_on_mask( + boundary_i=args.extra_frame_target, + boundary_j=args.extra_frame_target, + boundary_k=args.extra_frame_target, + unit="mm", + ) # Reconstruction space was provided by user else: @@ -127,22 +147,36 @@ def main(): recon0 = recon0.get_stack_multiplied_with_mask() if args.reconstruction_type in ["TVL2", "HuberL2"]: - ph.print_title("Compute Initial value for %s" % args.reconstruction_type) - SRR0 = tk.TikhonovSolver( - stacks=stacks, - reconstruction=recon0, - alpha=args.alpha, - iter_max=args.iter_max, - reg_type="TK1", - minimizer=args.minimizer, - data_loss=args.data_loss, - data_loss_scale=args.data_loss_scale, - # verbose=args.verbose, - ) + ph.print_title("Compute Initial value for %s" % + args.reconstruction_type) + SRR0 = tk.TikhonovSolver( + stacks=stacks, + reconstruction=recon0, + alpha=args.alpha, + iter_max=np.min([5, args.iter_max]), + reg_type="TK1", + minimizer="lsmr", + data_loss="linear", + use_masks=args.use_masks_srr, + # verbose=args.verbose, + ) + else: + SRR0 = tk.TikhonovSolver( + stacks=stacks, + reconstruction=recon0, + alpha=args.alpha, + iter_max=args.iter_max, + reg_type="TK1", + minimizer=args.minimizer, + data_loss=args.data_loss, + data_loss_scale=args.data_loss_scale, + use_masks=args.use_masks_srr, + # verbose=args.verbose, + ) SRR0.run() recon = SRR0.get_reconstruction() - recon.set_filename(SRR0.get_setting_specific_filename(args.prefix_output)) + recon.set_filename(SRR0.get_setting_specific_filename(args.prefix_output)) recon.write(args.dir_output) # List to store SRRs @@ -163,6 +197,7 @@ def main(): rho=args.rho, data_loss=args.data_loss, iterations=args.iterations, + use_masks=args.use_masks_srr, verbose=args.verbose, ) SRR.run() @@ -174,7 +209,6 @@ def main(): recon.write(args.dir_output) else: - SRR = pd.PrimalDualSolver( stacks=stacks, reconstruction=st.Stack.from_stack(SRR0.get_reconstruction()), diff --git a/niftymic/application/register_image.py b/niftymic/application/register_image.py index 91bd3367..ec0d55f1 100755 --- a/niftymic/application/register_image.py +++ b/niftymic/application/register_image.py @@ -7,8 +7,12 @@ # \date October 2017 # -import numpy as np +import re import os +import numpy as np +import SimpleITK as sitk +# import nipype.interfaces.fsl +# import nipype.interfaces.niftyreg import pysitk.python_helper as ph import pysitk.simple_itk_helper as sitkh @@ -17,8 +21,52 @@ import niftymic.base.stack as st import niftymic.registration.flirt as regflirt import niftymic.registration.niftyreg as niftyreg +import niftymic.validation.image_similarity_evaluator as ise from niftymic.utilities.input_arparser import InputArgparser +from niftymic.definitions import REGEX_FILENAMES, DIR_TMP + + +## +# Gets the ap flip transform in case the 'stack' physical coordinate is aligned +# with the voxel space +# \date 2018-02-07 19:59:39+0000 +# +# \param stack The stack +# \param initializer_type The initializer type +# +# \return The ap flip transform. +# +def get_ap_flip_transform(path_to_image, initializer_type="GEOMETRY"): + image_sitk = sitk.ReadImage(path_to_image) + initial_transform = sitk.CenteredTransformInitializer( + image_sitk, + image_sitk, + sitk.Euler3DTransform(), + eval("sitk.CenteredTransformInitializerFilter.%s" % ( + initializer_type))) + + translation = np.array(initial_transform.GetFixedParameters()[0:3]) + + transform_sitk1 = sitk.Euler3DTransform() + transform_sitk1.SetTranslation(-translation) + + transform_sitk2 = sitk.Euler3DTransform() + transform_sitk2.SetRotation(0, 0, np.pi) + + transform_sitk3 = sitk.Euler3DTransform(transform_sitk1.GetInverse()) + + transform_sitk = sitkh.get_composite_sitk_euler_transform( + transform_sitk2, transform_sitk1) + transform_sitk = sitkh.get_composite_sitk_euler_transform( + transform_sitk3, transform_sitk) + + # sitkh.show_sitk_image(( + # [stack.sitk, + # sitkh.get_transformed_sitk_image(stack.sitk, transform_sitk)])) + + return transform_sitk + def main(): @@ -35,139 +83,295 @@ def main(): "defined by the fixed.", ) input_parser.add_fixed(required=True) - input_parser.add_moving( - required=True, - nargs="+", - help="Specify moving image to be warped to fixed space. " - "If multiple images are provided, all images will be transformed " - "uniformly according to the registration obtained for the first one." - ) + input_parser.add_moving(required=True) + input_parser.add_fixed_mask() + input_parser.add_moving_mask() input_parser.add_dir_output(required=True) - input_parser.add_dir_input() - input_parser.add_suffix_mask(default="_mask") + input_parser.add_dir_input_mc() input_parser.add_search_angle(default=180) input_parser.add_option( - option_string="--transform-only", + option_string="--initial-transform", + type=str, + help="Set initial transform to be used.", + default=None) + input_parser.add_option( + option_string="--test-ap-flip", type=int, - help="Turn on/off functionality to transform moving image(s) to fixed " - "image only, i.e. no resampling to fixed image space", - default=0) + help="Turn on/off functionality to run an additional registration " + "after an AP-flip. Seems to be more robust to find a better " + "registration outcome in general.", + default=1) input_parser.add_option( - option_string="--write-transform", + option_string="--use-flirt", type=int, - help="Turn on/off functionality to write registration transform", - default=0) + help="Turn on/off functionality to use FLIRT for the registration.", + default=1) + input_parser.add_option( + option_string="--use-regaladin", + type=int, + help="Turn on/off functionality to use RegAladin for the " + "registration.", + default=1) input_parser.add_verbose(default=0) + input_parser.add_log_config(default=1) args = input_parser.parse_args() input_parser.print_arguments(args) - use_reg_aladin_for_refinement = True + debug = 0 + + if args.log_config: + input_parser.log_config(os.path.abspath(__file__)) + + if not args.use_regaladin and not args.use_flirt: + raise IOError("Either RegAladin or FLIRT must be activated.") # --------------------------------Read Data-------------------------------- ph.print_title("Read Data") - data_reader = dr.MultipleImagesReader(args.moving, suffix_mask="_mask") - data_reader.read_data() - moving = data_reader.get_data() + fixed = st.Stack.from_filename( + file_path=args.fixed, + file_path_mask=args.fixed_mask, + extract_slices=False) + moving = st.Stack.from_filename( + file_path=args.moving, + file_path_mask=args.moving_mask, + extract_slices=False) + + if args.fixed_mask is not None: + use_fixed_mask = True + if args.moving_mask is not None: + use_moving_mask = True - data_reader = dr.MultipleImagesReader([args.fixed], suffix_mask="_mask") - data_reader.read_data() - fixed = data_reader.get_data()[0] + path_to_transform = os.path.join( + args.dir_output, "registration_transform_sitk.txt") + if args.initial_transform is not None: + transform_sitk = sitkh.read_transform_sitk(args.initial_transform) + else: + transform_sitk = sitk.AffineTransform(fixed.sitk.GetDimension()) + sitk.WriteTransform(transform_sitk, path_to_transform) + + path_to_output = os.path.join( + args.dir_output, + ph.append_to_filename(os.path.basename(args.moving), + "ResamplingToTemplateSpace")) # -------------------Register Reconstruction to Template------------------- ph.print_title("Register Reconstruction to Template") - # Define search angle ranges for FLIRT in all three dimensions - search_angles = ["-searchr%s -%d %d" % - (x, args.search_angle, args.search_angle) - for x in ["x", "y", "z"]] - search_angles = (" ").join(search_angles) - options_args = [] - options_args.append(search_angles) - # cost = "mutualinfo" - # options_args.append("-searchcost %s -cost %s" % (cost, cost)) - registration = regflirt.FLIRT( - fixed=moving[0], - moving=fixed, - # use_fixed_mask=True, - # use_moving_mask=True, # moving mask only seems to work for SB cases - registration_type="Rigid", - use_verbose=False, - options=(" ").join(options_args), - ) - ph.print_info("Run Registration (FLIRT) ... ", newline=False) - registration.run() - print("done") - transform_sitk = registration.get_registration_transform_sitk() - - if args.write_transform: - path_to_transform = os.path.join( - args.dir_output, "registration_transform_sitk.txt") - sitk.WriteTransform(transform_sitk, path_to_transform) - - # Apply rigidly transform to align reconstruction (moving) with template - # (fixed) - for m in moving: - m.update_motion_correction(transform_sitk) - - # Additionally, use RegAladin for more accurate alignment - # Rationale: FLIRT has better capture range, but RegAladin seems to - # find better alignment once it is within its capture range. - if use_reg_aladin_for_refinement: - registration = niftyreg.RegAladin( - fixed=m, - use_fixed_mask=True, - moving=fixed, - registration_type="Rigid", - use_verbose=False, - ) - ph.print_info("Run Registration (RegAladin) ... ", newline=False) - registration.run() - print("done") - transform2_sitk = registration.get_registration_transform_sitk() - m.update_motion_correction(transform2_sitk) - transform_sitk = sitkh.get_composite_sitk_affine_transform( - transform2_sitk, transform_sitk) - - if args.transform_only: - for m in moving: - m.write(args.dir_output, write_mask=False) - ph.exit() - - # Resample reconstruction (moving) to template space (fixed) - warped_moving = [m.get_resampled_stack(fixed.sitk, interpolator="Linear") - for m in moving] - - for wm in warped_moving: - wm.set_filename( - wm.get_filename() + "ResamplingToTemplateSpace") - - if args.verbose: - sitkh.show_stacks([fixed, wm], segmentation=fixed) - - # Write resampled reconstruction (moving) - wm.write(args.dir_output, write_mask=False) - - if args.dir_input is not None: - data_reader = dr.ImageSlicesDirectoryReader( - path_to_directory=args.dir_input, - suffix_mask=args.suffix_mask) - data_reader.read_data() - stacks = data_reader.get_data() - - for i, stack in enumerate(stacks): - stack.update_motion_correction(transform_sitk) - ph.print_info("Stack %d/%d: All slice transforms updated" % - (i+1, len(stacks))) - - # Write transformed slices - stack.write( - os.path.join(args.dir_output, "motion_correction"), - write_mask=True, - write_slices=True, - write_transforms=True, - suffix_mask=args.suffix_mask, - ) + if args.use_flirt: + path_to_transform_flirt = os.path.join(DIR_TMP, "transform_flirt.txt") + + # Convert SimpleITK into FLIRT transform + cmd = "simplereg_transform -sitk2flirt %s %s %s %s" % ( + path_to_transform, args.fixed, args.moving, path_to_transform_flirt) + ph.execute_command(cmd, verbose=False) + + # Define search angle ranges for FLIRT in all three dimensions + search_angles = ["-searchr%s -%d %d" % + (x, args.search_angle, args.search_angle) + for x in ["x", "y", "z"]] + + # flt = nipype.interfaces.fsl.FLIRT() + # flt.inputs.in_file = args.moving + # flt.inputs.reference = args.fixed + # if args.initial_transform is not None: + # flt.inputs.in_matrix_file = path_to_transform_flirt + # flt.inputs.out_matrix_file = path_to_transform_flirt + # # flt.inputs.output_type = "NIFTI_GZ" + # flt.inputs.out_file = path_to_output + # flt.inputs.args = "-dof 6" + # flt.inputs.args += " %s" % " ".join(search_angles) + # if args.moving_mask is not None: + # flt.inputs.in_weight = args.moving_mask + # if args.fixed_mask is not None: + # flt.inputs.ref_weight = args.fixed_mask + # ph.print_info("Run Registration (FLIRT) ... ", newline=False) + # flt.run() + # print("done") + + cmd_args = ["flirt"] + cmd_args.append("-in %s" % args.moving) + cmd_args.append("-ref %s" % args.fixed) + if args.initial_transform is not None: + cmd_args.append("-init %s" % path_to_transform_flirt) + cmd_args.append("-omat %s" % path_to_transform_flirt) + cmd_args.append("-out %s" % path_to_output) + cmd_args.append("-dof 6") + cmd_args.append((" ").join(search_angles)) + if args.moving_mask is not None: + cmd_args.append("-inweight %s" % args.moving_mask) + if args.fixed_mask is not None: + cmd_args.append("-refweight %s" % args.fixed_mask) + ph.print_info("Run Registration (FLIRT) ... ", newline=False) + ph.execute_command(" ".join(cmd_args), verbose=False) + print("done") + + # Convert FLIRT to SimpleITK transform + cmd = "simplereg_transform -flirt2sitk %s %s %s %s" % ( + path_to_transform_flirt, args.fixed, args.moving, path_to_transform) + ph.execute_command(cmd, verbose=False) + + if debug: + ph.show_niftis([args.fixed, path_to_output]) + + # Additionally, use RegAladin for more accurate alignment + # Rationale: FLIRT has better capture range, but RegAladin seems to + # find better alignment once it is within its capture range. + if args.use_regaladin: + path_to_transform_regaladin = os.path.join( + DIR_TMP, "transform_regaladin.txt") + + # Convert SimpleITK to RegAladin transform + cmd = "simplereg_transform -sitk2nreg %s %s" % ( + path_to_transform, path_to_transform_regaladin) + ph.execute_command(cmd, verbose=False) + + # nreg = nipype.interfaces.niftyreg.RegAladin() + # nreg.inputs.ref_file = args.fixed + # nreg.inputs.flo_file = args.moving + # nreg.inputs.res_file = path_to_output + # nreg.inputs.in_aff_file = path_to_transform_regaladin + # nreg.inputs.aff_file = path_to_transform_regaladin + # nreg.inputs.args = "-rigOnly -voff" + # if args.moving_mask is not None: + # nreg.inputs.fmask_file = args.moving_mask + # if args.fixed_mask is not None: + # nreg.inputs.rmask_file = args.fixed_mask + # ph.print_info("Run Registration (RegAladin) ... ", newline=False) + # nreg.run() + # print("done") + + cmd_args = ["reg_aladin"] + cmd_args.append("-ref %s" % args.fixed) + cmd_args.append("-flo %s" % args.moving) + cmd_args.append("-res %s" % path_to_output) + cmd_args.append("-inaff %s" % path_to_transform_regaladin) + cmd_args.append("-aff %s" % path_to_transform_regaladin) + cmd_args.append("-rigOnly") + cmd_args.append("-voff") + if args.moving_mask is not None: + cmd_args.append("-fmask %s" % args.moving_mask) + if args.fixed_mask is not None: + cmd_args.append("-rmask %s" % args.fixed_mask) + ph.print_info("Run Registration (RegAladin) ... ", newline=False) + ph.execute_command(" ".join(cmd_args), verbose=False) + print("done") + + # Convert RegAladin to SimpleITK transform + cmd = "simplereg_transform -nreg2sitk %s %s" % ( + path_to_transform_regaladin, path_to_transform) + ph.execute_command(cmd, verbose=False) + + if debug: + ph.show_niftis([args.fixed, path_to_output]) + + if args.test_ap_flip: + path_to_transform_flip = os.path.join(DIR_TMP, "transform_flip.txt") + path_to_output_flip = os.path.join(DIR_TMP, "output_flip.nii.gz") + + # Get AP-flip transform + transform_ap_flip_sitk = get_ap_flip_transform(args.fixed) + path_to_transform_flip_regaladin = os.path.join( + DIR_TMP, "transform_flip_regaladin.txt") + sitk.WriteTransform(transform_ap_flip_sitk, path_to_transform_flip) + + # Compose current transform with AP flip transform + cmd = "simplereg_transform -c %s %s %s" % ( + path_to_transform, path_to_transform_flip, path_to_transform_flip) + ph.execute_command(cmd, verbose=False) + + # Convert SimpleITK to RegAladin transform + cmd = "simplereg_transform -sitk2nreg %s %s" % ( + path_to_transform_flip, path_to_transform_flip_regaladin) + ph.execute_command(cmd, verbose=False) + + # nreg = nipype.interfaces.niftyreg.RegAladin() + # nreg.inputs.ref_file = args.fixed + # nreg.inputs.flo_file = args.moving + # nreg.inputs.res_file = path_to_output_flip + # nreg.inputs.in_aff_file = path_to_transform_flip_regaladin + # nreg.inputs.aff_file = path_to_transform_flip_regaladin + # nreg.inputs.args = "-rigOnly -voff" + # if args.moving_mask is not None: + # nreg.inputs.fmask_file = args.moving_mask + # if args.fixed_mask is not None: + # nreg.inputs.rmask_file = args.fixed_mask + # ph.print_info("Run Registration AP-flipped (RegAladin) ... ", + # newline=False) + # nreg.run() + # print("done") + + cmd_args = ["reg_aladin"] + cmd_args.append("-ref %s" % args.fixed) + cmd_args.append("-flo %s" % args.moving) + cmd_args.append("-res %s" % path_to_output_flip) + cmd_args.append("-inaff %s" % path_to_transform_flip_regaladin) + cmd_args.append("-aff %s" % path_to_transform_flip_regaladin) + cmd_args.append("-rigOnly") + cmd_args.append("-voff") + if args.moving_mask is not None: + cmd_args.append("-fmask %s" % args.moving_mask) + if args.fixed_mask is not None: + cmd_args.append("-rmask %s" % args.fixed_mask) + ph.print_info("Run Registration AP-flipped (RegAladin) ... ", + newline=False) + ph.execute_command(" ".join(cmd_args), verbose=False) + print("done") + + if debug: + ph.show_niftis([args.fixed, path_to_output, path_to_output_flip]) + + warped_moving = st.Stack.from_filename( + path_to_output, extract_slices=False) + warped_moving_flip = st.Stack.from_filename( + path_to_output_flip, extract_slices=False) + fixed = st.Stack.from_filename(args.fixed, args.fixed_mask) + + stacks = [warped_moving, warped_moving_flip] + image_similarity_evaluator = ise.ImageSimilarityEvaluator( + stacks=stacks, reference=fixed) + image_similarity_evaluator.compute_similarities() + similarities = image_similarity_evaluator.get_similarities() + + if similarities["NMI"][1] > similarities["NMI"][0]: + ph.print_info("AP-flipped outcome better") + + # Convert RegAladin to SimpleITK transform + cmd = "simplereg_transform -nreg2sitk %s %s" % ( + path_to_transform_flip_regaladin, path_to_transform) + ph.execute_command(cmd, verbose=False) + + # Copy better outcome + cmd = "cp -p %s %s" % (path_to_output_flip, path_to_output) + ph.execute_command(cmd, verbose=False) + + else: + ph.print_info("AP-flip does not improve outcome") + + if args.dir_input_mc is not None: + transform_sitk = sitkh.read_transform_sitk( + path_to_transform, inverse=1) + + if args.dir_input_mc.endswith("/"): + subdir_mc = args.dir_input_mc.split("/")[-2] + else: + subdir_mc = args.dir_input_mc.split("/")[-1] + dir_output_mc = os.path.join(args.dir_output, subdir_mc) + + ph.create_directory(dir_output_mc) + pattern = REGEX_FILENAMES + "[.]tfm" + p = re.compile(pattern) + trafos = [t for t in os.listdir(args.dir_input_mc) if p.match(t)] + for t in trafos: + path_to_input_transform = os.path.join(args.dir_input_mc, t) + path_to_output_transform = os.path.join(dir_output_mc, t) + t_sitk = sitkh.read_transform_sitk(path_to_input_transform) + t_sitk = sitkh.get_composite_sitk_affine_transform( + transform_sitk, t_sitk) + sitk.WriteTransform(t_sitk, path_to_output_transform) + + if args.verbose: + ph.show_niftis([args.fixed, path_to_output]) elapsed_time_total = ph.stop_timing(time_start) @@ -177,5 +381,6 @@ def main(): return 0 + if __name__ == '__main__': main() diff --git a/niftymic/application/run_reconstruction_parameter_study.py b/niftymic/application/run_reconstruction_parameter_study.py index 14a0bfce..cd2e064f 100755 --- a/niftymic/application/run_reconstruction_parameter_study.py +++ b/niftymic/application/run_reconstruction_parameter_study.py @@ -34,11 +34,11 @@ def main(): description="Script to study reconstruction parameters and their " "impact on the volumetric reconstruction quality.", ) - input_parser.add_dir_input() - input_parser.add_filenames() - input_parser.add_image_selection() - input_parser.add_dir_output(required=True) + input_parser.add_filenames(required=True) + input_parser.add_filenames_masks() input_parser.add_suffix_mask(default="_mask") + input_parser.add_dir_input_mc() + input_parser.add_dir_output(required=True) input_parser.add_reconstruction_space() input_parser.add_reference( help="Path to reference NIfTI image file. If given the volumetric " @@ -57,7 +57,7 @@ def main(): input_parser.add_alpha(default=0.01) input_parser.add_data_loss(default="linear") input_parser.add_data_loss_scale(default=1) - input_parser.add_log_script_execution(default=1) + input_parser.add_log_config(default=1) input_parser.add_verbose(default=1) # Range for parameter sweeps @@ -77,35 +77,18 @@ def main(): raise IOError("Either reference (--reference) or reconstruction space " "(--reconstruction-space) must be provided.") - # Write script execution call - if args.log_script_execution: - input_parser.write_performed_script_execution( - os.path.abspath(__file__)) + if args.log_config: + input_parser.log_config(os.path.abspath(__file__)) # --------------------------------Read Data-------------------------------- ph.print_title("Read Data") - # Neither '--dir-input' nor '--filenames' was specified - if args.filenames is not None and args.dir_input is not None: - raise IOError( - "Provide input by either '--dir-input' or '--filenames' " - "but not both together") - - # '--dir-input' specified - elif args.dir_input is not None: - data_reader = dr.ImageSlicesDirectoryReader( - path_to_directory=args.dir_input, - suffix_mask=args.suffix_mask, - image_selection=args.image_selection) - - # '--filenames' specified - elif args.filenames is not None: - data_reader = dr.MultipleImagesReader( - args.filenames, suffix_mask=args.suffix_mask) - - else: - raise IOError( - "Provide input by either '--dir-input' or '--filenames'") + data_reader = dr.MultipleImagesReader( + file_paths=args.filenames, + file_paths_masks=args.filenames_masks, + suffix_mask=args.suffix_mask, + dir_motion_correction=args.dir_input_mc, + ) data_reader.read_data() stacks = data_reader.get_data() diff --git a/niftymic/application/run_reconstruction_pipeline.py b/niftymic/application/run_reconstruction_pipeline.py new file mode 100644 index 00000000..687ab2f4 --- /dev/null +++ b/niftymic/application/run_reconstruction_pipeline.py @@ -0,0 +1,360 @@ +## +# \file run_reconstruction_pipeline.py +# \brief Script to execute entire reconstruction pipeline +# +# \author Michael Ebner (michael.ebner.14@ucl.ac.uk) +# \date October 2017 +# + +import os +import re +import numpy as np + +import pysitk.python_helper as ph + +import niftymic.validation.simulate_stacks_from_reconstruction as \ + simulate_stacks_from_reconstruction +import niftymic.validation.evaluate_simulated_stack_similarity as \ + evaluate_simulated_stack_similarity +import niftymic.validation.show_evaluated_simulated_stack_similarity as \ + show_evaluated_simulated_stack_similarity +import niftymic.validation.export_side_by_side_simulated_vs_original_slice_comparison as \ + export_side_by_side_simulated_vs_original_slice_comparison +from niftymic.utilities.input_arparser import InputArgparser + +from niftymic.definitions import DIR_TEMPLATES + + + +def main(): + + time_start = ph.start_timing() + + np.set_printoptions(precision=3) + + input_parser = InputArgparser( + description="Run reconstruction pipeline including " + "(i) preprocessing (bias field correction + intensity correction), " + "(ii) volumetric reconstruction in subject space, " + "and (iii) volumetric reconstruction in template space.", + ) + input_parser.add_filenames(required=True) + input_parser.add_target_stack(required=False) + input_parser.add_suffix_mask(default="_mask") + input_parser.add_dir_output(required=True) + input_parser.add_alpha(default=0.01) + input_parser.add_verbose(default=0) + input_parser.add_gestational_age(required=False) + input_parser.add_prefix_output(default="") + input_parser.add_search_angle(default=180) + input_parser.add_multiresolution(default=0) + input_parser.add_log_config(default=1) + input_parser.add_dir_input_templates(default=DIR_TEMPLATES) + input_parser.add_isotropic_resolution() + input_parser.add_reference() + input_parser.add_reference_mask() + input_parser.add_bias_field_correction(default=1) + input_parser.add_intensity_correction(default=1) + input_parser.add_iter_max(default=10) + input_parser.add_two_step_cycles(default=3) + input_parser.add_option( + option_string="--run-bias-field-correction", + type=int, + help="Turn on/off bias field correction. " + "If off, it is assumed that this step was already performed", + default=1) + input_parser.add_option( + option_string="--run-recon-subject-space", + type=int, + help="Turn on/off reconstruction in subject space. " + "If off, it is assumed that this step was already performed", + default=1) + input_parser.add_option( + option_string="--run-recon-template-space", + type=int, + help="Turn on/off reconstruction in template space. " + "If off, it is assumed that this step was already performed", + default=1) + input_parser.add_option( + option_string="--run-data-vs-simulated-data", + type=int, + help="Turn on/off comparison of data vs data simulated from the " + "obtained volumetric reconstruction. " + "If off, it is assumed that this step was already performed", + default=1) + input_parser.add_outlier_rejection(default=1) + + args = input_parser.parse_args() + input_parser.print_arguments(args) + + if args.log_config: + input_parser.log_config(os.path.abspath(__file__)) + + prefix_bias = "N4ITK_" + prefix_ic = "IC_" + dir_output_preprocessing = os.path.join( + args.dir_output, "preprocessing") + dir_output_recon_subject_space = os.path.join( + args.dir_output, "recon_subject_space") + dir_output_recon_template_space = os.path.join( + args.dir_output, "recon_template_space") + dir_output_data_vs_simulatd_data = os.path.join( + args.dir_output, "data_vs_simulated_data") + + if args.run_recon_template_space and args.gestational_age is None: + raise IOError("Gestational age must be set in order to pick the " + "right template") + + if args.target_stack is None: + target_stack = args.filenames[0] + else: + target_stack = args.target_stack + + if args.bias_field_correction and args.run_bias_field_correction: + cmd_args = [] + cmd_args.append("--filenames %s" % (" ").join(args.filenames)) + cmd_args.append("--dir-output %s" % dir_output_preprocessing) + cmd_args.append("--prefix-output %s" % prefix_bias) + cmd_args.append("--suffix-mask %s" % args.suffix_mask) + # cmd_args.append("--verbose %d" % args.verbose) + cmd = "niftymic_correct_bias_field %s" % (" ").join(cmd_args) + time_start_bias = ph.start_timing() + exit_code = ph.execute_command(cmd) + if exit_code != 0: + raise RuntimeError("Bias field correction failed") + elapsed_time_bias = ph.stop_timing(time_start_bias) + filenames = [os.path.join(dir_output_preprocessing, "%s%s" % ( + prefix_bias, os.path.basename(f))) + for f in args.filenames] + elif args.bias_field_correction and not args.run_bias_field_correction: + elapsed_time_bias = ph.get_zero_time() + filenames = [os.path.join(dir_output_preprocessing, "%s%s" % ( + prefix_bias, os.path.basename(f))) + for f in args.filenames] + else: + elapsed_time_bias = ph.get_zero_time() + filenames = args.filenames + + if args.run_recon_subject_space: + + target_stack_index = args.filenames.index(target_stack) + + cmd_args = [] + cmd_args.append("--filenames %s" % (" ").join(filenames)) + cmd_args.append("--multiresolution %d" % args.multiresolution) + cmd_args.append("--target-stack-index %d" % target_stack_index) + cmd_args.append("--dir-output %s" % dir_output_recon_subject_space) + cmd_args.append("--suffix-mask %s" % args.suffix_mask) + cmd_args.append("--intensity-correction %d" % + args.intensity_correction) + cmd_args.append("--alpha %s" % args.alpha) + cmd_args.append("--iter-max %d" % args.iter_max) + cmd_args.append("--two-step-cycles %d" % args.two_step_cycles) + cmd_args.append("--outlier-rejection %d" % + args.outlier_rejection) + cmd_args.append("--verbose %d" % args.verbose) + if args.isotropic_resolution is not None: + cmd_args.append("--isotropic-resolution %f" % + args.isotropic_resolution) + if args.reference is not None: + cmd_args.append("--reference %s" % args.reference) + if args.reference_mask is not None: + cmd_args.append("--reference-mask %s" % args.reference_mask) + cmd = "niftymic_reconstruct_volume %s" % (" ").join(cmd_args) + time_start_volrec = ph.start_timing() + exit_code = ph.execute_command(cmd) + if exit_code != 0: + raise RuntimeError("Reconstruction in subject space failed") + elapsed_time_volrec = ph.stop_timing(time_start_volrec) + else: + elapsed_time_volrec = ph.get_zero_time() + + if args.run_recon_template_space: + # register recon to template space + template = os.path.join( + args.dir_input_templates, + "STA%d.nii.gz" % args.gestational_age) + template_mask = os.path.join( + args.dir_input_templates, + "STA%d_mask.nii.gz" % args.gestational_age) + pattern = "[a-zA-Z0-9_]+(stacks)[a-zA-Z0-9_]+(.nii.gz)" + p = re.compile(pattern) + reconstruction = [ + os.path.join( + dir_output_recon_subject_space, p.match(f).group(0)) + for f in os.listdir(dir_output_recon_subject_space) + if p.match(f)][0] + reconstruction = re.sub( + "%s.nii.gz" % args.suffix_mask, ".nii.gz", reconstruction) + cmd_args = [] + cmd_args.append("--fixed %s" % template) + cmd_args.append("--moving %s" % reconstruction) + cmd_args.append("--fixed-mask %s" % template_mask) + cmd_args.append("--moving-mask %s" % + ph.append_to_filename(reconstruction, args.suffix_mask)) + cmd_args.append("--dir-input-mc %s" % os.path.join( + dir_output_recon_subject_space, + "motion_correction")) + cmd_args.append("--dir-output %s" % dir_output_recon_template_space) + cmd_args.append("--verbose %s" % args.verbose) + cmd = "niftymic_register_image %s" % (" ").join(cmd_args) + exit_code = ph.execute_command(cmd) + if exit_code != 0: + raise RuntimeError("Registration to template space failed") + + # reconstruct volume in template space + # pattern = "[a-zA-Z0-9_.]+(ResamplingToTemplateSpace.nii.gz)" + # p = re.compile(pattern) + # reconstruction_space = [ + # os.path.join(dir_output_recon_template_space, p.match(f).group(0)) + # for f in os.listdir(dir_output_recon_template_space) + # if p.match(f)][0] + + dir_input_mc = os.path.join( + dir_output_recon_template_space, "motion_correction") + cmd_args = [] + cmd_args.append("--filenames %s" % (" ").join(filenames)) + cmd_args.append("--dir-input-mc %s" % dir_input_mc) + cmd_args.append("--dir-output %s" % dir_output_recon_template_space) + cmd_args.append("--reconstruction-space %s" % template) + cmd_args.append("--iter-max %d" % args.iter_max) + cmd_args.append("--alpha %s" % args.alpha) + cmd_args.append("--suffix-mask %s" % args.suffix_mask) + + cmd = "niftymic_reconstruct_volume_from_slices %s" % \ + (" ").join(cmd_args) + exit_code = ph.execute_command(cmd) + if exit_code != 0: + raise RuntimeError("Reconstruction in template space failed") + + pattern = "[a-zA-Z0-9_.]+(stacks[0-9]+).*(.nii.gz)" + p = re.compile(pattern) + reconstruction = { + p.match(f).group(1): + os.path.join( + dir_output_recon_template_space, p.match(f).group(0)) + for f in os.listdir(dir_output_recon_template_space) + if p.match(f) and not p.match(f).group(0).endswith( + "ResamplingToTemplateSpace.nii.gz")} + key = list(reconstruction.keys())[0] + path_to_recon = reconstruction[key] + path_to_recon = re.sub( + "%s.nii.gz" % args.suffix_mask, ".nii.gz", path_to_recon) + + # Copy SRR to output directory + output = "%sSRR_%s_GW%d.nii.gz" % ( + args.prefix_output, key, args.gestational_age) + path_to_output = os.path.join(args.dir_output, output) + cmd = "cp -p %s %s" % (path_to_recon, path_to_output) + exit_code = ph.execute_command(cmd) + if exit_code != 0: + raise RuntimeError("Copy of SRR to output directory failed") + + # Multiply template mask with reconstruction + cmd_args = [] + cmd_args.append("--filename %s" % path_to_output) + cmd_args.append("--gestational-age %s" % args.gestational_age) + cmd_args.append("--verbose %s" % args.verbose) + cmd_args.append("--dir-input-templates %s " % args.dir_input_templates) + cmd = "niftymic_multiply_stack_with_mask %s" % ( + " ").join(cmd_args) + exit_code = ph.execute_command(cmd) + if exit_code != 0: + raise RuntimeError("SRR brain masking failed") + + else: + elapsed_time_template = ph.get_zero_time() + + if args.run_data_vs_simulated_data: + + dir_input_mc = os.path.join( + dir_output_recon_template_space, "motion_correction") + + pattern = "[a-zA-Z0-9_.]+(stacks[0-9]+).*(.nii.gz)" + # pattern = "Masked_[a-zA-Z0-9_.]+(stacks[0-9]+).*(.nii.gz)" + p = re.compile(pattern) + reconstruction = { + p.match(f).group(1): + os.path.join( + dir_output_recon_template_space, p.match(f).group(0)) + for f in os.listdir(dir_output_recon_template_space) + if p.match(f) and not p.match(f).group(0).endswith( + "ResamplingToTemplateSpace.nii.gz")} + key = list(reconstruction.keys())[0] + path_to_recon = reconstruction[key] + path_to_recon = re.sub( + "%s.nii.gz" % args.suffix_mask, ".nii.gz", path_to_recon) + + # Get simulated/projected slices + cmd_args = [] + cmd_args.append("--filenames %s" % (" ").join(filenames)) + cmd_args.append("--dir-input-mc %s" % dir_input_mc) + cmd_args.append("--dir-output %s" % dir_output_data_vs_simulatd_data) + cmd_args.append("--reconstruction %s" % path_to_recon) + cmd_args.append("--copy-data 1") + cmd_args.append("--suffix-mask %s" % args.suffix_mask) + # cmd_args.append("--verbose %s" % args.verbose) + exe = os.path.abspath(simulate_stacks_from_reconstruction.__file__) + cmd = "python %s %s" % (exe, (" ").join(cmd_args)) + exit_code = ph.execute_command(cmd) + if exit_code != 0: + raise RuntimeError("SRR slice projections failed") + + filenames_simulated = [ + os.path.join(dir_output_data_vs_simulatd_data, os.path.basename(f)) + for f in filenames] + + dir_output_evaluation = os.path.join( + dir_output_data_vs_simulatd_data, "evaluation") + dir_output_figures = os.path.join( + dir_output_data_vs_simulatd_data, "figures") + dir_output_side_by_side = os.path.join( + dir_output_figures, "side-by-side") + + # Evaluate slice similarities to ground truth + cmd_args = [] + cmd_args.append("--filenames %s" % (" ").join(filenames_simulated)) + cmd_args.append("--suffix-mask %s" % args.suffix_mask) + cmd_args.append("--measures NCC SSIM") + cmd_args.append("--dir-output %s" % dir_output_evaluation) + exe = os.path.abspath(evaluate_simulated_stack_similarity.__file__) + cmd = "python %s %s" % (exe, (" ").join(cmd_args)) + exit_code = ph.execute_command(cmd) + if exit_code != 0: + raise RuntimeError("Evaluation of slice similarities failed") + + # Generate figures showing the quantitative comparison + cmd_args = [] + cmd_args.append("--dir-input %s" % dir_output_evaluation) + cmd_args.append("--dir-output %s" % dir_output_figures) + exe = os.path.abspath( + show_evaluated_simulated_stack_similarity.__file__) + cmd = "python %s %s" % (exe, (" ").join(cmd_args)) + exit_code = ph.execute_command(cmd) + if exit_code != 0: + ph.print_warning("Visualization of slice similarities failed") + + # Generate pdfs showing all the side-by-side comparisons + cmd_args = [] + cmd_args.append("--filenames %s" % (" ").join(filenames_simulated)) + cmd_args.append("--dir-output %s" % dir_output_side_by_side) + exe = os.path.abspath( + export_side_by_side_simulated_vs_original_slice_comparison.__file__) + cmd = "python %s %s" % (exe, (" ").join(cmd_args)) + exit_code = ph.execute_command(cmd) + if exit_code != 0: + raise RuntimeError("Generation of PDF overview failed") + + ph.print_title("Summary") + print("Computational Time for Bias Field Correction: %s" % + elapsed_time_bias) + print("Computational Time for Volumetric Reconstruction: %s" % + elapsed_time_volrec) + print("Computational Time for Pipeline: %s" % + ph.stop_timing(time_start)) + + return 0 + + +if __name__ == '__main__': + main() diff --git a/niftymic/base/data_reader.py b/niftymic/base/data_reader.py index 1f88a04b..b416b959 100644 --- a/niftymic/base/data_reader.py +++ b/niftymic/base/data_reader.py @@ -12,7 +12,7 @@ import numpy as np import os import re -# Import libraries +import six from abc import ABCMeta, abstractmethod import niftymic.base.stack as st @@ -23,12 +23,11 @@ from niftymic.definitions import REGEX_FILENAMES from niftymic.definitions import REGEX_FILENAME_EXTENSIONS + ## # DataReader is an abstract class to read data. # \date 2017-07-12 11:38:07+0100 # - - class DataReader(object): __metaclass__ = ABCMeta @@ -119,9 +118,11 @@ def read_data(self): p = re.compile(pattern) p_mask = re.compile(pattern_mask) - # Exclude potential mask filenames - # TODO: If folder contains A.nii and A.nii.gz that ambiguity will not - # be detected + # TODO: + # - If folder contains A.nii and A.nii.gz that ambiguity will not + # be detected + # - exclude potential mask filenames + # - hidden files are not excluded dic_filenames = {p.match(f).group(1): p.match(f).group(0) for f in os.listdir(abs_path_to_directory) if p.match(f) and not p_mask.match(f)} @@ -167,24 +168,38 @@ class MultipleImagesReader(ImageDataReader): # masks. # \date 2017-07-11 19:04:25+0100 # - # \param self The object - # \param file_paths The paths to filenames as list of strings, - # e.g. ["A.nii.gz", "B.nii", "C.nii.gz"] - # \param suffix_mask extension of stack filename as string - # indicating associated mask, e.g. "_mask" for - # "A_mask.nii". - # \param extract_slices Boolean to indicate whether given 3D image - # shall be split into its slices along the - # k-direction. + # \param self The object + # \param file_paths The paths to filenames as list of strings, + # e.g. ["A.nii.gz", "B.nii", "C.nii.gz"] + # \param file_paths_masks The paths to the filename masks as list of + # strings. If given 'suffix_mask' is ignored. + # It is assumed the the sequence matches the + # filename order. + # \param suffix_mask extension of stack filename as string + # indicating associated mask, e.g. "_mask" + # for "A_mask.nii". + # \param extract_slices Boolean to indicate whether given 3D image + # shall be split into its slices along the + # k-direction. # - def __init__(self, file_paths, suffix_mask="_mask", extract_slices=True): + def __init__(self, + file_paths, + file_paths_masks=None, + suffix_mask="_mask", + extract_slices=True, + dir_motion_correction=None, + prefix_slice="_slice", + ): super(self.__class__, self).__init__() # Get list of paths to image self._file_paths = file_paths + self._file_paths_masks = file_paths_masks self._suffix_mask = suffix_mask + self._dir_motion_correction = dir_motion_correction self._extract_slices = extract_slices + self._prefix_slice = prefix_slice ## # Reads the data of multiple images. @@ -192,48 +207,130 @@ def __init__(self, file_paths, suffix_mask="_mask", extract_slices=True): # def read_data(self): - self._stacks = [None] * len(self._file_paths) - - for i, file_path in enumerate(self._file_paths): + self._check_input() - # Build absolute path to directory of image - path_to_directory = os.path.dirname(file_path) - filename = os.path.basename(file_path) + stacks = [None] * len(self._file_paths) - if not ph.directory_exists(path_to_directory): - raise exceptions.DirectoryNotExistent(path_to_directory) - abs_path_to_directory = os.path.abspath(path_to_directory) + for i, file_path in enumerate(self._file_paths): - # Get absolute path mask to image - pattern = "(" + REGEX_FILENAMES + \ - ")[.]" + REGEX_FILENAME_EXTENSIONS - p = re.compile(pattern) - # filename = [p.match(f).group(1) if p.match(file_path)][0] - if not file_path.endswith(tuple(ALLOWED_EXTENSIONS)): - raise IOError("Input image type not correct. Allowed types %s" - % "(" + (", or ").join(ALLOWED_EXTENSIONS) + ")") - - # Strip extension from filename to find associated mask - filename = [re.sub("." + ext, "", filename) - for ext in ALLOWED_EXTENSIONS - if file_path.endswith(ext)][0] - pattern_mask = filename + self._suffix_mask + "[.]" + \ - REGEX_FILENAME_EXTENSIONS - p_mask = re.compile(pattern_mask) - filename_mask = [p_mask.match(f).group(0) - for f in os.listdir(abs_path_to_directory) - if p_mask.match(f)] - - if len(filename_mask) == 0: - abs_path_mask = None + if self._file_paths_masks is None: + file_path_mask = self._get_path_to_potential_mask(file_path) else: - abs_path_mask = os.path.join(abs_path_to_directory, - filename_mask[0]) - self._stacks[i] = st.Stack.from_filename( + if i < len(self._file_paths_masks): + file_path_mask = self._file_paths_masks[i] + else: + file_path_mask = None + + stacks[i] = st.Stack.from_filename( file_path, - abs_path_mask, + file_path_mask, extract_slices=self._extract_slices) + if self._dir_motion_correction is not None: + if not ph.directory_exists(self._dir_motion_correction): + raise exceptions.DirectoryNotExistent( + self._dir_motion_correction) + abs_path_to_directory = os.path.abspath( + self._dir_motion_correction) + stack_name = ph.strip_filename_extension( + os.path.basename(file_path))[0] + + path_to_stack_transform = os.path.join( + abs_path_to_directory, "%s.tfm" % stack_name) + if ph.file_exists(path_to_stack_transform): + transform_stack_sitk = sitkh.read_transform_sitk( + path_to_stack_transform) + transform_stack_sitk_inv = sitkh.read_transform_sitk( + path_to_stack_transform, inverse=True) + stacks[i].update_motion_correction(transform_stack_sitk) + ph.print_info( + "Stack '%s': Stack position updated" % stack_name) + else: + transform_stack_sitk_inv = sitk.Euler3DTransform() + + pattern_trafo_slices = stack_name + self._prefix_slice + \ + "([0-9]+)[.]tfm" + p = re.compile(pattern_trafo_slices) + dic_slice_transforms = { + int(p.match(f).group(1)): os.path.join( + abs_path_to_directory, p.match(f).group(0)) + for f in os.listdir(abs_path_to_directory) if p.match(f) + } + slices = stacks[i].get_slices() + for i_slice in range(stacks[i].get_number_of_slices()): + if i_slice in dic_slice_transforms.keys(): + transform_slice_sitk = sitkh.read_transform_sitk( + dic_slice_transforms[i_slice]) + transform_slice_sitk = sitkh.get_composite_sitk_affine_transform( + transform_slice_sitk, transform_stack_sitk_inv) + slices[i_slice].update_motion_correction( + transform_slice_sitk) + else: + stacks[i].delete_slice(i_slice) + if stacks[i].get_number_of_slices() == 0: + ph.print_info( + "Stack '%s' removed as all slices were deleted" % + stack_name) + stacks[i] = None + ph.print_info("Stack '%s': Slice positions updated" % stack_name) + + self._stacks = [s for s in stacks if s is not None] + + def _check_input(self): + if type(self._file_paths) is not list: + raise IOError("file_paths must be provided as list") + + if self._file_paths_masks is not None: + if type(self._file_paths_masks) is not list: + raise IOError("file_paths_masks must be provided as list") + + ## + # Gets the path to potential mask for a given file_path. + # \date 2018-01-25 12:34:23+0000 + # + # \param self The object + # \param file_path The file path + # + # \return The path the mask associated with the file or None in case no + # mask was found. + # + def _get_path_to_potential_mask(self, file_path): + # Build absolute path to directory of image + path_to_directory = os.path.dirname(file_path) + filename = os.path.basename(file_path) + + if not ph.directory_exists(path_to_directory): + raise exceptions.DirectoryNotExistent(path_to_directory) + abs_path_to_directory = os.path.abspath(path_to_directory) + + # Get absolute path mask to image + pattern = "(" + REGEX_FILENAMES + \ + ")[.]" + REGEX_FILENAME_EXTENSIONS + p = re.compile(pattern) + # filename = [p.match(f).group(1) if p.match(file_path)][0] + if not file_path.endswith(tuple(ALLOWED_EXTENSIONS)): + raise IOError("Input image type not correct. Allowed types %s" + % "(" + (", or ").join(ALLOWED_EXTENSIONS) + ")") + + # Strip extension from filename to find associated mask + filename = [re.sub("." + ext, "", filename) + for ext in ALLOWED_EXTENSIONS + if file_path.endswith(ext)][0] + pattern_mask = filename + self._suffix_mask + "[.]" + \ + REGEX_FILENAME_EXTENSIONS + p_mask = re.compile(pattern_mask) + filename_mask = [p_mask.match(f).group(0) + for f in os.listdir(abs_path_to_directory) + if p_mask.match(f)] + + if len(filename_mask) == 0: + abs_path_mask = None + else: + abs_path_mask = os.path.join(abs_path_to_directory, + filename_mask[0]) + + return abs_path_mask + ## # ImageSlicesDirectoryReader reads multiple stacks and their associated @@ -393,16 +490,91 @@ def __init__(self): DataReader.__init__(self) self._transforms_sitk = None + # Third line in *.tfm file contains information on the transform type + self._transform_type = { + "Euler3DTransform_double_3_3": sitk.Euler3DTransform, + "AffineTransform_double_3_3": sitk.AffineTransform, + } + def get_data(self): return self._transforms_sitk + def _get_sitk_transform_from_filepath(self, path_to_sitk_transform): + # Read transform as type sitk.Transform + transform_sitk = sitk.ReadTransform(path_to_sitk_transform) + + # Convert transform to respective type, e.g. Euler, Affine etc + # Third line in *.tfm file contains information on the transform type + transform_type = open(path_to_sitk_transform).readlines()[2] + transform_type = re.sub("\n", "", transform_type) + transform_type = transform_type.split(" ")[1] + transform_sitk = self._transform_type[transform_type](transform_sitk) + + return transform_sitk + +## +# Reads slice transformations stored in the format 'filename_slice#.tfm'. +# +# Rationale: Read only slice transformations associated with +# 'motion_correction' export achieved by the volumetric reconstruction +# algorithm +# \date 2018-01-31 19:16:00+0000 +# +class SliceTransformationDirectoryReader(TransformationDataReader): + + def __init__(self, directory, suffix_slice="_slice"): + TransformationDataReader.__init__(self) + self._directory = directory + self._suffix_slice = suffix_slice + + def read_data(self): + + if not ph.directory_exists(self._directory): + raise exceptions.DirectoryNotExistent(self._directory) + + # Create absolute path for directory + directory = os.path.abspath(self._directory) + + pattern = "(" + REGEX_FILENAMES + \ + ")%s([0-9]+)[.]tfm" % self._suffix_slice + p = re.compile(pattern) + + dic_tmp = { + (p.match(f).group(1), int(p.match(f).group(2))): + os.path.join(directory, p.match(f).group(0)) + for f in os.listdir(directory) if p.match(f) + } + fnames = list(set([k[0] for k in dic_tmp.keys()])) + self._transforms_sitk = {fname: {} for fname in fnames} + for (fname, slice_number), path in six.iteritems(dic_tmp): + self._transforms_sitk[fname][slice_number] = \ + self._get_sitk_transform_from_filepath(path) + + +## +# Reads all transformations in a given directory and stores them in an ordered +# list +# \date 2018-01-31 19:34:52+0000 +# class TransformationDirectoryReader(TransformationDataReader): def __init__(self, directory): TransformationDataReader.__init__(self) self._directory = directory + ## + # Reads all transformations in a given directory and stores them in an + # ordered list + # \date 2018-01-31 19:32:45+0000 + # + # \param self The object + # \param extension The extension + # \post self._transforms_sitk contains transformations as list of + # sitk.Transformation objects + # + # \return { description_of_the_return_value } + # def read_data(self, extension="tfm"): pattern = REGEX_FILENAMES + "[.]" + extension p = re.compile(pattern) @@ -418,28 +590,26 @@ def read_data(self, extension="tfm"): self._transforms_sitk = transforms_reader.get_data() +## +# Reads multiple transformations and store them as lists +# \date 2018-01-31 19:33:51+0000 +# class MultipleTransformationsReader(TransformationDataReader): def __init__(self, file_paths): super(self.__class__, self).__init__() self._file_paths = file_paths - # Third line in *.tfm file contains information on the transform type - self._transform_type = { - "Euler3DTransform_double_3_3": sitk.Euler3DTransform, - "AffineTransform_double_3_3": sitk.AffineTransform, - } - + ## + # Reads multiple transformations and store them as lists + # \date 2018-01-31 19:32:45+0000 + # + # \param self The object + # \post self._transforms_sitk contains transformations as list of + # sitk.Transformation objects + # def read_data(self): - self._transforms_sitk = [None] * len(self._file_paths) - - for i in range(len(self._file_paths)): - # Read transform as type sitk.Transform - transform_sitk = sitk.ReadTransform(self._file_paths[i]) - - # Convert transform to respective type, e.g. Euler, Affine etc - transform_type = open(self._file_paths[i]).readlines()[2] - transform_type = re.sub("\n", "", transform_type) - transform_type = transform_type.split(" ")[1] - self._transforms_sitk[i] = self._transform_type[transform_type]( - transform_sitk) + self._transforms_sitk = [ + self._get_sitk_transform_from_filepath(self._file_paths[i]) + for i in range(len(self._file_paths)) + ] diff --git a/niftymic/base/slice.py b/niftymic/base/slice.py index 803aac0d..0d7303f9 100644 --- a/niftymic/base/slice.py +++ b/niftymic/base/slice.py @@ -51,7 +51,14 @@ def from_sitk_image(cls, slice_sitk, slice_number, filename="unknown", slice_sit # Append masks (if provided) if slice_sitk_mask is not None: slice.sitk_mask = slice_sitk_mask - slice.itk_mask = sitkh.get_itk_from_sitk_image(slice_sitk_mask) + try: + # ensure mask occupies the same physical space + slice.sitk_mask.CopyInformation(slice.sitk) + except RuntimeError as e: + raise IOError( + "Given image and its mask do not occupy the same space: %s" % + e.message) + slice.itk_mask = sitkh.get_itk_from_sitk_image(slice.sitk_mask) else: slice.sitk_mask = slice._generate_identity_mask() slice.itk_mask = sitkh.get_itk_from_sitk_image(slice.sitk_mask) @@ -61,7 +68,8 @@ def from_sitk_image(cls, slice_sitk, slice_number, filename="unknown", slice_sit # HACK (for current Slice-to-Volume Registration) # See class SliceToVolumeRegistration # slice._sitk_upsampled = slice._get_upsampled_isotropic_resolution_slice(slice_sitk) - # slice._itk_upsampled = sitkh.get_itk_from_sitk_image(slice._sitk_upsampled) + # slice._itk_upsampled = + # sitkh.get_itk_from_sitk_image(slice._sitk_upsampled) # if slice_sitk_mask is not None: # slice._sitk_mask_upsampled = slice._get_upsampled_isotropic_resolution_slice(slice_sitk_mask) @@ -107,7 +115,7 @@ def from_filename(cls, file_path, slice_number, file_path_mask=None, verbose=Fal slice._slice_number = slice_number # Append stacks as SimpleITK and ITK Image objects - slice.sitk = sitk.ReadImage(file_path, sitk.sitkFloat64) + slice.sitk = sitkh.read_nifti_image_sitk(file_path, sitk.sitkFloat64) slice.itk = sitkh.get_itk_from_sitk_image(slice.sitk) # Append masks (if provided) @@ -120,7 +128,15 @@ def from_filename(cls, file_path, slice_number, file_path_mask=None, verbose=Fal else: if not ph.file_exists(file_path_mask): raise exceptions.FileNotExistent(file_path_mask) - slice.sitk_mask = sitk.ReadImage(file_path_mask, sitk.sitkUInt8) + slice.sitk_mask = sitkh.read_nifti_image_sitk( + file_path_mask, sitk.sitkUInt8) + try: + # ensure mask occupies the same physical space + slice.sitk_mask.CopyInformation(slice.sitk) + except RuntimeError as e: + raise IOError( + "Given image and its mask do not occupy the same space: %s" % + e.message) slice.itk_mask = sitkh.get_itk_from_sitk_image(slice.sitk_mask) @@ -162,7 +178,8 @@ def from_slice(cls, slice_to_copy): slice._slice_number = slice_to_copy.get_slice_number() slice._dir_input = slice_to_copy.get_directory() - # slice._history_affine_transforms, slice._history_motion_corrections = slice_to_copy.get_registration_history() + # slice._history_affine_transforms, slice._history_motion_corrections = + # slice_to_copy.get_registration_history() # Store current affine transform of image slice._affine_transform_sitk = sitkh.get_sitk_affine_transform_from_sitk_image( @@ -214,7 +231,9 @@ def update_motion_correction(self, affine_transform_sitk): # self._history_rigid_motion_estimates.append(current_rigid_motion_estimate) # ## New affine transform of slice after rigid motion correction - # affine_transform = sitkh.get_composite_sitk_affine_transform(rigid_transform_sitk, self._affine_transform_sitk) + # affine_transform = + # sitkh.get_composite_sitk_affine_transform(rigid_transform_sitk, + # self._affine_transform_sitk) # ## Update affine transform of slice, i.e. change image origin and direction in physical space # self._update_affine_transform(affine_transform) @@ -260,6 +279,13 @@ def get_registration_history(self): motion_corrections = list(self._history_motion_corrections) return affine_transforms, motion_corrections + def set_registration_history(self, registration_history): + affine_transform_sitk = registration_history[0][-1] + self._update_affine_transform(affine_transform_sitk) + + self._history_affine_transforms = [a for a in registration_history[0]] + self._history_motion_corrections = [t for t in registration_history[1]] + # Display slice with external viewer (ITK-Snap) # \param[in] show_segmentation display slice with or without associated segmentation (default=0) def show(self, show_segmentation=0, label=None, viewer=VIEWER, verbose=True): @@ -284,7 +310,15 @@ def show(self, show_segmentation=0, label=None, viewer=VIEWER, verbose=True): # - affine transformation describing physical space position of slice # \param[in] directory string specifying where the output will be written to (default="/tmp/") # \param[in] filename string specifyig the filename. If not given, filename of parent stack is used - def write(self, directory, filename=None, write_transform=False, suffix_mask="_mask", prefix_slice="_slice"): + def write(self, + directory, + filename=None, + write_slice=True, + write_transform=True, + suffix_mask="_mask", + prefix_slice="_slice", + use_float32=True, + ): # Create directory if not existing ph.create_directory(directory) @@ -299,24 +333,31 @@ def write(self, directory, filename=None, write_transform=False, suffix_mask="_m full_file_name = os.path.join(directory, filename_out) # Write slice and affine transform - sitkh.write_nifti_image_sitk(self.sitk, full_file_name + ".nii.gz") + if write_slice: + if use_float32: + image_sitk = sitk.Cast(self.sitk, sitk.sitkFloat32) + else: + image_sitk = self.sitk + sitkh.write_nifti_image_sitk( + image_sitk, full_file_name + ".nii.gz") + + # Write mask to specified location if given + if self.sitk_mask is not None: + nda = sitk.GetArrayFromImage(self.sitk_mask) + + # Write mask if it does not consist of only ones + if not np.all(nda): + sitkh.write_nifti_image_sitk(self.sitk_mask, full_file_name + + "%s.nii.gz" % (suffix_mask)) if write_transform: sitk.WriteTransform( # self.get_affine_transform(), self.get_motion_correction_transform(), full_file_name + ".tfm") - # Write mask to specified location if given - if self.sitk_mask is not None: - nda = sitk.GetArrayFromImage(self.sitk_mask) - - # Write mask if it does not consist of only ones - if not np.all(nda): - sitkh.write_nifti_image_sitk(self.sitk_mask, full_file_name + - "%s.nii.gz" % (suffix_mask)) - # print("Slice %r of stack %s was successfully written to %s" %(self._slice_number, self._filename, full_file_name)) - # print("Transformation of slice %r of stack %s was successfully written to %s" %(self._slice_number, self._filename, full_file_name)) + # print("Transformation of slice %r of stack %s was successfully + # written to %s" %(self._slice_number, self._filename, full_file_name)) # Update slice with new affine transform, specifying updated spatial # position of slice in physical space. The transform is obtained via diff --git a/niftymic/base/stack.py b/niftymic/base/stack.py index 21b208ea..5895c1b5 100644 --- a/niftymic/base/stack.py +++ b/niftymic/base/stack.py @@ -7,11 +7,12 @@ # -import SimpleITK as sitk -import numpy as np -# Import libraries import os import re +import numpy as np +import SimpleITK as sitk + +import simplereg.resampler import niftymic.base.slice as sl import pysitk.python_helper as ph @@ -29,6 +30,9 @@ class Stack: def __init__(self): self._is_unity_mask = True + self._deleted_slices = [] + self._history_affine_transforms = [] + self._history_motion_corrections = [] ## # Create Stack instance from file and add corresponding mask. Mask is @@ -64,7 +68,7 @@ def from_filename(cls, stack._filename = filename # Append stacks as SimpleITK and ITK Image objects - stack.sitk = sitk.ReadImage(file_path, sitk.sitkFloat64) + stack.sitk = sitkh.read_nifti_image_sitk(file_path, sitk.sitkFloat64) stack.itk = sitkh.get_itk_from_sitk_image(stack.sitk) # Append masks (either provided or binary mask) @@ -77,12 +81,33 @@ def from_filename(cls, else: if not ph.file_exists(file_path_mask): raise exceptions.FileNotExistent(file_path_mask) - stack.sitk_mask = sitk.ReadImage(file_path_mask, sitk.sitkUInt8) + stack.sitk_mask = sitkh.read_nifti_image_sitk( + file_path_mask, sitk.sitkUInt8) + try: + # ensure masks occupy same physical space + stack.sitk_mask.CopyInformation(stack.sitk) + except RuntimeError as e: + raise IOError( + "Given image and its mask do not occupy the same space: %s" % + e.message) stack._is_unity_mask = False # Append itk object stack.itk_mask = sitkh.get_itk_from_sitk_image(stack.sitk_mask) + # Store current affine transform of image + stack._affine_transform_sitk = sitkh.get_sitk_affine_transform_from_sitk_image( + stack.sitk) + + # Prepare history of affine transforms, i.e. encoded spatial + # position+orientation of stack, and motion estimates of stack + # obtained in the course of the registration/reconstruction process + stack._history_affine_transforms = [] + stack._history_affine_transforms.append(stack._affine_transform_sitk) + + stack._history_motion_corrections = [] + stack._history_motion_corrections.append(sitk.Euler3DTransform()) + # Extract all slices and their masks from the stack and store them if extract_slices: dimenson = stack.sitk.GetDimension() @@ -141,7 +166,7 @@ def from_slice_filenames(cls, stack._filename = prefix_stack # Get 3D images - stack.sitk = sitk.ReadImage( + stack.sitk = sitkh.read_nifti_image_sitk( dir_input + prefix_stack + ".nii.gz", sitk.sitkFloat64) stack.itk = sitkh.get_itk_from_sitk_image(stack.sitk) @@ -149,7 +174,7 @@ def from_slice_filenames(cls, if suffix_mask is not None and \ os.path.isfile(dir_input + prefix_stack + suffix_mask + ".nii.gz"): - stack.sitk_mask = sitk.ReadImage( + stack.sitk_mask = sitkh.read_nifti_image_sitk( dir_input + prefix_stack + suffix_mask + ".nii.gz", sitk.sitkUInt8) stack.itk_mask = sitkh.get_itk_from_sitk_image(stack.sitk_mask) @@ -226,6 +251,17 @@ def from_slices(cls, slices, stack_sitk=None, mask_sitk=None): stack.sitk = stack_sitk stack.itk = sitkh.get_itk_from_sitk_image(stack.sitk) + # Store current affine transform of image + stack._affine_transform_sitk = sitkh.get_sitk_affine_transform_from_sitk_image( + stack.sitk) + + stack._history_affine_transforms = [] + stack._history_affine_transforms.append( + stack._affine_transform_sitk) + + stack._history_motion_corrections = [] + stack._history_motion_corrections.append(sitk.Euler3DTransform()) + stack._N_slices = len(slices) stack._slices = slices @@ -254,7 +290,9 @@ def from_sitk_image(cls, image_sitk, filename="unknown", image_sitk_mask=None, - extract_slices=True): + extract_slices=True, + slice_numbers=None, + ): stack = cls() stack.sitk = sitk.Image(image_sitk) @@ -266,6 +304,13 @@ def from_sitk_image(cls, # Append masks (if provided) if image_sitk_mask is not None: stack.sitk_mask = image_sitk_mask + try: + # ensure mask occupies the same physical space + stack.sitk_mask.CopyInformation(stack.sitk) + except RuntimeError as e: + raise IOError( + "Given image and its mask do not occupy the same space: %s" % + e.message) stack.itk_mask = sitkh.get_itk_from_sitk_image(stack.sitk_mask) if sitk.GetArrayFromImage(stack.sitk_mask).prod() == 1: stack._is_unity_mask = True @@ -279,17 +324,27 @@ def from_sitk_image(cls, # Extract all slices and their masks from the stack and store them if extract_slices: stack._N_slices = stack.sitk.GetSize()[-1] - stack._slices = stack._extract_slices() + stack._slices = stack._extract_slices(slice_numbers=slice_numbers) else: stack._N_slices = 0 stack._slices = None + # Store current affine transform of image + stack._affine_transform_sitk = sitkh.get_sitk_affine_transform_from_sitk_image( + stack.sitk) + + stack._history_affine_transforms = [] + stack._history_affine_transforms.append(stack._affine_transform_sitk) + + stack._history_motion_corrections = [] + stack._history_motion_corrections.append(sitk.Euler3DTransform()) + return stack # Copy constructor # \param[in] stack_to_copy Stack object to be copied # \return copied Stack object - # TODO: That's not really well done! + # TODO: That's not really well done @classmethod def from_stack(cls, stack_to_copy, filename=None): stack = cls() @@ -311,16 +366,21 @@ def from_stack(cls, stack_to_copy, filename=None): else: stack._filename = filename stack._dir = stack_to_copy.get_directory() + stack._deleted_slices = stack_to_copy.get_deleted_slice_numbers() + + # Store current affine transform of image + stack.set_registration_history( + stack_to_copy.get_registration_history()) # Extract all slices and their masks from the stack and store them if # given if stack_to_copy.get_slices() is not None: stack._N_slices = stack_to_copy.get_number_of_slices() - stack._slices = [None]*stack._N_slices + stack._slices = [None] * stack._N_slices slices_to_copy = stack_to_copy.get_slices() - for i in range(0, stack._N_slices): - stack._slices[i] = sl.Slice.from_slice(slices_to_copy[i]) + for j, slice_j in enumerate(slices_to_copy): + stack._slices[j] = sl.Slice.from_slice(slice_j) else: stack._N_slices = 0 stack._slices = None @@ -342,36 +402,85 @@ def from_stack(cls, stack_to_copy, filename=None): # Get all slices of current stack # \return Array of sitk.Images containing slices in 3D space def get_slices(self): - return self._slices + return [s for s in self._slices if s is not None] + ## # Get one particular slice of current stack - # \return requested 3D slice of stack as Slice object + # \date 2018-04-18 22:06:38-0600 + # + # \param self The object + # \param index slice index as integer + # + # \return requested 3D slice of stack as Slice object + # def get_slice(self, index): index = int(index) if abs(index) > self._N_slices - 1: - raise ValueError("Enter a valid index between -%s and %s. Tried: %s" % - (self._N_slices-1, self._N_slices-1, index)) + raise ValueError( + "Enter a valid index between -%s and %s. Tried: %s" % + (self._N_slices - 1, self._N_slices - 1, index)) return self._slices[index] + ## + # Gets the deleted slice numbers, i.e. misregistered slice numbers detected + # by robust outlier algorithm. Indices refer to slice numbers within + # original stack + # \date 2018-07-08 23:06:24-0600 + # + # \param self The object + # + # \return The deleted slice numbers as list of integers. + # + def get_deleted_slice_numbers(self): + return list(self._deleted_slices) + + ## + # Sets the slice. + # \date 2018-04-18 22:05:28-0600 + # + # \param self The object + # \param slice slice as Slice object + # \param index slice index as integer + # + def set_slice(self, slice, index): + if not isinstance(slice, sl.Slice): + raise IOError("Input must be of type Slice") + + index = int(index) + if abs(index) > self._N_slices - 1: + raise ValueError( + "Enter a valid index between -%s and %s. Tried: %s" % + (self._N_slices - 1, self._N_slices - 1, index)) + + self._slices[index] = slice + ## # Delete slice at given index # \date 2017-12-01 00:38:56+0000 # + # Note that index refers to list index of slices (0 ... N_slices) whereas + # "deleted slice index" refers to actual slice number within original stack + # # \param self The object # \param index The index # def delete_slice(self, index): - if self._N_slices > 0: - # check whether slice exists - self.get_slice(index) - - # delete slice at given index - del self._slices[index] - self._N_slices -= 1 + # delete slice at given index + if index in range(self._N_slices): + slice_number = self._slices[index].get_slice_number() + self._deleted_slices.append(slice_number) + self._deleted_slices = sorted(list(set(self._deleted_slices))) + self._slices[index] = None else: - raise RuntimeError("No slice available anymore") + raise ValueError( + "Slice number must be between 0 and %d" % self._N_slices) + + + def get_deleted_slice_numbers(self): + return list(self._deleted_slices) + # Get name of directory where nifti was read from # \return string of directory wher nifti was read from @@ -388,10 +497,26 @@ def set_filename(self, filename): def get_filename(self): return self._filename + # Get history history of affine transforms, i.e. encoded spatial + # position+orientation of slice, and rigid motion estimates of slice + # obtained in the course of the registration/reconstruction process + # \return list of sitk.AffineTransform and sitk.Euler3DTransform objects + def get_registration_history(self): + affine_transforms = list(self._history_affine_transforms) + motion_corrections = list(self._history_motion_corrections) + return affine_transforms, motion_corrections + + def set_registration_history(self, registration_history): + affine_transform_sitk = registration_history[0][-1] + self._update_affine_transform(affine_transform_sitk) + + self._history_affine_transforms = [a for a in registration_history[0]] + self._history_motion_corrections = [t for t in registration_history[1]] + # Get number of slices of stack # \return number of slices of stack def get_number_of_slices(self): - return self._N_slices + return len(self.get_slices()) def is_unity_mask(self): return self._is_unity_mask @@ -425,7 +550,15 @@ def show_slices(self): # \param[in] directory string specifying where the output will be written to (default="/tmp/") # \param[in] filename string specifying the filename. If not given the assigned one within Stack will be chosen. # \param[in] write_slices boolean indicating whether each Slice of the stack shall be written (default=False) - def write(self, directory, filename=None, write_mask=False, write_slices=False, write_transforms=False, suffix_mask="_mask"): + def write(self, + directory, + filename=None, + write_stack=True, + write_mask=False, + write_slices=False, + write_transforms=False, + suffix_mask="_mask", + use_float32=True): # Create directory if not existing ph.create_directory(directory) @@ -437,10 +570,16 @@ def write(self, directory, filename=None, write_mask=False, write_slices=False, full_file_name = os.path.join(directory, filename) # Write file to specified location - ph.print_info("Write image stack to %s.nii.gz ... " % - (full_file_name), newline=False) - sitkh.write_nifti_image_sitk(self.sitk, full_file_name + ".nii.gz") - print("done") + if write_stack: + ph.print_info("Write image stack to %s.nii.gz ... " % + (full_file_name), newline=False) + if use_float32: + image_sitk = sitk.Cast(self.sitk, sitk.sitkFloat32) + else: + image_sitk = self.sitk + sitkh.write_nifti_image_sitk( + image_sitk, full_file_name + ".nii.gz") + print("done") # Write mask to specified location if given if self.sitk_mask is not None: @@ -455,11 +594,15 @@ def write(self, directory, filename=None, write_mask=False, write_slices=False, suffix_mask)) print("done") - # print("Stack was successfully written to %s.nii.gz" - # %(full_file_name)) + if write_transforms: + stack_transform_sitk = self._history_motion_corrections[-1] + sitk.WriteTransform( + stack_transform_sitk, + os.path.join(directory, self.get_filename() + ".tfm") + ) # Write each separate Slice of stack (if they exist) - if write_slices: + if write_slices or write_transforms: try: # Check whether variable exists # if 'self._slices' not in locals() or all(i is None for i in @@ -472,9 +615,13 @@ def write(self, directory, filename=None, write_mask=False, write_slices=False, # Write slices else: - if write_transforms: + if write_transforms and write_slices: + ph.print_info( + "Write image slices and slice transforms to %s ... " % + directory, newline=False) + elif write_transforms and not write_slices: ph.print_info( - "Write image slices + transforms to %s ... " % + "Write slice transforms to %s ... " % directory, newline=False) else: ph.print_info( @@ -485,6 +632,7 @@ def write(self, directory, filename=None, write_mask=False, write_slices=False, directory=directory, filename=filename, write_transform=write_transforms, + write_slice=write_slices, suffix_mask=suffix_mask) print("done") @@ -500,8 +648,18 @@ def write(self, directory, filename=None, write_mask=False, write_slices=False, # def update_motion_correction(self, affine_transform_sitk): - # Update sitk and itk stack position - self._update_sitk_and_itk_stack_position(affine_transform_sitk) + # Update rigid motion estimate + current_rigid_motion_estimate = sitkh.get_composite_sitk_affine_transform( + affine_transform_sitk, self._history_motion_corrections[-1]) + self._history_motion_corrections.append(current_rigid_motion_estimate) + + # New affine transform of slice after rigid motion correction + affine_transform = sitkh.get_composite_sitk_affine_transform( + affine_transform_sitk, self._affine_transform_sitk) + + # Update affine transform of stack, i.e. change image origin and + # direction in physical space + self._update_affine_transform(affine_transform) # Update slices if self.get_slices() is not None: @@ -527,15 +685,25 @@ def update_motion_correction_of_slices(self, affine_transforms_sitk): raise ValueError("Number of affine transforms does not match the " "number of slices") - def _update_sitk_and_itk_stack_position(self, affine_transform_sitk): + def _update_affine_transform(self, affine_transform_sitk): - # Apply transform to 3D image / stack of slices - self.sitk = sitkh.get_transformed_sitk_image( - self.sitk, affine_transform_sitk) + # Ensure correct object type + self._affine_transform_sitk = sitk.AffineTransform( + affine_transform_sitk) - # Update header information of other associated images - origin = self.sitk.GetOrigin() - direction = self.sitk.GetDirection() + # Append transform to registration history + self._history_affine_transforms.append(affine_transform_sitk) + + # Get origin and direction of transformed 3D slice given the new + # spatial transform + origin = sitkh.get_sitk_image_origin_from_sitk_affine_transform( + affine_transform_sitk, self.sitk) + direction = sitkh.get_sitk_image_direction_from_sitk_affine_transform( + affine_transform_sitk, self.sitk) + + # Update image objects + self.sitk.SetOrigin(origin) + self.sitk.SetDirection(direction) self.sitk_mask.SetOrigin(origin) self.sitk_mask.SetDirection(direction) @@ -708,35 +876,27 @@ def get_resampled_stack(self, resampling_grid=None, spacing=None, interpolator=" self.sitk_mask, resampling_grid, sitk.Euler3DTransform(), - interpolator, + sitk.sitkNearestNeighbor, 0, self.sitk_mask.GetPixelIDValue()) else: - spacing0 = np.array(self.sitk.GetSpacing()) - size0 = np.array(self.sitk.GetSize()) - size = np.round(size0 * spacing0 / spacing).astype("int") - - resampled_stack_sitk = sitk.Resample( - self.sitk, - size, - sitk.Euler3DTransform(), - interpolator, - self.sitk.GetOrigin(), - spacing, - self.sitk.GetDirection(), - default_pixel_value, - self.sitk.GetPixelIDValue()) - - resampled_stack_sitk_mask = sitk.Resample( - self.sitk_mask, - size, - sitk.Euler3DTransform(), - interpolator, - self.sitk.GetOrigin(), - spacing, - self.sitk.GetDirection(), - 0, - self.sitk_mask.GetPixelIDValue()) + resampler = simplereg.resampler.Resampler + resampled_stack_sitk = resampler.get_resampled_image_sitk( + image_sitk=self.sitk, + spacing=spacing, + interpolator=interpolator, + padding=default_pixel_value, + add_to_grid=extra_frame, + add_to_grid_unit="mm", + ) + resampled_stack_sitk_mask = resampler.get_resampled_image_sitk( + image_sitk=self.sitk_mask, + spacing=spacing, + interpolator=sitk.sitkNearestNeighbor, + padding=0, + add_to_grid=extra_frame, + add_to_grid_unit="mm", + ) # Create Stack instance if filename is None: @@ -793,11 +953,15 @@ def get_isotropically_resampled_stack_from_slices(self, resolution=None, interpo size_new = size spacing_new = spacing # Update information according to isotropic resolution - size_new[2] = np.round(spacing[2]/spacing[0]*size[2]).astype("int") + size_new[2] = np.round( + spacing[2] / spacing[0] * size[2]).astype("int") spacing_new[2] = spacing[0] else: - spacing_new = np.ones(3)*resolution - size_new = np.round(spacing/spacing_new*size).astype("int") + spacing_new = np.ones(3) * resolution + size_new = np.round(spacing / spacing_new * size).astype("int") + + # For Python3: sitk.Resample in Python3 does not like np.int types! + size_new = [int(i) for i in size_new] # Resample image and its mask to isotropic grid isotropic_resampled_stack_sitk = sitk.Resample( @@ -856,69 +1020,29 @@ def get_isotropically_resampled_stack(self, resolution=None, interpolator="Linea except: raise ValueError("Error: interpolator is not known") - # Read original spacing (voxel dimension) and size of target stack: - spacing = np.array(self.sitk.GetSpacing()) - size = np.array(self.sitk.GetSize()).astype("int") - origin = np.array(self.sitk.GetOrigin()) - direction = self.sitk.GetDirection() - if resolution is None: - size_new = size - spacing_new = spacing - # Update information according to isotropic resolution - size_new[2] = np.round(spacing[2]/spacing[0]*size[2]).astype("int") - spacing_new[2] = spacing[0] + spacing = self.sitk.GetSpacing()[0] else: - spacing_new = np.ones(3)*resolution - size_new = np.round(spacing/spacing_new*size).astype("int") - - if extra_frame is not 0: - - # Get extra_frame in voxel space - extra_frame_vox = np.round( - extra_frame/spacing_new[0]).astype("int") - - # Compute size of resampled stack by considering additional - # extra_frame - size_new = size_new + 2*extra_frame_vox - - # Compute origin of resampled stack by considering additional - # extra_frame - a_x = self.sitk.TransformIndexToPhysicalPoint((1, 0, 0)) - origin - a_y = self.sitk.TransformIndexToPhysicalPoint((0, 1, 0)) - origin - a_z = self.sitk.TransformIndexToPhysicalPoint((0, 0, 1)) - origin - e_x = a_x/np.linalg.norm(a_x) - e_y = a_y/np.linalg.norm(a_y) - e_z = a_z/np.linalg.norm(a_z) - - translation = (e_x + e_y + e_z)*extra_frame_vox*spacing_new - - origin = origin - translation + spacing = resolution # Resample image and its mask to isotropic grid - default_pixel_value = 0.0 - - isotropic_resampled_stack_sitk = sitk.Resample( - self.sitk, - size_new, - sitk.Euler3DTransform(), - interpolator, - origin, - spacing_new, - direction, - default_pixel_value, - self.sitk.GetPixelIDValue()) - - isotropic_resampled_stack_sitk_mask = sitk.Resample( - self.sitk_mask, - size_new, - sitk.Euler3DTransform(), - sitk.sitkNearestNeighbor, - origin, - spacing_new, - direction, - 0, - self.sitk_mask.GetPixelIDValue()) + resampler = simplereg.resampler.Resampler + isotropic_resampled_stack_sitk = resampler.get_resampled_image_sitk( + image_sitk=self.sitk, + spacing=spacing, + interpolator=interpolator, + padding=0.0, + add_to_grid=extra_frame, + add_to_grid_unit="mm", + ) + isotropic_resampled_stack_sitk_mask = resampler.get_resampled_image_sitk( + image_sitk=self.sitk_mask, + spacing=spacing, + interpolator=sitk.sitkNearestNeighbor, + padding=0, + add_to_grid=extra_frame, + add_to_grid_unit="mm", + ) if mask_dilation_radius > 0: dilater = sitk.BinaryDilateImageFilter() @@ -981,7 +1105,7 @@ def get_increased_stack(self, extra_slices_z=0): # Create Stack instance stack = self.from_sitk_image( - isotropic_resampled_stack_sitk, "zincreased_"+self._filename, isotropic_resampled_stack_sitk_mask) + isotropic_resampled_stack_sitk, "zincreased_" + self._filename, isotropic_resampled_stack_sitk_mask) return stack @@ -994,22 +1118,32 @@ def get_cropped_stack_based_on_mask(self, boundary_i=0, boundary_j=0, boundary_k if x_range is None: return None - # Crop to image region defined by rectangular mask - stack_crop_sitk = self._crop_image_to_region( - self.sitk, x_range, y_range, z_range) + if unit == "mm": + spacing = self.sitk.GetSpacing() + boundary_i = np.round(boundary_i / float(spacing[0])) + boundary_j = np.round(boundary_j / float(spacing[1])) + boundary_k = np.round(boundary_k / float(spacing[2])) - # Increase image region - stack_crop_sitk = sitkh.get_altered_field_of_view_sitk_image( - stack_crop_sitk, boundary_i, boundary_j, boundary_k, unit=unit) + shape = self.sitk.GetSize() + x_range[0] = np.max([0, x_range[0] - boundary_i]) + x_range[1] = np.min([shape[0], x_range[1] + boundary_i]) - # Resample original image and mask to specified image region - image_crop_sitk = sitk.Resample(self.sitk, stack_crop_sitk, sitk.Euler3DTransform( - ), sitk.sitkNearestNeighbor, 0, self.sitk.GetPixelIDValue()) - mask_crop_sitk = sitk.Resample(self.sitk_mask, stack_crop_sitk, sitk.Euler3DTransform( - ), sitk.sitkNearestNeighbor, 0, self.sitk_mask.GetPixelIDValue()) + y_range[0] = np.max([0, y_range[0] - boundary_j]) + y_range[1] = np.min([shape[1], y_range[1] + boundary_j]) + z_range[0] = np.max([0, z_range[0] - boundary_k]) + z_range[1] = np.min([shape[2], z_range[1] + boundary_k]) + + # Crop to image region defined by rectangular mask + image_crop_sitk = self._crop_image_to_region( + self.sitk, x_range, y_range, z_range) + mask_crop_sitk = self._crop_image_to_region( + self.sitk_mask, x_range, y_range, z_range) + + slice_numbers = range(z_range[0], z_range[1]) stack = self.from_sitk_image( - image_crop_sitk, self._filename, mask_crop_sitk) + image_crop_sitk, self._filename, mask_crop_sitk, + slice_numbers=slice_numbers) return stack @@ -1045,17 +1179,17 @@ def _get_rectangular_masked_region(self, mask_sitk): # Non-zero elements of numpy array nda defining x_range ran = np.nonzero(sum_yz)[0] range_x[0] = np.max([0, ran[0]]) - range_x[1] = np.min([shape[0], ran[-1]+1]) + range_x[1] = np.min([shape[0], ran[-1] + 1]) # Non-zero elements of numpy array nda defining y_range ran = np.nonzero(sum_xz)[0] range_y[0] = np.max([0, ran[0]]) - range_y[1] = np.min([shape[1], ran[-1]+1]) + range_y[1] = np.min([shape[1], ran[-1] + 1]) # Non-zero elements of numpy array nda defining z_range ran = np.nonzero(sum_xy)[0] range_z[0] = np.max([0, ran[0]]) - range_z[1] = np.min([shape[2], ran[-1]+1]) + range_z[1] = np.min([shape[2], ran[-1] + 1]) # Numpy reads the array as z,y,x coordinates! So swap them accordingly return range_z.astype(int), range_y.astype(int), range_x.astype(int) @@ -1078,17 +1212,20 @@ def _crop_image_to_region(self, image_sitk, range_x, range_y, range_z): # Burst the stack into its slices and return all slices of the stack # return list of Slice objects - def _extract_slices(self): + def _extract_slices(self, slice_numbers=None): + + slices = [None] * self._N_slices - slices = [None]*self._N_slices + if slice_numbers is None: + slice_numbers = range(0, self._N_slices) # Extract slices and add masks for i in range(0, self._N_slices): slices[i] = sl.Slice.from_sitk_image( - slice_sitk=self.sitk[:, :, i:i+1], + slice_sitk=self.sitk[:, :, i:i + 1], filename=self._filename, - slice_number=i, - slice_sitk_mask=self.sitk_mask[:, :, i:i+1]) + slice_number=slice_numbers[i], + slice_sitk_mask=self.sitk_mask[:, :, i:i + 1]) return slices diff --git a/niftymic/definitions.py b/niftymic/definitions.py index 41773a5f..248a48fc 100644 --- a/niftymic/definitions.py +++ b/niftymic/definitions.py @@ -17,5 +17,7 @@ REGEX_FILENAME_EXTENSIONS = "(" + "|".join(ALLOWED_EXTENSIONS) + ")" ALLOWED_INTERPOLATORS = ["NearestNeighbor", "Linear", "BSpline"] +TEMPLATES_INFO = "templates_info.json" + # Set default viewer VIEWER = ITKSNAP_EXE diff --git a/niftymic/reconstruction/admm_solver.py b/niftymic/reconstruction/admm_solver.py index a8e1de85..d51d56dc 100644 --- a/niftymic/reconstruction/admm_solver.py +++ b/niftymic/reconstruction/admm_solver.py @@ -65,6 +65,7 @@ def __init__(self, predefined_covariance=None, rho=0.5, iterations=10, + use_masks=1, verbose=1, ): @@ -82,6 +83,7 @@ def __init__(self, huber_gamma=huber_gamma, deconvolution_mode=deconvolution_mode, predefined_covariance=predefined_covariance, + use_masks=use_masks, verbose=verbose, ) diff --git a/niftymic/reconstruction/primal_dual_solver.py b/niftymic/reconstruction/primal_dual_solver.py index f3a01ff7..2c7b0b89 100644 --- a/niftymic/reconstruction/primal_dual_solver.py +++ b/niftymic/reconstruction/primal_dual_solver.py @@ -50,6 +50,7 @@ def __init__(self, reg_huber_gamma=0.05, iterations=10, alg_type="ALG2", + use_masks=1, verbose=0, ): @@ -66,6 +67,7 @@ def __init__(self, huber_gamma=huber_gamma, deconvolution_mode=deconvolution_mode, predefined_covariance=predefined_covariance, + use_masks=use_masks, verbose=verbose, ) @@ -142,6 +144,7 @@ def get_solver(self): x_scale=x_scale, data_loss=self._data_loss, data_loss_scale=self._data_loss_scale, + minimizer=self._minimizer, verbose=self._verbose) if self._reg_type == "TV": diff --git a/niftymic/reconstruction/solver.py b/niftymic/reconstruction/solver.py index ed74053f..caef756a 100644 --- a/niftymic/reconstruction/solver.py +++ b/niftymic/reconstruction/solver.py @@ -360,28 +360,24 @@ def _get_M_y(self): # Define index for first voxel of first slice within array i_min = 0 - for i in range(0, self._N_stacks): - stack = self._stacks[i] + for i, stack in enumerate(self._stacks): slices = stack.get_slices() # Get number of voxels of each slice in current stack N_slice_voxels = np.array(slices[0].sitk.GetSize()).prod() - for j in range(0, stack.get_number_of_slices()): + for j, slice_j in enumerate(slices): # Define index for last voxel to specify current slice # (exlusive) i_max = i_min + N_slice_voxels - # Get current slice - slice_k = slices[j] - # Apply M_k y_k if self._use_masks: slice_itk = self._linear_operators.M_itk( - slice_k.itk, slice_k.itk_mask) + slice_j.itk, slice_j.itk_mask) else: - slice_itk = slice_k.itk + slice_itk = slice_j.itk slice_nda_vec = self._itk2np.GetArrayFromImage( slice_itk).flatten() @@ -472,21 +468,20 @@ def _MA(self, reconstruction_nda_vec): # Define index for first voxel of first slice within array i_min = 0 - for i in range(0, self._N_stacks): - stack = self._stacks[i] + for i, stack in enumerate(self._stacks): slices = stack.get_slices() # Get number of voxels of each slice in current stack N_slice_voxels = np.array(slices[0].sitk.GetSize()).prod() - for j in range(0, stack.get_number_of_slices()): + for j, slice_j in enumerate(slices): # Define index for last voxel to specify current slice # (exclusive) i_max = i_min + N_slice_voxels # Compute M_k A_k y_k - slice_itk = self._Mk_Ak(x_itk, slices[j]) + slice_itk = self._Mk_Ak(x_itk, slice_j) slice_nda = self._itk2np.GetArrayFromImage(slice_itk) # Fill corresponding elements @@ -519,29 +514,25 @@ def _A_adj_M(self, stacked_slices_nda_vec): # Define index for first voxel of first slice within array i_min = 0 - for i in range(0, self._N_stacks): - stack = self._stacks[i] + for i, stack in enumerate(self._stacks): slices = stack.get_slices() # Get number of voxels of each slice in current stack N_slice_voxels = np.array(slices[0].sitk.GetSize()).prod() - for j in range(0, stack.get_number_of_slices()): + for j, slice_j in enumerate(slices): # Define index for last voxel to specify current slice # (exlusive) i_max = i_min + N_slice_voxels - # Get current slice - slice_k = slices[j] - # Extract 1D corresponding to current slice and convert it to # itk.Object slice_itk = self._get_itk_image_from_array_vec( - stacked_slices_nda_vec[i_min:i_max], slice_k.itk) + stacked_slices_nda_vec[i_min:i_max], slice_j.itk) # Apply A_k' M_k on current slice - Ak_adj_Mk_slice_itk = self._Ak_adj_Mk(slice_itk, slice_k) + Ak_adj_Mk_slice_itk = self._Ak_adj_Mk(slice_itk, slice_j) Ak_adj_Mk_slice_nda_vec = self._itk2np.GetArrayFromImage( Ak_adj_Mk_slice_itk).flatten() diff --git a/niftymic/reconstruction/tikhonov_solver.py b/niftymic/reconstruction/tikhonov_solver.py index b07e9909..6b322db6 100644 --- a/niftymic/reconstruction/tikhonov_solver.py +++ b/niftymic/reconstruction/tikhonov_solver.py @@ -102,8 +102,8 @@ def __init__(self, data_loss_scale=1, huber_gamma=1.345, predefined_covariance=None, - verbose=1, use_masks=True, + verbose=1, ): # Run constructor of superclass diff --git a/niftymic/registration/flirt.py b/niftymic/registration/flirt.py index 93e6284b..783e6d6d 100644 --- a/niftymic/registration/flirt.py +++ b/niftymic/registration/flirt.py @@ -90,9 +90,6 @@ def _run(self): elif self.get_registration_type() == "Affine": options += " -dof 12" - if self._use_verbose: - options += " -verbose 1" - self._registration_method = simplereg.flirt.FLIRT( fixed_sitk=self._fixed.sitk, moving_sitk=self._moving.sitk, diff --git a/niftymic/registration/niftyreg.py b/niftymic/registration/niftyreg.py index 2ce12b19..c1e90b25 100644 --- a/niftymic/registration/niftyreg.py +++ b/niftymic/registration/niftyreg.py @@ -31,7 +31,7 @@ def __init__(self, use_fixed_mask=False, use_moving_mask=False, use_verbose=False, - options="", + options="-voff", registration_type="Rigid", ): @@ -86,15 +86,13 @@ def _run(self): if self.get_registration_type() == "Rigid": options += " -rigOnly" - if not self._use_verbose: - options += " -voff" - self._registration_method = simplereg.niftyreg.RegAladin( fixed_sitk=self._fixed.sitk, moving_sitk=self._moving.sitk, fixed_sitk_mask=fixed_sitk_mask, moving_sitk_mask=moving_sitk_mask, options=options, + verbose=self._use_verbose, ) self._registration_method.run() @@ -113,7 +111,7 @@ def __init__(self, use_fixed_mask=False, use_moving_mask=False, use_verbose=False, - options="", + options="-voff", ): RegistrationMethod.__init__(self, @@ -160,15 +158,13 @@ def _run(self): options = self._options - if not self._use_verbose: - options += " -voff" - self._registration_method = simplereg.niftyreg.RegF3D( fixed_sitk=self._fixed.sitk, moving_sitk=self._moving.sitk, fixed_sitk_mask=fixed_sitk_mask, moving_sitk_mask=moving_sitk_mask, - options=options + options=options, + verbose=self._use_verbose, ) self._registration_method.run() diff --git a/niftymic/utilities/brain_stripping.py b/niftymic/utilities/brain_stripping.py index 1daee21f..5c7aef4c 100644 --- a/niftymic/utilities/brain_stripping.py +++ b/niftymic/utilities/brain_stripping.py @@ -18,6 +18,7 @@ import pysitk.simple_itk_helper as sitkh import pysitk.python_helper as ph +import niftymic.base.stack as st from niftymic.definitions import DIR_TMP @@ -58,6 +59,8 @@ def __init__(self, self._sitk_brain_mask = None self._sitk_skull_image = None + self._stack = None + ## # Initialize brain stripping class based on image to be read # \date 2016-10-12 12:19:18+0100 @@ -87,8 +90,9 @@ def from_filename(cls, compute_brain_mask=compute_brain_mask, compute_skull_image=compute_skull_image, dir_tmp=dir_tmp) - self._sitk = sitk.ReadImage( - dir_input + filename + ".nii.gz", sitk.sitkFloat64) + self._sitk = sitkh.read_nifti_image_sitk( + os.path.join(dir_input, "%s.nii.gz" % filename), + sitk.sitkFloat64) return self @@ -123,6 +127,39 @@ def from_sitk_image(cls, return self + ## + # Initialize brain stripping class based on given Stack object + # \date 2018-01-18 00:38:53+0000 + # + # \param cls The cls + # \param stack image as Stack object + # \param compute_brain_image Boolean flag for computing brain image + # \param compute_brain_mask Boolean flag for computing brain image + # mask + # \param compute_skull_image Boolean flag for computing skull mask + # \param dir_tmp Directory where temporary results are + # written to, string + # + # \return object + # + @classmethod + def from_stack(cls, + stack, + compute_brain_image=False, + compute_brain_mask=True, + compute_skull_image=False, + dir_tmp=os.path.join(DIR_TMP, "BrainExtractionTool")): + + self = cls(compute_brain_image=compute_brain_image, + compute_brain_mask=compute_brain_mask, + compute_skull_image=compute_skull_image, + dir_tmp=dir_tmp) + + self._stack = stack + self._sitk = sitk.Image(stack.sitk) + + return self + ## # Sets the sitk image for brain stripping # \date 2016-10-12 15:46:20+0100 @@ -188,6 +225,32 @@ def get_input_image_sitk(self): return sitk.Image(self._sitk) + ## + # Gets the brain masked stack. + # \date 2018-01-18 00:44:49+0000 + # + # \param self The object + # \param filename The filename + # \param extract_slices Extract slices of stack; boolean + # + # \return Returns image as Stack object holding obtained brain mask + # + def get_brain_masked_stack(self, filename="Unknown", extract_slices=False): + if self._sitk_brain_mask is None: + raise ValueError("Brain mask was not asked for. " + "Set option '-m' and run again.") + + if self._stack is not None: + filename = self._stack.get_filename() + + stack = st.Stack.from_sitk_image( + image_sitk=self._sitk, + image_sitk_mask=self._sitk_brain_mask, + filename=filename, + extract_slices=extract_slices + ) + return stack + ## # Get computed brain image # \date 2016-10-12 14:33:53+0100 @@ -317,7 +380,7 @@ def get_mask_around_skull(self, # Go slice by slice for i in range(0, shape[0]): - slice_mask_sitk = mask_sitk[:, :, i:i+1] + slice_mask_sitk = mask_sitk[:, :, i:i + 1] # Dilate mask of slice if dilate_radius > 0: @@ -397,12 +460,13 @@ def _run_bet_for_brain_stripping(self, debug=0): bet.run() if self._compute_brain_image: - self._sitk_brain_image = sitk.ReadImage( + self._sitk_brain_image = sitkh.read_nifti_image_sitk( path_to_res, sitk.sitkFloat64) if self._compute_brain_mask: - self._sitk_brain_mask = sitk.ReadImage( + self._sitk_brain_mask = sitkh.read_nifti_image_sitk( path_to_res_mask, sitk.sitkUInt8) if self._compute_skull_image: - self._sitk_skull_image = sitk.ReadImage(path_to_res_skull) + self._sitk_skull_image = sitkh.read_nifti_image_sitk( + path_to_res_skull) diff --git a/niftymic/utilities/data_preprocessing.py b/niftymic/utilities/data_preprocessing.py index dd02ca61..81839909 100644 --- a/niftymic/utilities/data_preprocessing.py +++ b/niftymic/utilities/data_preprocessing.py @@ -6,8 +6,6 @@ # \date May 2017 # -# Import libraries - import niftymic.base.stack as st import niftymic.utilities.intensity_correction as ic import niftymic.utilities.n4_bias_field_correction as n4bfc diff --git a/niftymic/utilities/input_arparser.py b/niftymic/utilities/input_arparser.py index c4959158..31502501 100644 --- a/niftymic/utilities/input_arparser.py +++ b/niftymic/utilities/input_arparser.py @@ -9,8 +9,8 @@ import os import re +import six import sys -import json import argparse import pysitk.python_helper as ph @@ -78,7 +78,20 @@ def print_arguments(self, args, title="Input Parameters:"): ph.print_title(title) for arg in sorted(vars(args)): ph.print_info("%s: " % (arg), newline=False) - print(getattr(args, arg)) + vals = getattr(args, arg) + + if type(vals) is list: + # print list element in new lines, unless only one entry in list + # if len(vals) == 1: + # print(vals[0]) + # else: + print("") + for val in vals: + print("\t%s" % val) + else: + print(vals) + ph.print_line_separator(add_newline=False) + print("") ## # Writes a performed script execution. @@ -89,15 +102,18 @@ def print_arguments(self, args, title="Input Parameters:"): # os.path.abspath(__file__) # \param prefix filename prefix # - def write_performed_script_execution(self, - file, - prefix="config"): + def log_config(self, + file, + prefix="config"): # parser returns options with underscores, e.g. 'dir_output' dic_with_underscores = vars(self._parser.parse_args()) # get output directory to write log/config file - dir_output = dic_with_underscores["dir_output"] + try: + dir_output = dic_with_underscores["dir_output"] + except KeyError: + dir_output = os.path.dirname(dic_with_underscores["output"]) # build output file name name = os.path.basename(file).split(".")[0] @@ -113,14 +129,10 @@ def write_performed_script_execution(self, # e.g. "dir-output" instead of "dir_output" dic = { re.sub("_", "-", k): v - for k, v in dic_with_underscores.iteritems()} + for k, v in six.iteritems(dic_with_underscores)} # write config file to output - ph.create_directory(dir_output) - with open(path_to_config_file, "w") as fp: - json.dump(dic, fp, sort_keys=True, indent=4) - ph.print_info( - "Configuration written to '%s'." % path_to_config_file) + ph.write_dictionary_to_json(dic, path_to_config_file, verbose=True) def add_filename( self, @@ -152,6 +164,17 @@ def add_dir_input( ): self._add_argument(dict(locals())) + def add_dir_input_mc( + self, + option_string="--dir-input-mc", + type=str, + help="Input directory where transformation files (.tfm) for " + "motion-corrected slices are stored.", + default=None, + required=False, + ): + self._add_argument(dict(locals())) + def add_subfolder_motion_correction( self, option_string="--subfolder-motion-correction", @@ -189,7 +212,17 @@ def add_filenames( self, option_string="--filenames", nargs="+", - help="Filenames.", + help="Paths to NIfTI file images %s." % (IMAGE_TYPES), + default=None, + required=False, + ): + self._add_argument(dict(locals())) + + def add_filenames_masks( + self, + option_string="--filenames-masks", + nargs="+", + help="Paths to NIfTI file image masks %s." % (IMAGE_TYPES), default=None, required=False, ): @@ -256,6 +289,16 @@ def add_metric( ): self._add_argument(dict(locals())) + def add_metric_radius( + self, + option_string="--metric-radius", + type=int, + help="Radius in case metric 'ANTSNeighborhoodCorrelation' is chosen.", + default=10, + required=False, + ): + self._add_argument(dict(locals())) + def add_labels( self, option_string="--labels", @@ -297,6 +340,16 @@ def add_reference_mask( ): self._add_argument(dict(locals())) + def add_output( + self, + option_string="--output", + type=str, + help="Path to output image file %s." % (IMAGE_TYPES), + default=None, + required=False, + ): + self._add_argument(dict(locals())) + def add_dir_output( self, option_string="--dir-output", @@ -332,6 +385,27 @@ def add_use_masks_srr( ): self._add_argument(dict(locals())) + def add_outlier_rejection( + self, + option_string="--outlier-rejection", + type=int, + help="Turn on/off use of outlier rejection mechanism to eliminate " + "misregistered slices.", + default=0, + required=False, + ): + self._add_argument(dict(locals())) + + def add_use_robust_registration( + self, + option_string="--use-robust-registration", + type=int, + help="Turn on/off use of robust slice-to-volume registration.", + default=0, + required=False, + ): + self._add_argument(dict(locals())) + def add_boundary_stacks( self, option_string="--boundary-stacks", @@ -513,6 +587,44 @@ def add_alpha_first( ): self._add_argument(dict(locals())) + def add_threshold( + self, + option_string="--threshold", + type=float, + help="Threshold between 0 and 1 to detect misregistered slices based " + "on NCC in final cycle.", + default=0.7, + ): + self._add_argument(dict(locals())) + + def add_threshold_first( + self, + option_string="--threshold-first", + type=float, + help="Threshold between 0 and 1 to detect misregistered slices based " + "on NCC in first cycle.", + default=0.6, + ): + self._add_argument(dict(locals())) + + def add_s2v_smoothing( + self, + option_string="--s2v-smoothing", + type=float, + help="Value for Gaussian process parameter smoothing.", + default=0.5, + ): + self._add_argument(dict(locals())) + + def add_interleave( + self, + option_string="--interleave", + type=int, + help="Interleave used for slice acquisition", + default=2, + ): + self._add_argument(dict(locals())) + def add_iter_max( self, option_string="--iter-max", @@ -586,7 +698,7 @@ def add_data_loss_scale( def add_pd_alg_type( self, - option_string="-pd_alg_type", + option_string="--pd-alg-type", type=str, help="Algorithm used to dynamically update parameters for each " "iteration of the dual algorithm. " @@ -648,11 +760,11 @@ def add_isotropic_resolution( ): self._add_argument(dict(locals())) - def add_log_script_execution( + def add_log_config( self, - option_string="--log-script-execution", + option_string="--log-config", type=int, - help="Turn on/off log for script execution.", + help="Turn on/off configuration of executed script.", default=0, ): self._add_argument(dict(locals())) @@ -706,6 +818,13 @@ def add_option( ): self._add_argument(dict(locals())) + def add_argument( + self, + *a, + **k + ): + self._parser.add_argument(*a, **k) + def add_psf_aware( self, option_string='--psf-aware', @@ -823,10 +942,10 @@ def _parse_config_file(self): # Read config file and insert all config entries into sys.argv (read by # argparse later) with open(path_to_config_file) as json_file: - dic = json.load(json_file) + dic = ph.read_dictionary_from_json(json_file) # Insert all config entries into sys.argv - for k, v in dic.iteritems(): + for k, v in six.iteritems(dic): # A 'None' entry should be ignored if v is None: @@ -864,13 +983,17 @@ def _add_argument(self, allvars): # Build dictionary for additional, optional parameters kwargs = {} - for key, value in allvars.iteritems(): + for key, value in six.iteritems(allvars): kwargs[key] = value # Add information on default value in case provided if 'default' in kwargs.keys(): - txt_default = " [default: %s]" % (str(kwargs['default'])) + if type(kwargs['default']) == list: + txt = " ".join([str(i) for i in kwargs['default']]) + else: + txt = str(kwargs['default']) + txt_default = " [default: %s]" % txt # Case where 'required' key is given: if 'required' in kwargs.keys(): diff --git a/niftymic/utilities/intensity_correction.py b/niftymic/utilities/intensity_correction.py index 39156ca3..5691e752 100644 --- a/niftymic/utilities/intensity_correction.py +++ b/niftymic/utilities/intensity_correction.py @@ -51,6 +51,7 @@ class IntensityCorrection(object): def __init__(self, stack=None, reference=None, + use_stack_mask=True, use_reference_mask=True, use_individual_slice_correction=False, use_verbose=False, @@ -87,6 +88,7 @@ def __init__(self, self._use_verbose = use_verbose self._use_reference_mask = use_reference_mask + self._use_stack_mask = use_stack_mask self._use_individual_slice_correction = use_individual_slice_correction self._prefix_corrected = prefix_corrected @@ -135,6 +137,9 @@ def use_verbose(self, verbose): def use_reference_mask(self, use_reference_mask): self._use_reference_mask = use_reference_mask + def use_stack_mask(self, use_stack_mask): + self._use_stack_mask = use_stack_mask + ## # Sets the use individual slice correction. # \date 2016-11-22 22:47:47+0000 @@ -262,7 +267,7 @@ def _run_intensity_correction(self, correction_model): for i in range(0, N_slices): if self._use_verbose: sys.stdout.write("Slice %2d/%d: " % - (i, self._stack.get_number_of_slices()-1)) + (i, self._stack.get_number_of_slices() - 1)) sys.stdout.flush() if self._additional_stack is None: nda[i, :, :], correction_coefficients[i, :] = self._apply_intensity_correction[ @@ -270,7 +275,6 @@ def _run_intensity_correction(self, correction_model): else: nda[i, :, :], correction_coefficients[i, :], nda_additional_stack[i, :, :] = self._apply_intensity_correction[ correction_model](nda[i, :, :], nda_reference[i, :, :], nda_mask[i, :, :], nda_additional_stack[i, :, :]) - correction_coefficients[:, ] = np.tile(cc, (N_slices, 1)) else: if self._use_verbose: ph.print_info("Run " + correction_model + @@ -322,9 +326,9 @@ def _apply_affine_intensity_correction(self, nda, nda_reference, nda_mask, nda_a ph.print_info("(c1, c0) = (%.3f, %.3f)" % (c1, c0)) if nda_additional_stack is None: - return nda*c1 + c0, np.array([c1, c0]) + return nda * c1 + c0, np.array([c1, c0]) else: - return nda*c1 + c0, np.array([c1, c0]), nda_additional_stack*c1 + c0 + return nda * c1 + c0, np.array([c1, c0]), nda_additional_stack * c1 + c0 ## # Perform linear intensity correction via normal equations @@ -349,15 +353,15 @@ def _apply_linear_intensity_correction(self, nda, nda_reference, nda_mask, nda_a # ph.show_2D_array_list([nda, nda_reference]) # Solve via normal equations: c1 = x'y/(x'x) - c1 = x.dot(y)/x.dot(x) + c1 = x.dot(y) / x.dot(x) if self._use_verbose: ph.print_info("c1 = %.3f" % (c1)) if nda_additional_stack is None: - return nda*c1, c1 + return nda * c1, c1 else: - return nda*c1, c1, nda_additional_stack*c1 + return nda * c1, c1, nda_additional_stack * c1 ## # Gets the data arrays prior to intensity correction. @@ -374,9 +378,15 @@ def _get_data_arrays_prior_to_intensity_correction(self): nda_reference = sitk.GetArrayFromImage(self._reference.sitk) if self._use_reference_mask: - nda_mask = sitk.GetArrayFromImage(self._reference.sitk_mask) + nda_mask_ref = sitk.GetArrayFromImage(self._reference.sitk_mask) + else: + nda_mask_ref = np.ones_like(nda) + + if self._use_stack_mask: + nda_mask_stack = sitk.GetArrayFromImage(self._stack.sitk_mask) else: - nda_mask = np.ones_like(nda) + nda_mask_stack = np.ones_like(nda) + nda_mask = nda_mask_ref * nda_mask_stack if self._additional_stack is None: nda_additional_stack = None @@ -403,4 +413,25 @@ def _create_stack_from_corrected_intensity_array(self, nda, stack): image_sitk = sitk.GetImageFromArray(nda) image_sitk.CopyInformation(stack.sitk) - return st.Stack.from_sitk_image(image_sitk, stack.get_filename(), stack.sitk_mask) + # Update registration history of stack + stack_ic = st.Stack.from_sitk_image( + image_sitk, stack.get_filename(), stack.sitk_mask) + + stack_ic.set_registration_history( + stack.get_registration_history()) + + # Update registration history of slices + slices = stack.get_slices() + slices_ic = stack_ic.get_slices() + kept_slice_numbers = [s.get_slice_number() for s in slices] + for i in range(len(slices_ic)): + # Update registration of kept slice + if i in kept_slice_numbers: + index = kept_slice_numbers.index(i) + slices_ic[i].set_registration_history( + slices[index].get_registration_history()) + # Otherwise, delete slices + else: + stack_ic.delete_slice(i) + + return stack_ic diff --git a/niftymic/utilities/robust_motion_estimator.py b/niftymic/utilities/robust_motion_estimator.py new file mode 100644 index 00000000..eb71686c --- /dev/null +++ b/niftymic/utilities/robust_motion_estimator.py @@ -0,0 +1,252 @@ +## +# \file robust_motion_estimator.py +# \brief Class to estimate motion parameters from estimated +# transformations. +# +# Regularisation of estimated motion parameters for robust slice-motion +# estimates +# +# \author Michael Ebner (michael.ebner.14@ucl.ac.uk) +# \date Feb 2018 +# + +import os +import scipy +import pymc3 +import theano +import numpy as np +import SimpleITK as sitk +import matplotlib.pyplot as plt + +import pysitk.python_helper as ph +import pysitk.simple_itk_helper as sitkh + + +## +# Class to compute robust slice-motion estimates +# \date 2018-03-26 16:47:47-0600 +# +class RobustMotionEstimator(object): + + def __init__(self, transforms_sitk, interleave=2): + self._transforms_sitk = transforms_sitk + self._interleave = interleave + + self._robust_transforms_sitk = [None] * len(self._transforms_sitk) + + ## + # Gets the robust transforms sitk. + # \date 2018-03-26 16:46:00-0600 + # + # \param self The object + # + # \return The robust transforms sitk as list of sitk.Transforms + # + def get_robust_transforms_sitk(self): + robust_transforms_sitk = [ + sitkh.copy_transform_sitk(t) for t in self._robust_transforms_sitk] + return robust_transforms_sitk + + ## + # Run Gaussian process smoothing for each dof individually + # \date 2018-03-13 19:31:02+0000 + # \see http://docs.pymc.io/notebooks/GP-smoothing.html + # + # \param self The object + # + def run_gaussian_process_smoothing(self, smoothing=0.5): + + params_nda = self._get_transformation_params_nda(self._transforms_sitk) + + # Iterate over each dof + for i_dof in range(params_nda.shape[0]): + + # Smooth each interleave package separately + for i in range(self._interleave): + indices = np.arange(i, params_nda.shape[1], self._interleave) + y = params_nda[i_dof, indices] + y_smoothed = self._run_gaussian_process_smoothing( + y=y, smoothing=smoothing) + params_nda[i_dof, indices] = y_smoothed + + self._update_robust_transforms_sitk_from_parameters(params_nda) + + def _run_gaussian_process_smoothing(self, y, smoothing): + LARGE_NUMBER = 1e5 + model = pymc3.Model() + with model: + smoothing_param = theano.shared(smoothing) + mu = pymc3.Normal("mu", sd=LARGE_NUMBER) + tau = pymc3.Exponential("tau", 1.0 / LARGE_NUMBER) + z = pymc3.distributions.timeseries.GaussianRandomWalk( + "z", + mu=mu, + tau=tau / (1.0 - smoothing_param), + shape=y.shape, + ) + obs = pymc3.Normal( + "obs", + mu=z, + tau=tau / smoothing_param, + observed=y, + ) + res = pymc3.find_MAP(vars=[z], method="L-BFGS-B") + return res['z'] + + # ## + # # { function_description } + # # \date 2018-03-13 19:23:39+0000 + # # \see http://docs.pymc.io/notebooks/GP-slice-sampling.html + # # + # # \param self The object + # # + # # \return { description_of_the_return_value } + # # + # def run_gaussian_process_regression(self): + # params_nda = self._get_transformation_params_nda(self._transforms_sitk) + + # # number of training points + # n = params_nda.shape[1] + # X0 = np.arange(0, params_nda.shape[1])[:, None] + + # # Number of points at which to interpolate + # X = np.arange(0, params_nda.shape[1])[:, None] + + # # Covariance kernel parameters + # noise = 0.1 + # lengthscale = 0.3 + # f_scale = 1 + + # cov = f_scale * pymc3.gp.cov.ExpQuad(1, lengthscale) + + # K = cov(X0) + # K_s = cov(X0, X) + # K_noise = K + noise * theano.tensor.eye(n) + + # # Add very slight perturbation to the covariance matrix diagonal to + # # improve numerical stability + # K_stable = K + 1e-12 * theano.tensor.eye(n) + + # # Observed data + # f = np.random.multivariate_normal(mean=np.zeros(n), cov=K_noise.eval()) + + # fig, ax = plt.subplots(figsize=(14, 6)) + # ax.scatter(X0, f, s=40, color='b', label='True points') + # ax.set_xticks(X0) + + # # Analytically compute posterior mean + # L = np.linalg.cholesky(K_noise.eval()) + # alpha = np.linalg.solve(L.T, np.linalg.solve(L, f)) + # post_mean = np.dot(K_s.T.eval(), alpha) + + # ax.plot(X, post_mean, color='g', alpha=0.8, label='Posterior mean') + + # ax.legend() + + # plt.show(True) + + ## + # Shows the estimated transform parameters. + # \date 2018-03-26 16:45:27-0600 + # + # \param self The object + # \param title The title + # \param fullscreen The fullscreen + # + def show_estimated_transform_parameters( + self, dir_output=None, title="RobustMotionEstimator", fullscreen=1): + params_nda = self._get_transformation_params_nda(self._transforms_sitk) + robust_params_nda = self._get_transformation_params_nda( + self.get_robust_transforms_sitk()) + + dof = params_nda.shape[0] + + N_rows = np.ceil(dof / 2.) + i_ref_marker = 0 + + fig = plt.figure(title) + fig.clf() + for i_dof in range(dof): + x = np.arange(params_nda.shape[1]) + y1 = params_nda[i_dof, :] + y2 = robust_params_nda[i_dof, :] + + ax = plt.subplot(N_rows, 2, i_dof + 1) + ax.plot(x, y1, + marker=ph.MARKERS[i_ref_marker], + color=ph.COLORS_TABLEAU20[0], + linestyle=":", + label="original", + markerfacecolor="w", + ) + ax.plot(x, y2, + marker=ph.MARKERS[i_ref_marker], + color=ph.COLORS_TABLEAU20[2], + linestyle="-.", + label="robust", + ) + ax.set_xticks(x) + plt.ylabel(sitkh.TRANSFORM_SITK_DOF_LABELS_LONG[dof][i_dof]) + plt.legend(loc="best") + plt.xlabel('Slice') + plt.suptitle(title) + + if fullscreen: + try: + # Open windows (and also save them) in full screen + manager = plt.get_current_fig_manager() + manager.full_screen_toggle() + except: + pass + + plt.show(block=False) + + if dir_output is not None: + filename = "%s.pdf" % title + ph.save_fig(fig, dir_output, filename) + + ## + # Get transformation parameters from sitk transform + # \date 2018-03-26 16:43:20-0600 + # + # \param self The object + # \param transforms_sitk List of sitk.Transforms + # + # \return The transformation parameters as (dof x #slices)-numpy array + # + def _get_transformation_params_nda(self, transforms_sitk): + dof = len(transforms_sitk[0].GetParameters()) + N_transformations = len(transforms_sitk) + + params_nda = np.zeros((dof, N_transformations)) + for i, transform_sitk in enumerate(transforms_sitk): + params_nda[:, i] = np.array(transform_sitk.GetParameters()) + + # params_nda = self._apply_interleave(params_nda) + + return params_nda + + ## + # Update robust transformations given parameter estimates + # \date 2018-03-26 15:52:19-0600 + # + # \param self The object + # \param params_nda The parameters as (dof x #slices)-numpy array + # + def _update_robust_transforms_sitk_from_parameters(self, params_nda): + + # params_nda = self._undo_interleave(params_nda) + + for i, transform_sitk in enumerate(self._transforms_sitk): + robust_transforms_sitk = sitkh.copy_transform_sitk(transform_sitk) + robust_transforms_sitk.SetParameters(params_nda[:, i]) + + self._robust_transforms_sitk[i] = robust_transforms_sitk + + # def _apply_interleave(self, params_nda): + + # indices = [] + # for i in range(self._interleave): + # indices.append(np.arange(i, self._interleave, params_nda.shape[1])) + + # def _undo_interleave(self, params_nda): diff --git a/niftymic/utilities/stack_mask_morphological_operations.py b/niftymic/utilities/stack_mask_morphological_operations.py index d64c84e5..9436962b 100644 --- a/niftymic/utilities/stack_mask_morphological_operations.py +++ b/niftymic/utilities/stack_mask_morphological_operations.py @@ -64,6 +64,9 @@ def from_stack(cls, stack=None, dilation_radius=0, dilation_kernel="Ball", use_d def set_mask_sitk(self, mask_sitk): self._mask_sitk = mask_sitk + def get_mask_sitk(self): + return self._mask_sitk + def get_stack(self): return st.Stack.from_stack(self._stack) diff --git a/niftymic/utilities/target_stack_estimator.py b/niftymic/utilities/target_stack_estimator.py new file mode 100644 index 00000000..9ce857ba --- /dev/null +++ b/niftymic/utilities/target_stack_estimator.py @@ -0,0 +1,54 @@ +## +# \file target_stack_estimator.py +# \brief Class to estimate target stack automatically +# +# \author Michael Ebner (michael.ebner.14@ucl.ac.uk) +# \date January 2018 +# + + +# Import libraries +import os +import re +import numpy as np +import SimpleITK as sitk + +import pysitk.simple_itk_helper as sitkh + +## +# Class to estimate target stack automatically +# \date 2018-01-26 16:32:11+0000 +# +class TargetStackEstimator(object): + + def __init__(self): + self._target_stack_index = None + + def get_target_stack_index(self): + return self._target_stack_index + + ## + # Use stack with largest mask volume as target stack + # \date 2018-01-26 16:52:39+0000 + # + # \param cls The cls + # \param file_paths_masks paths to image masks as list of strings + # + @classmethod + def from_masks(cls, file_paths_masks): + target_stack_estimator = cls() + + masks_sitk = [sitkh.read_nifti_image_sitk(str(f), sitk.sitkUInt8) + for f in file_paths_masks] + + # Compute volumes of all masks + volumes = np.zeros(len(masks_sitk)) + for i, mask_sitk in enumerate(masks_sitk): + mask_nda = sitk.GetArrayFromImage(mask_sitk) + spacing = np.array(mask_sitk.GetSpacing()) + volumes[i] = np.sum(mask_nda) * spacing.prod() + + # Get index corresponding to maximum volume stack mask + target_stack_estimator._target_stack_index = np.argmax(volumes) + + return target_stack_estimator \ No newline at end of file diff --git a/niftymic/utilities/template_stack_estimator.py b/niftymic/utilities/template_stack_estimator.py new file mode 100644 index 00000000..4b75c640 --- /dev/null +++ b/niftymic/utilities/template_stack_estimator.py @@ -0,0 +1,94 @@ +## +# \file template_stack_estimator.py +# \brief Class to estimate template stack automatically +# +# \author Michael Ebner (michael.ebner.14@ucl.ac.uk) +# \date January 2018 +# + + +import os +import re +import numpy as np +import json +import SimpleITK as sitk + +import pysitk.simple_itk_helper as sitkh + +from niftymic.definitions import DIR_TEMPLATES, TEMPLATES_INFO + + +## +# Class to estimate template stack automatically +# \date 2018-01-26 16:32:11+0000 +# +class TemplateStackEstimator(object): + + def __init__(self): + self._template_path = None + + ## + # Gets the path to estimated template. + # \date 2018-01-27 02:14:53+0000 + # + # \param self The object + # + # \return The path to template. + # + def get_path_to_template(self): + return self._template_path + + ## + # Select template with similar brain volume + # \date 2018-01-26 16:52:39+0000 + # + # \param cls The cls + # \param file_paths_masks paths to image masks as list of strings + # + @classmethod + def from_mask(cls, file_path_mask): + template_stack_estimator = cls() + + mask_sitk = sitkh.read_nifti_image_sitk( + file_path_mask, sitk.sitkUInt8) + mask_nda = sitk.GetArrayFromImage(mask_sitk) + spacing = np.array(mask_sitk.GetSpacing()) + volume = len(np.where(mask_nda > 0)[0]) * spacing.prod() + + # Read in template info + path_to_template_info = os.path.join(DIR_TEMPLATES, TEMPLATES_INFO) + with open(path_to_template_info) as json_file: + dic = json.load(json_file) + + # Get gestational ages as list of integers + gestational_ages = sorted([int(gw) for gw in dic.keys()]) + + # # Get matching gestational age + # template_volumes = np.array([dic[str(k)]["volume_mask"] + # for k in gestational_ages]) + # index = np.argmin(np.abs(template_volumes - volume)) + + # # Ensure valid index after correction + # index = np.max([0, index - 1]) + # # index = np.min([index + 1, len(template_volumes)-1]) + + # # Matching gestational age/week + # gw_match = str(gestational_ages[index]) + + # template_stack_estimator._template_path = os.path.join( + # DIR_TEMPLATES, dic[gw_match]["image"]) + + # return template_stack_estimator + + # Find template which has slightly smaller mask volume + for k in gestational_ages: + if dic[str(k)]["volume_mask_dil"] > volume: + key = str(np.max([gestational_ages[0], k - 1])) + template_stack_estimator._template_path = os.path.join( + DIR_TEMPLATES, dic[key]["image"]) + return template_stack_estimator + + # Otherwise, return path to oldest template image available + template_stack_estimator._template_path = os.path.join( + DIR_TEMPLATES, dic[str(gestational_ages[-1])]["image"]) + return template_stack_estimator diff --git a/niftymic/utilities/toolkit_executor.py b/niftymic/utilities/toolkit_executor.py new file mode 100644 index 00000000..eadc2af9 --- /dev/null +++ b/niftymic/utilities/toolkit_executor.py @@ -0,0 +1,166 @@ +## +# \file toolkit_executor.py +# \brief generates function calls to execute other reconstruction toolkits +# +# \author Michael Ebner (michael.ebner.14@ucl.ac.uk) +# \date Jan 2018 +# + +import os +import pysitk.python_helper as ph + +EXE_IRTK = { + "workstation": "/home/mebner/Development/VolumetricReconstruction_ImperialCollege/source/bin/SVRreconstructionGPU" +} + + +## +# Class to generate functions calls to execute other reconstruction toolkits +# \date 2018-01-27 02:12:00+0000 +# +class ToolkitExecuter(object): + + def __init__(self, paths_to_images, paths_to_masks, dir_output): + self._paths_to_images = paths_to_images + self._paths_to_masks = paths_to_masks + self._dir_output = dir_output + + # separator for command line export + self._sep = " \\\n" + + self._subdir_temp = "./temp" + + ## + # Gets the function call for fetalReconstruction toolkit provided by + # Bernhard Kainz. + # \date 2018-01-27 02:12:26+0000 + # + # \param self The object + # \param option_args The option arguments + # \param exe The executable + # \param output_name The output name + # + # \return The function call irtk. + # + def get_function_call_irtk( + self, + option_args=['-d 0', '--useCPU', '--resolution 1'], + exe=None, + output_name="IRTK_SRR.nii.gz", + kernel_mask_dilation=None, + ): + if exe is None: + exe = EXE_IRTK["workstation"] + + cmd_args = [] + + # store pwd + cmd_args.append("PWD=$(pwd)") + + # change to output directory + cmd_args.append("\necho 'Change to output directory'") + cmd_args.append("mkdir -p %s" % self._dir_output) + cmd_args.append("cd %s" % self._dir_output) + + # create temp directory if required + cmd_args.append("\necho 'Create temp directory'") + cmd_args.append("mkdir -p %s" % self._subdir_temp) + + # dilate masks + if kernel_mask_dilation is not None: + cmd_args.append("\necho 'Dilate masks'") + cmd_args.append(self._exe_dilate_masks( + kernel_mask_dilation, self._paths_to_masks)) + + # exe to determine slice thickness for toolkit + cmd_args.append("\necho 'Fetch slice thickness for all stacks'") + cmd_args.append(self._exe_to_fetch_slice_thickness( + self._paths_to_images)) + + # toolkit execution + cmd_args.append("\necho 'IRTK Toolkit Execution'") + exe_args = [exe] + exe_args.append("-o %s" % output_name) + exe_args.append("-i %s%s" % + (self._sep, self._sep.join(self._paths_to_images))) + # exe_args.append("--manualMask %s" % self._paths_to_masks[0]) # + # causes cuda sync error!? + exe_args.append("-m %s%s" % + (self._sep, self._sep.join(self._paths_to_masks))) + exe_args.append("--thickness `printf \"%s\" \"${thickness}\"`") + exe_args.extend(option_args) + toolkit_execution = "%s" % self._sep.join(exe_args) + cmd_args.append(toolkit_execution) + + cmd_args.append("\necho 'Delete temp directory'") + cmd_args.append("rm -rf %s" % self._subdir_temp) + + cmd_args.append("\necho 'Change back to original directory'") + cmd_args.append("cd ${PWD}") + cmd_args.append("\n") + + cmd = (" \n").join(cmd_args) + return cmd + + @staticmethod + def write_function_call_to_file(function_call, path_to_file): + text = "#!/bin/zsh\n\n%s" % function_call + ph.write_to_file(path_to_file, text, verbose=False) + ph.execute_command("chmod +x %s" % path_to_file, verbose=False) + + ## + # Provide bash-commands to read out slice thickness on-the-fly + # + # Rationale: IRTK recon toolkit assumes otherwise a thickness of twice the + # voxel spacing by default + # \date 2018-01-27 02:12:52+0000 + # + # \param paths_to_images The paths to images + # + # \return bash command as string + # + def _exe_to_fetch_slice_thickness(self, paths_to_images): + cmd_args = [] + cmd_args.append("args=()") + cmd_args.append("for i in %s" % (self._sep.join(paths_to_images))) + cmd_args.append("do") + cmd_args.append( + "t=$(fslhd ${i} | grep pixdim3 | awk -F ' ' '{print $2}')") + cmd_args.append("args+=(\" ${t}\")") + cmd_args.append("done") + cmd_args.append("thickness=${args[@]}") + cmd_args.append("") + cmd = ("\n").join(cmd_args) + return cmd + + ## + # Provide bash-commands to read out slice thickness on-the-fly + # + # Rationale: IRTK recon toolkit assumes otherwise a thickness of twice the + # voxel spacing by default + # \date 2018-01-27 02:12:52+0000 + # + # \param paths_to_images The paths to images + # + # \return bash command as string + # # + def _exe_dilate_masks(self, kernel, paths_to_masks, label=1): + cmd_loop = [] + kernel_str = [str(k) for k in kernel] + + # Export dilated mask to temp directory + for i_mask, path_to_mask in enumerate(paths_to_masks): + directory = os.path.dirname(path_to_mask) + mask_filename = os.path.basename(path_to_mask) + path_to_mask_dilated = os.path.join( + self._subdir_temp, ph.append_to_filename(mask_filename, "_dil")) + + cmd_args = ["c3d"] + cmd_args.append(path_to_mask) + cmd_args.append("-dilate %s %smm" % (label, "x".join(kernel_str))) + cmd_args.append("-o %s" % path_to_mask_dilated) + cmd = self._sep.join(cmd_args) + + cmd_loop.append(cmd) + paths_to_masks[i_mask] = path_to_mask_dilated + return "\n".join(cmd_loop) diff --git a/niftymic/utilities/volumetric_reconstruction_pipeline.py b/niftymic/utilities/volumetric_reconstruction_pipeline.py index 2044b851..ec079b3c 100644 --- a/niftymic/utilities/volumetric_reconstruction_pipeline.py +++ b/niftymic/utilities/volumetric_reconstruction_pipeline.py @@ -10,6 +10,7 @@ # \date Aug 2017 # +import six import numpy as np import SimpleITK as sitk from abc import ABCMeta, abstractmethod @@ -18,6 +19,9 @@ import pysitk.simple_itk_helper as sitkh import niftymic.base.stack as st +import niftymic.validation.motion_evaluator as me +import niftymic.validation.residual_evaluator as re +import niftymic.utilities.robust_motion_estimator as rme ## @@ -134,7 +138,7 @@ def _run(self): for i in range(0, len(self._stacks)): txt = "Volume-to-Volume Registration -- " \ - "Stack %d/%d" % (i+1, len(self._stacks)) + "Stack %d/%d" % (i + 1, len(self._stacks)) if self._verbose: ph.print_subtitle(txt) else: @@ -168,6 +172,8 @@ class SliceToVolumeRegistration(RegistrationPipeline): # \param verbose The verbose # \param print_prefix Print at each iteration at the # beginning, string + # \param threshold The threshold + # \param threshold_measure The threshold measure # def __init__(self, stacks, @@ -175,6 +181,10 @@ def __init__(self, registration_method, verbose=1, print_prefix="", + threshold=None, + threshold_measure="NCC", + s2v_smoothing=None, + interleave=2, ): RegistrationPipeline.__init__( self, @@ -183,39 +193,146 @@ def __init__(self, registration_method=registration_method, verbose=verbose) self._print_prefix = print_prefix + self._threshold = threshold + self._threshold_measure = threshold_measure + self._s2v_smoothing = s2v_smoothing + self._interleave = interleave def set_print_prefix(self, print_prefix): self._print_prefix = print_prefix + def set_threshold(self, threshold): + self._threshold = threshold + + def get_threshold(self): + return self._threshold + + def set_s2v_smoothing(self, s2v_smoothing): + self._s2v_smoothing = s2v_smoothing + + def get_s2v_smoothing(self): + return self._s2v_smoothing + + def set_threshold_measure(self, threshold_measure): + self._threshold_measure = threshold_measure + + def get_threshold_measure(self): + return self._threshold_measure + def _run(self): ph.print_title("Slice-to-Volume Registration") self._registration_method.set_moving(self._reference) - for i in range(0, len(self._stacks)): - stack = self._stacks[i] + for i, stack in enumerate(self._stacks): slices = stack.get_slices() - for j in range(0, len(slices)): + transforms_sitk = [None] * len(slices) + + for j, slice_j in enumerate(slices): txt = "%sSlice-to-Volume Registration -- " \ "Stack %d/%d -- Slice %d/%d" % ( self._print_prefix, - i+1, len(self._stacks), - j+1, len(slices)) + i + 1, len(self._stacks), + j + 1, len(slices)) if self._verbose: ph.print_subtitle(txt) else: ph.print_info(txt) - self._registration_method.set_fixed(slices[j]) + self._registration_method.set_fixed(slice_j) self._registration_method.run() - # Update position of slice + # Store information on registration transform transform_sitk = \ self._registration_method.\ get_registration_transform_sitk() - slices[j].update_motion_correction(transform_sitk) + transforms_sitk[j] = transform_sitk + + # Avoid slice misregistrations + if self._s2v_smoothing is not None: + ph.print_subtitle( + "Robust slice motion estimation " + "(GP smoothing = %g, interleave = %d)" % ( + self._s2v_smoothing, self._interleave)) + robust_motion_estimator = rme.RobustMotionEstimator( + transforms_sitk=transforms_sitk, + interleave=self._interleave) + robust_motion_estimator.run_gaussian_process_smoothing( + self._s2v_smoothing) + transforms_sitk = \ + robust_motion_estimator.get_robust_transforms_sitk() + + # Export figures + title = "%s_Stack%d%s" % ( + self._print_prefix, i, stack.get_filename()) + title = ph.replace_string_for_print(title) + robust_motion_estimator.show_estimated_transform_parameters( + dir_output="/tmp/fetal_brain/figs", title=title) + + # dir_output = "/tmp/fetal/figs" + # motion_evaluator = me.MotionEvaluator(transforms_sitk) + # motion_evaluator.run() + # motion_evaluator.display(dir_output=dir_output, title=title) + # motion_evaluator.show(dir_output=dir_output, title=title) + + # Update position of slice + for j, slice_j in enumerate(slices): + slice_j.update_motion_correction(transforms_sitk[j]) + + # Reject misregistered slices + if self._threshold is not None: + ph.print_subtitle( + "Slice Outlier Rejection (Threshold = %g @ %s)" % ( + self._threshold, self._threshold_measure)) + residual_evaluator = re.ResidualEvaluator( + stacks=self._stacks, + reference=self._reference, + use_slice_masks=False, + use_reference_mask=True, + verbose=False, + ) + residual_evaluator.compute_slice_projections() + residual_evaluator.evaluate_slice_similarities() + slice_sim = residual_evaluator.get_slice_similarities() + # residual_evaluator.show_slice_similarities( + # threshold=self._threshold, + # measures=[self._threshold_measure], + # directory="/tmp/spina/figs%s" % self._print_prefix[0:7], + # ) + + remove_stacks = [] + for i, stack in enumerate(self._stacks): + nda_sim = np.nan_to_num( + slice_sim[stack.get_filename()][self._threshold_measure]) + indices = np.where(nda_sim < self._threshold)[0] + N_slices = len(stack.get_slices()) + for j in indices: + stack.delete_slice(j) + + ph.print_info("Stack %d/%d: %d/%d slices deleted (%s)" % ( + i + 1, + len(self._stacks), + len(indices), + N_slices, + stack.get_filename(), + )) + + # Log stack where all slices were rejected + if stack.get_number_of_slices() == 0: + remove_stacks.append(stack) + + # Remove stacks where all slices where rejected + for stack in remove_stacks: + self._stacks.remove(stack) + ph.print_info("Stack '%s' removed entirely." % + stack.get_filename()) + + if len(self._stacks) == 0: + raise RuntimeError( + "All slices of all stacks were rejected " + "as outliers. Volumetric reconstruction is aborted.") ## @@ -264,7 +381,7 @@ def _run(self, debug=0): txt = "%sSliceSet-to-Volume Registration -- " \ "Stack %d/%d -- Slices %s" % ( self._print_prefix, - i+1, len(self._stacks), + i + 1, len(self._stacks), str(indices)) if self._verbose: ph.print_subtitle(txt) @@ -444,6 +561,12 @@ def __init__(self, alpha_range, cycles, verbose=1, + outlier_rejection=False, + threshold_measure="NCC", + threshold_range=[0.6, 0.7], + use_robust_registration=False, + s2v_smoothing=0.5, + interleave=2, ): ReconstructionRegistrationPipeline.__init__( @@ -456,6 +579,12 @@ def __init__(self, verbose=verbose) self._cycles = cycles + self._outlier_rejection = outlier_rejection + self._threshold_measure = threshold_measure + self._threshold_range = threshold_range + self._use_robust_registration = use_robust_registration + self._s2v_smoothing = s2v_smoothing + self._interleave = interleave def _run(self): @@ -467,11 +596,18 @@ def _run(self): self._alpha_range[0], self._alpha_range[1], self._cycles) alphas = alphas[0:self._cycles] + thresholds = np.linspace( + self._threshold_range[0], self._threshold_range[1], self._cycles) + thresholds = thresholds[0:self._cycles] + s2vreg = SliceToVolumeRegistration( stacks=self._stacks, reference=self._reference, registration_method=self._registration_method, - verbose=self._verbose) + verbose=self._verbose, + threshold_measure=self._threshold_measure, + interleave=self._interleave, + ) reference = self._reference @@ -479,7 +615,14 @@ def _run(self): # Slice-to-volume registration step s2vreg.set_reference(reference) - s2vreg.set_print_prefix("Cycle %d/%d: " % (cycle+1, self._cycles)) + s2vreg.set_print_prefix("Cycle %d/%d: " % + (cycle + 1, self._cycles)) + if self._outlier_rejection: + s2vreg.set_threshold(thresholds[cycle]) + if self._use_robust_registration and cycle == 0: + s2vreg.set_s2v_smoothing(self._s2v_smoothing) + else: + s2vreg.set_s2v_smoothing(None) s2vreg.run() self._computational_time_registration += \ @@ -497,7 +640,7 @@ def _run(self): # Store SRR filename = "Iter%d_%s" % ( - cycle+1, + cycle + 1, self._reconstruction_method.get_setting_specific_filename() ) self._reconstructions.insert(0, st.Stack.from_stack( @@ -569,9 +712,9 @@ def _run(self, debug=1): # Debug if debug: for i, stack in enumerate(self._stacks): - print("Stack %d/%d:" % (i+1, N_stacks)) - for k, v in slice_sets_indices[i].iteritems(): - print("\tCycle %d: arrays = %s" % (k+1, str(v))) + print("Stack %d/%d:" % (i + 1, N_stacks)) + for k, v in six.iteritems(slice_sets_indices[i]): + print("\tCycle %d: arrays = %s" % (k + 1, str(v))) N_cycles = np.max([len(slice_sets_indices[i]) for i in range(N_stacks)]) @@ -592,7 +735,7 @@ def _run(self, debug=1): } ss2vreg = SliceSetToVolumeRegistration( - print_prefix="Cycle %d/%d -- " % (i_cycle+1, N_cycles), + print_prefix="Cycle %d/%d -- " % (i_cycle + 1, N_cycles), stacks=self._stacks, reference=reference, registration_method=self._registration_method, @@ -675,11 +818,11 @@ def _get_slice_set_indices_per_cycle(self, stack, N_min): # Split into smaller subpackages while not finished: - i = i+1 + i = i + 1 # Get list of indices based on interleaved acquisition interleaved_acquisitions[i] = self._get_array_list_split( - interleaved_acquisitions[i-1], N_min) + interleaved_acquisitions[i - 1], N_min) # Stop if number of elements smaller than N_min. Remark, single # index splits can still occur. E.g. [1,3,5] is split into [1,3] @@ -763,7 +906,7 @@ def _run(self): for i in range(0, len(self._stacks)): ph.print_subtitle("Multi-Component Reconstruction -- " - "Stack %d/%d" % (i+1, len(self._stacks))) + "Stack %d/%d" % (i + 1, len(self._stacks))) stack = self._stacks[i] self._reconstruction_method.set_stacks([stack]) self._reconstruction_method.run() diff --git a/niftymic/validation/evaluate_image_similarity.py b/niftymic/validation/evaluate_image_similarity.py new file mode 100644 index 00000000..92f2243f --- /dev/null +++ b/niftymic/validation/evaluate_image_similarity.py @@ -0,0 +1,129 @@ +## +# \file evaluate_image_similarity.py +# \brief Evaluate similarity to a reference of one or more images +# +# \author Michael Ebner (michael.ebner.14@ucl.ac.uk) +# \date Feb 1 +# + +import os +import numpy as np +import pandas as pd +import SimpleITK as sitk + +import nsol.observer as obs +from nsol.similarity_measures import SimilarityMeasures as \ + SimilarityMeasures + +import niftymic.base.stack as st +import niftymic.base.data_reader as dr +import niftymic.registration.niftyreg as regniftyreg +from niftymic.utilities.input_arparser import InputArgparser + +import pysitk.python_helper as ph +import pysitk.simple_itk_helper as sitkh + + +def main(): + + # Set print options + np.set_printoptions(precision=3) + pd.set_option('display.width', 1000) + + input_parser = InputArgparser( + description=".", + ) + input_parser.add_filenames(required=True) + input_parser.add_reference(required=True) + input_parser.add_reference_mask() + input_parser.add_dir_output(required=False) + input_parser.add_measures(default=["PSNR", "RMSE", "SSIM", "NCC", "NMI"]) + input_parser.add_verbose(default=0) + args = input_parser.parse_args() + input_parser.print_arguments(args) + + ph.print_title("Image similarity") + data_reader = dr.MultipleImagesReader(args.filenames) + data_reader.read_data() + stacks = data_reader.get_data() + + reference = st.Stack.from_filename(args.reference, args.reference_mask) + + for stack in stacks: + try: + stack.sitk - reference.sitk + except RuntimeError as e: + raise IOError( + "All provided images must be at the same image space") + + x_ref = sitk.GetArrayFromImage(reference.sitk) + + if args.reference_mask is None: + indices = np.where(x_ref != np.inf) + else: + x_ref_mask = sitk.GetArrayFromImage(reference.sitk_mask) + indices = np.where(x_ref_mask > 0) + + measures_dic = { + m: lambda x, m=m: + SimilarityMeasures.similarity_measures[m]( + x[indices], x_ref[indices]) + # SimilarityMeasures.similarity_measures[m](x, x_ref) + for m in args.measures + } + + observer = obs.Observer() + observer.set_measures(measures_dic) + for stack in stacks: + nda = sitk.GetArrayFromImage(stack.sitk) + observer.add_x(nda) + + if args.verbose: + stacks_comparison = [s for s in stacks] + stacks_comparison.insert(0, reference) + sitkh.show_stacks( + stacks_comparison, + segmentation=reference, + ) + + observer.compute_measures() + measures = observer.get_measures() + + # Store information in array + error = np.zeros((len(stacks), len(measures))) + cols = measures + rows = [] + for i_stack, stack in enumerate(stacks): + error[i_stack, :] = np.array([measures[m][i_stack] for m in measures]) + rows.append(stack.get_filename()) + + header = "# Ref: %s, Ref-Mask: %d, %s \n" % ( + reference.get_filename(), + args.reference_mask is None, + ph.get_time_stamp(), + ) + header += "# %s\n" % ("\t").join(measures) + + path_to_file_filenames = os.path.join( + args.dir_output, "filenames.txt") + path_to_file_similarities = os.path.join( + args.dir_output, "similarities.txt") + + # Write to files + ph.write_to_file(path_to_file_similarities, header) + ph.write_array_to_file( + path_to_file_similarities, error, verbose=False) + text = header + text += "%s\n" %"\n".join(rows) + ph.write_to_file(path_to_file_filenames, text) + + # Print to screen + ph.print_subtitle("Computed Similarities") + df = pd.DataFrame(error, rows, cols) + print(df) + + return 0 + + +if __name__ == '__main__': + main() diff --git a/niftymic/validation/evaluate_simulated_stack_similarity.py b/niftymic/validation/evaluate_simulated_stack_similarity.py index 055eb9fa..7d34c650 100644 --- a/niftymic/validation/evaluate_simulated_stack_similarity.py +++ b/niftymic/validation/evaluate_simulated_stack_similarity.py @@ -106,9 +106,13 @@ def main(): text += "\n" ph.write_to_file(path_to_file, text, "w") for k in range(nda_3D_original.shape[0]): - x_2D_original = nda_3D_original[k, :, :].flatten() - x_2D_simulated = nda_3D_simulated[k, :, :].flatten() - x_2D_mask = nda_3D_mask[k, :, :].flatten() + x_2D_original = nda_3D_original[k, :, :] + x_2D_simulated = nda_3D_simulated[k, :, :] + + # zero slice, i.e. rejected during motion correction + if np.abs(x_2D_simulated).sum() < 1e-6: + x_2D_simulated[:] = np.nan + x_2D_mask = nda_3D_mask[k, :, :] indices = np.where(x_2D_mask > 0) @@ -117,7 +121,7 @@ def main(): similarities[m] = similarity_measures[measure]( x_2D_original[indices], x_2D_simulated[indices]) else: - similarities[m] = np.inf + similarities[m] = np.nan ph.write_array_to_file(path_to_file, similarities.reshape(1, -1)) return 0 diff --git a/niftymic/validation/evaluate_slice_residual_similarity.py b/niftymic/validation/evaluate_slice_residual_similarity.py new file mode 100644 index 00000000..e133c22a --- /dev/null +++ b/niftymic/validation/evaluate_slice_residual_similarity.py @@ -0,0 +1,100 @@ +## +# \file evaluate_slice_residual_similarity.py +# \brief Evaluate similarity to a reference of one or more images +# +# \author Michael Ebner (michael.ebner.14@ucl.ac.uk) +# \date Feb 1 +# + +import os +import numpy as np +import pandas as pd +import SimpleITK as sitk + +import niftymic.base.stack as st +import niftymic.base.data_reader as dr +import niftymic.registration.niftyreg as regniftyreg +from niftymic.utilities.input_arparser import InputArgparser +import niftymic.validation.residual_evaluator as res_ev + +import pysitk.python_helper as ph +import pysitk.simple_itk_helper as sitkh + + +def main(): + + time_start = ph.start_timing() + + # Set print options + np.set_printoptions(precision=3) + pd.set_option('display.width', 1000) + + input_parser = InputArgparser( + description=".", + ) + input_parser.add_dir_input() + input_parser.add_filenames() + input_parser.add_filenames_masks() + input_parser.add_suffix_mask(default="_mask") + input_parser.add_reference(required=True) + input_parser.add_reference_mask() + input_parser.add_dir_output(required=False) + input_parser.add_measures(default=["PSNR", "RMSE", "SSIM", "NCC", "NMI"]) + input_parser.add_verbose(default=0) + input_parser.add_option( + option_string="--use-reference-mask", type=int, default=1) + input_parser.add_option( + option_string="--use-slice-masks", type=int, default=1) + + args = input_parser.parse_args() + input_parser.print_arguments(args) + + # Neither '--dir-input' nor '--filenames' was specified + if args.filenames is not None and args.dir_input is not None: + raise IOError( + "Provide input by either '--dir-input' or '--filenames' " + "but not both together") + + # '--dir-input' specified + elif args.dir_input is not None: + data_reader = dr.ImageSlicesDirectoryReader( + args.dir_input, suffix_mask=args.suffix_mask) + + # '--filenames' specified + elif args.filenames is not None: + data_reader = dr.MultipleImagesReader( + file_paths=args.filenames, + file_paths_masks=args.filenames_masks, + suffix_mask=args.suffix_mask) + + else: + raise IOError( + "Provide input by either '--dir-input' or '--filenames'") + + ph.print_title("Slice Residual Similarity") + data_reader.read_data() + stacks = data_reader.get_data() + + reference = st.Stack.from_filename(args.reference, args.reference_mask) + + residual_evaluator = res_ev.ResidualEvaluator( + stacks=stacks, + reference=reference, + measures=args.measures, + use_reference_mask=args.use_reference_mask, + use_slice_masks=args.use_slice_masks, + ) + residual_evaluator.compute_slice_projections() + residual_evaluator.evaluate_slice_similarities() + residual_evaluator.write_slice_similarities(args.dir_output) + + elapsed_time = ph.stop_timing(time_start) + ph.print_title("Summary") + print("Computational Time for Slice Residual Evaluation: %s" % + (elapsed_time)) + + return 0 + + +if __name__ == '__main__': + main() diff --git a/niftymic/validation/export_side_by_side_simulated_vs_original_slice_comparison.py b/niftymic/validation/export_side_by_side_simulated_vs_original_slice_comparison.py index 7e2ad9cc..2c9248bc 100644 --- a/niftymic/validation/export_side_by_side_simulated_vs_original_slice_comparison.py +++ b/niftymic/validation/export_side_by_side_simulated_vs_original_slice_comparison.py @@ -43,7 +43,7 @@ def export_comparison_to_file(nda_original, dir_tmp = os.path.join(DIR_TMP, "ImageMagick") ph.clear_directory(dir_tmp, verbose=False) for k in range(nda_original.shape[0]): - ctr = k+1 + ctr = k + 1 # Export as individual image side-by-side _export_image_side_by_side( @@ -94,6 +94,8 @@ def _export_image_side_by_side( path_to_left = os.path.join(dir_output, "left.%s" % extension) path_to_right = os.path.join(dir_output, "right.%s" % extension) + nda_left = np.round(np.array(nda_left)).astype(np.uint8) + nda_right = np.round(np.array(nda_right)).astype(np.uint8) ph.write_image(nda_left, path_to_left, verbose=False) ph.write_image(nda_right, path_to_right, verbose=False) @@ -224,7 +226,7 @@ def main(): intensity_max = 255 intensity_min = 0 for i in range(len(stacks_original)): - ph.print_subtitle("Stack %d/%d" % (i+1, len(stacks_original))) + ph.print_subtitle("Stack %d/%d" % (i + 1, len(stacks_original))) nda_3D_original = sitk.GetArrayFromImage(stacks_original[i].sitk) nda_3D_simulated = sitk.GetArrayFromImage(stacks_simulated[i].sitk) diff --git a/niftymic/validation/image_similarity_evaluator.py b/niftymic/validation/image_similarity_evaluator.py new file mode 100644 index 00000000..f0c76216 --- /dev/null +++ b/niftymic/validation/image_similarity_evaluator.py @@ -0,0 +1,234 @@ +## +# \file image_similarity_evaluator.py +# \brief Class to evaluate image similarity between stacks and a reference +# +# \author Michael Ebner (michael.ebner.14@ucl.ac.uk) +# \date February 2018 +# + + +# Import libraries +import os +import re +import numpy as np +import SimpleITK as sitk + +import nsol.observer as obs +from nsol.similarity_measures import SimilarityMeasures as \ + SimilarityMeasures +import pysitk.python_helper as ph + +import niftymic.reconstruction.linear_operators as lin_op +import niftymic.base.exceptions as exceptions + + +## +# Class to evaluate image similarity between stacks and a reference +# \date 2018-02-08 16:16:08+0000 +# +class ImageSimilarityEvaluator(object): + + ## + # { constructor_description } + # \date 2018-02-08 14:13:19+0000 + # + # \param self The object + # \param stacks List of Stack objects + # \param reference Reference as Stack object + # \param use_reference_mask The use reference mask + # \param measures Similarity measures as given in + # nsol.similarity_measures, list of strings + # \param verbose The verbose + # + def __init__( + self, + stacks=None, + reference=None, + use_reference_mask=True, + measures=["NCC", "NMI", "PSNR", "SSIM", "RMSE"], + verbose=True, + ): + self._stacks = stacks + self._reference = reference + self._measures = measures + self._use_reference_mask = use_reference_mask + self._verbose = verbose + + self._similarities = None + + self._filename_filenames = "filenames.txt" + self._filename_similarities = "similarities.txt" + + ## + # Sets the stacks. + # \date 2018-02-08 14:13:27+0000 + # + # \param self The object + # \param stacks List of Stack objects + # + def set_stacks(self, stacks): + self._stacks = stacks + + ## + # Sets the reference from which the slices shall be simulated/projected. + # \date 2018-01-19 17:26:14+0000 + # + # \param self The object + # \param reference The reference + # + def set_reference(self, reference): + self._reference = reference + + ## + # Gets the computed similarities. + # \date 2018-02-08 14:36:15+0000 + # + # \param self The object + # + # \return The similarities as dictionary. E.g. { "NCC": np.array()} + # + def get_similarities(self): + return self._similarities + + def get_measures(self): + return self._measures + + ## + # Calculates the similarities. Outcome can be fetched using + # 'get_similarities' + # \date 2018-02-08 16:16:41+0000 + # + # \param self The object + # + # \post self._similarities updated. + # + def compute_similarities(self): + for stack in self._stacks: + try: + stack.sitk - self._reference.sitk + except RuntimeError as e: + raise IOError( + "All provided images must be at the same image space") + + x_ref = sitk.GetArrayFromImage(self._reference.sitk) + + x_ref_mask = np.ones_like(x_ref) + if self._use_reference_mask: + x_ref_mask *= sitk.GetArrayFromImage(self._reference.sitk_mask) + indices = np.where(x_ref_mask > 0) + + if len(indices[0]) == 0: + raise RuntimeError( + "Support to evaluate similarity measures is zero") + + # Define similarity measures as dic + measures_dic = { + m: lambda x, m=m: + SimilarityMeasures.similarity_measures[m]( + x[indices], x_ref[indices]) + # SimilarityMeasures.similarity_measures[m](x, x_ref) + for m in self._measures + } + + # Compute similarities + observer = obs.Observer() + observer.set_measures(measures_dic) + for stack in self._stacks: + nda = sitk.GetArrayFromImage(stack.sitk) + observer.add_x(nda) + observer.compute_measures() + self._similarities = observer.get_measures() + + # Add filenames to dictionary + image_names = [s.get_filename() for s in self._stacks] + self._similarities["filenames"] = image_names + + ## + # Writes the evaluated similarities to two files; one containing the + # similarity information, the other the filename information. + # \date 2018-02-08 14:58:29+0000 + # + # \param self The object + # \param directory The directory + # + def write_similarities(self, directory): + + # Store information in array + similarities_nda = np.zeros((len(self._stacks), len(self._measures))) + filenames = [] + for i_stack, stack in enumerate(self._stacks): + similarities_nda[i_stack, :] = np.array( + [self._similarities[m][i_stack] for m in self._measures]) + filenames.append(stack.get_filename()) + + # Build header of files + header = "# Ref: %s, Ref-Mask: %d, %s \n" % ( + self._reference.get_filename(), + self._use_reference_mask, + ph.get_time_stamp(), + ) + header += "# %s\n" % ("\t").join(self._measures) + + # Get filename paths + path_to_file_filenames, path_to_file_similarities = self._get_filename_paths( + directory) + + # Write similarities + ph.write_to_file(path_to_file_similarities, header) + ph.write_array_to_file( + path_to_file_similarities, similarities_nda, verbose=self._verbose) + + # Write stack filenames + text = header + text += "%s\n" % "\n".join(filenames) + ph.write_to_file(path_to_file_filenames, text, verbose=self._verbose) + + ## + # Reads similarities. + # \date 2018-02-08 15:32:04+0000 + # + # \param self The object + # \param directory The directory + # + def read_similarities(self, directory): + + if not ph.directory_exists(directory): + raise IOError("Directory '%s' does not exist." % directory) + + # Get filename paths + path_to_file_filenames, path_to_file_similarities = self._get_filename_paths( + directory) + + for f in [path_to_file_filenames, path_to_file_similarities]: + if not ph.file_exists(path_to_file_filenames): + raise IOError("File '%s' does not exist" % f) + + lines = ph.read_file_line_by_line(path_to_file_filenames) + + # Get image filenames + image_names = [re.sub("\n", "", f) for f in lines[2:]] + + # Get computed measures + measures = lines[1] + measures = re.sub("# ", "", measures) + measures = re.sub("\n", "", measures) + self._measures = measures.split("\t") + + # Get computed similarities + similarities_nda = np.loadtxt(path_to_file_similarities, skiprows=2) + + # Reconstruct similarity dictionary + self._similarities = {} + self._similarities["filenames"] = image_names + for i_m, m in enumerate(self._measures): + self._similarities[m] = similarities_nda[:, i_m] + + def _get_filename_paths(self, directory): + + # Define filename paths + path_to_file_filenames = os.path.join( + directory, self._filename_filenames) + path_to_file_similarities = os.path.join( + directory, self._filename_similarities) + + return path_to_file_filenames, path_to_file_similarities diff --git a/niftymic/validation/motion_evaluator.py b/niftymic/validation/motion_evaluator.py new file mode 100644 index 00000000..a5519644 --- /dev/null +++ b/niftymic/validation/motion_evaluator.py @@ -0,0 +1,176 @@ +## +# \file motion_evaluator.py +# \brief Class to evaluate computed motions +# +# Should help to assess the registration accuracy. +# \author Michael Ebner (michael.ebner.14@ucl.ac.uk) +# \date January 2018 +# + + +# Import libraries +import os +import re +import numpy as np +import pandas as pd +import SimpleITK as sitk +import matplotlib.pyplot as plt + +import pysitk.python_helper as ph +import pysitk.simple_itk_helper as sitkh + +import niftymic.base.exceptions as exceptions + + +## +# Class to evaluate computed motions +# \date 2018-01-25 22:53:37+0000 +# +class MotionEvaluator(object): + + def __init__(self, transforms_sitk): + self._transforms_sitk = transforms_sitk + self._transform_params = None + + self._scale = { + # convert radiant to degree for angles + 6: np.array([180. / np.pi, 180. / np.pi, 180. / np.pi, 1, 1, 1]), + } + self._labels_long = { + 6: ["angle_x [deg]", + "angle_y [deg]", + "angle_z [deg]", + "t_x [mm]", + "t_y [mm]", + "t_z [mm]"], + } + self._labels_short = { + 6: ["angle_x", + "angle_y", + "angle_z", + "t_x", + "t_y", + "t_z"], + } + + def run(self): + + # # Eliminate center information + # if self._transforms_sitk[0].GetName() in \ + # ["Euler2DTransform", "Euler3DTransform"]: + # identity = eval("sitk.Euler%dDTransform()" + # transforms_sitk = [ + # sitkh.get_composite_sitk_euler_transform() + # ] + + # Create (#transforms x DOF) numpy array + self._transform_params = np.zeros(( + len(self._transforms_sitk), + len(self._transforms_sitk[0].GetParameters()) + )) + + for j in range(self._transform_params.shape[0]): + self._transform_params[j, :] = \ + self._transforms_sitk[j].GetParameters() + + def display(self, title=None, dir_output=None): + pd.set_option('display.width', 1000) + N_trafos, dof = self._transform_params.shape + if dof == 6: + params = self._get_scaled_params(self._transform_params) + + # add mean value + params = np.concatenate( + (params, + np.mean(params, axis=0).reshape(1, -1), + np.std(params, axis=0).reshape(1, -1) + )) + + cols = self._labels_long[dof] + + else: + params = self._transform_params + cols = ["a%d" % (d + 1) for d in range(0, dof)] + + rows = ["Trafo %d" % (d + 1) for d in range(0, N_trafos)] + rows.append("Mean") + rows.append("Std") + + df = pd.DataFrame(params, rows, cols) + print(df) + + if dir_output is not None: + title = self._replace_string(title) + filename = "%s.csv" % title + ph.create_directory(dir_output) + df.to_csv(os.path.join(dir_output, filename)) + + ## + # Plot figure to show parameter distribution. + # Only works for 3D rigid transforms for now. + # \date 2018-01-25 23:30:45+0000 + # + # \param self The object + # + def show(self, title=None, dir_output=None): + params = self._get_scaled_params(self._transform_params) + + N_trafos, dof = self._transform_params.shape + + fig = plt.figure(title) + fig.clf() + + x = range(1, N_trafos+1) + ax = plt.subplot(2, 1, 1) + for i_param in range(0, 3): + ax.plot( + x, params[:, i_param], + marker=ph.MARKERS[i_param], + color=ph.COLORS_TABLEAU20[i_param*2], + linestyle=":", + label=self._labels_short[dof][i_param], + markerfacecolor="w", + ) + ax.set_xticks(x) + plt.ylabel('Rotation [deg]') + plt.legend(loc="best") + + ax = plt.subplot(2, 1, 2) + for i_param in range(0, 3): + ax.plot( + x, params[:, 3+i_param], + marker=ph.MARKERS[i_param], + color=ph.COLORS_TABLEAU20[i_param*2], + linestyle=":", + label=self._labels_short[dof][3+i_param], + markerfacecolor="w", + ) + ax.set_xticks(x) + plt.xlabel('Slice') + plt.ylabel('Translation [mm]') + plt.legend(loc="best") + plt.suptitle(title) + + try: + # Open windows (and also save them) in full screen + manager = plt.get_current_fig_manager() + manager.full_screen_toggle() + except: + pass + + plt.show(block=False) + + if dir_output is not None: + title = self._replace_string(title) + filename = "%s.pdf" % title + ph.save_fig(fig, dir_output, filename) + + def _get_scaled_params(self, transform_params): + dof = self._transform_params.shape[1] + return self._transform_params * self._scale[dof] + + def _replace_string(self, string): + string = re.sub(" ", "_", string) + string = re.sub(":", "", string) + string = re.sub("/", "_", string) + return string diff --git a/niftymic/validation/residual_evaluator.py b/niftymic/validation/residual_evaluator.py new file mode 100644 index 00000000..9abb1dba --- /dev/null +++ b/niftymic/validation/residual_evaluator.py @@ -0,0 +1,369 @@ +## +# \file residual_evaluator.py +# \brief Class to evaluate computed residuals between a +# simulated/projected and original/acquired slices of stacks +# +# Should help to assess the registration accuracy. +# \author Michael Ebner (michael.ebner.14@ucl.ac.uk) +# \date January 2018 +# + + +# Import libraries +import os +import re +import numpy as np +import SimpleITK as sitk +import matplotlib.pyplot as plt + +from nsol.similarity_measures import SimilarityMeasures as \ + SimilarityMeasures +import pysitk.python_helper as ph +import pysitk.statistics_helper as sh + +import niftymic.reconstruction.linear_operators as lin_op +import niftymic.base.exceptions as exceptions + + +## +# Class to evaluate computed residuals between a simulated/projected and +# original/acquired slices of stacks +# +# \date 2018-01-19 17:24:35+0000 +# +class ResidualEvaluator(object): + + ## + # { constructor_description } + # \date 2018-01-19 17:24:46+0000 + # + # \param self The object + # \param stacks List of Stack objects + # \param reference Reference as Stack object. Used to simulate slices + # at the position of the slices in stacks + # \param use_masks Turn on/off using masks for the residual + # evaluation + # \param measures Similarity measures as given in + # nsol.similarity_measures, list of strings + # + def __init__( + self, + stacks=None, + reference=None, + use_slice_masks=True, + use_reference_mask=True, + measures=["NCC", "NMI", "PSNR", "SSIM", "RMSE"], + verbose=True, + ): + self._stacks = stacks + self._reference = reference + self._measures = measures + self._use_slice_masks = use_slice_masks + self._use_reference_mask = use_reference_mask + self._verbose = verbose + + self._slice_projections = None + self._similarities = None + self._slice_similarities = None + + ## + # Sets the stacks. + # \date 2018-01-19 17:26:04+0000 + # + # \param self The object + # \param stacks List of Stack objects + # + def set_stacks(self, stacks): + self._stacks = stacks + + ## + # Sets the reference from which the slices shall be simulated/projected. + # \date 2018-01-19 17:26:14+0000 + # + # \param self The object + # \param reference The reference + # + def set_reference(self, reference): + self._reference = reference + + def get_measures(self): + return self._measures + + ## + # Gets the slice similarities computed between simulated/projected and + # original/acquired slices. + # \date 2018-01-19 17:26:44+0000 + # + # \param self The object + # + # \return The slice similarities for all slices and measures as + # dictionary. E.g. { + # fetal_brain_1: {'NCC': 1D-array[...], 'NMI': 1D-array[..]}, + # ... + # fetal_brain_N: {'NCC': 1D-array[...], 'NMI': 1D-array[..]} + # } + # + def get_slice_similarities(self): + return self._slice_similarities + + ## + # Gets the slice projections. + # \date 2018-01-19 17:27:41+0000 + # + # \param self The object + # + # \return The slice projections as list of lists. E.g. [ + # [stack1_slice1_sim, stack1_slice2_sim, ...], ... + # [stackN_slice1_sim, stackN_slice2_sim, ...] + # ] + # + def get_slice_projections(self): + return self._slice_projections + + ## + # Calculates the slice simulations/projections from the reference given the + # assumed slice acquisition protocol. + # \date 2018-01-19 17:29:20+0000 + # + # \param self The object + # + # \return The slice projections. + # + def compute_slice_projections(self): + + linear_operators = lin_op.LinearOperators() + self._slice_projections = [None] * len(self._stacks) + + for i_stack, stack in enumerate(self._stacks): + slices = stack.get_slices() + self._slice_projections[i_stack] = [None] * len(slices) + + if self._verbose: + ph.print_info( + "Stack %d/%d: Compute slice projections ... " % ( + i_stack + 1, len(self._stacks)), + newline=False) + + # Compute slice projections based on assumed slice acquisition + # protocol + for i_slice, slice in enumerate(slices): + self._slice_projections[i_stack][i_slice] = linear_operators.A( + self._reference, slice) + + if self._verbose: + print("done") + + ## + # Evaluate slice similarities for all simulated slices of all stacks for + # all similarity measures. + # \date 2018-01-19 17:30:37+0000 + # + # \param self The object + # + def evaluate_slice_similarities(self): + if self._slice_projections is None: + raise exceptions.ObjectNotCreated("compute_slice_projections") + + self._slice_similarities = { + stack.get_filename(): {} for stack in self._stacks + } + + similarity_measures = { + m: SimilarityMeasures.similarity_measures[m] + for m in self._measures + } + + for i_stack, stack in enumerate(self._stacks): + slices = stack.get_slices() + stack_name = stack.get_filename() + self._slice_similarities[stack_name] = { + m: np.zeros(len(slices)) for m in self._measures + } + if self._verbose: + ph.print_info( + "Stack %d/%d: Compute similarity measures ... " % ( + i_stack + 1, len(self._stacks)), + newline=False) + + for i_slice, slice in enumerate(slices): + + slice_nda = np.squeeze(sitk.GetArrayFromImage(slice.sitk)) + slice_proj_nda = np.squeeze(sitk.GetArrayFromImage( + self._slice_projections[i_stack][i_slice].sitk)) + + mask_nda = np.ones_like(slice_nda) + + if self._use_slice_masks: + mask_nda *= np.squeeze( + sitk.GetArrayFromImage(slice.sitk_mask)) + if self._use_reference_mask: + mask_nda *= np.squeeze( + sitk.GetArrayFromImage( + self._slice_projections[i_stack][i_slice].sitk_mask)) + indices = np.where(mask_nda > 0) + + if len(indices[0]) > 0: + for m in self._measures: + try: + self._slice_similarities[stack_name][m][i_slice] = \ + similarity_measures[m]( + slice_nda[indices], slice_proj_nda[indices]) + except ValueError as e: + # Error in case only a few/to less non-zero entries + # exist + if m == "SSIM": + self._slice_similarities[ + stack_name][m][i_slice] = \ + SimilarityMeasures.UNDEF[m] + else: + raise ValueError(e.message) + else: + for m in self._measures: + self._slice_similarities[ + stack_name][m][i_slice] = \ + SimilarityMeasures.UNDEF[m] + if self._verbose: + print("done") + + ## + # Writes the computed slice similarities for all stacks to output directory + # \date 2018-01-19 17:42:27+0000 + # + # \param self The object + # \param directory path to output directory, string + # + def write_slice_similarities(self, directory): + for i_stack, stack in enumerate(self._stacks): + stack_name = stack.get_filename() + path_to_file = os.path.join( + directory, "%s.txt" % stack_name) + + # Write header info + header = "# %s, %s\n" % (stack.get_filename(), ph.get_time_stamp()) + header += "# %s\n" % ("\t").join(self._measures) + ph.write_to_file(path_to_file, header, verbose=self._verbose) + + # Write array information + array = np.zeros( + (stack.get_number_of_slices(), len(self._measures))) + for i_m, m in enumerate(self._measures): + array[:, i_m] = self._slice_similarities[stack_name][m] + ph.write_array_to_file(path_to_file, array, verbose=self._verbose) + + ## + # Reads computed slice similarities for all files in directory. + # \date 2018-01-19 17:42:54+0000 + # + # \param self The object + # \param directory The directory + # \param ext The extent + # \post self._slice_similarities updated + # + def read_slice_similarities(self, directory, ext="txt"): + + if not ph.directory_exists(directory): + raise IOError("Given directory '%s' does not exist" % ( + directory)) + + pattern = "([a-zA-Z0-9_\+\-]+)[.]%s" % ext + p = re.compile(pattern) + + stack_names = [ + p.match(f).group(1) + for f in os.listdir(directory) if p.match(f) + ] + + self._slice_similarities = { + stack_name: {} for stack_name in stack_names + } + + for stack_name in stack_names: + path_to_file = os.path.join(directory, "%s.%s" % (stack_name, ext)) + + # Read computed measures + self._measures = ph.read_file_line_by_line(path_to_file)[1] + self._measures = re.sub("# ", "", self._measures) + self._measures = re.sub("\n", "", self._measures) + self._measures = self._measures.split("\t") + + # Read array + array = np.loadtxt(path_to_file, skiprows=2) + + # Ensure correct shape in case only a single slice available + array = array.reshape(-1, len(self._measures)) + + if array.ndim == 1: + array = array.reshape(len(array), 1) + + for i_m, m in enumerate(self._measures): + self._slice_similarities[stack_name][m] = array[:, i_m] + + ## + # Shows the slice similarities in plots. + # \date 2018-02-09 18:28:45+0000 + # + # \param self The object + # \param directory The directory + # \param title The title + # \param measures The measures + # \param threshold The threshold + # + def show_slice_similarities( + self, + directory=None, + title=None, + measures=["NCC"], + threshold=0.8, + ): + + for i_m, measure in enumerate(measures): + fig = plt.figure(measure) + fig.clf() + if title is not None: + title = "%s: %s" % (title, measure) + else: + title = measure + plt.suptitle(title) + + stack_names = self._slice_similarities.keys() + for i_name, stack_name in enumerate(stack_names): + ax = plt.subplot(np.ceil(len(stack_names) / 2.), 2, i_name + 1) + nda = self._slice_similarities[stack_name][measure] + + nda = np.nan_to_num(nda) + x = np.arange(nda.size) + + # indices_in = np.where(nda >= threshold) + indices_out = np.where(nda < threshold) + + plt.plot(x, nda, + color=ph.COLORS_TABLEAU20[0], + markerfacecolor="w", + marker=ph.MARKERS[0], + linestyle=":", + ) + plt.plot(indices_out[0], nda[indices_out], + color=ph.COLORS_TABLEAU20[6], + # markerfacecolor="w", + marker=ph.MARKERS[0], + linestyle="", + ) + plt.plot([x[0], x[-1]], np.ones(2) * threshold, + color=ph.COLORS_TABLEAU20[6], + linestyle="-.", + ) + + plt.xlabel("Slice") + # plt.ylabel(measure) + plt.title(stack_name) + + ax.set_xticks(x) + # ax.set_xticklabels(x + 1) + ax.set_ylim([0, 1]) + + sh.make_figure_fullscreen() + plt.show(block=False) + + if directory is not None: + filename = "slice_similarities_%s.pdf" % measure + ph.save_fig(fig, directory, filename) diff --git a/niftymic/validation/simulate_stacks_from_reconstruction.py b/niftymic/validation/simulate_stacks_from_reconstruction.py index a4f4dd57..32d7bf0d 100644 --- a/niftymic/validation/simulate_stacks_from_reconstruction.py +++ b/niftymic/validation/simulate_stacks_from_reconstruction.py @@ -41,7 +41,9 @@ def main(): "correct, the resulting stack of such obtained projected slices, " "corresponds to the originally acquired (motion corrupted) data.", ) - input_parser.add_dir_input(required=True) + input_parser.add_filenames(required=True) + input_parser.add_filenames_masks() + input_parser.add_dir_input_mc(required=True) input_parser.add_reconstruction(required=True) input_parser.add_dir_output(required=True) input_parser.add_suffix_mask(default="_mask") @@ -63,6 +65,7 @@ def main(): help="Choose the interpolator type to propagate the reconstruction " "mask (%s)." % (INTERPOLATOR_TYPES), default="NearestNeighbor") + input_parser.add_log_config(default=1) input_parser.add_verbose(default=0) args = input_parser.parse_args() @@ -73,10 +76,16 @@ def main(): "Unknown interpolator provided. Possible choices are %s" % ( INTERPOLATOR_TYPES)) + if args.log_config: + input_parser.log_config(os.path.abspath(__file__)) + # Read motion corrected data - data_reader = dr.ImageSlicesDirectoryReader( - path_to_directory=args.dir_input, - suffix_mask=args.suffix_mask) + data_reader = dr.MultipleImagesReader( + file_paths=args.filenames, + file_paths_masks=args.filenames_masks, + suffix_mask=args.suffix_mask, + dir_motion_correction=args.dir_input_mc, + ) data_reader.read_data() stacks = data_reader.get_data() @@ -89,22 +98,26 @@ def main(): # initialize image data array(s) nda = np.zeros_like(sitk.GetArrayFromImage(stack.sitk)) + nda[:] = np.nan if args.reconstruction_mask: nda_mask = np.zeros_like(sitk.GetArrayFromImage(stack.sitk_mask)) - # Simulate slices at estimated positions within reconstructed volume - simulated_slices = [ - linear_operators.A( - reconstruction, s, interpolator_mask=args.interpolator_mask) - for s in stack.get_slices() - ] + slices = stack.get_slices() + kept_indices = [s.get_slice_number() for s in slices] # Fill stack information "as if slice was acquired consecutively" # Therefore, simulated stack slices correspond to acquired slices # (in case motion correction was correct) - for j, simulated_slice in enumerate(simulated_slices): - nda[j, :, :] = sitk.GetArrayFromImage(simulated_slice.sitk) + for j in range(nda.shape[0]): + if j in kept_indices: + index = kept_indices.index(j) + simulated_slice = linear_operators.A( + reconstruction, + slices[index], + interpolator_mask=args.interpolator_mask + ) + nda[j, :, :] = sitk.GetArrayFromImage(simulated_slice.sitk) if args.reconstruction_mask: nda_mask[j, :, :] = sitk.GetArrayFromImage( @@ -129,8 +142,7 @@ def main(): if args.verbose: sitkh.show_stacks([ stack, simulated_stack], - segmentation=simulated_stack - if args.reconstruction_mask else None) + segmentation=stack) simulated_stack.write( args.dir_output, diff --git a/niftymic_correct_intensities.py b/niftymic_correct_intensities.py deleted file mode 100644 index acef40ec..00000000 --- a/niftymic_correct_intensities.py +++ /dev/null @@ -1,7 +0,0 @@ -# -*- coding: utf-8 -*- -import sys - -from niftymic.application.correct_intensities import main - -if __name__ == "__main__": - sys.exit(main()) diff --git a/requirements.txt b/requirements.txt index d4b94de0..730339c2 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,3 @@ --e git+https://github.com/gift-surg/PySiTK.git@v0.1.1#egg=PySiTK-0.1.1 --e git+https://github.com/gift-surg/NSoL.git@v0.1.1#egg=NSoL-0.1.1 --e git+https://github.com/gift-surg/SimpleReg.git@v0.1.1#egg=SimpleReg-0.1.1 +-e git+https://github.com/gift-surg/PySiTK.git@v0.2#egg=PySiTK-0.2 +-e git+https://github.com/gift-surg/NSoL.git@v0.1.4#egg=NSoL-0.1.4 +-e git+https://github.com/gift-surg/SimpleReg.git@v0.2#egg=SimpleReg-0.2 diff --git a/setup.py b/setup.py index 68d80e1e..84694c3c 100644 --- a/setup.py +++ b/setup.py @@ -61,7 +61,7 @@ def run(self): setup(name='NiftyMIC', - version='0.2.2', + version='0.3', description=description, long_description=long_description, url='https://github.com/gift-surg/NiftyMIC', @@ -70,14 +70,17 @@ def run(self): license='BSD-3-Clause', packages=['niftymic'], install_requires=[ - 'pysitk>=0.1', - 'nsol>=0.1', - 'simplereg>=0.1', + 'pysitk>=0.2', + 'nsol>=0.1.4', + 'simplereg>=0.2', 'scikit_image>=0.12.3', 'scipy>=0.19.1', 'natsort>=5.0.3', 'numpy>=1.13.1', 'SimpleITK>=1.0.1', + 'pymc3>=3.3', + 'theano>=1.0.1', + "six>=1.10.0", ], zip_safe=False, keywords='development numericalsolver convexoptimisation', @@ -102,7 +105,6 @@ def run(self): entry_points={ 'console_scripts': [ 'niftymic_correct_bias_field = niftymic.application.correct_bias_field:main', - 'niftymic_correct_intensities = niftymic.application.correct_intensities:main', 'niftymic_reconstruct_volume = niftymic.application.reconstruct_volume:main', 'niftymic_reconstruct_volume_from_slices = niftymic.application.reconstruct_volume_from_slices:main', 'niftymic_register_image = niftymic.application.register_image:main', diff --git a/tests/brain_stripping_test.py b/tests/brain_stripping_test.py index 7e367597..d1047e2f 100644 --- a/tests/brain_stripping_test.py +++ b/tests/brain_stripping_test.py @@ -5,12 +5,14 @@ # \date December 2015 -# Import libraries +import os import unittest +import numpy as np +import SimpleITK as sitk -import niftymic.utilities.brain_stripping as bs -# Import modules from src-folder import pysitk.simple_itk_helper as sitkh + +import niftymic.utilities.brain_stripping as bs from niftymic.definitions import DIR_TEST @@ -22,13 +24,15 @@ class BrainStrippingTest(unittest.TestCase): accuracy = 7 def setUp(self): - pass + self.precision = 7 + self.dir_data = os.path.join( + DIR_TEST, "case-studies", "fetal-brain", "input-data") + self.filename = "axial" def test_01_input_output(self): - filename = "stack0" brain_stripping = bs.BrainStripping.from_filename( - self.dir_test_data, filename) + self.dir_data, self.filename) brain_stripping.compute_brain_image(0) brain_stripping.compute_brain_mask(0) brain_stripping.compute_skull_image(0) @@ -53,10 +57,11 @@ def test_01_input_output(self): str(ve.exception)) def test_02_brain_mask(self): - filename = "stack0" + path_to_reference = os.path.join( + DIR_TEST, "case-studies", "fetal-brain", "brain_stripping", "axial_seg.nii.gz") brain_stripping = bs.BrainStripping.from_filename( - self.dir_test_data, filename) + self.dir_data, self.filename) brain_stripping.compute_brain_image(0) brain_stripping.compute_brain_mask(1) brain_stripping.compute_skull_image(0) @@ -64,5 +69,10 @@ def test_02_brain_mask(self): brain_stripping.run() original_sitk = brain_stripping.get_input_image_sitk() - brain_mask_sitk = brain_stripping.get_brain_mask_sitk() - sitkh.show_sitk_image([original_sitk], segmentation=brain_mask_sitk) + res_sitk = brain_stripping.get_brain_mask_sitk() + + ref_sitk = sitkh.read_nifti_image_sitk(path_to_reference) + + diff_sitk = res_sitk - ref_sitk + error = np.linalg.norm(sitk.GetArrayFromImage(diff_sitk)) + self.assertAlmostEqual(error, 0, places=self.precision) diff --git a/tests/case_study_fetal_brain_test.py b/tests/case_study_fetal_brain_test.py index ab7c3ecf..03dc5000 100644 --- a/tests/case_study_fetal_brain_test.py +++ b/tests/case_study_fetal_brain_test.py @@ -13,8 +13,9 @@ import SimpleITK as sitk import pysitk.python_helper as ph +import pysitk.simple_itk_helper as sitkh -from niftymic.definitions import DIR_TMP, DIR_TEST +from niftymic.definitions import DIR_TMP, DIR_TEST, REGEX_FILENAMES, DIR_TEMPLATES class CaseStudyFetalBrainTest(unittest.TestCase): @@ -22,65 +23,152 @@ class CaseStudyFetalBrainTest(unittest.TestCase): def setUp(self): self.precision = 7 self.dir_data = os.path.join(DIR_TEST, "case-studies", "fetal-brain") + self.filenames = [ + os.path.join(self.dir_data, + "input-data", + "%s.nii.gz" % f) + for f in ["axial", "coronal", "sagittal"]] self.dir_output = os.path.join(DIR_TMP, "case-studies", "fetal-brain") + self.suffix_mask = "_mask" - def test_reconstruct_volume_from_slices(self): - dir_root = os.path.join( - self.dir_data, "reconstruct_volume_from_slices") - dir_input = os.path.join(dir_root, "input-data") - dir_reference = os.path.join(dir_root, "result-comparison") - filename_reference = "SRR_stacks3_TK1_lsmr_alpha0p02_itermax10.nii.gz" + def test_reconstruct_volume(self): + dir_reference = os.path.join(self.dir_data, "reconstruct_volume") + dir_reference_mc = os.path.join(dir_reference, "motion_correction") + filename_reference = "SRR_stacks3_TK1_lsmr_alpha0p02_itermax5.nii.gz" path_to_reference = os.path.join(dir_reference, filename_reference) + path_to_reference_mask = ph.append_to_filename( + os.path.join(dir_reference, filename_reference), self.suffix_mask) + + two_step_cycles = 1 + iter_max = 5 + threshold = 0.8 + alpha = 0.02 cmd_args = [] - cmd_args.append("--dir-input %s" % dir_input) + cmd_args.append("--filenames %s" % " ".join(self.filenames)) cmd_args.append("--dir-output %s" % self.dir_output) - cmd_args.append("--reconstruction-space %s" % path_to_reference) + cmd_args.append("--suffix-mask %s" % self.suffix_mask) + cmd_args.append("--two-step-cycles %s" % two_step_cycles) + cmd_args.append("--iter-max %d" % iter_max) + cmd_args.append("--threshold-first %f" % threshold) + cmd_args.append("--threshold %f" % threshold) + cmd_args.append("--alpha %f" % alpha) - cmd = "niftymic_reconstruct_volume_from_slices %s" % ( + cmd = "niftymic_reconstruct_volume %s" % ( " ").join(cmd_args) self.assertEqual(ph.execute_command(cmd), 0) - # Check whether identical reconstruction has been created - path_to_reconstruction = os.path.join( - self.dir_output, filename_reference) - reconstruction_sitk = sitk.ReadImage(path_to_reconstruction) - reference_sitk = sitk.ReadImage(path_to_reference) + # Check SRR volume + res_sitk = sitkh.read_nifti_image_sitk( + os.path.join(self.dir_output, filename_reference)) + ref_sitk = sitkh.read_nifti_image_sitk(path_to_reference) - difference_sitk = reconstruction_sitk - reference_sitk - error = np.linalg.norm(sitk.GetArrayFromImage(difference_sitk)) + diff_sitk = res_sitk - ref_sitk + error = np.linalg.norm(sitk.GetArrayFromImage(diff_sitk)) + self.assertAlmostEqual(error, 0, places=self.precision) + + # Check SRR mask volume + res_sitk = sitkh.read_nifti_image_sitk( + ph.append_to_filename( + os.path.join(self.dir_output, filename_reference), + self.suffix_mask)) + ref_sitk = sitkh.read_nifti_image_sitk(path_to_reference_mask) + diff_sitk = res_sitk - ref_sitk + error = np.linalg.norm(sitk.GetArrayFromImage(diff_sitk)) self.assertAlmostEqual(error, 0, places=self.precision) - def test_reconstruct_volume(self): - dir_root = os.path.join(self.dir_data, "reconstruct_volume") - dir_input = os.path.join(dir_root, "input-data") - dir_reference = os.path.join(dir_root, "result-comparison") + # Check transforms + pattern = REGEX_FILENAMES + "[.]tfm" + p = re.compile(pattern) + dir_res_mc = os.path.join(self.dir_output, "motion_correction") + trafos_res = sorted( + [os.path.join(dir_res_mc, t) + for t in os.listdir(dir_res_mc) if p.match(t)]) + trafos_ref = sorted( + [os.path.join(dir_reference_mc, t) + for t in os.listdir(dir_reference_mc) if p.match(t)]) + self.assertEqual(len(trafos_res), len(trafos_ref)) + for i in range(len(trafos_ref)): + nda_res = sitkh.read_transform_sitk(trafos_res[i]).GetParameters() + nda_ref = sitkh.read_transform_sitk(trafos_ref[i]).GetParameters() + nda_diff = np.linalg.norm(np.array(nda_res) - nda_ref) + self.assertAlmostEqual(nda_diff, 0, places=self.precision) + + def test_reconstruct_volume_from_slices(self): + dir_reference = os.path.join( + self.dir_data, "reconstruct_volume_from_slices") + dir_input_mc = os.path.join( + self.dir_data, "reconstruct_volume", "motion_correction") filename_reference = "SRR_stacks3_TK1_lsmr_alpha0p02_itermax5.nii.gz" path_to_reference = os.path.join(dir_reference, filename_reference) - two_step_cycles = 1 iter_max = 5 + alpha = 0.02 cmd_args = [] - cmd_args.append("--dir-input %s" % dir_input) + cmd_args.append("--filenames %s" % " ".join(self.filenames)) + cmd_args.append("--dir-input-mc %s" % dir_input_mc) cmd_args.append("--dir-output %s" % self.dir_output) - cmd_args.append("--two-step-cycles %s" % two_step_cycles) - cmd_args.append("--iter-max %s" % iter_max) + cmd_args.append("--iter-max %d" % iter_max) + cmd_args.append("--alpha %f" % alpha) + cmd_args.append("--reconstruction-space %s" % path_to_reference) - cmd = "niftymic_reconstruct_volume %s" % ( + cmd = "niftymic_reconstruct_volume_from_slices %s" % ( " ").join(cmd_args) self.assertEqual(ph.execute_command(cmd), 0) # Check whether identical reconstruction has been created path_to_reconstruction = os.path.join( self.dir_output, filename_reference) - reconstruction_sitk = sitk.ReadImage(path_to_reconstruction) - reference_sitk = sitk.ReadImage(path_to_reference) + reconstruction_sitk = sitkh.read_nifti_image_sitk( + path_to_reconstruction) + reference_sitk = sitkh.read_nifti_image_sitk(path_to_reference) difference_sitk = reconstruction_sitk - reference_sitk error = np.linalg.norm(sitk.GetArrayFromImage(difference_sitk)) self.assertAlmostEqual(error, 0, places=self.precision) - # Obtained reconstructions could be tested too + def test_register_image(self): + filename_reference = "SRR_stacks3_TK1_lsmr_alpha0p02_itermax5.nii.gz" + path_to_recon = os.path.join( + self.dir_data, "reconstruct_volume", filename_reference) + dir_input_mc = os.path.join( + self.dir_data, "reconstruct_volume", "motion_correction") + gestational_age = 28 + + path_to_transform_ref = os.path.join( + self.dir_data, "register_image", "registration_transform_sitk.txt") + path_to_transform_res = os.path.join( + self.dir_output, "registration_transform_sitk.txt") + + template = os.path.join( + DIR_TEMPLATES, + "STA%d.nii.gz" % gestational_age) + template_mask = os.path.join( + DIR_TEMPLATES, + "STA%d_mask.nii.gz" % gestational_age) + + cmd_args = ["niftymic_register_image"] + cmd_args.append("--fixed %s" % template) + cmd_args.append("--moving %s" % path_to_recon) + cmd_args.append("--fixed-mask %s" % template_mask) + cmd_args.append("--moving-mask %s" % + ph.append_to_filename(path_to_recon, self.suffix_mask)) + cmd_args.append("--dir-input-mc %s" % dir_input_mc) + cmd_args.append("--dir-output %s" % self.dir_output) + cmd_args.append("--use-flirt 1") + cmd_args.append("--use-regaladin 1") + cmd_args.append("--test-ap-flip 1") + self.assertEqual(ph.execute_command(" ".join(cmd_args)), 0) + + res_sitk = sitkh.read_transform_sitk(path_to_transform_res) + ref_sitk = sitkh.read_transform_sitk(path_to_transform_ref) + + res_nda = res_sitk.GetParameters() + ref_nda = ref_sitk.GetParameters() + diff_nda = np.array(res_nda) - ref_nda + + self.assertAlmostEqual( + np.linalg.norm(diff_nda), 0, places=self.precision) diff --git a/tests/data_reader_test.py b/tests/data_reader_test.py new file mode 100644 index 00000000..1ca3163f --- /dev/null +++ b/tests/data_reader_test.py @@ -0,0 +1,68 @@ +## +# \file data_reader_test.py +# \brief Unit tests for data reader +# +# \author Michael Ebner (michael.ebner.14@ucl.ac.uk) +# \date January 2018 + + +import os +import unittest +import numpy as np +import re +import SimpleITK as sitk + +import pysitk.python_helper as ph +import niftymic.base.data_reader as dr + +from niftymic.definitions import DIR_TMP, DIR_TEST + + +class DataReaderTest(unittest.TestCase): + + def setUp(self): + self.precision = 7 + self.dir_data = os.path.join(DIR_TEST, "case-studies", "fetal-brain") + self.filenames = [ + os.path.join(self.dir_data, + "input-data", + "%s.nii.gz" % f) + for f in ["axial", "coronal", "sagittal"]] + self.dir_output = os.path.join(DIR_TMP, "case-studies", "fetal-brain") + self.suffix_mask = "_mask" + + ## + # Check that the same number of stacks (and slices therein) are read + # \date 2018-01-31 23:03:52+0000 + # + # \param self The object + # + def test_read_transformations(self): + + directory_motion_correction = os.path.join( + DIR_TEST, + "case-studies", + "fetal-brain", + "reconstruct_volume", + "motion_correction", + ) + + data_reader = dr.MultipleImagesReader( + file_paths=self.filenames, + dir_motion_correction=directory_motion_correction) + data_reader.read_data() + stacks = data_reader.get_data() + + data_reader = dr.SliceTransformationDirectoryReader( + directory_motion_correction) + data_reader.read_data() + transformations_dic = data_reader.get_data() + + self.assertEqual(len(stacks) - len(transformations_dic.keys()), 0) + + for stack in stacks: + N_slices = stack.get_number_of_slices() + N_slices2 = len(transformations_dic[stack.get_filename()].keys()) + self.assertEqual(N_slices - N_slices2, 0) + + diff --git a/tests/image_similarity_evaluator_test.py b/tests/image_similarity_evaluator_test.py new file mode 100644 index 00000000..0c4ebc82 --- /dev/null +++ b/tests/image_similarity_evaluator_test.py @@ -0,0 +1,72 @@ +## +# \file image_similarity_evaluator_test.py +# \brief Test ImageSimilarityEvaluator class +# +# \author Michael Ebner (michael.ebner.14@ucl.ac.uk) +# \date February 2018 + + +import os +import unittest +import numpy as np +import re +import SimpleITK as sitk + +import pysitk.python_helper as ph + +import niftymic.base.stack as st +import niftymic.validation.image_similarity_evaluator as ise +import niftymic.base.exceptions as exceptions +from niftymic.definitions import DIR_TMP, DIR_TEST + + +class ImageSimilarityEvaluatorTest(unittest.TestCase): + + def setUp(self): + self.precision = 7 + + def test_compute_write_read_similarities(self): + + paths_to_stacks = [ + os.path.join( + DIR_TEST, "fetal_brain_%d.nii.gz" % d) for d in range(0, 3) + ] + path_to_reference = os.path.join( + DIR_TEST, "FetalBrain_reconstruction_3stacks_myAlg.nii.gz") + + reference = st.Stack.from_filename( + path_to_reference, extract_slices=False) + + stacks = [ + st.Stack.from_filename(p, ph.append_to_filename(p, "_mask")) + for p in paths_to_stacks + ] + stacks = [s.get_resampled_stack(reference.sitk) for s in stacks] + + residual_evaluator = ise.ImageSimilarityEvaluator(stacks, reference) + residual_evaluator.compute_similarities() + residual_evaluator.write_similarities(DIR_TMP) + similarities = residual_evaluator.get_similarities() + + similarities1 = ise.ImageSimilarityEvaluator() + similarities1.read_similarities(DIR_TMP) + similarities1 = similarities1.get_similarities() + + for m in residual_evaluator.get_measures(): + rho_res = similarities[m] + rho_res1 = similarities1[m] + error = np.linalg.norm(rho_res - rho_res1) + self.assertAlmostEqual(error, 0, places=self.precision) + + def test_results_not_created(self): + residual_evaluator = ise.ImageSimilarityEvaluator() + + # Directory does not exist + self.assertRaises( + IOError, lambda: + residual_evaluator.read_similarities( + os.path.join(DIR_TMP, "whatevertestasdfsfasdasf"))) + + # Directory does exist but does not contain similarity result files + self.assertRaises(IOError, lambda: + residual_evaluator.read_similarities(DIR_TEST)) diff --git a/tests/intensity_correction_test.py b/tests/intensity_correction_test.py index 9aafd5ef..9edbd243 100644 --- a/tests/intensity_correction_test.py +++ b/tests/intensity_correction_test.py @@ -6,15 +6,13 @@ # \author Michael Ebner (michael.ebner.14@ucl.ac.uk) # \date October 2016 - +import os import unittest - -# Import libraries -import SimpleITK as sitk import numpy as np -from scipy.ndimage import imread +import SimpleITK as sitk + +import pysitk.python_helper as ph -# Import modules import niftymic.base.stack as st import niftymic.utilities.intensity_correction as ic from niftymic.definitions import DIR_TEST @@ -37,7 +35,9 @@ def test_linear_intensity_correction(self): shape_z = 15 # Original stack - nda_2D = imread(self.dir_test_data + "2D_Lena_256.png", flatten=True) + nda_2D = ph.read_image( + os.path.join(self.dir_test_data, "2D_Lena_256.png")) + nda_3D = np.tile(nda_2D, (shape_z, 1, 1)).astype('double') stack_sitk = sitk.GetImageFromArray(nda_3D) stack = st.Stack.from_sitk_image(stack_sitk, "Lena") @@ -45,7 +45,7 @@ def test_linear_intensity_correction(self): # 1) Create linearly corrupted intensity stack nda_3D_corruped = np.zeros_like(nda_3D) for i in range(0, shape_z): - nda_3D_corruped[i, :, :] = nda_3D[i, :, :]/(i+1.) + nda_3D_corruped[i, :, :] = nda_3D[i, :, :] / (i + 1.) stack_corrupted_sitk = sitk.GetImageFromArray(nda_3D_corruped) stack_corrupted = st.Stack.from_sitk_image( stack_corrupted_sitk, "stack_corrupted") @@ -55,10 +55,13 @@ def test_linear_intensity_correction(self): # Ground truth-parameter: ic_values = np.zeros((shape_z, 1)) for i in range(0, shape_z): - ic_values[i, :] = (i+1.) + ic_values[i, :] = (i + 1.) intensity_correction = ic.IntensityCorrection( - stack=stack_corrupted, reference=stack, use_individual_slice_correction=True, use_verbose=self.use_verbose) + stack=stack_corrupted, + reference=stack, + use_individual_slice_correction=True, + use_verbose=self.use_verbose) intensity_correction.run_linear_intensity_correction() ic_values_est = intensity_correction.get_intensity_correction_coefficients() @@ -72,7 +75,8 @@ def test_affine_intensity_correction(self): shape_z = 15 # Original stack - nda_2D = imread(self.dir_test_data + "2D_Lena_256.png", flatten=True) + nda_2D = ph.read_image( + os.path.join(self.dir_test_data, "2D_Lena_256.png")) nda_3D = np.tile(nda_2D, (shape_z, 1, 1)).astype('double') stack_sitk = sitk.GetImageFromArray(nda_3D) stack = st.Stack.from_sitk_image(stack_sitk, "Lena") @@ -80,7 +84,7 @@ def test_affine_intensity_correction(self): # 1) Create linearly corrupted intensity stack nda_3D_corruped = np.zeros_like(nda_3D) for i in range(0, shape_z): - nda_3D_corruped[i, :, :] = (nda_3D[i, :, :]-10*i)/(i+1.) + nda_3D_corruped[i, :, :] = (nda_3D[i, :, :] - 10 * i) / (i + 1.) stack_corrupted_sitk = sitk.GetImageFromArray(nda_3D_corruped) stack_corrupted = st.Stack.from_sitk_image( stack_corrupted_sitk, "stack_corrupted") @@ -90,13 +94,16 @@ def test_affine_intensity_correction(self): # Ground truth-parameter: ic_values = np.zeros((shape_z, 2)) for i in range(0, shape_z): - ic_values[i, :] = (i+1, 10*i) + ic_values[i, :] = (i + 1, 10 * i) intensity_correction = ic.IntensityCorrection( - stack=stack_corrupted, reference=stack, use_individual_slice_correction=True, use_verbose=self.use_verbose) + stack=stack_corrupted, + reference=stack, + use_individual_slice_correction=True, + use_verbose=self.use_verbose) intensity_correction.run_affine_intensity_correction() ic_values_est = intensity_correction.get_intensity_correction_coefficients() - + nda_diff = ic_values - ic_values_est self.assertEqual(np.round( np.linalg.norm(nda_diff), decimals=self.accuracy), 0) diff --git a/tests/linear_operators_test.py b/tests/linear_operators_test.py index b6eaf39d..f83048a6 100644 --- a/tests/linear_operators_test.py +++ b/tests/linear_operators_test.py @@ -22,31 +22,23 @@ simulate_stacks_from_reconstruction from niftymic.definitions import DIR_TMP, DIR_TEST + class LinearOperatorsTest(unittest.TestCase): def setUp(self): self.precision = 7 self.dir_output = os.path.join(DIR_TMP, "reconstruction") - self.dir_data = os.path.join(DIR_TEST, "reconstruction") - self.filenames = [ - "IC_N4ITK_HASTE_exam_3.5mm_800ms_3", - # "IC_N4ITK_HASTE_exam_3.5mm_800ms_4", - # "IC_N4ITK_HASTE_exam_3.5mm_800ms_5", - # "IC_N4ITK_HASTE_exam_3.5mm_800ms_6", - # "IC_N4ITK_HASTE_exam_3.5mm_800ms_7", - ] - self.filename_recon = "SRR_stacks5_alpha0p01" - - self.suffix_mask = "_brain" - self.paths_to_filenames = [ - os.path.join(self.dir_data, "motion_correction", f + ".nii.gz") - for f in self.filenames] - + self.dir_data = os.path.join(DIR_TEST, "case-studies", "fetal-brain") + self.filename = "axial" + self.suffix_mask = "_mask" + self.path_to_file = os.path.join( + self.dir_data, "input-data", "%s.nii.gz" % self.filename) + self.filename_recon = "SRR_stacks3_TK1_lsmr_alpha0p02_itermax5.nii.gz" self.path_to_recon = os.path.join( - self.dir_data, self.filename_recon + ".nii.gz") - self.path_to_recon_mask = os.path.join( - self.dir_data, self.filename_recon + self.suffix_mask + ".nii.gz") + self.dir_data, "recon_projections", self.filename_recon) + self.path_to_recon_mask = ph.append_to_filename( + self.path_to_recon, self.suffix_mask) ## # Test forward simulation of stack and associated propagation of @@ -55,12 +47,7 @@ def setUp(self): # def test_forward_operator_stack(self): - data_reader = dr.MultipleImagesReader( - self.paths_to_filenames, suffix_mask=self.suffix_mask) - data_reader.read_data() - stacks = data_reader.get_data() - stack = stacks[0] - + stack = st.Stack.from_filename(self.path_to_file) reconstruction = st.Stack.from_filename( self.path_to_recon, self.path_to_recon_mask) @@ -69,22 +56,22 @@ def test_forward_operator_stack(self): reconstruction, stack, interpolator_mask="Linear") simulated_stack.set_filename(stack.get_filename() + "_sim") - # sitkh.show_stacks([stack, simulated_stack]) - # simulated_stack.show(1) - # reconstruction.show(1) - # stack.show(1) + # sitkh.show_stacks( + # [stack, simulated_stack], segmentation=simulated_stack) + + filename_reference = os.path.join( + self.dir_data, + "recon_projections", + "stack", + "%s_sim.nii.gz" % self.filename) + filename_reference_mask = os.path.join( + self.dir_data, + "recon_projections", + "stack", + "%s_sim%s.nii.gz" % (self.filename, self.suffix_mask)) - filename_reference = "IC_N4ITK_HASTE_exam_3.5mm_800ms_3_simulated" reference_simulated_stack = st.Stack.from_filename( - os.path.join( - self.dir_data, - "result-comparison", - filename_reference + ".nii.gz"), - os.path.join( - self.dir_data, - "result-comparison", - filename_reference + self.suffix_mask + ".nii.gz") - ) + filename_reference, filename_reference_mask) # Error simulated stack difference_sitk = simulated_stack.sitk - \ @@ -107,15 +94,49 @@ def test_forward_operator_stack(self): def test_simulate_stacks_from_slices(self): cmd_args = [] - cmd_args.append("--dir-input %s" % - os.path.join(self.dir_data, "motion_correction")) + cmd_args.append("--filenames %s" % self.path_to_file) + cmd_args.append("--dir-input-mc %s" % + os.path.join( + self.dir_data, + "recon_projections", + "motion_correction")) cmd_args.append("--reconstruction %s" % self.path_to_recon) cmd_args.append("--reconstruction-mask %s" % self.path_to_recon_mask) cmd_args.append("--copy-data 1") - cmd_args.append("--suffix-mask _brain") - # cmd_args.append("--verbose 1") + cmd_args.append("--suffix-mask %s" % self.suffix_mask) cmd_args.append("--dir-output %s" % self.dir_output) exe = os.path.abspath(simulate_stacks_from_reconstruction.__file__) cmd = "python %s %s" % (exe, (" ").join(cmd_args)) self.assertEqual(ph.execute_command(cmd), 0) + + path_orig = os.path.join(self.dir_output, "%s.nii.gz" % self.filename) + path_sim = os.path.join( + self.dir_output, "Simulated_%s.nii.gz" % self.filename) + path_sim_mask = os.path.join( + self.dir_output, "Simulated_%s%s.nii.gz" % (self.filename, self.suffix_mask)) + + path_orig_ref = os.path.join(self.dir_data, + "recon_projections", + "slices", + "%s.nii.gz" % self.filename) + path_sim_ref = os.path.join(self.dir_data, + "recon_projections", + "slices", + "Simulated_%s.nii.gz" % self.filename) + path_sim_mask_ref = os.path.join(self.dir_data, + "recon_projections", + "slices", + "Simulated_%s%s.nii.gz" % ( + self.filename, self.suffix_mask)) + + for res, ref in zip( + [path_orig, path_sim, path_sim_mask], + [path_orig_ref, path_sim_ref, path_sim_mask_ref]): + res_sitk = sitk.ReadImage(res) + ref_sitk = sitk.ReadImage(ref) + + nda_diff = np.nan_to_num( + sitk.GetArrayFromImage(res_sitk - ref_sitk)) + self.assertAlmostEqual(np.linalg.norm( + nda_diff), 0, places=self.precision) diff --git a/tests/niftyreg_test.py b/tests/niftyreg_test.py index 499bc3d7..99146b67 100644 --- a/tests/niftyreg_test.py +++ b/tests/niftyreg_test.py @@ -12,7 +12,7 @@ import sys import os -# Import modules +import pysitk.python_helper as ph import pysitk.simple_itk_helper as sitkh import niftymic.registration.niftyreg as nreg @@ -37,9 +37,10 @@ def test_affine_transform_reg_aladin(self): filename_fixed = "stack1_rotated_angle_z_is_pi_over_10.nii.gz" filename_moving = "FetalBrain_reconstruction_3stacks_myAlg.nii.gz" + diff_ref = os.path.join( + DIR_TEST, "stack1_rotated_angle_z_is_pi_over_10_nreg_diff.nii.gz") moving = st.Stack.from_filename( os.path.join(self.dir_test_data, filename_moving), - # os.path.join(self.dir_test_data, filename_moving + "_mask.nii.gz") ) fixed = st.Stack.from_filename( os.path.join(self.dir_test_data, filename_fixed) @@ -57,33 +58,18 @@ def test_affine_transform_reg_aladin(self): # Get associated results affine_transform_sitk = nifty_reg.get_registration_transform_sitk() - moving_warped = nifty_reg.get_transformed_fixed() + moving_warped = nifty_reg.get_warped_moving() # Get SimpleITK result with "similar" interpolator (NiftyReg does not # state what interpolator is used but it seems to be BSpline) moving_warped_sitk = sitk.Resample( moving.sitk, fixed.sitk, affine_transform_sitk, sitk.sitkBSpline, 0.0, moving.sitk.GetPixelIDValue()) - # Check alignment of images - nda_NiftyReg = sitk.GetArrayFromImage(moving_warped.sitk) - nda_SimpleITK = sitk.GetArrayFromImage(moving_warped_sitk) - diff = nda_NiftyReg - nda_SimpleITK - abs_diff = abs(diff) - - try: - self.assertEqual(np.round( - np.linalg.norm(diff), decimals=self.accuracy), 0) - - except Exception as e: - print("FAIL: " + self.id() + " failed given norm of difference = %.2e > 1e-%s" % - (np.linalg.norm(diff), self.accuracy)) - print( - "\tCheck statistics of difference: (Maximum absolute difference per voxel might be acceptable)") - print("\tMaximum absolute difference per voxel: %s" % abs_diff.max()) - print("\tMean absolute difference per voxel: %s" % abs_diff.mean()) - print("\tMinimum absolute difference per voxel: %s" % abs_diff.min()) - - # Show results (difficult to compare directly given the different interpolators of NiftyReg and SimpleITK) - # sitkh.show_sitk_image(moving_warped.sitk, overlay=fixed.sitk, title="warpedMoving_fixed") - # sitkh.show_sitk_image(moving_warped.sitk, overlay=moving_warped_sitk, title="warpedMovingNiftyReg_warpedMovingSimpleITK") - # sitkh.show_sitk_image(moving_warped.sitk-moving_warped_sitk, title="difference_NiftyReg_SimpleITK") + diff_res_sitk = moving_warped.sitk - moving_warped_sitk + sitkh.write_nifti_image_sitk(diff_res_sitk, diff_ref) + diff_ref_sitk = sitk.ReadImage(diff_ref) + + res_diff_nda = sitk.GetArrayFromImage(diff_res_sitk - diff_ref_sitk) + + self.assertAlmostEqual( + np.linalg.norm(res_diff_nda), 0, places=self.accuracy) diff --git a/tests/parameter_normalization_test.py b/tests/parameter_normalization_test.py index fd0dd1cc..398d4859 100644 --- a/tests/parameter_normalization_test.py +++ b/tests/parameter_normalization_test.py @@ -6,13 +6,11 @@ # \date Nov 2016 -# Import libraries -import SimpleITK as sitk -import itk -import numpy as np -import unittest +import os import sys -from scipy.ndimage import imread +import unittest +import numpy as np +import SimpleITK as sitk # Import modules import pysitk.simple_itk_helper as sitkh @@ -43,9 +41,9 @@ def test_parameter_normalization(self): filename_stack_corrupted = "FetalBrain_reconstruction_3stacks_myAlg_corrupted_inplane" stack_sitk = sitk.ReadImage( - self.dir_test_data + filename_stack + ".nii.gz") + os.path.join(self.dir_test_data, filename_stack + ".nii.gz")) stack_corrupted_sitk = sitk.ReadImage( - self.dir_test_data + filename_stack_corrupted + ".nii.gz") + os.path.join(self.dir_test_data, filename_stack_corrupted + ".nii.gz")) stack_corrupted = st.Stack.from_sitk_image( stack_corrupted_sitk, "stack_corrupted") @@ -108,7 +106,7 @@ def test_parameter_normalization(self): # Check mean self.assertEqual(np.round( - abs(mean-coefficients[0, i]), decimals=self.accuracy), 0) + abs(mean - coefficients[0, i]), decimals=self.accuracy), 0) # Check standard deviation if abs(std) > 1e-8: @@ -117,4 +115,4 @@ def test_parameter_normalization(self): # Check parameter values self.assertEqual(np.round( - np.linalg.norm(parameters_tmp-parameters), decimals=self.accuracy), 0) + np.linalg.norm(parameters_tmp - parameters), decimals=self.accuracy), 0) diff --git a/tests/residual_evaluator_test.py b/tests/residual_evaluator_test.py new file mode 100644 index 00000000..b329af71 --- /dev/null +++ b/tests/residual_evaluator_test.py @@ -0,0 +1,79 @@ +## +# \file residual_evaluator_test.py +# \brief Test ResidualEvaluator class +# +# \author Michael Ebner (michael.ebner.14@ucl.ac.uk) +# \date November 2017 + + +import os +import unittest +import numpy as np +import re +import SimpleITK as sitk + +import pysitk.python_helper as ph + +import niftymic.base.stack as st +import niftymic.validation.residual_evaluator as res_ev +import niftymic.base.exceptions as exceptions +from niftymic.definitions import DIR_TMP, DIR_TEST + + +class ResidualEvaluatorTest(unittest.TestCase): + + def setUp(self): + self.precision = 7 + self.dir_tmp = os.path.join(DIR_TMP, "residual_evaluator") + + def test_compute_write_read_slice_similarities(self): + + paths_to_stacks = [ + os.path.join( + DIR_TEST, "fetal_brain_%d.nii.gz" % d) for d in range(0, 3) + ] + path_to_reference = os.path.join( + DIR_TEST, "FetalBrain_reconstruction_3stacks_myAlg.nii.gz") + + stacks = [ + st.Stack.from_filename(p, ph.append_to_filename(p, "_mask")) + for p in paths_to_stacks + ] + reference = st.Stack.from_filename( + path_to_reference, extract_slices=False) + + residual_evaluator = res_ev.ResidualEvaluator(stacks, reference) + residual_evaluator.compute_slice_projections() + residual_evaluator.evaluate_slice_similarities() + residual_evaluator.write_slice_similarities(self.dir_tmp) + slice_similarities = residual_evaluator.get_slice_similarities() + + residual_evaluator1 = res_ev.ResidualEvaluator() + residual_evaluator1.read_slice_similarities(self.dir_tmp) + slice_similarities1 = residual_evaluator1.get_slice_similarities() + + for stack_name in slice_similarities.keys(): + for m in slice_similarities[stack_name].keys(): + rho_res = slice_similarities[stack_name][m] + rho_res1 = slice_similarities1[stack_name][m] + error = np.linalg.norm(rho_res - rho_res1) + self.assertAlmostEqual(error, 0, places=self.precision) + + def test_slice_projections_not_created(self): + paths_to_stacks = [ + os.path.join( + DIR_TEST, "fetal_brain_%d.nii.gz" % d) for d in range(0, 1) + ] + path_to_reference = os.path.join( + DIR_TEST, "FetalBrain_reconstruction_3stacks_myAlg.nii.gz") + + stacks = [ + st.Stack.from_filename(p, ph.append_to_filename(p, "_mask")) + for p in paths_to_stacks + ] + reference = st.Stack.from_filename( + path_to_reference, extract_slices=False) + + residual_evaluator = res_ev.ResidualEvaluator(stacks, reference) + self.assertRaises(exceptions.ObjectNotCreated, lambda: + residual_evaluator.evaluate_slice_similarities()) \ No newline at end of file diff --git a/tests/run_tests.py b/tests/run_tests.py index 1927e45d..0cc549d1 100755 --- a/tests/run_tests.py +++ b/tests/run_tests.py @@ -9,22 +9,22 @@ # -# Import libraries import unittest -import sys -import os -# Import modules for unit testing from brain_stripping_test import * +from case_study_fetal_brain_test import * from cpp_itk_registration_test import * +from data_reader_test import * +from image_similarity_evaluator_test import * from intensity_correction_test import * -from intra_stack_registration_test import * +# from intra_stack_registration_test import * # TBC from linear_operators_test import * from niftyreg_test import * -from parameter_normalization_test import * -from registration_test import * +# from parameter_normalization_test import * +# from registration_test import * # TBC +from residual_evaluator_test import * from segmentation_propagation_test import * -from simulator_slice_acquisition_test import * +# from simulator_slice_acquisition_test import * from stack_test import * if __name__ == '__main__': diff --git a/tests/segmentation_propagation_test.py b/tests/segmentation_propagation_test.py index b1668f88..8f6ebc79 100644 --- a/tests/segmentation_propagation_test.py +++ b/tests/segmentation_propagation_test.py @@ -15,6 +15,7 @@ import pysitk.simple_itk_helper as sitkh import niftymic.base.stack as st import niftymic.registration.simple_itk_registration as regsitk +import niftymic.registration.niftyreg as nreg import niftymic.utilities.segmentation_propagation as segprop from niftymic.definitions import DIR_TEST @@ -33,7 +34,8 @@ def test_registration(self): filename = "fetal_brain_0" - parameters_gd = (0.1, 0.2, -0.3, 0, -4, 10) + parameters_gd = (0.1, 0.2, -0.3, 0, 0, 0) + # parameters_gd = np.zeros(6) template = st.Stack.from_filename( os.path.join(self.dir_test_data, filename + ".nii.gz"), @@ -49,8 +51,6 @@ def test_registration(self): stack = st.Stack.from_sitk_image(stack_sitk, filename="stack") - # sitkh.show_sitk_image([template.sitk, stack_sitk]) - optimizer = "RegularStepGradientDescent" optimizer_params = { 'learningRate': 1, @@ -58,19 +58,14 @@ def test_registration(self): 'numberOfIterations': 300 } - # optimizer="ConjugateGradientLineSearch" - # optimizer_params="{'learningRate': 1, 'numberOfIterations': 100}" - registration = regsitk.SimpleItkRegistration( - initializer_type="MOMENTS", + initializer_type="SelfGEOMETRY", use_verbose=True, metric="MeanSquares", optimizer=optimizer, optimizer_params=optimizer_params, - use_multiresolution_framework=True, + use_multiresolution_framework=False, ) - # registration = regitk.CppItkRegistration() - # registration = regniftyreg.RegAladin() segmentation_propagation = segprop.SegmentationPropagation( stack=stack, @@ -80,14 +75,13 @@ def test_registration(self): ) segmentation_propagation.run_segmentation_propagation() foo = segmentation_propagation.get_segmented_stack() - # sitkh.show_stacks( - # [template, foo, stack], - # label=["template", "stack_prop", "stack_orig"] - # ) + # Get transform and force center = 0 transform = segmentation_propagation.get_registration_transform_sitk() + transform = sitkh.get_composite_sitk_euler_transform( + transform, sitk.Euler3DTransform()) parameters = sitk.Euler3DTransform( transform.GetInverse()).GetParameters() self.assertEqual(np.round( - np.linalg.norm(np.array(parameters) - parameters_gd), decimals=0), 0) + np.linalg.norm(np.array(parameters) - parameters_gd), decimals=4), 0) diff --git a/tests/simulator_slice_acqusition_test.py b/tests/simulator_slice_acquisition_test.py similarity index 99% rename from tests/simulator_slice_acqusition_test.py rename to tests/simulator_slice_acquisition_test.py index 5cc4e8eb..2a0d3e97 100644 --- a/tests/simulator_slice_acqusition_test.py +++ b/tests/simulator_slice_acquisition_test.py @@ -5,19 +5,18 @@ # \date May 2016 -import unittest -# Import libraries -import SimpleITK as sitk +import os import itk +import unittest import numpy as np -import os +import SimpleITK as sitk + +import pysitk.simple_itk_helper as sitkh import niftymic.base.psf as psf -# Import modules import niftymic.base.stack as st import niftymic.prototyping.simulator_slice_acqusition as sa -import pysitk.simple_itk_helper as sitkh from niftymic.definitions import DIR_TEST # Pixel type of used 3D ITK image diff --git a/tests/stack_test.py b/tests/stack_test.py index 3b78ea0f..50b3b1e5 100644 --- a/tests/stack_test.py +++ b/tests/stack_test.py @@ -12,8 +12,11 @@ import os import niftymic.base.stack as st -from niftymic.definitions import DIR_TEST +import niftymic.base.data_reader as dr import niftymic.base.exceptions as exceptions +import niftymic.validation.motion_simulator as ms + +from niftymic.definitions import DIR_TEST, DIR_TMP class StackTest(unittest.TestCase): @@ -140,16 +143,94 @@ def test_delete_slices(self): for i in range(stack.get_number_of_slices()): # print ("----") slice_numbers = [s.get_slice_number() for s in stack.get_slices()] - # print slice_numbers + # print(slice_numbers) indices = np.arange(len(slice_numbers)) random.shuffle(indices) - index = indices[0] - # print index + index = slice_numbers[indices[0]] + # print(index) stack.delete_slice(index) - # print stack.get_number_of_slices() + # print(stack.get_number_of_slices()) # No slice left at the end of the loop self.assertEqual(stack.get_number_of_slices(), 0) # No slice left for deletion - self.assertRaises(RuntimeError, lambda: stack.delete_slice(-1)) + self.assertRaises(ValueError, lambda: stack.delete_slice(-1)) + + def test_update_write_transform(self): + + motion_simulator = ms.RandomRigidMotionSimulator( + dimension=3, + angle_max_deg=20, + translation_max=30) + + filenames = ["fetal_brain_%d" % d for d in range(3)] + stacks = [st.Stack.from_filename( + os.path.join(self.dir_test_data, "%s.nii.gz" % f)) + for f in filenames] + + # Generate random motions for all slices of each stack + motions_sitk = {f: {} for f in filenames} + for i, stack in enumerate(stacks): + motion_simulator.simulate_motion( + seed=i, simulations=stack.get_number_of_slices()) + motions_sitk[stack.get_filename()] = \ + motion_simulator.get_transforms_sitk() + + # Apply random motion to all slices of all stacks + dir_output = os.path.join(DIR_TMP, "test_update_write_transform") + for i, stack in enumerate(stacks): + for j, slice in enumerate(stack.get_slices()): + slice.update_motion_correction( + motions_sitk[stack.get_filename()][j]) + + # Write stacks to directory + stack.write(dir_output, write_slices=True, write_transforms=True) + + # Read written stacks/slices/transformations + data_reader = dr.ImageSlicesDirectoryReader( + dir_output) + data_reader.read_data() + stacks_2 = data_reader.get_data() + + data_reader = dr.SliceTransformationDirectoryReader( + dir_output) + data_reader.read_data() + transformations_dic = data_reader.get_data() + + filenames_2 = [s.get_filename() for s in stacks_2] + for i, stack in enumerate(stacks): + stack_2 = stacks_2[filenames_2.index(stack.get_filename())] + slices = stack.get_slices() + slices_2 = stack_2.get_slices() + + # test number of slices match + self.assertEqual(len(slices), len(slices_2)) + + # Test whether header of written slice coincides with transformed + # slice + for j in range(stack.get_number_of_slices()): + + # Check Spacing + self.assertAlmostEqual( + np.max(np.abs(np.array(slices[j].sitk.GetSpacing()) - + np.array(slices_2[j].sitk.GetSpacing()))), + 0, places=10) + # Check Origin + self.assertAlmostEqual( + np.max(np.abs(np.array(slices[j].sitk.GetOrigin()) - + np.array(slices_2[j].sitk.GetOrigin()))), + 0, places=4) + # Check Direction + self.assertAlmostEqual( + np.max(np.abs(np.array(slices[j].sitk.GetDirection()) - + np.array(slices_2[j].sitk.GetDirection()))), + 0, places=4) + + # Test whether parameters of written slice transforms match + params = np.array( + motions_sitk[stack.get_filename()][j].GetParameters()) + params_2 = np.array( + transformations_dic[stack.get_filename()][j].GetParameters()) + self.assertAlmostEqual( + np.max(np.abs(params - params_2)), 0, places=16)