From 3db925bfa0014b4af1c682eaf157754d4ca287cf Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 2 Oct 2024 13:17:27 +0100 Subject: [PATCH 1/3] optionally disable MessageRedirector --- petric.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/petric.py b/petric.py index 31e2d20..a53ffc0 100755 --- a/petric.py +++ b/petric.py @@ -212,14 +212,15 @@ class Dataset: def get_data(srcdir=".", outdir=OUTDIR, sirf_verbosity=0): """ Load data from `srcdir`, constructs prior and return as a `Dataset`. - Also redirects sirf.STIR log output to `outdir`. + Also redirects sirf.STIR log output to `outdir`, unless that's set to None """ srcdir = Path(srcdir) - outdir = Path(outdir) STIR.set_verbosity(sirf_verbosity) # set to higher value to diagnose problems STIR.AcquisitionData.set_storage_scheme('memory') # needed for get_subsets() - _ = STIR.MessageRedirector(str(outdir / 'info.txt'), str(outdir / 'warnings.txt'), str(outdir / 'errors.txt')) + if outdir is not None: + outdir = Path(outdir) + _ = STIR.MessageRedirector(str(outdir / 'info.txt'), str(outdir / 'warnings.txt'), str(outdir / 'errors.txt')) acquired_data = STIR.AcquisitionData(str(srcdir / 'prompts.hs')) additive_term = STIR.AcquisitionData(str(srcdir / 'additive_term.hs')) mult_factors = STIR.AcquisitionData(str(srcdir / 'mult_factors.hs')) From a530df225053d9eda97b0cb663100dfd6b66935e Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 2 Oct 2024 13:19:29 +0100 Subject: [PATCH 2/3] add get_penalisation_factor --- .../get_penalisation_factor.py | 88 +++++++++++++++++++ 1 file changed, 88 insertions(+) create mode 100755 SIRF_data_preparation/get_penalisation_factor.py diff --git a/SIRF_data_preparation/get_penalisation_factor.py b/SIRF_data_preparation/get_penalisation_factor.py new file mode 100755 index 0000000..c2ccaee --- /dev/null +++ b/SIRF_data_preparation/get_penalisation_factor.py @@ -0,0 +1,88 @@ +#!/usr/bin/env python +"""Find penalisation factor for one dataset based on another + +Usage: + get_penalisation_factor.py [--help | options] + +Options: + -h, --help + --dataset= dataset name (required) + --ref_dataset= reference dataset name (required) + -w, --write_penalisation_factor write in data//penalisation_factor.txt +""" + +# Copyright 2024 University College London +# Licence: Apache-2.0 + +# %% imports +import sirf.STIR +import numpy as np +import sys +from docopt import docopt +from pathlib import Path +from petric import get_data, Dataset +from SIRF_data_preparation.data_utilities import the_data_path + +# %% +__version__ = "0.1.0" + +write_penalisation_factor = False +if ( + not "ipykernel" in sys.argv[0] +): # clunky way to be able to set variables from within jupyter/VScode without docopt + args = docopt(__doc__, argv=None, version=__version__) + # logging.basicConfig(level=logging.INFO) + dataset = args["--dataset"] + ref_dataset = args["--ref_dataset"] + if dataset is None or ref_dataset is None: + print("Need to set the --dataset arguments") + exit(1) + if args["--write_penalisation_factor"] is not None: + write_penalisation_factor = False +else: + # set it by hand, e.g. + ref_dataset = "NeuroLF_Hoffman_Dataset" + dataset = "Siemens_mMR_NEMA_IQ" + + write_penalisation_factor = True + + +# %% +def VOImean(im: sirf.STIR.ImageData, background_mask: sirf.STIR.ImageData) -> float: + background_indices = np.where(background_mask.as_array()) + return np.mean(im.as_array()[background_indices]) / len(background_indices) + + +def backgroundVOImean(dataset: Path) -> float: + im = sirf.STIR.ImageData(str(dataset / "OSEM_image.hv")) + VOI = sirf.STIR.ImageData(str(dataset / "PETRIC" / "VOI_background.hv")) + return VOImean(im, VOI) + + +def get_penalisation_factor(ref_data: Dataset, cur_data: Dataset) -> float: + ref_mean = VOImean(ref_data.OSEM_image, ref_data.background_mask) + cur_mean = VOImean(cur_data.OSEM_image, cur_data.background_mask) + print( + f"ref_mean={ref_mean}, cur_mean={cur_mean}, c/r={cur_mean / ref_mean}, r/c={ref_mean / cur_mean}" + ) + ref_penalisation_factor = ref_data.prior.get_penalisation_factor() + penalisation_factor = ref_penalisation_factor * cur_mean / ref_mean + print(f"penalisation_factor={penalisation_factor}") + return penalisation_factor + + +# %% +refdir = Path(the_data_path(ref_dataset)) +curdir = Path(the_data_path(dataset)) +# %% +ref_data = get_data(refdir, outdir=None) +cur_data = get_data(curdir, outdir=None) +# %% +penalisation_factor = get_penalisation_factor(ref_data, cur_data) + +# %% +if write_penalisation_factor: + filename = curdir / "penalisation_factor.txt" + print(f"Writing it to {filename}") + with open(filename, "w") as file: + file.write(str(penalisation_factor)) From 52272ff9847219bb30b215c16c2fb51275006f37 Mon Sep 17 00:00:00 2001 From: Casper da Costa-Luis Date: Wed, 2 Oct 2024 14:26:04 +0100 Subject: [PATCH 3/3] lint --- .../get_penalisation_factor.py | 26 +++++++++---------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/SIRF_data_preparation/get_penalisation_factor.py b/SIRF_data_preparation/get_penalisation_factor.py index c2ccaee..a6ce158 100755 --- a/SIRF_data_preparation/get_penalisation_factor.py +++ b/SIRF_data_preparation/get_penalisation_factor.py @@ -14,24 +14,26 @@ # Copyright 2024 University College London # Licence: Apache-2.0 -# %% imports -import sirf.STIR -import numpy as np import sys -from docopt import docopt from pathlib import Path -from petric import get_data, Dataset + +import numpy as np +from docopt import docopt + +# %% imports +import sirf.STIR +from petric import Dataset, get_data from SIRF_data_preparation.data_utilities import the_data_path # %% __version__ = "0.1.0" write_penalisation_factor = False -if ( - not "ipykernel" in sys.argv[0] -): # clunky way to be able to set variables from within jupyter/VScode without docopt +if "ipykernel" not in sys.argv[0]: # clunky way to be able to set variables from within jupyter/VScode without docopt args = docopt(__doc__, argv=None, version=__version__) + # logging.basicConfig(level=logging.INFO) + dataset = args["--dataset"] ref_dataset = args["--ref_dataset"] if dataset is None or ref_dataset is None: @@ -39,11 +41,9 @@ exit(1) if args["--write_penalisation_factor"] is not None: write_penalisation_factor = False -else: - # set it by hand, e.g. +else: # set it by hand, e.g. ref_dataset = "NeuroLF_Hoffman_Dataset" dataset = "Siemens_mMR_NEMA_IQ" - write_penalisation_factor = True @@ -62,9 +62,7 @@ def backgroundVOImean(dataset: Path) -> float: def get_penalisation_factor(ref_data: Dataset, cur_data: Dataset) -> float: ref_mean = VOImean(ref_data.OSEM_image, ref_data.background_mask) cur_mean = VOImean(cur_data.OSEM_image, cur_data.background_mask) - print( - f"ref_mean={ref_mean}, cur_mean={cur_mean}, c/r={cur_mean / ref_mean}, r/c={ref_mean / cur_mean}" - ) + print(f"ref_mean={ref_mean}, cur_mean={cur_mean}, c/r={cur_mean / ref_mean}, r/c={ref_mean / cur_mean}") ref_penalisation_factor = ref_data.prior.get_penalisation_factor() penalisation_factor = ref_penalisation_factor * cur_mean / ref_mean print(f"penalisation_factor={penalisation_factor}")