Newer
Older
"""
Compute reconstructions with CTAC and DLAC (attenuation maps created by `compute_mu_maps`).
This script implements the (post-reconstructions attenuation correction) PRAC algorithm
by first projection a non-attenuation-corrected reconstruction and then reconstructing this
projection with an attenuation map.
"""
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
import argparse
import os
import pydicom
import stir
stir.Verbosity_set(0)
from mu_map.file.dicom import change_uid, load_dcm, update_dcm
from mu_map.file.util import load_as_interfile
from mu_map.logging import add_logging_args, get_logger_by_args
from mu_map.recon.filter import GaussianFilter
from mu_map.recon.osem import reconstruct
from mu_map.recon.project import forward_project
parser = argparse.ArgumentParser()
parser.add_argument(
"--mu_maps",
type=str,
required=True,
help="directory of mu maps computed with `mu_map.scripts.compute_mu_maps`",
)
parser.add_argument(
"--out_dir",
type=str,
required=True,
help="directory where resulting reconstructions are written to",
)
parser.add_argument(
"--description_postfix",
type=str,
help="postfix for the DICOM SeriesDescription",
)
parser.add_argument(
"--dataset_dir",
type=str,
default="data/second/",
help="directory of the dataset with reconstructions and real mu maps",
)
parser.add_argument(
"--images_dir",
type=str,
default="images",
help="directory under `--dataset_dir` actually containing the images",
)
parser.add_argument(
"--n_subsets",
type=int,
default=3,
help="number of subsets used for the OSEM reconstruction",
)
parser.add_argument(
"--n_iterations",
type=int,
default=10,
help="number of iterations performed by the OSEM reconstruction",
)
add_logging_args(parser, defaults={"--logfile": "compute_recons.log"})
args = parser.parse_args()
args.images_dir = os.path.join(args.dataset_dir, args.images_dir)
args.description_postfix = (
"" if args.description_postfix is None else f" - {args.description_postfix}"
)
args.logfile = os.path.join(args.out_dir, args.logfile)
logger = get_logger_by_args(args)
logger.info(args)
files_mu_maps_dl = os.listdir(args.mu_maps)
files_mu_maps_dl = filter(lambda f: f.split(".")[-1] == "dcm", files_mu_maps_dl)
files_mu_maps_dl = list(files_mu_maps_dl)
files_mu_maps_dl.sort()
for file_mu_map_dl in files_mu_maps_dl:
logger.info(f"Process {file_mu_map_dl}")
file_recon_ac = file_mu_map_dl.replace("mu_map_dl", "recon_ac_nsc")
file_recon_ac = os.path.join(args.images_dir, file_recon_ac)
recon_ac_dcm, _ = load_dcm(file_recon_ac, direction=1)
logger.info(f"Found AC recon: {file_recon_ac}")
file_recon_nac = file_mu_map_dl.replace("mu_map_dl", "recon_nac_nsc")
file_recon_nac = os.path.join(args.images_dir, file_recon_nac)
recon_nac = load_as_interfile(file_recon_nac, direction=1)
logger.info(f"Found nAC recon: {file_recon_nac}")
file_mu_map_ct = file_mu_map_dl.replace("mu_map_dl", "mu_map")
file_mu_map_ct = os.path.join(args.images_dir, file_mu_map_ct)
mu_map_ct = load_as_interfile(file_mu_map_ct, direction=1)
logger.info(f"Found CT-based mu map: {file_mu_map_ct}")
mu_map_dl = load_as_interfile(
os.path.join(args.mu_maps, file_mu_map_dl), direction=1
)
for mu_map, label in zip([mu_map_ct, mu_map_dl], ["CTAC", "DLAC"]):
logger.info(f"Compute {label} reconstruction")
projection = forward_project(recon_nac, n_slices=mu_map.image.shape[0])
_recon_ac = reconstruct(
projection,
mu_map=mu_map,
n_subsets=args.n_subsets,
n_iterations=args.n_iterations,
postfilter=GaussianFilter(projection),
)
# for some reason STIR reconstruction are inverted in the x-direction
# this is reverted here
_recon_ac.image = _recon_ac.image[:, :, ::-1]
_recon_ac_dcm = change_uid(update_dcm(recon_ac_dcm.copy(), _recon_ac.image))
_recon_ac_dcm.SeriesDescription = _recon_ac_dcm.SeriesDescription.replace(
"NoSC - AC ", f"{label}{args.description_postfix}"
)
_file_recon_ac = file_mu_map_dl.replace("mu_map_dl", f"recon_{label.lower()}")
_file_recon_ac = os.path.join(args.out_dir, _file_recon_ac)
pydicom.dcmwrite(_file_recon_ac, _recon_ac_dcm)
logger.info(f"Store {label} recon at {_file_recon_ac}")