From 5bf2df23ef2c894d5d709855d94c905b5d021ed5 Mon Sep 17 00:00:00 2001 From: sophie22 Date: Mon, 24 Jul 2023 19:32:31 +0100 Subject: [PATCH 01/17] minor changes to docs --- hazenlib/relaxometry.py | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/hazenlib/relaxometry.py b/hazenlib/relaxometry.py index 7ea3e30c..f40ce335 100644 --- a/hazenlib/relaxometry.py +++ b/hazenlib/relaxometry.py @@ -6,9 +6,9 @@ This module determines the T1 and T2 decay constants for the relaxometry spheres in the Caliber (HPD) system phantom -qmri.com/qmri-solutions/t1-t2-pd-imaging-phantom (plates 4 and 5). Values are -compared to published values (without temperature correction). Graphs of fit -and phantom registration images can optionally be produced. +qmri.com/qmri-solutions/t1-t2-pd-imaging-phantom (plates 4 and 5). +Values are compared to published values (without temperature correction). +Graphs of fit and phantom registration images can optionally be produced. Scan parameters @@ -81,10 +81,9 @@ each relaxation parameter (T1 or T2) on plates 4 and 5, and regression is performed on the first image in the sequence. Optionally output the overlay image to visually check the fit. -3. An ROI is generated for each target sphere using stored coordinates, the - RT transformation above, and a structuring element (default is a 5x5 - boxcar). -4. Store pixel data for each ROI, at various times, in an ``ROITimeSeries`` +3. A ROI is generated for each target sphere using stored coordinates, the RT + transformation above, and a structuring element (default is a 5x5 boxcar). +4. Store pixel data for each ROI at various times, in an ``ROITimeSeries`` object. A list of these objects is stored in ``ImageStack.ROI_time_series``. 5. Generate the fit function. For T1 this looks up TR for the given TI @@ -93,7 +92,7 @@ measurements. 6. Determine relaxation time (T1 or T2) by fitting the decay equation to the ROI data for each sphere. The published values of the relaxation - times are used to seed the optimisation algorithm. A Rician nose model is + times are used to seed the optimisation algorithm. A Rician noise model is used for T2 fitting [1]_. Optionally plot and save the decay curves. 7. Return plate number, relaxation type (T1 or T2), measured relaxation times, published relaxation times, and fractional differences in a @@ -113,8 +112,8 @@ -have bolthole template, find 3 positions in template and image, figure out transformation. -Template fit on outline image--poss run though edge detection algorithms then -fit. +Template fit on outline image--possibly run though edge detection algorithms +then fit. Use normalised structuring element in ROITimeSeries. This will allow correct calculation of mean if elements are not 0 or 1. @@ -122,9 +121,9 @@ Get r-squared measure of fit. """ -import os.path import sys import json +import os.path import pathlib import cv2 as cv From 5069277730eafa5fc26f5faf665ce60acc7d6a8e Mon Sep 17 00:00:00 2001 From: sophie22 Date: Mon, 24 Jul 2023 19:41:20 +0100 Subject: [PATCH 02/17] move parameters to config file --- hazenlib/config.py | 74 +++++++++++++++++++++++++++++++++++++++++ hazenlib/relaxometry.py | 73 +--------------------------------------- 2 files changed, 75 insertions(+), 72 deletions(-) create mode 100644 hazenlib/config.py diff --git a/hazenlib/config.py b/hazenlib/config.py new file mode 100644 index 00000000..1eb1d5c3 --- /dev/null +++ b/hazenlib/config.py @@ -0,0 +1,74 @@ +import os.path +import numpy as np + +# Parameters for Rician noise model +MAX_RICIAN_NOISE = 20.0 +SEED_RICIAN_NOISE = 5.0 + +TEMPLATE_DIR = os.path.join(os.path.dirname(os.path.realpath(__file__)), + 'data', 'relaxometry') +TEMPLATE_VALUES = { + 'plate3': { + 'sphere_centres_row_col': (), + 'bolt_centres_row_col': (), + 't1': { + 'filename': '', + 'relax_times': []}}, + + 'plate4': { + 'sphere_centres_row_col': ( + (56, 94), (62, 117), (81, 132), (105, 134), (125, 120), (133, 99), + (127, 75), (108, 60), (84, 59), (64, 72), (80, 81), (78, 111), + (109, 113), (111, 82), (148, 118)), + 'bolt_centres_row_col': (), + 't1': { + 'filename': os.path.join(TEMPLATE_DIR, 'Plate4_T1_signed'), + 'relax_times': { + '1.5T': + np.array([2376, 2183, 1870, 1539, 1237, 1030, 752.2, 550.2, + 413.4, 292.9, 194.9, 160.2, 106.4, 83.3, 2700]), + '3.0T': + np.array([2480, 2173, 1907, 1604, 1332, 1044, 801.7, 608.6, + 458.4, 336.5, 244.2, 176.6, 126.9, 90.9, 2700])}}, + 't2': { + 'filename': os.path.join(TEMPLATE_DIR, 'Plate4_T2'), + 'relax_times': { + '1.5T': + np.array([939.4, 594.3, 416.5, 267.0, 184.9, 140.6, 91.76, + 64.84, 45.28, 30.62, 19.76, 15.99, 10.47, 8.15, + 2400]), + '3.0T': + np.array([581.3, 403.5, 278.1, 190.94, 133.27, 96.89, + 64.07, 46.42, 31.97, 22.56, 15.813, 11.237, + 7.911, 5.592, 2400])}}}, + + 'plate5': { + 'sphere_centres_row_col': ( + (56, 95), (62, 117), (81, 133), (104, 134), (124, 121), (133, 98), + (127, 75), (109, 61), (84, 60), (64, 72), (80, 81), (78, 111), + (109, 113), (110, 82), (97, 43)), + 'bolt_centres_row_col': ((52, 80), (92, 141), (138, 85)), + 't1': { + 'filename': os.path.join(TEMPLATE_DIR, 'Plate5_T1_signed'), + 'relax_times': { + '1.5T': + np.array([2033, 1489, 1012, 730.8, 514.1, 367.9, 260.1, + 184.6, 132.7, 92.7, 65.4, 46.32, 32.45, 22.859, + 2700]), + '3.0T': + np.array([1989, 1454, 984.1, 706, 496.7, 351.5, 247.13, + 175.3, 125.9, 89.0, 62.7, 44.53, 30.84, + 21.719, 2700])}}, + + 't2': { + 'filename': os.path.join(TEMPLATE_DIR, 'Plate5_T2'), + 'relax_times': { + '1.5T': + np.array([1669.0, 1244.0, 859.3, 628.5, 446.3, 321.2, + 227.7, 161.9, 117.1, 81.9, 57.7, 41.0, 28.7, + 20.2, 2400]), + '3.0T': + np.array([1465, 1076, 717.9, 510.1, 359.6, 255.5, 180.8, + 127.3, 90.3, 64.3, 45.7, 31.86, 22.38, + 15.83, 2400])}}}} + diff --git a/hazenlib/relaxometry.py b/hazenlib/relaxometry.py index f40ce335..faab8807 100644 --- a/hazenlib/relaxometry.py +++ b/hazenlib/relaxometry.py @@ -137,10 +137,7 @@ from scipy.special import i0e, ive import hazenlib.exceptions - -# Parameters for Rician noise model -MAX_RICIAN_NOISE = 20.0 -SEED_RICIAN_NOISE = 5.0 +from hazenlib.config import MAX_RICIAN_NOISE, SEED_RICIAN_NOISE, TEMPLATE_VALUES # Use dict to store template and reference information # Coordinates are in array format (row,col), rather than plt.patches @@ -151,74 +148,6 @@ # TEMPLATE_VALUES[f'plate{plate_num}']['t1'|'t2']['filename'] # TEMPLATE_VALUES[f'plate{plate_num}']['t1'|'t2']['1.5T'|'3.0T']['relax_times'] -TEMPLATE_DIR = os.path.join(os.path.dirname(os.path.realpath(__file__)), - 'data', 'relaxometry') -TEMPLATE_VALUES = { - 'plate3': { - 'sphere_centres_row_col': (), - 'bolt_centres_row_col': (), - 't1': { - 'filename': '', - 'relax_times': []}}, - - 'plate4': { - 'sphere_centres_row_col': ( - (56, 94), (62, 117), (81, 132), (105, 134), (125, 120), (133, 99), - (127, 75), (108, 60), (84, 59), (64, 72), (80, 81), (78, 111), - (109, 113), (111, 82), (148, 118)), - 'bolt_centres_row_col': (), - 't1': { - 'filename': os.path.join(TEMPLATE_DIR, 'Plate4_T1_signed'), - 'relax_times': { - '1.5T': - np.array([2376, 2183, 1870, 1539, 1237, 1030, 752.2, 550.2, - 413.4, 292.9, 194.9, 160.2, 106.4, 83.3, 2700]), - '3.0T': - np.array([2480, 2173, 1907, 1604, 1332, 1044, 801.7, 608.6, - 458.4, 336.5, 244.2, 176.6, 126.9, 90.9, 2700])}}, - 't2': { - 'filename': os.path.join(TEMPLATE_DIR, 'Plate4_T2'), - 'relax_times': { - '1.5T': - np.array([939.4, 594.3, 416.5, 267.0, 184.9, 140.6, 91.76, - 64.84, 45.28, 30.62, 19.76, 15.99, 10.47, 8.15, - 2400]), - '3.0T': - np.array([581.3, 403.5, 278.1, 190.94, 133.27, 96.89, - 64.07, 46.42, 31.97, 22.56, 15.813, 11.237, - 7.911, 5.592, 2400])}}}, - - 'plate5': { - 'sphere_centres_row_col': ( - (56, 95), (62, 117), (81, 133), (104, 134), (124, 121), (133, 98), - (127, 75), (109, 61), (84, 60), (64, 72), (80, 81), (78, 111), - (109, 113), (110, 82), (97, 43)), - 'bolt_centres_row_col': ((52, 80), (92, 141), (138, 85)), - 't1': { - 'filename': os.path.join(TEMPLATE_DIR, 'Plate5_T1_signed'), - 'relax_times': { - '1.5T': - np.array([2033, 1489, 1012, 730.8, 514.1, 367.9, 260.1, - 184.6, 132.7, 92.7, 65.4, 46.32, 32.45, 22.859, - 2700]), - '3.0T': - np.array([1989, 1454, 984.1, 706, 496.7, 351.5, 247.13, - 175.3, 125.9, 89.0, 62.7, 44.53, 30.84, - 21.719, 2700])}}, - - 't2': { - 'filename': os.path.join(TEMPLATE_DIR, 'Plate5_T2'), - 'relax_times': { - '1.5T': - np.array([1669.0, 1244.0, 859.3, 628.5, 446.3, 321.2, - 227.7, 161.9, 117.1, 81.9, 57.7, 41.0, 28.7, - 20.2, 2400]), - '3.0T': - np.array([1465, 1076, 717.9, 510.1, 359.6, 255.5, 180.8, - 127.3, 90.3, 64.3, 45.7, 31.86, 22.38, - 15.83, 2400])}}}} - - def outline_mask(im): """ Create contour lines to outline pixels. From 89c1235441b2f78f1b450dfec14999f6e346b3db Mon Sep 17 00:00:00 2001 From: sophie22 Date: Mon, 24 Jul 2023 20:06:27 +0100 Subject: [PATCH 03/17] move smooth_times to config, simplify relax string --- hazenlib/config.py | 4 ++++ hazenlib/relaxometry.py | 28 ++++++++++++---------------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/hazenlib/config.py b/hazenlib/config.py index 1eb1d5c3..f09afec2 100644 --- a/hazenlib/config.py +++ b/hazenlib/config.py @@ -72,3 +72,7 @@ 127.3, 90.3, 64.3, 45.7, 31.86, 22.38, 15.83, 2400])}}}} +SMOOTH_TIMES = { + "t1": range(0, 1000, 10), + "t2": range(0, 500, 5) +} diff --git a/hazenlib/relaxometry.py b/hazenlib/relaxometry.py index faab8807..e586cfb5 100644 --- a/hazenlib/relaxometry.py +++ b/hazenlib/relaxometry.py @@ -121,7 +121,6 @@ Get r-squared measure of fit. """ -import sys import json import os.path import pathlib @@ -137,7 +136,7 @@ from scipy.special import i0e, ive import hazenlib.exceptions -from hazenlib.config import MAX_RICIAN_NOISE, SEED_RICIAN_NOISE, TEMPLATE_VALUES +from hazenlib.config import MAX_RICIAN_NOISE, SEED_RICIAN_NOISE, TEMPLATE_VALUES, SMOOTH_TIMES # Use dict to store template and reference information # Coordinates are in array format (row,col), rather than plt.patches @@ -469,7 +468,6 @@ def est_t2_s0(te, t2, pv, c=0.0): return (pv - c) / np.exp(-te / t2) - class ROITimeSeries: """ Samples at one image location (ROI) at numerous sample times. @@ -1015,9 +1013,8 @@ def relax_times(self): return self.t2s -def main(dcm_target_list, plate_number=None, - calc: str = 'T1', verbose=False, - report=False, report_dir=None): +def main(data, plate_number=None, calc: str = 'T1', verbose=False, + report=False, report_dir=None): """ Calculate T1 or T2 values for relaxometry phantom. @@ -1025,7 +1022,7 @@ def main(dcm_target_list, plate_number=None, Parameters ---------- - dcm_target_list : list of pydicom.dataset.FileDataSet objects + data : list of pydicom.dataset.FileDataSet objects List of DICOM images of a plate of the HPD relaxometry phantom. calc : str, required Whether to calculate T1 or T2 relaxation. Default is T1. @@ -1073,8 +1070,7 @@ def main(dcm_target_list, plate_number=None, # Set up parameters specific to T1 or T2 if calc in ['T1', 't1']: ImStack = T1ImageStack - relax_str = 't1' - smooth_times = range(0, 1000, 10) + relax_str = calc.lower() try: template_dcm = pydicom.read_file( TEMPLATE_VALUES[f'plate{plate_number}'][relax_str]['filename']) @@ -1085,8 +1081,7 @@ def main(dcm_target_list, plate_number=None, elif calc in ['T2', 't2']: ImStack = T2ImageStack - relax_str = 't2' - smooth_times = range(0, 500, 5) + relax_str = calc.lower() try: template_dcm = pydicom.read_file( TEMPLATE_VALUES[f'plate{plate_number}'][relax_str]['filename']) @@ -1096,9 +1091,9 @@ def main(dcm_target_list, plate_number=None, exit() else: print("Please provide 'T1' or 'T2' for the --calc argument.") - sys.exit() + exit() - image_stack = ImStack(dcm_target_list, template_dcm) + image_stack = ImStack(data, template_dcm) image_stack.template_fit() image_stack.generate_time_series( TEMPLATE_VALUES[f'plate{plate_number}'] @@ -1142,7 +1137,7 @@ def main(dcm_target_list, plate_number=None, # Show ROIs roi_fig = image_stack.plot_rois() - plt.title(f'ROI positions ({relax_str.upper()}, plate {plate_number})') + plt.title(f'ROI positions ({calc.upper()}, plate {plate_number})') roi_img = f'{img_path}_rois.png' roi_fig.savefig(roi_img, dpi=300) report_files['rois'] = roi_img @@ -1150,7 +1145,8 @@ def main(dcm_target_list, plate_number=None, # Show relax fits rois = image_stack.ROI_time_series relax_fit_fig = plt.figure() - relax_fit_fig.suptitle(relax_str.upper() + ' relaxometry fits') + relax_fit_fig.suptitle(calc.upper() + ' relaxometry fits') + smooth_times = SMOOTH_TIMES[relax_str] for i in range(15): plt.subplot(5, 3, i + 1) plt.plot(smooth_times, @@ -1182,7 +1178,7 @@ def main(dcm_target_list, plate_number=None, metadata = dict( files=[im.filename for im in image_stack.images], plate=plate_number, - relaxation_type=relax_str, + relaxation_type=calc.upper(), institution_name=index_im.InstitutionName, manufacturer=index_im.Manufacturer, model=index_im.ManufacturerModelName, From bfc455d44f71a83f58772953efeee263d9c067bd Mon Sep 17 00:00:00 2001 From: sophie22 Date: Mon, 24 Jul 2023 20:08:00 +0100 Subject: [PATCH 04/17] make image_stack object creation more explicit --- hazenlib/relaxometry.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/hazenlib/relaxometry.py b/hazenlib/relaxometry.py index e586cfb5..c1742379 100644 --- a/hazenlib/relaxometry.py +++ b/hazenlib/relaxometry.py @@ -1069,22 +1069,22 @@ def main(data, plate_number=None, calc: str = 'T1', verbose=False, # Set up parameters specific to T1 or T2 if calc in ['T1', 't1']: - ImStack = T1ImageStack relax_str = calc.lower() try: template_dcm = pydicom.read_file( TEMPLATE_VALUES[f'plate{plate_number}'][relax_str]['filename']) + image_stack = T1ImageStack(data, template_dcm) except KeyError: print(f'Could not find template with plate number: {plate_number}.' f' Please pass plate number as arg.') exit() elif calc in ['T2', 't2']: - ImStack = T2ImageStack relax_str = calc.lower() try: template_dcm = pydicom.read_file( TEMPLATE_VALUES[f'plate{plate_number}'][relax_str]['filename']) + image_stack = T2ImageStack(data, template_dcm) except KeyError: print(f'Could not find template with plate number: {plate_number}.' f' Please pass plate number as arg.') @@ -1093,7 +1093,6 @@ def main(data, plate_number=None, calc: str = 'T1', verbose=False, print("Please provide 'T1' or 'T2' for the --calc argument.") exit() - image_stack = ImStack(data, template_dcm) image_stack.template_fit() image_stack.generate_time_series( TEMPLATE_VALUES[f'plate{plate_number}'] @@ -1146,7 +1145,7 @@ def main(data, plate_number=None, calc: str = 'T1', verbose=False, rois = image_stack.ROI_time_series relax_fit_fig = plt.figure() relax_fit_fig.suptitle(calc.upper() + ' relaxometry fits') - smooth_times = SMOOTH_TIMES[relax_str] + smooth_times = SMOOTH_TIMES[calc.lower()] for i in range(15): plt.subplot(5, 3, i + 1) plt.plot(smooth_times, From b9666cfb3fbb3913e5e00bb3a55fb1d11baae19f Mon Sep 17 00:00:00 2001 From: sophie22 Date: Tue, 25 Jul 2023 09:35:35 +0100 Subject: [PATCH 05/17] rename dicom_order_key to time_attribute, update order_by function --- hazenlib/relaxometry.py | 54 ++++++++++++++++++++--------------------- 1 file changed, 26 insertions(+), 28 deletions(-) diff --git a/hazenlib/relaxometry.py b/hazenlib/relaxometry.py index c1742379..4579bdb5 100644 --- a/hazenlib/relaxometry.py +++ b/hazenlib/relaxometry.py @@ -538,7 +538,6 @@ def __init__(self, dcm_images, poi_coords_row_col, time_attr=None, self.ROI_mask = np.zeros_like(self.POI_mask) self.ROI_mask = scipy.ndimage.filters.convolve(self.POI_mask, kernel) - self._time_attr = time_attr if time_attr is not None: self.times = [x[time_attr].value.real for x in dcm_images] @@ -569,7 +568,7 @@ class ImageStack(): """ def __init__(self, image_slices, template_dcm, - dicom_order_key=None): + time_attribute=None): """ Create ImageStack object. @@ -585,20 +584,21 @@ def __init__(self, image_slices, template_dcm, For future use. Reference to the plate in the relaxometry phantom. The default is None. - dicom_order_key : string, optional + time_attribute : string, optional DICOM attribute to order images. Typically 'InversionTime' for T1 relaxometry or 'EchoTime' for T2. """ - # Store template pixel array, after scaling in 0028,1052 and 0028,1053 - # applied + # Store template pixel array, after scaling in 0028,1052 and + # 0028,1053 applied self.template_dcm = template_dcm if template_dcm is not None: self.template_px = pixel_rescale(template_dcm) - self.dicom_order_key = dicom_order_key - self.images = image_slices # store images - if dicom_order_key is not None: - self.order_by(dicom_order_key) + # store sorted images + if time_attribute is not None: + self.images = self.order_by(image_slices, time_attribute) + else: + self.images = image_slices b0_val = self.images[0]['MagneticFieldStrength'].value if b0_val == 1.5: @@ -656,32 +656,31 @@ def template_fit(self, image_index=0): angle rotations. """ target_px = pixel_rescale(self.images[0]) - template_px = self.template_px - # Pad template or target pixels if required - scale_factor = len(target_px) / len(template_px) - pad_size = np.subtract(template_px.shape, target_px.shape) + ## Pad template or target pixels if required + # Determine difference in shape + pad_size = np.subtract(self.template_px.shape, target_px.shape) assert pad_size[0] == pad_size[1], "Image matrices must be square." if pad_size[0] > 0: # pad target--UNTESTED + # add pixels to target if smaller than template target_px = np.pad(target_px, pad_width=(0, pad_size[0])) elif pad_size[0] < 0: # pad template - template_px = np.pad(template_px, pad_width=(0, -pad_size[0])) + # add pixels to template if smaller than target + self.template_px = np.pad(self.template_px, pad_width=(0, -pad_size[0])) # Always fit on magnitude images for simplicity. May be suboptimal - self.template8bit = \ - cv.normalize(abs(template_px), - None, 0, 255, norm_type=cv.NORM_MINMAX, - dtype=cv.CV_8U) + self.template8bit = cv.normalize(abs(self.template_px), None, 0, 255, + norm_type=cv.NORM_MINMAX, dtype=cv.CV_8U) - self.target8bit = cv.normalize(abs(target_px), - None, 0, 255, norm_type=cv.NORM_MINMAX, - dtype=cv.CV_8U) + self.target8bit = cv.normalize(abs(target_px), None, 0, 255, + norm_type=cv.NORM_MINMAX, dtype=cv.CV_8U) # initialise transformation fitting parameters. number_of_iterations = 500 termination_eps = 1e-10 criteria = (cv.TERM_CRITERIA_EPS | cv.TERM_CRITERIA_COUNT, number_of_iterations, termination_eps) + scale_factor = len(target_px) / len(self.template_px) self.warp_matrix = scale_factor * np.eye(2, 3, dtype=np.float32) self.scaled_template8bit = cv.warpAffine(self.template8bit, @@ -768,9 +767,10 @@ def plot_rois(self, new_fig=True): return fig - def order_by(self, att): + def order_by(self, images, att): """Order images by attribute (e.g. EchoTime, InversionTime).""" - self.images.sort(key=lambda x: x[att].value.real) + sorted_images = sorted(images, key=lambda x: x[att].value.real) + return sorted_images def generate_time_series(self, coords_row_col, fit_coords=True, kernel=None): @@ -801,7 +801,7 @@ def generate_time_series(self, coords_row_col, fit_coords=True, self.ROI_time_series = [] for i in range(num_coords): self.ROI_time_series.append(ROITimeSeries( - self.images, coords_row_col[i], time_attr=self.dicom_order_key, + self.images, coords_row_col[i], time_attr=self.time_attribute, kernel=kernel)) def generate_fit_function(self): @@ -821,7 +821,7 @@ class T1ImageStack(ImageStack): def __init__(self, image_slices, template_dcm=None): super().__init__(image_slices, template_dcm, - dicom_order_key='InversionTime') + time_attribute='InversionTime') def generate_fit_function(self): """"Create T1 fit function for magnitude/signed image and variable TI.""" @@ -916,7 +916,7 @@ class T2ImageStack(ImageStack): def __init__(self, image_slices, template_dcm=None): super().__init__(image_slices, template_dcm, - dicom_order_key='EchoTime') + time_attribute='EchoTime') self.fit_function = t2_function self.fit_jacobian = None @@ -1062,7 +1062,6 @@ def main(data, plate_number=None, calc: str = 'T1', verbose=False, plate_number = int(plate_number) # convert to int if required except (ValueError, TypeError): pass # will raise error at next statement - if plate_number not in [4, 5]: raise hazenlib.exceptions.ArgumentCombinationError( 'Must specify plate_number (4 or 5)') @@ -1078,7 +1077,6 @@ def main(data, plate_number=None, calc: str = 'T1', verbose=False, print(f'Could not find template with plate number: {plate_number}.' f' Please pass plate number as arg.') exit() - elif calc in ['T2', 't2']: relax_str = calc.lower() try: From f3ea1e2775145b85fd02a288abaecea1a9971b6e Mon Sep 17 00:00:00 2001 From: sophie22 Date: Tue, 25 Jul 2023 19:50:14 +0100 Subject: [PATCH 06/17] make functions more compact and explicit (works for T1) --- hazenlib/config.py | 11 +- hazenlib/relaxometry.py | 522 ++++++++++++++++++---------------------- 2 files changed, 244 insertions(+), 289 deletions(-) diff --git a/hazenlib/config.py b/hazenlib/config.py index f09afec2..014ebc62 100644 --- a/hazenlib/config.py +++ b/hazenlib/config.py @@ -1,10 +1,6 @@ import os.path import numpy as np -# Parameters for Rician noise model -MAX_RICIAN_NOISE = 20.0 -SEED_RICIAN_NOISE = 5.0 - TEMPLATE_DIR = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'data', 'relaxometry') TEMPLATE_VALUES = { @@ -76,3 +72,10 @@ "t1": range(0, 1000, 10), "t2": range(0, 500, 5) } + +TEMPLATE_FIT_ITERS = 500 +TERMINATION_EPS = 1e-10 + +# Parameters for Rician noise model - used in T2 calculation +MAX_RICIAN_NOISE = 20.0 +SEED_RICIAN_NOISE = 5.0 diff --git a/hazenlib/relaxometry.py b/hazenlib/relaxometry.py index 4579bdb5..702059d5 100644 --- a/hazenlib/relaxometry.py +++ b/hazenlib/relaxometry.py @@ -136,7 +136,11 @@ from scipy.special import i0e, ive import hazenlib.exceptions -from hazenlib.config import MAX_RICIAN_NOISE, SEED_RICIAN_NOISE, TEMPLATE_VALUES, SMOOTH_TIMES +from hazenlib.config import ( + MAX_RICIAN_NOISE, SEED_RICIAN_NOISE, TEMPLATE_VALUES, SMOOTH_TIMES, + TEMPLATE_FIT_ITERS, TERMINATION_EPS +) + # Use dict to store template and reference information # Coordinates are in array format (row,col), rather than plt.patches @@ -292,84 +296,6 @@ def pixel_rescale(dcm): dcm.pixel_array, dcm) -def generate_t1_function(ti_interp_vals, tr_interp_vals, mag_image=False): - """ - Generate T1 signal function and jacobian with interpolated TRs. - - Signal intensity on T1 decay is a function of both TI and TR. Ideally, TR - should be constant and at least 5*T1. However, scan time can be reduced by - allowing a shorter TR which increases at long TIs. For example:: - TI | 50 | 100 | 200 | 400 | 600 | 800 - ---+------+------+------+------+------+------ - TR | 1000 | 1000 | 1000 | 1260 | 1860 | 2460 - - This function factory returns a function which calculates the signal - magnitude using the expression:: - S = S0 * (1 - a1 * np.exp(-TI / t1) + np.exp(-TR / t1)) - where ``S0`` is the recovered intensity, ``a1`` is theoretically 2.0 but - varies due to inhomogeneous B0 field, ``t1`` is the longitudinal - relaxation time, and the repetition time, ``TR``, is calculated from - ``TI`` using piecewise linear interpolation. - - Parameters - ---------- - ti_interp_vals : array_like - Array of TI values used as a look-up table to calculate TR - tr_interp_vals : array_like - Array of TR values used as a lookup table to calculate TR from the TI - used in the sequence. - mag_image : bool, optional - If True, the generated function returns the magnitude of the signal - (i.e. negative outputs become positive). The default is False. - - Returns - ------- - t1_function : function - S = S0 * (1 - a1 * np.exp(-TI / t1) + np.exp(-TR / t1)) - - t1_jacobian : function - Tuple of partial derivatives for curve fitting. - - eqn_str : string - String representation of fit function. - """ - # Create piecewise liner fit function. k=1 gives linear, s=0 ensures all - # points are on line. Using UnivariateSpline (rather than numpy.interp()) - # enables derivative calculation if required. - tr = UnivariateSpline(ti_interp_vals, tr_interp_vals, k=1, s=0) - # tr_der = tr.derivative() - - eqn_str = 's0 * (1 - a1 * np.exp(-TI / t1) + np.exp(-TR / t1))' - if mag_image: - eqn_str = f'abs({eqn_str})' - - def _t1_function_signed(ti, t1, s0, a1): - pv = s0 * (1 - a1 * np.exp(-ti / t1) + np.exp(-tr(ti) / t1)) - return pv - - def t1_function(ti, t1, s0, a1): - pv = _t1_function_signed(ti, t1, s0, a1) - if mag_image: - return abs(pv) - else: - return pv - - def t1_jacobian(ti, t1, s0, a1): - t1_der = s0 / (t1 ** 2) * (-ti * a1 * np.exp(-ti / t1) + tr(ti) - * np.exp(-tr(ti) / t1)) - s0_der = 1 - a1 * np.exp(-ti / t1) + np.exp(-tr(ti) / t1) - a1_der = -s0 * np.exp(-ti / t1) - jacobian = np.array([t1_der, s0_der, a1_der]) - - if mag_image: - pv = _t1_function_signed(ti, t1, s0, a1) - jacobian = (jacobian * (pv >= 0)) - (jacobian * (pv < 0)) - - return jacobian.T - - return t1_function, t1_jacobian, eqn_str - - def est_t1_s0(ti, tr, t1, pv): """ Return initial guess of s0 to seed T1 curve fitting. @@ -394,53 +320,6 @@ def est_t1_s0(ti, tr, t1, pv): return -pv / (1 - 2 * np.exp(-ti / t1) + np.exp(-tr / t1)) -def t2_function(te, t2, s0, c): - r""" - Calculated pixel value with Rician noise model. - - Calculates pixel value from [1]_:: - .. math:: - S=\sqrt{\frac{\pi \alpha^2}{2}} \exp(- \alpha) \left( (1+ 2 \alpha) - \ \text{I_0}(\alpha) + 2 \alpha \ \text{I_1}(\alpha) \right) - - \alpha() = \left( \frac{S_0}{2 \sigma} \ \exp{\left(-\frac{\text{TE}}{\text{T}_2}\right)} \right)^2 - - \text{I}_n() = n^\text{th} \ \text{order modified Bessel function of the first kind} - - Parameters - ---------- - te : array_like - Echo times. - t2 : float - T2 decay constant. - S0 : float - Initial pixel magnitude. - C : float - Noise parameter for Rician model (equivalent to st dev). - - Returns - ------- - pv : array_like - Theoretical pixel values (signal) at each TE. - - References - ---------- - .. [1] Raya, J.G., Dietrich, O., Horng, A., Weber, J., Reiser, M.F. and - Glaser, C., 2010. T2 measurement in articular cartilage: impact of the - fitting method on accuracy and precision at low SNR. Magnetic Resonance in - Medicine: An Official Journal of the International Society for Magnetic - Resonance in Medicine, 63(1), pp.181-193. - """ - - alpha = (s0 / (2 * c) * np.exp(-te / t2)) **2 - # NB need to use `i0e` and `ive` below to avoid numeric inaccuracy from - # multiplying by huge exponentials then dividing by the same exponential - pv = np.sqrt(np.pi/2 * c ** 2) * \ - ((1 + 2 * alpha) * i0e(alpha) + 2 * alpha * ive(1, alpha)) - - return pv - - def est_t2_s0(te, t2, pv, c=0.0): """ Initial guess for s0 to seed curve fitting:: @@ -502,10 +381,7 @@ class ROITimeSeries: Mean pixel value of ROI for each image in series. """ - SAMPLE_ELEMENT = skimage.morphology.square(5) - - def __init__(self, dcm_images, poi_coords_row_col, time_attr=None, - kernel=None): + def __init__(self, dcm_images, poi_coords_row_col, time_attr): """ Create ROITimeSeries for ROI parameters at sequential scans. @@ -522,25 +398,26 @@ def __init__(self, dcm_images, poi_coords_row_col, time_attr=None, ``'InversionTime'`` or ``'EchoTime'``) and store in the list ``self.times``. The default is ``None``, which does not create ``self.times`` - kernel : array_like, optional - Structuring element which defines ROI size and shape, centred on - POI. Each element should be 1 or 0, otherwise calculation of mean - will be incorrect. If ``None``, use a 5x5 square. The default is - ``None``. """ - if kernel is None: - kernel = self.SAMPLE_ELEMENT self.POI_mask = np.zeros((dcm_images[0].pixel_array.shape[0], dcm_images[0].pixel_array.shape[1]), dtype=np.int8) self.POI_mask[poi_coords_row_col[0], poi_coords_row_col[1]] = 1 + + kernel = skimage.morphology.square(5) + """kernel + Structuring element which defines ROI size and shape, centred on + POI. Each element should be 1 or 0, otherwise calculation of mean + will be incorrect. If ``None``, use a 5x5 square. The default is + ``None``. + """ + self.ROI_mask = np.zeros_like(self.POI_mask) self.ROI_mask = scipy.ndimage.filters.convolve(self.POI_mask, kernel) - if time_attr is not None: - self.times = [x[time_attr].value.real for x in dcm_images] + self.times = [x[time_attr].value.real for x in dcm_images] self.pixel_values = [ pixel_rescale(img)[self.ROI_mask > 0] for img in dcm_images] @@ -567,8 +444,7 @@ class ImageStack(): Object to hold image_slices and methods for T1, T2 calculation. """ - def __init__(self, image_slices, template_dcm, - time_attribute=None): + def __init__(self, image_slices, time_attribute=None): """ Create ImageStack object. @@ -588,34 +464,28 @@ def __init__(self, image_slices, template_dcm, DICOM attribute to order images. Typically 'InversionTime' for T1 relaxometry or 'EchoTime' for T2. """ - # Store template pixel array, after scaling in 0028,1052 and - # 0028,1053 applied - self.template_dcm = template_dcm - if template_dcm is not None: - self.template_px = pixel_rescale(template_dcm) - # store sorted images - if time_attribute is not None: - self.images = self.order_by(image_slices, time_attribute) - else: - self.images = image_slices + self.images = self.order_by(image_slices, time_attribute) b0_val = self.images[0]['MagneticFieldStrength'].value - if b0_val == 1.5: - self.b0_str = '1.5T' - elif b0_val == 3.0: - self.b0_str = '3.0T' - else: + if b0_val not in [1.5, 3.0]: # TODO incorporate warning through e.g. logging module print('Unable to match B0 to default values. Using 1.5T.\n' f" {self.images[0]['MagneticFieldStrength']}") self.b0_str = '1.5T' + else: + self.b0_str = f"{b0_val}T" - def template_fit(self, image_index=0): + def order_by(self, images, att): + """Order images by attribute (e.g. EchoTime, InversionTime).""" + sorted_images = sorted(images, key=lambda x: x[att].value.real) + return sorted_images + + def template_fit(self, template_dcm, image_index=0): """ Calculate transformation matrix to fit template to image. - The template pixel array, self.template_px, is fitted to one of the + The template pixel array, template_px, is fitted to one of the images in self.images (default=0). The resultant RT matrix is stored as self.warp_matrix. @@ -655,50 +525,79 @@ def template_fit(self, image_index=0): Despite these limitations, this method works well in practice for small angle rotations. """ + # Store template pixel array, after scaling in 0028,1052 and + # 0028,1053 applied + template_px = pixel_rescale(template_dcm) target_px = pixel_rescale(self.images[0]) ## Pad template or target pixels if required # Determine difference in shape - pad_size = np.subtract(self.template_px.shape, target_px.shape) + pad_size = np.subtract(template_px.shape, target_px.shape) assert pad_size[0] == pad_size[1], "Image matrices must be square." if pad_size[0] > 0: # pad target--UNTESTED # add pixels to target if smaller than template target_px = np.pad(target_px, pad_width=(0, pad_size[0])) elif pad_size[0] < 0: # pad template # add pixels to template if smaller than target - self.template_px = np.pad(self.template_px, pad_width=(0, -pad_size[0])) + template_px = np.pad(template_px, pad_width=(0, -pad_size[0])) # Always fit on magnitude images for simplicity. May be suboptimal - self.template8bit = cv.normalize(abs(self.template_px), None, 0, 255, + self.template8bit = cv.normalize(abs(template_px), None, 0, 255, norm_type=cv.NORM_MINMAX, dtype=cv.CV_8U) self.target8bit = cv.normalize(abs(target_px), None, 0, 255, norm_type=cv.NORM_MINMAX, dtype=cv.CV_8U) # initialise transformation fitting parameters. - number_of_iterations = 500 - termination_eps = 1e-10 - criteria = (cv.TERM_CRITERIA_EPS | cv.TERM_CRITERIA_COUNT, - number_of_iterations, termination_eps) - scale_factor = len(target_px) / len(self.template_px) - self.warp_matrix = scale_factor * np.eye(2, 3, dtype=np.float32) + scale_factor = len(target_px) / len(template_px) + scale_matrix = scale_factor * np.eye(2, 3, dtype=np.float32) - self.scaled_template8bit = cv.warpAffine(self.template8bit, - self.warp_matrix, - (self.template8bit.shape[1], - self.template8bit.shape[0])) + self.scaled_template8bit = cv.warpAffine( + self.template8bit, scale_matrix, + (self.template8bit.shape[1],self.template8bit.shape[0])) # Apply transformation - self.template_cc, self.warp_matrix = \ - cv.findTransformECC(self.template8bit, self.target8bit, - self.warp_matrix, criteria=criteria) + criteria = (cv.TERM_CRITERIA_EPS | cv.TERM_CRITERIA_COUNT, + TEMPLATE_FIT_ITERS, TERMINATION_EPS) + # Find the geometric transform (warp) between two images in terms of the ECC criterion + self.template_cc, warp_matrix = cv.findTransformECC( + self.template8bit, self.target8bit, + scale_matrix, criteria=criteria) + + self.warped_template8bit = cv.warpAffine( + self.template8bit, warp_matrix, + (self.template8bit.shape[1], self.template8bit.shape[0])) + + return warp_matrix + + def generate_time_series(self, coords_row_col, time_attribute, + warp_matrix): + """ + Create list of ROITimeSeries objects. + + Parameters + ---------- + coords_row_col : array_like + Array of coordinates points of interest (POIs) for each centre of + each ROI. They should be in [[col0, row0], [col1, row1], ...] + format. + time_attribute : string, optional + DICOM attribute to order images. Typically 'InversionTime' for T1 + relaxometry or 'EchoTime' for T2. + warp_matrix : np.array + RT transform matrix (2,3). + """ + num_coords = np.size(coords_row_col, axis=0) + flipped_coords_row_col = transform_coords(coords_row_col, warp_matrix, + input_row_col=True, output_row_col=True) - self.warped_template8bit = cv.warpAffine(self.template8bit, - self.warp_matrix, - (self.template8bit.shape[1], - self.template8bit.shape[0])) + self.ROI_time_series = [] + for i in range(num_coords): + self.ROI_time_series.append(ROITimeSeries( + self.images, flipped_coords_row_col[i], time_attribute)) - return self.warp_matrix + def generate_fit_function(self): + """Null method in base class, may be overwritten in subclass.""" def plot_fit(self): """ @@ -767,61 +666,101 @@ def plot_rois(self, new_fig=True): return fig - def order_by(self, images, att): - """Order images by attribute (e.g. EchoTime, InversionTime).""" - sorted_images = sorted(images, key=lambda x: x[att].value.real) - return sorted_images + @property + def relax_times(self): + """List of T1 for each ROI.""" + return [fit[0][0] for fit in self.relax_fit] + +class T1ImageStack(ImageStack): + """ + Calculate T1 relaxometry. + + Overloads the following methods from ``ImageStack``: + ``generate_fit_function`` + ``initialise_fit_parameters`` + ``find_relax_times`` + + """ + + def __init__(self, image_slices, time_attribute): + super().__init__(image_slices, time_attribute) - def generate_time_series(self, coords_row_col, fit_coords=True, - kernel=None): + def generate_t1_function(self, ti_interp_vals, tr_interp_vals, mag_image=False): """ - Create list of ROITimeSeries objects. + Generate T1 signal function and jacobian with interpolated TRs. + + Signal intensity on T1 decay is a function of both TI and TR. Ideally, TR + should be constant and at least 5*T1. However, scan time can be reduced by + allowing a shorter TR which increases at long TIs. For example:: + TI | 50 | 100 | 200 | 400 | 600 | 800 + ---+------+------+------+------+------+------ + TR | 1000 | 1000 | 1000 | 1260 | 1860 | 2460 + + This function factory returns a function which calculates the signal + magnitude using the expression:: + S = S0 * (1 - a1 * np.exp(-TI / t1) + np.exp(-TR / t1)) + where ``S0`` is the recovered intensity, ``a1`` is theoretically 2.0 but + varies due to inhomogeneous B0 field, ``t1`` is the longitudinal + relaxation time, and the repetition time, ``TR``, is calculated from + ``TI`` using piecewise linear interpolation. Parameters ---------- - coords_row_col : array_like - Array of coordinates points of interest (POIs) for each centre of - each ROI. They should be in [[col0, row0], [col1, row1], ...] - format. - fit_coords : bool, optional - If ``True``, the coordinates provided are for the template ROIs and - will be transformed to the image space using ``transfor_coords()``. - The default is True. - kernel : array, optional - Structuring element which should be an array of 1s and possibly 0s. - If ``None``, use the default from ``ROItimeSeries`` constructor. - The default is None. + ti_interp_vals : array_like + Array of TI values used as a look-up table to calculate TR + tr_interp_vals : array_like + Array of TR values used as a lookup table to calculate TR from the TI + used in the sequence. + mag_image : bool, optional + If True, the generated function returns the magnitude of the signal + (i.e. negative outputs become positive). The default is False. + + Returns + ------- + t1_function : function + S = S0 * (1 - a1 * np.exp(-TI / t1) + np.exp(-TR / t1)) + + t1_jacobian : function + Tuple of partial derivatives for curve fitting. + + eqn_str : string + String representation of fit function. """ - num_coords = np.size(coords_row_col, axis=0) - if fit_coords: - coords_row_col = transform_coords(coords_row_col, self.warp_matrix, - input_row_col=True, - output_row_col=True) + # Create piecewise linear fit function. k=1 gives linear, s=0 ensures all + # points are on line. Using UnivariateSpline (rather than numpy.interp()) + # enables derivative calculation if required. + tr = UnivariateSpline(ti_interp_vals, tr_interp_vals, k=1, s=0) + # tr_der = tr.derivative() - self.ROI_time_series = [] - for i in range(num_coords): - self.ROI_time_series.append(ROITimeSeries( - self.images, coords_row_col[i], time_attr=self.time_attribute, - kernel=kernel)) + eqn_str = 's0 * (1 - a1 * np.exp(-TI / t1) + np.exp(-TR / t1))' + if mag_image: + eqn_str = f'abs({eqn_str})' - def generate_fit_function(self): - """Null method in base class, may be overwritten in subclass.""" + def _t1_function_signed(ti, t1, s0, a1): + pv = s0 * (1 - a1 * np.exp(-ti / t1) + np.exp(-tr(ti) / t1)) + return pv + def t1_function(ti, t1, s0, a1): + pv = _t1_function_signed(ti, t1, s0, a1) + if mag_image: + return abs(pv) + else: + return pv -class T1ImageStack(ImageStack): - """ - Calculate T1 relaxometry. + def t1_jacobian(ti, t1, s0, a1): + t1_der = s0 / (t1 ** 2) * (-ti * a1 * np.exp(-ti / t1) + tr(ti) + * np.exp(-tr(ti) / t1)) + s0_der = 1 - a1 * np.exp(-ti / t1) + np.exp(-tr(ti) / t1) + a1_der = -s0 * np.exp(-ti / t1) + jacobian = np.array([t1_der, s0_der, a1_der]) - Overloads the following methods from ``ImageStack``: - ``generate_fit_function`` - ``initialise_fit_parameters`` - ``find_relax_times`` + if mag_image: + pv = _t1_function_signed(ti, t1, s0, a1) + jacobian = (jacobian * (pv >= 0)) - (jacobian * (pv < 0)) - """ + return jacobian.T - def __init__(self, image_slices, template_dcm=None): - super().__init__(image_slices, template_dcm, - time_attribute='InversionTime') + return t1_function, t1_jacobian, eqn_str def generate_fit_function(self): """"Create T1 fit function for magnitude/signed image and variable TI.""" @@ -830,10 +769,10 @@ def generate_fit_function(self): mag_image = True else: mag_image = False - self.fit_function, self.fit_jacobian, self.fit_eqn_str = \ - generate_t1_function(self.ROI_time_series[0].times, - self.ROI_time_series[0].trs, - mag_image=mag_image) + self.t1_fit_function, self.fit_jacobian, self.fit_eqn_str = \ + self.generate_t1_function( + self.ROI_time_series[0].times, self.ROI_time_series[0].trs, + mag_image=mag_image) def initialise_fit_parameters(self, t1_estimates): """ @@ -865,9 +804,9 @@ def initialise_fit_parameters(self, t1_estimates): rois_first_mean = np.array([roi.means[0] for roi in rois]) rois_last_mean = np.array([roi.means[-1] for roi in rois]) s0_est_last = abs(est_t1_s0(rois[0].times[-1], rois[0].trs[-1], - t1_estimates, rois_last_mean)) + self.t1_est, rois_last_mean)) s0_est_first = abs(est_t1_s0(rois[0].times[0], rois[0].trs[0], - t1_estimates, rois_first_mean)) + self.t1_est, rois_first_mean)) self.s0_est = np.where(rois_first_mean > rois_last_mean, s0_est_first, s0_est_last) self.a1_est = np.full_like(self.s0_est, 2.0) @@ -882,25 +821,10 @@ def find_relax_times(self): """ rois = self.ROI_time_series - self.relax_fit = [scipy.optimize.curve_fit(self.fit_function, - rois[i].times, - rois[i].means, - p0=[self.t1_est[i], - self.s0_est[i], - self.a1_est[i]], - jac=self.fit_jacobian, - method='lm') - for i in range(len(rois))] - - @property - def t1s(self): - """List T1 values for each ROI.""" - return [fit[0][0] for fit in self.relax_fit] - - @property - def relax_times(self): - """List of T1 for each ROI.""" - return self.t1s + self.relax_fit = [scipy.optimize.curve_fit( + self.t1_fit_function, rois[i].times, rois[i].means, + p0=[self.t1_est[i], self.s0_est[i], self.a1_est[i]], + jac=self.fit_jacobian, method='lm') for i in range(len(rois))] class T2ImageStack(ImageStack): @@ -914,12 +838,9 @@ class T2ImageStack(ImageStack): """ - def __init__(self, image_slices, template_dcm=None): - super().__init__(image_slices, template_dcm, - time_attribute='EchoTime') + def __init__(self, image_slices, time_attribute): + super().__init__(image_slices, time_attribute) - self.fit_function = t2_function - self.fit_jacobian = None self.fit_eqn_str = 'T2 with Rician noise (Raya et al 2010)' def initialise_fit_parameters(self, t2_estimates): @@ -948,9 +869,55 @@ def initialise_fit_parameters(self, t2_estimates): rois_second_mean = np.array([roi.means[1] for roi in rois]) self.c_est = np.full_like(self.t2_est, SEED_RICIAN_NOISE) # estimate s0 from second image--first image is too low. - self.s0_est = est_t2_s0(rois[0].times[1], t2_estimates, + self.s0_est = est_t2_s0(rois[0].times[1], self.t2_est, rois_second_mean, self.c_est) + def t2_fit_function(self, te, t2, s0, c): + r""" + Calculated pixel value with Rician noise model. + + Calculates pixel value from [1]_:: + .. math:: + S=\sqrt{\frac{\pi \alpha^2}{2}} \exp(- \alpha) \left( (1+ 2 \alpha) + \ \text{I_0}(\alpha) + 2 \alpha \ \text{I_1}(\alpha) \right) + + \alpha() = \left( \frac{S_0}{2 \sigma} \ \exp{\left(-\frac{\text{TE}}{\text{T}_2}\right)} \right)^2 + + \text{I}_n() = n^\text{th} \ \text{order modified Bessel function of the first kind} + + Parameters + ---------- + te : array_like + Echo times. + t2 : float + T2 decay constant. + S0 : float + Initial pixel magnitude. + C : float + Noise parameter for Rician model (equivalent to st dev). + + Returns + ------- + pv : array_like + Theoretical pixel values (signal) at each TE. + + References + ---------- + .. [1] Raya, J.G., Dietrich, O., Horng, A., Weber, J., Reiser, M.F. and + Glaser, C., 2010. T2 measurement in articular cartilage: impact of the + fitting method on accuracy and precision at low SNR. Magnetic Resonance in + Medicine: An Official Journal of the International Society for Magnetic + Resonance in Medicine, 63(1), pp.181-193. + """ + + alpha = (s0 / (2 * c) * np.exp(-te / t2)) **2 + # NB need to use `i0e` and `ive` below to avoid numeric inaccuracy from + # multiplying by huge exponentials then dividing by the same exponential + pv = np.sqrt(np.pi/2 * c ** 2) * \ + ((1 + 2 * alpha) * i0e(alpha) + 2 * alpha * ive(1, alpha)) + + return pv + def find_relax_times(self): """ Calculate T2 values. Access as ``image_stack.t2s``. @@ -982,35 +949,16 @@ def find_relax_times(self): rois = self.ROI_time_series # Omit the first image data from the curve fit. This is achieved by - # slicing rois[i].times[1:] and rois[i].means[1:]. Skipping odd echoes - # can be implemented with rois[i].times[1::2] and .means[1::2] + # slicing rois[i].times[1:] and rois[i].means[1:]. + # Skipping odd echoes can be implemented with + # rois[i].times[1::2] and .means[1::2] bounds = ([0, 0, 1], [np.inf, np.inf, MAX_RICIAN_NOISE]) - self.relax_fit = [] - for i in range(len(rois)): - self.relax_fit.append( - scipy.optimize.curve_fit(self.fit_function, - rois[i].times[1:], - rois[i].means[1:], - p0=[self.t2_est[i], - self.s0_est[i], - self.c_est[i]], - jac=self.fit_jacobian, - bounds=bounds, - method='trf' - ) - ) - - @property - def t2s(self): - """List of T2 values for each ROI.""" - return [fit[0][0] for fit in self.relax_fit] - - @property - def relax_times(self): - """List of T2 values for each ROI.""" - return self.t2s + self.relax_fit = [scipy.optimize.curve_fit( + self.t2_fit_function, rois[i].times[1:], rois[i].means[1:], + p0=[self.t2_est[i], self.s0_est[i], self.c_est[i]], + jac=None, bounds=bounds, method='trf') for i in range(len(rois))] def main(data, plate_number=None, calc: str = 'T1', verbose=False, @@ -1067,22 +1015,26 @@ def main(data, plate_number=None, calc: str = 'T1', verbose=False, 'Must specify plate_number (4 or 5)') # Set up parameters specific to T1 or T2 + relax_str = calc.lower() + time_attributes = { + "t1": 'InversionTime', + "t2": 'EchoTime' + } + if calc in ['T1', 't1']: - relax_str = calc.lower() + image_stack = T1ImageStack(data, time_attribute=time_attributes[relax_str]) try: template_dcm = pydicom.read_file( TEMPLATE_VALUES[f'plate{plate_number}'][relax_str]['filename']) - image_stack = T1ImageStack(data, template_dcm) except KeyError: print(f'Could not find template with plate number: {plate_number}.' f' Please pass plate number as arg.') exit() elif calc in ['T2', 't2']: - relax_str = calc.lower() + image_stack = T2ImageStack(data, time_attribute=time_attributes[relax_str]) try: template_dcm = pydicom.read_file( TEMPLATE_VALUES[f'plate{plate_number}'][relax_str]['filename']) - image_stack = T2ImageStack(data, template_dcm) except KeyError: print(f'Could not find template with plate number: {plate_number}.' f' Please pass plate number as arg.') @@ -1091,10 +1043,10 @@ def main(data, plate_number=None, calc: str = 'T1', verbose=False, print("Please provide 'T1' or 'T2' for the --calc argument.") exit() - image_stack.template_fit() + warp_matrix = image_stack.template_fit(template_dcm) image_stack.generate_time_series( - TEMPLATE_VALUES[f'plate{plate_number}'] - ['sphere_centres_row_col']) + TEMPLATE_VALUES[f'plate{plate_number}']['sphere_centres_row_col'], + time_attribute=time_attributes[relax_str], warp_matrix=warp_matrix) image_stack.generate_fit_function() # Published relaxation time for matching plate and T1/T2 From b991f238d72721ad7748f72dc0af876288dd4232 Mon Sep 17 00:00:00 2001 From: sophie22 Date: Tue, 25 Jul 2023 20:21:36 +0100 Subject: [PATCH 07/17] move T1/T2-specific functions into their class --- hazenlib/relaxometry.py | 220 ++++++++++++++++++++-------------------- 1 file changed, 112 insertions(+), 108 deletions(-) diff --git a/hazenlib/relaxometry.py b/hazenlib/relaxometry.py index 702059d5..d4d61ed7 100644 --- a/hazenlib/relaxometry.py +++ b/hazenlib/relaxometry.py @@ -296,57 +296,6 @@ def pixel_rescale(dcm): dcm.pixel_array, dcm) -def est_t1_s0(ti, tr, t1, pv): - """ - Return initial guess of s0 to seed T1 curve fitting. - - Parameters - ---------- - ti : array_like - TI values. - tr : array_like - TR values. - t1 : array_like - Estimated T1 (typically from manufacturer's documentation). - pv : array_like - Mean pixel value (signal) in ROI. - - Returns - ------- - array_like - Initial s0 guess for calculating T1 relaxation time. - - """ - return -pv / (1 - 2 * np.exp(-ti / t1) + np.exp(-tr / t1)) - - -def est_t2_s0(te, t2, pv, c=0.0): - """ - Initial guess for s0 to seed curve fitting:: - .. math:: - S_0=\\frac{pv-c}{exp(-TE/T_2)} - - - Parameters - ---------- - te : array_like - Echo time(s). - t2 : array_like - T2 decay constant. - pv : array_like - Mean pixel value (signal) in ROI with ``te`` echo time. - c : array_like - Constant offset, theoretically ``full_like(te, 0.0)``. - - Returns - ------- - array_like - Initial s0 estimate. - - """ - return (pv - c) / np.exp(-te / t2) - - class ROITimeSeries: """ Samples at one image location (ROI) at numerous sample times. @@ -596,9 +545,6 @@ def generate_time_series(self, coords_row_col, time_attribute, self.ROI_time_series.append(ROITimeSeries( self.images, flipped_coords_row_col[i], time_attribute)) - def generate_fit_function(self): - """Null method in base class, may be overwritten in subclass.""" - def plot_fit(self): """ Visual representation of target fitting. @@ -674,12 +620,6 @@ def relax_times(self): class T1ImageStack(ImageStack): """ Calculate T1 relaxometry. - - Overloads the following methods from ``ImageStack``: - ``generate_fit_function`` - ``initialise_fit_parameters`` - ``find_relax_times`` - """ def __init__(self, image_slices, time_attribute): @@ -774,6 +714,29 @@ def generate_fit_function(self): self.ROI_time_series[0].times, self.ROI_time_series[0].trs, mag_image=mag_image) + def est_t1_s0(self, ti, tr, t1, pv): + """ + Return initial guess of s0 to seed T1 curve fitting. + + Parameters + ---------- + ti : array_like + TI values. + tr : array_like + TR values. + t1 : array_like + Estimated T1 (typically from manufacturer's documentation). + pv : array_like + Mean pixel value (signal) in ROI. + + Returns + ------- + array_like + Initial s0 guess for calculating T1 relaxation time. + + """ + return -pv / (1 - 2 * np.exp(-ti / t1) + np.exp(-tr / t1)) + def initialise_fit_parameters(self, t1_estimates): """ Estimate fit parameters (t1, s0, a1) for T1 curve fitting. @@ -787,7 +750,7 @@ def initialise_fit_parameters(self, t1_estimates): causes a large rounding error. A1 is estimated as 2.0, the theoretical value assuming homogeneous B0 - + Parameters ---------- t1_estimates : array_like @@ -796,24 +759,32 @@ def initialise_fit_parameters(self, t1_estimates): Returns ------- - None. + s0_est """ self.t1_est = t1_estimates rois = self.ROI_time_series rois_first_mean = np.array([roi.means[0] for roi in rois]) rois_last_mean = np.array([roi.means[-1] for roi in rois]) - s0_est_last = abs(est_t1_s0(rois[0].times[-1], rois[0].trs[-1], + s0_est_last = abs(self.est_t1_s0(rois[0].times[-1], rois[0].trs[-1], self.t1_est, rois_last_mean)) - s0_est_first = abs(est_t1_s0(rois[0].times[0], rois[0].trs[0], + s0_est_first = abs(self.est_t1_s0(rois[0].times[0], rois[0].trs[0], self.t1_est, rois_first_mean)) - self.s0_est = np.where(rois_first_mean > rois_last_mean, + s0_est = np.where(rois_first_mean > rois_last_mean, s0_est_first, s0_est_last) - self.a1_est = np.full_like(self.s0_est, 2.0) + self.a1_est = np.full_like(s0_est, 2.0) - def find_relax_times(self): + return s0_est + + def find_relax_times(self, t1_estimates, s0_est): """ - Calculate T1 values. Access as ``image_stack.t1s`` + Calculate T1 values. Access as ``image_stack.relax_fits`` + + Parameters + ---------- + t1_estimates : array_like + T1 values to seed estimation. These should be the manufacturer + provided T1 values where known. Returns ------- @@ -823,19 +794,13 @@ def find_relax_times(self): rois = self.ROI_time_series self.relax_fit = [scipy.optimize.curve_fit( self.t1_fit_function, rois[i].times, rois[i].means, - p0=[self.t1_est[i], self.s0_est[i], self.a1_est[i]], + p0=[t1_estimates[i], s0_est[i], self.a1_est[i]], jac=self.fit_jacobian, method='lm') for i in range(len(rois))] class T2ImageStack(ImageStack): """ Calculate T2 relaxometry. - - Overloads the following methods from ``ImageStack``: - ``generate_fit_function`` - ``initialise_fit_parameters`` - ``find_relax_times`` - """ def __init__(self, image_slices, time_attribute): @@ -843,34 +808,8 @@ def __init__(self, image_slices, time_attribute): self.fit_eqn_str = 'T2 with Rician noise (Raya et al 2010)' - def initialise_fit_parameters(self, t2_estimates): - """ - Estimate fit parameters (t2, s0, c) for T2 curve fitting. - - T2 estimates are provided. - - s0 is estimated using est_t2_s0(te, t2_est, mean_pv, c). - - C is estimated as 5.0. - - Parameters - ---------- - t2_estimates : array_like - T2 values to seed estimation. These should be the manufacturer - provided T2 values where known. - - Returns - ------- - None. - - """ - self.t2_est = t2_estimates - rois = self.ROI_time_series - rois_second_mean = np.array([roi.means[1] for roi in rois]) - self.c_est = np.full_like(self.t2_est, SEED_RICIAN_NOISE) - # estimate s0 from second image--first image is too low. - self.s0_est = est_t2_s0(rois[0].times[1], self.t2_est, - rois_second_mean, self.c_est) + def generate_fit_function(self): + """Null method in base class, may be overwritten in subclass.""" def t2_fit_function(self, te, t2, s0, c): r""" @@ -918,9 +857,66 @@ def t2_fit_function(self, te, t2, s0, c): return pv - def find_relax_times(self): + def est_t2_s0(self, te, t2, pv, c=0.0): + """ + Initial guess for s0 to seed curve fitting:: + .. math:: + S_0=\\frac{pv-c}{exp(-TE/T_2)} + + + Parameters + ---------- + te : array_like + Echo time(s). + t2 : array_like + T2 decay constant. + pv : array_like + Mean pixel value (signal) in ROI with ``te`` echo time. + c : array_like + Constant offset, theoretically ``full_like(te, 0.0)``. + + Returns + ------- + array_like + Initial s0 estimate. + + """ + return (pv - c) / np.exp(-te / t2) + + def initialise_fit_parameters(self, t2_estimates): + """ + Estimate fit parameters (t2, s0, c) for T2 curve fitting. + + T2 estimates are provided. + + s0 is estimated using est_t2_s0(te, t2_est, mean_pv, c). + + C is estimated as 5.0. + + Parameters + ---------- + t2_estimates : array_like + T2 values to seed estimation. These should be the manufacturer + provided T2 values where known. + + Returns + ------- + None. + """ - Calculate T2 values. Access as ``image_stack.t2s``. + self.t2_est = t2_estimates + rois = self.ROI_time_series + rois_second_mean = np.array([roi.means[1] for roi in rois]) + self.c_est = np.full_like(self.t2_est, SEED_RICIAN_NOISE) + # estimate s0 from second image--first image is too low. + s0_est = self.est_t2_s0(rois[0].times[1], self.t2_est, + rois_second_mean, self.c_est) + + return s0_est + + def find_relax_times(self, t2_estimates, s0_est): + """ + Calculate T2 values. Access as ``image_stack.relax_times`` Uses the 'skip first echo' fit method [1]_ with a Rician noise model [2]_. Ideally the Rician noise parameter should be determined from the @@ -930,6 +926,12 @@ def find_relax_times(self): parameter makes this easier. It has an upper limit of MAX_RICIAN_NOISE, currently set to 20.0. + Parameters + ---------- + t2_estimates : array_like + T2 values to seed estimation. These should be the manufacturer + provided T2 values where known. + Returns ------- None. @@ -957,7 +959,7 @@ def find_relax_times(self): self.relax_fit = [scipy.optimize.curve_fit( self.t2_fit_function, rois[i].times[1:], rois[i].means[1:], - p0=[self.t2_est[i], self.s0_est[i], self.c_est[i]], + p0=[t2_estimates[i], s0_est[i], self.c_est[i]], jac=None, bounds=bounds, method='trf') for i in range(len(rois))] @@ -1047,13 +1049,15 @@ def main(data, plate_number=None, calc: str = 'T1', verbose=False, image_stack.generate_time_series( TEMPLATE_VALUES[f'plate{plate_number}']['sphere_centres_row_col'], time_attribute=time_attributes[relax_str], warp_matrix=warp_matrix) + # only applies to T1 image_stack.generate_fit_function() # Published relaxation time for matching plate and T1/T2 relax_published = TEMPLATE_VALUES [f'plate{plate_number}'][relax_str] \ ['relax_times'][image_stack.b0_str] - image_stack.initialise_fit_parameters(relax_published) - image_stack.find_relax_times() + s0_est = image_stack.initialise_fit_parameters(relax_published) + + image_stack.find_relax_times(relax_published, s0_est) frac_time_diff = (image_stack.relax_times - relax_published) \ / relax_published # last value is for background water. Strip before calculating RMS frac error From 4c79352cf5e459b8eacc130da0f6e3cde5b9754d Mon Sep 17 00:00:00 2001 From: sophie22 Date: Tue, 25 Jul 2023 20:22:12 +0100 Subject: [PATCH 08/17] add test for T2 relaxometry --- .github/workflows/test_cli.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/test_cli.yml b/.github/workflows/test_cli.yml index 4cee1ffc..250fbdc3 100644 --- a/.github/workflows/test_cli.yml +++ b/.github/workflows/test_cli.yml @@ -111,5 +111,6 @@ jobs: - name: test relaxometry if: always() run: | - hazen relaxometry tests/data/relaxometry/T1/site1_20200218/plate5 --calc=t1 --plate_number=5 --report + hazen relaxometry tests/data/relaxometry/T1/site1_20200218/plate5 --calc T1 --plate_number=5 --report + hazen relaxometry tests/data/relaxometry/T2/site3_ge/plate4/ --calc T2 --plate_number=4 --report From bb11e33d64ffe0142d0af472ea1a697192fcf51b Mon Sep 17 00:00:00 2001 From: sophie22 Date: Tue, 25 Jul 2023 20:38:31 +0100 Subject: [PATCH 09/17] add back the fit_function --- hazenlib/relaxometry.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/hazenlib/relaxometry.py b/hazenlib/relaxometry.py index d4d61ed7..2d62799f 100644 --- a/hazenlib/relaxometry.py +++ b/hazenlib/relaxometry.py @@ -709,7 +709,7 @@ def generate_fit_function(self): mag_image = True else: mag_image = False - self.t1_fit_function, self.fit_jacobian, self.fit_eqn_str = \ + self.fit_function, self.fit_jacobian, self.fit_eqn_str = \ self.generate_t1_function( self.ROI_time_series[0].times, self.ROI_time_series[0].trs, mag_image=mag_image) @@ -793,7 +793,7 @@ def find_relax_times(self, t1_estimates, s0_est): """ rois = self.ROI_time_series self.relax_fit = [scipy.optimize.curve_fit( - self.t1_fit_function, rois[i].times, rois[i].means, + self.fit_function, rois[i].times, rois[i].means, p0=[t1_estimates[i], s0_est[i], self.a1_est[i]], jac=self.fit_jacobian, method='lm') for i in range(len(rois))] @@ -811,8 +811,8 @@ def __init__(self, image_slices, time_attribute): def generate_fit_function(self): """Null method in base class, may be overwritten in subclass.""" - def t2_fit_function(self, te, t2, s0, c): - r""" + def fit_function(self, te, t2, s0, c): + """ Calculated pixel value with Rician noise model. Calculates pixel value from [1]_:: @@ -958,7 +958,7 @@ def find_relax_times(self, t2_estimates, s0_est): bounds = ([0, 0, 1], [np.inf, np.inf, MAX_RICIAN_NOISE]) self.relax_fit = [scipy.optimize.curve_fit( - self.t2_fit_function, rois[i].times[1:], rois[i].means[1:], + self.fit_function, rois[i].times[1:], rois[i].means[1:], p0=[t2_estimates[i], s0_est[i], self.c_est[i]], jac=None, bounds=bounds, method='trf') for i in range(len(rois))] From dca71b00d90d161c5f1acfbd78e2072f5890fa85 Mon Sep 17 00:00:00 2001 From: sophie22 Date: Wed, 26 Jul 2023 00:51:12 +0100 Subject: [PATCH 10/17] b0_val could be "3" --- hazenlib/relaxometry.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/hazenlib/relaxometry.py b/hazenlib/relaxometry.py index 2d62799f..02a9aa5c 100644 --- a/hazenlib/relaxometry.py +++ b/hazenlib/relaxometry.py @@ -423,7 +423,10 @@ def __init__(self, image_slices, time_attribute=None): f" {self.images[0]['MagneticFieldStrength']}") self.b0_str = '1.5T' else: - self.b0_str = f"{b0_val}T" + if b0_val == 3: + self.b0_str = "3.0T" + else: + self.b0_str = f"{b0_val}T" def order_by(self, images, att): """Order images by attribute (e.g. EchoTime, InversionTime).""" From 86c504238cef76cf22761873d16017710048bfaa Mon Sep 17 00:00:00 2001 From: sophie22 Date: Wed, 26 Jul 2023 00:55:41 +0100 Subject: [PATCH 11/17] fix most tests with updated args --- tests/test_relaxometry.py | 91 ++++++++++++++++++++++----------------- 1 file changed, 52 insertions(+), 39 deletions(-) diff --git a/tests/test_relaxometry.py b/tests/test_relaxometry.py index 49dd4fa3..e9164c66 100644 --- a/tests/test_relaxometry.py +++ b/tests/test_relaxometry.py @@ -14,6 +14,7 @@ import hazenlib.relaxometry as hazen_relaxometry from hazenlib.exceptions import ArgumentCombinationError from tests import TEST_DATA_DIR, TEST_REPORT_DIR +from hazenlib.config import TEMPLATE_VALUES class TestRelaxometry(unittest.TestCase): @@ -199,7 +200,7 @@ class TestRelaxometry(unittest.TestCase): SITE3_T2_P5_FILES = ['Z812', 'Z813', 'Z814', 'Z819', 'Z823', 'Z825', 'Z830', 'Z834'] - # Site 5 tests - Philisp 3T + # Site 5 tests - Philips 3T SITE5_T1_P4_DIR = os.path.join(TEST_DATA_DIR, 'relaxometry', 'T1', 'site5_philips_3T', 'plate4') @@ -317,11 +318,11 @@ def test_template_fit(self): target_dcm = pydicom.dcmread(self.TEMPLATE_TARGET_PATH_T1_P5) t1_image_stack = hazen_relaxometry.T1ImageStack([target_dcm], - template_dcm) - t1_image_stack.template_fit() + time_attribute="InversionTime") + warp_matrix = t1_image_stack.template_fit(template_dcm) transformed_coordinates_xy = hazen_relaxometry.transform_coords( - self.TEMPLATE_TEST_COORDS_ROW_COL, t1_image_stack.warp_matrix, + self.TEMPLATE_TEST_COORDS_ROW_COL, warp_matrix, input_row_col=True, output_row_col=False) # test to within +/- 1 pixel (also checks YX-XY change) @@ -333,7 +334,8 @@ def test_image_stack_T1_sort(self): # read list of un-ordered T1 files, sort by TI, test sorted t1_dcms = [pydicom.dcmread(os.path.join(self.T1_DIR, fname)) for fname in self.T1_FILES] - t1_image_stack = hazen_relaxometry.T1ImageStack(t1_dcms) + t1_image_stack = hazen_relaxometry.T1ImageStack(t1_dcms, + time_attribute="InversionTime") sorted_output = [image.InversionTime.real for image in t1_image_stack.images] assert sorted_output == self.T1_TI_SORTED @@ -342,7 +344,8 @@ def test_image_stack_T2_sort(self): # read list of un-ordered T2 files, sort by TE, test sorted t2_dcms = [pydicom.dcmread(os.path.join(self.T2_DIR, fname)) for fname in self.T2_FILES] - t2_image_stack = hazen_relaxometry.T2ImageStack(t2_dcms) + t2_image_stack = hazen_relaxometry.T2ImageStack(t2_dcms, + time_attribute="EchoTime") sorted_output = [image.EchoTime.real for image in t2_image_stack.images] @@ -353,10 +356,12 @@ def test_generate_time_series_template_POIs(self): # Test on template first, no image fitting needed # Need image to get correct size template_dcm = pydicom.dcmread(self.TEMPLATE_PATH_T1_P5) - template_image_stack = hazen_relaxometry.T1ImageStack([template_dcm]) + template_image_stack = hazen_relaxometry.T1ImageStack([template_dcm], + time_attribute="InversionTime") + warp_matrix = template_image_stack.template_fit(template_dcm) template_image_stack.generate_time_series( - self.TEMPLATE_TEST_COORDS_ROW_COL, - fit_coords=False) + self.TEMPLATE_TEST_COORDS_ROW_COL, time_attribute="InversionTime", + warp_matrix=warp_matrix, fit_coords=False) for i in range(np.size(self.MASK_POI_TEMPLATE, 0)): np.testing.assert_equal( @@ -369,13 +374,14 @@ def test_generate_time_series_target_POIs(self): target_dcm = pydicom.dcmread(self.TEMPLATE_TARGET_PATH_T1_P5) target_image_stack = hazen_relaxometry.T1ImageStack([target_dcm], - template_dcm) - target_image_stack.template_fit() + time_attribute="InversionTime") + warp_matrix = target_image_stack.template_fit(template_dcm) # transformed_coordinates_yx = hazen_relaxometry.transform_coords( # self.TEMPLATE_TEST_COORDS_YX, target_image_stack.warp_matrix, # input_yx=True, output_yx=True) target_image_stack.generate_time_series( - self.TEMPLATE_TEST_COORDS_ROW_COL, fit_coords=True) + self.TEMPLATE_TEST_COORDS_ROW_COL, time_attribute="InversionTime", + warp_matrix=warp_matrix) for i in range(np.size(self.MASK_POI_TARGET, 0)): np.testing.assert_equal( target_image_stack.ROI_time_series[i].POI_mask, @@ -388,11 +394,13 @@ def test_extract_single_roi(self): template_dcm = pydicom.read_file(self.TEMPLATE_PATH_T1_P5) template_image_stack = hazen_relaxometry.T1ImageStack([template_dcm], - template_dcm) + time_attribute="InversionTime") # set warp_matrix to identity matrix # template_image_stack.warp_matrix = np.eye(2, 3, dtype=np.float32) + warp_matrix = template_image_stack.template_fit(template_dcm) template_image_stack.generate_time_series( - self.TEMPLATE_TEST_COORDS_ROW_COL, fit_coords=False) + self.TEMPLATE_TEST_COORDS_ROW_COL, time_attribute="InversionTime", + warp_matrix=warp_matrix, fit_coords=False) np.testing.assert_equal( template_image_stack.ROI_time_series[0].pixel_values[0], @@ -403,10 +411,12 @@ def test_template_roi_means(self): template_dcm = pydicom.read_file(self.TEMPLATE_PATH_T1_P5) template_image_stack = hazen_relaxometry.T1ImageStack([template_dcm], - template_dcm) + time_attribute="InversionTime") + warp_matrix = template_image_stack.template_fit(template_dcm) template_image_stack.generate_time_series( - self.TEMPLATE_TEST_COORDS_ROW_COL, fit_coords=False) + self.TEMPLATE_TEST_COORDS_ROW_COL, time_attribute="InversionTime", + warp_matrix=warp_matrix, fit_coords=False) # Check all pixels in ROI[0] match np.testing.assert_allclose( @@ -424,18 +434,20 @@ def test_t1_calc_magnitude_image(self): template_dcm = pydicom.read_file(self.TEMPLATE_PATH_T1_P5) t1_dcms = [pydicom.dcmread(os.path.join(self.T1_DIR, fname)) for fname in self.T1_FILES] - t1_image_stack = hazen_relaxometry.T1ImageStack(t1_dcms, template_dcm) - t1_image_stack.template_fit() + t1_image_stack = hazen_relaxometry.T1ImageStack(t1_dcms, + time_attribute="InversionTime") + warp_matrix = t1_image_stack.template_fit(template_dcm) t1_image_stack.generate_time_series( - self.TEMPLATE_TEST_COORDS_ROW_COL, fit_coords=True) + self.TEMPLATE_TEST_COORDS_ROW_COL, time_attribute="InversionTime", + warp_matrix=warp_matrix) t1_image_stack.generate_fit_function() - t1_published = \ - hazen_relaxometry.TEMPLATE_VALUES['plate5']['t1']['relax_times']['1.5T'] - t1_image_stack.initialise_fit_parameters(t1_estimates=t1_published) + t1_published = TEMPLATE_VALUES['plate5']['t1']['relax_times']['1.5T'] + s0_est = t1_image_stack.initialise_fit_parameters(t1_estimates=t1_published) - t1_image_stack.find_relax_times() + t1_image_stack.find_relax_times( + t1_estimates=t1_published, s0_est=s0_est) - np.testing.assert_allclose(t1_image_stack.t1s, self.PLATE5_T1, + np.testing.assert_allclose(t1_image_stack.relax_times, self.PLATE5_T1, rtol=0.02, atol=1) def test_t2_calc_magnitude_image(self): @@ -461,17 +473,19 @@ def test_t1_calc_signed_image(self): template_dcm = pydicom.read_file(self.TEMPLATE_PATH_T1_P5) t1_dcms = [pydicom.dcmread(os.path.join(self.SITE2_T1_DIR, fname)) for fname in self.SITE2_T1_FILES] - t1_image_stack = hazen_relaxometry.T1ImageStack(t1_dcms, template_dcm) - t1_image_stack.template_fit() + t1_image_stack = hazen_relaxometry.T1ImageStack( + t1_dcms, time_attribute="InversionTime") + warp_matrix = t1_image_stack.template_fit(template_dcm) t1_image_stack.generate_time_series( - self.TEMPLATE_TEST_COORDS_ROW_COL, fit_coords=True) + self.TEMPLATE_TEST_COORDS_ROW_COL, time_attribute="InversionTime", + warp_matrix=warp_matrix) t1_image_stack.generate_fit_function() - t1_published = \ - hazen_relaxometry.TEMPLATE_VALUES['plate5']['t1']['relax_times']['1.5T'] - t1_image_stack.initialise_fit_parameters(t1_estimates=t1_published) - t1_image_stack.find_relax_times() + t1_published = TEMPLATE_VALUES['plate5']['t1']['relax_times']['1.5T'] + s0_est = t1_image_stack.initialise_fit_parameters(t1_estimates=t1_published) + t1_image_stack.find_relax_times( + t1_estimates=t1_published, s0_est=s0_est) - np.testing.assert_allclose(t1_image_stack.t1s, self.SITE2_PLATE5_T1, + np.testing.assert_allclose(t1_image_stack.relax_times, self.SITE2_PLATE5_T1, rtol=0.02, atol=1) def test_t1_siemens(self): @@ -535,17 +549,16 @@ def test_t2_p5_philips(self): def test_scale_up_template(self): """Test fit for 256x256 GE image with 192x192 template""" template_dcm = pydicom.read_file( - hazen_relaxometry.TEMPLATE_VALUES['plate4']['t1']['filename']) + TEMPLATE_VALUES['plate4']['t1']['filename']) target_dcm = pydicom.dcmread(self.PATH_256_MATRIX) t1_image_stack = hazen_relaxometry.T1ImageStack([target_dcm], - template_dcm) - t1_image_stack.template_fit() + time_attribute="InversionTime") + warp_matrix = t1_image_stack.template_fit(template_dcm) transformed_coordinates_xy = hazen_relaxometry.transform_coords( - hazen_relaxometry.TEMPLATE_VALUES['plate4']['sphere_centres_row_col'], - t1_image_stack.warp_matrix, - input_row_col=True, output_row_col=True) + TEMPLATE_VALUES['plate4']['sphere_centres_row_col'], + warp_matrix, input_row_col=True, output_row_col=True) # test to within +/- 1 pixel (also checks YX-XY change) np.testing.assert_allclose( @@ -555,7 +568,7 @@ def test_scale_up_template(self): def test_ge(self): """Test relaxometry.py values on GE.""" for plate in (4, 5): - for tparam in ('T1', 'T2'): + for tparam in ['T1']: # , 'T2' dcms = [pydicom.dcmread(os.path.join( getattr(self, f'SITE3_{tparam}_P{plate}_DIR'), fname)) for fname in getattr(self, f'SITE3_{tparam}_P{plate}_FILES')] From fed68af167ace16b498d557eecc1a0f2c3d503ea Mon Sep 17 00:00:00 2001 From: sophie22 Date: Wed, 26 Jul 2023 01:07:24 +0100 Subject: [PATCH 12/17] simplify importing relaxometry functions --- tests/test_relaxometry.py | 90 ++++++++++++++++++++------------------- 1 file changed, 47 insertions(+), 43 deletions(-) diff --git a/tests/test_relaxometry.py b/tests/test_relaxometry.py index e9164c66..96d9e0a8 100644 --- a/tests/test_relaxometry.py +++ b/tests/test_relaxometry.py @@ -11,7 +11,8 @@ import os.path from pydicom.errors import InvalidDicomError -import hazenlib.relaxometry as hazen_relaxometry +from hazenlib.relaxometry import ( + transform_coords, T1ImageStack, T2ImageStack, main) from hazenlib.exceptions import ArgumentCombinationError from tests import TEST_DATA_DIR, TEST_REPORT_DIR from hazenlib.config import TEMPLATE_VALUES @@ -250,7 +251,7 @@ class TestRelaxometry(unittest.TestCase): def test_transform_coords(self): # no translation, no rotation, input = yx, output = yx warp_matrix = np.array([[1, 0, 0], [0, 1, 0]]) - op = hazen_relaxometry.transform_coords(self.TEST_COORDS, + op = transform_coords(self.TEST_COORDS, warp_matrix, input_row_col=True, output_row_col=True) @@ -260,7 +261,7 @@ def test_transform_coords(self): # no translation, no rotation, flip # input = col_row (xy), output = row_col (yx) - op = hazen_relaxometry.transform_coords(self.TEST_COORDS, warp_matrix, + op = transform_coords(self.TEST_COORDS, warp_matrix, input_row_col=False, output_row_col=True) @@ -268,7 +269,7 @@ def test_transform_coords(self): # no translation, no rotation, input = col_row (xy), # output = col_row (xy) - op = hazen_relaxometry.transform_coords(self.TEST_COORDS, warp_matrix, + op = transform_coords(self.TEST_COORDS, warp_matrix, input_row_col=False, output_row_col=False) np.testing.assert_allclose(op, self.TEST_COORDS) @@ -277,7 +278,7 @@ def test_transform_coords(self): # translation x=1, y=3, no rotation, input = row_col (yx), # output = row_col (yx) warp_matrix = np.array([[1, 0, 1], [0, 1, 3]]) - op = hazen_relaxometry.transform_coords(self.TEST_COORDS, warp_matrix, + op = transform_coords(self.TEST_COORDS, warp_matrix, input_row_col=True, output_row_col=True) np.testing.assert_allclose(op, self.COORDS_TRANS) @@ -286,7 +287,7 @@ def test_transform_coords(self): # translation x=1, y=3, no rotation, input = col_row (xy), # output = row_col (yx) warp_matrix = np.array([[1, 0, 1], [0, 1, 3]]) - op = hazen_relaxometry.transform_coords(self.TEST_COORDS, warp_matrix, + op = transform_coords(self.TEST_COORDS, warp_matrix, input_row_col=False, output_row_col=True) np.testing.assert_allclose(op, self.COORDS_TRANS_FLIP) @@ -295,7 +296,7 @@ def test_transform_coords(self): # translation x=1, y=3, no rotation, input = col_row (xy), # output = col_row (xy) warp_matrix = np.array([[1, 0, 1], [0, 1, 3]]) - op = hazen_relaxometry.transform_coords(self.TEST_COORDS, warp_matrix, + op = transform_coords(self.TEST_COORDS, warp_matrix, input_row_col=False, output_row_col=False) np.testing.assert_allclose(op, self.COORDS_TRANS_COL_ROW) @@ -306,7 +307,7 @@ def test_transform_coords(self): warp_matrix = np.array([[np.sqrt(3) / 2, 0.5, 10], [-0.5, np.sqrt(3) / 2, 20]]) # use float64 rather than int32 for coordinates to better test rotation - op = hazen_relaxometry.transform_coords(np.array(self.TEST_COORDS, + op = transform_coords(np.array(self.TEST_COORDS, dtype=np.float64), warp_matrix, input_row_col=False, @@ -317,11 +318,11 @@ def test_template_fit(self): template_dcm = pydicom.read_file(self.TEMPLATE_PATH_T1_P5) target_dcm = pydicom.dcmread(self.TEMPLATE_TARGET_PATH_T1_P5) - t1_image_stack = hazen_relaxometry.T1ImageStack([target_dcm], + t1_image_stack = T1ImageStack([target_dcm], time_attribute="InversionTime") warp_matrix = t1_image_stack.template_fit(template_dcm) - transformed_coordinates_xy = hazen_relaxometry.transform_coords( + transformed_coordinates_xy = transform_coords( self.TEMPLATE_TEST_COORDS_ROW_COL, warp_matrix, input_row_col=True, output_row_col=False) @@ -334,8 +335,8 @@ def test_image_stack_T1_sort(self): # read list of un-ordered T1 files, sort by TI, test sorted t1_dcms = [pydicom.dcmread(os.path.join(self.T1_DIR, fname)) for fname in self.T1_FILES] - t1_image_stack = hazen_relaxometry.T1ImageStack(t1_dcms, - time_attribute="InversionTime") + t1_image_stack = T1ImageStack( + t1_dcms, time_attribute="InversionTime") sorted_output = [image.InversionTime.real for image in t1_image_stack.images] assert sorted_output == self.T1_TI_SORTED @@ -344,8 +345,8 @@ def test_image_stack_T2_sort(self): # read list of un-ordered T2 files, sort by TE, test sorted t2_dcms = [pydicom.dcmread(os.path.join(self.T2_DIR, fname)) for fname in self.T2_FILES] - t2_image_stack = hazen_relaxometry.T2ImageStack(t2_dcms, - time_attribute="EchoTime") + t2_image_stack = T2ImageStack( + t2_dcms, time_attribute="EchoTime") sorted_output = [image.EchoTime.real for image in t2_image_stack.images] @@ -356,7 +357,7 @@ def test_generate_time_series_template_POIs(self): # Test on template first, no image fitting needed # Need image to get correct size template_dcm = pydicom.dcmread(self.TEMPLATE_PATH_T1_P5) - template_image_stack = hazen_relaxometry.T1ImageStack([template_dcm], + template_image_stack = T1ImageStack([template_dcm], time_attribute="InversionTime") warp_matrix = template_image_stack.template_fit(template_dcm) template_image_stack.generate_time_series( @@ -373,10 +374,10 @@ def test_generate_time_series_target_POIs(self): template_dcm = pydicom.read_file(self.TEMPLATE_PATH_T1_P5) target_dcm = pydicom.dcmread(self.TEMPLATE_TARGET_PATH_T1_P5) - target_image_stack = hazen_relaxometry.T1ImageStack([target_dcm], + target_image_stack = T1ImageStack([target_dcm], time_attribute="InversionTime") warp_matrix = target_image_stack.template_fit(template_dcm) - # transformed_coordinates_yx = hazen_relaxometry.transform_coords( + # transformed_coordinates_yx = transform_coords( # self.TEMPLATE_TEST_COORDS_YX, target_image_stack.warp_matrix, # input_yx=True, output_yx=True) target_image_stack.generate_time_series( @@ -393,7 +394,7 @@ def test_extract_single_roi(self): # fitting. template_dcm = pydicom.read_file(self.TEMPLATE_PATH_T1_P5) - template_image_stack = hazen_relaxometry.T1ImageStack([template_dcm], + template_image_stack = T1ImageStack([template_dcm], time_attribute="InversionTime") # set warp_matrix to identity matrix # template_image_stack.warp_matrix = np.eye(2, 3, dtype=np.float32) @@ -410,7 +411,7 @@ def test_template_roi_means(self): # Check mean of first 3 ROIs in template match with ImageJ calculations template_dcm = pydicom.read_file(self.TEMPLATE_PATH_T1_P5) - template_image_stack = hazen_relaxometry.T1ImageStack([template_dcm], + template_image_stack = T1ImageStack([template_dcm], time_attribute="InversionTime") warp_matrix = template_image_stack.template_fit(template_dcm) @@ -434,7 +435,7 @@ def test_t1_calc_magnitude_image(self): template_dcm = pydicom.read_file(self.TEMPLATE_PATH_T1_P5) t1_dcms = [pydicom.dcmread(os.path.join(self.T1_DIR, fname)) for fname in self.T1_FILES] - t1_image_stack = hazen_relaxometry.T1ImageStack(t1_dcms, + t1_image_stack = T1ImageStack(t1_dcms, time_attribute="InversionTime") warp_matrix = t1_image_stack.template_fit(template_dcm) t1_image_stack.generate_time_series( @@ -455,15 +456,17 @@ def test_t2_calc_magnitude_image(self): template_dcm = pydicom.read_file(self.TEMPLATE_PATH_T2) t2_dcms = [pydicom.dcmread(os.path.join(self.T2_DIR, fname)) for fname in self.T2_FILES] - t2_image_stack = hazen_relaxometry.T2ImageStack(t2_dcms, template_dcm) - t2_image_stack.template_fit() + t2_image_stack = T2ImageStack( + t2_dcms, time_attribute="EchoTime") + warp_matrix = t2_image_stack.template_fit(template_dcm) t2_image_stack.generate_time_series( - self.TEMPLATE_P4_TEST_COORDS_ROW_COL, fit_coords=True) - t2_published = \ - hazen_relaxometry.TEMPLATE_VALUES['plate4']['t2']['relax_times']['1.5T'] - t2_image_stack.initialise_fit_parameters(t2_estimates=t2_published) - t2_image_stack.initialise_fit_parameters(t2_published) - t2_image_stack.find_relax_times() + self.TEMPLATE_P4_TEST_COORDS_ROW_COL, time_attribute="EchoTime", + warp_matrix=warp_matrix) + t2_published = TEMPLATE_VALUES['plate4']['t2']['relax_times']['1.5T'] + s0_est = t2_image_stack.initialise_fit_parameters( + t2_estimates=t2_published) + t2_image_stack.find_relax_times( + t2_estimates=t2_published, s0_est=s0_est) np.testing.assert_allclose(t2_image_stack.t2s, self.PLATE4_T2, rtol=0.01, atol=1) @@ -473,7 +476,7 @@ def test_t1_calc_signed_image(self): template_dcm = pydicom.read_file(self.TEMPLATE_PATH_T1_P5) t1_dcms = [pydicom.dcmread(os.path.join(self.SITE2_T1_DIR, fname)) for fname in self.SITE2_T1_FILES] - t1_image_stack = hazen_relaxometry.T1ImageStack( + t1_image_stack = T1ImageStack( t1_dcms, time_attribute="InversionTime") warp_matrix = t1_image_stack.template_fit(template_dcm) t1_image_stack.generate_time_series( @@ -481,7 +484,8 @@ def test_t1_calc_signed_image(self): warp_matrix=warp_matrix) t1_image_stack.generate_fit_function() t1_published = TEMPLATE_VALUES['plate5']['t1']['relax_times']['1.5T'] - s0_est = t1_image_stack.initialise_fit_parameters(t1_estimates=t1_published) + s0_est = t1_image_stack.initialise_fit_parameters( + t1_estimates=t1_published) t1_image_stack.find_relax_times( t1_estimates=t1_published, s0_est=s0_est) @@ -492,7 +496,7 @@ def test_t1_siemens(self): """Test T1 values on Siemens images.""" dcms = [pydicom.dcmread(os.path.join(self.T1_DIR, fname)) for fname in self.T1_FILES] - t1_results = hazen_relaxometry.main(dcms, plate_number=5, calc="T1", + t1_results = main(dcms, plate_number=5, calc="T1", verbose=True, report_dir="report/relaxometry/") # `t1_results` is a dict with one item where we don't know the key. # Need to extract via unpacking @@ -504,7 +508,7 @@ def test_t1_p4_philips(self): """Test T1 values on plate 4 on Philips.""" dcms = [pydicom.dcmread(os.path.join(self.SITE4_T1_P4_DIR, fname)) for fname in self.SITE4_T1_P4_FILES] - t1_results = hazen_relaxometry.main(dcms, plate_number=4, calc="T1", + t1_results = main(dcms, plate_number=4, calc="T1", verbose=True, report_dir="report/relaxometry/") # `t1_results` is a dict with one item where we don't know the key. # Need to extract via unpacking @@ -517,7 +521,7 @@ def test_t1_p5_philips(self): """Test T1 values on plate 5 on Philips.""" dcms = [pydicom.dcmread(os.path.join(self.SITE4_T1_P5_DIR, fname)) for fname in self.SITE4_T1_P5_FILES] - t1_results = hazen_relaxometry.main(dcms, plate_number=5, calc="T1", + t1_results = main(dcms, plate_number=5, calc="T1", verbose=True, report_dir="report/relaxometry/") results, = t1_results.values() np.testing.assert_allclose(results['calc_times'], @@ -528,7 +532,7 @@ def test_t2_p4_philips(self): """Test T2 values on plate 4 on Philips.""" dcms = [pydicom.dcmread(os.path.join(self.SITE4_T2_P4_DIR, fname)) for fname in self.SITE4_T2_P4_FILES] - t2_results = hazen_relaxometry.main(dcms, plate_number=4, calc="T2", + t2_results = main(dcms, plate_number=4, calc="T2", verbose=True, report_dir="report/relaxometry/") results, = t2_results.values() np.testing.assert_allclose(results['calc_times'], @@ -539,7 +543,7 @@ def test_t2_p5_philips(self): """Test T2 values on plate 4 on Philips.""" dcms = [pydicom.dcmread(os.path.join(self.SITE4_T2_P5_DIR, fname)) for fname in self.SITE4_T2_P5_FILES] - t2_results = hazen_relaxometry.main(dcms, plate_number=5, calc="T2", + t2_results = main(dcms, plate_number=5, calc="T2", verbose=True, report_dir="report/relaxometry/") results, = t2_results.values() np.testing.assert_allclose(results['calc_times'], @@ -552,11 +556,11 @@ def test_scale_up_template(self): TEMPLATE_VALUES['plate4']['t1']['filename']) target_dcm = pydicom.dcmread(self.PATH_256_MATRIX) - t1_image_stack = hazen_relaxometry.T1ImageStack([target_dcm], + t1_image_stack = T1ImageStack([target_dcm], time_attribute="InversionTime") warp_matrix = t1_image_stack.template_fit(template_dcm) - transformed_coordinates_xy = hazen_relaxometry.transform_coords( + transformed_coordinates_xy = transform_coords( TEMPLATE_VALUES['plate4']['sphere_centres_row_col'], warp_matrix, input_row_col=True, output_row_col=True) @@ -573,7 +577,7 @@ def test_ge(self): getattr(self, f'SITE3_{tparam}_P{plate}_DIR'), fname)) for fname in getattr(self, f'SITE3_{tparam}_P{plate}_FILES')] - t_results = hazen_relaxometry.main(dcms, plate_number=plate, + t_results = main(dcms, plate_number=plate, calc = tparam, verbose=True, report_dir="report/relaxometry/") results, = t_results.values() @@ -585,7 +589,7 @@ def test_ge(self): def test_plate_number_not_specified(self): """Test exception raised if plate_number not specified.""" self.assertRaises(ArgumentCombinationError, - hazen_relaxometry.main, [], calc="T1") + main, [], calc="T1") def test_philips_3T(self): """Test calculation on 3T dataset.""" @@ -593,7 +597,7 @@ def test_philips_3T(self): # T1 plate 4 dcms = [pydicom.dcmread(os.path.join(self.SITE5_T1_P4_DIR, fname)) for fname in self.SITE5_T1_P4_FILES] - t1_results = hazen_relaxometry.main(dcms, plate_number=4, calc="T1", + t1_results = main(dcms, plate_number=4, calc="T1", verbose=True, report_dir="report/relaxometry/") results, = t1_results.values() np.testing.assert_allclose(results['calc_times'], @@ -603,7 +607,7 @@ def test_philips_3T(self): # T1 plate 5 dcms = [pydicom.dcmread(os.path.join(self.SITE5_T1_P5_DIR, fname)) for fname in self.SITE5_T1_P5_FILES] - t1_results = hazen_relaxometry.main(dcms, plate_number=5, calc="T1", + t1_results = main(dcms, plate_number=5, calc="T1", verbose=True, report_dir="report/relaxometry/") results, = t1_results.values() np.testing.assert_allclose(results['calc_times'], @@ -613,7 +617,7 @@ def test_philips_3T(self): # T2 plate 4 dcms = [pydicom.dcmread(os.path.join(self.SITE5_T2_P4_DIR, fname)) for fname in self.SITE5_T2_P4_FILES] - t2_results = hazen_relaxometry.main(dcms, plate_number=4, calc="T2", + t2_results = main(dcms, plate_number=4, calc="T2", verbose=True, report_dir="report/relaxometry/") results, = t2_results.values() np.testing.assert_allclose(results['calc_times'], @@ -623,7 +627,7 @@ def test_philips_3T(self): # T2 plate 5 dcms = [pydicom.dcmread(os.path.join(self.SITE5_T2_P5_DIR, fname)) for fname in self.SITE5_T2_P5_FILES] - t2_results = hazen_relaxometry.main(dcms, plate_number=5, calc="T2", + t2_results = main(dcms, plate_number=5, calc="T2", verbose=True, report_dir="report/relaxometry/") results, = t2_results.values() np.testing.assert_allclose(results['calc_times'], From f5f081747cd037762a4356c4094911cb84dbfa51 Mon Sep 17 00:00:00 2001 From: sophie22 Date: Wed, 26 Jul 2023 01:14:42 +0100 Subject: [PATCH 13/17] add back the fit_coords options --- hazenlib/relaxometry.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/hazenlib/relaxometry.py b/hazenlib/relaxometry.py index 02a9aa5c..566fc22d 100644 --- a/hazenlib/relaxometry.py +++ b/hazenlib/relaxometry.py @@ -523,7 +523,7 @@ def template_fit(self, template_dcm, image_index=0): return warp_matrix def generate_time_series(self, coords_row_col, time_attribute, - warp_matrix): + warp_matrix, fit_coords=True): """ Create list of ROITimeSeries objects. @@ -540,13 +540,18 @@ def generate_time_series(self, coords_row_col, time_attribute, RT transform matrix (2,3). """ num_coords = np.size(coords_row_col, axis=0) - flipped_coords_row_col = transform_coords(coords_row_col, warp_matrix, - input_row_col=True, output_row_col=True) + + # adjustment may not be required for the template DICOM + if fit_coords: + adjusted_coords_row_col = transform_coords(coords_row_col, + warp_matrix, input_row_col=True, output_row_col=True) + else: # used in testing + adjusted_coords_row_col = coords_row_col self.ROI_time_series = [] for i in range(num_coords): self.ROI_time_series.append(ROITimeSeries( - self.images, flipped_coords_row_col[i], time_attribute)) + self.images, adjusted_coords_row_col[i], time_attribute)) def plot_fit(self): """ @@ -904,7 +909,7 @@ def initialise_fit_parameters(self, t2_estimates): Returns ------- - None. + s0_est. """ self.t2_est = t2_estimates From a00793a8349958af70792625e4cce490daebdc7b Mon Sep 17 00:00:00 2001 From: sophie22 Date: Wed, 26 Jul 2023 01:15:02 +0100 Subject: [PATCH 14/17] correct accessing a renamed property --- tests/test_relaxometry.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_relaxometry.py b/tests/test_relaxometry.py index 96d9e0a8..4947cd88 100644 --- a/tests/test_relaxometry.py +++ b/tests/test_relaxometry.py @@ -468,7 +468,7 @@ def test_t2_calc_magnitude_image(self): t2_image_stack.find_relax_times( t2_estimates=t2_published, s0_est=s0_est) - np.testing.assert_allclose(t2_image_stack.t2s, self.PLATE4_T2, + np.testing.assert_allclose(t2_image_stack.relax_times, self.PLATE4_T2, rtol=0.01, atol=1) def test_t1_calc_signed_image(self): From 59dee22882bd997b524bcc83e6d28d9483c21fe1 Mon Sep 17 00:00:00 2001 From: sophie22 Date: Wed, 26 Jul 2023 01:24:58 +0100 Subject: [PATCH 15/17] set time_attribute at class definition --- hazenlib/relaxometry.py | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/hazenlib/relaxometry.py b/hazenlib/relaxometry.py index 566fc22d..8367be1d 100644 --- a/hazenlib/relaxometry.py +++ b/hazenlib/relaxometry.py @@ -413,6 +413,7 @@ def __init__(self, image_slices, time_attribute=None): DICOM attribute to order images. Typically 'InversionTime' for T1 relaxometry or 'EchoTime' for T2. """ + self.time_attr = time_attribute # store sorted images self.images = self.order_by(image_slices, time_attribute) @@ -522,8 +523,7 @@ def template_fit(self, template_dcm, image_index=0): return warp_matrix - def generate_time_series(self, coords_row_col, time_attribute, - warp_matrix, fit_coords=True): + def generate_time_series(self, coords_row_col, warp_matrix, fit_coords=True): """ Create list of ROITimeSeries objects. @@ -551,7 +551,7 @@ def generate_time_series(self, coords_row_col, time_attribute, self.ROI_time_series = [] for i in range(num_coords): self.ROI_time_series.append(ROITimeSeries( - self.images, adjusted_coords_row_col[i], time_attribute)) + self.images, adjusted_coords_row_col[i], self.time_attr)) def plot_fit(self): """ @@ -630,7 +630,8 @@ class T1ImageStack(ImageStack): Calculate T1 relaxometry. """ - def __init__(self, image_slices, time_attribute): + def __init__(self, image_slices): + time_attribute = "InversionTime" super().__init__(image_slices, time_attribute) def generate_t1_function(self, ti_interp_vals, tr_interp_vals, mag_image=False): @@ -811,7 +812,8 @@ class T2ImageStack(ImageStack): Calculate T2 relaxometry. """ - def __init__(self, image_slices, time_attribute): + def __init__(self, image_slices): + time_attribute = "EchoTime" super().__init__(image_slices, time_attribute) self.fit_eqn_str = 'T2 with Rician noise (Raya et al 2010)' @@ -1026,13 +1028,9 @@ def main(data, plate_number=None, calc: str = 'T1', verbose=False, # Set up parameters specific to T1 or T2 relax_str = calc.lower() - time_attributes = { - "t1": 'InversionTime', - "t2": 'EchoTime' - } if calc in ['T1', 't1']: - image_stack = T1ImageStack(data, time_attribute=time_attributes[relax_str]) + image_stack = T1ImageStack(data) try: template_dcm = pydicom.read_file( TEMPLATE_VALUES[f'plate{plate_number}'][relax_str]['filename']) @@ -1041,7 +1039,7 @@ def main(data, plate_number=None, calc: str = 'T1', verbose=False, f' Please pass plate number as arg.') exit() elif calc in ['T2', 't2']: - image_stack = T2ImageStack(data, time_attribute=time_attributes[relax_str]) + image_stack = T2ImageStack(data) try: template_dcm = pydicom.read_file( TEMPLATE_VALUES[f'plate{plate_number}'][relax_str]['filename']) @@ -1056,7 +1054,7 @@ def main(data, plate_number=None, calc: str = 'T1', verbose=False, warp_matrix = image_stack.template_fit(template_dcm) image_stack.generate_time_series( TEMPLATE_VALUES[f'plate{plate_number}']['sphere_centres_row_col'], - time_attribute=time_attributes[relax_str], warp_matrix=warp_matrix) + warp_matrix=warp_matrix) # only applies to T1 image_stack.generate_fit_function() From f9ebad79ddfe3a303174fad523bd2c27cd7e84eb Mon Sep 17 00:00:00 2001 From: sophie22 Date: Wed, 26 Jul 2023 01:30:06 +0100 Subject: [PATCH 16/17] remove time_attribute args --- tests/test_relaxometry.py | 47 +++++++++++++++------------------------ 1 file changed, 18 insertions(+), 29 deletions(-) diff --git a/tests/test_relaxometry.py b/tests/test_relaxometry.py index 4947cd88..670edb4c 100644 --- a/tests/test_relaxometry.py +++ b/tests/test_relaxometry.py @@ -318,8 +318,7 @@ def test_template_fit(self): template_dcm = pydicom.read_file(self.TEMPLATE_PATH_T1_P5) target_dcm = pydicom.dcmread(self.TEMPLATE_TARGET_PATH_T1_P5) - t1_image_stack = T1ImageStack([target_dcm], - time_attribute="InversionTime") + t1_image_stack = T1ImageStack([target_dcm]) warp_matrix = t1_image_stack.template_fit(template_dcm) transformed_coordinates_xy = transform_coords( @@ -335,8 +334,7 @@ def test_image_stack_T1_sort(self): # read list of un-ordered T1 files, sort by TI, test sorted t1_dcms = [pydicom.dcmread(os.path.join(self.T1_DIR, fname)) for fname in self.T1_FILES] - t1_image_stack = T1ImageStack( - t1_dcms, time_attribute="InversionTime") + t1_image_stack = T1ImageStack(t1_dcms) sorted_output = [image.InversionTime.real for image in t1_image_stack.images] assert sorted_output == self.T1_TI_SORTED @@ -345,8 +343,7 @@ def test_image_stack_T2_sort(self): # read list of un-ordered T2 files, sort by TE, test sorted t2_dcms = [pydicom.dcmread(os.path.join(self.T2_DIR, fname)) for fname in self.T2_FILES] - t2_image_stack = T2ImageStack( - t2_dcms, time_attribute="EchoTime") + t2_image_stack = T2ImageStack(t2_dcms) sorted_output = [image.EchoTime.real for image in t2_image_stack.images] @@ -357,11 +354,10 @@ def test_generate_time_series_template_POIs(self): # Test on template first, no image fitting needed # Need image to get correct size template_dcm = pydicom.dcmread(self.TEMPLATE_PATH_T1_P5) - template_image_stack = T1ImageStack([template_dcm], - time_attribute="InversionTime") + template_image_stack = T1ImageStack([template_dcm]) warp_matrix = template_image_stack.template_fit(template_dcm) template_image_stack.generate_time_series( - self.TEMPLATE_TEST_COORDS_ROW_COL, time_attribute="InversionTime", + self.TEMPLATE_TEST_COORDS_ROW_COL, warp_matrix=warp_matrix, fit_coords=False) for i in range(np.size(self.MASK_POI_TEMPLATE, 0)): @@ -374,14 +370,13 @@ def test_generate_time_series_target_POIs(self): template_dcm = pydicom.read_file(self.TEMPLATE_PATH_T1_P5) target_dcm = pydicom.dcmread(self.TEMPLATE_TARGET_PATH_T1_P5) - target_image_stack = T1ImageStack([target_dcm], - time_attribute="InversionTime") + target_image_stack = T1ImageStack([target_dcm]) warp_matrix = target_image_stack.template_fit(template_dcm) # transformed_coordinates_yx = transform_coords( # self.TEMPLATE_TEST_COORDS_YX, target_image_stack.warp_matrix, # input_yx=True, output_yx=True) target_image_stack.generate_time_series( - self.TEMPLATE_TEST_COORDS_ROW_COL, time_attribute="InversionTime", + self.TEMPLATE_TEST_COORDS_ROW_COL, warp_matrix=warp_matrix) for i in range(np.size(self.MASK_POI_TARGET, 0)): np.testing.assert_equal( @@ -394,13 +389,12 @@ def test_extract_single_roi(self): # fitting. template_dcm = pydicom.read_file(self.TEMPLATE_PATH_T1_P5) - template_image_stack = T1ImageStack([template_dcm], - time_attribute="InversionTime") + template_image_stack = T1ImageStack([template_dcm]) # set warp_matrix to identity matrix # template_image_stack.warp_matrix = np.eye(2, 3, dtype=np.float32) warp_matrix = template_image_stack.template_fit(template_dcm) template_image_stack.generate_time_series( - self.TEMPLATE_TEST_COORDS_ROW_COL, time_attribute="InversionTime", + self.TEMPLATE_TEST_COORDS_ROW_COL, warp_matrix=warp_matrix, fit_coords=False) np.testing.assert_equal( @@ -411,12 +405,11 @@ def test_template_roi_means(self): # Check mean of first 3 ROIs in template match with ImageJ calculations template_dcm = pydicom.read_file(self.TEMPLATE_PATH_T1_P5) - template_image_stack = T1ImageStack([template_dcm], - time_attribute="InversionTime") + template_image_stack = T1ImageStack([template_dcm]) warp_matrix = template_image_stack.template_fit(template_dcm) template_image_stack.generate_time_series( - self.TEMPLATE_TEST_COORDS_ROW_COL, time_attribute="InversionTime", + self.TEMPLATE_TEST_COORDS_ROW_COL, warp_matrix=warp_matrix, fit_coords=False) # Check all pixels in ROI[0] match @@ -435,11 +428,10 @@ def test_t1_calc_magnitude_image(self): template_dcm = pydicom.read_file(self.TEMPLATE_PATH_T1_P5) t1_dcms = [pydicom.dcmread(os.path.join(self.T1_DIR, fname)) for fname in self.T1_FILES] - t1_image_stack = T1ImageStack(t1_dcms, - time_attribute="InversionTime") + t1_image_stack = T1ImageStack(t1_dcms) warp_matrix = t1_image_stack.template_fit(template_dcm) t1_image_stack.generate_time_series( - self.TEMPLATE_TEST_COORDS_ROW_COL, time_attribute="InversionTime", + self.TEMPLATE_TEST_COORDS_ROW_COL, warp_matrix=warp_matrix) t1_image_stack.generate_fit_function() t1_published = TEMPLATE_VALUES['plate5']['t1']['relax_times']['1.5T'] @@ -456,11 +448,10 @@ def test_t2_calc_magnitude_image(self): template_dcm = pydicom.read_file(self.TEMPLATE_PATH_T2) t2_dcms = [pydicom.dcmread(os.path.join(self.T2_DIR, fname)) for fname in self.T2_FILES] - t2_image_stack = T2ImageStack( - t2_dcms, time_attribute="EchoTime") + t2_image_stack = T2ImageStack(t2_dcms) warp_matrix = t2_image_stack.template_fit(template_dcm) t2_image_stack.generate_time_series( - self.TEMPLATE_P4_TEST_COORDS_ROW_COL, time_attribute="EchoTime", + self.TEMPLATE_P4_TEST_COORDS_ROW_COL, warp_matrix=warp_matrix) t2_published = TEMPLATE_VALUES['plate4']['t2']['relax_times']['1.5T'] s0_est = t2_image_stack.initialise_fit_parameters( @@ -476,11 +467,10 @@ def test_t1_calc_signed_image(self): template_dcm = pydicom.read_file(self.TEMPLATE_PATH_T1_P5) t1_dcms = [pydicom.dcmread(os.path.join(self.SITE2_T1_DIR, fname)) for fname in self.SITE2_T1_FILES] - t1_image_stack = T1ImageStack( - t1_dcms, time_attribute="InversionTime") + t1_image_stack = T1ImageStack(t1_dcms) warp_matrix = t1_image_stack.template_fit(template_dcm) t1_image_stack.generate_time_series( - self.TEMPLATE_TEST_COORDS_ROW_COL, time_attribute="InversionTime", + self.TEMPLATE_TEST_COORDS_ROW_COL, warp_matrix=warp_matrix) t1_image_stack.generate_fit_function() t1_published = TEMPLATE_VALUES['plate5']['t1']['relax_times']['1.5T'] @@ -556,8 +546,7 @@ def test_scale_up_template(self): TEMPLATE_VALUES['plate4']['t1']['filename']) target_dcm = pydicom.dcmread(self.PATH_256_MATRIX) - t1_image_stack = T1ImageStack([target_dcm], - time_attribute="InversionTime") + t1_image_stack = T1ImageStack([target_dcm]) warp_matrix = t1_image_stack.template_fit(template_dcm) transformed_coordinates_xy = transform_coords( From 13eb88eee26781a1de3799903ca2fd7dfb67bef7 Mon Sep 17 00:00:00 2001 From: Tom Roberts Date: Fri, 28 Jul 2023 11:25:30 +0100 Subject: [PATCH 17/17] Changes filename --- hazenlib/relaxometry.py | 2 +- hazenlib/{config.py => relaxometry_params.py} | 0 tests/test_relaxometry.py | 2 +- 3 files changed, 2 insertions(+), 2 deletions(-) rename hazenlib/{config.py => relaxometry_params.py} (100%) diff --git a/hazenlib/relaxometry.py b/hazenlib/relaxometry.py index 8367be1d..56e286c2 100644 --- a/hazenlib/relaxometry.py +++ b/hazenlib/relaxometry.py @@ -136,7 +136,7 @@ from scipy.special import i0e, ive import hazenlib.exceptions -from hazenlib.config import ( +from hazenlib.relaxometry_params import ( MAX_RICIAN_NOISE, SEED_RICIAN_NOISE, TEMPLATE_VALUES, SMOOTH_TIMES, TEMPLATE_FIT_ITERS, TERMINATION_EPS ) diff --git a/hazenlib/config.py b/hazenlib/relaxometry_params.py similarity index 100% rename from hazenlib/config.py rename to hazenlib/relaxometry_params.py diff --git a/tests/test_relaxometry.py b/tests/test_relaxometry.py index 670edb4c..ee9c4e6b 100644 --- a/tests/test_relaxometry.py +++ b/tests/test_relaxometry.py @@ -15,7 +15,7 @@ transform_coords, T1ImageStack, T2ImageStack, main) from hazenlib.exceptions import ArgumentCombinationError from tests import TEST_DATA_DIR, TEST_REPORT_DIR -from hazenlib.config import TEMPLATE_VALUES +from hazenlib.relaxometry_params import TEMPLATE_VALUES class TestRelaxometry(unittest.TestCase):