Merge pull request #74 from KrisThielemans/analysis_improvements

extra features for petsird_analysis

cpp/petsird_analysis.cpp

@@ -20,23 +20,73 @@ using petsird::binary::PETSIRDReader;
 #include <xtensor/xio.hpp>
 #include <iostream>
 #include <variant>
+#include <cstdlib>
+
+void
+print_usage_and_exit(char const* prog_name)
+{
+  std::cerr << "Usage:\n"
+            << prog_name << " \\\n"
+            << "    [--print_events] [--input petsird_filename]\n\n"
+            << "Options:\n"
+            << "    -e, --print-events: Print event info\n"
+            << "    -i, --input       : Filename to read\n\n"
+            << "Currently, the following (deprecated) usage is also allowed:\n"
+            << prog_name << " [options] [---] petsird_filename\n"
+            << "Use of '--' is then required if the filename starts with -\n";
+
+  std::exit(EXIT_FAILURE);
+}
 
 int
-main(int argc, char* argv[])
+main(int argc, char const* argv[])
 {
+  auto prog_name = argv[0];
+
+  bool print_events = false;
+  std::string filename;
+  // option processing
+  while (argc > 1 && (strncmp(argv[1], "-", 1) == 0))
+    {
+      if (strcmp(argv[1], "--print_events") == 0 || strcmp(argv[1], "-e") == 0)
+        {
+          print_events = true;
+        }
+      else if (strcmp(argv[1], "--input") == 0 || strcmp(argv[1], "-i") == 0)
+        {
+          filename = argv[2];
+          ++argv;
+          --argc;
+        }
+      else if (strcmp(argv[1], "--") == 0)
+        {
+          // next argument is filename
+          ++argv;
+          --argc;
+          break;
+        }
+      else
+        print_usage_and_exit(prog_name);
+      ++argv;
+      --argc;
+    }
   // Check if the user has provided a file
-  if (argc < 2)
+  if (!filename.empty() && argc > 1)
+    print_usage_and_exit(prog_name);
+  else if (filename.empty())
     {
-      std::cerr << "Please provide a file to read" << std::endl;
-      return 1;
+      if (argc < 2)
+        print_usage_and_exit(prog_name);
+      else
+        filename = argv[1];
     }
 
   // Open the file
-  PETSIRDReader reader(argv[1]);
+  PETSIRDReader reader(filename);
   petsird::Header header;
   reader.ReadHeader(header);
 
-  std::cout << "Processing file: " << argv[1] << std::endl;
+  std::cout << "Processing file: " << filename << std::endl;
   if (header.exam) // only do this if present
     std::cout << "Subject ID: " << header.exam->subject.id << std::endl;
   std::cout << "Types of modules: " << header.scanner.scanner_geometry.replicated_modules.size() << std::endl;
@@ -65,6 +115,7 @@ main(int argc, char* argv[])
   // Process events in batches of up to 100
   float energy_1 = 0, energy_2 = 0;
   std::size_t num_prompts = 0;
+  std::size_t num_delayeds = 0;
   float last_time = 0.F;
   while (reader.ReadTimeBlocks(time_block))
     {
@@ -73,29 +124,37 @@ main(int argc, char* argv[])
           auto& event_time_block = std::get<petsird::EventTimeBlock>(time_block);
           last_time = event_time_block.start;
           num_prompts += event_time_block.prompt_events.size();
-          std::cout << "=====================  Events in time block from " << last_time << " ==============\n";
+          if (event_time_block.delayed_events)
+            num_delayeds += event_time_block.delayed_events->size();
+          if (print_events)
+            std::cout << "=====================  Prompt events in time block from " << last_time << " ==============\n";
 
           for (auto& event : event_time_block.prompt_events)
             {
               energy_1 += energy_mid_points[event.energy_indices[0]];
               energy_2 += energy_mid_points[event.energy_indices[1]];
 
-              std::cout << "CoincidenceEvent(detectorIds=[" << event.detector_ids[0] << ", " << event.detector_ids[1]
-                        << "], tofIdx=" << event.tof_idx << ", energyIndices=[" << event.energy_indices[0] << ", "
-                        << event.energy_indices[1] << "])\n";
-              const auto module_and_elems
-                  = petsird_helpers::get_module_and_element(header.scanner.scanner_geometry, event.detector_ids);
-              std::cout << "    "
-                        << "[ModuleAndElement(module=" << module_and_elems[0].module << ", "
-                        << "el=" << module_and_elems[0].el << "), ModuleAndElement(module=" << module_and_elems[0].module << ", "
-                        << "el=" << module_and_elems[0].el << ")]\n";
-              std::cout << "    efficiency:" << petsird_helpers::get_detection_efficiency(header.scanner, event) << "\n";
+              if (print_events)
+                {
+                  std::cout << "CoincidenceEvent(detectorIds=[" << event.detector_ids[0] << ", " << event.detector_ids[1]
+                            << "], tofIdx=" << event.tof_idx << ", energyIndices=[" << event.energy_indices[0] << ", "
+                            << event.energy_indices[1] << "])\n";
+                  const auto module_and_elems
+                      = petsird_helpers::get_module_and_element(header.scanner.scanner_geometry, event.detector_ids);
+                  std::cout << "    "
+                            << "[ModuleAndElement(module=" << module_and_elems[0].module << ", "
+                            << "el=" << module_and_elems[0].el << "), ModuleAndElement(module=" << module_and_elems[0].module
+                            << ", "
+                            << "el=" << module_and_elems[0].el << ")]\n";
+                  std::cout << "    efficiency:" << petsird_helpers::get_detection_efficiency(header.scanner, event) << "\n";
+                }
             }
         }
     }
 
   std::cout << "Last time block at " << last_time << " ms\n";
   std::cout << "Number of prompt events: " << num_prompts << std::endl;
+  std::cout << "Number of delayed events: " << num_delayeds << std::endl;
   if (num_prompts > 0)
     {
       std::cout << "Average energy_1: " << energy_1 / num_prompts << std::endl;

python/petsird_analysis.py

@@ -3,14 +3,40 @@
 #
 #  SPDX-License-Identifier: Apache-2.0
 
+import argparse
 import sys
 
 import petsird
 from petsird_helpers import (get_detection_efficiency, get_module_and_element,
                              get_num_det_els)
 
+
+def parserCreator():
+    parser = argparse.ArgumentParser(
+        prog='petsird_analysis',
+        description='Example program that lists basic content of a PETSIRD file'
+    )
+    parser.add_argument('-e', '--print_events', action='store_true')
+    parser.add_argument(
+        "-i",
+        "--input",
+        type=str,
+        default=None,
+        help="File to read from, or stdin if omitted",
+    )
+    return parser.parse_args()
+
+
 if __name__ == "__main__":
-    with petsird.BinaryPETSIRDReader(sys.stdin.buffer) as reader:
+    args = parserCreator()
+    file = None
+    if args.input is None:
+        file = sys.stdin.buffer
+    else:
+        file = open(args.input, "rb")
+    print_events = args.print_events
+
+    with petsird.BinaryPETSIRDReader(file) as reader:
         header = reader.read_header()
         if header.exam is not None:
             print(f"Subject ID: {header.exam.subject.id}")
@@ -42,28 +68,33 @@
               header.scanner.detection_efficiencies.module_pair_sgidlut)
         energy_1, energy_2 = 0.0, 0.0
         num_prompts = 0
+        num_delayeds = 0
         last_time = 0
         for time_block in reader.read_time_blocks():
             if isinstance(time_block, petsird.TimeBlock.EventTimeBlock):
                 last_time = time_block.value.start
                 num_prompts += len(time_block.value.prompt_events)
+                if time_block.value.delayed_events is not None:
+                    num_delayeds += len(time_block.value.delayed_events)
                 print("=====================  Events in time block from ",
                       last_time, " ==============")
                 for event in time_block.value.prompt_events:
                     energy_1 += energy_mid_points[event.energy_indices[0]]
                     energy_2 += energy_mid_points[event.energy_indices[1]]
-
-                    print(event)
-                    print(
-                        "   ",
-                        get_module_and_element(header.scanner.scanner_geometry,
-                                               event.detector_ids),
-                    )
-                    print("    efficiency:",
-                          get_detection_efficiency(header.scanner, event))
+                    if print_events:
+                        print(event)
+                        print(
+                            "   ",
+                            get_module_and_element(
+                                header.scanner.scanner_geometry,
+                                event.detector_ids),
+                        )
+                        print("    efficiency:",
+                              get_detection_efficiency(header.scanner, event))
 
         print(f"Last time block at {last_time} ms")
         print(f"Number of prompt events: {num_prompts}")
+        print(f"Number of delayed events: {num_delayeds}")
         if num_prompts > 0:
             print(f"Average energy_1: {energy_1 / num_prompts}")
             print(f"Average energy_2: {energy_2 / num_prompts}")