Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Enabling PET reconstruction from list mode data #1103

Merged
merged 52 commits into from
Apr 8, 2024
Merged
Show file tree
Hide file tree
Changes from 49 commits
Commits
Show all changes
52 commits
Select commit Hold shift + click to select a range
51fed6c
implemented basic-functionality listmode data class in C++ and Python
evgueni-ovtchinnikov May 10, 2022
c4be9dd
added objective function type for lismode reconstruction
evgueni-ovtchinnikov May 11, 2022
0063f78
introduced abstract base class PETScanData for all raw data objects
evgueni-ovtchinnikov May 12, 2022
b143cfc
[ci skip] interfaced PETScanData into Python
evgueni-ovtchinnikov May 12, 2022
b209f33
Minimal backbone for LM recons
NikEfth May 17, 2022
31ef16c
New test for LM data (test6)
NikEfth May 19, 2022
b38b35b
fix
NikEfth May 19, 2022
036c7c2
Merge branch 'HEAD' into lm-recon
NikEfth May 19, 2022
65f9fe0
Working test6
NikEfth May 19, 2022
97f76cd
Generalize the paths in test6
NikEfth May 19, 2022
26d36fc
Minor renaming
NikEfth May 19, 2022
db4479d
simplify link library to Boost
paskino Feb 1, 2022
debcd5b
Merge branch 'master' into HEAD
paskino May 21, 2022
b59524c
remove dangling endif
paskino May 21, 2022
f65e1de
Merge pull request #1 from paskino/edo_lm-recon
NikEfth May 21, 2022
e4632db
[ci skip] moved data folder TBPET out of SIRF repo
evgueni-ovtchinnikov May 23, 2022
ea18b8d
removed unused stuff from cstir_test6.cpp
evgueni-ovtchinnikov May 23, 2022
7a33289
restored backward compatibility of run_test6.cpp
evgueni-ovtchinnikov May 23, 2022
7e0e299
implemented C interface for list mode reconstruction objective function
evgueni-ovtchinnikov May 23, 2022
47eb909
started interfacing reconstruction from listmode into Python
evgueni-ovtchinnikov May 23, 2022
08f0425
implemented minimal Python interface for reconstruction from listmode…
evgueni-ovtchinnikov May 24, 2022
4dd2de7
resolved some issues raised by Kris
evgueni-ovtchinnikov May 25, 2022
50ef6d4
small amendments in STIR.py
evgueni-ovtchinnikov May 25, 2022
eb2ae63
removed cstir_test6 from ctest tests
evgueni-ovtchinnikov May 25, 2022
27e5d15
minor fixes to listmode classes
KrisThielemans May 26, 2022
ef4c9eb
removed some unneeded data_sptr() methods
evgueni-ovtchinnikov May 26, 2022
608f39a
replaced cache_path quick fix with a more proper handling
evgueni-ovtchinnikov May 26, 2022
8b96919
Merge remote-tracking branch 'origin/master' into lm-recon
KrisThielemans May 27, 2022
05083e8
[ci skip] corrected copyrights in cSTIR/tests/test6.cpp
evgueni-ovtchinnikov May 30, 2022
ed30e6b
make listmode recon compatible with STIR 5.1.0
KrisThielemans Jan 15, 2023
4fa6220
Merge branch 'master' into lm-recon
KrisThielemans Mar 2, 2023
cd0c423
rename PETScanData to ScanData(Python) or STIRScanData (C++)
KrisThielemans Mar 2, 2023
ab26634
Merge remote-tracking branch 'origin/master' into lm-recon
KrisThielemans Jan 15, 2024
dd4e33e
No longer derive ListmodeData from DataContainer
KrisThielemans Jan 17, 2024
a66667b
update listmode test (WIP)
KrisThielemans Jan 17, 2024
5698103
Derive DataContainer from ContainerBase
KrisThielemans Jan 18, 2024
61e692d
speed-up listmode recon test and use SIRF example data
KrisThielemans Jan 18, 2024
27eec3d
attended to Codacy issues
evgueni-ovtchinnikov Jan 19, 2024
81f1389
attended to further Codacy issues
evgueni-ovtchinnikov Jan 19, 2024
6ed55d8
fix ListmodeData hierarchy and add get_info to STIR containers
KrisThielemans Jan 21, 2024
b98c6d0
renamed ListmodeData to STIRListmodeData for consistency
KrisThielemans Jan 21, 2024
713ad10
add set_acquisition_model to listmode obj-fun
KrisThielemans Jan 21, 2024
835b499
minor doc corrections
KrisThielemans Mar 28, 2024
9ce3b51
listmode recon improvements
KrisThielemans Mar 28, 2024
f009fbf
Merge remote-tracking branch 'origin/master' into lm-recon
KrisThielemans Mar 28, 2024
bb42564
minor fix in demo
KrisThielemans Mar 28, 2024
56fbace
fix storing filename of listmode file
KrisThielemans Mar 28, 2024
b384409
add SIRF_DATA_PATH to STIR C++ test
KrisThielemans Apr 3, 2024
95b3a24
fixed bugs and issues found in lmrecontest2
evgueni-ovtchinnikov Apr 4, 2024
bf361ff
[ci skip] got rid of gpu stuff in reconstruct_from_listmode.py
evgueni-ovtchinnikov Apr 5, 2024
752c077
[ci skip] removed osem_lm_reconstruction.py, superseded by listmode_r…
evgueni-ovtchinnikov Apr 5, 2024
e4ac001
[ci skip] updated CHANGES.md
evgueni-ovtchinnikov Apr 8, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 11 additions & 0 deletions doc/UserGuide.md
Original file line number Diff line number Diff line change
Expand Up @@ -245,6 +245,16 @@ In what follows, we mark by `PET` classes defined in `sirf.STIR` only and by `MR

We remind that every derived class inherits all methods of its base class.

##### ListmodeData (PET)

An STIR-specitic acquisition data container class for list-mode data objects. Inherits from `sirf.SIRF.ContainerBase`.

###### Methods:

ListmodeData Constructor. If no arguments are present, creates an
empty object, otherwise specifies the file containing raw data
get_info (PET) Returns information on the acquisition data as a string.

##### AcquisitionData

An engine-specific acquisition data container class for acquisition data objects. Inherits from `sirf.SIRF.DataContainer`.
Expand Down Expand Up @@ -288,6 +298,7 @@ An engine-specific image data container class for data representing 3D objects.
fill Replaces the object data with user-supplied data.
as_array Returns the object data as an array.
read_from_file Reads the image data from file.
get_info (PET) Returns information on the data as a string.
get_ISMRMRD_info
(MR) Returns information on the image data as an array.
get_uniform_copy
Expand Down
267 changes: 267 additions & 0 deletions examples/Python/PET/listmode_reconstruction.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,267 @@
'''Listmode reconstruction demo

Usage:
listmode_reconstruction [--help | options]

Options:
-p <path>, --path=<path> path to data files, defaults to data/examples/PET/mMR
subfolder of SIRF root folder
-l <list>, --list=<list> listmode file [default: list.l.hdr]
-g <sino>, --sino=<sino> output file prefix [default: sinograms]
-a <attn>, --attn=<attn> attenuation image file file [default: mu_map.hv]
-n <norm>, --norm=<norm> ECAT8 bin normalization file [default: norm.n.hdr]
-I <int>, --interval=<int> scanning time interval to convert as string '(a,b)'
(no space after comma) [default: (0,50)]
-d <nxny>, --nxny=<nxny> image x and y dimensions as string '(nx,ny)'
(no space after comma) [default: (127,127)]
-S <subs>, --subs=<subs> number of subsets [default: 7]
-i <siter>, --subiter=<siter> number of sub-iterations [default: 2]
-o <outp>, --outp=<outp> output file prefix [default: recon]
-e <engn>, --engine=<engn> reconstruction engine [default: STIR]
-s <stsc>, --storage=<stsc> acquisition data storage scheme [default: file]
-C <cnts>, --counts=<cnts> account for delay between injection and acquisition start by shifting interval to start when counts exceed given threshold.
--visualisations show visualisations
--nifti save output as nifti
--gpu use gpu
--non-interactive do not show plots
'''

## SyneRBI Synergistic Image Reconstruction Framework (SIRF)
## Copyright 2018 - 2020 Rutherford Appleton Laboratory STFC
## Copyright 2018 - 2021, 2024 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.

__version__ = '1.1.0'
from docopt import docopt
args = docopt(__doc__, version=__version__)

from ast import literal_eval
import os

import PET_plot_functions

from sirf.Utilities import error, examples_data_path, existing_filepath
from sirf.Utilities import show_2D_array

# import engine module
import importlib
engine = args['--engine']
pet = importlib.import_module('sirf.' + engine)


# process command-line options
data_path = args['--path']
if data_path is None:
# default to data/examples/PET/mMR
# Note: seem to need / even on Windows
#data_path = os.path.join(examples_data_path('PET'), 'mMR')
data_path = examples_data_path('PET') + '/mMR'
print('Finding files in %s' % data_path)

list_file = args['--list']
sino_file = args['--sino']
norm_file = args['--norm']
attn_file = args['--attn']
outp_file = args['--outp']
# Check file exists (e.g., absolute path). Else prepend data_path
if not os.path.isfile(list_file):
list_file = existing_filepath(data_path, list_file)
if not os.path.isfile(norm_file):
norm_file = existing_filepath(data_path, norm_file)
if not os.path.isfile(attn_file):
attn_file = existing_filepath(data_path, attn_file)
nxny = literal_eval(args['--nxny'])
input_interval = literal_eval(args['--interval'])
num_subsets = int(args['--subs'])
num_subiterations = int(args['--subiter'])
storage = args['--storage']
count_threshold = args['--counts']

if args['--visualisations']:
visualisations = True
else:
visualisations = False
if args['--non-interactive']:
visualisations = False

if args['--gpu']:
use_gpu = True
# import sirf.Reg
else:
use_gpu = False


def main():

# engine's messages go to files, except error messages, which go to stdout
_ = pet.MessageRedirector('info.txt', 'warn.txt')

# select acquisition data storage scheme
pet.AcquisitionData.set_storage_scheme(storage)

listmode_data = pet.ListmodeData(list_file)

# First step is to create AcquisitionData ("sinograms") from the
# listmode file.
# See the listmode_to_sinograms demo for some more information on this step.

# create listmode-to-sinograms converter object
# See also the listmode_to_sinograms demo
lm2sino = pet.ListmodeToSinograms()

# set input, output and template files
lm2sino.set_input(listmode_data)
lm2sino.set_output_prefix(sino_file)
# need to be at maximum resolution to work in listmode reconstruction
acq_data_template = listmode_data.acquisition_data_template()
#print(acq_data_template.get_info())
lm2sino.set_template(acq_data_template)

if count_threshold is None:
interval = input_interval
else:
time_shift = lm2sino.get_time_at_which_num_prompts_exceeds_threshold(count_threshold)
if time_shift < 0:
print("No time found at which count rate exceeds " + str(time_shift) + ", not modifying interval")
interval = (input_interval[0]+time_shift, input_interval[1]+time_shift)
print("Time at which count rate exceeds " + str(count_threshold) + " = " + str(time_shift) + " s.")
print("Input intervals: " + str(input_interval[0]) + ", " + str(input_interval[1]))
print("Modified intervals: " + str(interval[0]) + ", " + str(interval[1]))

# set interval
lm2sino.set_time_interval(interval[0], interval[1])

# set up the converter
lm2sino.set_up()

# convert (need it for the scatter estimate)
lm2sino.process()

acq_data = lm2sino.get_output()

# Get the randoms
randoms = lm2sino.estimate_randoms()

# create initial image estimate of dimensions and voxel sizes
# compatible with the scanner geometry (included in the AcquisitionData
# object acq_data) and initialize each voxel to 1.0
image = acq_data.create_uniform_image(1.0, nxny)

# read attenuation image
attn_image = pet.ImageData(attn_file)
if visualisations:
z = attn_image.shape[0]//2
attn_image_as_array = attn_image.as_array()
show_2D_array('Attenuation image', attn_image_as_array[z,:,:])

# select acquisition model that implements the geometric
# forward projection by a ray tracing matrix multiplication
acq_model = pet.AcquisitionModelUsingRayTracingMatrix()
acq_model.set_num_tangential_LORs(10)

# create acquisition sensitivity model from normalisation data
asm_norm = pet.AcquisitionSensitivityModel(norm_file)

# create attenuation factors
asm_attn = pet.AcquisitionSensitivityModel(attn_image, acq_model)
# converting attenuation image into attenuation factors (one for every bin)
asm_attn.set_up(acq_data)
ac_factors = acq_data.get_uniform_copy(value=1)
print('applying attenuation (please wait, may take a while)...')
asm_attn.unnormalise(ac_factors)
asm_attn = pet.AcquisitionSensitivityModel(ac_factors)

# scatter estimation
print('estimating scatter (this will take a while!)')
scatter_estimator = pet.ScatterEstimator()

scatter_estimator.set_input(acq_data)
scatter_estimator.set_attenuation_image(attn_image)
scatter_estimator.set_randoms(randoms)
scatter_estimator.set_asm(asm_norm)
# invert attenuation factors to get the correction factors,
# as this is unfortunately what a ScatterEstimator needs
acf_factors=acq_data.get_uniform_copy()
acf_factors.fill(1/ac_factors.as_array())
scatter_estimator.set_attenuation_correction_factors(acf_factors)
scatter_estimator.set_output_prefix(sino_file + '_scatter')
scatter_estimator.set_num_iterations(3)
scatter_estimator.set_up()
scatter_estimator.process()
scatter_estimate = scatter_estimator.get_output()
if visualisations:
scatter_estimate_as_array = scatter_estimate.as_array()
show_2D_array('Scatter estimate (first sinogram)', scatter_estimate_as_array[0, 0, :, :])
PET_plot_functions.plot_sinogram_profile(acq_data, randoms=randoms, scatter=scatter_estimate)

# chain attenuation and ECAT8 normalisation
asm = pet.AcquisitionSensitivityModel(asm_norm, asm_attn)
asm.set_up(acq_data)

acq_model.set_acquisition_sensitivity(asm)
acq_model.set_background_term(randoms + scatter_estimate)

# define objective function to be maximized as
# Poisson logarithmic likelihood (with linear model for mean)
obj_fun = pet.PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin();
obj_fun.set_acquisition_model(acq_model)
obj_fun.set_acquisition_data(listmode_data)

# select Ordered Subsets Maximum A-Posteriori One Step Late as the
# reconstruction algorithm (since we are not using a penalty, or prior, in
# this example, we actually run OSEM);
# this algorithm does not converge to the maximum of the objective function
# but is used in practice to speed-up calculations
# See the reconstruction demos for more complicated examples
recon = pet.OSMAPOSLReconstructor()
recon.set_objective_function(obj_fun)
recon.set_num_subsets(num_subsets)
recon.set_num_subiterations(num_subiterations)

# set up the reconstructor based on a sample image
# (checks the validity of parameters, sets up objective function
# and other objects involved in the reconstruction, which involves
# computing/reading sensitivity image etc etc.)
print('setting up, please wait...')
recon.set_up(image)

# set the initial image estimate
recon.set_current_estimate(image)

# reconstruct
print('reconstructing, please wait...')
recon.process()
out = recon.get_output()
if not args['--nifti']:
out.write(outp_file)
else:
import sirf.Reg as reg
reg.NiftiImageData(out).write(outp_file)

if visualisations:
# show reconstructed image
z = out.shape[0]//2
image_array = out.as_array()
show_2D_array('Reconstructed image', image_array[z,:,:])
# pylab.show()


try:
main()
print('\n=== done with %s' % __file__)

except error as err:
print('%s' % err.value)
86 changes: 86 additions & 0 deletions examples/Python/PET/osem_lm_reconstruction.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,86 @@
'''OSEM reconstruction from listmode demo.


Usage:
osem_lm_reconstruction [--help | options] <data_path>

Options:
-e <engn>, --engine=<engn> reconstruction engine [default: STIR]
'''


__version__ = '0.1.0'
from docopt import docopt
args = docopt(__doc__, version=__version__)


from sirf.Utilities import error
#, examples_data_path, existing_filepath
#from sirf.Utilities import show_2D_array

# import engine module
import importlib
engine = args['--engine']
pet = importlib.import_module('sirf.' + engine)


data_path = args['<data_path>'] + '/'


def main():

# direct all engine's information and warnings printing to files
_ = pet.MessageRedirector('info.txt', 'warn.txt')

cache_path = data_path;
sens_filename = cache_path + "sens_0.hv";
#print(sens_filename)
tmpl_projdata_filename = data_path + "tmpl_scanner.hs";
#print(tmpl_projdata_filename)

acq_tmpl = pet.AcquisitionData(tmpl_projdata_filename)
#print(acq_tmpl.dimensions())

img_data = pet.ImageData(sens_filename)
#print(img_data.dimensions())

#obj_fun = PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin()
obj_fun = pet.make_Poisson_loglikelihood( \
likelihood_type='LinearModelForMeanAndListModeDataWithProjMatrixByBin')
obj_fun.set_cache_path(cache_path, True)
obj_fun.set_skip_lm_input_file(True);
obj_fun.set_skip_balanced_subsets(True);
obj_fun.set_acquisition_data(acq_tmpl);
obj_fun.set_skip_balanced_subsets(True);
obj_fun.set_max_ring_difference(60);
obj_fun.set_cache_max_size(1500000000);

#print(obj_fun.get_cache_path())
#print(obj_fun.get_cache_max_size())
#print(obj_fun.get_subsensitivity_filenames())

num_subiterations = 11
recon = pet.OSMAPOSLReconstructor()
recon.set_objective_function(obj_fun)
recon.set_num_subiterations(num_subiterations)
recon.set_output_filename_prefix(data_path + "/my_output");
recon.set_save_interval(1);

print('setting up, please wait...')
recon.set_up(img_data)

print('reconstructing, please wait...')
recon.reconstruct(img_data)

img_out = recon.get_output()
print(img_out.dimensions())

# if anything goes wrong, an exception will be thrown
# (cf. Error Handling section in the spec)
try:
main()
print('\n=== done with %s' % __file__)

except error as err:
# display error information
print('%s' % err.value)
Loading
Loading