From b9688883f1c0eb9058d49fb4b939e434477eced2 Mon Sep 17 00:00:00 2001 From: Claudia Date: Tue, 21 Jul 2026 12:59:38 +0200 Subject: [PATCH 1/4] Initial commit --- assets/Data.zip | 4 +- configs/config_test_wholebody_reference.yaml | 22 + .../scanner/Bruker_PET_MRI_scanner.yaml | 39 + .../params/scanner/Neuro_LF/NeuroLF_BPET.yaml | 38 + .../Neuro_LF/NeuroLF_scanner_10mm.yaml | 38 + .../Neuro_LF/NeuroLF_scanner_20mm.yaml | 38 + .../Neuro_LF/NeuroLF_scanner_25mm.yaml | 38 + .../Neuro_LF/NeuroLF_scanner_30mm.yaml | 38 + ...NeuroLF_scanner_3mm_crystal_thickness.yaml | 38 + .../NeuroLF_scanner_8bloques_axial_FOV.yaml | 38 + .../Neuro_LF/NeuroLF_scanner_cristal_2mm.yaml | 38 + .../Neuro_LF/NeuroLF_scanner_final.yaml | 38 + .../params/scanner/Vereos_PET_CT_scanner.yaml | 47 + configs/params/scanner/discovery_st.yaml | 2 +- configs/params/scanner/ls | 0 configs/params/test.yaml | 2 +- configs/params/test_wholebody_reference.yaml | 38 + makefile | 256 ---- requirements.txt | 2 +- scripts/experiment_wholebody.py | 3 +- src/simset/simset_sim.py | 49 +- src/simset/simset_tool.py | 0 src/simset/simset_tools.py | 9 +- src/stir/stir_sim.py | 3 + utils/Quantification_externo_Claudia.py | 407 +++++++ utils/wb_tools.py | 1083 ++++++++++++++++- wholebody.py | 237 +++- 27 files changed, 2215 insertions(+), 330 deletions(-) mode change 100755 => 100644 assets/Data.zip create mode 100644 configs/config_test_wholebody_reference.yaml create mode 100644 configs/params/scanner/Bruker_PET_MRI_scanner.yaml create mode 100644 configs/params/scanner/Neuro_LF/NeuroLF_BPET.yaml create mode 100644 configs/params/scanner/Neuro_LF/NeuroLF_scanner_10mm.yaml create mode 100644 configs/params/scanner/Neuro_LF/NeuroLF_scanner_20mm.yaml create mode 100644 configs/params/scanner/Neuro_LF/NeuroLF_scanner_25mm.yaml create mode 100644 configs/params/scanner/Neuro_LF/NeuroLF_scanner_30mm.yaml create mode 100644 configs/params/scanner/Neuro_LF/NeuroLF_scanner_3mm_crystal_thickness.yaml create mode 100644 configs/params/scanner/Neuro_LF/NeuroLF_scanner_8bloques_axial_FOV.yaml create mode 100644 configs/params/scanner/Neuro_LF/NeuroLF_scanner_cristal_2mm.yaml create mode 100644 configs/params/scanner/Neuro_LF/NeuroLF_scanner_final.yaml create mode 100644 configs/params/scanner/Vereos_PET_CT_scanner.yaml create mode 100644 configs/params/scanner/ls create mode 100644 configs/params/test_wholebody_reference.yaml delete mode 100644 makefile create mode 100644 src/simset/simset_tool.py create mode 100644 utils/Quantification_externo_Claudia.py diff --git a/assets/Data.zip b/assets/Data.zip old mode 100755 new mode 100644 index c9566ccc..5e898552 --- a/assets/Data.zip +++ b/assets/Data.zip @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:a7519a12bc1a21f260f0790d8eb3d1a2829b7422f972855f85efa255ef9b1658 -size 2353666 +oid sha256:27ccb5963664e7fa69b71d6db631a55597dc97447a0564eb65b4907cdb124167 +size 327152590 diff --git a/configs/config_test_wholebody_reference.yaml b/configs/config_test_wholebody_reference.yaml new file mode 100644 index 00000000..7000c2f7 --- /dev/null +++ b/configs/config_test_wholebody_reference.yaml @@ -0,0 +1,22 @@ +defaults: + - params: test_wholebody_reference +interactive_mode: 0 +dir_stir: ${root_path:root}/include/STIR/install +dir_simset: ${root_path:root}/include/SimSET/2.9.2 +matlab_mcr_path: "" +spm_path: "" +dir_data_path: ${root_path:root}/Data/Subjets +dir_results_path: ${root_path:root}/Data/Results +dir_stir_sino_corrections: ${root_path:root}/Normalization_sino_correction_Bruker_PET-MRI +stratification: "true" +forced_detection: "true" +forced_non_absortion: "true" +acceptance_angle: 90.0 +positron_range: "true" +isotope: "f18" +non_colinearity: "true" +minimum_energy: 350.0 +weight_window_ratio: 1.0 +point_source_voxels: "false" +coherent_scatter_object: "false" +coherent_scatter_detector: "false" diff --git a/configs/params/scanner/Bruker_PET_MRI_scanner.yaml b/configs/params/scanner/Bruker_PET_MRI_scanner.yaml new file mode 100644 index 00000000..915d3c01 --- /dev/null +++ b/configs/params/scanner/Bruker_PET_MRI_scanner.yaml @@ -0,0 +1,39 @@ +scanner_name: "Bruker_PET_MRI_scanner" +simset_material: 29 +average_doi: 0.84 +scanner_radius: 5.25 +num_rings: 100 +axial_fov: 15 +z_crystal_size: 0.1497 +transaxial_crystal_size: 0.241 +crystal_thickness: 1 +energy_resolution: 17 +num_aa_bins: 68 +num_td_bins: 136 +min_energy_window: 409 #20% +max_energy_window: 613 +coincidence_window: 5 #5 +numberOfSubsets: 1 +numberOfIterations: 48 +savingInterval: 16 +analytical_att_correction: 0 +stir_recons_att_corr: 1 +analytic_scatt_corr_factor: 0.15 +stir_scatt_corr_smoothing: 0 +stir_scatt_simulation: 0 +analytic_randoms_corr_factor: 0.2 +stir_randoms_corr_smoothing: 0 +recons_type: "OSEM3D" +inter_iteration_filter: 1 +subiteration_interval: 48 +x_dir_filter_FWHM: 0.7 +y_dir_filter_FWHM: 1 +z_dir_filter_FWHM: 0 +psf_value: 0 +add_noise: 0 +max_segment: 99 +zoomFactor: 1.0294 +xyOutputSize: 120 +zOutputSize: 200 #200 +zOutputVoxelSize: 3.27 + diff --git a/configs/params/scanner/Neuro_LF/NeuroLF_BPET.yaml b/configs/params/scanner/Neuro_LF/NeuroLF_BPET.yaml new file mode 100644 index 00000000..23260ae6 --- /dev/null +++ b/configs/params/scanner/Neuro_LF/NeuroLF_BPET.yaml @@ -0,0 +1,38 @@ +scanner_name: "NeuroLF_BPET" +simset_material: 29 +average_doi: 0.84 +scanner_radius: 13.4 +num_rings: 48 +axial_fov: 16.318 +z_crystal_size: 0.319 +transaxial_crystal_size: 0.319 +crystal_thickness: 1 +energy_resolution: 14 +num_aa_bins: 128 +num_td_bins: 256 +min_energy_window: 425 +max_energy_window: 650 +coincidence_window: 8 #5 #falta por revisar +numberOfSubsets: 8 +numberOfIterations: 128 #48 +savingInterval: 2 #8 +analytical_att_correction: 0 #falta por revisar +stir_recons_att_corr: 1 #1 esta corrigiendo por attenuacion; 0 no esta corrigiendo +analytic_scatt_corr_factor: 0.15 #falta por revisar +stir_scatt_corr_smoothing: 0 #falta por revisar +stir_scatt_simulation: 0 #falta por revisar +analytic_randoms_corr_factor: 0.2 #falta por revisar +stir_randoms_corr_smoothing: 0 #falta por revisar +recons_type: "OSEM3D" +inter_iteration_filter: 1 #1 se está filtrando +subiteration_interval: 4 +x_dir_filter_FWHM: 1.30 #valor inicial 1.5 #rango aceptable 1.20 a 1.45 +y_dir_filter_FWHM: 1.30 #valor inicial 1.5 #rango aceptable 1.20 a 1.45 +z_dir_filter_FWHM: 0.50 #valor inicial 1.5 +psf_value: 0 +add_noise: 0 +max_segment: 47 +zoomFactor: 1.264 +xyOutputSize: 325 +zOutputSize: 96 +zOutputVoxelSize: 1.64 #No hace nada diff --git a/configs/params/scanner/Neuro_LF/NeuroLF_scanner_10mm.yaml b/configs/params/scanner/Neuro_LF/NeuroLF_scanner_10mm.yaml new file mode 100644 index 00000000..d2bac01a --- /dev/null +++ b/configs/params/scanner/Neuro_LF/NeuroLF_scanner_10mm.yaml @@ -0,0 +1,38 @@ +scanner_name: "NeuroLF_BPET_claudia_scanner_10mm" +simset_material: 29 +average_doi: 0.285 #0.285 # 0.6 +scanner_radius: 13.4 +num_rings: 48 +axial_fov: 16.318 +z_crystal_size: 0.319 +transaxial_crystal_size: 0.319 +crystal_thickness: 1 #*modif de 1 +energy_resolution: 14 #revizar por que nosotros entramos un valor y reporta otro en stir_sinogram.hs :Energy_resolution := 0.0273972602739726 +num_aa_bins: 128 +num_td_bins: 256 +min_energy_window: 425 +max_energy_window: 650 +coincidence_window: 5 +numberOfSubsets: 8 +numberOfIterations: 48 +savingInterval: 16 +analytical_att_correction: 0 #1 esta corrigiendo por atenuación; 0 no esta corrigiendo. Si stir_recons_att_corr: 1 entonces analytical_att_correction debe ser 0 +stir_recons_att_corr: 1 # mas novedoso #1 esta corrigiendo por atenuación; 0 no esta corrigiendo. Si analytical_att_correction:1 entonces stir_recons_att_corr debe ser 0 +analytic_scatt_corr_factor: 0.15 #falta por revisar +stir_scatt_corr_smoothing: 0 #falta por revisar +stir_scatt_simulation: 0 #falta por revisar +analytic_randoms_corr_factor: 0.2 #falta por revisar +stir_randoms_corr_smoothing: 0 #falta por revisar +recons_type: "OSEM3D" #OSEM3D +inter_iteration_filter: 1 #1 se está filtrando +subiteration_interval: 4 +x_dir_filter_FWHM: 1.45 # 1.30 #valor inicial 1.5 #rango aceptable 1.20 a 1.45 +y_dir_filter_FWHM: 1.45 # 1.30 #valor inicial 1.5 #rango aceptable 1.20 a 1.45 +z_dir_filter_FWHM: 0.60 # 0.5 #valor inicial 1.5 +psf_value: 0 +add_noise: 0 +max_segment: 47 #47 3D #1 PARA 2D +zoomFactor: 1.264 #1.264 +xyOutputSize: 325 +zOutputSize: 96 +zOutputVoxelSize: 1.64 #No hace nada diff --git a/configs/params/scanner/Neuro_LF/NeuroLF_scanner_20mm.yaml b/configs/params/scanner/Neuro_LF/NeuroLF_scanner_20mm.yaml new file mode 100644 index 00000000..be67518c --- /dev/null +++ b/configs/params/scanner/Neuro_LF/NeuroLF_scanner_20mm.yaml @@ -0,0 +1,38 @@ +scanner_name: "NeuroLF_BPET_claudia_scanner_20mm" +simset_material: 29 +average_doi: 0.285 #0.285 # 0.6 +scanner_radius: 13.4 +num_rings: 48 +axial_fov: 16.318 +z_crystal_size: 0.319 +transaxial_crystal_size: 0.319 +crystal_thickness: 2 #*modif de 1 +energy_resolution: 14 #revizar por que nosotros entramos un valor y reporta otro en stir_sinogram.hs :Energy_resolution := 0.0273972602739726 +num_aa_bins: 128 +num_td_bins: 256 +min_energy_window: 425 +max_energy_window: 650 +coincidence_window: 5 +numberOfSubsets: 8 +numberOfIterations: 48 +savingInterval: 16 +analytical_att_correction: 0 #1 esta corrigiendo por atenuación; 0 no esta corrigiendo. Si stir_recons_att_corr: 1 entonces analytical_att_correction debe ser 0 +stir_recons_att_corr: 1 # mas novedoso #1 esta corrigiendo por atenuación; 0 no esta corrigiendo. Si analytical_att_correction:1 entonces stir_recons_att_corr debe ser 0 +analytic_scatt_corr_factor: 0.15 #falta por revisar +stir_scatt_corr_smoothing: 0 #falta por revisar +stir_scatt_simulation: 0 #falta por revisar +analytic_randoms_corr_factor: 0.2 #falta por revisar +stir_randoms_corr_smoothing: 0 #falta por revisar +recons_type: "OSEM3D" #OSEM3D +inter_iteration_filter: 1 #1 se está filtrando +subiteration_interval: 4 +x_dir_filter_FWHM: 1.45 # 1.30 #valor inicial 1.5 #rango aceptable 1.20 a 1.45 +y_dir_filter_FWHM: 1.45 # 1.30 #valor inicial 1.5 #rango aceptable 1.20 a 1.45 +z_dir_filter_FWHM: 0.60 # 0.5 #valor inicial 1.5 +psf_value: 0 +add_noise: 0 +max_segment: 47 #47 3D #1 PARA 2D +zoomFactor: 1.264 #1.264 +xyOutputSize: 325 +zOutputSize: 96 +zOutputVoxelSize: 1.64 #No hace nada diff --git a/configs/params/scanner/Neuro_LF/NeuroLF_scanner_25mm.yaml b/configs/params/scanner/Neuro_LF/NeuroLF_scanner_25mm.yaml new file mode 100644 index 00000000..58cfd05f --- /dev/null +++ b/configs/params/scanner/Neuro_LF/NeuroLF_scanner_25mm.yaml @@ -0,0 +1,38 @@ +scanner_name: "NeuroLF_BPET_claudia_scanner_25mm" +simset_material: 29 +average_doi: 0.285 #0.285 # 0.6 +scanner_radius: 13.4 +num_rings: 48 +axial_fov: 16.318 +z_crystal_size: 0.319 +transaxial_crystal_size: 0.319 +crystal_thickness: 2.5 #*modif de 1 +energy_resolution: 14 #revizar por que nosotros entramos un valor y reporta otro en stir_sinogram.hs :Energy_resolution := 0.0273972602739726 +num_aa_bins: 128 +num_td_bins: 256 +min_energy_window: 425 +max_energy_window: 650 +coincidence_window: 5 +numberOfSubsets: 8 +numberOfIterations: 48 +savingInterval: 16 +analytical_att_correction: 0 #1 esta corrigiendo por atenuación; 0 no esta corrigiendo. Si stir_recons_att_corr: 1 entonces analytical_att_correction debe ser 0 +stir_recons_att_corr: 1 # mas novedoso #1 esta corrigiendo por atenuación; 0 no esta corrigiendo. Si analytical_att_correction:1 entonces stir_recons_att_corr debe ser 0 +analytic_scatt_corr_factor: 0.15 #falta por revisar +stir_scatt_corr_smoothing: 0 #falta por revisar +stir_scatt_simulation: 0 #falta por revisar +analytic_randoms_corr_factor: 0.2 #falta por revisar +stir_randoms_corr_smoothing: 0 #falta por revisar +recons_type: "OSEM3D" #OSEM3D +inter_iteration_filter: 1 #1 se está filtrando +subiteration_interval: 4 +x_dir_filter_FWHM: 1.45 # 1.30 #valor inicial 1.5 #rango aceptable 1.20 a 1.45 +y_dir_filter_FWHM: 1.45 # 1.30 #valor inicial 1.5 #rango aceptable 1.20 a 1.45 +z_dir_filter_FWHM: 0.60 # 0.5 #valor inicial 1.5 +psf_value: 0 +add_noise: 0 +max_segment: 47 #47 3D #1 PARA 2D +zoomFactor: 1.264 #1.264 +xyOutputSize: 325 +zOutputSize: 96 +zOutputVoxelSize: 1.64 #No hace nada diff --git a/configs/params/scanner/Neuro_LF/NeuroLF_scanner_30mm.yaml b/configs/params/scanner/Neuro_LF/NeuroLF_scanner_30mm.yaml new file mode 100644 index 00000000..4ff00324 --- /dev/null +++ b/configs/params/scanner/Neuro_LF/NeuroLF_scanner_30mm.yaml @@ -0,0 +1,38 @@ +scanner_name: "NeuroLF_scanner_30mm" +simset_material: 29 +average_doi: 0.285 #0.285 # 0.6 +scanner_radius: 13.4 +num_rings: 48 +axial_fov: 16.318 +z_crystal_size: 0.319 +transaxial_crystal_size: 0.319 +crystal_thickness: 3 #*modif de 1 +energy_resolution: 14 #revizar por que nosotros entramos un valor y reporta otro en stir_sinogram.hs :Energy_resolution := 0.0273972602739726 +num_aa_bins: 128 +num_td_bins: 256 +min_energy_window: 425 +max_energy_window: 650 +coincidence_window: 5 +numberOfSubsets: 8 +numberOfIterations: 48 +savingInterval: 16 +analytical_att_correction: 0 #1 esta corrigiendo por atenuación; 0 no esta corrigiendo. Si stir_recons_att_corr: 1 entonces analytical_att_correction debe ser 0 +stir_recons_att_corr: 1 # mas novedoso #1 esta corrigiendo por atenuación; 0 no esta corrigiendo. Si analytical_att_correction:1 entonces stir_recons_att_corr debe ser 0 +analytic_scatt_corr_factor: 0.15 #falta por revisar +stir_scatt_corr_smoothing: 0 #falta por revisar +stir_scatt_simulation: 0 #falta por revisar +analytic_randoms_corr_factor: 0.2 #falta por revisar +stir_randoms_corr_smoothing: 0 #falta por revisar +recons_type: "FBP2D" #OSEM3D +inter_iteration_filter: 1 #1 se está filtrando +subiteration_interval: 4 +x_dir_filter_FWHM: 1.45 # 1.30 #valor inicial 1.5 #rango aceptable 1.20 a 1.45 +y_dir_filter_FWHM: 1.45 # 1.30 #valor inicial 1.5 #rango aceptable 1.20 a 1.45 +z_dir_filter_FWHM: 0.60 # 0.5 #valor inicial 1.5 +psf_value: 0 +add_noise: 0 +max_segment: 47 #47 3D #1 PARA 2D +zoomFactor: 1.264 #1.264 +xyOutputSize: 325 +zOutputSize: 96 +zOutputVoxelSize: 1.64 #No hace nada diff --git a/configs/params/scanner/Neuro_LF/NeuroLF_scanner_3mm_crystal_thickness.yaml b/configs/params/scanner/Neuro_LF/NeuroLF_scanner_3mm_crystal_thickness.yaml new file mode 100644 index 00000000..0cc46943 --- /dev/null +++ b/configs/params/scanner/Neuro_LF/NeuroLF_scanner_3mm_crystal_thickness.yaml @@ -0,0 +1,38 @@ +scanner_name: "NeuroLF_BPET_claudia_scanner_3mm_crystal_thickness" +simset_material: 29 +average_doi: 0.285 +scanner_radius: 13.4 +num_rings: 48 +axial_fov: 16.318 +z_crystal_size: 0.319 +transaxial_crystal_size: 0.319 +crystal_thickness: 0.3 +energy_resolution: 14 #revizar por que nosotros entramos un valor y reporta otro en stir_sinogram.hs :Energy_resolution := 0.0273972602739726 +num_aa_bins: 128 +num_td_bins: 256 +min_energy_window: 425 +max_energy_window: 650 +coincidence_window: 5 +numberOfSubsets: 8 +numberOfIterations: 48 +savingInterval: 16 +analytical_att_correction: 0 #1 esta corrigiendo por atenuación; 0 no esta corrigiendo. Si stir_recons_att_corr: 1 entonces analytical_att_correction debe ser 0 +stir_recons_att_corr: 1 #1 # mas novedoso #1 esta corrigiendo por atenuación; 0 no esta corrigiendo. Si analytical_att_correction:1 entonces stir_recons_att_corr debe ser 0 +analytic_scatt_corr_factor: 0.15 #falta por revisar +stir_scatt_corr_smoothing: 0 #falta por revisar +stir_scatt_simulation: 0 #falta por revisar +analytic_randoms_corr_factor: 0.2 #falta por revisar +stir_randoms_corr_smoothing: 0 #falta por revisar +recons_type: "OSEM3D" #OSEM3D +inter_iteration_filter: 1 #1 se está filtrando +subiteration_interval: 4 +x_dir_filter_FWHM: 1.45 # 1.30 #valor inicial 1.5 #rango aceptable 1.20 a 1.45 +y_dir_filter_FWHM: 1.45 # 1.30 #valor inicial 1.5 #rango aceptable 1.20 a 1.45 +z_dir_filter_FWHM: 0.60 # 0.5 #valor inicial 1.5 +psf_value: 0 +add_noise: 0 +max_segment: 47 #47 3D #1 PARA 2D +zoomFactor: 1 #1.264 +xyOutputSize: 260 #325 +zOutputSize: 96 +zOutputVoxelSize: 1.64 #No hace nada diff --git a/configs/params/scanner/Neuro_LF/NeuroLF_scanner_8bloques_axial_FOV.yaml b/configs/params/scanner/Neuro_LF/NeuroLF_scanner_8bloques_axial_FOV.yaml new file mode 100644 index 00000000..bbe00bbc --- /dev/null +++ b/configs/params/scanner/Neuro_LF/NeuroLF_scanner_8bloques_axial_FOV.yaml @@ -0,0 +1,38 @@ +scanner_name: "NeuroLF_scanner_8bloques_axial_FOV" +simset_material: 29 +average_doi: 0.285 +scanner_radius: 13.4 +num_rings: 64 +axial_fov: 20.5 +z_crystal_size: 0.319 +transaxial_crystal_size: 0.319 +crystal_thickness: 1 +energy_resolution: 14 #revizar por que nosotros entramos un valor y reporta otro en stir_sinogram.hs :Energy_resolution := 0.0273972602739726 +num_aa_bins: 128 +num_td_bins: 256 +min_energy_window: 425 +max_energy_window: 650 +coincidence_window: 5 +numberOfSubsets: 8 +numberOfIterations: 48 +savingInterval: 16 +analytical_att_correction: 0 #1 esta corrigiendo por atenuación; 0 no esta corrigiendo. Si stir_recons_att_corr: 1 entonces analytical_att_correction debe ser 0 +stir_recons_att_corr: 1 # mas novedoso #1 esta corrigiendo por atenuación; 0 no esta corrigiendo. Si analytical_att_correction:1 entonces stir_recons_att_corr debe ser 0 +analytic_scatt_corr_factor: 0.15 #falta por revisar +stir_scatt_corr_smoothing: 0 #falta por revisar +stir_scatt_simulation: 0 #falta por revisar +analytic_randoms_corr_factor: 0.2 #falta por revisar +stir_randoms_corr_smoothing: 0 #falta por revisar +recons_type: "OSEM3D" #OSEM3D +inter_iteration_filter: 1 #1 se está filtrando +subiteration_interval: 4 +x_dir_filter_FWHM: 1.45 # 1.30 #valor inicial 1.5 #rango aceptable 1.20 a 1.45 +y_dir_filter_FWHM: 1.45 # 1.30 #valor inicial 1.5 #rango aceptable 1.20 a 1.45 +z_dir_filter_FWHM: 0.60 # 0.5 #valor inicial 1.5 +psf_value: 0 +add_noise: 0 +max_segment: 63 #47 3D #1 PARA 2D +zoomFactor: 1.264 #1.264 +xyOutputSize: 325 +zOutputSize: 100 +zOutputVoxelSize: 1.64 #No hace nada diff --git a/configs/params/scanner/Neuro_LF/NeuroLF_scanner_cristal_2mm.yaml b/configs/params/scanner/Neuro_LF/NeuroLF_scanner_cristal_2mm.yaml new file mode 100644 index 00000000..18338325 --- /dev/null +++ b/configs/params/scanner/Neuro_LF/NeuroLF_scanner_cristal_2mm.yaml @@ -0,0 +1,38 @@ +scanner_name: "NeuroLF_BPET_claudia_scanner_cristal_2mm" +simset_material: 29 +average_doi: 0.285 +scanner_radius: 13.4 +num_rings: 76 +axial_fov: 16.318 +z_crystal_size: 0.2 +transaxial_crystal_size: 0.2 +crystal_thickness: 2.5 +energy_resolution: 14 #revizar por que nosotros entramos un valor y reporta otro en stir_sinogram.hs :Energy_resolution := 0.0273972602739726 +num_aa_bins: 128 #128 +num_td_bins: 360 #256 +bin_energy_window: 425 +max_energy_window: 650 +coincidence_window: 5 +numberOfSubsets: 8 +numberOfIterations: 48 +savingInterval: 16 +analytical_att_correction: 0 #1 esta corrigiendo por atenuación; 0 no esta corrigiendo. Si stir_recons_att_corr: 1 entonces analytical_att_correction debe ser 0 +stir_recons_att_corr: 1 #1 #mas novedoso,1 esta corrigiendo por atenuación; 0 no esta corrigiendo. Si analytical_att_correction:1 entonces stir_recons_att_corr debe ser 0 +analytic_scatt_corr_factor: 0.15 #falta por revisar +stir_scatt_corr_smoothing: 0 #falta por revisar +stir_scatt_simulation: 0 #falta por revisar +analytic_randoms_corr_factor: 0.2 #falta por revisar +stir_randoms_corr_smoothing: 0 #falta por revisar +recons_type: "OSEM3D" #OSEM3D +inter_iteration_filter: 1 #1 se está filtrando +subiteration_interval: 4 +x_dir_filter_FWHM: 1.45 # 1.30 #valor inicial 1.5 #rango aceptable 1.20 a 1.45 +y_dir_filter_FWHM: 1.45 # 1.30 #valor inicial 1.5 #rango aceptable 1.20 a 1.45 +z_dir_filter_FWHM: 0.60 # 0.5 #valor inicial 1.5 +psf_value: 0 +add_noise: 0 +max_segment: 75 # 75 3D #1 PARA 2D +zoomFactor: 1.37 #1.264 +xyOutputSize: 498 +zOutputSize: 152 #96 +zOutputVoxelSize: 1.64 #No hace nada diff --git a/configs/params/scanner/Neuro_LF/NeuroLF_scanner_final.yaml b/configs/params/scanner/Neuro_LF/NeuroLF_scanner_final.yaml new file mode 100644 index 00000000..560d157c --- /dev/null +++ b/configs/params/scanner/Neuro_LF/NeuroLF_scanner_final.yaml @@ -0,0 +1,38 @@ +scanner_name: "NeuroLF_scanner_final" +simset_material: 29 +average_doi: 0.285 +scanner_radius: 13.4 +num_rings: 48 +axial_fov: 16.318 #16.318 original +z_crystal_size: 0.319 +transaxial_crystal_size: 0.319 +crystal_thickness: 2.5 #2.5 +energy_resolution: 14 #revizar por que nosotros entramos un valor y reporta otro en stir_sinogram.hs :Energy_resolution := 0.0273972602739726 +num_aa_bins: 128 +num_td_bins: 256 +min_energy_window: 425 +max_energy_window: 650 +coincidence_window: 5 +numberOfSubsets: 8 +numberOfIterations: 48 +savingInterval: 16 +analytical_att_correction: 0 #1 esta corrigiendo por atenuación; 0 no esta corrigiendo. Si stir_recons_att_corr: 1 entonces analytical_att_correction debe ser 0 +stir_recons_att_corr: 1 #1 # mas novedoso #1 esta corrigiendo por atenuación; 0 no esta corrigiendo. Si analytical_att_correction:1 entonces stir_recons_att_corr debe ser 0 +analytic_scatt_corr_factor: 0.15 #falta por revisar +stir_scatt_corr_smoothing: 0 #falta por revisar +stir_scatt_simulation: 0 #falta por revisar +analytic_randoms_corr_factor: 0.2 #falta por revisar +stir_randoms_corr_smoothing: 0 #falta por revisar +recons_type: "OSEM3D" #OSEM3D +inter_iteration_filter: 1 #1 se está filtrando +subiteration_interval: 4 +x_dir_filter_FWHM: 1.45 # 1.30 #valor inicial 1.5 #rango aceptable 1.20 a 1.45 +y_dir_filter_FWHM: 1.45 # 1.30 #valor inicial 1.5 #rango aceptable 1.20 a 1.45 +z_dir_filter_FWHM: 0.60 # 0.5 #valor inicial 1.5 +psf_value: 0 +add_noise: 0 +max_segment: 47 #47 3D #1 PARA 2D +zoomFactor: 1.264 #1.264 +xyOutputSize: 325 +zOutputSize: 100 +zOutputVoxelSize: 1.64 #No hace nada diff --git a/configs/params/scanner/Vereos_PET_CT_scanner.yaml b/configs/params/scanner/Vereos_PET_CT_scanner.yaml new file mode 100644 index 00000000..b93b56f3 --- /dev/null +++ b/configs/params/scanner/Vereos_PET_CT_scanner.yaml @@ -0,0 +1,47 @@ +scanner_name: "Vereos_PET_CT_scanner" +simset_material: 29 +average_doi: 0.84 +scanner_radius: 38.2 +num_rings: 40 +axial_fov: 16.4 +z_crystal_size: 0.4 +transaxial_crystal_size: 0.4 +crystal_thickness: 1.9 +energy_resolution: 11.2 +num_aa_bins: 297 +num_td_bins: 333 +min_energy_window: 450 +max_energy_window: 613 +coincidence_window: 4 +numberOfSubsets: 11 +numberOfIterations: 33 #48 +savingInterval: 11 +analytical_att_correction: 1 #cambiar a 0 +stir_recons_att_corr: 0 #cambiar a 1 +analytic_scatt_corr_factor: 0.15 +stir_scatt_corr_smoothing: 0 +stir_scatt_simulation: 0 +analytic_randoms_corr_factor: 0.2 +stir_randoms_corr_smoothing: 0 +recons_type: "OSEM3D" +inter_iteration_filter: 0 # 0 no se está filtrando +subiteration_interval: 48 +x_dir_filter_FWHM: 0.7 +y_dir_filter_FWHM: 1 +z_dir_filter_FWHM: 0 +psf_value: 0 #Hace una convolución con psf para cuantificar la incertidumbre de efectos que no simulamos (efectos dentro del cristal). La ponemos a 0 para estudiar sin esta corrección. +add_noise: 0 #Como en el anterior filtrado se elimina el ruido, añadimos este analíticamente. +max_segment: 39 +zoomFactor: 1 # 0.795 para pixel de 0.75mm +xyOutputSize: 231 +zOutputSize: 79 +zOutputVoxelSize: 3.27 # Esto no se usa per ahora por defecto es la mitad del tamaño del cristal o sea 0.22/2 = 1.1 mm + +# NOTAS + +# El pixel en z corresponde a la mitad del tamaño del cristal. +# El pixel en x,y viene asociado al zoom factor. + +# bin_size=2*scanner_radius/num_td_bins +# xyVoxelsize=10*bin_size/zoom +# zVoxelsize=z_crystal_size/2 diff --git a/configs/params/scanner/discovery_st.yaml b/configs/params/scanner/discovery_st.yaml index e7f02dfd..3e0fb533 100644 --- a/configs/params/scanner/discovery_st.yaml +++ b/configs/params/scanner/discovery_st.yaml @@ -34,5 +34,5 @@ add_noise: 0 max_segment: 23 zoomFactor: 1.55 xyOutputSize: 128 -zOutputSize: 47 +zOutputSize: 160 #47 zOutputVoxelSize: 3.27 diff --git a/configs/params/scanner/ls b/configs/params/scanner/ls new file mode 100644 index 00000000..e69de29b diff --git a/configs/params/test.yaml b/configs/params/test.yaml index f47b6f29..2af9e7ad 100644 --- a/configs/params/test.yaml +++ b/configs/params/test.yaml @@ -10,7 +10,7 @@ model_type: "cylindrical" patient_dirname: "test_image" act_map: "act.hdr" att_map: "att.hdr" -output_dir: "test_image" +output_dir: "test_image_160slices" center_slice: 7 total_dose: 0.1 simulation_time: 30 diff --git a/configs/params/test_wholebody_reference.yaml b/configs/params/test_wholebody_reference.yaml new file mode 100644 index 00000000..f50f1b45 --- /dev/null +++ b/configs/params/test_wholebody_reference.yaml @@ -0,0 +1,38 @@ +defaults: + - scanner: Bruker_PET_MRI_scanner + +simulation_environment: 0 +sim_type: "SimSET" +do_simulation: 1 +do_reconstruction: 1 +divisions: 2 +model_type: "cylindrical" +patient_dirname: "Rat_Phantom_Phantech" +mask_dirname: "Rat_Phantom_Phantech" +act_map: "PET_study_9_real.hdr" +att_map: "RatPhantom_Phantech_att.hdr" +mask_map: "Phantom_Phantech_Reference_mask.hdr" +output_dir: "Rat_Phantom_Phantech_simulation" +center_slice: 0 +total_dose: 0.03592 +simulation_time: 600 +sampling_photons: 2000000 +photons: 0 +add_randoms: 0 +phglistmode: 0 +detlistmode: 0 +maximumIteration: 1 + +#Corrections: +correction_for_NECR: 1 # 0-NO 1-SI +correction_for_normalization: 1 # 0-No 1-Si +correction_total_fov: 1 # 0-No 1-Si + +#Quantification +quantification_info: 1 # 0-NO 1-SI + +# Variables needed for whole body simulation +whole_body_simulation : 1 # 0-NO 1-SI +z_min: 0 +z_max: 730 +joints_beds: 1 # 0-NO 1-SI diff --git a/makefile b/makefile deleted file mode 100644 index 8080b1c9..00000000 --- a/makefile +++ /dev/null @@ -1,256 +0,0 @@ -SHELL := /bin/bash -TMPDIR := /tmp -ROOT_DIR := $(dir $(abspath $(lastword $(MAKEFILE_LIST)))) -ASSETS_DIR := ${ROOT_DIR}assets -INCLUDE_DIR := ${ROOT_DIR}include -SUBMODULES_DIR := ${ROOT_DIR}submodules - -.PHONY: deps install-simset check-simset clean-simset -.PHONY: install-stir check-stir clean-stir -.PHONY: install-resources check-resources clean-resources -.PHONY: config-git clean-git smudge-filer config-paths clean-paths -.PHONY: dummy-data install clean help - -deps: - sudo apt-get -y -q update ;\ - sudo apt-get install -y -q \ - software-properties-common \ - git[all] \ - git-lfs \ - wget \ - unzip \ - sshpass \ - libboost-dev \ - libboost-all-dev \ - libpcre3 \ - libpcre3-dev \ - libncurses-dev \ - cmake \ - g++ \ - swig ;\ - sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-keys CC86BB64 ;\ - sudo add-apt-repository -y ppa:rmescandon/yq ;\ - sudo apt-get -y -q update ;\ - sudo apt-get install -y -q yq - -SIMSET_SUBMODULE_DIR = ${ROOT_DIR}/submodules/simset_simpet -SIMSET_DEST_DIR = ${INCLUDE_DIR}/SimSET -SIMSET_PATH = ${SIMSET_DEST_DIR}/2.9.2 -SIMSET_BIN = ${SIMSET_PATH}/bin -SIMSET_LIB = ${SIMSET_PATH}/lib -SIMSET_MKALL = ${SIMSET_PATH}/make_all.sh -SIMSET_MKFILE = ${SIMSET_PATH}/make.files/simset.make - -install-simset: - if [ ! -d ${SIMSET_DEST_DIR} ]; then \ - mkdir -p ${SIMSET_PATH} ;\ - cp -r ${SIMSET_SUBMODULE_DIR}/* ${SIMSET_PATH} ;\ - sed -i 's/^\(SIMSET_PATH = \).*$$/\1$(subst /,\/,${SIMSET_PATH})/' ${SIMSET_MKFILE} ;\ - cd ${SIMSET_PATH} && mkdir -p ${SIMSET_LIB} && bash ${SIMSET_MKALL} ;\ - else \ - echo "${SIMSET_DEST_DIR} already exists, run clean-simset if you really want to remove it (you will have to intall SimSET again.)" ;\ - fi - -SIMSET_TAR = ${ASSETS_DIR}/phg.2.9.2.tar.Z -SIMSET_STIR_PATCH = ${ASSETS_DIR}/simset_for_stir.patch - -install-canonic-simset: ${SIMSET_TAR} - if [ ! -d ${SIMSET_DEST_DIR} ]; then \ - mkdir -p ${SIMSET_DEST_DIR} && tar -xvf "${SIMSET_TAR}" --directory=${SIMSET_DEST_DIR} ;\ - cd ${SIMSET_DEST_DIR} && patch -s -p0 < ${SIMSET_STIR_PATCH} ;\ - sed -i 's/^\(SIMSET_PATH = \).*$$/\1$(subst /,\/,${SIMSET_PATH})/' ${SIMSET_MKFILE} ;\ - cd ${SIMSET_PATH} && mkdir -p ${SIMSET_LIB} && bash ${SIMSET_MKALL} ;\ - else \ - echo "${SIMSET_DEST_DIR} already exists, run clean-simset if you really want to remove it (you will have to intall SimSET again.)" ;\ - fi - -check-simset: - declare -a simset_files=(addrandoms bin calcattenuation combinehist makeindexfile phg timesort) ;\ - for file in "$${simset_files[@]}"; do \ - if [ ! -f "${SIMSET_BIN}/$${file}" ]; then \ - echo "${SIMSET_BIN}/$${file} does not exists, check your installation." ;\ - else \ - echo "${SIMSET_BIN}/$${file} exists." ;\ - fi ;\ - done - -clean-simset: - rm -rf ${SIMSET_DEST_DIR} - -SIMSET_SRC = ${SIMSET_PATH}/src -SIMSET_LIBSIMSET = ${SIMSET_LIB}/libsimset.so -STIR_DIR = ${SUBMODULES_DIR}/STIR/STIR -STIR_BUILD_DIR = ${SUBMODULES_DIR}/STIR/build -STIR_INSTALL_DIR = ${SUBMODULES_DIR}/STIR/install -STIR_INSTALL_BIN = ${STIR_INSTALL_DIR}/bin -STIR_MKFILE = ${STIR_BUILD_DIR}/CMakeCache.txt -STIR_FINAL_DEST_DIR = ${INCLUDE_DIR}/STIR -NPROC = $(shell nproc) - -install-stir: install-simset ${STIR_DIR} - if [ ! -d ${STIR_FINAL_DEST_DIR} ]; then \ - mkdir -p ${STIR_BUILD_DIR} ${STIR_INSTALL_DIR} ${STIR_FINAL_DEST_DIR} ;\ - cd ${STIR_BUILD_DIR} && cmake ${STIR_DIR} ;\ - sed -i \ - -e 's/^\(BUILD_SWIG_PYTHON\).*$$/\1:BOOL=OFF/' \ - -e 's/^\(CMAKE_INSTALL_PREFIX\).*$$/\1:PATH=$(subst /,\/,${STIR_INSTALL_DIR})/' \ - -e 's/^\(SIMSET_INCLUDE_DIRS\).*$$/\1:PATH=$(subst /,\/,${SIMSET_SRC})/' \ - -e 's/^\(SIMSET_LIBRARY\).*$$/\1:FILEPATH=$(subst /,\/,${SIMSET_LIBSIMSET})/' \ - -e 's/^\(STIR_OPENMP\).*$$/\1:BOOL=ON/' \ - ${STIR_MKFILE} ;\ - cmake ${STIR_DIR} ;\ - make -s -j${NPROC} ;\ - sudo make install ;\ - mv ${STIR_BUILD_DIR} ${STIR_INSTALL_DIR} ${STIR_FINAL_DEST_DIR} ;\ - else \ - echo "${STIR_FINAL_DEST_DIR} already exists, run clean-stir if you really want to remove it (you will have to intall STIR again)." ;\ - fi - -check-stir: - declare -a stir_files=(FBP2D FBP3DRP forward_project lm_to_projdata OSMAPOSL zoom_image) ;\ - for file in "$${stir_files[@]}"; do \ - if [ ! -f "${STIR_FINAL_DEST_DIR}/install/bin/$${file}" ]; then \ - echo "${STIR_FINAL_DEST_DIR}/install/bin/$${file} does not exists, check your installation." ;\ - else \ - echo "${STIR_FINAL_DEST_DIR}/install/bin/$${file} exists." ;\ - fi ;\ - done - -clean-stir: - rm -rf ${STIR_FINAL_DEST_DIR} - -RESOURCES_ZIP = ${ASSETS_DIR}/fruitcake.zip -RESOURCES_TMP = ${TMPDIR}/resources - -install-resources: ${RESOURCES_ZIP} - mkdir -p ${RESOURCES_TMP} ${INCLUDE_DIR} && unzip -o ${RESOURCES_ZIP} -d ${RESOURCES_TMP} ;\ - declare -a resources=(fruitcake format_converters) ;\ - for rce in "$${resources[@]}"; do \ - include_path="${INCLUDE_DIR}/$${rce}" ;\ - tmp_path="${RESOURCES_TMP}/$${rce}" ;\ - if [ ! -d $${include_path} ]; then \ - mv $${tmp_path} ${INCLUDE_DIR} ;\ - chmod -R +x $${include_path} ;\ - else \ - echo "$${include_path} already exists, run make clean-resources in order to clean resources installation (fruitcake will be removed as well)." ;\ - fi ;\ - done ;\ - rm -rf ${RESOURCES_TMP} - -check-resources: - declare -a resources_paths=(${FRUITCAKE_PATH} ${FORMAT_CONVERTERS_PATH}) ;\ - for path in "$${resources_paths[@]}"; do \ - if [ ! -d "$${path}" ]; then \ - echo "$${path} does not exists, clean fruitcake and format_converters with clean-resources and run install-resources." ;\ - else \ - echo "$${path} exists." ;\ - fi ;\ - done - -clean-resources: - rm -rf ${INCLUDE_DIR}/fruitcake ${INCLUDE_DIR}/format_converters - -config-git: - git config --local filter.config-test.smudge ${ROOT_DIR}scripts/smudge-config_test.sh - git config --local filter.config-test.clean ${ROOT_DIR}scripts/clean-config_test.sh - chmod +x -R ${ROOT_DIR}scripts - -clean-git: - git config --local --unset filter.config-test.smudge - git config --local --unset filter.config-test.clean - -smudge-filter: - rm -f ${ROOT_DIR}configs/config.yaml 2> /dev/null - git checkout HEAD -- ${ROOT_DIR}configs/config.yaml - -FRUITCAKE_PATH = ${INCLUDE_DIR}/fruitcake -FRUITCAKE_BIN = ${FRUITCAKE_PATH}/bin -FRUITCAKE_LIB = ${FRUITCAKE_PATH}/book/lib -FORMAT_CONVERTERS_PATH = ${INCLUDE_DIR}/format_converters - -config-paths: - touch "$${HOME}"/.bashrc ;\ - declare -a simpet_paths=( \ - 'export PATH=${FRUITCAKE_BIN}:$$PATH' \ - 'export LD_LIBRARY_PATH=${FRUITCAKE_LIB}:$$LD_LIBRARY_PATH' \ - 'export PATH=${FORMAT_CONVERTERS_PATH}:$$PATH' \ - ) ;\ - for path in "$${simpet_paths[@]}"; do \ - grep -qxF "$${path}" "$${HOME}"/.bashrc || echo "$${path}" >> "$${HOME}"/.bashrc ;\ - done - -clean-paths: - touch "$${HOME}/.bashrc" - sed -i \ - -e '/export PATH=$(subst /,\/,${FRUITCAKE_BIN}):$$PATH/d' \ - -e '/export LD_LIBRARY_PATH=$(subst /,\/,${FRUITCAKE_LIB}):$$LD_LIBRARY_PATH/d' \ - -e '/export PATH=$(subst /,\/,${FORMAT_CONVERTERS_PATH}):$$PATH/d' \ - "$${HOME}/.bashrc" - -DATA_DIR = ${ROOT_DIR}Data -DATA_ZIP = ${ASSETS_DIR}/Data.zip - -dummy-data: ${DATA_ZIP} - if [ ! -d "${DATA_DIR}" ]; then \ - mkdir -p ${DATA_DIR} && unzip -o ${DATA_ZIP} -d ${ROOT_DIR} ;\ - else \ - echo "${DATA_DIR} already exists, remove it manually." ;\ - fi - -install: - ${MAKE} \ - deps \ - install-simset \ - check-simset \ - install-stir \ - check-stir \ - install-resources \ - check-resources \ - config-paths \ - -clean: - ${MAKE} \ - clean-simset \ - clean-stir \ - clean-resources - -help: - @echo "Help:" - @echo " - deps: Install the dependencies of the projects via apt." - @echo "" - @echo " - install-simset: Install tweaked version of SimSET. - @echo "" - @echo " - install-canonic-simset: Install non-tweaked version of SimSET with STIR patch at directory ${SIMSET_DEST_DIR}." - @echo "" - @echo " - check-simset: Check that SimSET binaries exist." - @echo "" - @echo " - clean-simset: Clean SimSET installation (removes ${SIMSET_DEST_DIR} directory)." - @echo "" - @echo " - install-stir: Install STIR at directory $(shell dirname ${STIR_DIR}), install-simset is a prerequisite. Set NPROC=n to use n CPU cores in compilation." - @echo "" - @echo " - check-stir: Check that STIR binaries exist." - @echo "" - @echo " - clean-stir: Clean STIR installation removing ${STIR_INSTALL_DIR} and ${STIR_BUILD_DIR} directories." - @echo "" - @echo " - install-resources: Decompress ${RESOURCES_ZIP} and move fruitcake and format_converters to ${INCLUDE_DIR}" - @echo "" - @echo " - check-resources: Checks that fruitcake and format_converters are in ${INCLUDE_DIR}." - @echo "" - @echo " - clean-resources: Remove fruitcake and format_converters from ${INCLUDE_DIR}." - @echo "" - @echo " - config-git: Add project filter drivers to local git configuration and make them executable." - @echo "" - @echo " - clean-git: Remove project filter drivers from local git configuration." - @echo "" - @echo " - smudge-filer: Apply smudge filter to config file (updates paths)." - @echo "" - @echo " - config-paths: Iff not present in .bashrc, the paths of the projects will be appended to the file. If ~/.bashrc does not exists it will be created." - @echo "" - @echo " - clean-paths: If the paths of the project are present in .bashrc file, they will be deleted. If ~/.bashrc does not exists it will be created." - @echo "" - @echo " - dummy-data: Uncompress ${DATA_ZIP} at ${DATA_DIR} if ${DATA_DIR} does not exists." - @echo "" - @echo " - install: Run all install, config and check recipes." - @echo "" - @echo " - clean: Run all clean recipes." - diff --git a/requirements.txt b/requirements.txt index 40fb7a42..5dbd8935 100644 --- a/requirements.txt +++ b/requirements.txt @@ -9,5 +9,5 @@ nilearn==0.10.0 pexpect==4.8.0 hydra-core==1.3 pyprojroot==0.2.0 -black==26.3.1 +black==24.3.0 diff --git a/scripts/experiment_wholebody.py b/scripts/experiment_wholebody.py index 1cd1b8fb..a8c78fd9 100644 --- a/scripts/experiment_wholebody.py +++ b/scripts/experiment_wholebody.py @@ -7,7 +7,6 @@ sys.path.append(str(here())) import wholebody - try: OmegaConf.register_new_resolver("root_path", lambda s: str(here())) OmegaConf.register_new_resolver("home_path", lambda s: str(Path.home())) @@ -20,7 +19,7 @@ def wholebody_simulation(cfg: DictConfig) -> None: OmegaConf.resolve(cfg) wb = wholebody.WholebodySimulation(cfg) wb.run() - + if __name__ == "__main__": wholebody_simulation() diff --git a/src/simset/simset_sim.py b/src/simset/simset_sim.py index f06ddfcc..25f4cf87 100644 --- a/src/simset/simset_sim.py +++ b/src/simset/simset_sim.py @@ -6,12 +6,14 @@ import numpy as np import nibabel as nib import warnings +import math #añadido para Bruker PET/MRI PRECLINICA from multiprocessing import Process from pathlib import Path from os import PathLike from os.path import join, dirname, abspath, exists import src.simset.simset_tools as simset_tools from utils import tools +from utils import wb_tools def read_ws_from_simset_log(simset_log: PathLike) -> float: @@ -76,11 +78,16 @@ def __init__( self.s_photons = params.get("sampling_photons") self.photons = params.get("photons") self.sim_time = params.get("simulation_time") + self.divisions = params.get("divisions") self.detlistmode = params.get("detlistmode") self.phglistmode = params.get("phglistmode") self.add_randoms = params.get("add_randoms") + + #TODO #correcton remove to the code + + def run(self): processes = [] @@ -90,10 +97,10 @@ def run(self): print('Scanner: %s' % self.scanner.get("scanner_name")) print('Activity map: %s' % self.act_map) print('Attenuation map: %s' % self.att_map) - print('Dose: %s mCi' % self.sim_dose) + print('Total Dose: %s mCi' % self.sim_dose) print('Acquisition time: %s seconds' % self.sim_time) print('------------------------------------------------------------') - + for division in range(self.divisions): division_dir = join(self.output_dir, "division_" + str(division)) os.makedirs(division_dir) @@ -137,6 +144,28 @@ def run_simset_simulation(self, sim_dir): act_table_factor = self.sim_dose * 1000 / phantom_dose else: act_table_factor = 1 + + #Borrer a partir de aqui + print(f"act_table_factor_value:{act_table_factor}") + + img = nib.load(self.act_map) + data = img.get_fdata() + # Obtener los labels únicos (excluyendo el 0 si es fondo) + labels, counts = np.unique(data, return_counts=True) + + # Eliminar el fondo (label 0) si lo hay + mask = labels != 0 + labels = labels[mask] + counts = counts[mask] + + # Mostrar resultado + for label, count in zip(labels, counts): + print(f"Label {int(label)}: {int(count)} voxeles") + + # Mostrar cantidad total de labels + print(f"\nTotal de labels distintos: {len(labels)}") + + #Hasta la linea de alante borrar # Creates the data files from the simulation maps act_img = self.act_map[0:-3] + "img" @@ -225,6 +254,7 @@ def run_simset_simulation(self, sim_dir): w_quotient = read_ws_from_simset_log(my_log) sim_photons = int(self.s_photons * w_quotient) + print(f"sim_photoms : {sim_photons}") else: # If the user stated photons, the provided value will be used @@ -551,9 +581,22 @@ def run_recons(self): print("Starting STIR reconstruction") recons_algorithm = self.scanner.get("recons_type") - sinogram_stir = join(self.output_dir, "stir_sinogram.hs") + + #sinogram_stir = join(self.output_dir, "stir_sinogram.hs") + #TODO correction for sinogram using whole FOV adquisition + if self.params.get("correction_total_fov") == 1: + recons_dir = self.output_dir + print(f"Recons Dir = {recons_dir}") + sinogram_stir_corrected = wb_tools.total_fov_correction(self, recons_dir) + sinogram_stir = join(self.output_dir, sinogram_stir_corrected) + else: + sinogram_stir = join(self.output_dir, "stir_sinogram.hs") + + additive_sino_stir = join(self.output_dir, "stir_additivesino.hs") att_stir = join(self.output_dir, "stir_att.hs") + + if any( exists(i) == False for i in [sinogram_stir, additive_sino_stir, att_stir] diff --git a/src/simset/simset_tool.py b/src/simset/simset_tool.py new file mode 100644 index 00000000..e69de29b diff --git a/src/simset/simset_tools.py b/src/simset/simset_tools.py index a92b6b7f..f9648d45 100644 --- a/src/simset/simset_tools.py +++ b/src/simset/simset_tools.py @@ -84,6 +84,9 @@ def make_simset_phg( xbins = act.shape[0] ybins = act.shape[1] + + + act_fov = abs(act.shape * act.affine) xMin, xMax = round(-act_fov[0, 0] / 20, 3), round( act_fov[0, 0] / 20, 3 @@ -95,7 +98,7 @@ def make_simset_phg( ) # cm zMin, zMax = round(-z_offset, 3), round(act_fov[2, 2] / 10 - z_offset, 3) - dz = round((zMax - zMin) / nslices, 2) + dz = round((zMax - zMin) / nslices, 3) max_z_target = scanner_axial_fov / 2 min_z_target = -scanner_axial_fov / 2 @@ -152,8 +155,8 @@ def make_simset_phg( f.write("\n INT num_slices = %s" % str(nslices)) for i in range(nslices): - zMin_value = round(zMin + i * dz, 2) - zMax_value = round(zMin + (i + 1) * dz, 2) + zMin_value = round(zMin + i * dz, 3) + zMax_value = round(zMin + (i + 1) * dz, 3) f.write( "\n NUM_ELEMENTS_IN_LIST slice = 9 " + "\n INT slice_number = %s" % str(i) diff --git a/src/stir/stir_sim.py b/src/stir/stir_sim.py index ae9c1745..6c60c7a5 100644 --- a/src/stir/stir_sim.py +++ b/src/stir/stir_sim.py @@ -33,6 +33,9 @@ def __init__( self.target_size = scanner.scanner_target_size self.scanner_template = scanner.scanner_template self.proyector = scanner.proyector + + #TODO borrar + print(f"Scanner_template:{self.scanner_template} \n Scanner_target_size:{self.target_size} \n Scanner_proyector: {self.proyector}") # Configuring required resources self.stir_dir = join(simpet_dir, "include", "stir", "bin") diff --git a/utils/Quantification_externo_Claudia.py b/utils/Quantification_externo_Claudia.py new file mode 100644 index 00000000..6227bdfc --- /dev/null +++ b/utils/Quantification_externo_Claudia.py @@ -0,0 +1,407 @@ + +import sys +import math #new +from os.path import join, dirname +import nibabel as nib +import numpy as np +import os +from scipy.ndimage import zoom +from scipy.ndimage import gaussian_filter +from scipy.ndimage import median_filter +from pyprojroot import here + +sys.path.append(str(here())) +from utils import resources as rsc +from utils import spm_tools as spm +from utils import tools +from src.stir import stir_tools + +#Real Distribution +act_map = + +#Organs diferences for region with differents values +act_map_dif_region = +#Changed of dimensions +ct_image = join(output_dir, "Ct_image_act.img") +ct_image_2 = join(output_dir, "CT_image_act.img") +#address of map +rec_OSEM3D_48_norm.hdr +rec_OSEM3D_48_norm_wholeBody.hdr + +total_quantification = join(output_dir, 'rec_%s_%s_norm.hdr' % (recons_algorithm, recons_it)) +quantification_file = join(output_dir, "Quantification_data.txt") +info_act_map = + + +def rotate_image(act_map, ct_image_act): + # Upload original image (corresponding to the activity map) + act_img = nib.load(act_map) + #act_data = act_img.get_fdata() + act_data = act_img.dataobj[:] # conserva enteros y evita decimales + + + print(f"Image_before_rotate:{act_img.affine}") + + act_header = act_img.header + act_shape = np.array(act_data.shape[:3]) + act_voxel_sizes = act_img.header.get_zooms()[:3] + + # Rotation and flip + rot_xy = np.rot90(act_data, k=2, axes=(0,1)) + final = np.flip(rot_xy, axis=0) + + + # Save final image + final_img = nib.Nifti1Image(final, act_img.affine, act_img.header) + + + output_base = os.path.splitext(ct_image_act)[0] + if not output_base.endswith(".img") and not output_base.endswith(".hdr"): + output_base = output_base # base name without extension + + output_path = output_base + ".img" # nibabel will also generate the .hdr + nib.save(final_img, output_path) + + return final + +def change_act_dimensions(ct_image_act_2, ct_image_act, joint_beds): + """ + Adjust a CT image to match reference voxel size and shape, padding with zeros if necessary. + """ + # --- Load original image --- + img = nib.load(ct_image_act) + data = img.dataobj[:] # conserva enteros y evita decimales + + if data.ndim == 4 and data.shape[3] == 1: + data = np.squeeze(data, axis=3) + + + shape = np.array(data.shape[:3]) + voxel_size = np.array(img.header.get_zooms()[:3]) + size_mm = shape * voxel_size + + # --- Load reference image --- + recons_img = nib.load(joint_beds) + recons_shape = np.array(recons_img.shape[:3]) + recons_voxel_sizes = np.array(recons_img.header.get_zooms()[:3]) + + + # --- Adjust resolution --- + new_voxel_size = recons_voxel_sizes + new_shape = np.round(size_mm / new_voxel_size).astype(int) + scale_factor = new_shape / shape + + + print("Original Dimensions:", data.shape) + print("Scale Factor:", scale_factor) + + # --- Reescalar volumen (sin interpolación continua) --- + new_data = zoom(data, zoom=scale_factor, order=0, prefilter=False) + + # --- Corrige origen geométrico para que el centro quede alineado --- + orig_center = (shape * voxel_size) / 2 + new_center = (new_shape * new_voxel_size) / 2 + shift = orig_center - new_center # desplazamiento en mm + + # --- Nueva cabecera y affine --- + new_affine = recons_img.affine.copy() + new_affine[:3, :3] = np.diag(new_voxel_size) + new_affine[:3, 3] = 0 # centra la imagen correctamente + + # --- Create header and affine --- + nuevo_header = nib.Nifti1Header() + nuevo_header.set_data_shape(new_data.shape) + nuevo_header.set_zooms(new_voxel_size) + + # --- Create intermediate image --- + final_image = nib.Nifti1Image(new_data, affine=new_affine, header=nuevo_header) + + # --- Symmetric padding to match reference shape --- + final_shape = np.array(final_image.shape) + diff = recons_shape - final_shape[:3] + + pad_width = [] + cropped_data = final_image.get_fdata() + + for i in range(3): + if diff[i] >= 0: + # La imagen es más pequeña → agregamos ceros + pad_width.append((diff[i] // 2, diff[i] - diff[i] // 2)) + else: + # La imagen es más grande → recortamos + start = abs(diff[i]) // 2 + end = start + recons_shape[i] + cropped_data = np.take(cropped_data, indices=range(start, end), axis=i) + pad_width.append((0, 0)) + + padded_data = np.pad(final_image.get_fdata(), pad_width, mode='constant', constant_values=0) + + # --- Ajuste final de tamaño exacto --- + # En algunos casos, el padding o recorte previo deja 1 voxel de diferencia por redondeos + final_data = padded_data + for axis in range(3): + if final_data.shape[axis] > recons_shape[axis]: + start = (final_data.shape[axis] - recons_shape[axis]) // 2 + end = start + recons_shape[axis] + final_data = np.take(final_data, indices=range(start, end), axis=axis) + elif final_data.shape[axis] < recons_shape[axis]: + pad_before = (recons_shape[axis] - final_data.shape[axis]) // 2 + pad_after = recons_shape[axis] - final_data.shape[axis] - pad_before + final_data = np.pad( + final_data, + [(pad_before, pad_after) if ax == axis else (0, 0) for ax in range(3)], + mode='constant', + constant_values=0 + ) + # --- Create final NIfTI image --- + final_padded_image = nib.Nifti1Image(final_data, affine=new_affine, header=nuevo_header) + final_padded_image.set_qform(new_affine, code=1) + final_padded_image.set_sform(new_affine, code=1) + + + # --- Save final image --- + nib.save(final_padded_image, ct_image_act_2) + + # --- Delete original image --- + base_name = os.path.splitext(ct_image_act)[0] # quita la extensión .img + for f in [base_name + ".img", base_name + ".hdr"]: + if os.path.exists(f): + os.remove(f) + + return final_padded_image + +def total_quantification(ct_image_act_2, joint_norm_beds, quantification_file): + """ + Compute information of target image per labeled region in reference image and save to a TXT file. + Obtains information about the simulated distribution of activity in the phantom and saves it in a TXT file. + """ + # --- Load images --- + target_img = nib.load(joint_norm_beds) + roi_img = nib.load(ct_image_act_2) + + + target_data = target_img.get_fdata() + roi_data = roi_img.get_fdata() + + + # --- Ensure same shape --- + if target_data.shape != roi_data.shape: + raise ValueError("Target and reference images must have the same shape.") + + # --- Get all labels except 0 (background) --- + labels = np.unique(roi_data) + labels = labels[labels != 0] + + # --- Calculate mean per label --- + mean_values = [target_data[roi_data == label].mean() for label in labels] #KBq/cc + + #--- Calculate Volumen per label --- + voxel_counts = [np.sum(roi_data == label) for label in labels] + dx, dy, dz = np.array(roi_img.header.get_zooms()[:3]) + voxel_volume = dx * dy * dz + + volumes = [(count * voxel_volume)/1000 for count in voxel_counts] + + # --- Calculate activity per region --- + activity_region_KBq = [m * v for m, v in zip(mean_values, volumes)] + activity_region_mCi = [a / (3.7 * 10**4) for a in activity_region_KBq] + mean_values_mCi = [b / c for b, c in zip(activity_region_mCi, volumes)] + + # --- Calculate TOTAL activity across all labels --- + total_activity_KBq = float(np.sum(activity_region_KBq)) + total_activity_mCi = float(np.sum(activity_region_mCi)) + + # --- Save report as TXT --- + with open(quantification_file, 'w') as f: + # Encabezado + f.write( + f"{'Label':<10}{'MeanValue (mCi/cc)':>20}{'MeanValue (KBq/cc)':>20}{'Pixels':>15}" + f"{'Vol (ccm³)':>15}{'Act (mCi)':>15}{'Act (KBq)':>15}\n" + ) + + # Filas por cada label + for label, mean_val_mCi, mean_val_KBq, count, vol, act_mCi, act_KBq in zip( + labels, mean_values_mCi, mean_values, voxel_counts, volumes, activity_region_mCi, activity_region_KBq + ): + f.write( + f"{int(label):<10}{mean_val_mCi:>20.4f}{mean_val_KBq:>20.4f}{count:>15.0f}" + f"{vol:>15.4f}{act_mCi:>15.4f}{act_KBq:>15.4f}\n" + ) + # Linea separadora + f.write("="*90 + "\n") + + # Linea TOTAL + f.write( + f"{'TOTAL':<10}{'':>20}{'':>15}{'':>15}" + f"{total_activity_mCi:>15.4f}{total_activity_KBq:>15.4f}\n" + ) + + return {"labels": labels, "mean_values_mCi": mean_values_mCi, "mean_values_KBq": mean_values, "pixel_counts": voxel_counts, + "volumes": volumes, "activity_mCi": activity_region_mCi, "activity_KBq": activity_region_KBq, + "total_activity_mCi": total_activity_mCi, "total_activity_KBq": total_activity_KBq + } + + + """ + Compute information of target image per labeled region in reference image and save to a TXT file. + Obtains information about the simulated distribution of activity in the phantom and saves it in a TXT file. + """ + # --- Load images --- + target_img = nib.load(recons_norm_wholeBody_file) + roi_img = nib.load(ct_image_act_2) + + + target_data = target_img.get_fdata() + roi_data = roi_img.get_fdata() + + + # --- Ensure same shape --- + if target_data.shape != roi_data.shape: + raise ValueError("Target and reference images must have the same shape.") + + # --- Get all labels except 0 (background) --- + labels = np.unique(roi_data) + labels = labels[labels != 0] + + # --- Calculate mean per label --- + mean_values_KBq = [target_data[roi_data == label].mean() for label in labels] + mean_values_mCi = mean_values_mCi = [v * 0.00002703 for v in mean_values_KBq] #mCi + + #--- Calculate Volumen per label --- + voxel_counts = [np.sum(roi_data == label) for label in labels] + dx, dy, dz = np.array(roi_img.header.get_zooms()[:3]) + voxel_volume = dx * dy * dz + + volumes = [(count * voxel_volume)/1000 for count in voxel_counts] + + # --- Calculate activity per region --- + activity_region_KBq = [m * v for m, v in zip(mean_values_KBq, volumes)] + activity_region_mCi = [a / (3.7 * 10**4) for a in activity_region_KBq] + + # --- Calculate TOTAL activity across all labels --- + total_activity_KBq = float(np.sum(activity_region_KBq)) + total_activity_mCi = float(np.sum(activity_region_mCi)) + + + + # --- Save report as TXT --- + with open(quantification_file, 'w') as f: + # Encabezado + f.write( + f"{'Label':<10}{'MeanValue (mCi/cc)':>20}{'MeanValue (KBq/cc)':>20}{'Pixels':>15}" + f"{'Vol (ccm³)':>15}{'Act (mCi)':>15}{'Act (KBq)':>15}\n" + ) + + # Filas por cada label + for label, mean_val_mCi, mean_val_KBq, count, vol, act_mCi, act_KBq in zip( + labels,mean_values_KBq, mean_values_mCi, voxel_counts, volumes, activity_region_mCi, activity_region_KBq + ): + f.write( + f"{int(label):<10}{mean_val_mCi:>20.4f}{mean_val_KBq:>20.4f}{count:>15.0f}" + f"{vol:>15.4f}{act_mCi:>15.4f}{act_KBq:>15.4f}\n" + ) + # Linea separadora + f.write("="*90 + "\n") + + # Linea TOTAL + f.write( + f"{'TOTAL':<10}{'':>20}{'':>15}{'':>15}" + f"{total_activity_mCi:>15.4f}{total_activity_KBq:>15.4f}\n" + ) + + + return {"labels": labels, "mean_values_mCi": mean_values_mCi, "mean_values_KBq": mean_values_mCi, "pixel_counts": voxel_counts, + "volumes": volumes, "activity_mCi": activity_region_mCi, "activity_KBq": activity_region_KBq, + "total_activity_mCi": total_activity_mCi, "total_activity_KBq": total_activity_KBq + } + +def distribution_of_dose_into_phantom(act_map, info_act_map, self): + """ + Computes the actual distribution of activity by region labeled in act_map. + Obtains information about the actual distribution of activity in the phantom and saves it in a TXT file. + """ + + #Upload the act image + phantom_act = nib.load(act_map) + phantom_act_data = phantom_act.get_fdata() + phantom_voxel_volumen = abs(np.prod(phantom_act.header.get_zooms()[:3]) / 1000) #cm^3 + total_phantom_counts = np.sum(phantom_act_data) + + #Calculate the factor to convert to real phantom activity + if self.sim_dose != 0: + phantom_dose = abs(total_phantom_counts * phantom_voxel_volumen) # uCi + print(f"phantom_dose: {phantom_dose} au") + + act_table_factor = self.sim_dose * 1000 / phantom_dose + else: + act_table_factor = 1 + + print(f"Factor_actividad: {act_table_factor}") + print(f"self.sim_dose: {self.sim_dose} mCi") + print(f"total_phantom_counts: {total_phantom_counts}") + + phantom_real_act = phantom_act_data * (act_table_factor/1000) #se multilica por mil para llevar los valores a mCi + + # Save the file into the same Folder of act_map with different name + base_dir = os.path.dirname(act_map) + name, ext = os.path.splitext(os.path.basename(act_map)) + + if ext.lower() in [".hdr", ".img"]: + output_path = os.path.join(base_dir, f"{name}_real.hdr") + else: + output_path = os.path.join(base_dir, f"{name}_real.nii") + + #New image with the real activity to simulate + new_img = nib.Nifti1Image(phantom_real_act, affine=phantom_act.affine, header=phantom_act.header) + # Save new image + nib.save(new_img, output_path) + + + # Get information about the new image + + # --- Calculate average values per label --- + labels = np.unique(phantom_act_data) + labels = labels[labels != 0] + + mean_values_mCi = [phantom_real_act[phantom_act_data == label].mean() for label in labels] + voxel_counts = [np.sum(phantom_act_data == label) for label in labels] + volumes_region = [count * phantom_voxel_volumen for count in voxel_counts] + + print(f"voxel_counts: {voxel_counts}") + print(f"volumes_region: {volumes_region}") + print(f"phantom_volumen_voxel:{phantom_voxel_volumen}") + + # --- Calculate activity per region --- + activity_region_mCi = [m * v for m, v in zip(mean_values_mCi, volumes_region)] + activity_region_KBq = [a * 3.7 * 10**4 for a in activity_region_mCi] + mean_values_KBq = [b / c for b, c in zip(activity_region_KBq , volumes_region)] + + # --- Calculate TOTAL activity across all labels --- + total_activity_KBq = float(np.sum(activity_region_KBq)) + total_activity_mCi = float(np.sum(activity_region_mCi)) + + print(f"activity_region_KBq: {activity_region_KBq}") + # --- Save report as TXT --- + + with open(info_act_map, 'w') as f: + f.write(f"{'Label':<10}{'MeanValue (mCi/cc)':>20}{'MeanValue (KBq/cc)':>20}{'Pixels':>15}{'Vol (ccm³)':>15}{'Act (mCi)':>15}{'Act (KBq)':>15}\n") + for label, mean_val_mCi, mean_val_KBq, count, vol, act_mCi, act_KBq in zip( + labels, mean_values_mCi, mean_values_KBq, voxel_counts, volumes_region, activity_region_mCi, activity_region_KBq + ): + f.write(f"{int(label):<10}{mean_val_mCi:>20.4f}{mean_val_KBq:>20.4f}{count:>15.0f}{vol:>15.4f}{act_mCi:>15.4f}{act_KBq:>15.4f}\n") + + + # Linea separadora + f.write("="*90 + "\n") + + # Linea TOTAL + f.write( + f"{'TOTAL':<10}{'':>20}{'':>15}{'':>15}" + f"{total_activity_mCi:>15.4f}{total_activity_KBq:>15.4f}\n" + ) + + return { + "image_path": output_path, "factor": act_table_factor, + "labels": labels, "mean_values_mCi": mean_values_mCi,"mean_values_KBq": mean_values_KBq,"volumes": volumes_region, + "activity_mCi": act_mCi, "activity_KBq": act_KBq, "total_activity_mCi": total_activity_mCi, "total_activity_KBq": total_activity_KBq + } diff --git a/utils/wb_tools.py b/utils/wb_tools.py index 9c16b953..2bfaf6ae 100644 --- a/utils/wb_tools.py +++ b/utils/wb_tools.py @@ -1,15 +1,21 @@ import sys +import math #new from os.path import join, dirname import nibabel as nib import numpy as np +import os +from nibabel.processing import resample_from_to +from scipy.ndimage import zoom from scipy.ndimage import gaussian_filter from scipy.ndimage import median_filter +from itertools import product from pyprojroot import here sys.path.append(str(here())) from utils import resources as rsc from utils import spm_tools as spm from utils import tools +from src.stir import stir_tools change_format = rsc.get_rsc('change_format', 'fruitcake') @@ -193,64 +199,60 @@ def pet_to_actmap(self): tools.osrun(rcommand, self.log) -def calculate_center_slices(act_map, scanner, zmin, zmax, overlapping=0.1): +def calculate_center_slices(self, act_map, scanner, zmin, zmax, overlapping=0.1): + """Calculate the center slices of beds for the given axial range.""" + + # Load activity map to get voxel size act_img = nib.load(act_map) - z_voxsize = act_img.affine[2, 2] + self.z_voxsize = abs(act_img.affine[2, 2]) # Axial voxel size (mm) + # Compute axial range z_nvoxels = zmax - zmin + act_z_length = z_nvoxels * self.z_voxsize # mm - act_z_length = z_nvoxels * z_voxsize - - eff_aFOV = scanner.get("axial_fov") * 10 * (1 - 2 * overlapping) - + axial_fov = scanner.get("axial_fov") * 10 # mm + whole_body_simulation = self.params.get("whole_body_simulation") + beds_cs = [] - voxels_per_bed = eff_aFOV / z_voxsize - bed_0_center = zmin + np.rint(voxels_per_bed / 2) - beds_cs.append(bed_0_center) - bed_center = bed_0_center + # Single bed case + if whole_body_simulation == 0 or act_z_length <= axial_fov: + beds_cs = [self.params.get("center_slice")] - while (bed_center < zmax - voxels_per_bed): - bed_center = np.rint(bed_center + voxels_per_bed) - beds_cs.append(bed_center) + print( + "A single bed position is used because:\n" + "- The map length is smaller than the axial FOV, or\n" + "- Whole-body simulation is disabled." + ) - return beds_cs + else: + # Effective FOV considering overlap + eff_aFOV = axial_fov * (1 - 2 * overlapping) + if eff_aFOV <= 0: + raise ValueError("Overlapping too large. Effective axial FOV <= 0.") -def join_beds_wb(recons_beds, joint_beds): - bed_0 = nib.load(recons_beds[0]) - bed_0_data = tools.fix_4d_data(bed_0.get_fdata()) + # Determine number of beds + self.Bed_to_do = max(1, math.ceil(act_z_length / eff_aFOV)) - num_slices = bed_0_data.shape[2] - slices_to_remove = int(np.rint(num_slices / 10)) + # Number of voxels per bed + self.voxels_per_bed = max(1, int(np.ceil(z_nvoxels / self.Bed_to_do))) - bed_0_data = bed_0_data[:, :, :-slices_to_remove] + # First bed center + bed_center = zmin + self.voxels_per_bed / 2 + beds_cs.append(int(np.rint(bed_center))) - for i in range(1, len(recons_beds) - 1): - print(f"Bed: {i}") + # Additional beds + while len(beds_cs) < self.Bed_to_do: + bed_center += self.voxels_per_bed - bed = nib.load(recons_beds[i]) - bed_data = tools.fix_4d_data(bed.get_fdata()) - bed_data = bed_data[:, :, slices_to_remove:-slices_to_remove] - bed_0_data = np.append(bed_0_data, bed_data, axis=2) - - bed_last = nib.load(recons_beds[len(recons_beds) - 1]) - bed_data = tools.fix_4d_data(bed_last.get_fdata()) - bed_data = bed_data[:, :, slices_to_remove:] - bed_0_data = np.append(bed_0_data, bed_data, axis=2) - - bed_0_data = np.flipud(bed_0_data) - bed_0_data = np.fliplr(bed_0_data) + if bed_center + self.voxels_per_bed / 2 > zmax: + break - hdr1 = nib.AnalyzeHeader() - dtype = bed_0.get_data_dtype() - hdr1.set_data_dtype(dtype) - hdr1.set_data_shape(bed_0_data.shape) - affine = bed_0.get_affine() - hdr1.set_zooms((abs(affine[0, 0]), abs(affine[1, 1]), abs(affine[2, 2]))) + beds_cs.append(int(np.rint(bed_center))) - analyze_img = nib.AnalyzeImage(bed_0_data, hdr1.get_base_affine(), hdr1) - nib.save(analyze_img, joint_beds) + self.beds_cs = beds_cs + return beds_cs def update_act_map(spmrun, act_map, att_map, orig_pet, simu_pet, output): @@ -359,3 +361,1000 @@ def cut_image_min_max_slices(input_img, min_slice, max_slice, output): nib.save(analyze_img, output) return output + +def calculate_map_into_fov(self, act_map): + """ + Limits of the FOV per bed (map_into_FOV_start/end). + It is necessary to apply correction_for_NECR and join_beds_wb + """ + act = nib.load(act_map) + beds_cs = calculate_center_slices(self, act_map, self.scanner, self.zmin, self.zmax) + + zmin = int(self.zmin) + zmax = int(self.zmax) + + axial_fov = int(self.scanner["axial_fov"] * 10) # mm + map_voxel_size = float(round(act.affine[2, 2], 3)) # axial voxels size + half_fov_slices = (axial_fov / 2) / map_voxel_size + + self.map_into_FOV_start = [] + self.map_into_FOV_end = [] + + for j, cs in enumerate(beds_cs, start=1): + + # Slices into the FOV + map_into_FOV_start = int(round(cs - half_fov_slices)) + map_into_FOV_end = int(round(cs + half_fov_slices)) + + #Eliminar + #print(f"Before correction: start={map_into_FOV_start}, end={map_into_FOV_end}") + + if map_into_FOV_start < zmin: + map_into_FOV_start = zmin + else: + map_into_FOV_start = int(round(cs - half_fov_slices)) + + if map_into_FOV_end > zmax: + map_into_FOV_end = zmax + else: + map_into_FOV_end = int(round(cs + half_fov_slices)) + + # Save start and end slice for each bed position + self.map_into_FOV_start.append(map_into_FOV_start) + self.map_into_FOV_end.append(map_into_FOV_end) + + print(f"\nNumbers of Slices inside of FOV: [{map_into_FOV_start} , {map_into_FOV_end}]") + +def correction_for_NECR(self, act_map, sim_time_original): + + #Activity correction according to the bed (Bruker PET/MRI) + act = nib.load(act_map) + act_data = act.get_fdata() + + self.params = self.cfg["params"] + self.sim_dose = float(self.params.get("total_dose", 0)) + self.add_randoms = int(self.params.get("add_randoms", 0)) + + #Map into FOV + calculate_map_into_fov(self, act_map) + + self.phantom_doses_FOV = [] + sim_times_per_bed = [] + + zmin = int(self.zmin) + zmax = int(self.zmax) + + for j, (map_start, map_end) in enumerate(zip(self.map_into_FOV_start, self.map_into_FOV_end), start=1): + + #Activity at the time of acquisition according to the bed + phantom_counts_FOV = np.sum(act_data[:, :, int(map_start): int(map_end)]) + phantom_counts_total = np.sum(act_data[:, :, :]) + self.ratio_act_fov = abs (phantom_counts_FOV / phantom_counts_total) + + self.phantom_dose_FOV = float(round(abs(self.sim_dose * self.ratio_act_fov), 3)) + self.phantom_doses_FOV.append(self.phantom_dose_FOV) + + #TODO Use the fitted equation. + if self.add_randoms ==1: + y = abs(4.55 * self.phantom_dose_FOV + 1.9775) + else: + y = abs (4.569 * self.phantom_dose_FOV + 1.959) + + + print(f"\nBeds_to_perform: {j}") + print(f"Simulation Dose inside FOV: {float(round(self.phantom_dose_FOV, 3))} mCi") + + sim_time_bed = float(sim_time_original / y) + sim_time_bed = float(round(sim_time_bed, 1)) + sim_times_per_bed.append(sim_time_bed) + + print(f"Corrected simulation time = {sim_time_bed} seg") + + + # ####BORRAR + # #labels into each beds + # # Extract only the portion of the map within the FOV (Z-axis) + # # Assuming that axis 2 (index 2) is the Z-axis + act_data_FOV = act_data[:, :, map_start:map_end] + # # Get the labels within that range (excluding the background label, 0) + labels = np.unique(act_data_FOV) + labels = labels[labels != 0] + + # Contar voxeles por label dentro del FOV + voxel_counts = [np.sum(act_data_FOV == label) for label in labels] + phantom_voxel_volumen = abs(np.prod(act.header.get_zooms()[:3]) / 1000) #cm^3 + volumes_region = [count * phantom_voxel_volumen for count in voxel_counts] + + phantom_act = nib.load(act_map) + phantom_act_data = phantom_act.get_fdata() + phantom_voxel_volumen = abs(np.prod(phantom_act.header.get_zooms()[:3]) / 1000) #cm^3 + total_phantom_counts = np.sum(phantom_act_data) + + #Calculate the factor to convert to real phantom activity + if self.sim_dose != 0: + phantom_dose = abs(total_phantom_counts * phantom_voxel_volumen) # uCi + print(f"phantom_dose: {phantom_dose} au") + self.act_table_factor = self.sim_dose * 1000 / phantom_dose + else: + self.act_table_factor = 1 + + + new_concent_label = [a * self.act_table_factor for a in labels] + activity_region_uCi = [m * v for m, v in zip(new_concent_label, volumes_region)] + actividad_total_mCi = sum(activity_region_uCi)/1000 + + print(f"Actividad Total: {actividad_total_mCi} mCi") + + # Mostrar resumen por cama + print(f"\nBed {j}:") + for lbl, vox , volum, concent, act_mCi in zip(labels, voxel_counts, volumes_region, new_concent_label, activity_region_uCi): + print(f" Label {lbl}: {vox} voxeles dentro del FOV, Volumen de region {volum}, New_concentration: {concent}uci/cc, Activity uCi: {act_mCi}") + + return sim_times_per_bed, self.phantom_doses_FOV + +def join_beds_wb(self, act_map, recons_beds, joint_beds): + """ + Joins the reconstructed beds correcting overlap dynamically. + The trimming is split half to the previous bed and half to the next. + Saves the final image in Analyze format (.hdr/.img) as in the original function. + """ + #Map into FOV + calculate_map_into_fov(self, act_map) + + # Load activity map to get voxel size + act = nib.load(act_map) + act_voxel_size = float(round(act.affine[2, 2], 3)) # mm + + # Load first bed + bed_prev = nib.load(recons_beds[0]) + bed_prev_data = tools.fix_4d_data(bed_prev.get_fdata()) + img_recons_voxel_size = float(round(bed_prev.affine[2, 2], 3)) #mm + + for i in range(len(recons_beds) - 1): + # Load next bed + bed_next = nib.load(recons_beds[i + 1]) + bed_next_data = tools.fix_4d_data(bed_next.get_fdata()) + + # Calculate total dynamic trimming (only if there is overlap) + end_i = self.map_into_FOV_end[i] + start_next = self.map_into_FOV_start[i + 1] + diff = end_i - start_next + + if diff > 0: + slices_to_remove = int(round((diff * act_voxel_size) / img_recons_voxel_size)) + else: + slices_to_remove = 0 # No overlap + + # Split the trimming evenly + remove_prev = slices_to_remove // 2 + remove_next = slices_to_remove - remove_prev # covers the rest if odd + + print(f"Bed {i} end={end_i}, Bed {i+1} start={start_next} → diff={diff}, total remove={slices_to_remove}, remove_prev={remove_prev}, remove_next={remove_next}") + + # Trim previous bed (last slices) + if remove_prev > 0 and remove_prev < bed_prev_data.shape[2]: + bed_prev_data = bed_prev_data[:, :, :-remove_prev] + + # Trim next bed (first slices) + if remove_next > 0 and remove_next < bed_next_data.shape[2]: + bed_next_data = bed_next_data[:, :, remove_next:] + + # Concatenate beds + bed_prev_data = np.append(bed_prev_data, bed_next_data, axis=2) + + # --- Apply final flips as in the original function --- + bed_prev_data = np.flipud(bed_prev_data) + #bed_prev_data = np.fliplr(bed_prev_data) + + # --- Save the concatenated image as in the original function --- + hdr1 = nib.AnalyzeHeader() + dtype = bed_prev.get_data_dtype() + hdr1.set_data_dtype(dtype) + hdr1.set_data_shape(bed_prev_data.shape) + affine = bed_prev.get_affine() + hdr1.set_zooms((abs(affine[0, 0]), abs(affine[1, 1]), abs(affine[2, 2]))) + + analyze_img = nib.AnalyzeImage(bed_prev_data, hdr1.get_base_affine(), hdr1) + nib.save(analyze_img, joint_beds) + + print(f"\nConcatenated image saved at: {joint_beds}") + print(f"Data Type: {dtype}") + print(f"Matrix Size: {bed_prev_data.shape}") + print(f"Voxel Size: {hdr1.get_zooms()}") + +def normalization_factor_correction(self): + #factor_Q_norm = 400.13 #480.43 + results = [] + + print(f"Phantom dose (FOV):{self.phantom_doses_FOV}") + + for dose_mCi in self.phantom_doses_FOV: + # Convert mCi to kBq + phantom_dose_KBq = dose_mCi * 3.7e4 + sim_time_original_global = float(self.params.get("simulation_time", 0)) + print(f"Phantom dose in KBq: {phantom_dose_KBq}") #borrar + + # Linear adjustment equation + #lineal_ecuac = 5.9e-5 * phantom_dose_KBq + 0.978 + #value_Q_norm = float(lineal_ecuac * factor_Q_norm) + #value_Q_norm = float(0.061 * phantom_dose_KBq + 904.2) + concentrac = (phantom_dose_KBq / 98.96) + print(f"Phantom conc in KBq/cc: {concentrac}") + + + #voi 30 x 140 Imagen Simulada vs Concentracion teorica experimental + value_Q_norm = float(0.062 * (phantom_dose_KBq) + 986.39) + + value_Q_norm_corregido_time = value_Q_norm / (sim_time_original_global/300) # 300 seg es el tiempo de adq cyl calibración + print(f"Normalization value per bed: {value_Q_norm}") #borrar + #results.append(value_Q_norm) + results.append(value_Q_norm_corregido_time) + + + if len(results) == 0: + print(f"There is no FOV dose. Using factor = 1.0 by default") + return 1.0 + + print("Normalization factors by FOV:", [round(x, 3) for x in results]) + return results + +def normalization_factor_correction_whole_body(self, joint_beds): + factor_Q_norm = 940.34 + results = [] + sim_dose_original_global = float(self.params.get("total_dose", 0)) + print(f"Phantom dose (FOV):{sim_dose_original_global}") + + image_sim = nib.load(joint_beds) + image_sim_data = image_sim.get_fdata() + + # Reemplazar NaN o inf por 0 (solo si hay) + image_sim_data = np.nan_to_num(image_sim_data, nan=0.0, posinf=0.0, neginf=0.0) + + # Calcular volumen del voxel + dx, dy, dz = np.array(image_sim.header.get_zooms()[:3]) + voxel_volume = (dx * dy * dz)/1000 #cm^3 + + + # Calcular sumas + total_voxel_sum = np.sum(image_sim_data) + Activity_image = (voxel_volume * total_voxel_sum)/10 + + # Calcular valor medio + # Numero total de voxeles + #num_voxels = image_sim_data.size + + # Valor medio de todos los voxeles + #mean_voxel_value = total_voxel_sum / num_voxels + + # Creamos una máscara para voxeles distintos de cero + mask = image_sim_data != 0 + # Sumamos solo los voxeles distintos de cero + sum_nonzero = np.sum(image_sim_data[mask]) + # Contamos cuántos voxeles son distintos de cero + num_nonzero_voxels = np.count_nonzero(image_sim_data) + volumen_non_zero = voxel_volume * num_nonzero_voxels + print(f"Volumen_non_zero:{volumen_non_zero} cc") + # Valor medio ignorando ceros + #mean_nonzero = sum_nonzero / num_nonzero_voxels + + #Calcular la concentracion estimada segun el total sum de la imagen + + #conc_estimate_image = 1.34*10**-5 * sum_nonzero - 8.96*10**-9 + act_estimate_image = float(4.22*10**-4 * sum_nonzero - 4.65*10**-8) #volumen de la voi en cmm³ + conc_estimate_image = float(act_estimate_image/98.96) #volumen de la voi en cmm³ + #print(mean_nonzero) + print(f"Total_suma:{sum_nonzero} au/cc") + print(f"Concentracion estimada:{conc_estimate_image} au/cc Actividad estimada: {act_estimate_image}") + + print(f"Voxel volume: {voxel_volume:.3f} mm³") + print(f"Total_sum : {total_voxel_sum}") + print(f"Activity_image : {Activity_image}") + + + # Convert mCi to kBq + phantom_dose_KBq = sim_dose_original_global * 3.7*10**4 + #mean_nonzero = phantom_dose_KBq / num_nonzero_voxels + + #print(f"Phantom mean dose in KBq: {mean_nonzero}") + print(f"Phantom dose in KBq: {phantom_dose_KBq}") + + # Linear adjustment equation + #lineal_ecuac = (0.0038 * (phantom_dose_KBq / 31.43) + 0.9725) + lineal_ecuac = (0.0000606 * phantom_dose_KBq + 0.975) + print(f"Lineal Ecuation: {lineal_ecuac}") + + #value_Q_norm = float(lineal_ecuac * conc_estimate_image) + #value_Q_norm = float(lineal_ecuac * conc_estimate_image * factor_Q_norm) + value_Q_norm = float((lineal_ecuac * factor_Q_norm)/1.47) + #value_Q_norm = float(factor_Q_norm * lineal_ecuac ) + #value_Q_norm = float(lineal_ecuac) + + print(f"Normalization value per total: {value_Q_norm}") + results.append(value_Q_norm) + + + + if len(results) == 0: + print(f"There is no FOV dose. Using factor = 1.0 by default") + return 1.0 + + print("Normalization factors by FOV:", [round(x, 3) for x in results]) + return results + +def rotate_and_flip_mask(act_map, mask_map, rotated_mask_file, joint_beds, zmin, zmax): + """ + Rotates and flips a NIfTI mask (mask_map) based on the best alignment + calculated from act_map compared to a reference (joint_beds). + + PARAMETERS: + - act_map: path to activity map NIfTI + - mask_map: path to mask NIfTI to be rotated/flipped + - rotated_mask_file: output path for transformed mask + - joint_beds: path to reference NIfTI + - zmin, zmax: axial range of interest (inclusive, MRIcro-style) + + IMPORTANT: + - Original voxel intensities in mask_map are preserved. + - Binarization is performed ONLY on a copy of act_map for Dice calculation. + - Final transformation is applied to mask_map. + """ + + def overlap_coefficient(a, b): + """Computes overlap similarity between two binary masks.""" + a = a > 0 + b = b > 0 + intersection = np.sum(a & b) + return 2 * intersection / (np.sum(a) + np.sum(b) + 1e-8) + + def crop_z(data, zmin, zmax): + """ + Crops a volume along the axial axis including zmin and zmax (MRIcro-style, inclusive) + """ + zmin = int(zmin) + zmax = int(zmax) + + zmax = zmax -1 + if zmax is None: + zmax = data.shape[2] + zmin_adj = max(0, zmin - 1) # Ajuste inicio 0-indexed Python + zmax_adj = min(zmax, data.shape[2] - 1) # Ajuste fin + return data[:, :, zmin_adj:zmax_adj + 1] # +1 para incluir zmax + + try: + # --- Load mask to transform --- + mask_img = nib.load(mask_map) + mask_data_original = mask_img.get_fdata() + if mask_data_original.ndim == 4 and mask_data_original.shape[3] == 1: + mask_data_original = np.squeeze(mask_data_original, axis=3) + mask_data_original = crop_z(mask_data_original, zmin, zmax) + + # --- Load act_map for calculating best transform --- + act_img = nib.load(act_map) + act_data = act_img.get_fdata() + if act_data.ndim == 4 and act_data.shape[3] == 1: + act_data = np.squeeze(act_data, axis=3) + act_data = crop_z(act_data, zmin, zmax) + act_data_binary = (act_data > 0).astype(np.uint8) + + # --- Load reference --- + ref_img = nib.load(joint_beds) + ref_data = ref_img.get_fdata() + if ref_data.ndim == 4 and ref_data.shape[3] == 1: + ref_data = np.squeeze(ref_data, axis=3) + ref_data = crop_z(ref_data, zmin, zmax) + ref_data_binary = (ref_data > 0).astype(np.uint8) + + # --- Rescale act_map binary for overlap computation --- + act_voxel_size = act_img.header.get_zooms()[:3] + ref_voxel_size = ref_img.header.get_zooms()[:3] + scale_factors = np.array(act_voxel_size) / np.array(ref_voxel_size) + act_rescaled = zoom(act_data_binary, zoom=scale_factors, order=0, prefilter=False) + + # --- Pad volumes for overlap computation --- + max_shape = np.maximum(act_rescaled.shape, ref_data_binary.shape) + + def pad_to_shape(data, target_shape): + pad_width = [] + for s, t in zip(data.shape, target_shape): + total = t - s + before = total // 2 + after = total - before + pad_width.append((before, after)) + return np.pad(data, pad_width, mode='constant', constant_values=0) + + act_padded = pad_to_shape(act_rescaled, max_shape) + ref_padded = pad_to_shape(ref_data_binary, max_shape) + + # --- Search best rotation/flip using act_map --- + best_score = -1 + best_transform = (0, False, False) + + for k in range(4): + rotated = np.rot90(act_padded, k=k, axes=(0, 1)) + for flip_x, flip_y in product([False, True], repeat=2): + candidate = rotated.copy() + if flip_x: + candidate = np.flip(candidate, axis=0) + if flip_y: + candidate = np.flip(candidate, axis=1) + + score = overlap_coefficient(candidate, ref_padded) + if score > best_score: + best_score = score + best_transform = (k, flip_x, flip_y) + + # --- Apply best transform to ORIGINAL mask_map data --- + transformed_mask = np.rot90(mask_data_original, k=best_transform[0], axes=(0, 1)) + if best_transform[1]: + transformed_mask = np.flip(transformed_mask, axis=0) + if best_transform[2]: + transformed_mask = np.flip(transformed_mask, axis=1) + + # --- Save result --- + new_affine = mask_img.affine.copy() + final_img = nib.Nifti1Image(transformed_mask, affine=new_affine, header=mask_img.header.copy()) + final_img.set_qform(new_affine, code=1) + final_img.set_sform(new_affine, code=1) + + if not rotated_mask_file.endswith(".nii"): + rotated_mask_file = os.path.splitext(rotated_mask_file)[0] + ".nii" + + nib.save(final_img, rotated_mask_file) + return rotated_mask_file + + except Exception as e: + print("\n=== ERROR during rotate_and_flip_mask ===") + print("Error type:", type(e)) + print("Message:", e) + raise + +def change_act_dimensions(mask, rotated_mask, joint_beds): + """ + Adjust a mask image to match reference voxel size and shape, padding with zeros if necessary. + """ + # --- Load original image --- + img = nib.load(rotated_mask) + data = img.dataobj[:] # preserve integers and avoid decimals + + if data.ndim == 4 and data.shape[3] == 1: + data = np.squeeze(data, axis=3) + + shape = np.array(data.shape[:3]) + voxel_size = np.array(img.header.get_zooms()[:3]) + size_mm = shape * voxel_size + + # --- Load reference image --- + recons_img = nib.load(joint_beds) + recons_shape = np.array(recons_img.shape[:3]) + recons_voxel_sizes = np.array(recons_img.header.get_zooms()[:3]) + + # --- Adjust resolution --- + new_voxel_size = recons_voxel_sizes + new_shape = np.round(size_mm / new_voxel_size).astype(int) + scale_factor = new_shape / shape + #print("Original Dimensions:", data.shape) + #print("Scale Factor:", scale_factor) + + # --- Rescale volume (without continuous interpolation) --- + new_data = zoom(data, zoom=scale_factor, order=0, prefilter=False) + + # --- Correct geometric origin to align center --- + orig_center = (shape * voxel_size) / 2 + new_center = (new_shape * new_voxel_size) / 2 + shift = orig_center - new_center # shift in mm + + # --- New header and affine --- + new_affine = recons_img.affine.copy() + new_affine = np.eye(4) # complete 4x4 matrix + new_affine[:3, :3] = np.diag(new_voxel_size) + new_affine[:3, 3] = np.zeros(3) + + # --- Create header and affine --- + nuevo_header = nib.Nifti1Header() + nuevo_header.set_data_shape(new_data.shape) + nuevo_header.set_zooms(new_voxel_size) + + # --- Create intermediate image --- + final_image = nib.Nifti1Image(new_data, affine=new_affine, header=nuevo_header) + + # --- Symmetric padding to match reference shape --- + final_shape = np.array(final_image.shape) + diff = recons_shape - final_shape[:3] + pad_width = [] + cropped_data = final_image.get_fdata() + + for i in range(3): + if diff[i] >= 0: + # Image is smaller - add zeros + pad_width.append((diff[i] // 2, diff[i] - diff[i] // 2)) + else: + # Image is larger - crop + start = abs(diff[i]) // 2 + end = start + recons_shape[i] + cropped_data = np.take(cropped_data, indices=range(start, end), axis=i) + pad_width.append((0, 0)) + + padded_data = np.pad(cropped_data, pad_width, mode='constant', constant_values=0) + + # --- Final adjustment to exact size --- + final_data = padded_data + for axis in range(3): + if final_data.shape[axis] > recons_shape[axis]: + start = (final_data.shape[axis] - recons_shape[axis]) // 2 + end = start + recons_shape[axis] + final_data = np.take(final_data, indices=range(start, end), axis=axis) + elif final_data.shape[axis] < recons_shape[axis]: + pad_before = (recons_shape[axis] - final_data.shape[axis]) // 2 + pad_after = recons_shape[axis] - final_data.shape[axis] - pad_before + final_data = np.pad( + final_data, + [(pad_before, pad_after) if ax == axis else (0, 0) for ax in range(3)], + mode='constant', + constant_values=0 + ) + + # --- Create final NIfTI image --- + final_padded_image = nib.Nifti1Image(final_data, affine=new_affine, header=nuevo_header) + final_padded_image.set_qform(new_affine, code=1) + final_padded_image.set_sform(new_affine, code=1) + + # --- Save final image --- + nib.save(final_padded_image, mask) + + # --- Delete the .nii file --- + base_name = os.path.splitext(rotated_mask)[0] # remove the extension + for f in [base_name + ".nii"]: + if os.path.exists(f): + os.remove(f) + + return final_padded_image + +def total_quantification(mask_file, joint_norm_beds, quantification_file, label_file_mask): + """ + Compute information of target image per labeled region in reference image and save to a TXT file. + """ + + import nibabel as nib + import numpy as np + import os + + # --- Load images --- + target_img = nib.load(joint_norm_beds) + roi_img = nib.load(mask_file) + + target_data = target_img.get_fdata() + roi_data = roi_img.get_fdata().astype(int) + + # --- Ensure same shape --- + if target_data.shape != roi_data.shape: + raise ValueError("Target and reference images must have the same shape.") + + # --- Load label names --- + label_names = {} + + if os.path.exists(label_file_mask): + with open(label_file_mask, 'r') as f: + for line in f: + line = line.strip() + if not line or line.startswith("#"): + continue + parts = line.split() # debe estar dentro del for + label_id = int(parts[0]) + label_name = " ".join(parts[1:]) + label_names[label_id] = label_name + else: + print(f"INFO: Label file not found: {label_file_mask}. Using numeric labels from ROI.") + + + + # --- Get all labels in the mask except 0 (background) --- + labels = np.unique(roi_data) + labels = labels[labels != 0] + + # --- Calculate mean per label --- + mean_values = [target_data[roi_data == label].mean() for label in labels] # KBq/cc + + # --- Calculate volume per label --- + voxel_counts = [np.sum(roi_data == label) for label in labels] + dx, dy, dz = np.array(roi_img.header.get_zooms()[:3]) + voxel_volume = dx * dy * dz # mm³ + volumes = [(count * voxel_volume) / 1000 for count in voxel_counts] # cm³ + + # --- Calculate activity per region --- + activity_region_KBq = [m * v for m, v in zip(mean_values, volumes)] + activity_region_mCi = [a / (3.7 * 10**4) for a in activity_region_KBq] + mean_values_mCi = [b / c for b, c in zip(activity_region_mCi, volumes)] + + # --- Calculate TOTAL activity --- + total_activity_KBq = float(np.sum(activity_region_KBq)) + total_activity_mCi = float(np.sum(activity_region_mCi)) + + # --- Save report as TXT --- + with open(quantification_file, 'w') as f: + # Header + f.write( + f"{'Region':<15}{'MeanValue (mCi/cc)':>20}{'MeanValue (KBq/cc)':>20}" + f"{'Pixels':>15}{'Vol (cm³)':>15}{'Act (mCi)':>15}{'Act (KBq)':>15}\n" + ) + + # Rows + for label, mean_val_mCi, mean_val_KBq, count, vol, act_mCi, act_KBq in zip( + labels, mean_values_mCi, mean_values, + voxel_counts, volumes, activity_region_mCi, activity_region_KBq + ): + label_name = label_names.get(int(label), f"Label_{int(label)}") + f.write( + f"{label_name:<15}{mean_val_mCi:>20.4f}{mean_val_KBq:>20.4f}" + f"{count:>15.0f}{vol:>15.4f}{act_mCi:>15.4f}{act_KBq:>15.4f}\n" + ) + + # Separator + f.write("=" * 115 + "\n") + + # TOTAL + f.write( + f"{'TOTAL':<15}{'':>20}{'':>20}{'':>15}{'':>15}" + f"{total_activity_mCi:>15.4f}{total_activity_KBq:>15.4f}\n" + ) + + return { + "labels": labels, + "label_names": label_names, + "mean_values_mCi": mean_values_mCi, + "mean_values_KBq": mean_values, + "pixel_counts": voxel_counts, + "volumes": volumes, + "activity_mCi": activity_region_mCi, + "activity_KBq": activity_region_KBq, + "total_activity_mCi": total_activity_mCi, + "total_activity_KBq": total_activity_KBq + } + + + """ + Compute information of target image per labeled region in reference image and save to a TXT file. + Obtains information about the simulated distribution of activity in the phantom and saves it in a TXT file. + """ + # --- Load images --- + target_img = nib.load(joint_norm_beds) + roi_img = nib.load(ct_image_act_2) + + + target_data = target_img.get_fdata() + roi_data = roi_img.get_fdata() + + + # --- Ensure same shape --- + if target_data.shape != roi_data.shape: + raise ValueError("Target and reference images must have the same shape.") + + # --- Get all labels except 0 (background) --- + labels = np.unique(roi_data) + labels = labels[labels != 0] + + # --- Calculate mean per label --- + mean_values = [target_data[roi_data == label].mean() for label in labels] #KBq/cc + + #--- Calculate Volumen per label --- + voxel_counts = [np.sum(roi_data == label) for label in labels] + dx, dy, dz = np.array(roi_img.header.get_zooms()[:3]) + voxel_volume = dx * dy * dz + + volumes = [(count * voxel_volume)/1000 for count in voxel_counts] + + # --- Calculate activity per region --- + activity_region_KBq = [m * v for m, v in zip(mean_values, volumes)] + activity_region_mCi = [a / (3.7 * 10**4) for a in activity_region_KBq] + mean_values_mCi = [b / c for b, c in zip(activity_region_mCi, volumes)] + + # --- Calculate TOTAL activity across all labels --- + total_activity_KBq = float(np.sum(activity_region_KBq)) + total_activity_mCi = float(np.sum(activity_region_mCi)) + + # --- Save report as TXT --- + with open(quantification_file, 'w') as f: + # Encabezado + f.write( + f"{'Label':<10}{'MeanValue (mCi/cc)':>20}{'MeanValue (KBq/cc)':>20}{'Pixels':>15}" + f"{'Vol (ccm³)':>15}{'Act (mCi)':>15}{'Act (KBq)':>15}\n" + ) + + # Filas por cada label + for label, mean_val_mCi, mean_val_KBq, count, vol, act_mCi, act_KBq in zip( + labels, mean_values_mCi, mean_values, voxel_counts, volumes, activity_region_mCi, activity_region_KBq + ): + f.write( + f"{int(label):<10}{mean_val_mCi:>20.4f}{mean_val_KBq:>20.4f}{count:>15.0f}" + f"{vol:>15.4f}{act_mCi:>15.4f}{act_KBq:>15.4f}\n" + ) + # Linea separadora + f.write("="*90 + "\n") + + # Linea TOTAL + f.write( + f"{'TOTAL':<10}{'':>20}{'':>15}{'':>15}" + f"{total_activity_mCi:>15.4f}{total_activity_KBq:>15.4f}\n" + ) + + return {"labels": labels, "mean_values_mCi": mean_values_mCi, "mean_values_KBq": mean_values, "pixel_counts": voxel_counts, + "volumes": volumes, "activity_mCi": activity_region_mCi, "activity_KBq": activity_region_KBq, + "total_activity_mCi": total_activity_mCi, "total_activity_KBq": total_activity_KBq + } + +def total_quantification_wholeBody(mask_file, recons_norm_wholeBody_file, quantification_file): #TODO + """ + Compute information of target image per labeled region in reference image and save to a TXT file. + Obtains information about the simulated distribution of activity in the phantom and saves it in a TXT file. + """ + # --- Load images --- + target_img = nib.load(recons_norm_wholeBody_file) + roi_img = nib.load(mask_file) + + + target_data = target_img.get_fdata() + roi_data = roi_img.get_fdata() + + + # --- Ensure same shape --- + if target_data.shape != roi_data.shape: + raise ValueError("Target and reference images must have the same shape.") + + # --- Get all labels except 0 (background) --- + labels = np.unique(roi_data) + labels = labels[labels != 0] + + # --- Calculate mean per label --- + mean_values_KBq = [target_data[roi_data == label].mean() for label in labels] + mean_values_mCi = mean_values_mCi = [v * 0.00002703 for v in mean_values_KBq] #mCi + + #--- Calculate Volumen per label --- + voxel_counts = [np.sum(roi_data == label) for label in labels] + dx, dy, dz = np.array(roi_img.header.get_zooms()[:3]) + voxel_volume = dx * dy * dz + + volumes = [(count * voxel_volume)/1000 for count in voxel_counts] + + # --- Calculate activity per region --- + activity_region_KBq = [m * v for m, v in zip(mean_values_KBq, volumes)] + activity_region_mCi = [a / (3.7 * 10**4) for a in activity_region_KBq] + + # --- Calculate TOTAL activity across all labels --- + total_activity_KBq = float(np.sum(activity_region_KBq)) + total_activity_mCi = float(np.sum(activity_region_mCi)) + + + + # --- Save report as TXT --- + with open(quantification_file, 'w') as f: + # Encabezado + f.write( + f"{'Label':<10}{'MeanValue (mCi/cc)':>20}{'MeanValue (KBq/cc)':>20}{'Pixels':>15}" + f"{'Vol (ccm³)':>15}{'Act (mCi)':>15}{'Act (KBq)':>15}\n" + ) + + # Filas por cada label + for label, mean_val_mCi, mean_val_KBq, count, vol, act_mCi, act_KBq in zip( + labels, mean_values_mCi, mean_values_KBq, voxel_counts, volumes, activity_region_mCi, activity_region_KBq + ): + f.write( + f"{int(label):<10}{mean_val_mCi:>20.4f}{mean_val_KBq:>20.4f}{count:>15.0f}" + f"{vol:>15.4f}{act_mCi:>15.4f}{act_KBq:>15.4f}\n" + ) + # Linea separadora + f.write("="*90 + "\n") + + # Linea TOTAL + f.write( + f"{'TOTAL':<10}{'':>20}{'':>15}{'':>15}" + f"{total_activity_mCi:>15.4f}{total_activity_KBq:>15.4f}\n" + ) + + + return {"labels": labels, "mean_values_mCi": mean_values_mCi, "mean_values_KBq": mean_values_mCi, "pixel_counts": voxel_counts, + "volumes": volumes, "activity_mCi": activity_region_mCi, "activity_KBq": activity_region_KBq, + "total_activity_mCi": total_activity_mCi, "total_activity_KBq": total_activity_KBq + } + +def distribution_of_dose_into_phantom(self, maps_dir, act_map, info_act_map, label_file_act_map): + """ + Computes the actual distribution of activity by region labeled in act_map. + Obtains information about the actual distribution of activity in the phantom and saves it in a TXT file. + """ + + #Upload the act image + phantom_act = nib.load(act_map) + phantom_act_data = phantom_act.get_fdata() + phantom_voxel_volumen = abs(np.prod(phantom_act.header.get_zooms()[:3]) / 1000) #cm^3 + total_phantom_counts = np.sum(phantom_act_data) + + #Calculate the factor to convert to real phantom activity + if self.sim_dose != 0: + phantom_dose = abs(total_phantom_counts * phantom_voxel_volumen) # uCi + #print(f"phantom_dose: {phantom_dose} au") + + act_table_factor = self.sim_dose * 1000 / phantom_dose + else: + act_table_factor = 1 + + #print(f"Factor_actividad: {act_table_factor}") + #print(f"self.sim_dose: {self.sim_dose} mCi") + #print(f"total_phantom_counts: {total_phantom_counts}") + + phantom_real_act = phantom_act_data * (act_table_factor/1000) #se multilica por mil para llevar los valores a mCi + + # Save the file into the same Folder of act_map with different name + base_dir = os.path.dirname(act_map) + name, ext = os.path.splitext(os.path.basename(act_map)) + + if ext.lower() in [".hdr", ".img"]: + output_path = os.path.join(base_dir, f"{name}_real.hdr") + else: + output_path = os.path.join(base_dir, f"{name}_real.nii") + + #New image with the real activity to simulate + new_img = nib.Nifti1Image(phantom_real_act, affine=phantom_act.affine, header=phantom_act.header) + # Save new image + nib.save(new_img, output_path) + + # --- Load label names --- + label_names = {} + + if label_file_act_map and os.path.exists(label_file_act_map): + with open(label_file_act_map, 'r') as f: + for line in f: + line = line.strip() + if not line or line.startswith("#"): + continue + parts = line.split() + if len(parts) < 2: + # Linea invalida, ignorar + continue + try: + label_id = int(parts[0]) + label_name = " ".join(parts[1:]) + label_names[label_id] = label_name + except ValueError: + continue # Si no se puede convertir a entero, ignorar + else: + print(f"No label file found: Using ROI numbers as labels") + + + + # --- Get labels --- + labels = np.unique(phantom_act_data) + labels = labels[labels != 0] + + # --- Calculate stats per label --- + mean_values_mCi = [phantom_real_act[phantom_act_data == label].mean() for label in labels] + voxel_counts = [np.sum(phantom_act_data == label) for label in labels] + volumes_region = [count * phantom_voxel_volumen for count in voxel_counts] + activity_region_mCi = [m * v for m, v in zip(mean_values_mCi, volumes_region)] + activity_region_KBq = [a * 3.7 * 10**4 for a in activity_region_mCi] + mean_values_KBq = [b / c for b, c in zip(activity_region_KBq, volumes_region)] + + total_activity_KBq = float(np.sum(activity_region_KBq)) + total_activity_mCi = float(np.sum(activity_region_mCi)) + + # --- Save report as TXT --- + with open(info_act_map, 'w') as f: + f.write(f"{'Region':<15}{'MeanValue (mCi/cc)':>20}{'MeanValue (KBq/cc)':>20}" + f"{'Pixels':>15}{'Vol (cm³)':>15}{'Act (mCi)':>15}{'Act (KBq)':>15}\n") + + for label, mean_val_mCi, mean_val_KBq, count, vol, act_mCi, act_KBq in zip( + labels, mean_values_mCi, mean_values_KBq, voxel_counts, volumes_region, activity_region_mCi, activity_region_KBq): + + # Sustituir número por nombre + label_name = label_names.get(int(label), str(int(label))) + f.write(f"{label_name:<15}{mean_val_mCi:>20.4f}{mean_val_KBq:>20.4f}" + f"{count:>15.0f}{vol:>15.4f}{act_mCi:>15.4f}{act_KBq:>15.4f}\n") + + # Separador y total + f.write("="*115 + "\n") + f.write(f"{'TOTAL':<15}{'':>20}{'':>20}{'':>15}{'':>15}" + f"{total_activity_mCi:>15.4f}{total_activity_KBq:>15.4f}\n") + + return { + "image_path": output_path, + "factor": act_table_factor, + "labels": labels, + "label_names": label_names, + "mean_values_mCi": mean_values_mCi, + "mean_values_KBq": mean_values_KBq, + "volumes": volumes_region, + "activity_mCi": activity_region_mCi, + "activity_KBq": activity_region_KBq, + "total_activity_mCi": total_activity_mCi, + "total_activity_KBq": total_activity_KBq + } + +def coincidencias_mask_vs_image_Claudia(mask_file, act_map): #ELIMINAR ES SOLO PARA COMPROBACIÓN + # --- Verificación de volúmenes por label --- + import numpy as np + import nibabel as nib + + try: + orig_img = nib.load(act_map) # imagen original base + new_img = nib.load(mask_file) + + orig_data = orig_img.get_fdata() + new_data = new_img.get_fdata() + + voxel_orig = np.prod(orig_img.header.get_zooms()[:3]) / 1000 # mm³ -> cm³ + voxel_new = np.prod(new_img.header.get_zooms()[:3]) / 1000 + + labels = np.unique(orig_data) + labels = labels[labels != 0] + + + print("\n--- Comparación de volúmenes por región ---") + for label in labels: + count_orig = np.sum(orig_data == label) + count_new = np.sum(new_data == label) + vol_orig = count_orig * voxel_orig + vol_new = count_new * voxel_new + ratio = vol_new / vol_orig if vol_orig > 0 else np.nan + + print(f"Label {int(label):3d}: original={vol_orig:.3f} cm³, new={vol_new:.3f} cm³, ratio={ratio:.3f}") + except Exception as e: + print(f"[Aviso] No se pudo comparar volúmenes por label: {e}") + +def total_fov_correction(self, recons_dir): + + """ + This function performs the STIR reconstruction sinogram correction, if desired, + using a sinogram obtained from a cylinder that fills the entire FOV. + IT IS RECOMMENDED TO NORMALIZE THE SINOGRAM USING ITS MEAN VALUE. + """ + import shutil + + whole_FOV_dir_path = self.config.get("dir_stir_sino_corrections") + whole_FOV_sino_file = os.path.join(whole_FOV_dir_path, "sinogram_TOTAL_FOV_divid_VM.img") + whole_FOV_stir_sino = nib.load(whole_FOV_sino_file) + #whole_FOV_stir_sino_data = whole_FOV_stir_sino.get_fdata() + whole_FOV_stir_sino_data = np.squeeze(whole_FOV_stir_sino.get_fdata()) + + if not whole_FOV_dir_path or not os.path.exists(whole_FOV_dir_path) or os.path.getsize(whole_FOV_dir_path) == 0: + raise FileNotFoundError(f"The directory or file '{whole_FOV_dir_path}' does not exist or is empty.") + + + stir_sino_path = os.path.join(recons_dir, "stir_sinogram") + stir_sino_img = nib.load(stir_sino_path + ".img") + stir_sinogram = stir_sino_img.get_fdata() + + # Division pixel per pixel, wherever where the pixel value is 0, the result will be 0 + corrected_sinogram = np.where(whole_FOV_stir_sino_data == 0, 0, stir_sinogram / whole_FOV_stir_sino_data) + corrected_img = nib.AnalyzeImage(corrected_sinogram, affine=stir_sino_img.affine, header=stir_sino_img.header) + + #Output dir to corrected sinogram + output_base = os.path.join(recons_dir, "stir_sinogram_corrected") + + #Save .hdr and .img + nib.save(corrected_img, output_base) + + #Copy .img to .s + img_file = output_base + ".img" + s_file = output_base + ".s" + shutil.copy(img_file, s_file) + + #Copy file strir_sinogram.hs + stir_sinogram_hs = stir_sino_path + ".hs" + hs_file_new = output_base + ".hs" + shutil.copy(stir_sinogram_hs, hs_file_new) + + #Open hs_file_new to modify + with open(hs_file_new, 'r') as f: + lines = f.readlines() + + #Modify line where appear stir_sinogram.s to stir_sinogram_corrected.s + lines = [line.replace("stir_sinogram.s", "stir_sinogram_corrected.s") for line in lines] + + #Save changes + with open(hs_file_new, 'w') as f: + f.writelines(lines) + + return (hs_file_new) + + + + + diff --git a/wholebody.py b/wholebody.py index 2c09902a..d0c8c1e5 100644 --- a/wholebody.py +++ b/wholebody.py @@ -1,6 +1,9 @@ import os import yaml import sys +import shutil +import nibabel as nib +import numpy as np from os.path import join, exists from omegaconf import DictConfig, OmegaConf from pyprojroot import here @@ -33,8 +36,8 @@ def __init__(self, cfg: DictConfig): # The following lines will read the general, scanner and config parameters self.sim_type = self.params.get("sim_type") - self.zmin = self.params.get("z_min") - self.zmax = self.params.get("z_max") + self.zmin = float(self.params.get("z_min")) + self.zmax = float(self.params.get("z_max")) # This will load the environment config self.cesga = self.config.get("cesga") @@ -69,37 +72,233 @@ def run(self): os.makedirs(output_dir) log_file = join(output_dir, "logfile.log") + - beds_cs = wb_tools.calculate_center_slices(act_map, self.scanner, self.zmin, self.zmax) + beds_cs = wb_tools.calculate_center_slices(self, act_map, self.scanner, self.zmin, self.zmax) - print("The number of beds to simulate is: %s" % len(beds_cs)) - print("Beds center slides: %s" % beds_cs) + print("\nNumber of simulation to be performed: %s" % len(beds_cs)) #num_beds = len(beds_cs) + print("Beds center slides: %s\n" % beds_cs) - j = 1 + # Global Params + sim_time_original_global = float(self.params.get("simulation_time", 0)) + sim_dose_original_global = float(self.params.get("total_dose", 0)) + self.apply_correction_for_NECR = self.params.get("correction_for_NECR", 0) - for cs in beds_cs: + + # Apply NECR correction if enabled - bed_dir = join(output_dir, "Bed_cs_%s" % cs) - if not exists(bed_dir): - os.makedirs(bed_dir) + # === Correction NECR === + if int(self.apply_correction_for_NECR) == 1: + print("\n>>> Applying NECR correction for this simulation...") + print(f"Original Time of Simulation: {sim_time_original_global} seg") - self.cfg_omega.params.center_slice = int(cs) - self.cfg_omega.params.output_dir = output_name + "/Bed_cs_%s" % cs + # Calculate corrected times (one per bed) + sim_times_per_bed, phantom_doses_FOV = wb_tools.correction_for_NECR(self, act_map, sim_time_original_global) - print("Simulating bed %s with center slice %s" % (j, cs)) + else: + print("\nNECR correction disabled _ keeping original simulation time.") + sim_times_per_bed = [sim_time_original_global] * len(beds_cs) + phantom_doses_FOV = [sim_dose_original_global] * len(beds_cs) + + + # === Run a simulation per bed === + for j, (cs, sim_time_bed, dose ) in enumerate(zip(beds_cs, sim_times_per_bed, phantom_doses_FOV), start=1): + + print(f"\n>>> Running simulation for bed {j} with center slice {cs}, total dose into the FOV: {dose} and time = {sim_time_bed} seg <<<") + + # Create folder for this bed + bed_dir = join(output_dir, f"Bed_{j}_CenterSlice_{cs}") + os.makedirs(bed_dir, exist_ok=True) - bed_simu = SimPET(self.cfg_omega) + # Create copy of configuration file + cfg_copy = OmegaConf.create(OmegaConf.to_container(self.cfg_omega)) + cfg_copy.params.center_slice = int(cs) + cfg_copy.params.simulation_time = sim_time_bed + cfg_copy.params.output_dir = join(output_name, f"Bed_{j}_CenterSlice_{cs}") + + # Run simulation + bed_simu = SimPET(cfg_copy) bed_simu.run() + print("\nAll simulations completed successfully.\n") + + # The following line will be ralated with the simulation of each beds: recons_beds = [] + recons_norm_beds = [] + recons_algorithm = self.scanner.get('recons_type') recons_it = self.scanner.get('numberOfIterations') - for cs in beds_cs: - recons_dir = join(output_dir, "Bed_cs_%s" % cs, "%s_Sim_%s" % (self.sim_type, self.scanner_model), - recons_algorithm) + for i, cs in enumerate(beds_cs, start=1): + recons_dir = join(output_dir, "Bed_%s_CenterSlice_%s" % (i, cs), "%s_Sim_%s" % (self.sim_type, self.scanner_model),recons_algorithm) recons_file = join(recons_dir, 'rec_%s_%s.hdr' % (recons_algorithm, recons_it)) recons_beds.append(recons_file) - joint_beds = join(output_dir, 'rec_%s_%s.hdr' % (recons_algorithm, recons_it)) - wb_tools.join_beds_wb(recons_beds, joint_beds) + + # === Correction for Normalization === + self.Normalization = self.params.get("correction_for_normalization", 0) + + if self.Normalization == 1 and len(recons_beds) > 0: + print("\n>>> Applying NORMALIZATION correction for Reconstructed Image...") + + # Normalization factors for each bed + factors = wb_tools.normalization_factor_correction(self) + + for i, recons_file in enumerate(recons_beds, start=1): + if not os.path.exists(recons_file): + print(f"File not found: {recons_file}") + continue + + if i - 1 >= len(factors): + raise RuntimeError("Normalization factors do not match number of beds") + + print(f"\nApplying normalization correction to bed {i}: {recons_file}") + + factor = factors[i - 1] + print(f">>> Bed {i} normalization factor applied: {factor}") + + # Load reconstructed image + img = nib.load(recons_file) + data = img.get_fdata() + + # Apply normalization + data_norm = data * factor + + # Save normalized image + recons_dir = os.path.dirname(recons_file) + recons_norm_file = join(recons_dir, f"rec_{recons_algorithm}_{recons_it}_norm.hdr") + + img_norm = nib.Nifti1Image(data_norm, img.affine, img.header) + nib.save(img_norm, recons_norm_file) + + recons_norm_beds.append(recons_norm_file) + print(f"Normalized image saved in: {recons_norm_file}") + + # === Joined Beds === + self.joints_beds = self.params.get("joints_beds", 0) + + # Paths for the final joined images + joint_beds = join(output_dir, f"rec_{recons_algorithm}_{recons_it}.hdr") + joint_norm_beds = join(output_dir, f"rec_{recons_algorithm}_{recons_it}_norm.hdr") + + n_beds = len(recons_beds) + n_norm_beds = len(recons_norm_beds) + print(f"Norm_bed = {n_norm_beds}") + + # ---------- NO JOIN (0 or single bed) ---------- + if self.joints_beds == 0 or n_beds == 1: + print("\n>>> No bed joining required (single bed or joints_beds = 0)") + + # Copy the reconstructed image + src_hdr = recons_beds[0] + src_img = src_hdr.replace('.hdr', '.img') # asociado .img + if not os.path.exists(src_hdr) or not os.path.exists(src_img): + raise FileNotFoundError(f"Reconstructed bed files not found: {src_hdr} / {src_img}") + + shutil.copy(src_hdr, joint_beds) + shutil.copy(src_img, joint_beds.replace('.hdr', '.img')) + print(f"Copied single bed image to: {joint_beds}") + + # Copy normalized image if Normalization is enabled + if self.Normalization == 1 and n_beds == 1: + src_norm_hdr = recons_norm_beds[0] + src_norm_img = src_norm_hdr.replace('.hdr', '.img') + shutil.copy(src_norm_hdr, joint_norm_beds) + shutil.copy(src_norm_img, joint_norm_beds.replace('.hdr', '.img')) + print(f"Copied normalized bed image to: {joint_norm_beds}") + + + # ---------- ACTUAL JOIN more than one bed) ---------- + else: + try: + print("\n>>> Joining reconstructed beds...") + wb_tools.join_beds_wb(self, act_map, recons_beds, joint_beds) + print("Bed joining completed.") + except Exception as e: + print(f"Error while joining beds: {e}") + + if self.Normalization == 1 and n_beds > 1: + try: + joint_norm_beds = join(output_dir, f"rec_{recons_algorithm}_{recons_it}_norm.hdr") + wb_tools.join_beds_wb(self, act_map, recons_norm_beds, joint_norm_beds) + + print("Joining of normalized beds completed.") + except Exception as e: + print(f"Error while joining normalized beds: {e}") + + + # === Quantification === + self.quantification = self.params.get("quantification_info", 0) + joint_norm_beds = join(output_dir, f"rec_{recons_algorithm}_{recons_it}_norm.hdr") ## tratar de mencionarlo solo una vez!! #TODO + + try: + #Copy the mask in the same orientation of the simulated image + mask_dir = self.params.get("mask_dirname") + mask_map = join(self.dir_data, mask_dir, self.params.get("mask_map")) + + rotated_mask_file = join(output_dir, "rotated_mask.nii") + rotated_mask = wb_tools.rotate_and_flip_mask(act_map, mask_map, rotated_mask_file, joint_beds, self.zmin, self.zmax) + + + #Changed the dimension of the mash to the same of the simulated image + mask_file = join(output_dir, "mask_image.nii") + reshaped_mask_image = wb_tools.change_act_dimensions(mask_file, rotated_mask, joint_beds) + + #Doing Quantification of Final Image + if self.quantification == 1 and self.joints_beds == 1 and joint_norm_beds is not None: # es si existen los ficheros Arreglar!!!! + + quantification_file = join(output_dir, "Quantification_data.txt") + info_act_map = join(output_dir, "Activity_distribution_data.txt") + + label_file_mask = join(self.dir_data, mask_dir, "mask_map.txt") + label_file_act_map = join(self.dir_data, patient_dir, "act_map.txt") + + import inspect + print(inspect.signature(wb_tools.total_quantification)) + wb_tools.total_quantification(mask_file, joint_norm_beds, quantification_file, label_file_mask) + + wb_tools.distribution_of_dose_into_phantom(self, maps_dir, act_map, info_act_map, label_file_act_map) + + print("Quantification finished") + + except Exception as e: + print(f"Error: Quantification was not performed: {e}") + + + + # #Apply normalization to whole body + # if self.Normalization == 1 and self.joints_beds == 1 and os.path.exists(joint_beds): + # factor_wholeBody = wb_tools.normalization_factor_correction_whole_body(self, joint_beds) + + # # Load reconstructed image + # img = nib.load(joint_beds) + # data = img.get_fdata() + + # # Apply the factor corresponding to each bed + # data_norm = data * factor_wholeBody + + # # Save normalizated image + # recons_dir = os.path.dirname(output_dir) + # recons_norm_wholeBody_file = join(output_dir, 'rec_%s_%s_norm_wholeBody.hdr' % (recons_algorithm, recons_it)) + # img_norm = nib.Nifti1Image(data_norm, img.affine, img.header) + # nib.save(img_norm, recons_norm_wholeBody_file) + + # #Quantification + # quantification_file = join(output_dir, "Quantification_data_whole_body.txt") + + # wb_tools.total_quantification_wholeBody(mask_file, recons_norm_wholeBody_file, quantification_file) + + + # print("Normalization Whole Body completed") + # else: + # print("No Normalization Whole Body completed was performed") + + #wb_tools.coincidencias_mask_vs_image_Claudia(mask_file, act_map) + + + + + + + + From 4d10bfb6c8c27c33db79ceb7bc876e8318c3e4aa Mon Sep 17 00:00:00 2001 From: claudiadom <163831034+claudiadom@users.noreply.github.com> Date: Thu, 23 Jul 2026 13:04:14 +0200 Subject: [PATCH 2/4] Update zOutputSize in discovery_st.yaml reverted --- configs/params/scanner/discovery_st.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/configs/params/scanner/discovery_st.yaml b/configs/params/scanner/discovery_st.yaml index 3e0fb533..e7f02dfd 100644 --- a/configs/params/scanner/discovery_st.yaml +++ b/configs/params/scanner/discovery_st.yaml @@ -34,5 +34,5 @@ add_noise: 0 max_segment: 23 zoomFactor: 1.55 xyOutputSize: 128 -zOutputSize: 160 #47 +zOutputSize: 47 zOutputVoxelSize: 3.27 From bb2f5a76d3529a54bc994ccb029c301747b5a509 Mon Sep 17 00:00:00 2001 From: Claudia Date: Thu, 23 Jul 2026 15:02:30 +0200 Subject: [PATCH 3/4] . --- makefile | 256 +++++++++++++++ utils/Quantification_externo_Claudia.py | 407 ------------------------ utils/wb_tools.py | 86 +---- wholebody.py | 6 +- 4 files changed, 266 insertions(+), 489 deletions(-) create mode 100644 makefile delete mode 100644 utils/Quantification_externo_Claudia.py diff --git a/makefile b/makefile new file mode 100644 index 00000000..8080b1c9 --- /dev/null +++ b/makefile @@ -0,0 +1,256 @@ +SHELL := /bin/bash +TMPDIR := /tmp +ROOT_DIR := $(dir $(abspath $(lastword $(MAKEFILE_LIST)))) +ASSETS_DIR := ${ROOT_DIR}assets +INCLUDE_DIR := ${ROOT_DIR}include +SUBMODULES_DIR := ${ROOT_DIR}submodules + +.PHONY: deps install-simset check-simset clean-simset +.PHONY: install-stir check-stir clean-stir +.PHONY: install-resources check-resources clean-resources +.PHONY: config-git clean-git smudge-filer config-paths clean-paths +.PHONY: dummy-data install clean help + +deps: + sudo apt-get -y -q update ;\ + sudo apt-get install -y -q \ + software-properties-common \ + git[all] \ + git-lfs \ + wget \ + unzip \ + sshpass \ + libboost-dev \ + libboost-all-dev \ + libpcre3 \ + libpcre3-dev \ + libncurses-dev \ + cmake \ + g++ \ + swig ;\ + sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-keys CC86BB64 ;\ + sudo add-apt-repository -y ppa:rmescandon/yq ;\ + sudo apt-get -y -q update ;\ + sudo apt-get install -y -q yq + +SIMSET_SUBMODULE_DIR = ${ROOT_DIR}/submodules/simset_simpet +SIMSET_DEST_DIR = ${INCLUDE_DIR}/SimSET +SIMSET_PATH = ${SIMSET_DEST_DIR}/2.9.2 +SIMSET_BIN = ${SIMSET_PATH}/bin +SIMSET_LIB = ${SIMSET_PATH}/lib +SIMSET_MKALL = ${SIMSET_PATH}/make_all.sh +SIMSET_MKFILE = ${SIMSET_PATH}/make.files/simset.make + +install-simset: + if [ ! -d ${SIMSET_DEST_DIR} ]; then \ + mkdir -p ${SIMSET_PATH} ;\ + cp -r ${SIMSET_SUBMODULE_DIR}/* ${SIMSET_PATH} ;\ + sed -i 's/^\(SIMSET_PATH = \).*$$/\1$(subst /,\/,${SIMSET_PATH})/' ${SIMSET_MKFILE} ;\ + cd ${SIMSET_PATH} && mkdir -p ${SIMSET_LIB} && bash ${SIMSET_MKALL} ;\ + else \ + echo "${SIMSET_DEST_DIR} already exists, run clean-simset if you really want to remove it (you will have to intall SimSET again.)" ;\ + fi + +SIMSET_TAR = ${ASSETS_DIR}/phg.2.9.2.tar.Z +SIMSET_STIR_PATCH = ${ASSETS_DIR}/simset_for_stir.patch + +install-canonic-simset: ${SIMSET_TAR} + if [ ! -d ${SIMSET_DEST_DIR} ]; then \ + mkdir -p ${SIMSET_DEST_DIR} && tar -xvf "${SIMSET_TAR}" --directory=${SIMSET_DEST_DIR} ;\ + cd ${SIMSET_DEST_DIR} && patch -s -p0 < ${SIMSET_STIR_PATCH} ;\ + sed -i 's/^\(SIMSET_PATH = \).*$$/\1$(subst /,\/,${SIMSET_PATH})/' ${SIMSET_MKFILE} ;\ + cd ${SIMSET_PATH} && mkdir -p ${SIMSET_LIB} && bash ${SIMSET_MKALL} ;\ + else \ + echo "${SIMSET_DEST_DIR} already exists, run clean-simset if you really want to remove it (you will have to intall SimSET again.)" ;\ + fi + +check-simset: + declare -a simset_files=(addrandoms bin calcattenuation combinehist makeindexfile phg timesort) ;\ + for file in "$${simset_files[@]}"; do \ + if [ ! -f "${SIMSET_BIN}/$${file}" ]; then \ + echo "${SIMSET_BIN}/$${file} does not exists, check your installation." ;\ + else \ + echo "${SIMSET_BIN}/$${file} exists." ;\ + fi ;\ + done + +clean-simset: + rm -rf ${SIMSET_DEST_DIR} + +SIMSET_SRC = ${SIMSET_PATH}/src +SIMSET_LIBSIMSET = ${SIMSET_LIB}/libsimset.so +STIR_DIR = ${SUBMODULES_DIR}/STIR/STIR +STIR_BUILD_DIR = ${SUBMODULES_DIR}/STIR/build +STIR_INSTALL_DIR = ${SUBMODULES_DIR}/STIR/install +STIR_INSTALL_BIN = ${STIR_INSTALL_DIR}/bin +STIR_MKFILE = ${STIR_BUILD_DIR}/CMakeCache.txt +STIR_FINAL_DEST_DIR = ${INCLUDE_DIR}/STIR +NPROC = $(shell nproc) + +install-stir: install-simset ${STIR_DIR} + if [ ! -d ${STIR_FINAL_DEST_DIR} ]; then \ + mkdir -p ${STIR_BUILD_DIR} ${STIR_INSTALL_DIR} ${STIR_FINAL_DEST_DIR} ;\ + cd ${STIR_BUILD_DIR} && cmake ${STIR_DIR} ;\ + sed -i \ + -e 's/^\(BUILD_SWIG_PYTHON\).*$$/\1:BOOL=OFF/' \ + -e 's/^\(CMAKE_INSTALL_PREFIX\).*$$/\1:PATH=$(subst /,\/,${STIR_INSTALL_DIR})/' \ + -e 's/^\(SIMSET_INCLUDE_DIRS\).*$$/\1:PATH=$(subst /,\/,${SIMSET_SRC})/' \ + -e 's/^\(SIMSET_LIBRARY\).*$$/\1:FILEPATH=$(subst /,\/,${SIMSET_LIBSIMSET})/' \ + -e 's/^\(STIR_OPENMP\).*$$/\1:BOOL=ON/' \ + ${STIR_MKFILE} ;\ + cmake ${STIR_DIR} ;\ + make -s -j${NPROC} ;\ + sudo make install ;\ + mv ${STIR_BUILD_DIR} ${STIR_INSTALL_DIR} ${STIR_FINAL_DEST_DIR} ;\ + else \ + echo "${STIR_FINAL_DEST_DIR} already exists, run clean-stir if you really want to remove it (you will have to intall STIR again)." ;\ + fi + +check-stir: + declare -a stir_files=(FBP2D FBP3DRP forward_project lm_to_projdata OSMAPOSL zoom_image) ;\ + for file in "$${stir_files[@]}"; do \ + if [ ! -f "${STIR_FINAL_DEST_DIR}/install/bin/$${file}" ]; then \ + echo "${STIR_FINAL_DEST_DIR}/install/bin/$${file} does not exists, check your installation." ;\ + else \ + echo "${STIR_FINAL_DEST_DIR}/install/bin/$${file} exists." ;\ + fi ;\ + done + +clean-stir: + rm -rf ${STIR_FINAL_DEST_DIR} + +RESOURCES_ZIP = ${ASSETS_DIR}/fruitcake.zip +RESOURCES_TMP = ${TMPDIR}/resources + +install-resources: ${RESOURCES_ZIP} + mkdir -p ${RESOURCES_TMP} ${INCLUDE_DIR} && unzip -o ${RESOURCES_ZIP} -d ${RESOURCES_TMP} ;\ + declare -a resources=(fruitcake format_converters) ;\ + for rce in "$${resources[@]}"; do \ + include_path="${INCLUDE_DIR}/$${rce}" ;\ + tmp_path="${RESOURCES_TMP}/$${rce}" ;\ + if [ ! -d $${include_path} ]; then \ + mv $${tmp_path} ${INCLUDE_DIR} ;\ + chmod -R +x $${include_path} ;\ + else \ + echo "$${include_path} already exists, run make clean-resources in order to clean resources installation (fruitcake will be removed as well)." ;\ + fi ;\ + done ;\ + rm -rf ${RESOURCES_TMP} + +check-resources: + declare -a resources_paths=(${FRUITCAKE_PATH} ${FORMAT_CONVERTERS_PATH}) ;\ + for path in "$${resources_paths[@]}"; do \ + if [ ! -d "$${path}" ]; then \ + echo "$${path} does not exists, clean fruitcake and format_converters with clean-resources and run install-resources." ;\ + else \ + echo "$${path} exists." ;\ + fi ;\ + done + +clean-resources: + rm -rf ${INCLUDE_DIR}/fruitcake ${INCLUDE_DIR}/format_converters + +config-git: + git config --local filter.config-test.smudge ${ROOT_DIR}scripts/smudge-config_test.sh + git config --local filter.config-test.clean ${ROOT_DIR}scripts/clean-config_test.sh + chmod +x -R ${ROOT_DIR}scripts + +clean-git: + git config --local --unset filter.config-test.smudge + git config --local --unset filter.config-test.clean + +smudge-filter: + rm -f ${ROOT_DIR}configs/config.yaml 2> /dev/null + git checkout HEAD -- ${ROOT_DIR}configs/config.yaml + +FRUITCAKE_PATH = ${INCLUDE_DIR}/fruitcake +FRUITCAKE_BIN = ${FRUITCAKE_PATH}/bin +FRUITCAKE_LIB = ${FRUITCAKE_PATH}/book/lib +FORMAT_CONVERTERS_PATH = ${INCLUDE_DIR}/format_converters + +config-paths: + touch "$${HOME}"/.bashrc ;\ + declare -a simpet_paths=( \ + 'export PATH=${FRUITCAKE_BIN}:$$PATH' \ + 'export LD_LIBRARY_PATH=${FRUITCAKE_LIB}:$$LD_LIBRARY_PATH' \ + 'export PATH=${FORMAT_CONVERTERS_PATH}:$$PATH' \ + ) ;\ + for path in "$${simpet_paths[@]}"; do \ + grep -qxF "$${path}" "$${HOME}"/.bashrc || echo "$${path}" >> "$${HOME}"/.bashrc ;\ + done + +clean-paths: + touch "$${HOME}/.bashrc" + sed -i \ + -e '/export PATH=$(subst /,\/,${FRUITCAKE_BIN}):$$PATH/d' \ + -e '/export LD_LIBRARY_PATH=$(subst /,\/,${FRUITCAKE_LIB}):$$LD_LIBRARY_PATH/d' \ + -e '/export PATH=$(subst /,\/,${FORMAT_CONVERTERS_PATH}):$$PATH/d' \ + "$${HOME}/.bashrc" + +DATA_DIR = ${ROOT_DIR}Data +DATA_ZIP = ${ASSETS_DIR}/Data.zip + +dummy-data: ${DATA_ZIP} + if [ ! -d "${DATA_DIR}" ]; then \ + mkdir -p ${DATA_DIR} && unzip -o ${DATA_ZIP} -d ${ROOT_DIR} ;\ + else \ + echo "${DATA_DIR} already exists, remove it manually." ;\ + fi + +install: + ${MAKE} \ + deps \ + install-simset \ + check-simset \ + install-stir \ + check-stir \ + install-resources \ + check-resources \ + config-paths \ + +clean: + ${MAKE} \ + clean-simset \ + clean-stir \ + clean-resources + +help: + @echo "Help:" + @echo " - deps: Install the dependencies of the projects via apt." + @echo "" + @echo " - install-simset: Install tweaked version of SimSET. + @echo "" + @echo " - install-canonic-simset: Install non-tweaked version of SimSET with STIR patch at directory ${SIMSET_DEST_DIR}." + @echo "" + @echo " - check-simset: Check that SimSET binaries exist." + @echo "" + @echo " - clean-simset: Clean SimSET installation (removes ${SIMSET_DEST_DIR} directory)." + @echo "" + @echo " - install-stir: Install STIR at directory $(shell dirname ${STIR_DIR}), install-simset is a prerequisite. Set NPROC=n to use n CPU cores in compilation." + @echo "" + @echo " - check-stir: Check that STIR binaries exist." + @echo "" + @echo " - clean-stir: Clean STIR installation removing ${STIR_INSTALL_DIR} and ${STIR_BUILD_DIR} directories." + @echo "" + @echo " - install-resources: Decompress ${RESOURCES_ZIP} and move fruitcake and format_converters to ${INCLUDE_DIR}" + @echo "" + @echo " - check-resources: Checks that fruitcake and format_converters are in ${INCLUDE_DIR}." + @echo "" + @echo " - clean-resources: Remove fruitcake and format_converters from ${INCLUDE_DIR}." + @echo "" + @echo " - config-git: Add project filter drivers to local git configuration and make them executable." + @echo "" + @echo " - clean-git: Remove project filter drivers from local git configuration." + @echo "" + @echo " - smudge-filer: Apply smudge filter to config file (updates paths)." + @echo "" + @echo " - config-paths: Iff not present in .bashrc, the paths of the projects will be appended to the file. If ~/.bashrc does not exists it will be created." + @echo "" + @echo " - clean-paths: If the paths of the project are present in .bashrc file, they will be deleted. If ~/.bashrc does not exists it will be created." + @echo "" + @echo " - dummy-data: Uncompress ${DATA_ZIP} at ${DATA_DIR} if ${DATA_DIR} does not exists." + @echo "" + @echo " - install: Run all install, config and check recipes." + @echo "" + @echo " - clean: Run all clean recipes." + diff --git a/utils/Quantification_externo_Claudia.py b/utils/Quantification_externo_Claudia.py deleted file mode 100644 index 6227bdfc..00000000 --- a/utils/Quantification_externo_Claudia.py +++ /dev/null @@ -1,407 +0,0 @@ - -import sys -import math #new -from os.path import join, dirname -import nibabel as nib -import numpy as np -import os -from scipy.ndimage import zoom -from scipy.ndimage import gaussian_filter -from scipy.ndimage import median_filter -from pyprojroot import here - -sys.path.append(str(here())) -from utils import resources as rsc -from utils import spm_tools as spm -from utils import tools -from src.stir import stir_tools - -#Real Distribution -act_map = - -#Organs diferences for region with differents values -act_map_dif_region = -#Changed of dimensions -ct_image = join(output_dir, "Ct_image_act.img") -ct_image_2 = join(output_dir, "CT_image_act.img") -#address of map -rec_OSEM3D_48_norm.hdr -rec_OSEM3D_48_norm_wholeBody.hdr - -total_quantification = join(output_dir, 'rec_%s_%s_norm.hdr' % (recons_algorithm, recons_it)) -quantification_file = join(output_dir, "Quantification_data.txt") -info_act_map = - - -def rotate_image(act_map, ct_image_act): - # Upload original image (corresponding to the activity map) - act_img = nib.load(act_map) - #act_data = act_img.get_fdata() - act_data = act_img.dataobj[:] # conserva enteros y evita decimales - - - print(f"Image_before_rotate:{act_img.affine}") - - act_header = act_img.header - act_shape = np.array(act_data.shape[:3]) - act_voxel_sizes = act_img.header.get_zooms()[:3] - - # Rotation and flip - rot_xy = np.rot90(act_data, k=2, axes=(0,1)) - final = np.flip(rot_xy, axis=0) - - - # Save final image - final_img = nib.Nifti1Image(final, act_img.affine, act_img.header) - - - output_base = os.path.splitext(ct_image_act)[0] - if not output_base.endswith(".img") and not output_base.endswith(".hdr"): - output_base = output_base # base name without extension - - output_path = output_base + ".img" # nibabel will also generate the .hdr - nib.save(final_img, output_path) - - return final - -def change_act_dimensions(ct_image_act_2, ct_image_act, joint_beds): - """ - Adjust a CT image to match reference voxel size and shape, padding with zeros if necessary. - """ - # --- Load original image --- - img = nib.load(ct_image_act) - data = img.dataobj[:] # conserva enteros y evita decimales - - if data.ndim == 4 and data.shape[3] == 1: - data = np.squeeze(data, axis=3) - - - shape = np.array(data.shape[:3]) - voxel_size = np.array(img.header.get_zooms()[:3]) - size_mm = shape * voxel_size - - # --- Load reference image --- - recons_img = nib.load(joint_beds) - recons_shape = np.array(recons_img.shape[:3]) - recons_voxel_sizes = np.array(recons_img.header.get_zooms()[:3]) - - - # --- Adjust resolution --- - new_voxel_size = recons_voxel_sizes - new_shape = np.round(size_mm / new_voxel_size).astype(int) - scale_factor = new_shape / shape - - - print("Original Dimensions:", data.shape) - print("Scale Factor:", scale_factor) - - # --- Reescalar volumen (sin interpolación continua) --- - new_data = zoom(data, zoom=scale_factor, order=0, prefilter=False) - - # --- Corrige origen geométrico para que el centro quede alineado --- - orig_center = (shape * voxel_size) / 2 - new_center = (new_shape * new_voxel_size) / 2 - shift = orig_center - new_center # desplazamiento en mm - - # --- Nueva cabecera y affine --- - new_affine = recons_img.affine.copy() - new_affine[:3, :3] = np.diag(new_voxel_size) - new_affine[:3, 3] = 0 # centra la imagen correctamente - - # --- Create header and affine --- - nuevo_header = nib.Nifti1Header() - nuevo_header.set_data_shape(new_data.shape) - nuevo_header.set_zooms(new_voxel_size) - - # --- Create intermediate image --- - final_image = nib.Nifti1Image(new_data, affine=new_affine, header=nuevo_header) - - # --- Symmetric padding to match reference shape --- - final_shape = np.array(final_image.shape) - diff = recons_shape - final_shape[:3] - - pad_width = [] - cropped_data = final_image.get_fdata() - - for i in range(3): - if diff[i] >= 0: - # La imagen es más pequeña → agregamos ceros - pad_width.append((diff[i] // 2, diff[i] - diff[i] // 2)) - else: - # La imagen es más grande → recortamos - start = abs(diff[i]) // 2 - end = start + recons_shape[i] - cropped_data = np.take(cropped_data, indices=range(start, end), axis=i) - pad_width.append((0, 0)) - - padded_data = np.pad(final_image.get_fdata(), pad_width, mode='constant', constant_values=0) - - # --- Ajuste final de tamaño exacto --- - # En algunos casos, el padding o recorte previo deja 1 voxel de diferencia por redondeos - final_data = padded_data - for axis in range(3): - if final_data.shape[axis] > recons_shape[axis]: - start = (final_data.shape[axis] - recons_shape[axis]) // 2 - end = start + recons_shape[axis] - final_data = np.take(final_data, indices=range(start, end), axis=axis) - elif final_data.shape[axis] < recons_shape[axis]: - pad_before = (recons_shape[axis] - final_data.shape[axis]) // 2 - pad_after = recons_shape[axis] - final_data.shape[axis] - pad_before - final_data = np.pad( - final_data, - [(pad_before, pad_after) if ax == axis else (0, 0) for ax in range(3)], - mode='constant', - constant_values=0 - ) - # --- Create final NIfTI image --- - final_padded_image = nib.Nifti1Image(final_data, affine=new_affine, header=nuevo_header) - final_padded_image.set_qform(new_affine, code=1) - final_padded_image.set_sform(new_affine, code=1) - - - # --- Save final image --- - nib.save(final_padded_image, ct_image_act_2) - - # --- Delete original image --- - base_name = os.path.splitext(ct_image_act)[0] # quita la extensión .img - for f in [base_name + ".img", base_name + ".hdr"]: - if os.path.exists(f): - os.remove(f) - - return final_padded_image - -def total_quantification(ct_image_act_2, joint_norm_beds, quantification_file): - """ - Compute information of target image per labeled region in reference image and save to a TXT file. - Obtains information about the simulated distribution of activity in the phantom and saves it in a TXT file. - """ - # --- Load images --- - target_img = nib.load(joint_norm_beds) - roi_img = nib.load(ct_image_act_2) - - - target_data = target_img.get_fdata() - roi_data = roi_img.get_fdata() - - - # --- Ensure same shape --- - if target_data.shape != roi_data.shape: - raise ValueError("Target and reference images must have the same shape.") - - # --- Get all labels except 0 (background) --- - labels = np.unique(roi_data) - labels = labels[labels != 0] - - # --- Calculate mean per label --- - mean_values = [target_data[roi_data == label].mean() for label in labels] #KBq/cc - - #--- Calculate Volumen per label --- - voxel_counts = [np.sum(roi_data == label) for label in labels] - dx, dy, dz = np.array(roi_img.header.get_zooms()[:3]) - voxel_volume = dx * dy * dz - - volumes = [(count * voxel_volume)/1000 for count in voxel_counts] - - # --- Calculate activity per region --- - activity_region_KBq = [m * v for m, v in zip(mean_values, volumes)] - activity_region_mCi = [a / (3.7 * 10**4) for a in activity_region_KBq] - mean_values_mCi = [b / c for b, c in zip(activity_region_mCi, volumes)] - - # --- Calculate TOTAL activity across all labels --- - total_activity_KBq = float(np.sum(activity_region_KBq)) - total_activity_mCi = float(np.sum(activity_region_mCi)) - - # --- Save report as TXT --- - with open(quantification_file, 'w') as f: - # Encabezado - f.write( - f"{'Label':<10}{'MeanValue (mCi/cc)':>20}{'MeanValue (KBq/cc)':>20}{'Pixels':>15}" - f"{'Vol (ccm³)':>15}{'Act (mCi)':>15}{'Act (KBq)':>15}\n" - ) - - # Filas por cada label - for label, mean_val_mCi, mean_val_KBq, count, vol, act_mCi, act_KBq in zip( - labels, mean_values_mCi, mean_values, voxel_counts, volumes, activity_region_mCi, activity_region_KBq - ): - f.write( - f"{int(label):<10}{mean_val_mCi:>20.4f}{mean_val_KBq:>20.4f}{count:>15.0f}" - f"{vol:>15.4f}{act_mCi:>15.4f}{act_KBq:>15.4f}\n" - ) - # Linea separadora - f.write("="*90 + "\n") - - # Linea TOTAL - f.write( - f"{'TOTAL':<10}{'':>20}{'':>15}{'':>15}" - f"{total_activity_mCi:>15.4f}{total_activity_KBq:>15.4f}\n" - ) - - return {"labels": labels, "mean_values_mCi": mean_values_mCi, "mean_values_KBq": mean_values, "pixel_counts": voxel_counts, - "volumes": volumes, "activity_mCi": activity_region_mCi, "activity_KBq": activity_region_KBq, - "total_activity_mCi": total_activity_mCi, "total_activity_KBq": total_activity_KBq - } - - - """ - Compute information of target image per labeled region in reference image and save to a TXT file. - Obtains information about the simulated distribution of activity in the phantom and saves it in a TXT file. - """ - # --- Load images --- - target_img = nib.load(recons_norm_wholeBody_file) - roi_img = nib.load(ct_image_act_2) - - - target_data = target_img.get_fdata() - roi_data = roi_img.get_fdata() - - - # --- Ensure same shape --- - if target_data.shape != roi_data.shape: - raise ValueError("Target and reference images must have the same shape.") - - # --- Get all labels except 0 (background) --- - labels = np.unique(roi_data) - labels = labels[labels != 0] - - # --- Calculate mean per label --- - mean_values_KBq = [target_data[roi_data == label].mean() for label in labels] - mean_values_mCi = mean_values_mCi = [v * 0.00002703 for v in mean_values_KBq] #mCi - - #--- Calculate Volumen per label --- - voxel_counts = [np.sum(roi_data == label) for label in labels] - dx, dy, dz = np.array(roi_img.header.get_zooms()[:3]) - voxel_volume = dx * dy * dz - - volumes = [(count * voxel_volume)/1000 for count in voxel_counts] - - # --- Calculate activity per region --- - activity_region_KBq = [m * v for m, v in zip(mean_values_KBq, volumes)] - activity_region_mCi = [a / (3.7 * 10**4) for a in activity_region_KBq] - - # --- Calculate TOTAL activity across all labels --- - total_activity_KBq = float(np.sum(activity_region_KBq)) - total_activity_mCi = float(np.sum(activity_region_mCi)) - - - - # --- Save report as TXT --- - with open(quantification_file, 'w') as f: - # Encabezado - f.write( - f"{'Label':<10}{'MeanValue (mCi/cc)':>20}{'MeanValue (KBq/cc)':>20}{'Pixels':>15}" - f"{'Vol (ccm³)':>15}{'Act (mCi)':>15}{'Act (KBq)':>15}\n" - ) - - # Filas por cada label - for label, mean_val_mCi, mean_val_KBq, count, vol, act_mCi, act_KBq in zip( - labels,mean_values_KBq, mean_values_mCi, voxel_counts, volumes, activity_region_mCi, activity_region_KBq - ): - f.write( - f"{int(label):<10}{mean_val_mCi:>20.4f}{mean_val_KBq:>20.4f}{count:>15.0f}" - f"{vol:>15.4f}{act_mCi:>15.4f}{act_KBq:>15.4f}\n" - ) - # Linea separadora - f.write("="*90 + "\n") - - # Linea TOTAL - f.write( - f"{'TOTAL':<10}{'':>20}{'':>15}{'':>15}" - f"{total_activity_mCi:>15.4f}{total_activity_KBq:>15.4f}\n" - ) - - - return {"labels": labels, "mean_values_mCi": mean_values_mCi, "mean_values_KBq": mean_values_mCi, "pixel_counts": voxel_counts, - "volumes": volumes, "activity_mCi": activity_region_mCi, "activity_KBq": activity_region_KBq, - "total_activity_mCi": total_activity_mCi, "total_activity_KBq": total_activity_KBq - } - -def distribution_of_dose_into_phantom(act_map, info_act_map, self): - """ - Computes the actual distribution of activity by region labeled in act_map. - Obtains information about the actual distribution of activity in the phantom and saves it in a TXT file. - """ - - #Upload the act image - phantom_act = nib.load(act_map) - phantom_act_data = phantom_act.get_fdata() - phantom_voxel_volumen = abs(np.prod(phantom_act.header.get_zooms()[:3]) / 1000) #cm^3 - total_phantom_counts = np.sum(phantom_act_data) - - #Calculate the factor to convert to real phantom activity - if self.sim_dose != 0: - phantom_dose = abs(total_phantom_counts * phantom_voxel_volumen) # uCi - print(f"phantom_dose: {phantom_dose} au") - - act_table_factor = self.sim_dose * 1000 / phantom_dose - else: - act_table_factor = 1 - - print(f"Factor_actividad: {act_table_factor}") - print(f"self.sim_dose: {self.sim_dose} mCi") - print(f"total_phantom_counts: {total_phantom_counts}") - - phantom_real_act = phantom_act_data * (act_table_factor/1000) #se multilica por mil para llevar los valores a mCi - - # Save the file into the same Folder of act_map with different name - base_dir = os.path.dirname(act_map) - name, ext = os.path.splitext(os.path.basename(act_map)) - - if ext.lower() in [".hdr", ".img"]: - output_path = os.path.join(base_dir, f"{name}_real.hdr") - else: - output_path = os.path.join(base_dir, f"{name}_real.nii") - - #New image with the real activity to simulate - new_img = nib.Nifti1Image(phantom_real_act, affine=phantom_act.affine, header=phantom_act.header) - # Save new image - nib.save(new_img, output_path) - - - # Get information about the new image - - # --- Calculate average values per label --- - labels = np.unique(phantom_act_data) - labels = labels[labels != 0] - - mean_values_mCi = [phantom_real_act[phantom_act_data == label].mean() for label in labels] - voxel_counts = [np.sum(phantom_act_data == label) for label in labels] - volumes_region = [count * phantom_voxel_volumen for count in voxel_counts] - - print(f"voxel_counts: {voxel_counts}") - print(f"volumes_region: {volumes_region}") - print(f"phantom_volumen_voxel:{phantom_voxel_volumen}") - - # --- Calculate activity per region --- - activity_region_mCi = [m * v for m, v in zip(mean_values_mCi, volumes_region)] - activity_region_KBq = [a * 3.7 * 10**4 for a in activity_region_mCi] - mean_values_KBq = [b / c for b, c in zip(activity_region_KBq , volumes_region)] - - # --- Calculate TOTAL activity across all labels --- - total_activity_KBq = float(np.sum(activity_region_KBq)) - total_activity_mCi = float(np.sum(activity_region_mCi)) - - print(f"activity_region_KBq: {activity_region_KBq}") - # --- Save report as TXT --- - - with open(info_act_map, 'w') as f: - f.write(f"{'Label':<10}{'MeanValue (mCi/cc)':>20}{'MeanValue (KBq/cc)':>20}{'Pixels':>15}{'Vol (ccm³)':>15}{'Act (mCi)':>15}{'Act (KBq)':>15}\n") - for label, mean_val_mCi, mean_val_KBq, count, vol, act_mCi, act_KBq in zip( - labels, mean_values_mCi, mean_values_KBq, voxel_counts, volumes_region, activity_region_mCi, activity_region_KBq - ): - f.write(f"{int(label):<10}{mean_val_mCi:>20.4f}{mean_val_KBq:>20.4f}{count:>15.0f}{vol:>15.4f}{act_mCi:>15.4f}{act_KBq:>15.4f}\n") - - - # Linea separadora - f.write("="*90 + "\n") - - # Linea TOTAL - f.write( - f"{'TOTAL':<10}{'':>20}{'':>15}{'':>15}" - f"{total_activity_mCi:>15.4f}{total_activity_KBq:>15.4f}\n" - ) - - return { - "image_path": output_path, "factor": act_table_factor, - "labels": labels, "mean_values_mCi": mean_values_mCi,"mean_values_KBq": mean_values_KBq,"volumes": volumes_region, - "activity_mCi": act_mCi, "activity_KBq": act_KBq, "total_activity_mCi": total_activity_mCi, "total_activity_KBq": total_activity_KBq - } diff --git a/utils/wb_tools.py b/utils/wb_tools.py index 2bfaf6ae..28ad9977 100644 --- a/utils/wb_tools.py +++ b/utils/wb_tools.py @@ -386,7 +386,7 @@ def calculate_map_into_fov(self, act_map): map_into_FOV_start = int(round(cs - half_fov_slices)) map_into_FOV_end = int(round(cs + half_fov_slices)) - #Eliminar + #Remove #print(f"Before correction: start={map_into_FOV_start}, end={map_into_FOV_end}") if map_into_FOV_start < zmin: @@ -451,7 +451,7 @@ def correction_for_NECR(self, act_map, sim_time_original): print(f"Corrected simulation time = {sim_time_bed} seg") - # ####BORRAR + # ####Remove # #labels into each beds # # Extract only the portion of the map within the FOV (Z-axis) # # Assuming that axis 2 (index 2) is the Z-axis @@ -637,7 +637,7 @@ def normalization_factor_correction_whole_body(self, joint_beds): # Valor medio ignorando ceros #mean_nonzero = sum_nonzero / num_nonzero_voxels - #Calcular la concentracion estimada segun el total sum de la imagen + #Calculate the estimated concentration based on the total sum of the image #conc_estimate_image = 1.34*10**-5 * sum_nonzero - 8.96*10**-9 act_estimate_image = float(4.22*10**-4 * sum_nonzero - 4.65*10**-8) #volumen de la voi en cmm³ @@ -1011,76 +1011,6 @@ def total_quantification(mask_file, joint_norm_beds, quantification_file, label_ } - """ - Compute information of target image per labeled region in reference image and save to a TXT file. - Obtains information about the simulated distribution of activity in the phantom and saves it in a TXT file. - """ - # --- Load images --- - target_img = nib.load(joint_norm_beds) - roi_img = nib.load(ct_image_act_2) - - - target_data = target_img.get_fdata() - roi_data = roi_img.get_fdata() - - - # --- Ensure same shape --- - if target_data.shape != roi_data.shape: - raise ValueError("Target and reference images must have the same shape.") - - # --- Get all labels except 0 (background) --- - labels = np.unique(roi_data) - labels = labels[labels != 0] - - # --- Calculate mean per label --- - mean_values = [target_data[roi_data == label].mean() for label in labels] #KBq/cc - - #--- Calculate Volumen per label --- - voxel_counts = [np.sum(roi_data == label) for label in labels] - dx, dy, dz = np.array(roi_img.header.get_zooms()[:3]) - voxel_volume = dx * dy * dz - - volumes = [(count * voxel_volume)/1000 for count in voxel_counts] - - # --- Calculate activity per region --- - activity_region_KBq = [m * v for m, v in zip(mean_values, volumes)] - activity_region_mCi = [a / (3.7 * 10**4) for a in activity_region_KBq] - mean_values_mCi = [b / c for b, c in zip(activity_region_mCi, volumes)] - - # --- Calculate TOTAL activity across all labels --- - total_activity_KBq = float(np.sum(activity_region_KBq)) - total_activity_mCi = float(np.sum(activity_region_mCi)) - - # --- Save report as TXT --- - with open(quantification_file, 'w') as f: - # Encabezado - f.write( - f"{'Label':<10}{'MeanValue (mCi/cc)':>20}{'MeanValue (KBq/cc)':>20}{'Pixels':>15}" - f"{'Vol (ccm³)':>15}{'Act (mCi)':>15}{'Act (KBq)':>15}\n" - ) - - # Filas por cada label - for label, mean_val_mCi, mean_val_KBq, count, vol, act_mCi, act_KBq in zip( - labels, mean_values_mCi, mean_values, voxel_counts, volumes, activity_region_mCi, activity_region_KBq - ): - f.write( - f"{int(label):<10}{mean_val_mCi:>20.4f}{mean_val_KBq:>20.4f}{count:>15.0f}" - f"{vol:>15.4f}{act_mCi:>15.4f}{act_KBq:>15.4f}\n" - ) - # Linea separadora - f.write("="*90 + "\n") - - # Linea TOTAL - f.write( - f"{'TOTAL':<10}{'':>20}{'':>15}{'':>15}" - f"{total_activity_mCi:>15.4f}{total_activity_KBq:>15.4f}\n" - ) - - return {"labels": labels, "mean_values_mCi": mean_values_mCi, "mean_values_KBq": mean_values, "pixel_counts": voxel_counts, - "volumes": volumes, "activity_mCi": activity_region_mCi, "activity_KBq": activity_region_KBq, - "total_activity_mCi": total_activity_mCi, "total_activity_KBq": total_activity_KBq - } - def total_quantification_wholeBody(mask_file, recons_norm_wholeBody_file, quantification_file): #TODO """ Compute information of target image per labeled region in reference image and save to a TXT file. @@ -1267,13 +1197,13 @@ def distribution_of_dose_into_phantom(self, maps_dir, act_map, info_act_map, lab "total_activity_KBq": total_activity_KBq } -def coincidencias_mask_vs_image_Claudia(mask_file, act_map): #ELIMINAR ES SOLO PARA COMPROBACIÓN - # --- Verificación de volúmenes por label --- +def coincidences_mask_vs_image(mask_file, act_map): #Remove it is only a check + # --- Verification of the label volumen --- import numpy as np import nibabel as nib try: - orig_img = nib.load(act_map) # imagen original base + orig_img = nib.load(act_map) # original image new_img = nib.load(mask_file) orig_data = orig_img.get_fdata() @@ -1286,7 +1216,7 @@ def coincidencias_mask_vs_image_Claudia(mask_file, act_map): #ELIMINAR ES SOLO labels = labels[labels != 0] - print("\n--- Comparación de volúmenes por región ---") + print("\n--- Comparation of volumen per label ---") for label in labels: count_orig = np.sum(orig_data == label) count_new = np.sum(new_data == label) @@ -1296,7 +1226,7 @@ def coincidencias_mask_vs_image_Claudia(mask_file, act_map): #ELIMINAR ES SOLO print(f"Label {int(label):3d}: original={vol_orig:.3f} cm³, new={vol_new:.3f} cm³, ratio={ratio:.3f}") except Exception as e: - print(f"[Aviso] No se pudo comparar volúmenes por label: {e}") + print(f"Notice: Could not compare volumes by label: {e}") def total_fov_correction(self, recons_dir): diff --git a/wholebody.py b/wholebody.py index d0c8c1e5..4a940da8 100644 --- a/wholebody.py +++ b/wholebody.py @@ -1,6 +1,4 @@ -import os -import yaml -import sys +import os, yaml, sys import shutil import nibabel as nib import numpy as np @@ -293,7 +291,7 @@ def run(self): # else: # print("No Normalization Whole Body completed was performed") - #wb_tools.coincidencias_mask_vs_image_Claudia(mask_file, act_map) + #wb_tools.coincidences_mask_vs_image(mask_file, act_map) From ed85d46d2198466adbc97c41de48f86d498c1545 Mon Sep 17 00:00:00 2001 From: Claudia Date: Wed, 5 Aug 2026 18:36:46 +0200 Subject: [PATCH 4/4] Added: Save the configuration for the config entry at the specified path, as well as the configuration for each bed for whole-body simulation --- src/simset/simset_sim.py | 26 +----- utils/tools.py | 4 +- utils/wb_tools.py | 170 +-------------------------------------- wholebody.py | 54 +++---------- 4 files changed, 17 insertions(+), 237 deletions(-) diff --git a/src/simset/simset_sim.py b/src/simset/simset_sim.py index 25f4cf87..39d5d21c 100644 --- a/src/simset/simset_sim.py +++ b/src/simset/simset_sim.py @@ -6,7 +6,7 @@ import numpy as np import nibabel as nib import warnings -import math #añadido para Bruker PET/MRI PRECLINICA +import math from multiprocessing import Process from pathlib import Path from os import PathLike @@ -85,7 +85,7 @@ def __init__( self.phglistmode = params.get("phglistmode") self.add_randoms = params.get("add_randoms") - #TODO #correcton remove to the code + @@ -144,28 +144,6 @@ def run_simset_simulation(self, sim_dir): act_table_factor = self.sim_dose * 1000 / phantom_dose else: act_table_factor = 1 - - #Borrer a partir de aqui - print(f"act_table_factor_value:{act_table_factor}") - - img = nib.load(self.act_map) - data = img.get_fdata() - # Obtener los labels únicos (excluyendo el 0 si es fondo) - labels, counts = np.unique(data, return_counts=True) - - # Eliminar el fondo (label 0) si lo hay - mask = labels != 0 - labels = labels[mask] - counts = counts[mask] - - # Mostrar resultado - for label, count in zip(labels, counts): - print(f"Label {int(label)}: {int(count)} voxeles") - - # Mostrar cantidad total de labels - print(f"\nTotal de labels distintos: {len(labels)}") - - #Hasta la linea de alante borrar # Creates the data files from the simulation maps act_img = self.act_map[0:-3] + "img" diff --git a/utils/tools.py b/utils/tools.py index 6f411d04..1e0317c9 100644 --- a/utils/tools.py +++ b/utils/tools.py @@ -544,9 +544,7 @@ def reorient_dcmtonii(image_path): # shutil.copy(join(path,"f"+image), image_path) os.remove(join(path,"c"+image)) # os.remove(join(path,"f"+image)) - - - + def petmr2maps(pet_image, mri_image, ct_image, log_file, spm_run, output_dir, mode="SimSET"): """ It will create act and att maps from PET and MR images. diff --git a/utils/wb_tools.py b/utils/wb_tools.py index 28ad9977..f55f9a36 100644 --- a/utils/wb_tools.py +++ b/utils/wb_tools.py @@ -198,7 +198,6 @@ def pet_to_actmap(self): rcommand = '%s %s %s 1B >> %s' % (change_format, act_out, act_out, self.log) tools.osrun(rcommand, self.log) - def calculate_center_slices(self, act_map, scanner, zmin, zmax, overlapping=0.1): """Calculate the center slices of beds for the given axial range.""" @@ -254,7 +253,6 @@ def calculate_center_slices(self, act_map, scanner, zmin, zmax, overlapping=0.1) self.beds_cs = beds_cs return beds_cs - def update_act_map(spmrun, act_map, att_map, orig_pet, simu_pet, output): output_dir = dirname(output) mfile = join(dirname(output), "fusion.m") @@ -343,7 +341,6 @@ def update_act_map(spmrun, act_map, att_map, orig_pet, simu_pet, output): updated_act_img = nib.AnalyzeImage(updated_act, simpet.affine, simpet.header) nib.save(updated_act_img, output) - def cut_image_min_max_slices(input_img, min_slice, max_slice, output): img = nib.load(input_img) img_data = tools.fix_4d_data(img.get_fdata()) @@ -572,7 +569,7 @@ def normalization_factor_correction(self): # Convert mCi to kBq phantom_dose_KBq = dose_mCi * 3.7e4 sim_time_original_global = float(self.params.get("simulation_time", 0)) - print(f"Phantom dose in KBq: {phantom_dose_KBq}") #borrar + print(f"Phantom dose in KBq: {phantom_dose_KBq}") # Linear adjustment equation #lineal_ecuac = 5.9e-5 * phantom_dose_KBq + 0.978 @@ -586,7 +583,7 @@ def normalization_factor_correction(self): value_Q_norm = float(0.062 * (phantom_dose_KBq) + 986.39) value_Q_norm_corregido_time = value_Q_norm / (sim_time_original_global/300) # 300 seg es el tiempo de adq cyl calibración - print(f"Normalization value per bed: {value_Q_norm}") #borrar + #print(f"Normalization value per bed: {value_Q_norm}") #results.append(value_Q_norm) results.append(value_Q_norm_corregido_time) @@ -598,89 +595,6 @@ def normalization_factor_correction(self): print("Normalization factors by FOV:", [round(x, 3) for x in results]) return results -def normalization_factor_correction_whole_body(self, joint_beds): - factor_Q_norm = 940.34 - results = [] - sim_dose_original_global = float(self.params.get("total_dose", 0)) - print(f"Phantom dose (FOV):{sim_dose_original_global}") - - image_sim = nib.load(joint_beds) - image_sim_data = image_sim.get_fdata() - - # Reemplazar NaN o inf por 0 (solo si hay) - image_sim_data = np.nan_to_num(image_sim_data, nan=0.0, posinf=0.0, neginf=0.0) - - # Calcular volumen del voxel - dx, dy, dz = np.array(image_sim.header.get_zooms()[:3]) - voxel_volume = (dx * dy * dz)/1000 #cm^3 - - - # Calcular sumas - total_voxel_sum = np.sum(image_sim_data) - Activity_image = (voxel_volume * total_voxel_sum)/10 - - # Calcular valor medio - # Numero total de voxeles - #num_voxels = image_sim_data.size - - # Valor medio de todos los voxeles - #mean_voxel_value = total_voxel_sum / num_voxels - - # Creamos una máscara para voxeles distintos de cero - mask = image_sim_data != 0 - # Sumamos solo los voxeles distintos de cero - sum_nonzero = np.sum(image_sim_data[mask]) - # Contamos cuántos voxeles son distintos de cero - num_nonzero_voxels = np.count_nonzero(image_sim_data) - volumen_non_zero = voxel_volume * num_nonzero_voxels - print(f"Volumen_non_zero:{volumen_non_zero} cc") - # Valor medio ignorando ceros - #mean_nonzero = sum_nonzero / num_nonzero_voxels - - #Calculate the estimated concentration based on the total sum of the image - - #conc_estimate_image = 1.34*10**-5 * sum_nonzero - 8.96*10**-9 - act_estimate_image = float(4.22*10**-4 * sum_nonzero - 4.65*10**-8) #volumen de la voi en cmm³ - conc_estimate_image = float(act_estimate_image/98.96) #volumen de la voi en cmm³ - #print(mean_nonzero) - print(f"Total_suma:{sum_nonzero} au/cc") - print(f"Concentracion estimada:{conc_estimate_image} au/cc Actividad estimada: {act_estimate_image}") - - print(f"Voxel volume: {voxel_volume:.3f} mm³") - print(f"Total_sum : {total_voxel_sum}") - print(f"Activity_image : {Activity_image}") - - - # Convert mCi to kBq - phantom_dose_KBq = sim_dose_original_global * 3.7*10**4 - #mean_nonzero = phantom_dose_KBq / num_nonzero_voxels - - #print(f"Phantom mean dose in KBq: {mean_nonzero}") - print(f"Phantom dose in KBq: {phantom_dose_KBq}") - - # Linear adjustment equation - #lineal_ecuac = (0.0038 * (phantom_dose_KBq / 31.43) + 0.9725) - lineal_ecuac = (0.0000606 * phantom_dose_KBq + 0.975) - print(f"Lineal Ecuation: {lineal_ecuac}") - - #value_Q_norm = float(lineal_ecuac * conc_estimate_image) - #value_Q_norm = float(lineal_ecuac * conc_estimate_image * factor_Q_norm) - value_Q_norm = float((lineal_ecuac * factor_Q_norm)/1.47) - #value_Q_norm = float(factor_Q_norm * lineal_ecuac ) - #value_Q_norm = float(lineal_ecuac) - - print(f"Normalization value per total: {value_Q_norm}") - results.append(value_Q_norm) - - - - if len(results) == 0: - print(f"There is no FOV dose. Using factor = 1.0 by default") - return 1.0 - - print("Normalization factors by FOV:", [round(x, 3) for x in results]) - return results - def rotate_and_flip_mask(act_map, mask_map, rotated_mask_file, joint_beds, zmin, zmax): """ Rotates and flips a NIfTI mask (mask_map) based on the best alignment @@ -1010,81 +924,6 @@ def total_quantification(mask_file, joint_norm_beds, quantification_file, label_ "total_activity_KBq": total_activity_KBq } - -def total_quantification_wholeBody(mask_file, recons_norm_wholeBody_file, quantification_file): #TODO - """ - Compute information of target image per labeled region in reference image and save to a TXT file. - Obtains information about the simulated distribution of activity in the phantom and saves it in a TXT file. - """ - # --- Load images --- - target_img = nib.load(recons_norm_wholeBody_file) - roi_img = nib.load(mask_file) - - - target_data = target_img.get_fdata() - roi_data = roi_img.get_fdata() - - - # --- Ensure same shape --- - if target_data.shape != roi_data.shape: - raise ValueError("Target and reference images must have the same shape.") - - # --- Get all labels except 0 (background) --- - labels = np.unique(roi_data) - labels = labels[labels != 0] - - # --- Calculate mean per label --- - mean_values_KBq = [target_data[roi_data == label].mean() for label in labels] - mean_values_mCi = mean_values_mCi = [v * 0.00002703 for v in mean_values_KBq] #mCi - - #--- Calculate Volumen per label --- - voxel_counts = [np.sum(roi_data == label) for label in labels] - dx, dy, dz = np.array(roi_img.header.get_zooms()[:3]) - voxel_volume = dx * dy * dz - - volumes = [(count * voxel_volume)/1000 for count in voxel_counts] - - # --- Calculate activity per region --- - activity_region_KBq = [m * v for m, v in zip(mean_values_KBq, volumes)] - activity_region_mCi = [a / (3.7 * 10**4) for a in activity_region_KBq] - - # --- Calculate TOTAL activity across all labels --- - total_activity_KBq = float(np.sum(activity_region_KBq)) - total_activity_mCi = float(np.sum(activity_region_mCi)) - - - - # --- Save report as TXT --- - with open(quantification_file, 'w') as f: - # Encabezado - f.write( - f"{'Label':<10}{'MeanValue (mCi/cc)':>20}{'MeanValue (KBq/cc)':>20}{'Pixels':>15}" - f"{'Vol (ccm³)':>15}{'Act (mCi)':>15}{'Act (KBq)':>15}\n" - ) - - # Filas por cada label - for label, mean_val_mCi, mean_val_KBq, count, vol, act_mCi, act_KBq in zip( - labels, mean_values_mCi, mean_values_KBq, voxel_counts, volumes, activity_region_mCi, activity_region_KBq - ): - f.write( - f"{int(label):<10}{mean_val_mCi:>20.4f}{mean_val_KBq:>20.4f}{count:>15.0f}" - f"{vol:>15.4f}{act_mCi:>15.4f}{act_KBq:>15.4f}\n" - ) - # Linea separadora - f.write("="*90 + "\n") - - # Linea TOTAL - f.write( - f"{'TOTAL':<10}{'':>20}{'':>15}{'':>15}" - f"{total_activity_mCi:>15.4f}{total_activity_KBq:>15.4f}\n" - ) - - - return {"labels": labels, "mean_values_mCi": mean_values_mCi, "mean_values_KBq": mean_values_mCi, "pixel_counts": voxel_counts, - "volumes": volumes, "activity_mCi": activity_region_mCi, "activity_KBq": activity_region_KBq, - "total_activity_mCi": total_activity_mCi, "total_activity_KBq": total_activity_KBq - } - def distribution_of_dose_into_phantom(self, maps_dir, act_map, info_act_map, label_file_act_map): """ Computes the actual distribution of activity by region labeled in act_map. @@ -1283,8 +1122,3 @@ def total_fov_correction(self, recons_dir): f.writelines(lines) return (hs_file_new) - - - - - diff --git a/wholebody.py b/wholebody.py index 4a940da8..45617da0 100644 --- a/wholebody.py +++ b/wholebody.py @@ -8,6 +8,10 @@ from utils import tools from utils import wb_tools +def save_cfg(cfg, path): + with open(path, 'w') as f: + yaml.dump(cfg, f, default_flow_style=False) + sys.path.append(str(here())) from simpet import SimPET @@ -68,7 +72,10 @@ def run(self): # THIS WILL POP UP EVEN IF ONLY RECONSTRUCTION IS DONE. MOVE IT OUT if not exists(output_dir): os.makedirs(output_dir) - + + #Summary of General Parameters + save_cfg(self.cfg, join(output_dir, f"{patient_dir}_{self.scanner_model}_wholebody.yaml")) + log_file = join(output_dir, "logfile.log") @@ -113,7 +120,9 @@ def run(self): cfg_copy.params.center_slice = int(cs) cfg_copy.params.simulation_time = sim_time_bed cfg_copy.params.output_dir = join(output_name, f"Bed_{j}_CenterSlice_{cs}") - + + # Save the configuration for each bed + save_cfg(OmegaConf.to_container(cfg_copy), join(bed_dir, f"{patient_dir}_{self.scanner_model}_bed{j}.yaml")) # Run simulation bed_simu = SimPET(cfg_copy) bed_simu.run() @@ -260,43 +269,4 @@ def run(self): print("Quantification finished") except Exception as e: - print(f"Error: Quantification was not performed: {e}") - - - - # #Apply normalization to whole body - # if self.Normalization == 1 and self.joints_beds == 1 and os.path.exists(joint_beds): - # factor_wholeBody = wb_tools.normalization_factor_correction_whole_body(self, joint_beds) - - # # Load reconstructed image - # img = nib.load(joint_beds) - # data = img.get_fdata() - - # # Apply the factor corresponding to each bed - # data_norm = data * factor_wholeBody - - # # Save normalizated image - # recons_dir = os.path.dirname(output_dir) - # recons_norm_wholeBody_file = join(output_dir, 'rec_%s_%s_norm_wholeBody.hdr' % (recons_algorithm, recons_it)) - # img_norm = nib.Nifti1Image(data_norm, img.affine, img.header) - # nib.save(img_norm, recons_norm_wholeBody_file) - - # #Quantification - # quantification_file = join(output_dir, "Quantification_data_whole_body.txt") - - # wb_tools.total_quantification_wholeBody(mask_file, recons_norm_wholeBody_file, quantification_file) - - - # print("Normalization Whole Body completed") - # else: - # print("No Normalization Whole Body completed was performed") - - #wb_tools.coincidences_mask_vs_image(mask_file, act_map) - - - - - - - - + print(f"Error: Quantification was not performed: {e}") \ No newline at end of file