From ec7a95eb8901dbe60c44b7e6b7b779f457ec2372 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Tue, 18 Jun 2024 17:30:46 +0100 Subject: [PATCH] Nema data preparation (#3) create scripts for converting Siemens mMR data in the format required by PETRIC --------- Co-authored-by: Evgueni Ovtchinnikov Co-authored-by: Edoardo Pasca <14138589+paskino@users.noreply.github.com> --- .gitignore | 9 + lib/data_utilities/__init__.py | 183 ++++++++++++++++++ lib/data_utilities/data_path.py | 1 + .../PET_plot_functions.py | 53 +++++ .../prepare_mMR_NEMA_IQ_data.py | 57 ++++++ 5 files changed, 303 insertions(+) create mode 100644 .gitignore create mode 100644 lib/data_utilities/__init__.py create mode 100644 lib/data_utilities/data_path.py create mode 100644 src/SIRF_data_preparation/PET_plot_functions.py create mode 100644 src/SIRF_data_preparation/prepare_mMR_NEMA_IQ_data.py diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..ad4f10a --- /dev/null +++ b/.gitignore @@ -0,0 +1,9 @@ +data/ +tmp*/ +errr.txt +info.txt +warn.txt +__pycache__ +*.s +*.hs + diff --git a/lib/data_utilities/__init__.py b/lib/data_utilities/__init__.py new file mode 100644 index 0000000..6366d5c --- /dev/null +++ b/lib/data_utilities/__init__.py @@ -0,0 +1,183 @@ +'''Library of Siemens data preparation utilities.''' + +# Author: Ashley Gillman, Kris Thielemans, Evgueni Ovtchinnikov +# Copyright (C) 2021 Commonwealth Scientific and Industrial Research Organisation +# Copyright (C) 2024 University College London +# Copyright (C) 2024 STFC, UK Research and Innovation +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# http://www.apache.org/licenses/LICENSE-2.0 +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import os +import importlib +pet = importlib.import_module('sirf.STIR') +import logging + +logger = logging.getLogger("PETRIC") + + +def the_data_path(*data_type): + ''' + Returns the path to data. + + data_type: either 'PET', 'MR' or 'Synergistic', or use multiple arguments for + subdirectories like the_data_path('PET', 'mMR', 'NEMA_IQ'). + ''' + try: + from .data_path import data_path + except ImportError: + raise RuntimeError( + "Path to data not found.") + + return os.path.join(data_path, *data_type) + + +def fix_siemens_norm_EOL(in_filename, out_filename): + with open(in_filename, mode="rb") as f: + data = bytearray(f.read()) + for i in range(len(data)): + if data[i] == 13: #'\r' + data[i] = 10 # \n + with open(out_filename, mode="wb") as f: + f.write(data) + + +def prepare_challenge_Siemens_data(data_path, challenge_data_path, intermediate_data_path, + f_root, f_listmode, f_mumap, f_attn, f_norm, f_stir_norm, f_template, + f_prompts, f_multfactors, f_additive, f_randoms, f_af, f_acf, f_scatter, start, stop): + '''Prepares data for SyneRBI Challenge24 + + data_path: path to Siemens data + challenge_data_path: path to prepared data + data_path: path to intermediate data + f_root: common prefix for some data files' names (list-mode data, mu-maps etc.) + f_listmode: list-mode data file suffix + f_numap: mu-map file suffix + f_attn: mu-map header suffix + f_norm: Siemens normalisation data file suffix + f_stir_norm: STIR normalisation data file name + f_template: template for prompts file name + f_prompts: prompts file name + f_multfactors: multfactors file name + f_additive: additive term file name + f_randoms: estimated randoms file name + f_af: attenuation factor file name + f_acf: attenuation correction factor file name + f_scatter: scatter estimate file name + start: start time for data acquisition + stop: end time for data acquisition + ''' + + logging.info(f"Start time for data: {start} sec") + logging.info(f"End time for data: {stop} sec") + + f_listmode = os.path.join(data_path, f_root + f_listmode) + f_siemens_attn_image = os.path.join(data_path, f_root + f_mumap) + f_siemens_attn_header = f_siemens_attn_image + '.hdr' + f_siemens_norm = os.path.join(data_path, f_root + f_norm) + f_siemens_norm_header = f_siemens_norm + '.hdr' + f_stir_norm_header = os.path.join(challenge_data_path, f_stir_norm) + f_stir_attn_header = os.path.join(challenge_data_path, f_root + f_attn) + f_prompts = os.path.join(challenge_data_path, f_prompts) + f_multfactors = os.path.join(challenge_data_path, f_multfactors) + f_additive = os.path.join(challenge_data_path, f_additive) + f_randoms = os.path.join(intermediate_data_path, f_randoms) + f_af = os.path.join(intermediate_data_path, f_af) + f_acf = os.path.join(intermediate_data_path, f_acf) + f_scatter = os.path.join(intermediate_data_path, f_scatter) + f_info = os.path.join(intermediate_data_path, 'info.txt') + f_warn = os.path.join(intermediate_data_path, 'warn.txt') + + os.system('cp ' + f_siemens_attn_image + ' ' + challenge_data_path) + os.system('convertSiemensInterfileToSTIR.sh ' + f_siemens_attn_header + ' ' + f_stir_attn_header) + os.system('cp ' + f_siemens_norm + ' ' + challenge_data_path) + + fix_siemens_norm_EOL(f_siemens_norm_header, f_stir_norm_header) + + # engine's messages go to files, except error messages, which go to stdout + _ = pet.MessageRedirector(f_info, f_warn) + + # select acquisition data storage scheme + pet.AcquisitionData.set_storage_scheme('memory') + + # read acquisition data template + acq_data_template = pet.AcquisitionData(f_template) + logger.info(acq_data_template.norm()) + + output_prefix = os.path.join(challenge_data_path, "prompts") + + listmode_data = pet.ListmodeData(f_listmode) + + logger.info("Processing listmode data...") + logger.info(f"Start time: {start} sec. End time: {stop} sec.") + logger.info(f"Output prefix: {output_prefix}") + # create listmode-to-sinograms converter object + lm2sino = pet.ListmodeToSinograms() + lm2sino.set_input(listmode_data) + lm2sino.set_output_prefix(output_prefix) + lm2sino.set_template(acq_data_template) + lm2sino.set_time_interval(start, stop) + lm2sino.set_up() + lm2sino.process() + + prompts = lm2sino.get_output() + randoms = lm2sino.estimate_randoms() + logger.info('data shape: %s' % repr(prompts.shape)) + logger.info('prompts norm: %f' % prompts.norm()) + logger.info('randoms norm: %f' % randoms.norm()) + + logger.info(f'writing prompts to {f_prompts} and randoms to {f_randoms}') + prompts.write(f_prompts) + randoms.write(f_randoms) + + attn_image = pet.ImageData(f_stir_attn_header) + af, acf = pet.AcquisitionSensitivityModel.compute_attenuation_factors(prompts, attn_image) + logger.info('norm of the attenuation factor: %f' % af.norm()) + logger.info('norm of the attenuation correction factor: %f' % acf.norm()) + logger.info(f'writing intermediate attenuation factors to {f_af}') + logger.info(f'writing intermediate attenuation coefficient factors to {f_acf}') + af.write(f_af) + acf.write(f_acf) + + asm = pet.AcquisitionSensitivityModel(f_stir_norm_header) + se = pet.ScatterEstimator() + se.set_input(prompts) + se.set_attenuation_image(attn_image) + se.set_randoms(randoms) + se.set_asm(asm) + se.set_attenuation_correction_factors(acf) + se.set_num_iterations(4) + se.set_OSEM_num_subsets(7) + se.set_output_prefix(f_scatter) + se.set_up() + se.process() + scatter = se.get_output() + logger.info('norm of the scatter estimate: %f' % scatter.norm()) + + multfact = af.clone() + asm.set_up(af) + asm.unnormalise(multfact) + logger.info(multfact.norm()) + logger.info(f'writing multiplicative factors to {f_multfactors}') + multfact.write(f_multfactors) + + background = randoms + scatter + logger.info('norm of the background term: %f' % background.norm()) + + asm_mf = pet.AcquisitionSensitivityModel(multfact) + asm_mf.set_up(background) + asm_mf.normalise(background) + logger.info('norm of the additive term: %f' % background.norm()) + logger.info(f'writing additive term to {f_additive}') + + background.write(f_additive) + + return + diff --git a/lib/data_utilities/data_path.py b/lib/data_utilities/data_path.py new file mode 100644 index 0000000..7dad7bd --- /dev/null +++ b/lib/data_utilities/data_path.py @@ -0,0 +1 @@ +data_path = '/home/sirfuser/devel/SyneRBI-Challenge/data' diff --git a/src/SIRF_data_preparation/PET_plot_functions.py b/src/SIRF_data_preparation/PET_plot_functions.py new file mode 100644 index 0000000..0b6730b --- /dev/null +++ b/src/SIRF_data_preparation/PET_plot_functions.py @@ -0,0 +1,53 @@ +''' +Some functions for plotting used in the examples +''' + +## CCP SyneRBI Synergistic Image Reconstruction Framework (SIRF) +## Copyright 2020-2021 University College London +## +## This is software developed for the Collaborative Computational +## Project in Synergistic Reconstruction for Biomedical Imaging (formerly CCP PETMR) +## (http://www.ccpsynerbi.ac.uk/). +## +## Licensed under the Apache License, Version 2.0 (the "License"); +## you may not use this file except in compliance with the License. +## You may obtain a copy of the License at +## http://www.apache.org/licenses/LICENSE-2.0 +## Unless required by applicable law or agreed to in writing, software +## distributed under the License is distributed on an "AS IS" BASIS, +## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +## See the License for the specific language governing permissions and +## limitations under the License. + +import sirf.STIR as PET +import numpy as np +import matplotlib.pyplot as plt + +def plot_sinogram_profile(prompts, randoms=None, scatter=None, sumaxis=(0,1), select=0): + """ + Plot a profile through sirf.STIR.AcquisitionData + + sumaxis: axes to sum over (passed to numpy.sum(..., axis)) + select: element to select after summing + """ + if isinstance(prompts, str): + prompts = PET.AcquisitionData(prompts) + if isinstance(randoms, str): + randoms = PET.AcquisitionData(randoms) + if isinstance(scatter, str): + scatter = PET.AcquisitionData(scatter) + # we will average over all sinograms to reduce noise + plt.figure() + ax = plt.subplot(111) + plt.plot(np.sum(prompts.as_array(), axis=sumaxis)[select,:], label='prompts') + if randoms is None: + if scatter is not None: + plt.plot(np.sum(scatter.as_array(), axis=sumaxis)[select,:], label='scatter') + else: + randoms_as_array = randoms.as_array() + plt.plot(np.sum(randoms_as_array, axis=sumaxis)[select,:], label='randoms') + if scatter is not None: + plt.plot(np.sum(scatter.as_array() + randoms_as_array, axis=sumaxis)[select,:], label='randoms+scatter') + ax.legend() + plt.show() + diff --git a/src/SIRF_data_preparation/prepare_mMR_NEMA_IQ_data.py b/src/SIRF_data_preparation/prepare_mMR_NEMA_IQ_data.py new file mode 100644 index 0000000..542d699 --- /dev/null +++ b/src/SIRF_data_preparation/prepare_mMR_NEMA_IQ_data.py @@ -0,0 +1,57 @@ +import os +import sys +this_directory = os.path.dirname(__file__) +repo_directory = os.path.dirname(os.path.dirname(this_directory)) +sys.path.append(os.path.join(repo_directory, 'lib')) + +import importlib +pet = importlib.import_module('sirf.STIR') + +from sirf.Utilities import examples_data_path +from data_utilities import the_data_path +from data_utilities import prepare_challenge_Siemens_data + +import logging +import argparse + +if __name__ == '__main__': + parser = argparse.ArgumentParser(description='SyneRBI Challenge 2024 data preparation script.') + + parser.add_argument('--log', type=str, default='warning') + parser.add_argument('--start', type=float, default=0) + parser.add_argument('--end', type=float, default=600) + parser.add_argument('--raw_data_path', type=str, default=None) + args = parser.parse_args() + + if args.log in ['debug', 'info', 'warning', 'error', 'critical']: + level = eval(f'logging.{args.log.upper()}') + logging.basicConfig(level=level) + logging.info(f"Setting logging level to {args.log.upper()}") + + start = args.start + end = args.end + + if args.raw_data_path is None: + data_path = the_data_path('PET', 'mMR', 'NEMA_IQ') + else: + data_path = args.raw_data_path + + data_path = os.path.abspath(data_path) + logging.debug(f"Raw data path: {data_path}") + + + sirf_data_path = os.path.join(examples_data_path('PET'), 'mMR') + challenge_data_path = os.path.join(repo_directory, 'data') + intermediate_data_path = os.path.join(challenge_data_path, 'intermediate') + os.makedirs(challenge_data_path, exist_ok=True) + os.chdir(challenge_data_path) + os.makedirs(intermediate_data_path, exist_ok=True) + + f_template = os.path.join(sirf_data_path, 'mMR_template_span11') + os.system('cp ' + f_template + '* ' + challenge_data_path) + + prepare_challenge_Siemens_data(data_path, challenge_data_path, intermediate_data_path, '20170809_NEMA_', + '60min_UCL.l.hdr', 'MUMAP_UCL.v', 'MUMAP_UCL.hv', 'UCL.n', 'norm.n.hdr', f_template + '.hs', + 'prompts', 'multfactors', 'additive', 'randoms', + 'attenuation_factor', 'attenuation_correction_factor', 'scatter', start, end) +