From cd8544b7f9e7915f2ce22c6beae7d17d9e13f188 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Sun, 3 Nov 2024 19:11:30 +0000 Subject: [PATCH 01/13] update submodule --- .gitmodules | 3 ++- PETSIRD | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/.gitmodules b/.gitmodules index c1703ce..e392867 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,3 +1,4 @@ [submodule "PETSIRD"] path = PETSIRD - url = https://github.com/ETSInitiative/PRDdefinition + url = https://github.com/ETSInitiative/PETSIRD + branch=main diff --git a/PETSIRD b/PETSIRD index f168903..ae3a048 160000 --- a/PETSIRD +++ b/PETSIRD @@ -1 +1 @@ -Subproject commit f16890385facda7ac949d8ab622e4e8fee3332d2 +Subproject commit ae3a0484a128f805a2143cc0d010d01557cf902f From 08ef5aa453d932c222f2dcc7a93eab7ad9c5859b Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Sun, 3 Nov 2024 19:13:44 +0000 Subject: [PATCH 02/13] update devcontainer --- .devcontainer/Dockerfile | 22 +++++++++------------- .devcontainer/devcontainer.bashrc | 3 +-- .devcontainer/devcontainer.json | 21 +++++++++------------ 3 files changed, 19 insertions(+), 27 deletions(-) diff --git a/.devcontainer/Dockerfile b/.devcontainer/Dockerfile index 4e0fbee..8677bde 100644 --- a/.devcontainer/Dockerfile +++ b/.devcontainer/Dockerfile @@ -29,18 +29,16 @@ ARG USERNAME="vscode" ARG USER_UID=1000 ARG USER_GID=$USER_UID ARG CONDA_GID=900 -ARG CONDA_ENVIRONMENT_NAME=prd -ARG VSCODE_DEV_CONTAINERS_SCRIPT_LIBRARY_VERSION=v0.229.0 +ARG CONDA_ENVIRONMENT_NAME=petsird RUN apt-get update && apt-get install -y \ libc6-dbg \ && rm -rf /var/lib/apt/lists/* # Enable non-root Docker access in container -ARG ENABLE_NONROOT_DOCKER="true" -# Use the OSS Moby CLI instead of the licensed Docker CLI -ARG USE_MOBY="false" -RUN script=$(curl -fsSL "https://raw.githubusercontent.com/microsoft/vscode-dev-containers/${VSCODE_DEV_CONTAINERS_SCRIPT_LIBRARY_VERSION}/script-library/docker-debian.sh") && bash -c "$script" -- "${ENABLE_NONROOT_DOCKER}" "/var/run/docker-host.sock" "/var/run/docker.sock" "${USERNAME}" "${USE_MOBY}" +ARG DOCKER_VERSION="27.0.3" +RUN script=$(curl -fsSL "https://raw.githubusercontent.com/devcontainers/features/2951f0481a488ea43a6f2ea6f18a47f0a0bf744d/src/docker-outside-of-docker/install.sh") \ + && VERSION=${DOCKER_VERSION} DOCKERDASHCOMPOSEVERSION="none" bash -c "$script" # Setting the ENTRYPOINT to docker-init.sh will configure non-root access to # the Docker socket if "overrideCommand": false is set in devcontainer.json. @@ -48,10 +46,8 @@ RUN script=$(curl -fsSL "https://raw.githubusercontent.com/microsoft/vscode-dev- ENTRYPOINT [ "/usr/local/share/docker-init.sh" ] CMD [ "sleep", "infinity" ] -ARG MAMBAFORGE_VERSION=22.9.0-2 - # Based on https://github.com/conda-forge/miniforge-images/blob/master/ubuntu/Dockerfile -RUN wget --no-hsts --quiet https://github.com/conda-forge/miniforge/releases/download/${MAMBAFORGE_VERSION}/Mambaforge-${MAMBAFORGE_VERSION}-Linux-$(uname -m).sh -O /tmp/miniforge.sh \ +RUN wget --no-hsts --quiet https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-$(uname)-$(uname -m).sh -O /tmp/miniforge.sh \ && /bin/bash /tmp/miniforge.sh -b -p /opt/conda \ && rm /tmp/miniforge.sh \ && /opt/conda/bin/conda clean --tarballs --index-cache --packages --yes \ @@ -67,8 +63,8 @@ RUN wget --no-hsts --quiet https://github.com/conda-forge/miniforge/releases/dow # Create a conda environment from the environment file in the repo root. COPY --from=file-normalizer --chown=$USER_UID:conda /data/environment.yml /tmp/build/ RUN umask 0002 \ - && /opt/conda/bin/mamba env create -f /tmp/build/environment.yml \ - && /opt/conda/bin/mamba clean -fy \ + && /opt/conda/bin/conda env create -f /tmp/build/environment.yml \ + && /opt/conda/bin/conda clean -fy \ && sudo chown -R :conda /opt/conda/envs # Add a file that is to be sourced from .bashrc and from the devops pipeline stages @@ -87,11 +83,11 @@ ENV CMAKE_GENERATOR=Ninja # Create a kits file for the VSCode CMake Tools extension, so you are not prompted for which kit to select whenever you open VSCode RUN mkdir -p /home/vscode/.local/share/CMakeTools \ - && echo '[{"name":"GCC-10","compilers":{"C":"/opt/conda/envs/prd/bin/x86_64-conda_cos6-linux-gnu-gcc","CXX":"/opt/conda/envs/prd/bin/x86_64-conda_cos6-linux-gnu-g++"}}]' > /home/vscode/.local/share/CMakeTools/cmake-tools-kits.json \ + && echo '[{"name":"GCC-10","compilers":{"C":"/opt/conda/envs/petsird/bin/x86_64-conda_cos6-linux-gnu-gcc","CXX":"/opt/conda/envs/petsird/bin/x86_64-conda_cos6-linux-gnu-g++"}}]' > /home/vscode/.local/share/CMakeTools/cmake-tools-kits.json \ && chown vscode:conda /home/vscode/.local/share/CMakeTools/cmake-tools-kits.json # Install the yardl tool -ARG YARDL_VERSION=0.3.2 +ARG YARDL_VERSION=0.6.2 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/" \ diff --git a/.devcontainer/devcontainer.bashrc b/.devcontainer/devcontainer.bashrc index 2812baa..aa9919c 100644 --- a/.devcontainer/devcontainer.bashrc +++ b/.devcontainer/devcontainer.bashrc @@ -2,9 +2,8 @@ # shellcheck source=/dev/null source /opt/conda/etc/profile.d/conda.sh -conda activate prd +conda activate petsird source <(yardl completion bash) -source <(just --completions bash) if [[ "${BASH_ENV:-}" == "$(readlink -f "${BASH_SOURCE[0]:-}")" ]]; then # We don't want subshells to unnecessarily source this again. diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json index cc2ee7d..ed472ee 100644 --- a/.devcontainer/devcontainer.json +++ b/.devcontainer/devcontainer.json @@ -1,7 +1,6 @@ -// For format details, see https://aka.ms/devcontainer.json. For config options, see the README at: -// https://github.com/microsoft/vscode-dev-containerdevcotns/tree/v0.238.0/containers/go +// For format details, see https://aka.ms/devcontainer.json. { - "name": "prd", + "name": "petsird", "build": { "dockerfile": "Dockerfile", "context": ".." @@ -48,7 +47,7 @@ "cmake.buildDirectory": "${workspaceFolder}/cpp/build", "cmake.configureOnOpen": false, - "python.defaultInterpreterPath": "/opt/conda/envs/prd/bin/python", + "python.defaultInterpreterPath": "/opt/conda/envs/petsird/bin/python", "python.analysis.typeCheckingMode": "strict", "python.analysis.diagnosticMode": "workspace", "python.analysis.diagnosticSeverityOverrides": { @@ -88,7 +87,7 @@ "gitlens.showWhatsNewAfterUpgrades": false }, - "gcovViewer.gcovBinary": "/opt/conda/envs/prd/bin/x86_64-conda-linux-gnu-gcov", + "gcovViewer.gcovBinary": "/opt/conda/envs/petsird/bin/x86_64-conda-linux-gnu-gcov", "gcovViewer.buildDirectories": ["${workspaceFolder}/cpp/build"], "search.useIgnoreFiles": false, @@ -109,19 +108,17 @@ "ms-python.python", "ms-vscode.cmake-tools", "ms-vscode.cpptools", - "sclu1034.justfile", "timonwong.shellcheck", "twxs.cmake" ] } }, - // Use 'forwardPorts' to make a list of ports inside the container available locally. - // "forwardPorts": [], - - // Use 'postCreateCommand' to run commands after the container is created. - // "postCreateCommand": "go version", + "overrideCommand": false, + "mounts": [ + // Bind mount docker socket under an alias to support docker-from-docker + "source=/var/run/docker.sock,target=/var/run/docker-host.sock,type=bind" + ], - // Comment out to connect as root instead. More info: https://aka.ms/vscode-remote/containers/non-root. "remoteUser": "vscode" } From 1152d9450d964c79f2166ad4f49e2681ae2c9bcf Mon Sep 17 00:00:00 2001 From: nkarakatsanis Date: Sun, 3 Nov 2024 22:50:14 +0000 Subject: [PATCH 03/13] Intermediate changes to the converter --- cpp/main.cpp | 153 ++++++++++++++++++++++++++------ data/root/scanner_geometry.json | 1 + environment.yml | 2 +- 3 files changed, 127 insertions(+), 29 deletions(-) diff --git a/cpp/main.cpp b/cpp/main.cpp index bad19d2..75e5adf 100644 --- a/cpp/main.cpp +++ b/cpp/main.cpp @@ -101,6 +101,21 @@ using namespace std ; +#include "petsird_helpers.h" + +// these are constants for now +constexpr uint32_t NUMBER_OF_ENERGY_BINS = 3; +constexpr uint32_t NUMBER_OF_TOF_BINS = 300; +constexpr float RADIUS = 400.F; +constexpr std::array CRYSTAL_LENGTH{ 20.F, 4.F, 4.F }; +constexpr std::array NUM_CRYSTALS_PER_MODULE{ 2, 4, 5 }; +constexpr uint32_t NUM_MODULES_ALONG_RING{ 20 }; +constexpr uint32_t NUM_MODULES_ALONG_AXIS{ 2 }; +constexpr float MODULE_AXIS_SPACING{ (NUM_CRYSTALS_PER_MODULE[2] + 4) * CRYSTAL_LENGTH[2] }; + +constexpr uint32_t NUMBER_OF_TIME_BLOCKS = 6; +constexpr float COUNT_RATE = 500.F; + struct ScannerGeometry { int n_rings; @@ -118,6 +133,7 @@ struct ScannerGeometry int n_crystal; int n_cry_xy; int n_cry_z; + int n_cry_layer; int max_d_ring; int number_of_tof_bins; int number_of_energy_bins; @@ -202,6 +218,7 @@ ScannerGeometry ReadScannerGeometry(const std::string& filename) scanner_geometry.n_crystal = j["n_crystal"]; scanner_geometry.n_cry_xy = j["n_cry_xy"]; scanner_geometry.n_cry_z = j["n_cry_z"]; + scanner_geometry.n_cry_layer = j["n_cry_layer"]; scanner_geometry.max_d_ring = j["max_d_ring"]; scanner_geometry.number_of_tof_bins = j["number_of_tof_bins"]; scanner_geometry.number_of_energy_bins = j["number_of_energy_bins"]; @@ -250,6 +267,7 @@ int calculate_detector_id(int gantry_id, int rsector_id, int module_id, int subm int N_SMOD_z = scannerGeometry.n_smod_z; int N_CRY_xy = scannerGeometry.n_cry_xy; int N_CRY_z = scannerGeometry.n_cry_z; + int N_CRY_layer = scannerGeometry.n_cry_layer; int ring = (Int_t)(gantry_id)*N_RSEC_z*N_MOD_z*N_SMOD_z*N_CRY_z + (Int_t)(rsector_id/N_RSEC_xy)*N_MOD_z*N_SMOD_z*N_CRY_z @@ -266,7 +284,7 @@ int calculate_detector_id(int gantry_id, int rsector_id, int module_id, int subm } // single ring as example -prd::ScannerInformation +petsird::ScannerInformation get_scanner_info(ScannerGeometry& scannerGeometry) { float radius = scannerGeometry.radius; @@ -283,14 +301,92 @@ get_scanner_info(ScannerGeometry& scannerGeometry) angles.push_back(static_cast(2 * M_PI * (1.0f*i) / n_detectors)); } - std::vector detectors; +//! return a cuboid volume +petsird::BoxSolidVolume +get_crystal() +{ + using petsird::Coordinate; + petsird::BoxShape crystal_shape{ Coordinate{ { 0, 0, 0 } }, + Coordinate{ { 0, 0, detector_z_dim } }, + Coordinate{ { 0, detector_y_dim, detector_z_dim } }, + Coordinate{ { 0, detector_y_dim, 0 } }, + Coordinate{ { detector_x_dim, 0, 0 } }, + Coordinate{ { detector_x_dim, 0, detector_z_dim } }, + Coordinate{ { detector_x_dim, detector_y_dim, detector_z_dim } }, + Coordinate{ { detector_x_dim, detector_y_dim, 0 } } }; + + petsird::BoxSolidVolume crystal{ crystal_shape, /* material_id */ 1 }; + return crystal; +} + +//! return a module of NUM_CRYSTALS_PER_MODULE cuboids +petsird::DetectorModule +get_detector_module() +{ + petsird::ReplicatedBoxSolidVolume rep_volume; + { + rep_volume.object = get_crystal(); + for (int rep_mod_xy = 0; rep_mod_xy < n_mod_xy; ++rep_mod_xy) + for (int rep_mod_z = 0; rep_mod_z < n_mod_z; ++rep_mod_z) + for (int rep_cry_layer = 0; rep_cry_layer < n_cry_layer; ++rep_cry_layer) + for (int rep_cry_xy = 0; rep_cry_xy < n_cry_xy; ++rep_cry_xy) + for (int rep_cry_z = 0; rep_cry_z < n_cry_z; ++rep_cry_z) + { + petsird::RigidTransformation transform{ { { 1.0, 0.0, 0.0, radius + rep_cry_layer * detector_x_dim }, + { 0.0, 1.0, 0.0, (rep_modxy - n_mod_xy / 2) * (rep_cry_xy - n_cry_xy / 2) * detector_y_dim }, + { 0.0, 0.0, 1.0, (rep_mod_z - n_mod_z / 2) *(rep_cry_z - n_cry_z / 2) * detector_z_dim } } }; + rep_volume.transforms.push_back(transform); + rep_volume.ids.push_back(rep_cry_z + n_cry_z * (rep_cry_xy + n_cry_xy * rep_cry_layer)); + } + } + + petsird::DetectorModule detector_module; + detector_module.detecting_elements.push_back(rep_volume); + detector_module.detecting_element_ids.push_back(0); + + return detector_module; +} + + + +//! return scanner build by rotating a module around the (0,0,1) axis +petsird::ScannerGeometry +get_scanner_geometry() +{ + petsird::ReplicatedDetectorModule rep_module; + { + rep_module.object = get_detector_module(); + int module_id = 0; + std::vector angles; + for (unsigned int i = 0; i < NUM_MODULES_ALONG_RING; ++i) + { + angles.push_back(static_cast((2 * M_PI * i) / NUM_MODULES_ALONG_RING)); + } + for (auto angle : angles) + for (unsigned ax_mod = 0; ax_mod < NUM_MODULES_ALONG_AXIS; ++ax_mod) + { + petsird::RigidTransformation transform{ { { std::cos(angle), std::sin(angle), 0.F, 0.F }, + { -std::sin(angle), std::cos(angle), 0.F, 0.F }, + { 0.F, 0.F, 1.F, MODULE_AXIS_SPACING * ax_mod } } }; + rep_module.ids.push_back(module_id++); + rep_module.transforms.push_back(transform); + } + } + petsird::ScannerGeometry scanner_geometry; + scanner_geometry.replicated_modules.push_back(rep_module); + scanner_geometry.ids.push_back(0); + return scanner_geometry; +} + + + std::vector detectors; int detector_id = 0; for (int r =0; r < n_rings; r++) { for (auto angle : angles) { // Create a new detector - prd::Detector d; + petsird::Detector d; d.x = radius * std::cos(angle); d.y = radius * std::sin(angle); d.z = ((-n_rings/2.0f)*scannerGeometry.detector_z_dim) + scannerGeometry.detector_z_dim*r; @@ -311,8 +407,9 @@ get_scanner_info(ScannerGeometry& scannerGeometry) for (std::size_t i = 0; i < energy_bin_edges.size(); ++i) { energy_bin_edges[i] = energy_LLD + i * (energy_ULD - energy_LLD) / NUMBER_OF_ENERGY_BINS; } - prd::ScannerInformation scanner_info; - scanner_info.detectors = detectors; + petsird::ScannerInformation scanner_info; + scanner_info.scanner_geometry = get_scanner_geometry(); + //scanner_info.detectors = detectors; scanner_info.tof_bin_edges = tof_bin_edges; scanner_info.tof_resolution = scannerGeometry.TOF_resolution*0.3; // conversion from psec to mm (e.g. 200ps TOF is equivalent to 60mm uncertainty) scanner_info.energy_bin_edges = energy_bin_edges; @@ -321,7 +418,7 @@ get_scanner_info(ScannerGeometry& scannerGeometry) return scanner_info; } -uint32_t tofToIdx(double delta_time_psec, const prd::ScannerInformation& scanner_info) +uint32_t tofToIdx(double delta_time_psec, const petsird::ScannerInformation& scanner_info) { float tofPos_mm = delta_time_psec * 0.15; //conversion from time difference (in psec) to spatial position in LOR (in mm) DT*C/2 for (size_t i = 0; i < scanner_info.tof_bin_edges.size() - 1; ++i) @@ -337,7 +434,7 @@ uint32_t tofToIdx(double delta_time_psec, const prd::ScannerInformation& scanner throw std::runtime_error("TOF out of range"); } -uint32_t energyToIdx(float energy, const prd::ScannerInformation& scanner_info) +uint32_t energyToIdx(float energy, const petsird::ScannerInformation& scanner_info) { for (size_t i = 0; i < scanner_info.energy_bin_edges.size() - 1; ++i) { @@ -492,8 +589,8 @@ int main(int argc, char** argv) printf("Total Number of Coincidence Events in the ROOT file:= %llu \n",nentries ); // Output PETSIRD - prd::Header header; - prd::ScannerInformation scanner = get_scanner_info(scannerGeometry); + petsird::Header header; + petsird::ScannerInformation scanner = get_scanner_info(scannerGeometry); if (verbose) { // Print scanner information @@ -510,17 +607,17 @@ int main(int argc, char** argv) //} } - prd::ExamInformation exam; + petsird::ExamInformation exam; header.exam = exam; header.scanner = scanner; // Write PETSiRD file - prd::binary::PrdExperimentWriter writer(petsird_file); + petsird::binary::petsirdExperimentWriter writer(petsird_file); writer.WriteHeader(header); long current_time_block = -1; - prd::TimeBlock time_block; + petsird::EventTimeBlock time_block; unsigned long Counts_binned = 0; for (unsigned long long int i = 0 ; i < nentries ; i++) { @@ -532,41 +629,41 @@ int main(int argc, char** argv) if (eventID1 == eventID2) { if (comptonPhantom1 == 0 && comptonPhantom2 == 0) { - prd::CoincidenceEvent event; - event.detector_1_id = calculate_detector_id(gantryID1, rsectorID1, moduleID1, submoduleID1, crystalID1, scannerGeometry); - event.detector_2_id = calculate_detector_id(gantryID2, rsectorID2, moduleID2, submoduleID2, crystalID2, scannerGeometry); + petsird::CoincidenceEvent event; + event.detector_ids[0] = calculate_detector_id(gantryID1, rsectorID1, moduleID1, submoduleID1, crystalID1, scannerGeometry); + event.detector_ids[1] = calculate_detector_id(gantryID2, rsectorID2, moduleID2, submoduleID2, crystalID2, scannerGeometry); double dt_psec = 1.0e12f*(time1 - time2); //in psec if (abs(dt_psec) > scannerGeometry.TxFOV_TOF/0.3f) { continue; } event.tof_idx = static_cast(tofToIdx(dt_psec, scanner)); - event.energy_1_idx = static_cast(energyToIdx(1.0e3*energy1, scanner)); - event.energy_2_idx = static_cast(energyToIdx(1.0e3*energy2, scanner)); + event.energy_indices[0] = static_cast(energyToIdx(1.0e3*energy1, scanner)); + event.energy_indices[1] = static_cast(energyToIdx(1.0e3*energy2, scanner)); if (verbose && i%100000 == 0) { std::cout << "Event " << i << std::endl; - std::cout << " detector_1_id: " << event.detector_1_id << std::endl; - std::cout << " detector_2_id: " << event.detector_2_id << std::endl; + std::cout << " detector_1_id: " << event.detector_ids[0] << std::endl; + std::cout << " detector_2_id: " << event.detector_ids[1] << std::endl; std::cout << " tof_idx: " << event.tof_idx << std::endl; - std::cout << " energy_1_idx: " << event.energy_1_idx << std::endl; - std::cout << " energy_2_idx: " << event.energy_2_idx << std::endl; - std::cout << " detector 1 position: " << scanner.detectors[event.detector_1_id].x << ", " << scanner.detectors[event.detector_1_id].y << ", " << scanner.detectors[event.detector_1_id].z << std::endl; + std::cout << " energy_1_idx: " << event.energy_indices[0] << std::endl; + std::cout << " energy_2_idx: " << event.energy_indices[1] << std::endl; + std::cout << " detector 1 position: " << scanner.detectors[event.detector_ids[0]].x << ", " << scanner.detectors[event.detector_ids[0]].y << ", " << scanner.detectors[event.detector_ids[0]].z << std::endl; std::cout << " GlobalPosition 1: " << globalPosX1 << ", " << globalPosY1 << ", " << globalPosZ1 << std::endl; - float distance_1 = std::sqrt(std::pow(scanner.detectors[event.detector_1_id].x-globalPosX1, 2) + std::pow(scanner.detectors[event.detector_1_id].y-globalPosY1, 2) + std::pow(scanner.detectors[event.detector_1_id].z-globalPosZ1, 2)); + float distance_1 = std::sqrt(std::pow(scanner.detectors[event.detector_ids[0]].x-globalPosX1, 2) + std::pow(scanner.detectors[event.detector_ids[0]].y-globalPosY1, 2) + std::pow(scanner.detectors[event.detector_ids[0]].z-globalPosZ1, 2)); std::cout << " Distance 1: " << distance_1 << std::endl; - std::cout << " detector 2 position: " << scanner.detectors[event.detector_2_id].x << ", " << scanner.detectors[event.detector_2_id].y << ", " << scanner.detectors[event.detector_2_id].z << std::endl; + std::cout << " detector 2 position: " << scanner.detectors[event.detector_ids[1]].x << ", " << scanner.detectors[event.detector_ids[1]].y << ", " << scanner.detectors[event.detector_ids[1]].z << std::endl; std::cout << " GlobalPosition 2: " << globalPosX2 << ", " << globalPosY2 << ", " << globalPosZ2 << std::endl; - float distance_2 = std::sqrt(std::pow(scanner.detectors[event.detector_2_id].x-globalPosX2, 2) + std::pow(scanner.detectors[event.detector_2_id].y-globalPosY2, 2) + std::pow(scanner.detectors[event.detector_2_id].z-globalPosZ2, 2)); + float distance_2 = std::sqrt(std::pow(scanner.detectors[event.detector_ids[1]].x-globalPosX2, 2) + std::pow(scanner.detectors[event.detector_ids[1]].y-globalPosY2, 2) + std::pow(scanner.detectors[event.detector_ids[1]].z-globalPosZ2, 2)); std::cout << " Distance 2: " << distance_2 << std::endl; } - long this_time_block = static_cast(time1*1.0e3 / scanner.listmode_time_block_duration); + long this_time_block = static_cast(time1*1.0e3 / scanner.event_time_block_duration); if (this_time_block != current_time_block) { if (current_time_block != -1) { writer.WriteTimeBlocks(time_block); } current_time_block = this_time_block; - time_block = prd::TimeBlock(); - time_block.id = static_cast(current_time_block); + time_block = petsird::TimeBlock(); + time_block.start = time1*1.0e3; } time_block.prompt_events.push_back(event); Counts_binned++; diff --git a/data/root/scanner_geometry.json b/data/root/scanner_geometry.json index 80aa662..9b8aa61 100644 --- a/data/root/scanner_geometry.json +++ b/data/root/scanner_geometry.json @@ -2,6 +2,7 @@ "ax_phys_crystal_num": 5, "ax_virtual_crystal_num": 0, "max_d_ring": 39, + "n_cry_layer":1, "n_cry_xy": 5, "n_cry_z": 5, "n_crystal": 25, diff --git a/environment.yml b/environment.yml index 78bda2e..d4471a7 100644 --- a/environment.yml +++ b/environment.yml @@ -1,4 +1,4 @@ -name: prd +name: petsird channels: - conda-forge - defaults From 67db64421eed6f7f22a0029421acd122f3f13f27 Mon Sep 17 00:00:00 2001 From: nkarakatsanis Date: Mon, 4 Nov 2024 04:57:19 +0000 Subject: [PATCH 04/13] Added geometry info --- cpp/main.cpp | 100 +++++++++++++++++++++----------- data/root/scanner_geometry.json | 10 +++- 2 files changed, 76 insertions(+), 34 deletions(-) diff --git a/cpp/main.cpp b/cpp/main.cpp index 75e5adf..06b2f4f 100644 --- a/cpp/main.cpp +++ b/cpp/main.cpp @@ -133,7 +133,15 @@ struct ScannerGeometry int n_crystal; int n_cry_xy; int n_cry_z; - int n_cry_layer; + int n_cry_layers; + int cry_ax_gap; + int cry_tx_gap; + int smod_ax_gap; + int smod_tx_gap; + int mod_ax_gap; + int mod_tx_gap; + int rsec_ax_gap; + int rsec_tx_gap; int max_d_ring; int number_of_tof_bins; int number_of_energy_bins; @@ -170,6 +178,15 @@ void WriteScannerGeometry(const ScannerGeometry& scanner_geometry, const std::st j["n_crystal"] = scanner_geometry.n_crystal; j["n_cry_xy"] = scanner_geometry.n_cry_xy; j["n_cry_z"] = scanner_geometry.n_cry_z; + j["n_cry_layers"] = scanner_geometry.n_cry_layers; + j["cry_ax_gap"] = scanner_geometry.cry_ax_gap; + j["cry_tx_gap"] = scanner_geometry.cry_tx_gap; + j["smod_ax_gap"] = scanner_geometry.smod_ax_gap; + j["smod_tx_gap"] = scanner_geometry.smod_tx_gap; + j["mod_ax_gap"] = scanner_geometry.mod_ax_gap; + j["mod_tx_gap"] = scanner_geometry.mod_tx_gap; + j["rsec_ax_gap"] = scanner_geometry.rsec_ax_gap; + j["rsec_tx_gap"] = scanner_geometry.rsec_tx_gap; j["max_d_ring"] = scanner_geometry.max_d_ring; j["number_of_tof_bins"] = scanner_geometry.number_of_tof_bins; j["number_of_energy_bins"] = scanner_geometry.number_of_energy_bins; @@ -218,7 +235,15 @@ ScannerGeometry ReadScannerGeometry(const std::string& filename) scanner_geometry.n_crystal = j["n_crystal"]; scanner_geometry.n_cry_xy = j["n_cry_xy"]; scanner_geometry.n_cry_z = j["n_cry_z"]; - scanner_geometry.n_cry_layer = j["n_cry_layer"]; + scanner_geometry.n_cry_layers = j["n_cry_layers"]; + scanner_geometry.cry_ax_gap = j["cry_ax_gap"]; + scanner_geometry.cry_tx_gap = j["cry_tx_gap"]; + scanner_geometry.smod_ax_gap = j["smod_ax_gap"]; + scanner_geometry.smod_tx_gap = j["smod_tx_gap"]; + scanner_geometry.mod_ax_gap = j["mod_ax_gap"]; + scanner_geometry.mod_tx_gap = j["mod_tx_gap"]; + scanner_geometry.rsec_ax_gap = j["rsec_ax_gap"]; + scanner_geometry.rsec_tx_gap = j["rsec_tx_gap"]; scanner_geometry.max_d_ring = j["max_d_ring"]; scanner_geometry.number_of_tof_bins = j["number_of_tof_bins"]; scanner_geometry.number_of_energy_bins = j["number_of_energy_bins"]; @@ -239,6 +264,7 @@ ScannerGeometry ReadScannerGeometry(const std::string& filename) scanner_geometry.ArcLength = scanner_geometry.s_width * scanner_geometry.detector_y_dim / 2.0f; scanner_geometry.TxFOV = 2 * scanner_geometry.radius * sin (scanner_geometry.ArcLength / (2 * scanner_geometry.radius) ); scanner_geometry.TxFOV_TOF = scanner_geometry.TxFOV + 0.3*scanner_geometry.TOF_resolution; + scanner_geometry.module_axial_pitch = n_cry_z * detector_z_dim + (n_cry_z - 1) * cry_ax_gap; return scanner_geometry; } @@ -267,18 +293,18 @@ int calculate_detector_id(int gantry_id, int rsector_id, int module_id, int subm int N_SMOD_z = scannerGeometry.n_smod_z; int N_CRY_xy = scannerGeometry.n_cry_xy; int N_CRY_z = scannerGeometry.n_cry_z; - int N_CRY_layer = scannerGeometry.n_cry_layer; + int N_CRY_layers = scannerGeometry.n_cry_layers; int ring = (Int_t)(gantry_id)*N_RSEC_z*N_MOD_z*N_SMOD_z*N_CRY_z - + (Int_t)(rsector_id/N_RSEC_xy)*N_MOD_z*N_SMOD_z*N_CRY_z - + (Int_t)(module_id/N_MOD_xy)*N_SMOD_z*N_CRY_z - + (Int_t)(submodule_id/N_SMOD_xy)*N_CRY_z - + (Int_t)(crystal_id/N_CRY_xy); + + (Int_t)(rsector_id/N_RSEC_xy)*N_MOD_z*N_SMOD_z*N_CRY_z + + (Int_t)(module_id/N_MOD_xy)*N_SMOD_z*N_CRY_z + + (Int_t)(submodule_id/N_SMOD_xy)*N_CRY_z + + (Int_t)(crystal_id/N_CRY_xy); int crystal = (Int_t)(crystal_id%N_CRY_xy) - + (Int_t)(submodule_id%N_SMOD_xy)*N_CRY_xy - + (Int_t)(module_id%N_MOD_xy)*N_SMOD_xy*N_CRY_xy - + (Int_t)(rsector_id%N_RSEC_xy)*N_MOD_xy*N_SMOD_xy*N_CRY_xy; + + (Int_t)(submodule_id%N_SMOD_xy)*N_CRY_xy + + (Int_t)(module_id%N_MOD_xy)*N_SMOD_xy*N_CRY_xy + + (Int_t)(rsector_id%N_RSEC_xy)*N_MOD_xy*N_SMOD_xy*N_CRY_xy; return crystal + (ring-rmin)*N_DET; } @@ -326,18 +352,17 @@ get_detector_module() petsird::ReplicatedBoxSolidVolume rep_volume; { rep_volume.object = get_crystal(); - for (int rep_mod_xy = 0; rep_mod_xy < n_mod_xy; ++rep_mod_xy) - for (int rep_mod_z = 0; rep_mod_z < n_mod_z; ++rep_mod_z) - for (int rep_cry_layer = 0; rep_cry_layer < n_cry_layer; ++rep_cry_layer) - for (int rep_cry_xy = 0; rep_cry_xy < n_cry_xy; ++rep_cry_xy) - for (int rep_cry_z = 0; rep_cry_z < n_cry_z; ++rep_cry_z) - { - petsird::RigidTransformation transform{ { { 1.0, 0.0, 0.0, radius + rep_cry_layer * detector_x_dim }, - { 0.0, 1.0, 0.0, (rep_modxy - n_mod_xy / 2) * (rep_cry_xy - n_cry_xy / 2) * detector_y_dim }, - { 0.0, 0.0, 1.0, (rep_mod_z - n_mod_z / 2) *(rep_cry_z - n_cry_z / 2) * detector_z_dim } } }; - rep_volume.transforms.push_back(transform); - rep_volume.ids.push_back(rep_cry_z + n_cry_z * (rep_cry_xy + n_cry_xy * rep_cry_layer)); - } + for (int rep_cry_layer = 0; rep_cry_layer < n_cry_layers; ++rep_cry_layer) + for (int rep_cry_xy = 0; rep_cry_xy < n_cry_xy; ++rep_cry_xy) + for (int rep_cry_z = 0; rep_cry_z < n_cry_z; ++rep_cry_z) + { + petsird::RigidTransformation transform{ { { 1.0, 0.0, 0.0, radius + rep_cry_layer * detector_x_dim }, + { 0.0, 1.0, 0.0, (rep_cry_xy - n_cry_xy / 2) * detector_y_dim }, + { 0.0, 0.0, 1.0, (rep_cry_z - n_cry_z / 2) * detector_z_dim } } }; + rep_volume.transforms.push_back(transform); + //rep_volume.ids.push_back(rep_cry_layer + n_cry_layer * (rep_cry_xy + n_cry_xy * rep_cry_z)); + rep_volume.ids.push_back(rep_cry_z + n_cry_z * (rep_cry_xy + n_cry_xy * rep_cry_layer)); + } } petsird::DetectorModule detector_module; @@ -358,19 +383,28 @@ get_scanner_geometry() rep_module.object = get_detector_module(); int module_id = 0; std::vector angles; - for (unsigned int i = 0; i < NUM_MODULES_ALONG_RING; ++i) + for (unsigned int i = 0; i < n_rsec_xy; ++i) { - angles.push_back(static_cast((2 * M_PI * i) / NUM_MODULES_ALONG_RING)); + angles.push_back(static_cast((2 * M_PI * i) / n_rsec_xy)); } - for (auto angle : angles) - for (unsigned ax_mod = 0; ax_mod < NUM_MODULES_ALONG_AXIS; ++ax_mod) - { - petsird::RigidTransformation transform{ { { std::cos(angle), std::sin(angle), 0.F, 0.F }, - { -std::sin(angle), std::cos(angle), 0.F, 0.F }, - { 0.F, 0.F, 1.F, MODULE_AXIS_SPACING * ax_mod } } }; - rep_module.ids.push_back(module_id++); - rep_module.transforms.push_back(transform); - } + + for (int rep_smod_xy = 0; rep_smod_xy < n_smod_xy; ++rep_smod_xy) + for (int rep_smod_z = 0; rep_smod_z < n_smod_z; ++rep_smod_z) + for (int rep_mod_xy = 0; rep_mod_xy < n_mod_xy; ++rep_mod_xy) + for (int rep_mod_z = 0; rep_mod_z < n_mod_z; ++rep_mod_z) + for (auto angle : angles) + for (int rep_rsec_z = 0; rep_rsec_z < n_rsec_z; ++rep_rsec_z) + { + petsird::RigidTransformation transform{ { { std::cos(angle), std::sin(angle), 0.F, 0.F }, + { -std::sin(angle), std::cos(angle), 0.F, (rep_smod_xy - n_smod_xy / 2) * n_smod_xy * detector_y_dim + + (rep_mod_xy - n_mod_xy / 2) * n_mod_xy * n_smod_xy * detector_y_dim}, + { 0.F, 0.F, 1.F, (rep_smod_z - n_smod_z / 2) * n_smod_z * detector_z_dim + + (rep_mod_z - n_mod_z / 2) * n_mod_z * n_smod_z * detector_z_dim + + (rep_rsec_z - n_rsec_z / 2) * n_rsec_z * n_mod_z * n_smod_z * detector_z_dim} } }; + + rep_module.ids.push_back(module_id++); + rep_module.transforms.push_back(transform); + } } petsird::ScannerGeometry scanner_geometry; scanner_geometry.replicated_modules.push_back(rep_module); diff --git a/data/root/scanner_geometry.json b/data/root/scanner_geometry.json index 9b8aa61..bdb9b3f 100644 --- a/data/root/scanner_geometry.json +++ b/data/root/scanner_geometry.json @@ -2,7 +2,7 @@ "ax_phys_crystal_num": 5, "ax_virtual_crystal_num": 0, "max_d_ring": 39, - "n_cry_layer":1, + "n_cry_layers":1, "n_cry_xy": 5, "n_cry_z": 5, "n_crystal": 25, @@ -17,6 +17,14 @@ "n_smod_xy": 1, "n_smod_z": 1, "n_submod": 1, + "cry_ax_gap":0, + "cry_tx_gap":0, + "smod_ax_gap":0, + "smod_tx_gap":0, + "mod_ax_gap":0, + "mod_tx_gap":0, + "rsec_ax_gap":0, + "rsec_tx_gap":0, "number_of_energy_bins": 1, "number_of_tof_bins": 62, "radius": 313.0, From d97ebc4e6a3a30f09d091cbc7ba2a5714de731bc Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Mon, 4 Nov 2024 05:12:29 +0000 Subject: [PATCH 05/13] fix CMake file and include petsird_helpers.h --- cpp/CMakeLists.txt | 11 +++++++---- cpp/main.cpp | 1 + 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 33c5a7a..e44ec9e 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -1,5 +1,5 @@ cmake_minimum_required(VERSION 3.12.0) # older would work, but could give warnings on policy CMP0074 -project(PETSIRDUseCaseTemplate VERSION 0.1.0) +project(root_to_petsird VERSION 0.2.0) set(CMAKE_CXX_STANDARD 17) @@ -16,11 +16,14 @@ else() add_compile_options(-Wall -Wextra -pedantic) endif() -add_subdirectory(../PETSIRD/cpp/generated PETSIRD_generated) +set(PETSIRD_dir ../PETSIRD/cpp/generated) +add_subdirectory(${PETSIRD_dir} PETSIRD_generated) find_package(ROOT REQUIRED COMPONENTS Core Imt RIO Net Hist Graf Graf3d Gpad Tree TreePlayer Rint Postscript Matrix Physics MathCore Thread MultiProc ROOTDataFrame) add_executable(root_to_petsird main.cpp) -include_directories(root_to_petsird PUBLIC ../PETSIRD/cpp/generated) -target_link_libraries(root_to_petsird PUBLIC prd_generated) +target_include_directories(root_to_petsird PUBLIC ${PETSIRD_dir}) +# needed for helpers +target_include_directories(root_to_petsird PUBLIC ${PETSIRD_dir}/..) +target_link_libraries(root_to_petsird PUBLIC petsird_generated) target_link_libraries(root_to_petsird PUBLIC ROOT::Core ROOT::Imt ROOT::RIO ROOT::Net ROOT::Hist ROOT::Graf ROOT::Graf3d ROOT::Gpad ROOT::Tree ROOT::TreePlayer ROOT::Rint ROOT::Postscript ROOT::Matrix ROOT::Physics ROOT::MathCore ROOT::Thread ROOT::MultiProc ROOT::ROOTDataFrame) install(TARGETS root_to_petsird DESTINATION bin) diff --git a/cpp/main.cpp b/cpp/main.cpp index 06b2f4f..e99c885 100644 --- a/cpp/main.cpp +++ b/cpp/main.cpp @@ -98,6 +98,7 @@ #include "protocols.h" #include "types.h" #include "binary/protocols.h" +#include "petsird_helpers.h" using namespace std ; From d259e3f875c6901f627ad5ea2f3f0d25851c420b Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Mon, 4 Nov 2024 05:18:19 +0000 Subject: [PATCH 06/13] [CI] remove python from tests --- .github/workflows/ci.yml | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 5b7b603..2c7caef 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -23,7 +23,7 @@ jobs: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 with: submodules: recursive @@ -46,11 +46,6 @@ jobs: cd PETSIRD/model yardl generate - - name: Python - run: | - cd python - python start.py - - name: Cpp run: | cd cpp && mkdir -p build && cd build From 26db7563e041533860d6991b0f8e730ee37cfdd3 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Mon, 4 Nov 2024 05:28:38 +0000 Subject: [PATCH 07/13] remove duplicate petsird_helpers.h --- cpp/main.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/cpp/main.cpp b/cpp/main.cpp index e99c885..06b2f4f 100644 --- a/cpp/main.cpp +++ b/cpp/main.cpp @@ -98,7 +98,6 @@ #include "protocols.h" #include "types.h" #include "binary/protocols.h" -#include "petsird_helpers.h" using namespace std ; From 75730a8f568466b4e4ed695850a2d1c05b053d4b Mon Sep 17 00:00:00 2001 From: nkarakatsanis Date: Mon, 4 Nov 2024 14:21:05 +0000 Subject: [PATCH 08/13] update v.0.1 --- cpp/main.cpp | 154 ++++++++++++++++++++++----------------------------- 1 file changed, 66 insertions(+), 88 deletions(-) diff --git a/cpp/main.cpp b/cpp/main.cpp index 06b2f4f..510ec61 100644 --- a/cpp/main.cpp +++ b/cpp/main.cpp @@ -15,21 +15,21 @@ // ********************************************************************* // //###################################################################################### -//# Authors : Nicolas A Karakatsanis, Sadek A. Nehmeh, CR Schmidtlein # +//# Authors : Nicolas A Karakatsanis, # //# # -//# Program : Bin_GATE_v1.0.c 29-JUL-2010 # +//# Date 03-NOV-2024 # //# # //# Objective : To read the coincidences TTree from the .root file, and generates the # //# corresponding Michelogram and Projection files. # //# # -//# Input : Monte Carlo data from GATE and egsPET # +//# Input : Monte Carlo data from GATE # //# # //# Output : 1 Michelogram files according to various binning definitions # //# : 2 Projection files according to various binning definitions # //# # //###################################################################################### //# # -//# This file is last modified on Nov 12, 2023 by: N. Karakatsanis # +//# This file is last modified on Nov 03, 2024 by: N. Karakatsanis # //# # //# The data are input from a root file produced by Gate simulating extended FOV of # //# mCT scanner. This scanner will have 5xFOV thus the root file contains information # @@ -38,11 +38,10 @@ //# # //# The virtual rings between the blocks are taken into consideration here. # //# # -//# The central FOV is taken into consideration => N_RINGS = 55 # +//# The central FOV is taken into consideration # //# The maximum and minimum rings should be specified if the user wishes to change # //# the number or the order of gantries. # //# # -//# The odd rings are removed.... # //# # //# # //# # @@ -101,8 +100,7 @@ using namespace std ; -#include "petsird_helpers.h" - +//#include "petsird_helpers.h" // these are constants for now constexpr uint32_t NUMBER_OF_ENERGY_BINS = 3; constexpr uint32_t NUMBER_OF_TOF_BINS = 300; @@ -158,6 +156,7 @@ struct ScannerGeometry float ArcLength; float TxFOV; float TxFOV_TOF; + float module_axial_pitch; }; void WriteScannerGeometry(const ScannerGeometry& scanner_geometry, const std::string& filename) @@ -207,6 +206,8 @@ void WriteScannerGeometry(const ScannerGeometry& scanner_geometry, const std::st j["ArcLength"] = scanner_geometry.s_width * scanner_geometry.detector_y_dim / 2.0f; j["TxFOV"] = 2 * scanner_geometry.radius * sin (scanner_geometry.ArcLength / (2 * scanner_geometry.radius) ); j["TxFOV_TOF"] = scanner_geometry.TxFOV + 0.3*scanner_geometry.TOF_resolution; + j["module_axial_pitch"] = scanner_geometry.n_cry_z * scanner_geometry.detector_z_dim + (scanner_geometry.n_cry_z - 1) * scanner_geometry.cry_ax_gap; + std::ofstream o(filename); o << std::setw(4) << j << std::endl; @@ -264,7 +265,7 @@ ScannerGeometry ReadScannerGeometry(const std::string& filename) scanner_geometry.ArcLength = scanner_geometry.s_width * scanner_geometry.detector_y_dim / 2.0f; scanner_geometry.TxFOV = 2 * scanner_geometry.radius * sin (scanner_geometry.ArcLength / (2 * scanner_geometry.radius) ); scanner_geometry.TxFOV_TOF = scanner_geometry.TxFOV + 0.3*scanner_geometry.TOF_resolution; - scanner_geometry.module_axial_pitch = n_cry_z * detector_z_dim + (n_cry_z - 1) * cry_ax_gap; + scanner_geometry.module_axial_pitch = scanner_geometry.n_cry_z * scanner_geometry.detector_z_dim + (scanner_geometry.n_cry_z - 1) * scanner_geometry.cry_ax_gap; return scanner_geometry; } @@ -321,11 +322,28 @@ get_scanner_info(ScannerGeometry& scannerGeometry) float energy_LLD = scannerGeometry.energy_LLD; float energy_ULD =scannerGeometry.energy_ULD; - std::vector angles; - for (int i = 0; i < n_detectors; ++i) - { - angles.push_back(static_cast(2 * M_PI * (1.0f*i) / n_detectors)); + typedef yardl::NDArray FArray1D; + // TOF info (in mm) + FArray1D::shape_type tof_bin_edges_shape = { NUMBER_OF_TOF_BINS + 1 }; + FArray1D tof_bin_edges(tof_bin_edges_shape); + for (std::size_t i = 0; i < tof_bin_edges.size(); ++i) { + tof_bin_edges[i] = (i - NUMBER_OF_TOF_BINS / 2.F) / NUMBER_OF_TOF_BINS * scannerGeometry.TxFOV_TOF; + } + FArray1D::shape_type energy_bin_edges_shape = { NUMBER_OF_ENERGY_BINS + 1 }; + FArray1D energy_bin_edges(energy_bin_edges_shape); + for (std::size_t i = 0; i < energy_bin_edges.size(); ++i) { + energy_bin_edges[i] = energy_LLD + i * (energy_ULD - energy_LLD) / NUMBER_OF_ENERGY_BINS; } + petsird::ScannerInformation scanner_info; + scanner_info.scanner_geometry = get_scanner_geometry(); + scanner_info.detectors = detectors; + scanner_info.tof_bin_edges = tof_bin_edges; + scanner_info.tof_resolution = scannerGeometry.TOF_resolution*0.3; // conversion from psec to mm (e.g. 200ps TOF is equivalent to 60mm uncertainty) + scanner_info.energy_bin_edges = energy_bin_edges; + scanner_info.energy_resolution_at_511 = scannerGeometry.EnergyResolutionAt511; // as fraction of 511 (e.g. 0.11F) + scanner_info.listmode_time_block_duration = scannerGeometry.LM_TimeBlockDuration; // ms + return scanner_info; +} //! return a cuboid volume petsird::BoxSolidVolume @@ -333,13 +351,13 @@ get_crystal() { using petsird::Coordinate; petsird::BoxShape crystal_shape{ Coordinate{ { 0, 0, 0 } }, - Coordinate{ { 0, 0, detector_z_dim } }, - Coordinate{ { 0, detector_y_dim, detector_z_dim } }, - Coordinate{ { 0, detector_y_dim, 0 } }, - Coordinate{ { detector_x_dim, 0, 0 } }, - Coordinate{ { detector_x_dim, 0, detector_z_dim } }, - Coordinate{ { detector_x_dim, detector_y_dim, detector_z_dim } }, - Coordinate{ { detector_x_dim, detector_y_dim, 0 } } }; + Coordinate{ { 0, 0, scannerGeometry.detector_z_dim } }, + Coordinate{ { 0, scannerGeometry.detector_y_dim, scannerGeometry.detector_z_dim } }, + Coordinate{ { 0, scannerGeometry.detector_y_dim, 0 } }, + Coordinate{ { scannerGeometry.detector_x_dim, 0, 0 } }, + Coordinate{ { scannerGeometry.detector_x_dim, 0, scannerGeometry.detector_z_dim } }, + Coordinate{ { scannerGeometry.detector_x_dim, scannerGeometry.detector_y_dim, scannerGeometry.detector_z_dim } }, + Coordinate{ { scannerGeometry.detector_x_dim, scannerGeometry.detector_y_dim, 0 } } }; petsird::BoxSolidVolume crystal{ crystal_shape, /* material_id */ 1 }; return crystal; @@ -356,9 +374,9 @@ get_detector_module() for (int rep_cry_xy = 0; rep_cry_xy < n_cry_xy; ++rep_cry_xy) for (int rep_cry_z = 0; rep_cry_z < n_cry_z; ++rep_cry_z) { - petsird::RigidTransformation transform{ { { 1.0, 0.0, 0.0, radius + rep_cry_layer * detector_x_dim }, - { 0.0, 1.0, 0.0, (rep_cry_xy - n_cry_xy / 2) * detector_y_dim }, - { 0.0, 0.0, 1.0, (rep_cry_z - n_cry_z / 2) * detector_z_dim } } }; + petsird::RigidTransformation transform{ { { 1.0, 0.0, 0.0, radius + rep_cry_layer * sscannerGeometry.detector_x_dim }, + { 0.0, 1.0, 0.0, (rep_cry_xy - n_cry_xy / 2) * scannerGeometry.detector_y_dim }, + { 0.0, 0.0, 1.0, (rep_cry_z - n_cry_z / 2) * scannerGeometry.detector_z_dim } } }; rep_volume.transforms.push_back(transform); //rep_volume.ids.push_back(rep_cry_layer + n_cry_layer * (rep_cry_xy + n_cry_xy * rep_cry_z)); rep_volume.ids.push_back(rep_cry_z + n_cry_z * (rep_cry_xy + n_cry_xy * rep_cry_layer)); @@ -388,23 +406,23 @@ get_scanner_geometry() angles.push_back(static_cast((2 * M_PI * i) / n_rsec_xy)); } - for (int rep_smod_xy = 0; rep_smod_xy < n_smod_xy; ++rep_smod_xy) - for (int rep_smod_z = 0; rep_smod_z < n_smod_z; ++rep_smod_z) - for (int rep_mod_xy = 0; rep_mod_xy < n_mod_xy; ++rep_mod_xy) - for (int rep_mod_z = 0; rep_mod_z < n_mod_z; ++rep_mod_z) - for (auto angle : angles) - for (int rep_rsec_z = 0; rep_rsec_z < n_rsec_z; ++rep_rsec_z) - { - petsird::RigidTransformation transform{ { { std::cos(angle), std::sin(angle), 0.F, 0.F }, - { -std::sin(angle), std::cos(angle), 0.F, (rep_smod_xy - n_smod_xy / 2) * n_smod_xy * detector_y_dim - + (rep_mod_xy - n_mod_xy / 2) * n_mod_xy * n_smod_xy * detector_y_dim}, - { 0.F, 0.F, 1.F, (rep_smod_z - n_smod_z / 2) * n_smod_z * detector_z_dim - + (rep_mod_z - n_mod_z / 2) * n_mod_z * n_smod_z * detector_z_dim - + (rep_rsec_z - n_rsec_z / 2) * n_rsec_z * n_mod_z * n_smod_z * detector_z_dim} } }; - - rep_module.ids.push_back(module_id++); - rep_module.transforms.push_back(transform); - } + for (int rep_smod_xy = 0; rep_smod_xy < n_smod_xy; ++rep_smod_xy) + for (int rep_smod_z = 0; rep_smod_z < n_smod_z; ++rep_smod_z) + for (int rep_mod_xy = 0; rep_mod_xy < n_mod_xy; ++rep_mod_xy) + for (int rep_mod_z = 0; rep_mod_z < n_mod_z; ++rep_mod_z) + for (auto angle : angles) + for (int rep_rsec_z = 0; rep_rsec_z < n_rsec_z; ++rep_rsec_z) + { + petsird::RigidTransformation transform{ { { std::cos(angle), std::sin(angle), 0.F, 0.F }, + { -std::sin(angle), std::cos(angle), 0.F, (rep_smod_xy - n_smod_xy / 2) * n_smod_xy * scannerGeometry.detector_y_dim + + (rep_mod_xy - n_mod_xy / 2) * n_mod_xy * n_smod_xy * scannerGeometry.detector_y_dim}, + { 0.F, 0.F, 1.F, (rep_smod_z - n_smod_z / 2) * n_smod_z * scannerGeometry.detector_z_dim + + (rep_mod_z - n_mod_z / 2) * n_mod_z * n_smod_z * scannerGeometry.detector_z_dim + + (rep_rsec_z - n_rsec_z / 2) * n_rsec_z * n_mod_z * n_smod_z * scannerGeometry.detector_z_dim} } }; + + rep_module.ids.push_back(module_id++); + rep_module.transforms.push_back(transform); + } } petsird::ScannerGeometry scanner_geometry; scanner_geometry.replicated_modules.push_back(rep_module); @@ -412,46 +430,6 @@ get_scanner_geometry() return scanner_geometry; } - - std::vector detectors; - int detector_id = 0; - for (int r =0; r < n_rings; r++) - { - for (auto angle : angles) - { - // Create a new detector - petsird::Detector d; - d.x = radius * std::cos(angle); - d.y = radius * std::sin(angle); - d.z = ((-n_rings/2.0f)*scannerGeometry.detector_z_dim) + scannerGeometry.detector_z_dim*r; - d.id = detector_id++; - detectors.push_back(d); - } - } - - typedef yardl::NDArray FArray1D; - // TOF info (in mm) - FArray1D::shape_type tof_bin_edges_shape = { NUMBER_OF_TOF_BINS + 1 }; - FArray1D tof_bin_edges(tof_bin_edges_shape); - for (std::size_t i = 0; i < tof_bin_edges.size(); ++i) { - tof_bin_edges[i] = (i - NUMBER_OF_TOF_BINS / 2.F) / NUMBER_OF_TOF_BINS * scannerGeometry.TxFOV_TOF; - } - FArray1D::shape_type energy_bin_edges_shape = { NUMBER_OF_ENERGY_BINS + 1 }; - FArray1D energy_bin_edges(energy_bin_edges_shape); - for (std::size_t i = 0; i < energy_bin_edges.size(); ++i) { - energy_bin_edges[i] = energy_LLD + i * (energy_ULD - energy_LLD) / NUMBER_OF_ENERGY_BINS; - } - petsird::ScannerInformation scanner_info; - scanner_info.scanner_geometry = get_scanner_geometry(); - //scanner_info.detectors = detectors; - scanner_info.tof_bin_edges = tof_bin_edges; - scanner_info.tof_resolution = scannerGeometry.TOF_resolution*0.3; // conversion from psec to mm (e.g. 200ps TOF is equivalent to 60mm uncertainty) - scanner_info.energy_bin_edges = energy_bin_edges; - scanner_info.energy_resolution_at_511 = scannerGeometry.EnergyResolutionAt511; // as fraction of 511 (e.g. 0.11F) - scanner_info.listmode_time_block_duration = scannerGeometry.LM_TimeBlockDuration; // ms - return scanner_info; -} - uint32_t tofToIdx(double delta_time_psec, const petsird::ScannerInformation& scanner_info) { float tofPos_mm = delta_time_psec * 0.15; //conversion from time difference (in psec) to spatial position in LOR (in mm) DT*C/2 @@ -681,14 +659,14 @@ int main(int argc, char** argv) std::cout << " tof_idx: " << event.tof_idx << std::endl; std::cout << " energy_1_idx: " << event.energy_indices[0] << std::endl; std::cout << " energy_2_idx: " << event.energy_indices[1] << std::endl; - std::cout << " detector 1 position: " << scanner.detectors[event.detector_ids[0]].x << ", " << scanner.detectors[event.detector_ids[0]].y << ", " << scanner.detectors[event.detector_ids[0]].z << std::endl; - std::cout << " GlobalPosition 1: " << globalPosX1 << ", " << globalPosY1 << ", " << globalPosZ1 << std::endl; - float distance_1 = std::sqrt(std::pow(scanner.detectors[event.detector_ids[0]].x-globalPosX1, 2) + std::pow(scanner.detectors[event.detector_ids[0]].y-globalPosY1, 2) + std::pow(scanner.detectors[event.detector_ids[0]].z-globalPosZ1, 2)); - std::cout << " Distance 1: " << distance_1 << std::endl; - std::cout << " detector 2 position: " << scanner.detectors[event.detector_ids[1]].x << ", " << scanner.detectors[event.detector_ids[1]].y << ", " << scanner.detectors[event.detector_ids[1]].z << std::endl; - std::cout << " GlobalPosition 2: " << globalPosX2 << ", " << globalPosY2 << ", " << globalPosZ2 << std::endl; - float distance_2 = std::sqrt(std::pow(scanner.detectors[event.detector_ids[1]].x-globalPosX2, 2) + std::pow(scanner.detectors[event.detector_ids[1]].y-globalPosY2, 2) + std::pow(scanner.detectors[event.detector_ids[1]].z-globalPosZ2, 2)); - std::cout << " Distance 2: " << distance_2 << std::endl; + //std::cout << " detector 1 position: " << scanner.detectors[event.detector_ids[0]].x << ", " << scanner.detectors[event.detector_ids[0]].y << ", " << scanner.detectors[event.detector_ids[0]].z << std::endl; + //std::cout << " GlobalPosition 1: " << globalPosX1 << ", " << globalPosY1 << ", " << globalPosZ1 << std::endl; + //float distance_1 = std::sqrt(std::pow(scanner.detectors[event.detector_ids[0]].x-globalPosX1, 2) + std::pow(scanner.detectors[event.detector_ids[0]].y-globalPosY1, 2) + std::pow(scanner.detectors[event.detector_ids[0]].z-globalPosZ1, 2)); + //std::cout << " Distance 1: " << distance_1 << std::endl; + //std::cout << " detector 2 position: " << scanner.detectors[event.detector_ids[1]].x << ", " << scanner.detectors[event.detector_ids[1]].y << ", " << scanner.detectors[event.detector_ids[1]].z << std::endl; + //std::cout << " GlobalPosition 2: " << globalPosX2 << ", " << globalPosY2 << ", " << globalPosZ2 << std::endl; + //float distance_2 = std::sqrt(std::pow(scanner.detectors[event.detector_ids[1]].x-globalPosX2, 2) + std::pow(scanner.detectors[event.detector_ids[1]].y-globalPosY2, 2) + std::pow(scanner.detectors[event.detector_ids[1]].z-globalPosZ2, 2)); + //std::cout << " Distance 2: " << distance_2 << std::endl; } long this_time_block = static_cast(time1*1.0e3 / scanner.event_time_block_duration); if (this_time_block != current_time_block) { From 7efeb80090ccdce1b6fa63e61e6bc84f640ce933 Mon Sep 17 00:00:00 2001 From: Joe Naegele Date: Mon, 4 Nov 2024 15:24:32 +0000 Subject: [PATCH 09/13] Fix last bugs --- cpp/main.cpp | 128 +++++++++++++++++++++++++-------------------------- 1 file changed, 64 insertions(+), 64 deletions(-) diff --git a/cpp/main.cpp b/cpp/main.cpp index 510ec61..f085539 100644 --- a/cpp/main.cpp +++ b/cpp/main.cpp @@ -294,7 +294,7 @@ int calculate_detector_id(int gantry_id, int rsector_id, int module_id, int subm int N_SMOD_z = scannerGeometry.n_smod_z; int N_CRY_xy = scannerGeometry.n_cry_xy; int N_CRY_z = scannerGeometry.n_cry_z; - int N_CRY_layers = scannerGeometry.n_cry_layers; + // int N_CRY_layers = scannerGeometry.n_cry_layers; int ring = (Int_t)(gantry_id)*N_RSEC_z*N_MOD_z*N_SMOD_z*N_CRY_z + (Int_t)(rsector_id/N_RSEC_xy)*N_MOD_z*N_SMOD_z*N_CRY_z @@ -310,44 +310,9 @@ int calculate_detector_id(int gantry_id, int rsector_id, int module_id, int subm return crystal + (ring-rmin)*N_DET; } -// single ring as example -petsird::ScannerInformation -get_scanner_info(ScannerGeometry& scannerGeometry) -{ - float radius = scannerGeometry.radius; - int n_detectors = scannerGeometry.n_det; - int n_rings = scannerGeometry.n_rings; - unsigned long NUMBER_OF_TOF_BINS = static_cast(scannerGeometry.number_of_tof_bins); - unsigned long NUMBER_OF_ENERGY_BINS = static_cast(scannerGeometry.number_of_energy_bins); - float energy_LLD = scannerGeometry.energy_LLD; - float energy_ULD =scannerGeometry.energy_ULD; - - typedef yardl::NDArray FArray1D; - // TOF info (in mm) - FArray1D::shape_type tof_bin_edges_shape = { NUMBER_OF_TOF_BINS + 1 }; - FArray1D tof_bin_edges(tof_bin_edges_shape); - for (std::size_t i = 0; i < tof_bin_edges.size(); ++i) { - tof_bin_edges[i] = (i - NUMBER_OF_TOF_BINS / 2.F) / NUMBER_OF_TOF_BINS * scannerGeometry.TxFOV_TOF; - } - FArray1D::shape_type energy_bin_edges_shape = { NUMBER_OF_ENERGY_BINS + 1 }; - FArray1D energy_bin_edges(energy_bin_edges_shape); - for (std::size_t i = 0; i < energy_bin_edges.size(); ++i) { - energy_bin_edges[i] = energy_LLD + i * (energy_ULD - energy_LLD) / NUMBER_OF_ENERGY_BINS; - } - petsird::ScannerInformation scanner_info; - scanner_info.scanner_geometry = get_scanner_geometry(); - scanner_info.detectors = detectors; - scanner_info.tof_bin_edges = tof_bin_edges; - scanner_info.tof_resolution = scannerGeometry.TOF_resolution*0.3; // conversion from psec to mm (e.g. 200ps TOF is equivalent to 60mm uncertainty) - scanner_info.energy_bin_edges = energy_bin_edges; - scanner_info.energy_resolution_at_511 = scannerGeometry.EnergyResolutionAt511; // as fraction of 511 (e.g. 0.11F) - scanner_info.listmode_time_block_duration = scannerGeometry.LM_TimeBlockDuration; // ms - return scanner_info; -} - //! return a cuboid volume petsird::BoxSolidVolume -get_crystal() +get_crystal(ScannerGeometry& scannerGeometry) { using petsird::Coordinate; petsird::BoxShape crystal_shape{ Coordinate{ { 0, 0, 0 } }, @@ -365,21 +330,21 @@ get_crystal() //! return a module of NUM_CRYSTALS_PER_MODULE cuboids petsird::DetectorModule -get_detector_module() +get_detector_module(ScannerGeometry& scannerGeometry) { petsird::ReplicatedBoxSolidVolume rep_volume; { - rep_volume.object = get_crystal(); - for (int rep_cry_layer = 0; rep_cry_layer < n_cry_layers; ++rep_cry_layer) - for (int rep_cry_xy = 0; rep_cry_xy < n_cry_xy; ++rep_cry_xy) - for (int rep_cry_z = 0; rep_cry_z < n_cry_z; ++rep_cry_z) + rep_volume.object = get_crystal(scannerGeometry); + for (int rep_cry_layer = 0; rep_cry_layer < scannerGeometry.n_cry_layers; ++rep_cry_layer) + for (int rep_cry_xy = 0; rep_cry_xy < scannerGeometry.n_cry_xy; ++rep_cry_xy) + for (int rep_cry_z = 0; rep_cry_z < scannerGeometry.n_cry_z; ++rep_cry_z) { - petsird::RigidTransformation transform{ { { 1.0, 0.0, 0.0, radius + rep_cry_layer * sscannerGeometry.detector_x_dim }, - { 0.0, 1.0, 0.0, (rep_cry_xy - n_cry_xy / 2) * scannerGeometry.detector_y_dim }, - { 0.0, 0.0, 1.0, (rep_cry_z - n_cry_z / 2) * scannerGeometry.detector_z_dim } } }; + petsird::RigidTransformation transform{ { { 1.0, 0.0, 0.0, scannerGeometry.radius + rep_cry_layer * scannerGeometry.detector_x_dim }, + { 0.0, 1.0, 0.0, (rep_cry_xy - scannerGeometry.n_cry_xy / 2) * scannerGeometry.detector_y_dim }, + { 0.0, 0.0, 1.0, (rep_cry_z - scannerGeometry.n_cry_z / 2) * scannerGeometry.detector_z_dim } } }; rep_volume.transforms.push_back(transform); //rep_volume.ids.push_back(rep_cry_layer + n_cry_layer * (rep_cry_xy + n_cry_xy * rep_cry_z)); - rep_volume.ids.push_back(rep_cry_z + n_cry_z * (rep_cry_xy + n_cry_xy * rep_cry_layer)); + rep_volume.ids.push_back(rep_cry_z + scannerGeometry.n_cry_z * (rep_cry_xy + scannerGeometry.n_cry_xy * rep_cry_layer)); } } @@ -391,34 +356,33 @@ get_detector_module() } - //! return scanner build by rotating a module around the (0,0,1) axis petsird::ScannerGeometry -get_scanner_geometry() +get_scanner_geometry(ScannerGeometry& scannerGeometry) { petsird::ReplicatedDetectorModule rep_module; { - rep_module.object = get_detector_module(); + rep_module.object = get_detector_module(scannerGeometry); int module_id = 0; std::vector angles; - for (unsigned int i = 0; i < n_rsec_xy; ++i) + for (int i = 0; i < scannerGeometry.n_rsec_xy; ++i) { - angles.push_back(static_cast((2 * M_PI * i) / n_rsec_xy)); + angles.push_back(static_cast((2 * M_PI * i) / scannerGeometry.n_rsec_xy)); } - for (int rep_smod_xy = 0; rep_smod_xy < n_smod_xy; ++rep_smod_xy) - for (int rep_smod_z = 0; rep_smod_z < n_smod_z; ++rep_smod_z) - for (int rep_mod_xy = 0; rep_mod_xy < n_mod_xy; ++rep_mod_xy) - for (int rep_mod_z = 0; rep_mod_z < n_mod_z; ++rep_mod_z) + for (int rep_smod_xy = 0; rep_smod_xy < scannerGeometry.n_smod_xy; ++rep_smod_xy) + for (int rep_smod_z = 0; rep_smod_z < scannerGeometry.n_smod_z; ++rep_smod_z) + for (int rep_mod_xy = 0; rep_mod_xy < scannerGeometry.n_mod_xy; ++rep_mod_xy) + for (int rep_mod_z = 0; rep_mod_z < scannerGeometry.n_mod_z; ++rep_mod_z) for (auto angle : angles) - for (int rep_rsec_z = 0; rep_rsec_z < n_rsec_z; ++rep_rsec_z) + for (int rep_rsec_z = 0; rep_rsec_z < scannerGeometry.n_rsec_z; ++rep_rsec_z) { petsird::RigidTransformation transform{ { { std::cos(angle), std::sin(angle), 0.F, 0.F }, - { -std::sin(angle), std::cos(angle), 0.F, (rep_smod_xy - n_smod_xy / 2) * n_smod_xy * scannerGeometry.detector_y_dim - + (rep_mod_xy - n_mod_xy / 2) * n_mod_xy * n_smod_xy * scannerGeometry.detector_y_dim}, - { 0.F, 0.F, 1.F, (rep_smod_z - n_smod_z / 2) * n_smod_z * scannerGeometry.detector_z_dim - + (rep_mod_z - n_mod_z / 2) * n_mod_z * n_smod_z * scannerGeometry.detector_z_dim - + (rep_rsec_z - n_rsec_z / 2) * n_rsec_z * n_mod_z * n_smod_z * scannerGeometry.detector_z_dim} } }; + { -std::sin(angle), std::cos(angle), 0.F, (rep_smod_xy - scannerGeometry.n_smod_xy / 2) * scannerGeometry.n_smod_xy * scannerGeometry.detector_y_dim + + (rep_mod_xy - scannerGeometry.n_mod_xy / 2) * scannerGeometry.n_mod_xy * scannerGeometry.n_smod_xy * scannerGeometry.detector_y_dim}, + { 0.F, 0.F, 1.F, (rep_smod_z - scannerGeometry.n_smod_z / 2) * scannerGeometry.n_smod_z * scannerGeometry.detector_z_dim + + (rep_mod_z - scannerGeometry.n_mod_z / 2) * scannerGeometry.n_mod_z * scannerGeometry.n_smod_z * scannerGeometry.detector_z_dim + + (rep_rsec_z - scannerGeometry.n_rsec_z / 2) * scannerGeometry.n_rsec_z * scannerGeometry.n_mod_z * scannerGeometry.n_smod_z * scannerGeometry.detector_z_dim} } }; rep_module.ids.push_back(module_id++); rep_module.transforms.push_back(transform); @@ -430,6 +394,43 @@ get_scanner_geometry() return scanner_geometry; } +// single ring as example +petsird::ScannerInformation +get_scanner_info(ScannerGeometry& scannerGeometry) +{ + // float radius = scannerGeometry.radius; + // int n_detectors = scannerGeometry.n_det; + // int n_rings = scannerGeometry.n_rings; + unsigned long NUMBER_OF_TOF_BINS = static_cast(scannerGeometry.number_of_tof_bins); + unsigned long NUMBER_OF_ENERGY_BINS = static_cast(scannerGeometry.number_of_energy_bins); + float energy_LLD = scannerGeometry.energy_LLD; + float energy_ULD =scannerGeometry.energy_ULD; + + typedef yardl::NDArray FArray1D; + // TOF info (in mm) + FArray1D::shape_type tof_bin_edges_shape = { NUMBER_OF_TOF_BINS + 1 }; + FArray1D tof_bin_edges(tof_bin_edges_shape); + for (std::size_t i = 0; i < tof_bin_edges.size(); ++i) { + tof_bin_edges[i] = (i - NUMBER_OF_TOF_BINS / 2.F) / NUMBER_OF_TOF_BINS * scannerGeometry.TxFOV_TOF; + } + FArray1D::shape_type energy_bin_edges_shape = { NUMBER_OF_ENERGY_BINS + 1 }; + FArray1D energy_bin_edges(energy_bin_edges_shape); + for (std::size_t i = 0; i < energy_bin_edges.size(); ++i) { + energy_bin_edges[i] = energy_LLD + i * (energy_ULD - energy_LLD) / NUMBER_OF_ENERGY_BINS; + } + petsird::ScannerInformation scanner_info; + scanner_info.scanner_geometry = get_scanner_geometry(scannerGeometry); + // scanner_info.detectors = detectors; + scanner_info.tof_bin_edges = tof_bin_edges; + scanner_info.tof_resolution = scannerGeometry.TOF_resolution*0.3; // conversion from psec to mm (e.g. 200ps TOF is equivalent to 60mm uncertainty) + scanner_info.energy_bin_edges = energy_bin_edges; + scanner_info.energy_resolution_at_511 = scannerGeometry.EnergyResolutionAt511; // as fraction of 511 (e.g. 0.11F) + scanner_info.event_time_block_duration = scannerGeometry.LM_TimeBlockDuration; // ms + return scanner_info; +} + + + uint32_t tofToIdx(double delta_time_psec, const petsird::ScannerInformation& scanner_info) { float tofPos_mm = delta_time_psec * 0.15; //conversion from time difference (in psec) to spatial position in LOR (in mm) DT*C/2 @@ -607,7 +608,6 @@ int main(int argc, char** argv) if (verbose) { // Print scanner information std::cout << "Scanner information:" << std::endl; - std::cout << "Number of detectors: " << scanner.NumberOfDetectors() << std::endl; std::cout << "Number of TOF bins: " << scanner.NumberOfTOFBins() << std::endl; std::cout << "Number of energy bins: " << scanner.NumberOfEnergyBins() << std::endl; const auto& tof_bin_edges = scanner.tof_bin_edges; @@ -625,7 +625,7 @@ int main(int argc, char** argv) header.scanner = scanner; // Write PETSiRD file - petsird::binary::petsirdExperimentWriter writer(petsird_file); + petsird::binary::PETSIRDWriter writer(petsird_file); writer.WriteHeader(header); long current_time_block = -1; @@ -674,7 +674,7 @@ int main(int argc, char** argv) writer.WriteTimeBlocks(time_block); } current_time_block = this_time_block; - time_block = petsird::TimeBlock(); + time_block = petsird::EventTimeBlock{}; time_block.start = time1*1.0e3; } time_block.prompt_events.push_back(event); From b693dfc9cdba073697ae3e8f5991f310c93ec037 Mon Sep 17 00:00:00 2001 From: nkarakatsanis Date: Mon, 4 Nov 2024 18:37:06 +0000 Subject: [PATCH 10/13] updated v0.2 --- cpp/main.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/main.cpp b/cpp/main.cpp index f085539..a569dc2 100644 --- a/cpp/main.cpp +++ b/cpp/main.cpp @@ -378,8 +378,8 @@ get_scanner_geometry(ScannerGeometry& scannerGeometry) for (int rep_rsec_z = 0; rep_rsec_z < scannerGeometry.n_rsec_z; ++rep_rsec_z) { petsird::RigidTransformation transform{ { { std::cos(angle), std::sin(angle), 0.F, 0.F }, - { -std::sin(angle), std::cos(angle), 0.F, (rep_smod_xy - scannerGeometry.n_smod_xy / 2) * scannerGeometry.n_smod_xy * scannerGeometry.detector_y_dim - + (rep_mod_xy - scannerGeometry.n_mod_xy / 2) * scannerGeometry.n_mod_xy * scannerGeometry.n_smod_xy * scannerGeometry.detector_y_dim}, + { -std::sin(angle), std::cos(angle), 0.F, std::sin(angle)*((rep_smod_xy - scannerGeometry.n_smod_xy / 2) * scannerGeometry.n_smod_xy * scannerGeometry.detector_y_dim + + (rep_mod_xy - scannerGeometry.n_mod_xy / 2) * scannerGeometry.n_mod_xy * scannerGeometry.n_smod_xy * scannerGeometry.detector_y_dim)}, { 0.F, 0.F, 1.F, (rep_smod_z - scannerGeometry.n_smod_z / 2) * scannerGeometry.n_smod_z * scannerGeometry.detector_z_dim + (rep_mod_z - scannerGeometry.n_mod_z / 2) * scannerGeometry.n_mod_z * scannerGeometry.n_smod_z * scannerGeometry.detector_z_dim + (rep_rsec_z - scannerGeometry.n_rsec_z / 2) * scannerGeometry.n_rsec_z * scannerGeometry.n_mod_z * scannerGeometry.n_smod_z * scannerGeometry.detector_z_dim} } }; From 5952db10efbad24d333c997889adb6469b9d7d36 Mon Sep 17 00:00:00 2001 From: nkarakatsanis Date: Mon, 4 Nov 2024 20:32:33 +0000 Subject: [PATCH 11/13] update v0.3 --- justfile | 25 +++++++++++++++++++++---- 1 file changed, 21 insertions(+), 4 deletions(-) diff --git a/justfile b/justfile index f760871..8ec14cd 100644 --- a/justfile +++ b/justfile @@ -25,8 +25,25 @@ default: run export AZURE_STORAGE_SAS_TOKEN="sp=rl&st=2023-11-14T18:39:07Z&se=2024-11-15T02:39:07Z&spr=https&sv=2022-11-02&sr=c&sig=0KVD7ORBM7Mx1%2BhVrVbqYcQycshhvT2XvdmrWVetiQM%3D"; \ dvc pull -@run: build download-test-data - cpp/build/root_to_petsird --root-prefix data/root/ETSIPETscanner_mIEC_ --number-of-root-files 1 --petsird-file petsird.bin -s data/root/scanner_geometry.json +#@run: build download-test-data +# cpp/build/root_to_petsird --root-prefix data/root/ETSIPETscanner_mIEC_ --number-of-root-files 1 --petsird-file petsird.bin -s data/root/scanner_geometry.json -@run-full: build download-test-data-full - cpp/build/root_to_petsird --root-prefix data/root/ETSIPETscanner_mIEC_ --number-of-root-files 36 --petsird-file petsird-full.bin -s data/root/scanner_geometry.json +#@run-full: build download-test-data-full +# cpp/build/root_to_petsird --root-prefix data/root/ETSIPETscanner_mIEC_ --number-of-root-files 36 --petsird-file petsird-full.bin -s data/root/scanner_geometry.json + + + +#@run_mIEC: build +# cpp/build/root_to_petsird --root-prefix ROOT_DATA/IEC/ETSIPETscanner_mIEC_ --number-of-root-files 1 --petsird-file petsird_mIEC_1_v0.2.bin -s data/root/scanner_geometry.json + +#@run-full_mIEC: build +# cpp/build/root_to_petsird --root-prefix ROOT_DATA/IEC/ETSIPETscanner_mIEC_ --number-of-root-files 36 --petsird-file petsird-full_mIEC_10_v0.2.bin -s data/root/scanner_geometry.json + + +@run_voxBrain: build + cpp/build/root_to_petsird --root-prefix ROOT_DATA/voxBrain/ETSIPETscanner256mmAFOV_positronSource/ETSIPETscanner256mmAFOV_voxBrain_ --number-of-root-files 1 --petsird-file petsird_voxBrain_1_v0.2.bin -s data/root/scanner_geometry.json + +@run-full_voxBrain: build + cpp/build/root_to_petsird --root-prefix ROOT_DATA/voxBrain/ETSIPETscanner256mmAFOV_positronSource/ETSIPETscanner256mmAFOV_voxBrain_ --number-of-root-files 36 --petsird-file petsird-full_voxBrain_10_v0.2.bin -s data/root/scanner_geometry.json + +@run: run-full_voxBrain From a38837717a94a250e56a21cada7b5c0a798b034e Mon Sep 17 00:00:00 2001 From: nkarakatsanis Date: Tue, 5 Nov 2024 07:33:35 +0000 Subject: [PATCH 12/13] Geometry description updated v0.2 --- cpp/main.cpp | 71 +++++++++++++++++++++++----------------------------- justfile | 31 +++++++++++------------ 2 files changed, 46 insertions(+), 56 deletions(-) diff --git a/cpp/main.cpp b/cpp/main.cpp index a569dc2..3bbce99 100644 --- a/cpp/main.cpp +++ b/cpp/main.cpp @@ -101,18 +101,6 @@ using namespace std ; //#include "petsird_helpers.h" -// these are constants for now -constexpr uint32_t NUMBER_OF_ENERGY_BINS = 3; -constexpr uint32_t NUMBER_OF_TOF_BINS = 300; -constexpr float RADIUS = 400.F; -constexpr std::array CRYSTAL_LENGTH{ 20.F, 4.F, 4.F }; -constexpr std::array NUM_CRYSTALS_PER_MODULE{ 2, 4, 5 }; -constexpr uint32_t NUM_MODULES_ALONG_RING{ 20 }; -constexpr uint32_t NUM_MODULES_ALONG_AXIS{ 2 }; -constexpr float MODULE_AXIS_SPACING{ (NUM_CRYSTALS_PER_MODULE[2] + 4) * CRYSTAL_LENGTH[2] }; - -constexpr uint32_t NUMBER_OF_TIME_BLOCKS = 6; -constexpr float COUNT_RATE = 500.F; struct ScannerGeometry { @@ -335,17 +323,27 @@ get_detector_module(ScannerGeometry& scannerGeometry) petsird::ReplicatedBoxSolidVolume rep_volume; { rep_volume.object = get_crystal(scannerGeometry); - for (int rep_cry_layer = 0; rep_cry_layer < scannerGeometry.n_cry_layers; ++rep_cry_layer) - for (int rep_cry_xy = 0; rep_cry_xy < scannerGeometry.n_cry_xy; ++rep_cry_xy) - for (int rep_cry_z = 0; rep_cry_z < scannerGeometry.n_cry_z; ++rep_cry_z) - { - petsird::RigidTransformation transform{ { { 1.0, 0.0, 0.0, scannerGeometry.radius + rep_cry_layer * scannerGeometry.detector_x_dim }, - { 0.0, 1.0, 0.0, (rep_cry_xy - scannerGeometry.n_cry_xy / 2) * scannerGeometry.detector_y_dim }, - { 0.0, 0.0, 1.0, (rep_cry_z - scannerGeometry.n_cry_z / 2) * scannerGeometry.detector_z_dim } } }; - rep_volume.transforms.push_back(transform); - //rep_volume.ids.push_back(rep_cry_layer + n_cry_layer * (rep_cry_xy + n_cry_xy * rep_cry_z)); - rep_volume.ids.push_back(rep_cry_z + scannerGeometry.n_cry_z * (rep_cry_xy + scannerGeometry.n_cry_xy * rep_cry_layer)); - } + int crystal_id = 0; + for (int rep_mod_xy = 0; rep_mod_xy < scannerGeometry.n_mod_xy; ++rep_mod_xy) + for (int rep_mod_z = 0; rep_mod_z < scannerGeometry.n_mod_z; ++rep_mod_z) + for (int rep_smod_xy = 0; rep_smod_xy < scannerGeometry.n_smod_xy; ++rep_smod_xy) + for (int rep_smod_z = 0; rep_smod_z < scannerGeometry.n_smod_z; ++rep_smod_z) + for (int rep_cry_xy = 0; rep_cry_xy < scannerGeometry.n_cry_xy; ++rep_cry_xy) + for (int rep_cry_z = 0; rep_cry_z < scannerGeometry.n_cry_z; ++rep_cry_z) + for (int rep_cry_layer = 0; rep_cry_layer < scannerGeometry.n_cry_layers; ++rep_cry_layer) + { + petsird::RigidTransformation transform{ { { 1.0, 0.0, 0.0, scannerGeometry.radius + rep_cry_layer * scannerGeometry.detector_x_dim }, + { 0.0, 1.0, 0.0, (rep_mod_xy - scannerGeometry.n_mod_xy / 2) * scannerGeometry.n_smod_xy * scannerGeometry.n_cry_xy * scannerGeometry.detector_y_dim + + (rep_smod_xy - scannerGeometry.n_smod_xy / 2) * scannerGeometry.n_cry_xy * scannerGeometry.detector_y_dim + + (rep_cry_xy - scannerGeometry.n_cry_xy / 2) * scannerGeometry.detector_y_dim }, + { 0.0, 0.0, 1.0, (rep_mod_z - scannerGeometry.n_mod_z / 2) * scannerGeometry.n_smod_z * scannerGeometry.n_cry_z * scannerGeometry.detector_z_dim + + (rep_smod_z - scannerGeometry.n_smod_z / 2) * scannerGeometry.n_cry_z * scannerGeometry.detector_z_dim + + (rep_cry_z - scannerGeometry.n_cry_z / 2) * scannerGeometry.detector_z_dim } } }; + rep_volume.transforms.push_back(transform); + //rep_volume.ids.push_back(rep_cry_layer + n_cry_layer * (rep_cry_xy + n_cry_xy * rep_cry_z)); + //rep_volume.ids.push_back(rep_cry_z + scannerGeometry.n_cry_z * (rep_cry_xy + scannerGeometry.n_cry_xy * rep_cry_layer)); + rep_volume.ids.push_back(crystal_id++); + } } petsird::DetectorModule detector_module; @@ -370,23 +368,16 @@ get_scanner_geometry(ScannerGeometry& scannerGeometry) angles.push_back(static_cast((2 * M_PI * i) / scannerGeometry.n_rsec_xy)); } - for (int rep_smod_xy = 0; rep_smod_xy < scannerGeometry.n_smod_xy; ++rep_smod_xy) - for (int rep_smod_z = 0; rep_smod_z < scannerGeometry.n_smod_z; ++rep_smod_z) - for (int rep_mod_xy = 0; rep_mod_xy < scannerGeometry.n_mod_xy; ++rep_mod_xy) - for (int rep_mod_z = 0; rep_mod_z < scannerGeometry.n_mod_z; ++rep_mod_z) - for (auto angle : angles) - for (int rep_rsec_z = 0; rep_rsec_z < scannerGeometry.n_rsec_z; ++rep_rsec_z) - { - petsird::RigidTransformation transform{ { { std::cos(angle), std::sin(angle), 0.F, 0.F }, - { -std::sin(angle), std::cos(angle), 0.F, std::sin(angle)*((rep_smod_xy - scannerGeometry.n_smod_xy / 2) * scannerGeometry.n_smod_xy * scannerGeometry.detector_y_dim - + (rep_mod_xy - scannerGeometry.n_mod_xy / 2) * scannerGeometry.n_mod_xy * scannerGeometry.n_smod_xy * scannerGeometry.detector_y_dim)}, - { 0.F, 0.F, 1.F, (rep_smod_z - scannerGeometry.n_smod_z / 2) * scannerGeometry.n_smod_z * scannerGeometry.detector_z_dim - + (rep_mod_z - scannerGeometry.n_mod_z / 2) * scannerGeometry.n_mod_z * scannerGeometry.n_smod_z * scannerGeometry.detector_z_dim - + (rep_rsec_z - scannerGeometry.n_rsec_z / 2) * scannerGeometry.n_rsec_z * scannerGeometry.n_mod_z * scannerGeometry.n_smod_z * scannerGeometry.detector_z_dim} } }; - - rep_module.ids.push_back(module_id++); - rep_module.transforms.push_back(transform); - } + for (auto angle : angles) + for (int rep_rsec_z = 0; rep_rsec_z < scannerGeometry.n_rsec_z; ++rep_rsec_z) + { + petsird::RigidTransformation transform{ { { std::cos(angle), std::sin(angle), 0.F, 0.F }, + { -std::sin(angle), std::cos(angle), 0.F, 0.F}, + { 0.F, 0.F, 1.F, (rep_rsec_z - scannerGeometry.n_rsec_z / 2) * scannerGeometry.n_mod_z * scannerGeometry.n_smod_z * scannerGeometry.n_cry_z * scannerGeometry.detector_z_dim} } }; + + rep_module.ids.push_back(module_id++); + rep_module.transforms.push_back(transform); + } } petsird::ScannerGeometry scanner_geometry; scanner_geometry.replicated_modules.push_back(rep_module); diff --git a/justfile b/justfile index 8ec14cd..a3b9366 100644 --- a/justfile +++ b/justfile @@ -17,13 +17,6 @@ default: run cd cpp/build; \ ninja -@download-test-data: - export AZURE_STORAGE_SAS_TOKEN="sp=rl&st=2023-11-14T18:39:07Z&se=2024-11-15T02:39:07Z&spr=https&sv=2022-11-02&sr=c&sig=0KVD7ORBM7Mx1%2BhVrVbqYcQycshhvT2XvdmrWVetiQM%3D"; \ - dvc pull data/root/ETSIPETscanner_mIEC_1.root - -@download-test-data-full: - export AZURE_STORAGE_SAS_TOKEN="sp=rl&st=2023-11-14T18:39:07Z&se=2024-11-15T02:39:07Z&spr=https&sv=2022-11-02&sr=c&sig=0KVD7ORBM7Mx1%2BhVrVbqYcQycshhvT2XvdmrWVetiQM%3D"; \ - dvc pull #@run: build download-test-data # cpp/build/root_to_petsird --root-prefix data/root/ETSIPETscanner_mIEC_ --number-of-root-files 1 --petsird-file petsird.bin -s data/root/scanner_geometry.json @@ -33,17 +26,23 @@ default: run -#@run_mIEC: build -# cpp/build/root_to_petsird --root-prefix ROOT_DATA/IEC/ETSIPETscanner_mIEC_ --number-of-root-files 1 --petsird-file petsird_mIEC_1_v0.2.bin -s data/root/scanner_geometry.json +@run_0min_mIEC: build + cpp/build/root_to_petsird --root-prefix ROOT_DATA/IEC/ETSIPETscanner_mIEC_ --number-of-root-files 0 --petsird-file petsird_ETSIPETscanner_mIEC_0min_v0.2.bin -s data/root/scanner_geometry.json + +@run_10sec_mIEC: build + cpp/build/root_to_petsird --root-prefix ROOT_DATA/IEC/ETSIPETscanner_mIEC_ --number-of-root-files 1 --petsird-file petsird_ETSIPETscanner_mIEC_10sec_v0.2.bin -s data/root/scanner_geometry.json + +@run-6min_mIEC: build + cpp/build/root_to_petsird --root-prefix ROOT_DATA/IEC/ETSIPETscanner_mIEC_ --number-of-root-files 36 --petsird-file petsird_ETSIPETscanner_mIEC_6min_v0.2.bin -s data/root/scanner_geometry.json -#@run-full_mIEC: build -# cpp/build/root_to_petsird --root-prefix ROOT_DATA/IEC/ETSIPETscanner_mIEC_ --number-of-root-files 36 --petsird-file petsird-full_mIEC_10_v0.2.bin -s data/root/scanner_geometry.json +@run_0min_voxBrain: build + cpp/build/root_to_petsird --root-prefix ROOT_DATA/voxBrain/ETSIPETscanner256mmAFOV_positronSource/ETSIPETscanner256mmAFOV_voxBrain_ --number-of-root-files 0 --petsird-file petsird_ETSIPETscanner256mmAFOV_voxBrain_0min_v0.2.bin -s data/root/scanner_geometry.json -@run_voxBrain: build - cpp/build/root_to_petsird --root-prefix ROOT_DATA/voxBrain/ETSIPETscanner256mmAFOV_positronSource/ETSIPETscanner256mmAFOV_voxBrain_ --number-of-root-files 1 --petsird-file petsird_voxBrain_1_v0.2.bin -s data/root/scanner_geometry.json +@run_3min_voxBrain: build + cpp/build/root_to_petsird --root-prefix ROOT_DATA/voxBrain/ETSIPETscanner256mmAFOV_positronSource/ETSIPETscanner256mmAFOV_voxBrain_ --number-of-root-files 18 --petsird-file petsird_ETSIPETscanner256mmAFOV_voxBrain_3min_v0.2.bin -s data/root/scanner_geometry.json -@run-full_voxBrain: build - cpp/build/root_to_petsird --root-prefix ROOT_DATA/voxBrain/ETSIPETscanner256mmAFOV_positronSource/ETSIPETscanner256mmAFOV_voxBrain_ --number-of-root-files 36 --petsird-file petsird-full_voxBrain_10_v0.2.bin -s data/root/scanner_geometry.json +@run_6min_voxBrain: build + cpp/build/root_to_petsird --root-prefix ROOT_DATA/voxBrain/ETSIPETscanner256mmAFOV_positronSource/ETSIPETscanner256mmAFOV_voxBrain_ --number-of-root-files 36 --petsird-file petsird_ETSIPETscanner256mmAFOV_voxBrain_6min_v0.2.bin -s data/root/scanner_geometry.json -@run: run-full_voxBrain +@run: run_0min_voxBrain From 3a371201cfbbf239311ab0426d67a789392fe827 Mon Sep 17 00:00:00 2001 From: nkarakatsanis Date: Tue, 5 Nov 2024 17:41:27 +0000 Subject: [PATCH 13/13] v0.2 Fixed Geometry Description --- .../ETSIPETscanner256mmAFOV_geometry.json | 42 +++++++++++++++++++ data/root/ETSIPETscanner_geometry.json | 42 +++++++++++++++++++ justfile | 6 +-- 3 files changed, 87 insertions(+), 3 deletions(-) create mode 100644 data/root/ETSIPETscanner256mmAFOV_geometry.json create mode 100644 data/root/ETSIPETscanner_geometry.json diff --git a/data/root/ETSIPETscanner256mmAFOV_geometry.json b/data/root/ETSIPETscanner256mmAFOV_geometry.json new file mode 100644 index 0000000..4e18b22 --- /dev/null +++ b/data/root/ETSIPETscanner256mmAFOV_geometry.json @@ -0,0 +1,42 @@ +{ + "ax_phys_crystal_num": 5, + "ax_virtual_crystal_num": 0, + "max_d_ring": 79, + "n_cry_layers": 1, + "n_cry_xy": 5, + "n_cry_z": 5, + "n_crystal": 25, + "n_det": 600, + "n_mod_xy": 2, + "n_mod_z": 16, + "n_module": 32, + "n_rings": 80, + "n_rsec": 60, + "n_rsec_xy": 60, + "n_rsec_z": 1, + "n_smod_xy": 1, + "n_smod_z": 1, + "n_submod": 1, + "cry_ax_gap": 0, + "cry_tx_gap": 0, + "smod_ax_gap": 0, + "smod_tx_gap": 0, + "mod_ax_gap": 0, + "mod_tx_gap": 0, + "rsec_ax_gap": 0, + "rsec_tx_gap": 0, + "number_of_energy_bins": 1, + "number_of_tof_bins": 62, + "radius": 313.0, + "s_width": 380, + "tx_phys_crystal_num": 10, + "tx_virtual_crystal_num": 0, + "detector_x_dim": 20.0, + "detector_y_dim": 3.2, + "detector_z_dim": 3.2, + "energy_LLD": 430.0, + "energy_ULD": 650.0, + "EnergyResolutionAt511": 0.11, + "TOF_resolution": 200, + "LM_TimeBlockDuration": 1 +} diff --git a/data/root/ETSIPETscanner_geometry.json b/data/root/ETSIPETscanner_geometry.json new file mode 100644 index 0000000..e73d2b9 --- /dev/null +++ b/data/root/ETSIPETscanner_geometry.json @@ -0,0 +1,42 @@ +{ + "ax_phys_crystal_num": 5, + "ax_virtual_crystal_num": 0, + "max_d_ring": 39, + "n_cry_layers": 1, + "n_cry_xy": 5, + "n_cry_z": 5, + "n_crystal": 25, + "n_det": 600, + "n_mod_xy": 2, + "n_mod_z": 8, + "n_module": 16, + "n_rings": 40, + "n_rsec": 60, + "n_rsec_xy": 60, + "n_rsec_z": 1, + "n_smod_xy": 1, + "n_smod_z": 1, + "n_submod": 1, + "cry_ax_gap": 0, + "cry_tx_gap": 0, + "smod_ax_gap": 0, + "smod_tx_gap": 0, + "mod_ax_gap": 0, + "mod_tx_gap": 0, + "rsec_ax_gap": 0, + "rsec_tx_gap": 0, + "number_of_energy_bins": 1, + "number_of_tof_bins": 62, + "radius": 313.0, + "s_width": 380, + "tx_phys_crystal_num": 10, + "tx_virtual_crystal_num": 0, + "detector_x_dim": 20.0, + "detector_y_dim": 3.2, + "detector_z_dim": 3.2, + "energy_LLD": 430.0, + "energy_ULD": 650.0, + "EnergyResolutionAt511": 0.11, + "TOF_resolution": 200, + "LM_TimeBlockDuration": 1 +} diff --git a/justfile b/justfile index a3b9366..4560934 100644 --- a/justfile +++ b/justfile @@ -37,12 +37,12 @@ default: run @run_0min_voxBrain: build - cpp/build/root_to_petsird --root-prefix ROOT_DATA/voxBrain/ETSIPETscanner256mmAFOV_positronSource/ETSIPETscanner256mmAFOV_voxBrain_ --number-of-root-files 0 --petsird-file petsird_ETSIPETscanner256mmAFOV_voxBrain_0min_v0.2.bin -s data/root/scanner_geometry.json + cpp/build/root_to_petsird --root-prefix ROOT_DATA/voxBrain/ETSIPETscanner256mmAFOV_positronSource/ETSIPETscanner256mmAFOV_voxBrain_ --number-of-root-files 0 --petsird-file petsird_ETSIPETscanner256mmAFOV_voxBrain_0min_v0.2.bin -s data/root/ETSIPETscanner256mmAFOV_geometry.json @run_3min_voxBrain: build - cpp/build/root_to_petsird --root-prefix ROOT_DATA/voxBrain/ETSIPETscanner256mmAFOV_positronSource/ETSIPETscanner256mmAFOV_voxBrain_ --number-of-root-files 18 --petsird-file petsird_ETSIPETscanner256mmAFOV_voxBrain_3min_v0.2.bin -s data/root/scanner_geometry.json + cpp/build/root_to_petsird --root-prefix ROOT_DATA/voxBrain/ETSIPETscanner256mmAFOV_positronSource/ETSIPETscanner256mmAFOV_voxBrain_ --number-of-root-files 18 --petsird-file petsird_ETSIPETscanner256mmAFOV_voxBrain_3min_v0.2.bin -s data/root/ETSIPETscanner256mmAFOV_geometry.json @run_6min_voxBrain: build - cpp/build/root_to_petsird --root-prefix ROOT_DATA/voxBrain/ETSIPETscanner256mmAFOV_positronSource/ETSIPETscanner256mmAFOV_voxBrain_ --number-of-root-files 36 --petsird-file petsird_ETSIPETscanner256mmAFOV_voxBrain_6min_v0.2.bin -s data/root/scanner_geometry.json + cpp/build/root_to_petsird --root-prefix ROOT_DATA/voxBrain/ETSIPETscanner256mmAFOV_positronSource/ETSIPETscanner256mmAFOV_voxBrain_ --number-of-root-files 36 --petsird-file petsird_ETSIPETscanner256mmAFOV_voxBrain_6min_v0.2.bin -s data/root/ETSIPETscanner256mmAFOV_geometry.json @run: run_0min_voxBrain