Implement MATLAB generator and analysis tools

.devcontainer/Dockerfile

@@ -87,7 +87,7 @@ RUN mkdir -p /home/vscode/.local/share/CMakeTools \
     && chown vscode:conda /home/vscode/.local/share/CMakeTools/cmake-tools-kits.json
 
 # Install the yardl tool
-ARG YARDL_VERSION=0.6.2
+ARG YARDL_VERSION=0.6.3
 RUN wget --quiet "https://github.com/microsoft/yardl/releases/download/v${YARDL_VERSION}/yardl_${YARDL_VERSION}_linux_x86_64.tar.gz" \
     && tar -xzf "yardl_${YARDL_VERSION}_linux_x86_64.tar.gz" \
     && mv yardl "/opt/conda/envs/${CONDA_ENVIRONMENT_NAME}/bin/" \

.github/workflows/ci.yml

@@ -8,7 +8,7 @@ on:
   pull_request: {branches: [main]}
 defaults: {run: {shell: 'bash -el {0}'}} # https://github.com/marketplace/actions/setup-miniconda#important
 env:
-  YARDL_VERSION: 0.6.2
+  YARDL_VERSION: 0.6.3
 jobs:
   validate:
     strategy:

.pre-commit-config.yaml

@@ -56,3 +56,4 @@ repos:
   rev: v14.0.6
   hooks:
   - id: clang-format
+    types_or: [c++]

matlab/.gitignore

@@ -0,0 +1 @@
+release/

matlab/build-toolbox.sh

@@ -0,0 +1,3 @@
+#!/usr/bin/env bash
+
+matlab -batch buildtool

matlab/buildfile.m

@@ -0,0 +1,41 @@
+function plan = buildfile
+    plan = buildplan(localfunctions);
+    plan.DefaultTasks = ["packageToolbox"];
+end
+
+function packageToolboxTask(~)
+    buildToolbox("release");
+end
+
+function buildToolbox (outdir)
+    if ~exist(outdir, 'dir')
+        mkdir(outdir);
+    end
+
+    uuid = string(java.util.UUID.randomUUID());
+    toolboxFolder = "./toolbox/";
+    opts = matlab.addons.toolbox.ToolboxOptions(toolboxFolder, uuid);
+
+    toolboxDirs = unique(fileparts(opts.ToolboxFiles));
+    if ~all(contains(toolboxDirs, opts.ToolboxFolder))
+        error("No symbolic links allowed in toolboxFolder (BUG through at least R2023a)");
+    end
+
+    opts.ToolboxName = "PETSIRD";
+
+    opts.ToolboxVersion = "0.2.0";
+    opts.OutputFile = fullfile(outdir, sprintf("petsird-%s.mltbx", opts.ToolboxVersion));
+
+    opts.Description = "Positron Emission Tomography Standardization Initiative Raw Data (PETSIRD) toolbox for MATLAB";
+    % opts.Summary = "";
+    % opts.AuthorCompany = "";
+
+    opts.MinimumMatlabRelease = "R2022b";
+
+    % Must also specify which folders should be added to MATLAB path upon toolbox installation.
+    % Must also include at least one *file *in the toolbox folder.
+    % This seems to be bug on Linux, Matlab R2023a. On Windows, this isn't required.
+    opts.ToolboxMatlabPath = toolboxFolder;
+
+    matlab.addons.toolbox.packageToolbox(opts);
+end

matlab/toolbox/+petsird/+helpers/get_detection_efficiency.m

@@ -0,0 +1,37 @@
+function eff = get_detection_efficiency(scanner, event)
+% Compute the detection efficiency for a coincidence event
+
+arguments
+    scanner (1,1) petsird.ScannerInformation
+    event (1,1) petsird.CoincidenceEvent
+end
+
+eff = 1;
+
+% Per det_el efficiencies
+det_el_efficiencies = scanner.detection_efficiencies.det_el_efficiencies;
+if det_el_efficiencies ~= yardl.None
+    eff = eff .* (det_el_efficiencies(event.energy_indices(1)+1, event.detector_ids(1)+1) .* ...
+            det_el_efficiencies(event.energy_indices(2)+1, event.detector_ids(2)+1));
+end
+
+% Per module-pair efficiencies
+module_pair_efficiencies_vector = scanner.detection_efficiencies.module_pair_efficiencies_vector;
+if module_pair_efficiencies_vector ~= yardl.None
+    module_pair_SGID_LUT = scanner.detection_efficiencies.module_pair_sgidlut;
+    assert(module_pair_SGID_LUT ~= yardl.None);
+    mod_and_els = petsird.helpers.get_module_and_element(scanner.scanner_geometry, event.detector_ids);
+    assert(length(scanner.scanner_geometry.replicated_modules) == 1);
+    SGID = module_pair_SGID_LUT(mod_and_els(1).module+1, mod_and_els(2).module+1);
+    assert(SGID >= 0);
+    module_pair_efficiencies = module_pair_efficiencies_vector(SGID+1);
+    assert(module_pair_efficiencies.sgid == SGID);
+    eff = eff * module_pair_efficiencies.values( ...
+        event.energy_indices(2)+1, ...
+        mod_and_els(2).el+1, ...
+        event.energy_indices(1)+1, ...
+        mod_and_els(1).el+1 ...
+    );
+end
+
+end

matlab/toolbox/+petsird/+helpers/get_module_and_element.m

@@ -0,0 +1,16 @@
+function mod_and_el = get_module_and_element(scanner_geometry, scanner_det_ids)
+% Find ModuleAndElement for a list of detector_ids
+
+arguments
+    scanner_geometry (1,1) petsird.ScannerGeometry
+    scanner_det_ids (:,1) uint32
+end
+
+assert(length(scanner_geometry.replicated_modules) == 1);
+rep_module = scanner_geometry.replicated_modules(1);
+assert(length(rep_module.object.detecting_elements) == 1);
+num_el_per_module = length(rep_module.object.detecting_elements(1).ids);
+
+mod_and_el = arrayfun(@(id) struct('module', idivide(id, num_el_per_module), 'el', mod(id, num_el_per_module)), scanner_det_ids);
+
+end

matlab/toolbox/+petsird/+helpers/get_num_detecting_elements.m

@@ -0,0 +1,18 @@
+function num_det_els = get_num_detecting_elements(scanner_geometry)
+% Compute total number of detecting elements in the scanner
+
+arguments
+    scanner_geometry (1,1) petsird.ScannerGeometry
+end
+
+num_det_els = 0;
+for rep_mod_idx = 1:length(scanner_geometry.replicated_modules)
+    rep_module = scanner_geometry.replicated_modules(rep_mod_idx);
+    det_els = rep_module.object.detecting_elements;
+    for rep_volume_idx = 1:length(det_els)
+        rep_volume = det_els(rep_volume_idx);
+        num_det_els = num_det_els + length(rep_volume.transforms) .* length(rep_module.transforms);
+    end
+end
+
+end

matlab/toolbox/petsird_analysis.m

@@ -0,0 +1,80 @@
+function petsird_generator(input)
+
+arguments
+    input (1,1) string
+end
+
+reader = petsird.binary.PETSIRDReader(input);
+
+header = reader.read_header();
+
+if header.exam ~= yardl.None
+    fprintf("Subject ID: %s\n", header.exam.subject.id);
+end
+
+fprintf("Scanner name: %s\n", header.scanner.model_name);
+
+fprintf("Types of modules: %d\n", length(header.scanner.scanner_geometry.replicated_modules));
+fprintf("Number of modules of first type: %d\n", ...
+    length(header.scanner.scanner_geometry.replicated_modules(1).transforms));
+fprintf("Number of types of detecting elements in modules of first type: %d\n", ...
+    length(header.scanner.scanner_geometry.replicated_modules(1).object.detecting_elements));
+fprintf("Number of elements of first type in modules of first type: %d\n", ...
+    length(header.scanner.scanner_geometry.replicated_modules(1).object.detecting_elements(1).transforms));
+fprintf("Total number of 'crystals': %d\n", ...
+    petsird.helpers.get_num_detecting_elements(header.scanner.scanner_geometry));
+fprintf("Number of TOF bins: %d\n", header.scanner.number_of_tof_bins());
+fprintf("Number of energy bins: %d\n", header.scanner.number_of_energy_bins());
+
+energy_bin_edges = header.scanner.energy_bin_edges;
+fprintf("Energy bin edges: " + join(repelem("%f", length(energy_bin_edges))) + "\n", energy_bin_edges);
+
+energy_mid_points = (energy_bin_edges(1:end-1) + energy_bin_edges(2:end)) / 2;
+fprintf("Energy mid points: " + join(repelem("%f", length(energy_mid_points))) + "\n", energy_mid_points);
+
+% sgidlut = header.scanner.detection_efficiencies.module_pair_sgidlut;
+% fprintf("SGID LUT:\n");
+% disp(sgidlut);
+
+energy_1 = 0.0;
+energy_2 = 0.0;
+num_prompts = 0;
+num_delayeds = 0;
+last_time = 0;
+while reader.has_time_blocks()
+    time_block = reader.read_time_blocks();
+
+    if time_block.isEventTimeBlock()
+        last_time = time_block.value.start;
+        num_prompts = num_prompts + length(time_block.value.prompt_events);
+        if time_block.value.delayed_events ~= yardl.None
+            num_delayeds = num_delayeds + length(time_block.value.delayed_events.value);
+        end
+        fprintf("=====================  Events in time block from %f  ==============\n", last_time);
+        for event_idx = 1:length(time_block.value.prompt_events)
+            event = time_block.value.prompt_events(event_idx);
+            energy_1 = energy_1 + energy_mid_points(event.energy_indices(1) + 1);
+            energy_2 = energy_2 + energy_mid_points(event.energy_indices(2) + 1);
+
+            fprintf("CoincidenceEvent(detectorIds=[%d, %d], tofIdx=%d, energyIndices=[%d, %d])\n", ...
+                event.detector_ids(1), event.detector_ids(2), event.tof_idx, ...
+                event.energy_indices(1), event.energy_indices(2));
+            mod_and_el = petsird.helpers.get_module_and_element(header.scanner.scanner_geometry, event.detector_ids);
+            fprintf("    [ModuleAndElement(module=%d, el=%d), ModuleAndElement(module=%d, el=%d)]\n", ...
+                mod_and_el(1).module, mod_and_el(1).el, mod_and_el(2).module, mod_and_el(2).el);
+            fprintf("    efficiency: %f\n", petsird.helpers.get_detection_efficiency(header.scanner, event));
+        end
+    end
+end
+
+fprintf("Last time block at %f ms\n", last_time);
+fprintf("Number of prompt events: %d\n", num_prompts);
+fprintf("Number of delayed events: %d\n", num_delayeds);
+if num_prompts > 0
+    fprintf("Average energy_1: %f\n", energy_1 / num_prompts)
+    fprintf("Average energy_2: %f\n", energy_2 / num_prompts)
+end
+
+reader.close();
+
+end