Skip to content

Commit e007ee5

Browse files
KrisThielemansKrisThielemans
authored
Merge pull request #168 from ETSInitiative/different_module_types
generator/analysis helpers adaptations for multiple module types
2 parents f582758 + 7112ab5 commit e007ee5

File tree

4 files changed

+448
-169
lines changed

4 files changed

+448
-169
lines changed

‎cpp/helpers/petsird_analysis.cpp‎

Lines changed: 54 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@ using petsird::binary::PETSIRDReader;
2222
#include <iostream>
2323
#include <variant>
2424
#include <cstdlib>
25+
#include <vector>
2526

2627
void
2728
print_usage_and_exit(char const* prog_name)
@@ -107,53 +108,59 @@ main(int argc, char const* argv[])
107108
std::cout << "Subject ID: " << header.exam->subject.id << std::endl;
108109
const auto num_module_types = scanner.scanner_geometry.replicated_modules.size();
109110
std::cout << "Types of modules: " << num_module_types << std::endl;
110-
// TODO more than 1 type of module
111-
const petsird::TypeOfModule type_of_module{ 0 };
112-
std::cout << "Number of modules of first type: "
113-
<< scanner.scanner_geometry.replicated_modules[type_of_module].transforms.size() << std::endl;
114-
std::cout << "Number of elements in modules of first type: "
115-
<< scanner.scanner_geometry.replicated_modules[type_of_module].object.detecting_elements.transforms.size()
116-
<< std::endl;
117-
std::cout << "Total number of 'crystals' in modules of first type : "
118-
<< petsird_helpers::get_num_det_els(scanner, type_of_module) << std::endl;
119-
120-
const auto& tof_bin_edges = scanner.tof_bin_edges[type_of_module][type_of_module];
121-
const auto num_tof_bins = tof_bin_edges.NumberOfBins();
122-
std::cout << "Number of TOF bins: " << num_tof_bins << std::endl;
123-
std::cout << "TOF bin edges: " << tof_bin_edges.edges << std::endl;
124-
const auto& event_energy_bin_edges = scanner.event_energy_bin_edges[type_of_module];
125-
const auto num_event_energy_bins = event_energy_bin_edges.NumberOfBins();
126-
std::cout << "Number of energy bins: " << num_event_energy_bins << std::endl;
127-
std::cout << "Event energy bin edges: " << event_energy_bin_edges.edges << std::endl;
128-
const auto energy_mid_points = (xt::view(event_energy_bin_edges.edges, xt::range(0, event_energy_bin_edges.edges.size() - 1))
129-
+ xt::view(event_energy_bin_edges.edges, xt::range(1, event_energy_bin_edges.edges.size())))
130-
/ 2;
131-
std::cout << "Event energy mid points: " << energy_mid_points << std::endl;
132-
133-
std::cout << "Singles Histogram Level: ";
134-
switch (scanner.singles_histogram_level)
111+
std::vector<yardl::NDArray<float, 1>> all_energy_mid_points;
112+
for (petsird::TypeOfModule type_of_module = 0; type_of_module < num_module_types; ++type_of_module)
135113
{
136-
case petsird::SinglesHistogramLevelType::kNone:
137-
std::cout << "none\n";
138-
break;
139-
case petsird::SinglesHistogramLevelType::kModule:
140-
std::cout << "module\n";
141-
break;
142-
case petsird::SinglesHistogramLevelType::kAll:
143-
std::cout << "all\n";
144-
break;
145-
}
146-
if (scanner.singles_histogram_level != petsird::SinglesHistogramLevelType::kNone)
147-
{
148-
const auto& singles_histogram_energy_bin_edges = scanner.singles_histogram_energy_bin_edges[type_of_module];
149-
std::cout << "Singles Histogram Energy Bin Edges: " << singles_histogram_energy_bin_edges.edges << std::endl;
150-
std::cout << "Number of Singles Histogram Energy Windows: " << singles_histogram_energy_bin_edges.NumberOfBins()
114+
std::cout << "------ Module type " << type_of_module << std::endl;
115+
std::cout << "Number of modules of this type: "
116+
<< scanner.scanner_geometry.replicated_modules[type_of_module].transforms.size() << std::endl;
117+
std::cout << "Number of elements in modules of this type: "
118+
<< scanner.scanner_geometry.replicated_modules[type_of_module].object.detecting_elements.transforms.size()
151119
<< std::endl;
152-
}
120+
std::cout << "Total number of 'crystals' in modules of this type : "
121+
<< petsird_helpers::get_num_det_els(scanner, type_of_module) << std::endl;
153122

154-
petsird::TimeBlock time_block;
123+
const auto& tof_bin_edges = scanner.tof_bin_edges[type_of_module][type_of_module];
124+
const auto num_tof_bins = tof_bin_edges.NumberOfBins();
125+
std::cout << "Number of TOF bins: " << num_tof_bins << std::endl;
126+
std::cout << "TOF bin edges: " << tof_bin_edges.edges << std::endl;
127+
const auto& event_energy_bin_edges = scanner.event_energy_bin_edges[type_of_module];
128+
const auto num_event_energy_bins = event_energy_bin_edges.NumberOfBins();
129+
std::cout << "Number of energy bins: " << num_event_energy_bins << std::endl;
130+
std::cout << "Event energy bin edges: " << event_energy_bin_edges.edges << std::endl;
131+
const auto energy_mid_points
132+
= (xt::view(event_energy_bin_edges.edges, xt::range(0, event_energy_bin_edges.edges.size() - 1))
133+
+ xt::view(event_energy_bin_edges.edges, xt::range(1, event_energy_bin_edges.edges.size())))
134+
/ 2;
135+
std::cout << "Event energy mid points: " << energy_mid_points << std::endl;
136+
all_energy_mid_points.push_back(energy_mid_points);
137+
138+
std::cout << "Singles Histogram Level: ";
139+
switch (scanner.singles_histogram_level)
140+
{
141+
case petsird::SinglesHistogramLevelType::kNone:
142+
std::cout << "none\n";
143+
break;
144+
case petsird::SinglesHistogramLevelType::kModule:
145+
std::cout << "module\n";
146+
break;
147+
case petsird::SinglesHistogramLevelType::kAll:
148+
std::cout << "all\n";
149+
break;
150+
}
151+
if (scanner.singles_histogram_level != petsird::SinglesHistogramLevelType::kNone)
152+
{
153+
const auto& singles_histogram_energy_bin_edges = scanner.singles_histogram_energy_bin_edges[type_of_module];
154+
std::cout << "Singles Histogram Energy Bin Edges: " << singles_histogram_energy_bin_edges.edges << std::endl;
155+
std::cout << "Number of Singles Histogram Energy Windows: " << singles_histogram_energy_bin_edges.NumberOfBins()
156+
<< std::endl;
157+
}
158+
} // loop over type_of_module
159+
160+
std::cout << "------------------------- " << std::endl;
155161

156-
// Process events in batches of up to 100
162+
// Now read events and print some things
163+
petsird::TimeBlock time_block;
157164
float energy_1 = 0, energy_2 = 0;
158165
std::size_t num_prompts = 0;
159166
std::size_t num_delayeds = 0;
@@ -170,9 +177,11 @@ main(int argc, char const* argv[])
170177

171178
for (unsigned mtype0 = 0; mtype0 < num_module_types; ++mtype0)
172179
{
180+
const auto& energy_mid_points0 = all_energy_mid_points[mtype0];
173181
for (unsigned mtype1 = 0; mtype1 < num_module_types; ++mtype1)
174182
{
175183
const petsird::TypeOfModulePair mtype_pair{ mtype0, mtype1 };
184+
const auto& energy_mid_points1 = all_energy_mid_points[mtype1];
176185

177186
// This code would need work to be able to handle a list-mode file without prompts
178187
const auto& prompt_events = event_time_block.prompt_events[mtype0][mtype1];
@@ -198,8 +207,8 @@ main(int argc, char const* argv[])
198207
= petsird_helpers::expand_detection_bin(scanner, mtype1, event.detection_bins[1]);
199208

200209
// accumulate energies to print average below
201-
energy_1 += energy_mid_points[expanded_det_bin0.energy_index];
202-
energy_2 += energy_mid_points[expanded_det_bin1.energy_index];
210+
energy_1 += energy_mid_points0[expanded_det_bin0.energy_index];
211+
energy_2 += energy_mid_points1[expanded_det_bin1.energy_index];
203212

204213
if (print_events)
205214
{

‎python/petsird/helpers/analysis.py‎

Lines changed: 48 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66
import argparse
77
import sys
88

9+
import numpy as np
910
import petsird
1011
import petsird.helpers.geometry
1112
from petsird.helpers import (expand_detection_bin, get_detection_efficiency,
@@ -44,41 +45,57 @@ def parserCreator():
4445
print(f"Subject ID: {header.exam.subject.id}")
4546
print(f"Scanner name: {scanner.model_name}")
4647
num_module_types = len(scanner.scanner_geometry.replicated_modules)
48+
all_energy_mid_points = []
4749
print("Types of modules: ", num_module_types)
48-
type_of_module = 0
49-
print(
50-
"Number of modules of first type: ",
51-
len(scanner.scanner_geometry.replicated_modules[type_of_module].
52-
transforms))
53-
print(
54-
"Number of elements in modules of first type: ",
55-
len(scanner.scanner_geometry.replicated_modules[type_of_module].
56-
object.detecting_elements.transforms))
57-
print("Total number of 'crystals': ",
58-
get_num_det_els(scanner, type_of_module))
59-
60-
tof_bin_edges = scanner.tof_bin_edges[type_of_module][type_of_module]
61-
num_tof_bins = tof_bin_edges.number_of_bins()
62-
event_energy_bin_edges = scanner.event_energy_bin_edges[type_of_module]
63-
num_event_energy_bins = event_energy_bin_edges.number_of_bins()
64-
print("Number of TOF bins: ", num_tof_bins)
65-
print("Tof bin edges: ", tof_bin_edges)
66-
print("Number of energy bins: ", num_event_energy_bins)
67-
event_energy_bin_edges = event_energy_bin_edges.edges
68-
print("Event energy bin edges: ", event_energy_bin_edges)
69-
energy_mid_points = (event_energy_bin_edges[:-1] +
70-
event_energy_bin_edges[1:]) / 2
71-
print("Event energy mid points: ", energy_mid_points)
50+
for type_of_module in range(num_module_types):
51+
print("------ Module type ", type_of_module)
52+
print(
53+
"Number of modules of this type: ",
54+
len(scanner.scanner_geometry.
55+
replicated_modules[type_of_module].transforms))
56+
print(
57+
"Number of elements in modules of this type: ",
58+
len(scanner.scanner_geometry.replicated_modules[type_of_module]
59+
.object.detecting_elements.transforms))
60+
print("Total number of 'crystals': ",
61+
get_num_det_els(scanner, type_of_module))
62+
63+
tof_bin_edges = scanner.tof_bin_edges[type_of_module][
64+
type_of_module]
65+
num_tof_bins = tof_bin_edges.number_of_bins()
66+
event_energy_bin_edges = scanner.event_energy_bin_edges[
67+
type_of_module]
68+
num_event_energy_bins = event_energy_bin_edges.number_of_bins()
69+
print("Number of TOF bins: ", num_tof_bins)
70+
print("Tof bin edges: ", tof_bin_edges)
71+
print("Number of energy bins: ", num_event_energy_bins)
72+
event_energy_bin_edges = event_energy_bin_edges.edges
73+
print("Event energy bin edges: ", event_energy_bin_edges)
74+
energy_mid_points = (event_energy_bin_edges[:-1] +
75+
event_energy_bin_edges[1:]) / 2
76+
all_energy_mid_points.append(energy_mid_points)
77+
print("Event energy mid points: ", energy_mid_points)
78+
if scanner.detection_efficiencies.module_pair_sgidlut is not None:
79+
for type_of_module0 in range(num_module_types):
80+
for type_of_module1 in range(num_module_types):
81+
print("------ Module type pair ", type_of_module0,
82+
type_of_module1)
83+
with np.printoptions(threshold=2000, linewidth=140):
84+
print(
85+
"SGID LUT:\n",
86+
scanner.detection_efficiencies.module_pair_sgidlut[
87+
type_of_module0][type_of_module1])
7288
print("Singles histogram level: ", scanner.singles_histogram_level)
7389
if scanner.singles_histogram_level != petsird.SinglesHistogramLevelType.NONE:
7490
print(
75-
"Number of singles histograms energy windows for first module-type: ",
91+
"Number of singles histograms energy windows for this module-type: ",
7692
scanner.singles_histogram_energy_bin_edges[type_of_module].
7793
number_of_bins())
7894
print("Singles histogram energy bin edges: ",
7995
scanner.singles_histogram_energy_bin_edges)
80-
print("SGID LUT:\n",
81-
scanner.detection_efficiencies.module_pair_sgidlut)
96+
print("------------------------- ")
97+
98+
# Now read events and print some things
8299
energy_1, energy_2 = 0.0, 0.0
83100
num_prompts = 0
84101
num_delayeds = 0
@@ -87,7 +104,9 @@ def parserCreator():
87104
if isinstance(time_block, petsird.TimeBlock.EventTimeBlock):
88105
last_time = time_block.value.time_interval.stop
89106
for mtype0 in range(num_module_types):
107+
energy_mid_points0 = all_energy_mid_points[mtype0]
90108
for mtype1 in range(num_module_types):
109+
energy_mid_points1 = all_energy_mid_points[mtype1]
91110
mtype_pair = petsird.TypeOfModulePair((mtype0, mtype1))
92111

93112
# count events
@@ -111,9 +130,9 @@ def parserCreator():
111130
scanner, mtype1, event.detection_bins[1])
112131

113132
# accumulate energies to print average below
114-
energy_1 += float(energy_mid_points[
133+
energy_1 += float(energy_mid_points0[
115134
expanded_det_bin0.energy_index])
116-
energy_2 += float(energy_mid_points[
135+
energy_2 += float(energy_mid_points1[
117136
expanded_det_bin1.energy_index])
118137
if print_events:
119138
print(event)

0 commit comments

Comments
 (0)