From b043a13a50819866aadf4780f27024c670fb364a Mon Sep 17 00:00:00 2001 From: Casper da Costa-Luis Date: Tue, 9 Jul 2024 16:20:19 +0100 Subject: [PATCH 1/4] add metrics as per wiki --- petric.py | 76 +++++++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 60 insertions(+), 16 deletions(-) diff --git a/petric.py b/petric.py index c282e80..01c7c05 100755 --- a/petric.py +++ b/petric.py @@ -22,7 +22,7 @@ from traceback import print_exc import numpy as np -from skimage.metrics import mean_squared_error, peak_signal_noise_ratio +from skimage.metrics import mean_squared_error as mse from tensorboardX import SummaryWriter import sirf.STIR as STIR @@ -88,26 +88,52 @@ def __call__(self, algo: Algorithm): log.debug("...logged") +class QualityMetrics(ImageQualityCallback): + """From https://github.com/SyneRBI/PETRIC/wiki#metrics-and-thresholds""" + def __init__(self, reference_image, backround_mask, foreground_mask, **kwargs): + super().__init__(reference_image, **kwargs) + self.background = np.where(backround_mask == 1) + self.foreground = np.where(foreground_mask == 1) + + def __call__(self, algorithm): + iteration = algorithm.iteration + if iteration % algorithm.update_objective_interval != 0 and iteration != algorithm.max_iteration: + return + test_image = algorithm.x # CIL or SIRF ImageData + + test_im_arr, ref_im_arr = test_image.as_array(), self.reference_image.as_array() + + for filter_name, filter_func in self.filter.items(): + if filter_func is not None: + test_im, ref_im = map(filter_func, (test_im_arr, ref_im_arr)) + + # (1) global metrics & statistics + norm = ref_im[self.background].mean() + self.tb_summary_writer.add_scalar(f"RMSE_foreground{filter_name}", + np.sqrt(mse(test_im[self.foreground], ref_im[self.foreground])) / norm, + iteration) + self.tb_summary_writer.add_scalar(f"RMSE_background{filter_name}", + np.sqrt(mse(test_im[self.background], ref_im[self.background])) / norm, + iteration) + + # (2) local metrics & statistics + for roi_name, roi_inds in self.roi_indices.items(): + # AEM not to be confused with MAE + self.tb_summary_writer.add_scalar(f"AEM_VOI_{roi_name}{filter_name}", + np.abs(test_im[roi_inds].mean() - ref_im[roi_inds].mean()) / norm, + iteration) + + class MetricsWithTimeout(cbks.Callback): """Stops the algorithm after `seconds`""" - def __init__(self, seconds=300, outdir=OUTDIR, transverse_slice=None, coronal_slice=None, reference_image=None, - verbose=1): + def __init__(self, seconds=300, outdir=OUTDIR, transverse_slice=None, coronal_slice=None, verbose=1): super().__init__(verbose) self._seconds = seconds self.callbacks = [ cbks.ProgressCallback(), SaveIters(outdir=outdir), (tb_cbk := TensorBoard(logdir=outdir, transverse_slice=transverse_slice, coronal_slice=coronal_slice))] - - if reference_image: - roi_image_dict = {f'S{i}': STIR.ImageData(f'S{i}.hv') for i in range(1, 8)} - # NB: these metrics are for testing only. - # The final evaluation will use metrics described in https://github.com/SyneRBI/PETRIC/wiki - self.callbacks.append( - ImageQualityCallback( - reference_image, tb_cbk.tb, roi_mask_dict=roi_image_dict, metrics_dict={ - 'MSE': mean_squared_error, 'MAE': self.mean_absolute_error, 'PSNR': peak_signal_noise_ratio}, - statistics_dict={'MEAN': np.mean, 'STDDEV': np.std, 'MAX': np.max})) + self.tb = tb_cbk.tb self.reset() def reset(self, seconds=None): @@ -144,7 +170,9 @@ def construct_RDP(penalty_strength, initial_image, kappa, max_scaling=1e-3): return prior -Dataset = namedtuple('Dataset', ['acquired_data', 'additive_term', 'mult_factors', 'OSEM_image', 'prior', 'kappa']) +Dataset = namedtuple('Dataset', [ + 'acquired_data', 'additive_term', 'mult_factors', 'OSEM_image', 'prior', 'kappa', 'reference_image', + 'background_mask', 'foreground_mask', 'voi_masks']) def get_data(srcdir=".", outdir=OUTDIR, sirf_verbosity=0): @@ -165,7 +193,18 @@ def get_data(srcdir=".", outdir=OUTDIR, sirf_verbosity=0): penalty_strength = 1 / 700 # default choice prior = construct_RDP(penalty_strength, OSEM_image, kappa) - return Dataset(acquired_data, additive_term, mult_factors, OSEM_image, prior, kappa) + reference_image = STIR.ImageData(str(srcdir / 'reference_image.hv')) if (srcdir / + 'reference_image.hv').is_file() else None + background_mask = STIR.ImageData(str(srcdir / 'VOI_background.hv')) if (srcdir / + 'VOI_background.hv').is_file() else None + foreground_mask = STIR.ImageData(str(srcdir / 'VOI_foreground.hv')) if (srcdir / + 'VOI_foreground.hv').is_file() else None + voi_masks = { + voi.stem: STIR.ImageData(str(voi)) + for voi in srcdir.glob("VOI_*.hv") if voi.stem[4:] not in ('background', 'foreground')} + + return Dataset(acquired_data, additive_term, mult_factors, OSEM_image, prior, kappa, reference_image, + background_mask, foreground_mask, voi_masks) if SRCDIR.is_dir(): @@ -194,7 +233,12 @@ def get_data(srcdir=".", outdir=OUTDIR, sirf_verbosity=0): assert issubclass(Submission, Algorithm) for srcdir, outdir, metrics in data_dirs_metrics: data = get_data(srcdir=srcdir, outdir=outdir) - metrics[0].reset() # timeout from now + metrics_cbk = metrics[0] + if data.reference_image is not None: + metrics_cbk.callbacks.append( + QualityMetrics(data.reference_image, data.background_mask, data.foreground_mask, + tb_summary_writer=metrics_cbk.tb, roi_mask_dict=data.voi_masks)) + metrics_cbk.reset() # timeout from now algo = Submission(data) try: algo.run(np.inf, callbacks=metrics + submission_callbacks) From d397e0e162a67b48d79bd0dedb55191463cd3bc5 Mon Sep 17 00:00:00 2001 From: Casper da Costa-Luis Date: Tue, 9 Jul 2024 19:27:52 +0100 Subject: [PATCH 2/4] update after review --- petric.py | 58 +++++++++++++++++++++++++++---------------------------- 1 file changed, 29 insertions(+), 29 deletions(-) diff --git a/petric.py b/petric.py index 01c7c05..69d105c 100755 --- a/petric.py +++ b/petric.py @@ -90,38 +90,37 @@ def __call__(self, algo: Algorithm): class QualityMetrics(ImageQualityCallback): """From https://github.com/SyneRBI/PETRIC/wiki#metrics-and-thresholds""" - def __init__(self, reference_image, backround_mask, foreground_mask, **kwargs): + def __init__(self, reference_image, whole_object_mask, background_mask, **kwargs): super().__init__(reference_image, **kwargs) - self.background = np.where(backround_mask == 1) - self.foreground = np.where(foreground_mask == 1) + self.whole_object_indices = np.where(whole_object_mask == 1) + self.background_indices = np.where(background_mask == 1) + self.ref_im_arr = reference_image.as_array() def __call__(self, algorithm): iteration = algorithm.iteration if iteration % algorithm.update_objective_interval != 0 and iteration != algorithm.max_iteration: return - test_image = algorithm.x # CIL or SIRF ImageData - - test_im_arr, ref_im_arr = test_image.as_array(), self.reference_image.as_array() for filter_name, filter_func in self.filter.items(): - if filter_func is not None: - test_im, ref_im = map(filter_func, (test_im_arr, ref_im_arr)) + if filter_func is None: + filter_func = lambda x: x + test_im, ref_im = (filter_func(img_data).as_array() for img_data in (algorithm.x, self.reference_image)) # (1) global metrics & statistics - norm = ref_im[self.background].mean() - self.tb_summary_writer.add_scalar(f"RMSE_foreground{filter_name}", - np.sqrt(mse(test_im[self.foreground], ref_im[self.foreground])) / norm, - iteration) - self.tb_summary_writer.add_scalar(f"RMSE_background{filter_name}", - np.sqrt(mse(test_im[self.background], ref_im[self.background])) / norm, - iteration) + norm = ref_im[self.background_indices].mean() + self.tb_summary_writer.add_scalar( + f"RMSE_whole_object{filter_name}", + np.sqrt(mse(ref_im[self.whole_object_indices], test_im[self.whole_object_indices])) / norm, iteration) + self.tb_summary_writer.add_scalar( + f"RMSE_background{filter_name}", + np.sqrt(mse(ref_im[self.background_indices], test_im[self.background_indices])) / norm, iteration) # (2) local metrics & statistics - for roi_name, roi_inds in self.roi_indices.items(): + for voi_name, voi_indices in self.voi_indices.items(): # AEM not to be confused with MAE - self.tb_summary_writer.add_scalar(f"AEM_VOI_{roi_name}{filter_name}", - np.abs(test_im[roi_inds].mean() - ref_im[roi_inds].mean()) / norm, - iteration) + self.tb_summary_writer.add_scalar( + f"AEM_VOI_{voi_name}{filter_name}", + np.abs(test_im[voi_indices].mean() - ref_im[voi_indices].mean()) / norm, iteration) class MetricsWithTimeout(cbks.Callback): @@ -172,7 +171,7 @@ def construct_RDP(penalty_strength, initial_image, kappa, max_scaling=1e-3): Dataset = namedtuple('Dataset', [ 'acquired_data', 'additive_term', 'mult_factors', 'OSEM_image', 'prior', 'kappa', 'reference_image', - 'background_mask', 'foreground_mask', 'voi_masks']) + 'whole_object_mask', 'background_mask', 'voi_masks']) def get_data(srcdir=".", outdir=OUTDIR, sirf_verbosity=0): @@ -195,16 +194,17 @@ def get_data(srcdir=".", outdir=OUTDIR, sirf_verbosity=0): reference_image = STIR.ImageData(str(srcdir / 'reference_image.hv')) if (srcdir / 'reference_image.hv').is_file() else None + whole_object_mask = STIR.ImageData(str(srcdir / + 'VOI_whole_object.hv')) if (srcdir / + 'VOI_whole_object.hv').is_file() else None background_mask = STIR.ImageData(str(srcdir / 'VOI_background.hv')) if (srcdir / 'VOI_background.hv').is_file() else None - foreground_mask = STIR.ImageData(str(srcdir / 'VOI_foreground.hv')) if (srcdir / - 'VOI_foreground.hv').is_file() else None voi_masks = { voi.stem: STIR.ImageData(str(voi)) - for voi in srcdir.glob("VOI_*.hv") if voi.stem[4:] not in ('background', 'foreground')} + for voi in srcdir.glob("VOI_*.hv") if voi.stem[4:] not in ('background', 'whole_object')} return Dataset(acquired_data, additive_term, mult_factors, OSEM_image, prior, kappa, reference_image, - background_mask, foreground_mask, voi_masks) + whole_object_mask, background_mask, voi_masks) if SRCDIR.is_dir(): @@ -233,12 +233,12 @@ def get_data(srcdir=".", outdir=OUTDIR, sirf_verbosity=0): assert issubclass(Submission, Algorithm) for srcdir, outdir, metrics in data_dirs_metrics: data = get_data(srcdir=srcdir, outdir=outdir) - metrics_cbk = metrics[0] + metrics_with_timeout = metrics[0] if data.reference_image is not None: - metrics_cbk.callbacks.append( - QualityMetrics(data.reference_image, data.background_mask, data.foreground_mask, - tb_summary_writer=metrics_cbk.tb, roi_mask_dict=data.voi_masks)) - metrics_cbk.reset() # timeout from now + metrics_with_timeout.callbacks.append( + QualityMetrics(data.reference_image, data.whole_object_mask, data.background_mask, + tb_summary_writer=metrics_with_timeout.tb, roi_mask_dict=data.voi_masks)) + metrics_with_timeout.reset() # timeout from now algo = Submission(data) try: algo.run(np.inf, callbacks=metrics + submission_callbacks) From d7e6fe477524da9993a5dec7e62374996865d0b3 Mon Sep 17 00:00:00 2001 From: Casper da Costa-Luis Date: Wed, 10 Jul 2024 16:20:27 +0100 Subject: [PATCH 3/4] drop filters --- petric.py | 38 +++++++++++++++++++------------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/petric.py b/petric.py index 69d105c..0a0ba6e 100755 --- a/petric.py +++ b/petric.py @@ -95,32 +95,32 @@ def __init__(self, reference_image, whole_object_mask, background_mask, **kwargs self.whole_object_indices = np.where(whole_object_mask == 1) self.background_indices = np.where(background_mask == 1) self.ref_im_arr = reference_image.as_array() + self.norm = self.ref_im_arr[self.background_indices].mean() def __call__(self, algorithm): iteration = algorithm.iteration if iteration % algorithm.update_objective_interval != 0 and iteration != algorithm.max_iteration: return - for filter_name, filter_func in self.filter.items(): - if filter_func is None: - filter_func = lambda x: x - test_im, ref_im = (filter_func(img_data).as_array() for img_data in (algorithm.x, self.reference_image)) - - # (1) global metrics & statistics - norm = ref_im[self.background_indices].mean() - self.tb_summary_writer.add_scalar( - f"RMSE_whole_object{filter_name}", - np.sqrt(mse(ref_im[self.whole_object_indices], test_im[self.whole_object_indices])) / norm, iteration) + assert not any(self.filter.values()), "Filtering not implemented" + test_im_arr = algorithm.x.as_array() + + # (1) global metrics & statistics + self.tb_summary_writer.add_scalar( + "RMSE_whole_object", + np.sqrt(mse(self.ref_im_arr[self.whole_object_indices], test_im_arr[self.whole_object_indices])) / + self.norm, iteration) + self.tb_summary_writer.add_scalar( + "RMSE_background", + np.sqrt(mse(self.ref_im_arr[self.background_indices], test_im_arr[self.background_indices])) / self.norm, + iteration) + + # (2) local metrics & statistics + for voi_name, voi_indices in sorted(self.voi_indices.items()): + # AEM not to be confused with MAE self.tb_summary_writer.add_scalar( - f"RMSE_background{filter_name}", - np.sqrt(mse(ref_im[self.background_indices], test_im[self.background_indices])) / norm, iteration) - - # (2) local metrics & statistics - for voi_name, voi_indices in self.voi_indices.items(): - # AEM not to be confused with MAE - self.tb_summary_writer.add_scalar( - f"AEM_VOI_{voi_name}{filter_name}", - np.abs(test_im[voi_indices].mean() - ref_im[voi_indices].mean()) / norm, iteration) + f"AEM_VOI_{voi_name}", + np.abs(test_im_arr[voi_indices].mean() - self.ref_im_arr[voi_indices].mean()) / self.norm, iteration) class MetricsWithTimeout(cbks.Callback): From a5860b581d3b40f0ef41ad76ac1b409a2ebb4420 Mon Sep 17 00:00:00 2001 From: Casper da Costa-Luis Date: Wed, 10 Jul 2024 16:50:28 +0100 Subject: [PATCH 4/4] slight tidy --- petric.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/petric.py b/petric.py index 0a0ba6e..ec62dc7 100755 --- a/petric.py +++ b/petric.py @@ -192,16 +192,17 @@ def get_data(srcdir=".", outdir=OUTDIR, sirf_verbosity=0): penalty_strength = 1 / 700 # default choice prior = construct_RDP(penalty_strength, OSEM_image, kappa) - reference_image = STIR.ImageData(str(srcdir / 'reference_image.hv')) if (srcdir / - 'reference_image.hv').is_file() else None - whole_object_mask = STIR.ImageData(str(srcdir / - 'VOI_whole_object.hv')) if (srcdir / - 'VOI_whole_object.hv').is_file() else None - background_mask = STIR.ImageData(str(srcdir / 'VOI_background.hv')) if (srcdir / - 'VOI_background.hv').is_file() else None + def get_image(fname): + if (source := srcdir / 'PETRIC' / fname).is_file(): + return STIR.ImageData(str(source)) + return None # explicit to suppress linter warnings + + reference_image = get_image('reference_image.hv') + whole_object_mask = get_image('VOI_whole_object.hv') + background_mask = get_image('VOI_background.hv') voi_masks = { voi.stem: STIR.ImageData(str(voi)) - for voi in srcdir.glob("VOI_*.hv") if voi.stem[4:] not in ('background', 'whole_object')} + for voi in (srcdir / 'PETRIC').glob("VOI_*.hv") if voi.stem[4:] not in ('background', 'whole_object')} return Dataset(acquired_data, additive_term, mult_factors, OSEM_image, prior, kappa, reference_image, whole_object_mask, background_mask, voi_masks)