Skip to content

Commit cc0228c

Browse files
lukepolsonlukepolson
authored
Merge pull request #40 from PyTomography/main
Update branch with main
2 parents 3d98246 + 444982a commit cc0228c

File tree

49 files changed

+4488
-4288
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

49 files changed

+4488
-4288
lines changed

‎CMakeLists.txt‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,15 +1,15 @@
11
cmake_minimum_required(VERSION 3.16.3...3.19.7 FATAL_ERROR)
22

3-
project(Pytomography Slicer SPECT Reconstruction)
3+
project(SPECTRecon)
44

55
#-----------------------------------------------------------------------------
66
# Extension meta-information
7-
set(EXTENSION_HOMEPAGE "https://github.com/PyTomography/slicer_spect_recon")
7+
set(EXTENSION_HOMEPAGE "https://github.com/PyTomography/Slicer")
88
set(EXTENSION_CATEGORY "Tomographic Reconstruction")
99
set(EXTENSION_CONTRIBUTORS "Obed Dzikunu (QURIT), Luke Polson (QURIT), Maziar Sabouri (QURIT), Shadab Ahamed (QURIT)")
10-
set(EXTENSION_DESCRIPTION "This module implement GPU accelerated SPECT image reconstruction. User Manual: https://github.com/PyTomography/slicer_spect_recon/blob/main/User_Manual.md")
11-
set(EXTENSION_ICONURL "https://github.com/PyTomography/slicer_spect_recon/blob/main/Pytomography.png")
12-
set(EXTENSION_SCREENSHOTURLS "https://github.com/PyTomography/slicer_spect_recon/blob/main/images/screenshot.png")
10+
set(EXTENSION_DESCRIPTION "This module implement GPU accelerated SPECT image reconstruction. User Manual: https://github.com/PyTomography/Slicer/blob/main/User_Manual.md")
11+
set(EXTENSION_ICONURL "https://raw.githubusercontent.com/PyTomography/Slicer/main/Pytomography.png")
12+
set(EXTENSION_SCREENSHOTURLS "https://raw.githubusercontent.com/PyTomography/Slicer/main/images/screenshot.png")
1313
set(EXTENSION_DEPENDS "NA") # Specified as a list or "NA" if no dependencies
1414

1515
#-----------------------------------------------------------------------------
@@ -19,7 +19,7 @@ include(${Slicer_USE_FILE})
1919

2020
#-----------------------------------------------------------------------------
2121
# Extension modules
22-
add_subdirectory(SlicerSPECTRecon)
22+
add_subdirectory(SPECTRecon)
2323
## NEXT_MODULE
2424

2525
#-----------------------------------------------------------------------------

‎README.md‎

Lines changed: 6 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
# SPECT Tomographic Reconstruction 3D Slicer Extension
22

3-
This is the official repository for the `Slicer` extension `SlicerSPECTRecon`.
3+
This is the official repository for the `Slicer` extension `SlicerSPECTRecon`. **For details on how to use, see the [user manual](https://github.com/PyTomography/Slicer/blob/main/User_Manual.md) and the [associated youtube video tutorial](https://www.youtube.com/watch?v=DzV1soWTzEA)**
44

55
This module enables the reconstruction of raw SPECT projection data, providing customizable options for image modeling and image reconstruction. The module has a SIMIND to DICOM converter to permit reconstruction of SIMIND Monte Carlo data.
66

@@ -15,11 +15,13 @@ The module is divided into the following sections:
1515
Please refer to the `User_Manual.md` file for further information
1616

1717
## User interface
18+
Within Slicer, the "SPECT reconstruction" module is located under the parent category "Tomographic Reconstruction". For more details, see the [user manual](https://github.com/PyTomography/Slicer/blob/main/User_Manual.md)
1819

1920
- Inputs
20-
- Input volume: input SPECT/CT dicom files, simind file (convert to dicom using the data converter)
21+
- SPECT projection data
22+
- (optional) CT data for attenuation correction
2123
- Outputs
22-
- Reconstructed volume: The volume will be saved under the specified name (or as the dataset name appended with _reconstructed) and will be located within the Subject Hierarchy in the Data Module.
24+
- a 3D reconstruction SPECT image. The volume will be saved under the specified name (or as the dataset name appended with _reconstructed) and will be located within the Subject Hierarchy in the Data Module.
2325

2426
## Resources
2527

@@ -47,4 +49,4 @@ SlicerSPECTRecon is subject to the `MIT License`, which is in the project's root
4749

4850
## Contact
4951

50-
Please post any questions to the [Pytomography Discourse Forum](https://pytomography.discourse.group/).
52+
Please post any questions to the [Pytomography Discourse Forum](https://pytomography.discourse.group/).

SlicerSPECTRecon/CMakeLists.txt renamed to SPECTRecon/CMakeLists.txt

Lines changed: 45 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -1,45 +1,45 @@
1-
#-----------------------------------------------------------------------------
2-
set(MODULE_NAME SlicerSPECTRecon)
3-
4-
#-----------------------------------------------------------------------------
5-
set(MODULE_PYTHON_SCRIPTS
6-
${MODULE_NAME}.py
7-
Logic/Algorithms.py
8-
Logic/Callbacks.py
9-
Logic/MetadataUtils.py
10-
Logic/Priors.py
11-
Logic/SimindToDicom.py
12-
Logic/SlicerSPECTReconLogic.py
13-
Logic/VolumeUtils.py
14-
Logic/VtkkmrmlUtils.py
15-
Testing/Python/reconstruction.py
16-
Testing/Python/simind_to_dicom.py
17-
Testing/Python/utils.py
18-
)
19-
20-
set(MODULE_PYTHON_RESOURCES
21-
Resources/Icons/${MODULE_NAME}.png
22-
Resources/UI/${MODULE_NAME}.ui
23-
Resources/algorithmTestSettings.json
24-
Resources/psfMeta.json
25-
Resources/sampleDataMetaData.json
26-
)
27-
28-
#-----------------------------------------------------------------------------
29-
slicerMacroBuildScriptedModule(
30-
NAME ${MODULE_NAME}
31-
SCRIPTS ${MODULE_PYTHON_SCRIPTS}
32-
RESOURCES ${MODULE_PYTHON_RESOURCES}
33-
WITH_GENERIC_TESTS
34-
)
35-
36-
#-----------------------------------------------------------------------------
37-
if(BUILD_TESTING)
38-
39-
# Register the unittest subclass in the main script as a ctest.
40-
# Note that the test will also be available at runtime.
41-
# slicer_add_python_unittest(SCRIPT ${MODULE_NAME}.py)
42-
43-
# Additional build-time testing
44-
add_subdirectory(Logic)
45-
endif()
1+
#-----------------------------------------------------------------------------
2+
set(MODULE_NAME SPECTRecon)
3+
4+
#-----------------------------------------------------------------------------
5+
set(MODULE_PYTHON_SCRIPTS
6+
${MODULE_NAME}.py
7+
Logic/Algorithms.py
8+
Logic/Callbacks.py
9+
Logic/MetadataUtils.py
10+
Logic/Priors.py
11+
Logic/SimindToDicom.py
12+
Logic/SPECTReconLogic.py
13+
Logic/VolumeUtils.py
14+
Logic/VtkkmrmlUtils.py
15+
Testing/Python/reconstruction.py
16+
Testing/Python/simind_to_dicom.py
17+
Testing/Python/utils.py
18+
)
19+
20+
set(MODULE_PYTHON_RESOURCES
21+
Resources/Icons/${MODULE_NAME}.png
22+
Resources/UI/${MODULE_NAME}.ui
23+
Resources/algorithmTestSettings.json
24+
Resources/psfMeta.json
25+
Resources/sampleDataMetaData.json
26+
)
27+
28+
#-----------------------------------------------------------------------------
29+
slicerMacroBuildScriptedModule(
30+
NAME ${MODULE_NAME}
31+
SCRIPTS ${MODULE_PYTHON_SCRIPTS}
32+
RESOURCES ${MODULE_PYTHON_RESOURCES}
33+
WITH_GENERIC_TESTS
34+
)
35+
36+
#-----------------------------------------------------------------------------
37+
if(BUILD_TESTING)
38+
39+
# Register the unittest subclass in the main script as a ctest.
40+
# Note that the test will also be available at runtime.
41+
# slicer_add_python_unittest(SCRIPT ${MODULE_NAME}.py)
42+
43+
# Additional build-time testing
44+
add_subdirectory(Logic)
45+
endif()

‎SPECTRecon/Logic/AdvancedPSF.py‎

Lines changed: 80 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,80 @@
1+
import numpy as np
2+
import torch
3+
from pathlib import Path
4+
from spectpsftoolbox.kernel1d import ArbitraryKernel1D, FunctionKernel1D
5+
from spectpsftoolbox.operator2d import GaussianOperator, Rotate1DConvOperator, RotateSeperable2DConvOperator
6+
import pytomography
7+
device = pytomography.device
8+
9+
psf_operator_dict = {
10+
'GI-HEGP @ 440 keV': 'params_440keV_GIHEGP_256shift_sparse_DW_pyrex66_1DR',
11+
'SY-HE @ 440 keV': 'params_440keV_SYHE_256shift_sparse_DW_pyrex66_1DR'
12+
}
13+
14+
def getScriptPath():
15+
try:
16+
script_path = Path(__file__).resolve().parents[1]
17+
except:
18+
from inspect import getsourcefile
19+
script_path = Path(getsourcefile(lambda: 0)).resolve().parents[1]
20+
finally:
21+
return script_path
22+
23+
24+
def get_psf_operator(psf_name):
25+
script_path = getScriptPath()
26+
path = script_path / 'Resources' / 'Data' / 'PSF_models' / f'{psf_operator_dict[psf_name]}_'
27+
path = str(path)
28+
gauss_amplitude_params = torch.load(path+'gauss_amplitude.pt')
29+
gauss_sigma_params = torch.load(path+'gauss_sigma.pt')
30+
tail_amplitude_params = torch.load(path+'tail_amplitude.pt')
31+
tail_sigma_params = torch.load(path+'tail_sigma.pt')
32+
bkg_amplitude_params = torch.load(path+'bkg_amplitude.pt')
33+
bkg_sigma_params = torch.load(path+'bkg_sigma.pt')
34+
Nx0 = 255
35+
tail_kernel_decay = 0.5
36+
if "_SY" in path:
37+
rot = 90
38+
dx0 = 0.24
39+
elif "_GI" in path:
40+
rot = 0
41+
dx0 = 0.2208336
42+
a_min = 1.
43+
a_max = 55.
44+
# Gaussian component
45+
gauss_amplitude_fn = lambda a, bs: bs[0]*torch.exp(-a*bs[1]) + bs[2]*torch.exp(-a*bs[3])
46+
gauss_sigma_fn = lambda a, bs: bs[0] + bs[1]*(torch.sqrt(a**2 + bs[2]**2) - torch.abs(bs[2]))
47+
gaussian_operator = GaussianOperator(
48+
gauss_amplitude_fn,
49+
gauss_sigma_fn,
50+
gauss_amplitude_params,
51+
gauss_sigma_params,
52+
)
53+
# Tail component
54+
Nx_tail = round(np.sqrt(2)*Nx0)
55+
if Nx_tail%2==0:
56+
Nx_tail += 1
57+
x = torch.arange(-(Nx_tail-1)/2, (Nx_tail+1)/2, 1).to(device) * dx0
58+
tail_kernel = torch.tensor(torch.exp(-tail_kernel_decay*torch.abs(x)), requires_grad=True, device=device)
59+
tail_amplitude_fn = lambda a, bs: bs[0]*torch.exp(-a*bs[1]) + bs[2]*torch.exp(-a*bs[3])
60+
tail_sigma_fn = lambda a, bs, a_min=a_min: 1 + bs[0]*(torch.sqrt((a-a_min)**2 + bs[1]**2) - torch.abs(bs[1]))
61+
tail_kernel1D = ArbitraryKernel1D(tail_kernel, tail_amplitude_fn, tail_sigma_fn, tail_amplitude_params, tail_sigma_params, dx0, grid_sample_mode='bicubic')
62+
tail_operator = Rotate1DConvOperator(
63+
tail_kernel1D,
64+
N_angles = 3,
65+
additive=True,
66+
rot=rot,
67+
)
68+
# Bkg component
69+
bkg_amplitude_fn = lambda a, bs: bs[0]*torch.exp(-a*bs[1]) + bs[2]*torch.exp(-a*bs[3])
70+
bkg_sigma_fn = lambda a, bs: bs[0] + bs[1]*(torch.sqrt(a**2 + bs[2]**2) - torch.abs(bs[2]))
71+
kernel_fn = lambda x: torch.exp(-torch.abs(x))
72+
bkg_kernel1D = FunctionKernel1D(kernel_fn, bkg_amplitude_fn, bkg_sigma_fn, bkg_amplitude_params, bkg_sigma_params, a_min=a_min, a_max=a_max)
73+
bkg_operator = RotateSeperable2DConvOperator(
74+
bkg_kernel1D,
75+
N_angles = 1,
76+
additive=False
77+
)
78+
# Operator
79+
psf_operator = (tail_operator + bkg_operator) * gaussian_operator + gaussian_operator
80+
return psf_operator

SlicerSPECTRecon/Logic/Algorithms.py renamed to SPECTRecon/Logic/Algorithms.py

Lines changed: 12 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,12 @@
1-
from pytomography.algorithms import OSEM, BSREM, OSMAPOSL
2-
3-
def selectAlgorithm(algorithm, likelihood, prior):
4-
if algorithm == "OSEM":
5-
reconstruction_algorithm = OSEM(likelihood)
6-
return reconstruction_algorithm
7-
elif algorithm == "BSREM":
8-
reconstruction_algorithm = BSREM(likelihood, prior=prior)
9-
return reconstruction_algorithm
10-
elif algorithm == "OSMAPOSL":
11-
reconstruction_algorithm = OSMAPOSL(likelihood, prior=prior)
12-
return reconstruction_algorithm
1+
from pytomography.algorithms import OSEM, BSREM, OSMAPOSL
2+
3+
def selectAlgorithm(algorithm, likelihood, prior):
4+
if algorithm == "OSEM":
5+
reconstruction_algorithm = OSEM(likelihood)
6+
return reconstruction_algorithm
7+
elif algorithm == "BSREM":
8+
reconstruction_algorithm = BSREM(likelihood, prior=prior)
9+
return reconstruction_algorithm
10+
elif algorithm == "OSMAPOSL":
11+
reconstruction_algorithm = OSMAPOSL(likelihood, prior=prior)
12+
return reconstruction_algorithm

SlicerSPECTRecon/Logic/CMakeLists.txt renamed to SPECTRecon/Logic/CMakeLists.txt

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,2 +1,2 @@
1-
slicer_add_python_unittest(SCRIPT ./../Testing/Python/reconstruction.py SLICER_ARGS --no-main-window)
1+
slicer_add_python_unittest(SCRIPT ./../Testing/Python/reconstruction.py SLICER_ARGS --no-main-window)
22
slicer_add_python_unittest(SCRIPT ./../Testing/Python/simind_To_dicom.py SLICER_ARGS --no-main-window)

SlicerSPECTRecon/Logic/Callbacks.py renamed to SPECTRecon/Logic/Callbacks.py

Lines changed: 25 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -1,26 +1,26 @@
1-
import slicer
2-
from pytomography.callbacks import Callback, DataStorageCallback
3-
4-
class LoadingCallback(Callback):
5-
def __init__(self, progressDialog, n_iters, n_subsets):
6-
self.progressDialog = progressDialog
7-
self.total_subiterations = n_iters*n_subsets
8-
self.n_subsets = n_subsets
9-
10-
def run(self, object, n_iter, n_subset):
11-
self.progressDialog.value = int((n_iter*self.n_subsets+n_subset+1)/self.total_subiterations * 100)
12-
slicer.app.processEvents()
13-
return object
14-
15-
class DataStorageWithLoadingCallback(DataStorageCallback):
16-
def __init__(self, likelihood, object_prediction, progressDialog, n_iters, n_subsets):
17-
super().__init__(likelihood, object_prediction)
18-
self.progressDialog = progressDialog
19-
self.total_subiterations = n_iters*n_subsets
20-
self.n_subsets = n_subsets
21-
22-
def run(self, object, n_iter, n_subset):
23-
object = super().run(object, n_iter, n_subset)
24-
self.progressDialog.value = int((n_iter*self.n_subsets+n_subset+1)/self.total_subiterations * 100)
25-
slicer.app.processEvents()
1+
import slicer
2+
from pytomography.callbacks import Callback, DataStorageCallback
3+
4+
class LoadingCallback(Callback):
5+
def __init__(self, progressDialog, n_iters, n_subsets):
6+
self.progressDialog = progressDialog
7+
self.total_subiterations = n_iters*n_subsets
8+
self.n_subsets = n_subsets
9+
10+
def run(self, object, n_iter, n_subset):
11+
self.progressDialog.value = int((n_iter*self.n_subsets+n_subset+1)/self.total_subiterations * 100)
12+
slicer.app.processEvents()
13+
return object
14+
15+
class DataStorageWithLoadingCallback(DataStorageCallback):
16+
def __init__(self, likelihood, object_prediction, progressDialog, n_iters, n_subsets):
17+
super().__init__(likelihood, object_prediction)
18+
self.progressDialog = progressDialog
19+
self.total_subiterations = n_iters*n_subsets
20+
self.n_subsets = n_subsets
21+
22+
def run(self, object, n_iter, n_subset):
23+
object = super().run(object, n_iter, n_subset)
24+
self.progressDialog.value = int((n_iter*self.n_subsets+n_subset+1)/self.total_subiterations * 100)
25+
slicer.app.processEvents()
2626
return object

SlicerSPECTRecon/Logic/MetadataUtils.py renamed to SPECTRecon/Logic/MetadataUtils.py

Lines changed: 38 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -1,38 +1,39 @@
1-
from pytomography.io.SPECT import dicom
2-
import numpy as np
3-
import pydicom
4-
5-
def getEnergyWindow(directory):
6-
ds = pydicom.dcmread(directory)
7-
window_names = []
8-
mean_window_energies = []
9-
for energy_window_information in ds.EnergyWindowInformationSequence:
10-
lower_limit = energy_window_information.EnergyWindowRangeSequence[0].EnergyWindowLowerLimit
11-
upper_limit = energy_window_information.EnergyWindowRangeSequence[0].EnergyWindowUpperLimit
12-
energy_window_name = energy_window_information.EnergyWindowName
13-
mean_window_energies.append((lower_limit+upper_limit)/2)
14-
window_names.append(f'{energy_window_name} ({lower_limit:.2f}keV - {upper_limit:.2f}keV)')
15-
idx_sorted = list(np.argsort(mean_window_energies))
16-
window_names = list(np.array(window_names)[idx_sorted])
17-
mean_window_energies = list(np.array(mean_window_energies)[idx_sorted])
18-
# Insert None option to beginning
19-
window_names.insert(0, 'None')
20-
idx_sorted.insert(0, None)
21-
return window_names, mean_window_energies, idx_sorted
22-
23-
def get_object_meta_proj_meta(bed_idx, files_NM, index_peak):
24-
file_NM = files_NM[bed_idx]
25-
object_meta, proj_meta = dicom.get_metadata(file_NM, index_peak)
26-
return object_meta, proj_meta
27-
28-
def get_photopeak_scatter(bed_idx, files_NM, index_peak, index_lower=None, index_upper=None):
29-
projectionss = dicom.load_multibed_projections(files_NM)
30-
photopeak = projectionss[bed_idx][index_peak]
31-
# No scatter
32-
if (index_lower is None)*(index_upper is None):
33-
scatter = None
34-
# Dual or triple energy window
35-
else:
36-
file_NM = files_NM[bed_idx]
37-
scatter = dicom.get_energy_window_scatter_estimate_projections(file_NM, projectionss[bed_idx], index_peak, index_lower, index_upper)
1+
from pytomography.io.SPECT import dicom
2+
import numpy as np
3+
import pydicom
4+
5+
def getEnergyWindow(directory):
6+
ds = pydicom.dcmread(directory)
7+
window_names = []
8+
mean_window_energies = []
9+
for energy_window_information in ds.EnergyWindowInformationSequence:
10+
lower_limit = energy_window_information.EnergyWindowRangeSequence[0].EnergyWindowLowerLimit
11+
upper_limit = energy_window_information.EnergyWindowRangeSequence[0].EnergyWindowUpperLimit
12+
energy_window_name = energy_window_information.EnergyWindowName
13+
mean_window_energies.append((lower_limit+upper_limit)/2)
14+
window_names.append(f'{energy_window_name} ({lower_limit:.2f}keV - {upper_limit:.2f}keV)')
15+
idx_sorted = list(np.argsort(mean_window_energies))
16+
window_names = list(np.array(window_names)[idx_sorted])
17+
mean_window_energies = list(np.array(mean_window_energies)[idx_sorted])
18+
# Insert None option to beginning
19+
window_names.insert(0, 'None')
20+
idx_sorted.insert(0, None)
21+
return window_names, mean_window_energies, idx_sorted
22+
23+
def get_object_meta_proj_meta(bed_idx, files_NM, index_peak):
24+
file_NM = files_NM[bed_idx]
25+
object_meta, proj_meta = dicom.get_metadata(file_NM, index_peak)
26+
return object_meta, proj_meta
27+
28+
def get_photopeak_scatter(bed_idx, files_NM, index_peak, index_lower=None, index_upper=None, scatter_sigma=0.0):
29+
object_meta, proj_meta = dicom.get_metadata(files_NM[0], index_peak)
30+
projectionss = dicom.load_multibed_projections(files_NM)
31+
photopeak = projectionss[bed_idx][index_peak]
32+
# No scatter
33+
if (index_lower is None)*(index_upper is None):
34+
scatter = None
35+
# Dual or triple energy window
36+
else:
37+
file_NM = files_NM[bed_idx]
38+
scatter = dicom.get_energy_window_scatter_estimate_projections(file_NM, projectionss[bed_idx], index_peak, index_lower, index_upper, sigma_r=scatter_sigma, sigma_z=scatter_sigma, proj_meta=proj_meta)
3839
return photopeak, scatter

0 commit comments

Comments
 (0)