diff --git a/defaults.py b/defaults.py index 678157a..888be9b 100644 --- a/defaults.py +++ b/defaults.py @@ -1,3 +1,14 @@ -DEFAULT_BATCH_SIZE = 32 +# Default inference batch size - based on RAM of computer +# for 16 GB RAM , recommend 16 +DEFAULT_BATCH_SIZE = 16 -FIX_VISUALIZATION_BOUNDS = True \ No newline at end of file +# Fix the bounds for the visualization colorbar based on the quantitative measure (T2, T1rho, etc) +# if any value exceeds these bounds, automatically set this parameter to False +FIX_VISUALIZATION_BOUNDS = True + +# The dilation rate we use for dilating any mask before registration +DEFAULT_MASK_DIL_RATE = 8.0 +DEFAULT_MASK_DIL_THRESHOLD = 0.2 + +# The R^2 fit threshold to include when estimating quantitative values +DEFAULT_R2_THRESHOLD = 0.9 diff --git a/elastix_params/parameters-affine-interregister.txt b/elastix_params/parameters-affine-interregister.txt new file mode 100644 index 0000000..5b382d6 --- /dev/null +++ b/elastix_params/parameters-affine-interregister.txt @@ -0,0 +1,69 @@ +// ********************* +// * affine +// ********************* + +// ********************* +// * ImageTypes +// ********************* +(FixedInternalImagePixelType "float") +(MovingInternalImagePixelType "float") +(UseDirectionCosines "true") + +// ********************* +// * Components +// ********************* +(FixedImagePyramid "FixedSmoothingImagePyramid") +(MovingImagePyramid "MovingSmoothingImagePyramid") +(Registration "MultiResolutionRegistration") +(Interpolator "BSplineInterpolator") +(ResampleInterpolator "FinalBSplineInterpolator") +(Metric "AdvancedMattesMutualInformation") +(BSplineInterpolationOrder 1) +(Resampler "DefaultResampler") +(Optimizer "AdaptiveStochasticGradientDescent") +(Transform "AffineTransform") + +// ********************* +// * Mask settings +// ********************* +(ErodeMask "false") +(ErodeFixedMask "false") + +// ********************* +// * Optimizer settings +// ********************* +(NumberOfResolutions 1) +(MaximumNumberOfIterations 1000) +(ASGDParameterEstimationMethod "Original") +(AutomaticParameterEstimation "true") +(AutomaticTransformInitialization "true") +(AutomaticScalesEstimation "true") + +// ********************* +// * Transform settings +// ********************* +(HowToCombineTransforms "Compose") + +// ********************* +// * Pyramid settings +// ********************* +(NumberOfHistogramBins 32) + + +// ********************* +// * Sampler parameters +// ********************* +(NumberOfSpatialSamples 2000) +(ImageSampler "RandomSparseMask") +(CheckNumberOfSamples "true") +(NewSamplesEveryIteration "true") +(FinalBSplineInterpolationOrder 3) + +// ********************* +// * Output settings +// ********************* +(DefaultPixelValue 0) +(WriteTransformParametersEachIteration "false") +(WriteResultImage "true") +(ResultImageFormat "nii.gz") +(ResultImagePixelType "float") diff --git a/elastix_params/parameters-rigid-interregister.txt b/elastix_params/parameters-rigid-interregister.txt new file mode 100644 index 0000000..4bc850f --- /dev/null +++ b/elastix_params/parameters-rigid-interregister.txt @@ -0,0 +1,69 @@ +// ********************* +// * rigid +// ********************* + +// ********************* +// * ImageTypes +// ********************* +(FixedInternalImagePixelType "float") +(MovingInternalImagePixelType "float") +(UseDirectionCosines "true") + +// ********************* +// * Components +// ********************* +(FixedImagePyramid "FixedSmoothingImagePyramid") +(MovingImagePyramid "MovingSmoothingImagePyramid") +(Registration "MultiResolutionRegistration") +(Interpolator "BSplineInterpolator") +(ResampleInterpolator "FinalBSplineInterpolator") +(Metric "AdvancedMattesMutualInformation") +(BSplineInterpolationOrder 1) +(Resampler "DefaultResampler") +(Optimizer "AdaptiveStochasticGradientDescent") +(Transform "EulerTransform") + +// ********************* +// * Mask settings +// ********************* +(ErodeMask "false") +(ErodeFixedMask "false") + +// ********************* +// * Optimizer settings +// ********************* +(NumberOfResolutions 3) +(MaximumNumberOfIterations 1000) +(ASGDParameterEstimationMethod "Original") +(AutomaticParameterEstimation "true") +(AutomaticTransformInitialization "true") +(AutomaticScalesEstimation "true") + +// ********************* +// * Transform settings +// ********************* +(HowToCombineTransforms "Compose") + +// ********************* +// * Pyramid settings +// ********************* +(NumberOfHistogramBins 32) + + +// ********************* +// * Sampler parameters +// ********************* +(NumberOfSpatialSamples 2000) +(ImageSampler "RandomSparseMask") +(CheckNumberOfSamples "true") +(NewSamplesEveryIteration "true") +(FinalBSplineInterpolationOrder 3) + +// ********************* +// * Output settings +// ********************* +(DefaultPixelValue 0) +(WriteTransformParametersEachIteration "false") +(WriteResultImage "true") +(ResultImageFormat "nii.gz") +(ResultImagePixelType "float") diff --git a/environment.yml b/environment.yml index 1c72552..85b97eb 100644 --- a/environment.yml +++ b/environment.yml @@ -1,111 +1,111 @@ name: tf channels: - - conda-forge - - defaults +- conda-forge +- defaults dependencies: - - boost=1.68.0=py36h3e44d54_0 - - boost-cpp=1.68.0=h3a22d5f_0 - - bzip2=1.0.6=1 - - icu=58.2=hfc679d8_0 - - blas=1.0=mkl - - ca-certificates=2018.03.07=0 - - certifi=2018.8.24=py36_1 - - cycler=0.10.0=py36hfc81398_0 - - freetype=2.9.1=hb4e5f40_0 - - intel-openmp=2019.0=118 - - kiwisolver=1.0.1=py36h0a44026_0 - - libcxx=4.0.1=h579ed51_0 - - libcxxabi=4.0.1=hebd6815_0 - - libedit=3.1.20170329=hb402a30_2 - - libffi=3.2.1=h475c297_4 - - libgfortran=3.0.1=h93005f0_2 - - libopenblas=0.2.20=hdc02c5d_7 - - libpng=1.6.34=he12f830_0 - - matplotlib=3.0.0=py36h54f8f79_0 - - mkl=2019.0=118 - - mkl_fft=1.0.6=py36hb8a8100_0 - - mkl_random=1.0.1=py36h5d10147_1 - - ncurses=6.1=h0a44026_0 - - numpy=1.15.2=py36h6a91979_0 - - numpy-base=1.15.2=py36h8a80b8c_0 - - openssl=1.0.2p=h1de35cc_0 - - pip=10.0.1=py36_0 - - pyparsing=2.2.0=py36_1 - - python=3.6.6=hc167b69_0 - - python-dateutil=2.7.3=py36_0 - - pytz=2018.5=py36_0 - - readline=7.0=hc1231fa_4 - - setuptools=40.0.0=py36_0 - - six=1.11.0=py36_1 - - sqlite=3.24.0=ha441bb4_0 - - tk=8.6.7=h35a86e2_3 - - tornado=5.1=py36h1de35cc_0 - - wheel=0.31.1=py36_0 - - xz=5.2.4=h1de35cc_4 - - zlib=1.2.11=hf3cbc9b_2 - - pip: - - absl-py==0.4.0 - - apipkg==1.5 - - astor==0.7.1 - - atomicwrites==1.2.1 - - attrs==18.2.0 - - bleach==1.5.0 - - click==6.7 - - configparser==3.5.0 - - decorator==4.3.0 - - deprecated==1.2.1 - - dicom2nifti==2.0.3 - - et-xmlfile==1.0.1 - - execnet==1.5.0 - - funcsigs==1.0.2 - - future==0.16.0 - - gast==0.2.0 - - graphviz==0.9 - - grpcio==1.14.1 - - h5py==2.8.0 - - html5lib==0.9999999 - - isodate==0.6.0 - - jdcal==1.4 - - keras==2.1.6 - - lxml==4.2.5 - - markdown==2.6.11 - - mock==2.0.0 - - more-itertools==4.3.0 - - mpmath==1.0.0 - - natsort==5.3.3 - - networkx==2.1 - - nibabel==2.3.0 - - nipy==0.4.2 - - nipype==1.1.2 - - opencv-python==3.4.2.17 - - openpyxl==2.5.8 - - packaging==17.1 - - pandas==0.23.4 - - pbr==4.2.0 - - pillow==5.2.0 - - pluggy==0.7.1 - - protobuf==3.6.1 - - prov==1.5.0 - - psutil==5.4.7 - - py==1.6.0 - - pydicom==1.1.0 - - pydot==1.2.4 - - pydotplus==2.0.2 - - pytest==3.8.0 - - pytest-forked==0.2 - - pytest-xdist==1.23.0 - - pywavelets==0.5.2 - - pyyaml==3.13 - - rdflib==4.2.2 - - scikit-image==0.13.1 - - scipy==1.1.0 - - simpleitk==1.1.0 - - simplejson==3.16.0 - - sympy==1.3 - - tensorboard==1.8.0 - - tensorflow==1.8.0 - - termcolor==1.1.0 - - traits==4.6.0 - - werkzeug==0.14.1 - - wrapt==1.10.11 +- boost=1.68.0=py36h3e44d54_0 +- boost-cpp=1.68.0=h3a22d5f_0 +- bzip2=1.0.6=1 +- icu=58.2=hfc679d8_0 +- blas=1.0=mkl +- ca-certificates=2018.03.07=0 +- certifi=2018.8.24=py36_1 +- cycler=0.10.0=py36hfc81398_0 +- freetype=2.9.1=hb4e5f40_0 +- intel-openmp=2019.0=118 +- kiwisolver=1.0.1=py36h0a44026_0 +- libcxx=4.0.1=h579ed51_0 +- libcxxabi=4.0.1=hebd6815_0 +- libedit=3.1.20170329=hb402a30_2 +- libffi=3.2.1=h475c297_4 +- libgfortran=3.0.1=h93005f0_2 +- libopenblas=0.2.20=hdc02c5d_7 +- libpng=1.6.34=he12f830_0 +- matplotlib=3.0.0=py36h54f8f79_0 +- mkl=2019.0=118 +- mkl_fft=1.0.6=py36hb8a8100_0 +- mkl_random=1.0.1=py36h5d10147_1 +- ncurses=6.1=h0a44026_0 +- numpy=1.15.2=py36h6a91979_0 +- numpy-base=1.15.2=py36h8a80b8c_0 +- openssl=1.0.2p=h1de35cc_0 +- pip=10.0.1=py36_0 +- pyparsing=2.2.0=py36_1 +- python=3.6.6=hc167b69_0 +- python-dateutil=2.7.3=py36_0 +- pytz=2018.5=py36_0 +- readline=7.0=hc1231fa_4 +- setuptools=40.0.0=py36_0 +- six=1.11.0=py36_1 +- sqlite=3.24.0=ha441bb4_0 +- tk=8.6.7=h35a86e2_3 +- tornado=5.1=py36h1de35cc_0 +- wheel=0.31.1=py36_0 +- xz=5.2.4=h1de35cc_4 +- zlib=1.2.11=hf3cbc9b_2 +- pip: + - absl-py==0.4.0 + - apipkg==1.5 + - astor==0.7.1 + - atomicwrites==1.2.1 + - attrs==18.2.0 + - bleach==1.5.0 + - click==6.7 + - configparser==3.5.0 + - decorator==4.3.0 + - deprecated==1.2.1 + - dicom2nifti==2.0.3 + - et-xmlfile==1.0.1 + - execnet==1.5.0 + - funcsigs==1.0.2 + - future==0.16.0 + - gast==0.2.0 + - graphviz==0.9 + - grpcio==1.14.1 + - h5py==2.8.0 + - html5lib==0.9999999 + - isodate==0.6.0 + - jdcal==1.4 + - keras==2.1.6 + - lxml==4.2.5 + - markdown==2.6.11 + - mock==2.0.0 + - more-itertools==4.3.0 + - mpmath==1.0.0 + - natsort==5.3.3 + - networkx==2.1 + - nibabel==2.3.0 + - nipy==0.4.2 + - nipype==1.1.2 + - opencv-python==3.4.2.17 + - openpyxl==2.5.8 + - packaging==17.1 + - pandas==0.23.4 + - pbr==4.2.0 + - pillow==5.2.0 + - pluggy==0.7.1 + - protobuf==3.6.1 + - prov==1.5.0 + - psutil==5.4.7 + - py==1.6.0 + - pydicom==1.1.0 + - pydot==1.2.4 + - pydotplus==2.0.2 + - pytest==3.8.0 + - pytest-forked==0.2 + - pytest-xdist==1.23.0 + - pywavelets==0.5.2 + - pyyaml==3.13 + - rdflib==4.2.2 + - scikit-image==0.13.1 + - scipy==1.1.0 + - simpleitk==1.1.0 + - simplejson==3.16.0 + - sympy==1.3 + - tensorboard==1.8.0 + - tensorflow==1.8.0 + - termcolor==1.1.0 + - traits==4.6.0 + - werkzeug==0.14.1 + - wrapt==1.10.11 prefix: /Users/arjundesai/anaconda3/envs/tf diff --git a/file_constants.py b/file_constants.py index 1330d5c..27e9f1b 100644 --- a/file_constants.py +++ b/file_constants.py @@ -11,6 +11,10 @@ ELASTIX_BSPLINE_PARAMS_FILE = os.path.join(__PATH_TO_ELASTIX_FOLDER__, 'parameters-bspline.txt') ELASTIX_RIGID_PARAMS_FILE = os.path.join(__PATH_TO_ELASTIX_FOLDER__, 'parameters-rigid.txt') +ELASTIX_AFFINE_INTERREGISTER_PARAMS_FILE = os.path.join(__PATH_TO_ELASTIX_FOLDER__, + 'parameters-affine-interregister.txt') +ELASTIX_RIGID_INTERREGISTER_PARAMS_FILE = os.path.join(__PATH_TO_ELASTIX_FOLDER__, + 'parameters-rigid-intraregister.txt') # Temporary file path TEMP_FOLDER_PATH = os.path.join(__DIR__, 'temp') @@ -18,7 +22,8 @@ # TODO: nipype logging NIPYPE_LOGGING = 'none' + def set_debug(): global DEBUG, NIPYPE_LOGGING - DEBUG=1 - NIPYPE_LOGGING = 'stream' \ No newline at end of file + DEBUG = 1 + NIPYPE_LOGGING = 'stream' diff --git a/med_objects/__init__.py b/med_objects/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/med_objects/med_volume.py b/med_objects/med_volume.py new file mode 100644 index 0000000..f0b5bb3 --- /dev/null +++ b/med_objects/med_volume.py @@ -0,0 +1,10 @@ +from utils import io_utils + + +class MedicalVolume(): + def __init__(self, volume, pixel_spacing): + self.volume = volume + self.pixel_spacing = pixel_spacing + + def save_volume(self, filename): + io_utils.save_nifti(filename, self.volume, self.pixel_spacing) diff --git a/models/get_model.py b/models/get_model.py index 06f118a..aff03a7 100644 --- a/models/get_model.py +++ b/models/get_model.py @@ -2,8 +2,9 @@ SUPPORTED_MODELS = ['unet2d'] + def get_model(model_str, input_shape, weights_path): if model_str == 'unet2d': return Unet2D(input_shape, weights_path) - raise ValueError('%s model type not supported' % model_str) \ No newline at end of file + raise ValueError('%s model type not supported' % model_str) diff --git a/models/unet2d.py b/models/unet2d.py index 7d0bc51..c8d4150 100644 --- a/models/unet2d.py +++ b/models/unet2d.py @@ -1,8 +1,8 @@ +import keras.backend as K import numpy as np -from keras.models import Model -from keras.layers import Input, Conv2D, MaxPooling2D, Conv2DTranspose, concatenate, Dropout, Concatenate from keras.layers import BatchNormalization as BN -import keras.backend as K +from keras.layers import Input, Conv2D, MaxPooling2D, Conv2DTranspose, Dropout, Concatenate +from keras.models import Model from models.model import SegModel @@ -76,7 +76,7 @@ def __load_keras_model__(self, input_shape): up = Concatenate(axis=3)([Conv2DTranspose(nfeatures[depth_cnt], (3, 3), padding='same', strides=unpooling_size)(conv), - conv_ptr[depth_cnt]]) + conv_ptr[depth_cnt]]) conv = Conv2D(nfeatures[depth_cnt], (3, 3), padding='same', @@ -125,4 +125,3 @@ def generate_mask(self, volume): K.clear_session() return mask - diff --git a/msk/knee.py b/msk/knee.py index 11900e7..f6c9390 100644 --- a/msk/knee.py +++ b/msk/knee.py @@ -1,6 +1,6 @@ from tissues.femoral_cartilage import FemoralCartilage -from utils.quant_vals import QuantitativeValues as QV from utils import io_utils +from utils.quant_vals import QuantitativeValues as QV KNEE_KEY = 'knee' ORIENTATION_KEY = 'orientation' @@ -47,7 +47,7 @@ def handle_knee(vargin): tissue.load_data(load_path) print('') - print('=='*40) + print('==' * 40) print(tissue.FULL_NAME) print('==' * 40) @@ -55,9 +55,8 @@ def handle_knee(vargin): # load file print('Analyzing %s' % qv.name.lower()) filepath = find_filepath_with_qv(load_path, qv) - tmp = io_utils.load_h5(filepath) - qv_map = tmp['data'] - tissue.calc_quant_vals(qv_map, qv) + tmp = io_utils.load_nifti(filepath) + tissue.calc_quant_vals(tmp, qv) for tissue in tissues: tissue.save_data(vargin[SAVE_KEY]) @@ -67,12 +66,12 @@ def handle_knee(vargin): def find_filepath_with_qv(load_path, qv): import glob, os - dirlist = glob.glob(os.path.join(load_path, '*', '%s.h5' % qv.name.lower())) + dirlist = glob.glob(os.path.join(load_path, '*', '%s.nii.gz' % qv.name.lower())) name = qv.name.lower() if len(dirlist) == 0: - raise ValueError('No map for %s found. Must have name %s.h5' % (name, name)) + raise ValueError('No map for %s found. Must have name %s.nii.gz' % (name, name)) if len(dirlist) > 1: raise ValueError('Multiple %s maps found. Delete extra %s maps' % (name, name)) diff --git a/pipeline.py b/pipeline.py index a67c485..68cae14 100644 --- a/pipeline.py +++ b/pipeline.py @@ -2,21 +2,19 @@ Main file for scan pipeline - handle argparse """ import argparse -import os, time +import os +import time +import defaults +import file_constants as fc from models.get_model import SUPPORTED_MODELS -from scan_sequences.dess import Dess -from scan_sequences.cube_quant import CubeQuant -from scan_sequences.cones import Cones from models.get_model import get_model -from tissues.femoral_cartilage import FemoralCartilage - -from utils.quant_vals import QuantitativeValues as QV, get_qv -import file_constants as fc - from msk import knee - -import defaults +from scan_sequences.cones import Cones +from scan_sequences.cube_quant import CubeQuant +from scan_sequences.dess import Dess +from tissues.femoral_cartilage import FemoralCartilage +from utils.quant_vals import QuantitativeValues as QV SUPPORTED_TISSUES = [FemoralCartilage()] SUPPORTED_QUANTITATIVE_VALUES = [QV.T2, QV.T1_RHO, QV.T2_STAR] @@ -49,7 +47,7 @@ TARGET_MASK_KEY = 'tm' INTERREGISTERED_FILES_DIR_KEY = 'd' -ORIENTATION_KEY='orientation' +ORIENTATION_KEY = 'orientation' TISSUES_KEY = 'tissues' @@ -74,13 +72,13 @@ def add_segmentation_subparser(parser): parser_segment = parser.add_parser('segment') parser_segment.add_argument('--%s' % SEGMENTATION_MODEL_KEY, choices=SUPPORTED_MODELS, nargs='?', default='unet2d') parser_segment.add_argument('--%s' % SEGMENTATION_WEIGHTS_DIR_KEY, type=str, nargs=1, - help='path to directory with weights') + help='path to directory with weights') parser_segment.add_argument('--%s' % SEGMENTATION_BATCH_SIZE_KEY, metavar='B', type=int, default=defaults.DEFAULT_BATCH_SIZE, nargs='?', help='batch size for inference. Default: 32') for tissue in SUPPORTED_TISSUES: parser_segment.add_argument('-%s' % tissue.STR_ID, action='store_const', default=False, const=True, - help='handle %s' % tissue.FULL_NAME) + help='handle %s' % tissue.FULL_NAME) def add_interregister_subparser(parser): @@ -169,8 +167,8 @@ def handle_cubequant(vargin): def handle_cones(vargin): print('\nAnalyze cones') scan = Cones(dicom_path=vargin[DICOM_KEY], - dicom_ext=vargin[EXT_KEY], - load_path=vargin[LOAD_KEY]) + dicom_ext=vargin[EXT_KEY], + load_path=vargin[LOAD_KEY]) if vargin[FOCUSED_MASK_KEY]: scan.focused_mask_filepath = vargin[FOCUSED_MASK_KEY] @@ -249,11 +247,11 @@ def parse_args(): # Cones parser parser_cones = subparsers.add_parser(CONES_KEY, help='analyze cones sequence') parser_cones.add_argument('-%s' % T2_STAR_KEY, action='store_const', default=False, const=True, - help='do t2* analysis') + help='do t2* analysis') parser_cones.add_argument('-%s' % FOCUSED_MASK_KEY, - nargs='?', - default=None, - help='focused mask to speed up t1rho calculation') + nargs='?', + default=None, + help='focused mask to speed up t1rho calculation') subparsers_cones = parser_cones.add_subparsers(help='sub-command help', dest=ACTION_KEY) add_interregister_subparser(subparsers_cones) diff --git a/scan_sequences/cones.py b/scan_sequences/cones.py index 25f71c9..9147778 100644 --- a/scan_sequences/cones.py +++ b/scan_sequences/cones.py @@ -1,19 +1,17 @@ import os -from scan_sequences.scans import NonTargetSequence -from utils import dicom_utils, io_utils -from utils import quant_vals as qv -from nipype.interfaces.elastix import Registration, ApplyWarp import numpy as np -import file_constants as fc -import scipy.ndimage as sni from natsort import natsorted -import re + +import file_constants as fc +from scan_sequences.scans import NonTargetSequence +from utils import io_utils +from utils import quant_vals as qv +from utils.fits import MonoExponentialFit __EXPECTED_NUM_ECHO_TIMES__ = 4 -__R_SQUARED_THRESHOLD__ = 0.9 -__INITIAL_T2_STAR_VAL__ = 30.0 # ms +__INITIAL_T2_STAR_VAL__ = 30.0 # ms __T2_STAR_LOWER_BOUND__ = 0 __T2_STAR_UPPER_BOUND__ = np.inf @@ -27,6 +25,7 @@ def __init__(self, dicom_path=None, dicom_ext=None, load_path=None): super().__init__(dicom_path, dicom_ext) self.t2star_map = None + self.r2 = None self.subvolumes = None self.focused_mask_filepath = None self.echo_times = [] @@ -59,9 +58,7 @@ def interregister(self, target_path, mask_path=None): for i in range(len(echo_time_inds)): raw_filepath = os.path.join(temp_raw_dirpath, '%03d.nii.gz' % i) - io_utils.save_nifti(raw_filepath, - subvolumes[i], - self.pixel_spacing) + subvolumes[i].save_volume(raw_filepath) raw_filepaths[i] = raw_filepath # last echo should be base @@ -77,128 +74,83 @@ def interregister(self, target_path, mask_path=None): print('Mask: %s' % mask_path) print('==' * 40) - # get number of slices - volume_shape = self.subvolumes[base_echo_time].shape - - # Register base image to the target image - print('Registering %s (base image)' % base_image) - reg = Registration() - reg.inputs.fixed_image = target_path - reg.inputs.moving_image = base_image - reg.inputs.output_path = io_utils.check_dir(os.path.join(temp_interregistered_dirpath, - '%03d' % base_echo_time)) - reg.inputs.parameters = [fc.ELASTIX_RIGID_PARAMS_FILE, fc.ELASTIX_AFFINE_PARAMS_FILE] - - if mask_path is not None: - #raise ValueError('Moving mask not supported') - - mask, mask_spacing = io_utils.load_nifti(mask_path) - - fixed_mask = np.asarray(sni.gaussian_filter(np.asarray(mask, dtype=np.float32), sigma=3.0) > 0.2, - dtype=np.int8) - fixed_mask_filepath = os.path.join(temp_interregistered_dirpath, 'dilated-mask.nii.gz') - io_utils.save_nifti(fixed_mask_filepath, fixed_mask, mask_spacing) - reg.inputs.fixed_mask = fixed_mask_filepath - - reg.terminal_output = fc.NIPYPE_LOGGING - reg_output = reg.run() - reg_output = reg_output.outputs - transformation_files = reg_output.transform - warped_files = [(base_echo_time, reg_output.warped_file)] - - files = [] + files_to_warp = [] for echo_time_ind in raw_filepaths.keys(): + if echo_time_ind == base_echo_time: + continue filepath = raw_filepaths[echo_time_ind] - files.append((echo_time_ind, filepath)) + files_to_warp.append((echo_time_ind, filepath)) - # Load the transformation file. Apply same transform to the remaining images - for echo_time_ind, filename in files: - print('Applying transform %s' % filename) - warped_file = '' - for f in transformation_files: - reg = ApplyWarp() - reg.inputs.moving_image = filename - - reg.inputs.transform_file = f - reg.inputs.output_path = io_utils.check_dir(os.path.join(temp_interregistered_dirpath, - '%03d' % echo_time_ind)) - reg.terminal_output = fc.NIPYPE_LOGGING - reg_output = reg.run() - warped_file = reg_output.outputs.warped_file - - assert warped_file != '' + if not mask_path: + parameter_files = [fc.ELASTIX_RIGID_PARAMS_FILE, fc.ELASTIX_AFFINE_PARAMS_FILE] + else: + parameter_files = [fc.ELASTIX_RIGID_INTERREGISTER_PARAMS_FILE, fc.ELASTIX_AFFINE_INTERREGISTER_PARAMS_FILE] + warped_file, transformation_files = self.__interregister_base_file__((base_image, base_echo_time), + target_path, + temp_interregistered_dirpath, + mask_path=mask_path, + parameter_files=parameter_files) + warped_files = [(base_echo_time, warped_file)] + + # Load the transformation file. Apply same transform to the remaining images + for echo_time, filename in files_to_warp: + warped_file = self.__apply_transform__((filename, echo_time), + transformation_files, + temp_interregistered_dirpath) # append the last warped file - this has all the transforms applied - warped_files.append((echo_time_ind, warped_file)) + warped_files.append((echo_time, warped_file)) # copy each of the interregistered warped files to their own output subvolumes = dict() - for echo_time_ind, warped_file in warped_files: - subvolumes[echo_time_ind], _ = io_utils.load_nifti(warped_file) + for echo_time, warped_file in warped_files: + subvolumes[echo_time] = io_utils.load_nifti(warped_file) self.subvolumes = subvolumes - def save_data(self, save_dirpath): - super().save_data(save_dirpath) - save_dirpath = self.__save_dir__(save_dirpath) - - if self.t2star_map is not None: - data = {'data': self.t2star_map} - io_utils.save_h5(os.path.join(save_dirpath, '%s.h5' % qv.QuantitativeValues.T2_STAR.name.lower()), data) - io_utils.save_nifti(os.path.join(save_dirpath, '%s.nii.gz' % qv.QuantitativeValues.T2_STAR.name.lower()), - self.t2star_map, self.pixel_spacing) - - # Save interregistered files - interregistered_dirpath = os.path.join(save_dirpath, 'interregistered') - - for spin_lock_time in self.subvolumes.keys(): - filepath = os.path.join(interregistered_dirpath, '%03d.nii.gz' % spin_lock_time) - io_utils.save_nifti(filepath, self.subvolumes[spin_lock_time], self.pixel_spacing) - def generate_t2_star_map(self): - svs = [] + msk = None spin_lock_times = [] - original_shape = None + subvolumes_list = [] if self.focused_mask_filepath: print('Using focused mask: %s' % self.focused_mask_filepath) - msk, _ = io_utils.load_nifti(self.focused_mask_filepath) - msk = msk.reshape(1, -1) + msk = io_utils.load_nifti(self.focused_mask_filepath) - for spin_lock_time in self.subvolumes.keys(): - spin_lock_times.append(spin_lock_time) - sv = self.subvolumes[spin_lock_time] + for echo_time in self.subvolumes.keys(): + spin_lock_times.append(echo_time) + subvolumes_list.append(self.subvolumes[echo_time]) - if original_shape is None: - original_shape = sv.shape - else: - assert (sv.shape == original_shape) + mef = MonoExponentialFit(spin_lock_times, + subvolumes_list, + mask=msk, + bounds=(__T2_STAR_LOWER_BOUND__, __T2_STAR_UPPER_BOUND__), + tc0=__INITIAL_T2_STAR_VAL__, + decimal_precision=__T2_STAR_DECIMAL_PRECISION__) - svr = sv.reshape((1, -1)) - if self.focused_mask_filepath: - svr = svr * msk + self.t2star_map, self.r2 = mef.fit() - svs.append(svr) + return self.t2star_map - svs = np.concatenate(svs) - vals, r_squared = qv.fit_monoexp_tc(spin_lock_times, svs, __INITIAL_T2_STAR_VAL__) - - map_unfiltered = vals.reshape(original_shape) - r_squared = r_squared.reshape(original_shape) - - t2star_map = map_unfiltered * (r_squared > __R_SQUARED_THRESHOLD__) + def save_data(self, save_dirpath): + super().save_data(save_dirpath) + save_dirpath = self.__save_dir__(save_dirpath) - # Filter calculated T1-rho values that are below 0ms and over 100ms - t2star_map[t2star_map <= __T2_STAR_LOWER_BOUND__] = np.nan - t2star_map = np.nan_to_num(t2star_map) - t2star_map[t2star_map > __T2_STAR_UPPER_BOUND__] = np.nan - t2star_map = np.nan_to_num(t2star_map) + if self.t2star_map is not None: + assert self.r2 is not None + self.t2star_map.save_volume(os.path.join(save_dirpath, + '%s.nii.gz' % qv.QuantitativeValues.T2_STAR.name.lower())) - t2star_map = np.around(t2star_map, __T2_STAR_DECIMAL_PRECISION__) + t1rho_r2_map_filepath = os.path.join(save_dirpath, + '%s_r2.nii.gz' % qv.QuantitativeValues.T2_STAR.name.lower()) + self.r2.save_volume(t1rho_r2_map_filepath) - self.t2star_map = t2star_map + # Save interregistered files + interregistered_dirpath = os.path.join(save_dirpath, 'interregistered') - return t2star_map + for spin_lock_time in self.subvolumes.keys(): + filepath = os.path.join(interregistered_dirpath, '%03d.nii.gz' % spin_lock_time) + self.subvolumes[spin_lock_time].save_volume(filepath) def load_data(self, load_dirpath): super().load_data(load_dirpath) @@ -208,41 +160,6 @@ def load_data(self, load_dirpath): self.subvolumes = self.__load_interregistered_files__(interregistered_dirpath) - def __load_interregistered_files__(self, interregistered_dirpath): - if 'interregistered' not in interregistered_dirpath: - raise ValueError('Invalid path for loading %s interregistered files' % self.NAME) - - subfiles = os.listdir(interregistered_dirpath) - subfiles = natsorted(subfiles) - - if len(subfiles) == 0: - raise ValueError('No interregistered files found') - - spin_lock_times = [] - subvolumes = [] - for subfile in subfiles: - subfile_nums = re.findall(r"[-+]?\d*\.\d+|\d+", subfile) - if len(subfile_nums) == 0: - raise ValueError('%s is not an interregisterd \'.gz.nii\' file.' % subfile) - - subfile_num = float(subfile_nums[0]) - spin_lock_times.append(subfile_num) - - filepath = os.path.join(interregistered_dirpath, subfile) - subvolume_arr, self.pixel_spacing = io_utils.load_nifti(filepath) - subvolumes.append(subvolume_arr) - - assert len(spin_lock_times) == len(subvolumes), "Number of subvolumes mismatch" - - if len(subvolumes) == 0: - raise ValueError('No interregistered files found') - - subvolumes_dict = dict() - for i in range(len(spin_lock_times)): - subvolumes_dict[spin_lock_times[i]] = subvolumes[i] - - return subvolumes_dict - def __serializable_variables__(self): var_names = super().__serializable_variables__() var_names.extend(['echo_times']) @@ -251,4 +168,4 @@ def __serializable_variables__(self): if __name__ == '__main__': - cq = Cones('../dicoms/healthy07/009', 'dcm', './') \ No newline at end of file + cq = Cones('../dicoms/healthy07/009', 'dcm', './') diff --git a/scan_sequences/cube_quant.py b/scan_sequences/cube_quant.py index 5b9cc7b..00818b8 100644 --- a/scan_sequences/cube_quant.py +++ b/scan_sequences/cube_quant.py @@ -1,21 +1,20 @@ import os -from scan_sequences.scans import NonTargetSequence from natsort import natsorted -from nipype.interfaces.elastix import Registration, ApplyWarp -from utils import io_utils -import scipy.ndimage as sni -import numpy as np +from nipype.interfaces.elastix import Registration + import file_constants as fc +from scan_sequences.scans import NonTargetSequence +from utils import io_utils from utils import quant_vals as qv - +from utils.fits import MonoExponentialFit __EXPECTED_NUM_SPIN_LOCK_TIMES__ = 4 __R_SQUARED_THRESHOLD__ = 0.9 __INITIAL_T1_RHO_VAL__ = 70.0 __T1_RHO_LOWER_BOUND__ = 0 -__T1_RHO_UPPER_BOUND__ = np.inf +__T1_RHO_UPPER_BOUND__ = 500 __T1_RHO_DECIMAL_PRECISION__ = 3 @@ -26,6 +25,7 @@ def __init__(self, dicom_path=None, dicom_ext=None, load_path=None): super().__init__(dicom_path, dicom_ext, load_path=load_path) self.t1rho_map = None + self.r2 = None self.subvolumes = None self.focused_mask_filepath = None @@ -53,108 +53,55 @@ def interregister(self, target_path, mask_path=None): print('Mask: %s' % mask_path) print('==' * 40) - # get number of slices - volume_shape = self.subvolumes[base_spin_lock_time].shape + if not mask_path: + parameter_files = [fc.ELASTIX_RIGID_PARAMS_FILE, fc.ELASTIX_AFFINE_PARAMS_FILE] + else: + parameter_files = [fc.ELASTIX_RIGID_INTERREGISTER_PARAMS_FILE, fc.ELASTIX_AFFINE_INTERREGISTER_PARAMS_FILE] - # Register base image to the target image - print('Registering %s (base image)' % base_image) - reg = Registration() - reg.inputs.fixed_image = target_path - reg.inputs.moving_image = base_image - reg.inputs.output_path = io_utils.check_dir(os.path.join(temp_interregistered_dirpath, - '%03d' % base_spin_lock_time)) - reg.inputs.parameters = [fc.ELASTIX_RIGID_PARAMS_FILE, fc.ELASTIX_AFFINE_PARAMS_FILE] - - if mask_path is not None: - #raise ValueError('Moving mask not supported') - - mask, mask_spacing = io_utils.load_nifti(mask_path) - - fixed_mask = np.asarray(sni.gaussian_filter(np.asarray(mask, dtype=np.float32), sigma=3.0) > 0.2, dtype=np.int8) - fixed_mask_filepath = os.path.join(temp_interregistered_dirpath, 'dilated-mask.nii.gz') - io_utils.save_nifti(fixed_mask_filepath, fixed_mask, mask_spacing) - reg.inputs.fixed_mask = fixed_mask_filepath - - reg.terminal_output = fc.NIPYPE_LOGGING - reg_output = reg.run() - reg_output = reg_output.outputs - transformation_files = reg_output.transform - warped_files = [(base_spin_lock_time, reg_output.warped_file)] + warped_file, transformation_files = self.__interregister_base_file__((base_image, base_spin_lock_time), + target_path, + temp_interregistered_dirpath, + mask_path=mask_path, + parameter_files=parameter_files) + warped_files = [(base_spin_lock_time, warped_file)] # Load the transformation file. Apply same transform to the remaining images for spin_lock_time, filename in files: - print('Applying transform %s' % filename) - warped_file = '' - for f in transformation_files: - reg = ApplyWarp() - reg.inputs.moving_image = filename - - reg.inputs.transform_file = f - reg.inputs.output_path = io_utils.check_dir(os.path.join(temp_interregistered_dirpath, - '%03d' % spin_lock_time)) - reg.terminal_output = fc.NIPYPE_LOGGING - reg_output = reg.run() - warped_file = reg_output.outputs.warped_file - - assert warped_file != '' - + warped_file = self.__apply_transform__((filename, spin_lock_time), + transformation_files, + temp_interregistered_dirpath) # append the last warped file - this has all the transforms applied warped_files.append((spin_lock_time, warped_file)) # copy each of the interregistered warped files to their own output subvolumes = dict() for spin_lock_time, warped_file in warped_files: - subvolumes[spin_lock_time], _ = io_utils.load_nifti(warped_file) + subvolumes[spin_lock_time] = io_utils.load_nifti(warped_file) - #self.subvolumes = self.__load_interregistered_files__(interregistered_dirpath) self.subvolumes = subvolumes def generate_t1_rho_map(self): - svs = [] spin_lock_times = [] - original_shape = None - + subvolumes_list = [] + msk = None if self.focused_mask_filepath: print('Using focused mask: %s' % self.focused_mask_filepath) - msk, _ = io_utils.load_nifti(self.focused_mask_filepath) - msk = msk.reshape(1, -1) + msk = io_utils.load_nifti(self.focused_mask_filepath) sorted_keys = natsorted(list(self.subvolumes.keys())) for spin_lock_time_index in sorted_keys: + subvolumes_list.append(self.subvolumes[spin_lock_time_index]) spin_lock_times.append(self.spin_lock_times[spin_lock_time_index]) - sv = self.subvolumes[spin_lock_time_index] - - if original_shape is None: - original_shape = sv.shape - else: - assert(sv.shape == original_shape) - - svr = sv.reshape((1, -1)) - if self.focused_mask_filepath: - svr = svr * msk - - svs.append(svr) - - svs = np.concatenate(svs) - vals, r_squared = qv.fit_monoexp_tc(spin_lock_times, svs, __INITIAL_T1_RHO_VAL__) + mef = MonoExponentialFit(spin_lock_times, subvolumes_list, + mask=msk, + bounds=(__T1_RHO_LOWER_BOUND__, __T1_RHO_UPPER_BOUND__), + tc0=__INITIAL_T1_RHO_VAL__, + decimal_precision=__T1_RHO_DECIMAL_PRECISION__) - map_unfiltered = vals.reshape(original_shape) - r_squared = r_squared.reshape(original_shape) + self.t1rho_map, self.r2 = mef.fit() - t1rho_map = map_unfiltered * (r_squared > __R_SQUARED_THRESHOLD__) - - # Filter calculated T1-rho values that are below 0ms and over 100ms - t1rho_map[t1rho_map <= __T1_RHO_LOWER_BOUND__] = np.nan - t1rho_map = np.nan_to_num(t1rho_map) - t1rho_map[t1rho_map > __T1_RHO_UPPER_BOUND__] = np.nan - t1rho_map = np.nan_to_num(t1rho_map) - - t1rho_map = np.around(t1rho_map, __T1_RHO_DECIMAL_PRECISION__) - - self.t1rho_map = t1rho_map - - return t1rho_map + return self.t1rho_map def __intraregister__(self, subvolumes): """ @@ -180,12 +127,12 @@ def __intraregister__(self, subvolumes): filepath = os.path.join(raw_volumes_base_path, '%03d' % spin_lock_time_index + '.nii.gz') spin_lock_nii_files.append(filepath) - io_utils.save_nifti(filepath, subvolumes[spin_lock_time_index], self.pixel_spacing) + subvolumes[spin_lock_time_index].save_volume(filepath) target_filepath = spin_lock_nii_files[0] intraregistered_files = [] - for i in range(1,len(spin_lock_nii_files)): + for i in range(1, len(spin_lock_nii_files)): spin_file = spin_lock_nii_files[i] spin_lock_time_index = ordered_spin_lock_time_indices[i] @@ -211,17 +158,21 @@ def save_data(self, save_dirpath): save_dirpath = self.__save_dir__(save_dirpath) if self.t1rho_map is not None: - data = {'data': self.t1rho_map} - io_utils.save_h5(os.path.join(save_dirpath, '%s.h5' % qv.QuantitativeValues.T1_RHO.name.lower()), data) - io_utils.save_nifti(os.path.join(save_dirpath, '%s.nii.gz' % qv.QuantitativeValues.T1_RHO.name.lower()), self.t1rho_map, - self.pixel_spacing) + assert (self.r2 is not None) + + t1rho_map_filepath = os.path.join(save_dirpath, '%s.nii.gz' % qv.QuantitativeValues.T1_RHO.name.lower()) + self.t1rho_map.save_volume(t1rho_map_filepath) + + t1rho_r2_map_filepath = os.path.join(save_dirpath, + '%s_r2.nii.gz' % qv.QuantitativeValues.T1_RHO.name.lower()) + self.r2.save_volume(t1rho_r2_map_filepath) # Save interregistered files interregistered_dirpath = os.path.join(save_dirpath, 'interregistered') for spin_lock_time_index in self.subvolumes.keys(): filepath = os.path.join(interregistered_dirpath, '%03d.nii.gz' % spin_lock_time_index) - io_utils.save_nifti(filepath, self.subvolumes[spin_lock_time_index], self.pixel_spacing) + self.subvolumes[spin_lock_time_index].save_volume(filepath) def load_data(self, load_dirpath): super().load_data(load_dirpath) diff --git a/scan_sequences/dess.py b/scan_sequences/dess.py index dadb751..89c4a24 100644 --- a/scan_sequences/dess.py +++ b/scan_sequences/dess.py @@ -1,11 +1,15 @@ -import math, os +import math +import os + import numpy as np from pydicom.tag import Tag +from med_objects.med_volume import MedicalVolume from scan_sequences.scans import TargetSequence from utils import dicom_utils, im_utils, io_utils from utils.quant_vals import QuantitativeValues + class Dess(TargetSequence): NAME = 'dess' @@ -22,7 +26,7 @@ class Dess(TargetSequence): # Clipping bounds for t2 __T2_LOWER_BOUND__ = 0 __T2_UPPER_BOUND__ = 100 - __T2_DECIMAL_PRECISION__ = 1 # 0.1 ms + __T2_DECIMAL_PRECISION__ = 1 # 0.1 ms use_rms = False @@ -37,7 +41,7 @@ def __init__(self, dicom_path, dicom_ext=None, load_path=None): raise ValueError('dicoms in \'%s\' are not acquired from DESS sequence' % self.dicom_path) def split_volume(self): - volume = self.volume + volume = self.volume.volume echos = self.__NUM_ECHOS__ if len(volume.shape) != self.__VOLUME_DIMENSIONS__: @@ -52,7 +56,7 @@ def split_volume(self): sub_volumes = [] for i in range(echos): - sub_volumes.append(volume[:, :, i::echos]) + sub_volumes.append(MedicalVolume(volume[:, :, i::echos], self.volume.pixel_spacing)) return sub_volumes @@ -73,12 +77,15 @@ def segment(self, model, tissue): else: segmentation_volume = self.subvolumes[0] + pixel_spacing = segmentation_volume.pixel_spacing + segmentation_volume = segmentation_volume.volume + volume = dicom_utils.whiten_volume(segmentation_volume) # Segment tissue and add it to list mask = model.generate_mask(volume) - tissue.set_mask(mask, self.pixel_spacing) + tissue.set_mask(MedicalVolume(mask, pixel_spacing)) self.__add_tissue__(tissue) return mask @@ -99,12 +106,12 @@ def generate_t2_map(self): all invalid pixels are denoted by the value 0 """ - dicom_array = self.volume - ref_dicom = self.ref_dicom - if self.volume is None or self.ref_dicom is None: raise ValueError('volume and ref_dicom fields must be initialized') + dicom_array = self.volume.volume + ref_dicom = self.ref_dicom + if len(dicom_array.shape) != 3: raise ValueError("dicom_array must be 3D volume") @@ -112,8 +119,8 @@ def generate_t2_map(self): subvolumes = self.subvolumes # Split echos - echo_1 = subvolumes[0] - echo_2 = subvolumes[1] + echo_1 = subvolumes[0].volume + echo_2 = subvolumes[1].volume # All timing in seconds TR = float(ref_dicom.RepetitionTime) * 1e-3 @@ -130,7 +137,8 @@ def generate_t2_map(self): dkL = gamma * Gl * Tg # Simply math - k = math.pow((math.sin(alpha / 2)), 2) * (1 + math.exp(-TR / self.__T1__ - TR * math.pow(dkL, 2) * self.__D__)) / ( + k = math.pow((math.sin(alpha / 2)), 2) * ( + 1 + math.exp(-TR / self.__T1__ - TR * math.pow(dkL, 2) * self.__D__)) / ( 1 - math.cos(alpha) * math.exp(-TR / self.__T1__ - TR * math.pow(dkL, 2) * self.__D__)) c1 = (TR - Tg / 3) * (math.pow(dkL, 2)) * self.__D__ @@ -153,22 +161,22 @@ def generate_t2_map(self): t2map = np.around(t2map, self.__T2_DECIMAL_PRECISION__) - self.t2map = t2map + self.t2map = MedicalVolume(t2map, subvolumes[0].pixel_spacing) - return t2map + return self.t2map def save_data(self, save_dirpath): super().save_data(save_dirpath) save_dirpath = self.__save_dir__(save_dirpath) - data = {'data': self.t2map} - io_utils.save_h5(os.path.join(save_dirpath, '%s.h5' % QuantitativeValues.T2.name.lower()), data) - io_utils.save_nifti(os.path.join(save_dirpath, '%s.nii.gz' % QuantitativeValues.T2.name.lower()), self.t2map, self.pixel_spacing) + + if self.t2map is not None: + self.t2map.save_volume(os.path.join(save_dirpath, '%s.nii.gz' % QuantitativeValues.T2.name.lower())) # write echos for i in range(len(self.subvolumes)): - nii_registration_filepath = os.path.join(save_dirpath, 'echo%d.nii.gz' % (i+1)) - io_utils.save_nifti(nii_registration_filepath, self.subvolumes[i], self.pixel_spacing) + nii_registration_filepath = os.path.join(save_dirpath, 'echo%d.nii.gz' % (i + 1)) + self.subvolumes[i].save_volume(nii_registration_filepath) def load_data(self, load_dirpath): super().load_data(load_dirpath) @@ -176,8 +184,8 @@ def load_data(self, load_dirpath): self.subvolumes = [] # Load subvolumes from nifti file for i in range(self.__NUM_ECHOS__): - nii_registration_filepath = os.path.join(load_dirpath, 'echo%d.nii.gz' % (i+1)) - subvolume, _ = io_utils.load_nifti(nii_registration_filepath) + nii_registration_filepath = os.path.join(load_dirpath, 'echo%d.nii.gz' % (i + 1)) + subvolume = io_utils.load_nifti(nii_registration_filepath) self.subvolumes.append(subvolume) def calc_rms(self): @@ -186,8 +194,8 @@ def calc_rms(self): assert len(self.subvolumes) == 2, "2 Echos expected" - echo1 = np.asarray(self.subvolumes[0], dtype=np.float64) - echo2 = np.asarray(self.subvolumes[1], dtype=np.float64) + echo1 = np.asarray(self.subvolumes[0].volume, dtype=np.float64) + echo2 = np.asarray(self.subvolumes[1].volume, dtype=np.float64) assert (echo1 >= 0).all() assert (echo2 >= 0).all() @@ -200,4 +208,4 @@ def calc_rms(self): rms = np.sqrt(echo1 ** 2 + echo2 ** 2) - return rms + return MedicalVolume(rms, self.subvolumes[0].pixel_spacing) diff --git a/scan_sequences/scans.py b/scan_sequences/scans.py index c367117..ee20e92 100644 --- a/scan_sequences/scans.py +++ b/scan_sequences/scans.py @@ -1,12 +1,18 @@ -from abc import ABC, abstractmethod -from utils import dicom_utils -import numpy as np -from utils import io_utils import os -import file_constants as fc +import re +from abc import ABC, abstractmethod from time import gmtime, strftime + +import numpy as np +import scipy.ndimage as sni from natsort import natsorted -import re +from nipype.interfaces.elastix import Registration, ApplyWarp + +import defaults +import file_constants as fc +from med_objects.med_volume import MedicalVolume +from utils import dicom_utils +from utils import io_utils class ScanSequence(ABC): @@ -20,7 +26,6 @@ def __init__(self, dicom_path=None, dicom_ext=None, load_path=None): self.dicom_ext = dicom_ext self.volume = None - self.pixel_spacing = None self.ref_dicom = None if load_path: @@ -33,12 +38,12 @@ def __load_dicom__(self): dicom_path = self.dicom_path dicom_ext = self.dicom_ext - self.volume, self.refs_dicom, self.pixel_spacing = dicom_utils.load_dicom(dicom_path, dicom_ext) + self.volume, self.refs_dicom = dicom_utils.load_dicom(dicom_path, dicom_ext) self.ref_dicom = self.refs_dicom[0] def get_dimensions(self): - return self.volume.shape + return self.volume.volume.shape def __add_tissue__(self, new_tissue): contains_tissue = any([tissue.ID == new_tissue.ID for tissue in self.tissues]) @@ -90,8 +95,7 @@ def __save_dir__(self, dirpath, create_dir=True): return scan_dirpath def __serializable_variables__(self): - return ['volume', 'dicom_path', 'dicom_ext', 'pixel_spacing'] - + return ['volume', 'dicom_path', 'dicom_ext'] class TargetSequence(ScanSequence): @@ -129,10 +133,8 @@ def interregister(self, target): def __split_volumes__(self, expected_num_echos): refs_dicom = self.refs_dicom volume = self.volume - subvolumes_dict = dict() echo_times = [] - subvolumes = [] # get echo_time time for each refs_dicom for i in range(len(refs_dicom)): echo_time = float(refs_dicom[i].EchoTime) @@ -156,17 +158,20 @@ def __split_volumes__(self, expected_num_echos): # number of spin lock times should go evenly into subvolume if (len(refs_dicom) % num_echo_times) != 0: - raise ValueError('Uneven number of dicom files - %d dicoms, %d spin lock times/echos' % (len(refs_dicom), num_echo_times)) + raise ValueError('Uneven number of dicom files - %d dicoms, %d spin lock times/echos' % ( + len(refs_dicom), num_echo_times)) num_slices_per_subvolume = len(refs_dicom) / num_echo_times for key in subvolumes_dict.keys(): if (len(subvolumes_dict[key]) != num_slices_per_subvolume): - raise ValueError('subvolume \'%s\' (echo_time %s) has %d slices. Expected %d slices' % (key, echo_times[key-1], len(subvolumes_dict[key]), num_slices_per_subvolume)) + raise ValueError('subvolume \'%s\' (echo_time %s) has %d slices. Expected %d slices' % ( + key, echo_times[key - 1], len(subvolumes_dict[key]), num_slices_per_subvolume)) for key in subvolumes_dict.keys(): sv = subvolumes_dict[key] - sv = volume[..., sv] + sv = MedicalVolume(volume.volume[..., sv], volume.pixel_spacing) + subvolumes_dict[key] = sv return subvolumes_dict, echo_times @@ -193,8 +198,9 @@ def __load_interregistered_files__(self, interregistered_dirpath): indices.append(subfile_num) filepath = os.path.join(interregistered_dirpath, subfile) - subvolume_arr, self.pixel_spacing = io_utils.load_nifti(filepath) - subvolumes.append(subvolume_arr) + subvolume = io_utils.load_nifti(filepath) + + subvolumes.append(subvolume) assert len(indices) == len(subvolumes), "Number of subvolumes mismatch" @@ -206,3 +212,67 @@ def __load_interregistered_files__(self, interregistered_dirpath): subvolumes_dict[indices[i]] = subvolumes[i] return subvolumes_dict + + def __dilate_mask__(self, mask_path, temp_path, dil_rate=defaults.DEFAULT_MASK_DIL_RATE, + dil_threshold=defaults.DEFAULT_MASK_DIL_THRESHOLD): + mask = io_utils.load_nifti(mask_path) + mask = mask.volume + dilated_mask = sni.gaussian_filter(np.asarray(mask.volume, dtype=np.float32), + sigma=dil_rate) > dil_threshold + fixed_mask = np.asarray(dilated_mask, + dtype=np.int8) + fixed_mask_filepath = os.path.join(temp_path, 'dilated-mask.nii.gz') + + dilated_mask_volume = MedicalVolume(fixed_mask, mask.pixel_spacing) + dilated_mask_volume.save_volume(fixed_mask_filepath) + + return fixed_mask_filepath + + def __interregister_base_file__(self, base_image_info, target_path, temp_path, mask_path=None, + parameter_files=[fc.ELASTIX_RIGID_PARAMS_FILE, fc.ELASTIX_AFFINE_PARAMS_FILE]): + base_image_path, base_time_id = base_image_info + # Register base image to the target image + print('Registering %s (base image)' % base_image_path) + transformation_files = [] + + use_mask_arr = [False, True] + + for i in range(len(parameter_files)): + use_mask = use_mask_arr[i] + pfile = parameter_files[i] + reg = Registration() + reg.inputs.fixed_image = target_path + reg.inputs.moving_image = base_image_path + reg.inputs.output_path = io_utils.check_dir(os.path.join(temp_path, + '%03d_param%i' % (base_time_id, i))) + reg.inputs.parameters = pfile + + if use_mask and mask_path is not None: + fixed_mask_filepath = self.__dilate_mask__(mask_path, temp_path) + reg.inputs.fixed_mask = fixed_mask_filepath + + reg.terminal_output = fc.NIPYPE_LOGGING + reg_output = reg.run() + reg_output = reg_output.outputs + transformation_files.append(reg_output.transform[0]) + + return reg_output.warped_file, transformation_files + + def __apply_transform__(self, image_info, transformation_files, temp_path): + filename, image_id = image_info + print('Applying transform %s' % filename) + warped_file = '' + for f in transformation_files: + reg = ApplyWarp() + reg.inputs.moving_image = filename + + reg.inputs.transform_file = f + reg.inputs.output_path = io_utils.check_dir(os.path.join(temp_path, + '%03d' % image_id)) + reg.terminal_output = fc.NIPYPE_LOGGING + reg_output = reg.run() + warped_file = reg_output.outputs.warped_file + + assert warped_file != '' + + return warped_file diff --git a/scripts/test-script b/scripts/test-script index 0715fd6..21b571c 100755 --- a/scripts/test-script +++ b/scripts/test-script @@ -59,7 +59,7 @@ for i in $FILES; do # cones 3D t2_star map if [ ! -e $CONES_T2STAR ]; then - python -m pipeline -l $SAVE_DIRNAME cones -t2star -fm $MASK + python -m pipeline --debug -l $SAVE_DIRNAME cones -t2star -fm $MASK fi # analyze femoral cartilage diff --git a/test_pipeline.py b/test_pipeline.py index f0d73e6..6520111 100644 --- a/test_pipeline.py +++ b/test_pipeline.py @@ -1,18 +1,17 @@ +import os import unittest + import keras.backend as K import numpy as np import scipy.io as sio -import os +import file_constants as fc import pipeline -from tissues.femoral_cartilage import FemoralCartilage from scan_sequences.cube_quant import CubeQuant - +from tissues.femoral_cartilage import FemoralCartilage from utils import io_utils, dicom_utils from utils.quant_vals import QuantitativeValues -import file_constants as fc - DESS_003_DICOM_PATH = './dicoms/003' DESS_003_T2_MAP_PATH = './dicoms/003_t2_map.mat' @@ -83,13 +82,13 @@ def test_t2_map(self): # need to convert all np.nan values to 0 before comparing # np.nan does not equal np.nan, so we need the same values to compare mat_t2_map = np.nan_to_num(mat_t2_map) - #py_t2_map = np.nan_to_num(py_t2_map) + # py_t2_map = np.nan_to_num(py_t2_map) # Round to the nearest 1000th (0.001) mat_t2_map = np.round(mat_t2_map, decimals=DECIMAL_PRECISION) py_t2_map = np.round(py_t2_map, decimals=DECIMAL_PRECISION) - assert((mat_t2_map == py_t2_map).all()) + assert ((mat_t2_map == py_t2_map).all()) def test_loading_data(self): vargin = self.get_vargin() @@ -101,13 +100,12 @@ def test_loading_data(self): pipeline.save_info(vargin[pipeline.SAVE_KEY], scan) - fc = FemoralCartilage() fc.load_data(vargin[pipeline.SAVE_KEY]) mask2 = fc.mask - assert((mask == mask2).all()) + assert ((mask == mask2).all()) def test_t2_map_load(self): vargin = self.get_vargin() @@ -245,7 +243,7 @@ def test_simpleitk_io(self): def test_h5(self): filepath = './dicoms/sample_h5_data' - datas = [{'type': np.zeros([4,4,4]), 'type2': np.zeros([4,4,4])}] + datas = [{'type': np.zeros([4, 4, 4]), 'type2': np.zeros([4, 4, 4])}] for data in datas: io_utils.save_h5(filepath, data) @@ -285,9 +283,10 @@ def test_femoral_cartilage(self): if not os.path.isdir(SAVE_PATH): self.segment_dess() - vargin = {'dicom': None, 'save': 'dicoms/healthy07/data', 'load': 'dicoms/healthy07/data', 'ext': 'dcm', 'gpu': None, + vargin = {'dicom': None, 'save': 'dicoms/healthy07/data', 'load': 'dicoms/healthy07/data', 'ext': 'dcm', + 'gpu': None, 'scan': 'tissues', 'fc': True, 't2': 'dicoms/healthy07/data/dess_data/t2.h5', 't1_rho': None, - 't2_star': None, 'tissues':[FemoralCartilage()]} + 't2_star': None, 'tissues': [FemoralCartilage()]} tissues = pipeline.handle_tissues(vargin) @@ -313,8 +312,5 @@ def test_reinit_variables(self): assert on == 'stream' - - if __name__ == '__main__': unittest.main() - diff --git a/tissues/femoral_cartilage.py b/tissues/femoral_cartilage.py index f593225..ddef50b 100644 --- a/tissues/femoral_cartilage.py +++ b/tissues/femoral_cartilage.py @@ -1,24 +1,23 @@ -import numpy as np import os -from tissues.tissue import Tissue - -from utils import io_utils -from utils.geometry_utils import circle_fit, cart2pol -import scipy.ndimage as sni -import pandas as pd - -import nipy.labs.mask as nlm +import warnings -from utils.quant_vals import QuantitativeValues import matplotlib.pyplot as plt +import nipy.labs.mask as nlm +import numpy as np +import pandas as pd +import scipy.ndimage as sni import defaults -import warnings +from tissues.tissue import Tissue +from utils import io_utils +from utils.geometry_utils import circle_fit, cart2pol +from utils.quant_vals import QuantitativeValues BOUNDS = {QuantitativeValues.T2: 100.0, QuantitativeValues.T1_RHO: 150.0, QuantitativeValues.T2_STAR: 100.0} + class FemoralCartilage(Tissue): ID = 1 STR_ID = 'fc' @@ -37,7 +36,7 @@ class FemoralCartilage(Tissue): LATERAL_KEY = 1 SAGGITAL_KEYS = [MEDIAL_KEY, LATERAL_KEY] - def __init__(self, weights_dir = None): + def __init__(self, weights_dir=None): super().__init__(weights_dir=weights_dir) self.regions_mask = None @@ -65,12 +64,12 @@ def unroll(self, qv_map): # ...considering the DEEP layers ### - mask = self.mask + mask = self.mask.volume - if (qv_map.shape != mask.shape): + if qv_map.shape != mask.shape: raise ValueError('t2_map and mask must have same shape') - if (len(qv_map.shape) != 3): + if len(qv_map.shape) != 3: raise ValueError('t2_map and mask must be 3D') num_slices = qv_map.shape[-1] @@ -162,7 +161,6 @@ def unroll(self, qv_map): def split_regions(self, unrolled_quantitative_map): # WARNING: this method has to be called after the unroll method. - import matplotlib.pyplot as plt # create unrolled mask from unrolled map unrolled_mask_indexes = np.nonzero(unrolled_quantitative_map) @@ -212,10 +210,13 @@ def calc_quant_vals(self, quant_map, map_type): :param map_type: A QuantitativeValue instance :return: """ + + super().calc_quant_vals(quant_map, map_type) + if self.mask is None: raise ValueError('Please initialize mask') - total, superficial, deep = self.unroll(quant_map) + total, superficial, deep = self.unroll(quant_map.volume) assert total.shape == deep.shape assert deep.shape == superficial.shape @@ -238,7 +239,7 @@ def calc_quant_vals(self, quant_map, map_type): for sagital in [self.ANTERIOR_KEY, self.CENTRAL_KEY, self.POSTERIOR_KEY]: curr_region_mask = (coronal_region_mask == coronal) * (sagital_region_mask == sagital) * axial_map - curr_region_mask[curr_region_mask==0] = np.nan + curr_region_mask[curr_region_mask == 0] = np.nan # discard all values that are 0 c_mean = np.nanmean(curr_region_mask) c_std = np.nanstd(curr_region_mask) @@ -250,19 +251,24 @@ def calc_quant_vals(self, quant_map, map_type): depth_keys = np.array(['deep', 'deep', 'superficial', 'superficial', 'total', 'total']) coronal_keys = np.array(['medial', 'lateral'] * 3) sagital_keys = ['anterior', 'central', 'posterior'] - df = pd.DataFrame(data=np.transpose(tissue_values), index=sagital_keys, columns=pd.MultiIndex.from_tuples(zip(depth_keys, coronal_keys))) + df = pd.DataFrame(data=np.transpose(tissue_values), index=sagital_keys, + columns=pd.MultiIndex.from_tuples(zip(depth_keys, coronal_keys))) qv_name = map_type.name - maps = [{'title': '%s deep' % qv_name, 'data': deep, 'xlabel': 'Slice', 'ylabel': 'Angle (binned)', 'filename': '%s_deep.png' % qv_name}, - {'title': '%s superficial' % qv_name, 'data': superficial, 'xlabel': 'Slice', 'ylabel': 'Angle (binned)', 'filename': '%s_superficial.png' % qv_name}, - {'title': '%s total' % qv_name, 'data': total, 'xlabel': 'Slice', 'ylabel': 'Angle (binned)', 'filename': '%s_total.png' % qv_name}] + maps = [{'title': '%s deep' % qv_name, 'data': deep, 'xlabel': 'Slice', 'ylabel': 'Angle (binned)', + 'filename': '%s_deep.png' % qv_name}, + {'title': '%s superficial' % qv_name, 'data': superficial, 'xlabel': 'Slice', + 'ylabel': 'Angle (binned)', 'filename': '%s_superficial.png' % qv_name}, + {'title': '%s total' % qv_name, 'data': total, 'xlabel': 'Slice', 'ylabel': 'Angle (binned)', + 'filename': '%s_total.png' % qv_name}] self.__store_quant_vals__(maps, df, map_type) - def set_mask(self, mask, pixel_spacing): - mask = np.asarray(nlm.largest_cc(mask), dtype=np.uint8) + def set_mask(self, mask): + msk = np.asarray(nlm.largest_cc(mask.volume), dtype=np.uint8) + mask.volume = msk self.regions_mask = None - super().set_mask(mask, pixel_spacing) + super().set_mask(mask) def __save_quant_data__(self, dirpath): q_names = [] diff --git a/tissues/meniscus.py b/tissues/meniscus.py index 8977abc..cfa857f 100644 --- a/tissues/meniscus.py +++ b/tissues/meniscus.py @@ -1,6 +1,6 @@ -import numpy as np from tissues.tissue import Tissue + class Meniscus(Tissue): ID = 4 STR_ID = 'men' @@ -13,11 +13,10 @@ class Meniscus(Tissue): SUPERFICIAL_CENTRAL_REGION_KEY = 'superficial_central' SUPERFICIAL_LATERAL_REGION_KEY = 'superficial_lateral' - def split_regions(self, mask): - #TODO: implement spliting regions + # TODO: implement spliting regions pass def calc_quant_vals(self, quant_map, mask=None): # TODO: implement getting quantitative values for regions regions - pass \ No newline at end of file + pass diff --git a/tissues/patellar_cartilage.py b/tissues/patellar_cartilage.py index 3d54e80..f91febe 100644 --- a/tissues/patellar_cartilage.py +++ b/tissues/patellar_cartilage.py @@ -1,4 +1,3 @@ -import numpy as np from tissues.tissue import Tissue @@ -7,9 +6,9 @@ class PatellarCartilage(Tissue): STR_ID = 'pc' def split_regions(self, mask): - #TODO: implement spliting regions + # TODO: implement spliting regions pass def calc_quant_vals(self, quant_map, mask=None): # TODO: implement getting quantitative values for regions regions - pass \ No newline at end of file + pass diff --git a/tissues/tibial_cartilage.py b/tissues/tibial_cartilage.py index 5eadbbf..6b7ad3e 100644 --- a/tissues/tibial_cartilage.py +++ b/tissues/tibial_cartilage.py @@ -1,4 +1,3 @@ -import numpy as np from tissues.tissue import Tissue @@ -8,11 +7,10 @@ class TibialCartilage(Tissue): # TODO: implement region keys - def split_regions(self, mask): - #TODO: implement spliting regions + # TODO: implement spliting regions pass def calc_quant_vals(self, quant_map, mask=None): # TODO: implement getting quantitative values for regions regions - pass \ No newline at end of file + pass diff --git a/tissues/tissue.py b/tissues/tissue.py index 42341f1..0cd6027 100644 --- a/tissues/tissue.py +++ b/tissues/tissue.py @@ -1,25 +1,26 @@ -from abc import ABC, abstractmethod import os +from abc import ABC, abstractmethod + +from med_objects.med_volume import MedicalVolume from utils import io_utils from utils.quant_vals import QuantitativeValues -import cv2 -import numpy as np WEIGHTS_FILE_EXT = 'h5' + class Tissue(ABC): ID = -1 # should be unique to all tissues, and should not change STR_ID = '' FULL_NAME = '' ORIENTATION = '' - def __init__(self, weights_dir = None): + def __init__(self, weights_dir=None): self.regions = dict() self.mask = None self.quant_vals = dict() self.weights_filepath = None - self.pixel_spacing = None - if (weights_dir is not None): + + if weights_dir is not None: self.weights_filepath = self.find_weights(weights_dir) @abstractmethod @@ -39,6 +40,10 @@ def calc_quant_vals(self, quant_map, map_type): :param map_type: an enum instance of QuantitativeValue :return: a dictionary of quantitative values, save in quant_vals """ + + assert type(quant_map) is MedicalVolume + assert type(map_type) is QuantitativeValues + pass def __store_quant_vals__(self, quant_map, quant_df, map_type): @@ -73,7 +78,7 @@ def save_data(self, dirpath): if self.mask is not None: mask_filepath = os.path.join(dirpath, '%s.nii.gz' % self.STR_ID) - io_utils.save_nifti(mask_filepath, self.mask, self.pixel_spacing) + self.mask.save_volume(mask_filepath) self.__save_quant_data__(dirpath) @@ -88,11 +93,11 @@ def load_data(self, dirpath): raise FileNotFoundError('File \'%s\' does not exist' % mask_filepath) filepath = os.path.join(dirpath, '%s.nii.gz' % self.STR_ID) - self.mask, self.pixel_spacing = io_utils.load_nifti(filepath) + self.mask = io_utils.load_nifti(filepath) def __save_dirpath__(self, dirpath): return io_utils.check_dir(os.path.join(dirpath, '%s' % self.STR_ID)) - def set_mask(self, mask, pixel_spacing): + def set_mask(self, mask): + assert (type(mask) is MedicalVolume) self.mask = mask - self.pixel_spacing = pixel_spacing diff --git a/utils/dicom_utils.py b/utils/dicom_utils.py index 7f441b3..ce73e35 100644 --- a/utils/dicom_utils.py +++ b/utils/dicom_utils.py @@ -1,10 +1,10 @@ import os + import numpy as np -from natsort import natsorted import pydicom -import math -from pydicom.tag import Tag +from natsort import natsorted +from med_objects.med_volume import MedicalVolume __VOLUME_DIMENSIONS__ = 3 __EPSILON__ = 1e-8 @@ -59,7 +59,9 @@ def load_dicom(dicom_path, dicom_ext=None): refs_dicom.append(ds) dicom_array[:, :, lstFilesDCM.index(dicom_filename)] = ds.pixel_array - return dicom_array, refs_dicom, pixelSpacing + volume = MedicalVolume(dicom_array, pixelSpacing) + + return volume, refs_dicom def whiten_volume(x): @@ -71,4 +73,4 @@ def whiten_volume(x): raise ValueError("Dimension Error: input has %d dimensions. Expected %d" % (x.ndims, __VOLUME_DIMENSIONS__)) # Add epsilon to avoid dividing by 0 - return (x - np.mean(x)) / (np.std(x) + __EPSILON__) \ No newline at end of file + return (x - np.mean(x)) / (np.std(x) + __EPSILON__) diff --git a/utils/fits.py b/utils/fits.py new file mode 100644 index 0000000..dde21d9 --- /dev/null +++ b/utils/fits.py @@ -0,0 +1,129 @@ +import warnings +from abc import ABC, abstractmethod + +import numpy as np +from scipy import optimize as sop + +import defaults +from med_objects.med_volume import MedicalVolume + + +class Fit(ABC): + @abstractmethod + def fit(self): + pass + + +class MonoExponentialFit(Fit): + def __init__(self, ts, subvolumes, mask=None, bounds=(0, 100.0), tc0=30.0, decimal_precision=1): + assert (len(ts) == len(subvolumes)) + self.ts = ts + + assert (type(subvolumes) is list) + for sv in subvolumes: + assert type(sv) is MedicalVolume + + self.subvolumes = subvolumes + + if mask is not None: + assert (type(mask) is MedicalVolume) + self.mask = mask + + assert (type(bounds) is tuple and len(bounds) == 2) + self.bounds = bounds + + self.tc0 = tc0 + self.decimal_precision = decimal_precision + + def fit(self): + original_shape = None + svs = [] + msk = None + if self.mask: + msk = self.mask.volume + msk = msk.reshape(1, -1) + + pixel_spacing = self.subvolumes[0].pixel_spacing + + for i in range(len(self.ts)): + sv = self.subvolumes[i].volume + + if original_shape is None: + original_shape = sv.shape + else: + assert (sv.shape == original_shape) + + svr = sv.reshape((1, -1)) + if msk is not None: + svr = svr * msk + + svs.append(svr) + + svs = np.concatenate(svs) + + vals, r_squared = fit_monoexp_tc(self.ts, svs, self.tc0) + + map_unfiltered = vals.reshape(original_shape) + r_squared = r_squared.reshape(original_shape) + + tc_map = map_unfiltered * (r_squared > defaults.DEFAULT_R2_THRESHOLD) + + # Filter calculated T1-rho values that are below 0ms and over 100ms + tc_map[tc_map <= self.bounds[0]] = np.nan + tc_map = np.nan_to_num(tc_map) + tc_map[tc_map > self.bounds[1]] = np.nan + tc_map = np.nan_to_num(tc_map) + + tc_map = np.around(tc_map, self.decimal_precision) + + return MedicalVolume(tc_map, pixel_spacing), MedicalVolume(r_squared, pixel_spacing) + + +__EPSILON__ = 1e-8 + + +def __fit_mono_exp__(x, y, p0=None): + def func(t, a, b): + exp = np.exp(b * t) + return a * exp + + x = np.asarray(x) + y = np.asarray(y) + + popt, _ = sop.curve_fit(func, x, y, p0=p0, maxfev=1000) + + residuals = y - func(x, popt[0], popt[1]) + ss_res = np.sum(residuals ** 2) + ss_tot = np.sum((y - np.mean(y)) ** 2) + + r_squared = 1 - (ss_res / (ss_tot + __EPSILON__)) + + return popt, r_squared + + +def fit_monoexp_tc(x, ys, tc0): + p0 = (1.0, -1 / tc0) + time_constants = np.zeros([1, ys.shape[-1]]) + r_squared = np.zeros([1, ys.shape[-1]]) + + warned_negative = False + for i in range(ys.shape[1]): + y = ys[..., i] + if (y < 0).any() and not warned_negative: + warned_negative = True + warnings.warn("Negative values found. Failure in monoexponential fit will result in t1_rho=np.nan") + + # Skip any negative values or all values that are 0s + if (y < 0).any() or (y == 0).all(): + continue + + try: + params, r2 = __fit_mono_exp__(x, y, p0=p0) + tc = 1 / abs(params[-1]) + except RuntimeError: + tc, r2 = (np.nan, 0.0) + + time_constants[..., i] = tc + r_squared[..., i] = r2 + + return time_constants, r_squared diff --git a/utils/im_utils.py b/utils/im_utils.py index a9d357a..e427b95 100644 --- a/utils/im_utils.py +++ b/utils/im_utils.py @@ -1,7 +1,9 @@ -import cv2 import os -import numpy as np + import SimpleITK as sitk +import cv2 +import numpy as np + from utils import io_utils TIFF_EXTENSION = '.tiff' @@ -45,4 +47,4 @@ def write_3d(filepath, mask): visual_mask = np.asarray(mask * 255, np.uint8) visual_mask = np.transpose(visual_mask, permute_order) img = sitk.GetImageFromArray(visual_mask) - sitk.WriteImage(img, filepath) \ No newline at end of file + sitk.WriteImage(img, filepath) diff --git a/utils/io_utils.py b/utils/io_utils.py index 2970745..664c029 100644 --- a/utils/io_utils.py +++ b/utils/io_utils.py @@ -1,12 +1,11 @@ -import h5py -import pickle import os -import pandas as pd -import SimpleITK as sitk - +import pickle import warnings -import numpy as np +import SimpleITK as sitk +import h5py +import numpy as np +import pandas as pd DATA_EXT = 'data' INFO_EXT = 'info' @@ -24,7 +23,7 @@ def get_subdirs(dirpath): if not os.path.isdir(dirpath): raise NotADirectoryError('%s not a directory' % dirpath) - subdirs =[] + subdirs = [] for file in os.listdir(dirpath): possible_dir = os.path.join(dirpath, file) if os.path.isdir(possible_dir): @@ -96,7 +95,7 @@ def save_tables(filepath, data_frames, sheet_names=None): if sheet_names is None: sheet_names = [] for i in range(len(data_frames)): - sheet_names.append('Sheet%d' % (i+1)) + sheet_names.append('Sheet%d' % (i + 1)) if len(data_frames) != len(sheet_names): raise ValueError('Number of data_frames and sheet_frames should be the same') @@ -109,7 +108,6 @@ def save_tables(filepath, data_frames, sheet_names=None): def save_nifti(filepath, img_array, spacing): - assert filepath.endswith('.nii.gz') if img_array is None or len(img_array.shape) < 2: warnings.warn('%s not saved. Input array is None' % img_array) @@ -128,6 +126,8 @@ def save_nifti(filepath, img_array, spacing): def load_nifti(filepath): + from med_objects.med_volume import MedicalVolume + assert filepath.endswith('.nii.gz') image = sitk.ReadImage(filepath) spacing = image.GetSpacing() @@ -137,8 +137,4 @@ def load_nifti(filepath): # invert array for convention of SimpleITK - array is now in form (row, column, depth) img_array = np.transpose(img_array, [1, 2, 0]) - return img_array, spacing - - - - + return MedicalVolume(img_array, spacing) diff --git a/utils/quant_vals.py b/utils/quant_vals.py index d7bee80..11b0f63 100644 --- a/utils/quant_vals.py +++ b/utils/quant_vals.py @@ -1,12 +1,4 @@ from enum import Enum -import scipy.optimize as sop -import numpy as np -import glob -import warnings - - -__EPSILON__ = 1e-8 -__R_SQUARED_THRESHOLD__ = 0.9 class QuantitativeValues(Enum): @@ -21,53 +13,6 @@ def get_qv(id): return qv -def __fit_mono_exp__(x, y, p0=None): - def func(t, a, b): - exp = np.exp(b * t) - return a * exp - - x = np.asarray(x) - y = np.asarray(y) - - popt, _ = sop.curve_fit(func, x, y, p0=p0, maxfev=1000) - - residuals = y - func(x, popt[0], popt[1]) - ss_res = np.sum(residuals ** 2) - ss_tot = np.sum((y - np.mean(y)) ** 2) - - r_squared = 1 - (ss_res / (ss_tot + __EPSILON__)) - - return popt, r_squared - - -def fit_monoexp_tc(x, ys, tc0): - p0 = (1.0, -1/tc0) - time_constants = np.zeros([1, ys.shape[-1]]) - r_squared = np.zeros([1, ys.shape[-1]]) - - warned_negative = False - for i in range(ys.shape[1]): - y = ys[..., i] - if (y < 0).any() and not warned_negative: - warned_negative = True - warnings.warn("Negative values found. Failure in monoexponential fit will result in t1_rho=np.nan") - - # Skip any negative values or all values that are 0s - if (y < 0).any() or (y == 0).all(): - continue - - try: - params, r2 = __fit_mono_exp__(x, y, p0=p0) - tc = 1 / abs(params[-1]) - except RuntimeError: - tc, r2 = (np.nan, 0.0) - - time_constants[..., i] = tc - r_squared[..., i] = r2 - - return time_constants, r_squared - - if __name__ == '__main__': print(type(QuantitativeValues.T1_RHO.name)) - print(QuantitativeValues.T1_RHO.value == 1) \ No newline at end of file + print(QuantitativeValues.T1_RHO.value == 1)