Skip to content

Commit 7112ab5

Browse files
KrisThielemansKrisThielemans
committed
updated C++ analysis helper for multiple module-types
1 parent 29b7d10 commit 7112ab5

File tree

1 file changed

+54
-45
lines changed

1 file changed

+54
-45
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
{

0 commit comments

Comments
 (0)