diff --git a/.gitmodules b/.gitmodules index b432fb08..1727b476 100644 --- a/.gitmodules +++ b/.gitmodules @@ -3,4 +3,4 @@ url = https://github.com/UCL/STIR [submodule "submodules/simset_simpet"] path = submodules/simset_simpet - url = https://github.com/txusser/simset_simpet.git + url = https://github.com/arnaufme8/simset_simpet.git diff --git a/README.md b/README.md index c465a7d5..af1c1063 100644 --- a/README.md +++ b/README.md @@ -10,7 +10,7 @@ The SIMPET project is intended to allow to setup and launch MC simulation on a s - Apply the BrainViset procedure to obtain realistic Activity and Attenuation maps. - Run Analytic simulations using STIR simulation procedure and MC simulation using SimSET. -# Installtion +# Installation 1. Install [Git LFS](https://git-lfs.com/). 2. Clone the repository by adding the `--recurse-submodules` flag: @@ -52,6 +52,7 @@ Sometimes, even activating the virtual environemnt, the shell will use the wide ``` make install ``` +_NOTE:_ At the moment, to install SimPET you need to have sudo available in your computer. If you don't (e.g., when installing on a HPC), you can make use of the makefile that does not require sudo: ```makefile-withoutsudo```. To use it, simply remove or change the name of ```makefile```, and rename the file ```makefile-withoutsudo``` to ```makefile```. This will allow the installation without root privileges. 7. Decompress dummy data with: ``` @@ -277,8 +278,9 @@ See [SimSET](https://depts.washington.edu/simset/html/simset_main.html) document #### Attenuation correction -- **analytical_att_correction**: performed by SimSET calcattenuation. +- **analytical_att_correction**: performed by SimSET calcattenuation. (_NOTE: obsolete at the moment and will be removed in further updates_). - **stir_recons_att_corr**: performed in STIR by entering the att image as a normalization map. +- **attenuation_mode**: it has three options. 0 (no attenuation), 1 (attenuation by SimSET calcattenuation) or 2 (attenuation by STIR computation of the ACFs). #### Scatter Correction @@ -316,11 +318,13 @@ See [SimSET](https://depts.washington.edu/simset/html/simset_main.html) document - Jesús Silva-Rodríguez - Pablo Aguiar -- Aida Ninyerola-Baizan +- Aida Niñerola-Baizán - Jeremiah Poveda - Francisco Javier López-González - Nikos Efthimiou -- Arnau Farre +- Arnau Farré-Melero +- Claudia Dominguez-Borrero + # References diff --git a/configs/params/scanner/discovery_st.yaml b/configs/params/scanner/discovery_st.yaml index e7f02dfd..a586a990 100644 --- a/configs/params/scanner/discovery_st.yaml +++ b/configs/params/scanner/discovery_st.yaml @@ -23,7 +23,7 @@ stir_scatt_corr_smoothing: 0 stir_scatt_simulation: 0 analytic_randoms_corr_factor: 0.2 stir_randoms_corr_smoothing: 0 -recons_type: "OSEM3D" +recons_type: "OSEM3D" #"FBP3D" #"OSEM3D" inter_iteration_filter: 0 subiteration_interval: 4 x_dir_filter_FWHM: 1.5 @@ -32,7 +32,9 @@ z_dir_filter_FWHM: 3 psf_value: 0 add_noise: 0 max_segment: 23 -zoomFactor: 1.55 +zoomFactor: 1 #For brain: 1.55 xyOutputSize: 128 zOutputSize: 47 zOutputVoxelSize: 3.27 + +attenuation_mode: 2 #0: no attenuation, 1: attenuation with SimSET, 2: attenuation with STIR diff --git a/configs/params/scanner/siemens_quadra.yaml b/configs/params/scanner/siemens_quadra.yaml new file mode 100644 index 00000000..eac432e1 --- /dev/null +++ b/configs/params/scanner/siemens_quadra.yaml @@ -0,0 +1,42 @@ +scanner_name: "Siemens Quadra" +simset_material: 18 +average_doi: 0.7 +scanner_radius: 41.0 +num_rings: 320 #Real: 320, With module gaps: 323 +axial_fov: 106 +z_crystal_size: 0.32 #cm +transaxial_crystal_size: 0.32 #cm +crystal_thickness: 2 #cm +energy_resolution: 9 +num_aa_bins: 399 +num_td_bins: 520 +min_energy_window: 435 +max_energy_window: 585 +coincidence_window: 4.7 #ns + +numberOfSubsets: 7 +numberOfIterations: 14 +savingInterval: 7 +analytical_att_correction: 0 +stir_recons_att_corr: 0 +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" #Options: "FBP2D", "FBP3D", "OSEM2D", "OSEM3D" +inter_iteration_filter: 0 +subiteration_interval: 4 +x_dir_filter_FWHM: 1.5 +y_dir_filter_FWHM: 1.5 +z_dir_filter_FWHM: 3 +psf_value: 0 +add_noise: 0 + +max_segment: 85 #Reduced protocol: 85; Full protocol: 319 or 322 +zoomFactor: 1 +xyOutputSize: 220 +zOutputSize: 645 +zOutputVoxelSize: 1.65 + +attenuation_mode: 0 #0: no attenuation, 1: attenuation with SimSET, 2: attenuation with STIR diff --git a/configs/params/test.yaml b/configs/params/test.yaml index f47b6f29..42a624a2 100644 --- a/configs/params/test.yaml +++ b/configs/params/test.yaml @@ -11,12 +11,14 @@ patient_dirname: "test_image" act_map: "act.hdr" att_map: "att.hdr" output_dir: "test_image" -center_slice: 7 +center_slice: 0 total_dose: 0.1 simulation_time: 30 sampling_photons: 0 photons: 0 add_randoms: 1 -phglistmode: 0 -detlistmode: 1 -maximumIteration: 1 + +listmode: True # TODO: to unify both listmodes to only the necessary. +phglistmode: 1 #MUST BE ACTIVATED IF listmode. +detlistmode: 0 #At the moment unused, but may be used in the future. +maximumIteration: 1 #Unused at the moment. diff --git a/makefile-withoutsudo b/makefile-withoutsudo new file mode 100644 index 00000000..7a7dba7a --- /dev/null +++ b/makefile-withoutsudo @@ -0,0 +1,231 @@ +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: 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 + +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} ;\ + 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} \ + 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 " - 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/src/simset/simset_sim.py b/src/simset/simset_sim.py index f06ddfcc..67bd5633 100644 --- a/src/simset/simset_sim.py +++ b/src/simset/simset_sim.py @@ -81,6 +81,8 @@ def __init__( self.detlistmode = params.get("detlistmode") self.phglistmode = params.get("phglistmode") self.add_randoms = params.get("add_randoms") + + self.list_mode = params.get("listmode") def run(self): processes = [] @@ -115,7 +117,9 @@ def run(self): div_0_dir = join(self.output_dir, "division_0") for j in ['trues', 'scatter', 'randoms']: - hdr_ = join(div_0_dir, "%s.hdr" % j) + #hdr_ = join(div_0_dir, "%s.hdr" % j) + hdr_ = join(div_0_dir, "%s.nii" % j) + #hdr_ = join(div_0_dir, "%s.nii.gz" % j) if exists(hdr_): counts_ = tools.ncounts(hdr_) print("Number of %s in simulation: %s" % (j, counts_)) @@ -192,6 +196,7 @@ def run_simset_simulation(self, sim_dir): rec_weight = join(sim_dir, "rec.weight") det_hf = join(sim_dir, "det_hf.hist") phg_hf = join(sim_dir, "phg_hf.hist") + if self.s_photons != 0 and self.params.get("add_randoms") != 1: @@ -199,7 +204,8 @@ def run_simset_simulation(self, sim_dir): if self.photons == 0: # Removes counts for preparing for the next simulation - os.remove(rec_weight) + if exists(rec_weight): + os.remove(rec_weight) if exists(det_hf): os.remove(det_hf) if exists(phg_hf): @@ -222,6 +228,7 @@ def run_simset_simulation(self, sim_dir): command = "%s/bin/phg %s > %s" % (self.simset_dir, my_phg, my_log) tools.osrun(command, log_file) + w_quotient = read_ws_from_simset_log(my_log) sim_photons = int(self.s_photons * w_quotient) @@ -231,7 +238,8 @@ def run_simset_simulation(self, sim_dir): sim_photons = self.photons # Removes counts for preparing for the next simulation - os.remove(rec_weight) + if exists(rec_weight): + os.remove(rec_weight) if exists(det_hf): os.remove(det_hf) if exists(phg_hf): @@ -247,10 +255,11 @@ def run_simset_simulation(self, sim_dir): command = "%s/bin/phg %s > %s" % (self.simset_dir, my_phg, my_log) tools.osrun(command, log_file) + if self.add_randoms == 1: coincidence_window = self.scanner.get("coincidence_window") - + simset_tools.add_randoms( sim_dir, self.simset_dir, @@ -258,10 +267,12 @@ def run_simset_simulation(self, sim_dir): rebin=True, log_file=log_file, ) - - simset_tools.process_weights( - rec_weight, sim_dir, self.scanner, self.add_randoms - ) + + if not self.list_mode: + + simset_tools.process_weights( + rec_weight, sim_dir, self.scanner, self.add_randoms + ) print("Finished simulation for %s" % os.path.basename(sim_dir)) @@ -273,6 +284,7 @@ def prepare_simset_files( model_type = self.params.get("model_type") scanner_radius = self.scanner.get("scanner_radius") scanner_axial_fov = self.scanner.get("axial_fov") + list_mode = self.params.get("listmode") # We activate det_listmode if demanded by user or if add_randoms is on if self.add_randoms == 1: @@ -327,6 +339,7 @@ def prepare_simset_files( self.scanner, add_randoms, log_file=log_file, + list_mode=list_mode ) simset_tools.make_index_file(sim_dir, self.simset_dir, log_file=log_file) @@ -336,85 +349,151 @@ def prepare_simset_files( def simulation_postprocessing(self): print("Postprocessing simulation...") print(" ") - # All the parallel simulations are combined in division 0 + + + # All the parallel simulations are combined in division 0. NOTE: CHANGE TO simulation folder at some point. log_file = join(self.output_dir, "postprocessing.log") division_zero = join(self.output_dir, "division_0") + + #ATTENTION: We need to think a better place to do this. At the moment, we keep it here. + #This is generating rec.weight files when listmoding. NOTE: WE WILL USE phg_hf.hist at the moment. Ignoring det.hist. + + if self.list_mode: + + zero_hist = join(division_zero, "phg_hf.hist") + full_hist = join(division_zero, "full_phg_hf.hist") + rec_weight = join(division_zero, "rec.weight") + my_phg = join(division_zero, "phg.rec") + file_list = zero_hist + + print("Adding History Files for PHG") + + message = "Adding History Files for PHG" + tools.log_message(log_file, message) + + for division in range(1, self.divisions): + + division_dir = join(self.output_dir, "division_" + str(division)) + division_hist = join(division_dir, "phg_hf.hist") + file_list = file_list + " " + division_hist + + if self.divisions > 1: + simset_tools.combine_history_files( + self.simset_dir, file_list, full_hist, log_file + ) + + os.remove(zero_hist) + os.rename(full_hist, zero_hist) + + #ATTENTION: THIS IS ADDED ONLY FOR LIST MODE. MAY BECOME UNUSED OR PROBLEMATIC. + + with open(join(division_zero, "bin.rec"), "a") as f: + f.write('STR weight_image_path = "' + rec_weight + '"\n') + f.close() + + command = "%s/bin/bin -p %s" % (self.simset_dir, my_phg) + tools.osrun(command, log_file) + + + + simset_tools.process_weights( + rec_weight, division_zero, self.scanner, self.add_randoms + ) + + + + else: - for image in ["trues", "scatter", "randoms"]: - zero_image = join(division_zero, image + ".hdr") - - if exists(zero_image): - print("Adding sinograms for %s" % image) + for image in ["trues", "scatter", "randoms"]: #TODO: randoms needed here? + zero_image = join(division_zero, image + ".nii") - for division in range(1, self.divisions): - division_dir = join(self.output_dir, "division_" + str(division)) - division_image = join(division_dir, image + ".hdr") - message = "Adding %s from simulation %s" % (image, division) + if exists(zero_image): + print("Adding sinograms for %s" % image) + + #This is to create the postlog file in case there is only one simulation, to avoid double attenuation computing! + message = "Adding sinograms for %s" % image tools.log_message(log_file, message) - tools.operate_images_analyze( - zero_image, division_image, zero_image, "sum" - ) - os.remove(division_image) - os.remove(division_image[0:-3] + "img") - - for hist in ["phg_hf.hist", "det_hf.hist"]: - zero_hist = join(division_zero, hist) - - if exists(zero_hist): - print(" ") - print("Adding History Files for %s" % hist) - - output = join(division_zero, "tmp_" + hist) - - for division in range(1, self.divisions): - division_dir = join(self.output_dir, "division_" + str(division)) - division_hist = join(division_dir, hist) - file_list = zero_hist + " " + division_hist - simset_tools.combine_history_files( - self.simset_dir, file_list, output, log_file - ) - shutil.move(output, zero_hist) - # Once everything is combined in division_0, remove the other division - shutil.rmtree(division_dir) + + for division in range(1, self.divisions): + division_dir = join(self.output_dir, "division_" + str(division)) + division_image = join(division_dir, image + ".nii") + message = "Adding %s from simulation %s" % (image, division) + tools.log_message(log_file, message) + + tools.operate_sinograms_nii( + zero_image, division_image, zero_image, "sum" + ) + + #Remove unused divisions... + for division in range(1, self.divisions): + division_dir = join(self.output_dir, "division_" + str(division)) + shutil.rmtree(division_dir) # Once everything is combined in division_0, remove the other division + - if self.add_randoms == 1: + if self.add_randoms == 1: #TODO: THIS HAS TO BE TESTED IN NEW SETUP. # To have randoms in the final history file, we need to add randoms to the final det_hf.hist print("Adding randoms to the history file...") coincidence_window = self.scanner.get("coincidence_window") - - simset_tools.add_randoms( - division_zero, - self.simset_dir, - coincidence_window, - rebin=False, - log_file=log_file, - ) + + #ATTENTION: AT THE MOMENT RANDOMS + LIST MODE IS NOT WELL IMPLEMENTED, NEED TO HAVE A LOOK: + #NOTE: I THINK THIS PROCESS SHOULD BE DONE BEFORE REMOVING EACH DIVISION. THAT IS: EXTRACTING AND COMBINING FILES FOR EACH DIVISION. CHECK IF WORKS! + #WARNING: AT THE MOMENT THIS IS NOT NEEDED. RECOVER WHEN DETLISTMODE IS USED. + #simset_tools.add_randoms( + # division_zero, + # self.simset_dir, + # coincidence_window, + # rebin=False, + # log_file=log_file, + #) # os.remove(join(division_zero, "sorted_det_hf.hist")) - - det_hist = join(division_zero, "det_hf.hist") - randoms_hist = join(division_zero, "randoms.hist") - output = join(division_zero, "full_det_hf.hist") - - file_list = det_hist + " " + randoms_hist - - simset_tools.combine_history_files( - self.simset_dir, file_list, output, log_file + + #WARNING: AT THE MOMENT THIS IS NOT NEEDED. RECOVER WHEN DETLISTMODE IS USED. + #det_hist = join(division_zero, "det_hf.hist") + #randoms_hist = join(division_zero, "randoms.hist") + #output = join(division_zero, "full_det_hf.hist") + + #file_list = det_hist + " " + randoms_hist + + #simset_tools.combine_history_files( + # self.simset_dir, file_list, output, log_file + #) + + self.attenuation_mode = self.scanner.get("attenuation_mode") + + if self.attenuation_mode == 1: #self.stir_norm_from_att_map != 1: + + print("Calculating attenuation map...") + print(" ") + + output_atten = "attenuationsino" + hdr_to_copy = join("trues.nii") + + simset_tools.simset_calcattenuation( + self.simset_dir, division_zero, output_atten, hdr_to_copy, nrays=1, timeout=None ) - - print("Calculating attenuation map...") - print(" ") - - output_atten = "attenuationsino" - hdr_to_copy = join("trues.hdr") - - simset_tools.simset_calcattenuation( - self.simset_dir, division_zero, output_atten, hdr_to_copy, nrays=1 - ) - + + #TODO: Add a parameter that allows to remove unnecessary files, for saving space... + """ + if exists(join(division_zero, "rec.weight")): + os.remove(join(division_zero, "rec.weight")) + + if exists(join(division_zero, "rec.act_indexes")): + os.remove(join(division_zero, "rec.act_indexes")) + + if exists(join(division_zero, "rec.att_indexes")): + os.remove(join(division_zero, "rec.att_indexes")) + + if exists(join(division_zero, "rec.activity_image")): + os.remove(join(division_zero, "rec.activity_image")) + + if exists(join(division_zero, "rec.attenuation_image")): + os.remove(join(division_zero, "rec.attenuation_image")) + """ + class SimSET_Reconstruction(object): """This class provides functions to reconstruct a SimSET simulation.""" @@ -426,7 +505,7 @@ def __init__( self.simpet_dir = dirname(abspath(__file__)) self.dir_stir = config.get("dir_stir") - self.input_dir = join(projections_dir, "division_0") + self.input_dir = join(projections_dir, "division_0") #CHANGE TO "SIMULATION"!!! self.output_dir = reconstructions_dir self.params = params @@ -451,7 +530,8 @@ def __init__( def run(self): if not exists(self.output_dir): os.makedirs(self.output_dir) - + + self.prepare_recons() self.run_recons() @@ -459,71 +539,206 @@ def prepare_recons(self): from src.stir import stir_tools print("Preparing files for reconstruction") - - trues_sino = join(self.input_dir, "trues.hdr") - scatter_sino = join(self.input_dir, "scatter.hdr") - randoms_sino = join(self.input_dir, "randoms.hdr") - - corr_scatter_sino = join(self.input_dir, "corr_scatter.hdr") - corr_randoms_sino = join(self.input_dir, "corr_randoms.hdr") - my_simset_sino = join(self.input_dir, "my_sinogram.hdr") - additive_sinogram = join(self.input_dir, "additive_sinogram.hdr") - - tools.operate_single_image( - scatter_sino, - "mult", - self.scatt_corr_factor, - corr_scatter_sino, - self.log_file, - ) - tools.operate_images_analyze( - trues_sino, corr_scatter_sino, my_simset_sino, operation="sum" - ) + + self.simset_dir = self.config.get("dir_simset") + self.stir_dir = self.config.get("dir_stir") + + num_rings = self.scanner.get("num_rings") + max_segment = self.scanner.get("max_segment") + + self.attenuation_mode = self.scanner.get("attenuation_mode") + + if ((self.attenuation_mode == 1) and (not exists(join(self.input_dir, "attenuationsino.nii")))): + + print("Attenuation map was not computed: Calculating attenuation map...") + print(" ") + + output_atten = "attenuationsino" + hdr_to_copy = join("trues.nii") + + simset_tools.simset_calcattenuation( + self.simset_dir, self.input_dir, output_atten, hdr_to_copy, nrays=1, timeout=None + ) + + trues_sino = join(self.input_dir, "trues.nii") + scatter_sino = join(self.input_dir, "scatter.nii") + randoms_sino = join(self.input_dir, "randoms.nii") + + corr_scatter_sino = join(self.input_dir, "corr_scatter.nii") + corr_randoms_sino = join(self.input_dir, "corr_randoms.nii") + my_simset_sino = join(self.input_dir, "my_sinogram.nii") + additive_sinogram = join(self.input_dir, "additive_sinogram.nii") + att_sino = join(self.input_dir, "attenuationsino.nii") + + print("MULTIPLYING CORR_SCATTER_SINO") + + if not exists(corr_scatter_sino): + tools.operate_single_image_nii( + scatter_sino, + "mult", + self.scatt_corr_factor, + corr_scatter_sino, + self.log_file, + check_nans = False, + ) + + print("SUMMING TRUES AND CORR SCATTER") + + if not exists(my_simset_sino): + tools.operate_images_nii( + trues_sino, corr_scatter_sino, my_simset_sino, operation="sum", check_nans = False, + ) if self.add_randoms == 1: - tools.operate_single_image( + tools.operate_single_image_nii( randoms_sino, "mult", self.random_corr_factor, corr_randoms_sino, self.log_file, + #check_nans = False ) - tools.operate_images_analyze( - my_simset_sino, corr_randoms_sino, my_simset_sino, operation="sum" + + tools.operate_sinograms_nii( + my_simset_sino, corr_randoms_sino, my_simset_sino, operation="sum", #check_nans = False, ) if self.scanner.get("stir_randoms_corr_smoothing") == 1: - tools.operate_images_analyze( - scatter_sino, randoms_sino, additive_sinogram, operation="sum" + + tools.operate_sinograms_nii( + scatter_sino, randoms_sino, additive_sinogram, operation="sum", #check_nans = False, ) else: - tools.copy_analyze(scatter_sino, additive_sinogram) + tools.copy_nifti(scatter_sino, additive_sinogram) else: - tools.copy_analyze(scatter_sino, additive_sinogram) - - tools.smooth_analyze(additive_sinogram, 10, additive_sinogram) - - sinogram_stir = join(self.output_dir, "stir_sinogram.hdr") - tools.convert_simset_sino_to_stir(my_simset_sino, sinogram_stir) - shutil.copy(sinogram_stir[0:-3] + "img", sinogram_stir[0:-3] + "s") - stir_tools.create_stir_hs_from_detparams( - self.scanner, sinogram_stir[0:-3] + "hs" - ) - - additive_sino_stir = join(self.output_dir, "stir_additivesino.hdr") - tools.convert_simset_sino_to_stir(additive_sinogram, additive_sino_stir) - shutil.copy(additive_sino_stir[0:-3] + "img", additive_sino_stir[0:-3] + "s") + if self.scanner.get("stir_scatt_corr_smoothing") == 1: + print("COPYING ADDITIVE SINO") + if not exists(additive_sinogram): + tools.copy_nifti(scatter_sino, additive_sinogram) + print("SMOOTHING ADDITIVE SINO") + tools.smooth_nifti(additive_sinogram, 10, additive_sinogram) + + + print("GENERATING STIR SINOGRAM") + + sinogram_stir_nii = join(self.input_dir, "stir_sinogram.nii") + sinogram_stir_s = join(self.output_dir, "stir_sinogram.s") + sinogram_stir_hs = join(self.output_dir, "stir_sinogram.hs") + + if not exists(sinogram_stir_nii): + tools.convert_simset_sino_to_stir_nii(my_simset_sino, sinogram_stir_nii) + + tools.copy_sinogram_stir_to_output(sinogram_stir_nii, sinogram_stir_s) + #tools.copy_reduced_sinogram_stir_to_output(sinogram_stir_nii, sinogram_stir_s, num_rings, max_segment) #NOTE: Needed for FBP + stir_tools.create_stir_hs_from_detparams( - self.scanner, additive_sino_stir[0:-3] + "hs" + self.scanner, sinogram_stir_hs ) + + + if self.scanner.get("stir_scatt_corr_smoothing") == 1: + + print("GENERATING STIR ADDITIVE") #TODO: CHECK IF NECESSARY IN FUNCTION OF stir_scatt_corr_smoothing: + + additivesino_stir_nii = join(self.input_dir, "stir_additivesino.nii") + additivesino_stir_s = join(self.output_dir, "stir_additivesino.s") + additivesino_stir_hs = join(self.output_dir, "stir_additivesino.hs") + + if not exists(additivesino_stir_nii): + tools.convert_simset_sino_to_stir_nii(additive_sinogram, additivesino_stir_nii) + + tools.copy_sinogram_stir_to_output(additivesino_stir_nii, additivesino_stir_s) + #tools.copy_reduced_sinogram_stir_to_output(additivesino_stir_nii, additivesino_stir_s, num_rings, max_segment) #NOTE: Needed for FBP + + stir_tools.create_stir_hs_from_detparams( + self.scanner, additivesino_stir_hs + ) + + print("GENERATING STIR ATT") + att_stir_nii = join(self.input_dir, "stir_att.nii") + att_stir_s = join(self.output_dir, "stir_att.s") + att_stir_hs = join(self.output_dir, "stir_att.hs") + + + if self.attenuation_mode == 1: + if not exists(att_stir_nii): + tools.convert_simset_sino_to_stir_nii(att_sino, att_stir_nii) + + tools.copy_sinogram_stir_to_output(att_stir_nii, att_stir_s) + #tools.copy_reduced_sinogram_stir_to_output(att_stir_nii, att_stir_s, num_rings, max_segment) #NOTE: Needed for FBP + + stir_tools.create_stir_hs_from_detparams( + self.scanner, att_stir_hs + ) + + + #Obtain Mu-map directly from att map: + + if self.attenuation_mode == 2: + + if not exists(att_stir_nii): #Added to save time... Remove if this gives problems. + + nib.save(nib.load(join(self.input_dir, 'stir_sinogram.nii')), join(self.input_dir, 'stir_sinogram.hdr')) #DELETE, this is only for testing... + + attmap = nib.load(join('Data', self.params.get("patient_dirname"), self.params.get("att_map"))) + att_indexes = np.uint8(attmap.get_fdata()) + + attmap_dim = np.shape(attmap) + attmap_pixsize = attmap.header['pixdim'][1:4] + + mu_map = np.zeros(np.shape(att_indexes), dtype=np.float32) + + for attidx in np.unique(att_indexes): + mu_map[att_indexes == attidx] = tools.mu_coef_511keV(attidx) + + mu_map_img = nib.Nifti1Image(mu_map, attmap.affine) + nib.save(mu_map_img, join(self.input_dir, 'mu_map.hdr')) + nib.save(mu_map_img, join(self.input_dir, 'mu_map.nii')) + + shutil.copy(join(self.input_dir, 'mu_map.img'), join(self.output_dir, 'mu_map.v')) + tools.write_interfile_header_mu(join(self.output_dir, 'mu_map.hv'), attmap_dim[0], attmap_pixsize[0], + attmap_dim[1], attmap_pixsize[1], + attmap_dim[2], attmap_pixsize[2]) + + + tools.write_fwdproj_parfile(join(self.output_dir, "fwdproj_par.par")) + + + #Command to compute ACF + #NOTE: For total-body this is crashing at some point of the execution... Wait for further versions to be solved. + ACF_command = "%s --ACF %s %s %s %s 2>&1 | tee %s" % (join(self.stir_dir, "bin", "calculate_attenuation_coefficients"), join(self.output_dir, "stir_att.hs"), join(self.output_dir, 'mu_map.hv'), sinogram_stir_hs, join(self.output_dir, "fwdproj_par.par"), join(self.output_dir, "att_logging.log")) + + #TODO: Still have to think about this... Way to save the nifti to the simulation folder. + os.system(ACF_command) + + shutil.copyfile(join(self.output_dir, "stir_att.s"), join(self.input_dir, "stir_att.img")) + + os.rename(join(self.input_dir, 'stir_sinogram.hdr'), join(self.input_dir, 'stir_att.hdr')) + + stiratt_img = nib.load(join(self.input_dir, "stir_att.img")) + stiratt_data = stiratt_img.dataobj + + stiratt_nifti2 = nib.Nifti2Image(stiratt_data, stiratt_img.affine, stiratt_img.header) + nib.save(stiratt_nifti2, join(self.input_dir, "stir_att.nii")) + + #Remove not needed files anymore... + os.remove(join(self.input_dir, 'mu_map.hdr')) + os.remove(join(self.input_dir, 'mu_map.img')) + os.remove(join(self.input_dir, 'stir_att.hdr')) + os.remove(join(self.input_dir, 'stir_att.img')) + os.remove(join(self.input_dir, 'stir_sinogram.img')) + + else: + + tools.copy_sinogram_stir_to_output(att_stir_nii, att_stir_s) - att_sino = join(self.input_dir, "attenuationsino.hdr") - att_stir = join(self.output_dir, "stir_att.hdr") - tools.convert_simset_sino_to_stir(att_sino, att_stir) - shutil.copy(att_stir[0:-3] + "img", att_stir[0:-3] + "s") - stir_tools.create_stir_hs_from_detparams(self.scanner, att_stir[0:-3] + "hs") - + stir_tools.create_stir_hs_from_detparams( + self.scanner, att_stir_hs, output_format = "STIR" + ) + + + #WARNING: This has not been revised yet for the new version. if self.scanner.get("analytical_att_correction") == 1: catt_sino = join(self.output_dir, "catt_sinogram.hdr") tools.operate_images_analyze( @@ -538,7 +753,9 @@ def prepare_recons(self): shutil.copy(catt_add_sino[0:-3] + "img", additive_sino_stir[0:-3] + "s") if self.scanner.get("psf_value") != 0: - stir_tools.apply_psf(self.scanner, sinogram_stir, self.log_file) + #stir_tools.apply_psf(self.scanner, sinogram_stir, self.log_file) #Original, recover if not working. + stir_tools.apply_psf(self.scanner, join(self.output_dir, "stir_sinogram.nii"), self.log_file) + if self.scanner.get("add_noise") != 0: stir_tools.add_noise( @@ -547,20 +764,38 @@ def prepare_recons(self): def run_recons(self): from src.stir import stir_tools - + + start_recons = False + print("Starting STIR reconstruction") recons_algorithm = self.scanner.get("recons_type") 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] - ): - print("Something is not ready for the reconstruction") + + + if self.scanner.get("attenuation_mode") != 0: + + if any( + exists(i) == False for i in [sinogram_stir, att_stir] #[sinogram_stir, additive_sino_stir, att_stir] #NOTE: Integrate in the future... + ): + print("Something is not ready for the reconstruction") + else: + start_recons = True + print("Starting STIR reconstruction") + else: - print("Starting STIR reconstruction") + + if any( + exists(i) == False for i in [sinogram_stir] #[sinogram_stir, additive_sino_stir] #NOTE: Integrate in the future... + ): + print("Something is not ready for the reconstruction") + else: + start_recons = True + print("Starting STIR reconstruction") + + if start_recons: if recons_algorithm == "FBP2D": reconsFile_hdr = stir_tools.FBP2D_recons( diff --git a/src/simset/simset_tools.py b/src/simset/simset_tools.py index a92b6b7f..54f5c588 100644 --- a/src/simset/simset_tools.py +++ b/src/simset/simset_tools.py @@ -1,655 +1,684 @@ -import random -import os -import shutil -import subprocess as sp -from os.path import join, exists -import pexpect -from utils import tools - - -def make_simset_act_table(act_table_factor, my_act_table, log_file=False): - with open(my_act_table, "w") as f: - f.write("256\n") - for i in range(256): - act_value = i * 0.000001 * act_table_factor - f.write(str("{:.12f}".format(act_value)) + "\n") - f.close() - - if log_file: - message = "Wrote act_table with factor %s" % act_table_factor - tools.log_message(log_file, message, "info") - - -def make_simset_phg( - config, - output_file, - simulation_dir, - act, - scanner_radius, - scanner_axial_fov, - center_slice, - photons, - sim_time, - add_randoms=False, - phg_hf=False, - S=0, - log_file=False, -): - stir_identifier = "# Hello, I am a SimSET PHG file!\n" - - # define types to be used - boolean = "BOOL " - real = "REAL " - longl = "LONGLONG " - integer = "INT " - enum = "ENUM " - string = "STR " - - # Configuring the importance sampling parameters - if add_randoms == 1: - stratification = "false" - forced_detection = "false" - non_absortion = "false" - # We need to force coincidences + singles for randoms - sim_PET_coinc_only = "false" - sim_PET_coinc_plus_singles = "true" - else: - stratification = config.get("stratification") - forced_detection = config.get("forced_detection") - non_absortion = config.get("forced_non_absortion") - sim_PET_coinc_only = "false" - sim_PET_coinc_plus_singles = "true" - - # Configuration of photons an simulation time - - num_to_simulate = str(photons) - length_of_scan = str(sim_time) - - # Import other configurations form config.yml - acceptance_angle = str(config.get("acceptance_angle")) - positron_range = config.get("positron_range") - non_col = config.get("non_colinearity") - minimum_energy = str(config.get("minimum_energy")) - weight_window_ratio = str(config.get("weight_window_ratio")) - point_source_voxels = config.get("point_source_voxels") - coherent_scatter_object = config.get("coherent_scatter_object") - coherent_scatter_detector = config.get("coherent_scatter_detector") - simulated_isotope = config.get("isotope") - - # Generate a random seed for the MC - random_seed = str(random.randint(1, 1e12)) - - # Object stuff - nslices = act.shape[2] - 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 - ) # 20 is because it must be cm - yMin, yMax = round(-act_fov[1, 1] / 20, 3), round(act_fov[1, 1] / 20, 3) - - z_offset = abs( - act.affine[2, 2] * center_slice / 10 + 0.5 * act.affine[2, 2] / 10 - ) # cm - zMin, zMax = round(-z_offset, 3), round(act_fov[2, 2] / 10 - z_offset, 3) - - dz = round((zMax - zMin) / nslices, 2) - - max_z_target = scanner_axial_fov / 2 - min_z_target = -scanner_axial_fov / 2 - - # Directory variables - simset_dir = config.get("dir_simset") - simset_phgdata_dir = join(simset_dir, "phg.data") - - # Now we write the phg file - with open(output_file, "w") as f: - # First we write the configuration stuff for the simulation - f.write(stir_identifier) - f.write("\n\n# Runtime options\n") - f.write(boolean + "simulate_stratification = %s \n" % stratification) - f.write(boolean + "simulate_forced_detection = " + forced_detection + "\n") - f.write(boolean + "forced_non_absorbtion = " + non_absortion + "\n") - f.write(real + "acceptance_angle = " + acceptance_angle + "\n") - f.write(longl + "num_to_simulate = " + num_to_simulate + "\n") - f.write(real + "length_of_scan = " + length_of_scan + "\n") - f.write(boolean + "simulate_SPECT = false\n") - f.write( - boolean + "simulate_PET_coincidences_only = " + sim_PET_coinc_only + "\n" - ) - f.write( - boolean - + "simulate_PET_coincidences_plus_singles = " - + sim_PET_coinc_plus_singles - + "\n" - ) - f.write(boolean + "adjust_for_positron_range = " + positron_range + "\n") - f.write(boolean + "adjust_for_collinearity = " + non_col + "\n") - f.write(real + "minimum_energy = " + minimum_energy + "\n") - f.write(real + "photon_energy = 511.0\n") - f.write(real + "weight_window_ratio = " + weight_window_ratio + "\n") - f.write(boolean + "point_source_voxels = " + point_source_voxels + "\n") - f.write(integer + "random_seed = " + random_seed + "\n") - f.write( - boolean - + "model_coherent_scatter_in_obj = " - + coherent_scatter_object - + "\n" - ) - f.write( - boolean - + "model_coherent_scatter_in_tomo = " - + coherent_scatter_detector - + "\n" - ) - f.write(enum + "isotope = " + simulated_isotope + "\n") - - # Now comes the object stuff (old make_phg_simset from fruitcake) - f.write("\n\n# OBJECT GEOMETRY VALUES\n") - f.write("\nNUM_ELEMENTS_IN_LIST object = %s" % str(nslices + 1)) - 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) - f.write( - "\n NUM_ELEMENTS_IN_LIST slice = 9 " - + "\n INT slice_number = %s" % str(i) - + "\n REAL zMin = %s" % str(zMin_value) - + "\n REAL zMax = %s" % str(zMax_value) - + "\n REAL xMin = %s" % str(xMin) - + "\n REAL xMax = %s" % str(xMax) - + "\n REAL yMin = %s" % str(yMin) - + "\n REAL yMax = %s" % str(yMax) - + "\n INT num_X_bins = %s" % str(xbins) - + "\n INT num_Y_bins = %s" % str(ybins) - ) - - # Now we fill the values for the target cylinder - - f.write("\n\n# TARGET CYLINDER INFORMATION") - f.write("\nNUM_ELEMENTS_IN_LIST target_cylinder = 3") - f.write("\n REAL target_zMin = %s" % str(min_z_target)) - f.write("\n REAL target_zMax = %s" % str(max_z_target)) - f.write("\n REAL radius = %s\n\n" % str(scanner_radius)) - - # Now we need the directory stuff - f.write( - string - + 'coherent_scatter_table = "' - + join(simset_phgdata_dir, "coh.tables") - + '"\n' - ) - f.write( - string - + 'activity_indexes = "' - + join(simulation_dir, "rec.act_indexes") - + '"\n' - ) - f.write( - string - + 'activity_table = "' - + join(simulation_dir, "phg_act_table") - + '"\n' - ) - f.write( - string - + 'activity_index_trans = "' - + join(simset_phgdata_dir, "phg_act_index_trans") - + '"\n' - ) - f.write( - string - + 'activity_image = "' - + join(simulation_dir, "rec.activity_image") - + '"\n' - ) - f.write( - string - + 'attenuation_indexes = "' - + join(simulation_dir, "rec.att_indexes") - + '"\n' - ) - f.write( - string - + 'attenuation_table = "' - + join(simset_phgdata_dir, "phg_att_table") - + '"\n' - ) - f.write( - string - + 'attenuation_index_trans = "' - + join(simset_phgdata_dir, "phg_att_index_trans") - + '"\n' - ) - f.write( - string - + 'attenuation_image = "' - + join(simulation_dir, "rec.attenuation_image") - + '"\n' - ) - - # This part changes productivity tables between S=0 and S=1 adq - if S == 1: - f.write( - string - + 'productivity_input_table = "' - + join(simulation_dir, "sampling_rec") - + '"\n' - ) - f.write(string + 'productivity_output_table = ""\n') - else: - f.write(string + 'productivity_input_table = ""\n') - f.write( - string - + 'productivity_output_table = "' - + join(simulation_dir, "sampling_rec") - + '"\n' - ) - - f.write( - string + 'statistics_file = "' + join(simulation_dir, "rec.stat") + '"\n' - ) - f.write( - string - + 'isotope_data_file = "' - + join(simset_phgdata_dir, "isotope_positron_energy_data") - + '"\n' - ) - f.write( - string - + 'detector_params_file = "' - + join(simulation_dir, "det.rec") - + '"\n' - ) - f.write( - string + 'bin_params_file = "' + join(simulation_dir, "bin.rec") + '"\n' - ) - # Adds the phg history file if it is activated - if phg_hf == 1: - f.write( - string - + 'history_file = "' - + join(simulation_dir, "phg_hf.hist" + '"\n') - ) - - f.close() - - if log_file: - message = "Created phg_file in: %s" % output_file - tools.log_message(log_file, message, "info") - - -def make_simset_bin( - config, output_file, simulation_dir, scanner, add_randoms=False, log_file=False -): - stir_identifier = "# Hello, I am a SimSET BIN file!\n" - - # define types to be used - boolean = "BOOL " - real = "REAL " - integer = "INT " - string = "STR " - - # We will create an extra bin for randoms if the randoms simulation is activated. - if add_randoms: - accept_randoms = "true" - scatter_param = "6" - else: - accept_randoms = "false" - scatter_param = "1" - - min_s = "0" - max_s = "9" - - # Here we get the parameters from the parameter files. - num_z_bins = str(scanner.get("num_rings")) - axial_fov = scanner.get("axial_fov") - min_z, max_z = -axial_fov / 2, axial_fov / 2 - num_aa_bins = str(scanner.get("num_aa_bins")) - num_td_bins = str(scanner.get("num_td_bins")) - min_td = str(-scanner.get("scanner_radius")) - max_td = str(scanner.get("scanner_radius")) - min_e = str(scanner.get("min_energy_window")) - max_e = str(scanner.get("max_energy_window")) - rec_weight_file = join(simulation_dir, "rec.weight") - - with open(output_file, "w") as f: - f.write(stir_identifier) - f.write(boolean + " accept_randoms = " + accept_randoms + "\n") - f.write(integer + "scatter_param = " + scatter_param + "\n") - f.write(integer + "min_s = " + min_s + "\n") - f.write(integer + "max_s = " + max_s + "\n") - f.write(integer + "num_z_bins = " + num_z_bins + "\n") - f.write(real + "min_z = " + str(min_z) + "\n") - f.write(real + "max_z = " + str(max_z) + "\n") - f.write(integer + "num_aa_bins = " + num_aa_bins + "\n") - f.write(integer + "num_td_bins = " + num_td_bins + "\n") - f.write(real + "min_td = " + min_td + "\n") - f.write(real + "max_td = " + max_td + "\n") - f.write("INT num_e1_bins = 1\n") - f.write("INT num_e2_bins = 1\n") - f.write(real + "min_e = " + min_e + "\n") - f.write(real + "max_e = " + max_e + "\n") - f.write("INT weight_image_type = 2\n") - f.write("INT count_image_type = 2\n") - f.write("BOOL add_to_existing_img = false\n") - f.write(string + 'weight_image_path = "' + rec_weight_file + '"\n') - - f.close() - - if log_file: - message = "Created bin_file in: %s" % output_file - tools.log_message(log_file, message, "info") - - -def make_simset_simp_det(scanner_params, output, sim_dir, det_hf=0, log_file=False): - energy_resolution = scanner_params.get("energy_resolution") - new_file = open(output, "w") - new_file.write( - "ENUM detector_type = simple_pet \n\n" - + "REAL reference_energy_keV = 511.0 \n" - + "REAL energy_resolution_percentage = %s \n" % energy_resolution - ) - - if det_hf == 1: - new_file.write( - 'STR history_file = "' + join(sim_dir, "det_hf.hist" + '"\n') - ) - - new_file.close() - - if log_file: - message = ( - "Created det_file with:\n" + "Energy resolution: %s" % energy_resolution - ) - - tools.log_message(log_file, message, "info") - - -def make_simset_cyl_det(scanner_params, output, sim_dir, det_hf=0, log_file=False): - num_rings = scanner_params.get("num_rings") - z_crystal_size = scanner_params.get("z_crystal_size") - axial_fov = scanner_params.get("axial_fov") - max_z = axial_fov / 2 - min_z = -axial_fov / 2 - gap_z_size = (max_z - min_z - z_crystal_size * num_rings) / (num_rings - 1) - cyln_inner_radius = scanner_params.get("scanner_radius") - cyln_outer_radius = cyln_inner_radius + scanner_params.get("crystal_thickness") - energy_resolution = scanner_params.get("energy_resolution") - timing_resolution = scanner_params.get("timing_resolution") - material = scanner_params.get("simset_material") - - nrings_total = 2 * num_rings - 1 - - new_file = open(output, "w") - new_file.write( - "ENUM detector_type = cylindrical \n\n" - + "# This detector example has %s axial rings and %s gaps \n" - % (num_rings, num_rings - 1) - + "INT cyln_num_rings = %s \n\n" % nrings_total - ) - - for i in range(1, num_rings + 1): - ring_zmin = min_z + (i - 1) * z_crystal_size + (i - 1) * gap_z_size - ring_zmax = ring_zmin + z_crystal_size - gap_zmax = ring_zmin + z_crystal_size + gap_z_size - - new_file.write( - "# RING #%s \n" % i - + "# The following defines the ring parameters \n" - + "LIST cyln_ring_info_list = 5 \n" - + "INT cyln_num_layers = 1 \n" - + "LIST cyln_layer_info_list = 4 \n" - + "BOOL cyln_layer_is_active = TRUE \n" - + "INT cyln_layer_material = %s \n" % material - + "REAL cyln_layer_inner_radius = %s \n" % cyln_inner_radius - + "REAL cyln_layer_outer_radius = %s \n" % cyln_outer_radius - + "REAL cyln_min_z = %s \n" % ring_zmin - + "REAL cyln_max_z = %s \n\n" % ring_zmax - ) - - if not i == num_rings: - new_file.write( - "# GAP #%s \n" % i - + "# The following defines the gap parameters \n" - + "LIST cyln_ring_info_list = 5 \n" - + "INT cyln_num_layers = 1 \n" - + "LIST cyln_layer_info_list = 4 \n" - + "BOOL cyln_layer_is_active = FALSE \n" - + "INT cyln_layer_material = 0 \n" - + "REAL cyln_layer_inner_radius = %s \n" % cyln_inner_radius - + "REAL cyln_layer_outer_radius = %s \n" % cyln_outer_radius - + "REAL cyln_min_z = %s \n" % ring_zmax - + "REAL cyln_max_z = %s \n\n" % gap_zmax - ) - - new_file.write( - "REAL reference_energy_keV = 511.0 \n" - + "REAL energy_resolution_percentage = %s \n" % energy_resolution - + "REAL photon_time_fwhm_ns = %s \n" % timing_resolution - ) - if det_hf == 1: - new_file.write( - 'STR history_file = "' + join(sim_dir, "det_hf.hist" + '"\n') - ) - - new_file.close() - - if log_file: - message = ( - "Created det_file with:\n" - + "Cristal z: %s cm\n" % z_crystal_size - + "Gap size: %s cm\n" % gap_z_size - + "Ring thickness: %s cm" % (cyln_outer_radius - cyln_inner_radius) - + "Energy resolution: %s" % energy_resolution - + "Timing resolution: %s" % timing_resolution - ) - - tools.log_message(log_file, message, "info") - - -def make_index_file(simulation_dir, simset_dir, log_file=False): - output = join(simulation_dir, "index_file.log") - phg_file = join(simulation_dir, "phg.rec") - act_dat = join(simulation_dir, "act.dat") - att_dat = join(simulation_dir, "att.dat") - make_index_file = join(simset_dir, "bin", "makeindexfile") - - prompts = ( - phg_file - + "\n" - + "y\n" - + "y\n" - + "0\n" - + "y\n" - + act_dat - + "\n" - + "0\n" - + "0\n" - + "1\n" - + "n\n" - + "n\n" - + "y\n" - + "y\n" - + "0\n" - + "y\n" - + att_dat - + "\n" - + "0\n" - + "0\n" - + "1\n" - + "n\n" - + "n\n" - ) - - os.system('printf "%s" | %s > %s' % (prompts, make_index_file, output)) - - if log_file: - message = 'printf "%s" | %s > %s' % (prompts, make_index_file, output) - tools.log_message(log_file, message, "info") - - -def process_weights(weights_file, output_dir, scanner, add_randoms=0): - nbins = scanner.get("num_td_bins") - nangles = scanner.get("num_aa_bins") - nrings = scanner.get("num_rings") - nslices = nrings * nrings - - Simset_offset = 32768 - block_size = nbins * nangles * nslices * 4 - - trues_start = Simset_offset - trues_end = Simset_offset + block_size - trues_file = join(output_dir, "w1") - - with open(weights_file, "rb") as in_file: - with open(trues_file, "wb") as out_file: - out_file.write(in_file.read()[trues_start:trues_end]) - - output = join(output_dir, "trues.hdr") - tools.create_analyze_from_imgdata( - trues_file, output, nbins, nangles, nslices, 1, 1, 1, "fl" - ) - os.remove(trues_file) - - scatter_start = trues_end - scatter_end = trues_end + block_size - scatter_file = join(output_dir, "w2") - - with open(weights_file, "rb") as in_file: - with open(scatter_file, "wb") as out_file: - out_file.write(in_file.read()[scatter_start:scatter_end]) - - output = join(output_dir, "scatter.hdr") - tools.create_analyze_from_imgdata( - scatter_file, output, nbins, nangles, nslices, 1, 1, 1, "fl" - ) - os.remove(scatter_file) - - if add_randoms == 1: - randoms_start = scatter_end - randoms_end = scatter_end + block_size - randoms_file = join(output_dir, "w3") - - with open(weights_file, "rb") as in_file: - with open(randoms_file, "wb") as out_file: - out_file.write(in_file.read()[randoms_start:randoms_end]) - - output = join(output_dir, "randoms.hdr") - tools.create_analyze_from_imgdata( - randoms_file, output, nbins, nangles, nslices, 1, 1, 1, "fl" - ) - os.remove(randoms_file) - - -def add_randoms(sim_dir, simset_dir, coincidence_window, rebin=True, log_file=False): - string = "STR " - integer = "INT " - enum = "ENUM " - real = "REAL " - - sorted_file = join(sim_dir, "sorted_det_hf.hist") - template = join(sim_dir, "sort.params") - det_hf = join(sim_dir, "det_hf.hist") - - with open(template, "w") as f: - f.write(string + 'history_file = "' + det_hf + '"\n') - f.write(string + 'sorted_history_file = "' + sorted_file + '"\n') - f.write(integer + "buffer_size = 1024\n") - - f.close() - - sorting_bin = join(simset_dir, "bin", "timesort -d") - - command = command = "%s %s >> %s" % (sorting_bin, template, log_file) - tools.osrun(command, log_file) - - randoms_file = join(sim_dir, "randoms.hist") - template = join(sim_dir, "randoms.params") - - with open(template, "w") as f: - f.write(enum + "detector_type = cylindrical\n") - f.write(string + 'history_file = "' + sorted_file + '"\n') - f.write(string + 'randoms_history_file = "' + randoms_file + '"\n') - f.write( - real + "coincidence_timing_window_in_ns = " + str(coincidence_window) + "\n" - ) - - f.close() - - addrand_bin = join(simset_dir, "bin", "addrandoms") - - command = "%s %s >> %s" % (addrand_bin, template, log_file) - tools.osrun(command, log_file) - - if exists(sorted_file): - os.remove(sorted_file) - - # The following will replace the existing det-hist file with the new including randoms - - if rebin: - det_file = join(sim_dir, "det.rec") - - with open(det_file, "r") as f: - lines = f.readlines() - with open(det_file, "w") as f: - for line in lines: - line = line.replace("det_hf.hist", "randoms.hist") - f.write(line) - - binfile = join(sim_dir, "phg.rec") - phgbin = join(simset_dir, "bin", "bin -d") - - command = command = "%s %s >> %s" % (phgbin, binfile, log_file) - tools.osrun(command, log_file) - - -def combine_history_files(simset_dir, history_files, output, log_file): - combinehist = join(simset_dir, "bin", "combinehist") - - rcommand = "%s %s %s" % (combinehist, history_files, output) - - proc = sp.Popen( - rcommand, - universal_newlines=True, - shell=True, - stdin=sp.PIPE, - stdout=sp.PIPE, - stderr=sp.PIPE, - ).communicate("Yes\n") - - # tools.osrun(rcommand, log_file) - - -def simset_calcattenuation( - simset_dir, sim_dir, output, hdr_to_copy, nrays=1, timeout=36000 -): - calcattenuation = join(simset_dir, "bin", "calcattenuation") - - current_dir = os.getcwd() - os.chdir(sim_dir) - - child = pexpect.spawn(calcattenuation, timeout=timeout) - # child.logfile = sys.stdout.buffer #I comment it because it failed when I run the code on the spyder console - child.expect("Enter name of param file: ") - child.sendline("phg.rec") - child.expect("Enter name of output file: ") - child.sendline(output) - child.expect("Enter the number of sub-samples*") - child.sendline(str(nrays)) - child.wait() - - CHUNK_SIZE = os.path.getsize(hdr_to_copy[0:-3] + "img") - print(CHUNK_SIZE) - with open(output, "rb") as f: - chunk = f.read(CHUNK_SIZE) - with open(output + ".img", "wb") as chunk_file: - chunk_file.write(chunk) - chunk_file.close() - - shutil.copy(hdr_to_copy, output + ".hdr") - - os.chdir(current_dir) +import random +import os +import shutil +import subprocess as sp +from os.path import join, exists +import pexpect +from utils import tools + +import nibabel as nib + + +def make_simset_act_table(act_table_factor, my_act_table, log_file=False): + with open(my_act_table, "w") as f: + f.write("256\n") + for i in range(256): + act_value = i * 0.000001 * act_table_factor + f.write(str("{:.12f}".format(act_value)) + "\n") + f.close() + + if log_file: + message = "Wrote act_table with factor %s" % act_table_factor + tools.log_message(log_file, message, "info") + + +def make_simset_phg( + config, + output_file, + simulation_dir, + act, + scanner_radius, + scanner_axial_fov, + center_slice, + photons, + sim_time, + add_randoms=False, + phg_hf=False, + S=0, + log_file=False, +): + stir_identifier = "# Hello, I am a SimSET PHG file!\n" + + # define types to be used + boolean = "BOOL " + real = "REAL " + longl = "LONGLONG " + integer = "INT " + enum = "ENUM " + string = "STR " + + # Configuring the importance sampling parameters + if add_randoms == 1: + stratification = "false" + forced_detection = "false" + non_absortion = "false" + # We need to force coincidences + singles for randoms + sim_PET_coinc_only = "false" + sim_PET_coinc_plus_singles = "true" + else: + stratification = config.get("stratification") + forced_detection = config.get("forced_detection") + non_absortion = config.get("forced_non_absortion") + sim_PET_coinc_only = "false" + sim_PET_coinc_plus_singles = "true" + + # Configuration of photons an simulation time + + num_to_simulate = str(photons) + length_of_scan = str(sim_time) + + # Import other configurations form config.yml + acceptance_angle = str(config.get("acceptance_angle")) + positron_range = config.get("positron_range") + non_col = config.get("non_colinearity") + minimum_energy = str(config.get("minimum_energy")) + weight_window_ratio = str(config.get("weight_window_ratio")) + point_source_voxels = config.get("point_source_voxels") + coherent_scatter_object = config.get("coherent_scatter_object") + coherent_scatter_detector = config.get("coherent_scatter_detector") + simulated_isotope = config.get("isotope") + + # Generate a random seed for the MC + random_seed = str(random.randint(1, 1e12)) + + # Object stuff + nslices = act.shape[2] + 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 + ) # 20 is because it must be cm + yMin, yMax = round(-act_fov[1, 1] / 20, 3), round(act_fov[1, 1] / 20, 3) + + z_offset = abs( + act.affine[2, 2] * center_slice / 10 + 0.5 * act.affine[2, 2] / 10 + ) # cm + zMin, zMax = round(-z_offset, 3), round(act_fov[2, 2] / 10 - z_offset, 3) + + dz = round((zMax - zMin) / nslices, 2) + + max_z_target = scanner_axial_fov / 2 + min_z_target = -scanner_axial_fov / 2 + + # Directory variables + simset_dir = config.get("dir_simset") + simset_phgdata_dir = join(simset_dir, "phg.data") + + # Now we write the phg file + with open(output_file, "w") as f: + # First we write the configuration stuff for the simulation + f.write(stir_identifier) + f.write("\n\n# Runtime options\n") + f.write(boolean + "simulate_stratification = %s \n" % stratification) + f.write(boolean + "simulate_forced_detection = " + forced_detection + "\n") + f.write(boolean + "forced_non_absorbtion = " + non_absortion + "\n") + f.write(real + "acceptance_angle = " + acceptance_angle + "\n") + f.write(longl + "num_to_simulate = " + num_to_simulate + "\n") + f.write(real + "length_of_scan = " + length_of_scan + "\n") + f.write(boolean + "simulate_SPECT = false\n") + f.write( + boolean + "simulate_PET_coincidences_only = " + sim_PET_coinc_only + "\n" + ) + f.write( + boolean + + "simulate_PET_coincidences_plus_singles = " + + sim_PET_coinc_plus_singles + + "\n" + ) + f.write(boolean + "adjust_for_positron_range = " + positron_range + "\n") + f.write(boolean + "adjust_for_collinearity = " + non_col + "\n") + f.write(real + "minimum_energy = " + minimum_energy + "\n") + f.write(real + "photon_energy = 511.0\n") + f.write(real + "weight_window_ratio = " + weight_window_ratio + "\n") + f.write(boolean + "point_source_voxels = " + point_source_voxels + "\n") + f.write(integer + "random_seed = " + random_seed + "\n") + f.write( + boolean + + "model_coherent_scatter_in_obj = " + + coherent_scatter_object + + "\n" + ) + f.write( + boolean + + "model_coherent_scatter_in_tomo = " + + coherent_scatter_detector + + "\n" + ) + f.write(enum + "isotope = " + simulated_isotope + "\n") + + # Now comes the object stuff (old make_phg_simset from fruitcake) + f.write("\n\n# OBJECT GEOMETRY VALUES\n") + f.write("\nNUM_ELEMENTS_IN_LIST object = %s" % str(nslices + 1)) + 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) + f.write( + "\n NUM_ELEMENTS_IN_LIST slice = 9 " + + "\n INT slice_number = %s" % str(i) + + "\n REAL zMin = %s" % str(zMin_value) + + "\n REAL zMax = %s" % str(zMax_value) + + "\n REAL xMin = %s" % str(xMin) + + "\n REAL xMax = %s" % str(xMax) + + "\n REAL yMin = %s" % str(yMin) + + "\n REAL yMax = %s" % str(yMax) + + "\n INT num_X_bins = %s" % str(xbins) + + "\n INT num_Y_bins = %s" % str(ybins) + ) + + # Now we fill the values for the target cylinder + + f.write("\n\n# TARGET CYLINDER INFORMATION") + f.write("\nNUM_ELEMENTS_IN_LIST target_cylinder = 3") + f.write("\n REAL target_zMin = %s" % str(min_z_target)) + f.write("\n REAL target_zMax = %s" % str(max_z_target)) + f.write("\n REAL radius = %s\n\n" % str(scanner_radius)) + + # Now we need the directory stuff + f.write( + string + + 'coherent_scatter_table = "' + + join(simset_phgdata_dir, "coh.tables") + + '"\n' + ) + f.write( + string + + 'activity_indexes = "' + + join(simulation_dir, "rec.act_indexes") + + '"\n' + ) + f.write( + string + + 'activity_table = "' + + join(simulation_dir, "phg_act_table") + + '"\n' + ) + f.write( + string + + 'activity_index_trans = "' + + join(simset_phgdata_dir, "phg_act_index_trans") + + '"\n' + ) + f.write( + string + + 'activity_image = "' + + join(simulation_dir, "rec.activity_image") + + '"\n' + ) + f.write( + string + + 'attenuation_indexes = "' + + join(simulation_dir, "rec.att_indexes") + + '"\n' + ) + f.write( + string + + 'attenuation_table = "' + + join(simset_phgdata_dir, "phg_att_table") + + '"\n' + ) + f.write( + string + + 'attenuation_index_trans = "' + + join(simset_phgdata_dir, "phg_att_index_trans") + + '"\n' + ) + f.write( + string + + 'attenuation_image = "' + + join(simulation_dir, "rec.attenuation_image") + + '"\n' + ) + + # This part changes productivity tables between S=0 and S=1 adq + if S == 1: + f.write( + string + + 'productivity_input_table = "' + + join(simulation_dir, "sampling_rec") + + '"\n' + ) + f.write(string + 'productivity_output_table = ""\n') + else: + f.write(string + 'productivity_input_table = ""\n') + f.write( + string + + 'productivity_output_table = "' + + join(simulation_dir, "sampling_rec") + + '"\n' + ) + + f.write( + string + 'statistics_file = "' + join(simulation_dir, "rec.stat") + '"\n' + ) + f.write( + string + + 'isotope_data_file = "' + + join(simset_phgdata_dir, "isotope_positron_energy_data") + + '"\n' + ) + f.write( + string + + 'detector_params_file = "' + + join(simulation_dir, "det.rec") + + '"\n' + ) + f.write( + string + 'bin_params_file = "' + join(simulation_dir, "bin.rec") + '"\n' + ) + # Adds the phg history file if it is activated + if phg_hf == 1: + f.write( + string + + 'history_file = "' + + join(simulation_dir, "phg_hf.hist" + '"\n') + ) + + f.close() + + if log_file: + message = "Created phg_file in: %s" % output_file + tools.log_message(log_file, message, "info") + + +def make_simset_bin( + config, output_file, simulation_dir, scanner, add_randoms=False, log_file=False, list_mode=False +): + stir_identifier = "# Hello, I am a SimSET BIN file!\n" + + # define types to be used + boolean = "BOOL " + real = "REAL " + integer = "INT " + string = "STR " + + # We will create an extra bin for randoms if the randoms simulation is activated. + if add_randoms: + accept_randoms = "true" + scatter_param = "6" + else: + accept_randoms = "false" + scatter_param = "1" + + min_s = "0" + max_s = "9" + + # Here we get the parameters from the parameter files. + num_z_bins = str(scanner.get("num_rings")) + axial_fov = scanner.get("axial_fov") + min_z, max_z = -axial_fov / 2, axial_fov / 2 + num_aa_bins = str(scanner.get("num_aa_bins")) + num_td_bins = str(scanner.get("num_td_bins")) + min_td = str(-scanner.get("scanner_radius")) + max_td = str(scanner.get("scanner_radius")) + min_e = str(scanner.get("min_energy_window")) + max_e = str(scanner.get("max_energy_window")) + rec_weight_file = join(simulation_dir, "rec.weight") + + with open(output_file, "w") as f: + f.write(stir_identifier) + f.write(boolean + " accept_randoms = " + accept_randoms + "\n") + f.write(integer + "scatter_param = " + scatter_param + "\n") + f.write(integer + "min_s = " + min_s + "\n") + f.write(integer + "max_s = " + max_s + "\n") + f.write(integer + "num_z_bins = " + num_z_bins + "\n") + f.write(real + "min_z = " + str(min_z) + "\n") + f.write(real + "max_z = " + str(max_z) + "\n") + f.write(integer + "num_aa_bins = " + num_aa_bins + "\n") + f.write(integer + "num_td_bins = " + num_td_bins + "\n") + f.write(real + "min_td = " + min_td + "\n") + f.write(real + "max_td = " + max_td + "\n") + f.write("INT num_e1_bins = 1\n") + f.write("INT num_e2_bins = 1\n") + f.write(real + "min_e = " + min_e + "\n") + f.write(real + "max_e = " + max_e + "\n") + f.write("INT weight_image_type = 2\n") #Originally was 3. Change to 2 for optimisation. + f.write("INT count_image_type = 2\n") #Originally was 3. Change to 2 for consistence with weight image. + f.write("BOOL add_to_existing_img = false\n") + f.write("BOOL sum_according_to_type = true\n") #Originally was false. Change to true for optimisation. + + if list_mode == False: #Uses the param if no listmode is used. Otherwise does not write weights (RAM saving). + f.write(string + 'weight_image_path = "' + rec_weight_file + '"\n') + f.close() + + if log_file: + message = "Created bin_file in: %s" % output_file + tools.log_message(log_file, message, "info") + + +def make_simset_simp_det(scanner_params, output, sim_dir, det_hf=0, log_file=False): + energy_resolution = scanner_params.get("energy_resolution") + new_file = open(output, "w") + new_file.write( + "ENUM detector_type = simple_pet \n\n" + + "REAL reference_energy_keV = 511.0 \n" + + "REAL energy_resolution_percentage = %s \n" % energy_resolution + ) + + if det_hf == 1: + new_file.write( + 'STR history_file = "' + join(sim_dir, "det_hf.hist" + '"\n') + ) + + new_file.close() + + if log_file: + message = ( + "Created det_file with:\n" + "Energy resolution: %s" % energy_resolution + ) + + tools.log_message(log_file, message, "info") + + +def make_simset_cyl_det(scanner_params, output, sim_dir, det_hf=0, log_file=False): + num_rings = scanner_params.get("num_rings") + z_crystal_size = scanner_params.get("z_crystal_size") + axial_fov = scanner_params.get("axial_fov") + max_z = axial_fov / 2 + min_z = -axial_fov / 2 + gap_z_size = (max_z - min_z - z_crystal_size * num_rings) / (num_rings - 1) + cyln_inner_radius = scanner_params.get("scanner_radius") + cyln_outer_radius = cyln_inner_radius + scanner_params.get("crystal_thickness") + energy_resolution = scanner_params.get("energy_resolution") + timing_resolution = scanner_params.get("timing_resolution") + material = scanner_params.get("simset_material") + + nrings_total = 2 * num_rings - 1 + + new_file = open(output, "w") + new_file.write( + "ENUM detector_type = cylindrical \n\n" + + "# This detector example has %s axial rings and %s gaps \n" + % (num_rings, num_rings - 1) + + "INT cyln_num_rings = %s \n\n" % nrings_total + ) + + for i in range(1, num_rings + 1): + ring_zmin = min_z + (i - 1) * z_crystal_size + (i - 1) * gap_z_size + ring_zmax = ring_zmin + z_crystal_size + gap_zmax = ring_zmin + z_crystal_size + gap_z_size + + new_file.write( + "# RING #%s \n" % i + + "# The following defines the ring parameters \n" + + "LIST cyln_ring_info_list = 5 \n" + + "INT cyln_num_layers = 1 \n" + + "LIST cyln_layer_info_list = 4 \n" + + "BOOL cyln_layer_is_active = TRUE \n" + + "INT cyln_layer_material = %s \n" % material + + "REAL cyln_layer_inner_radius = %s \n" % cyln_inner_radius + + "REAL cyln_layer_outer_radius = %s \n" % cyln_outer_radius + + "REAL cyln_min_z = %s \n" % ring_zmin + + "REAL cyln_max_z = %s \n\n" % ring_zmax + ) + + if not i == num_rings: + new_file.write( + "# GAP #%s \n" % i + + "# The following defines the gap parameters \n" + + "LIST cyln_ring_info_list = 5 \n" + + "INT cyln_num_layers = 1 \n" + + "LIST cyln_layer_info_list = 4 \n" + + "BOOL cyln_layer_is_active = FALSE \n" + + "INT cyln_layer_material = 0 \n" #NOTE: if needed to change the gap material, change it here: e.g., Lead (8), Aluminium (20). This has to be to 0 by default (air). + + "REAL cyln_layer_inner_radius = %s \n" % cyln_inner_radius + + "REAL cyln_layer_outer_radius = %s \n" % cyln_outer_radius + + "REAL cyln_min_z = %s \n" % ring_zmax + + "REAL cyln_max_z = %s \n\n" % gap_zmax + ) + + new_file.write( + "REAL reference_energy_keV = 511.0 \n" + + "REAL energy_resolution_percentage = %s \n" % energy_resolution + + "REAL photon_time_fwhm_ns = %s \n" % timing_resolution + ) + if det_hf == 1: + new_file.write( + 'STR history_file = "' + join(sim_dir, "det_hf.hist" + '"\n') + ) + + new_file.close() + + if log_file: + message = ( + "Created det_file with:\n" + + "Cristal z: %s cm\n" % z_crystal_size + + "Gap size: %s cm\n" % gap_z_size + + "Ring thickness: %s cm" % (cyln_outer_radius - cyln_inner_radius) + + "Energy resolution: %s" % energy_resolution + + "Timing resolution: %s" % timing_resolution + ) + + tools.log_message(log_file, message, "info") + + + +def make_index_file(simulation_dir, simset_dir, log_file=False): + output = join(simulation_dir, "index_file.log") + phg_file = join(simulation_dir, "phg.rec") + act_dat = join(simulation_dir, "act.dat") + att_dat = join(simulation_dir, "att.dat") + make_index_file = join(simset_dir, "bin", "makeindexfile") + + prompts = ( + phg_file + + "\n" + + "y\n" + + "y\n" + + "0\n" + + "y\n" + + act_dat + + "\n" + + "0\n" + + "0\n" + + "1\n" + + "n\n" + + "n\n" + + "y\n" + + "y\n" + + "0\n" + + "y\n" + + att_dat + + "\n" + + "0\n" + + "0\n" + + "1\n" + + "n\n" + + "n\n" + ) + + os.system('printf "%s" | %s > %s' % (prompts, make_index_file, output)) + + if log_file: + message = 'printf "%s" | %s > %s' % (prompts, make_index_file, output) + tools.log_message(log_file, message, "info") + + +def process_weights(weights_file, output_dir, scanner, add_randoms=0): + nbins = scanner.get("num_td_bins") + nangles = scanner.get("num_aa_bins") + nrings = scanner.get("num_rings") + nslices = nrings * nrings + + Simset_offset = 32768 + block_size = nbins * nangles * nslices * 4 + + trues_start = Simset_offset + trues_end = Simset_offset + block_size + trues_file = join(output_dir, "w1") + + #This is much faster than originally opening and writing a file: + with open(weights_file, "rb") as in_file: + with open(trues_file, "wb") as out_file: + in_file.seek(trues_start) + out_file.write(in_file.read(block_size)) + + + output = join(output_dir, "trues.nii") + + tools.create_nifti_from_imgdata( + trues_file, output, nbins, nangles, nslices, 1, 1, 1, "fl" + ) + os.remove(trues_file) + + + scatter_start = trues_end + scatter_end = trues_end + block_size + scatter_file = join(output_dir, "w2") + + #This is much faster than originally opening and writing a file: + with open(weights_file, "rb") as in_file: + with open(scatter_file, "wb") as out_file: + in_file.seek(scatter_start) + out_file.write(in_file.read(block_size)) + + output = join(output_dir, "scatter.nii") + + tools.create_nifti_from_imgdata( + scatter_file, output, nbins, nangles, nslices, 1, 1, 1, "fl" + ) + os.remove(scatter_file) + + #TODO: Still few things need to be tested in new version... + if add_randoms == 1: + randoms_start = scatter_end + randoms_end = scatter_end + block_size + randoms_file = join(output_dir, "w3") + + + #This is much faster than originally opening and writing a file: + with open(weights_file, "rb") as in_file: + with open(randoms_file, "wb") as out_file: + in_file.seek(randoms_start) + out_file.write(in_file.read(block_size)) + + output = join(output_dir, "randoms.nii") + + tools.create_nifti_from_imgdata( + randoms_file, output, nbins, nangles, nslices, 1, 1, 1, "fl" + ) + os.remove(randoms_file) + +def add_randoms(sim_dir, simset_dir, coincidence_window, rebin=True, log_file=False): + string = "STR " + integer = "INT " + enum = "ENUM " + real = "REAL " + + sorted_file = join(sim_dir, "sorted_det_hf.hist") + template = join(sim_dir, "sort.params") + det_hf = join(sim_dir, "det_hf.hist") + + with open(template, "w") as f: + f.write(string + 'history_file = "' + det_hf + '"\n') + f.write(string + 'sorted_history_file = "' + sorted_file + '"\n') + f.write(integer + "buffer_size = 1024\n") + + f.close() + + sorting_bin = join(simset_dir, "bin", "timesort -d") + + command = command = "%s %s >> %s" % (sorting_bin, template, log_file) + tools.osrun(command, log_file) + + randoms_file = join(sim_dir, "randoms.hist") + template = join(sim_dir, "randoms.params") + + with open(template, "w") as f: + f.write(enum + "detector_type = cylindrical\n") + f.write(string + 'history_file = "' + sorted_file + '"\n') + f.write(string + 'randoms_history_file = "' + randoms_file + '"\n') + f.write( + real + "coincidence_timing_window_in_ns = " + str(coincidence_window) + "\n" + ) + + f.close() + + addrand_bin = join(simset_dir, "bin", "addrandoms") + + command = "%s %s >> %s" % (addrand_bin, template, log_file) + tools.osrun(command, log_file) + + if exists(sorted_file): + os.remove(sorted_file) + + # The following will replace the existing det-hist file with the new including randoms + + if rebin: + det_file = join(sim_dir, "det.rec") + + with open(det_file, "r") as f: + lines = f.readlines() + with open(det_file, "w") as f: + for line in lines: + line = line.replace("det_hf.hist", "randoms.hist") + f.write(line) + + binfile = join(sim_dir, "phg.rec") + phgbin = join(simset_dir, "bin", "bin -d") + + command = command = "%s %s >> %s" % (phgbin, binfile, log_file) + tools.osrun(command, log_file) + + +def combine_history_files(simset_dir, history_files, output, log_file): + combinehist = join(simset_dir, "bin", "combinehist") + + rcommand = "%s %s %s" % (combinehist, history_files, output) + + proc = sp.Popen( + rcommand, + universal_newlines=True, + shell=True, + stdin=sp.PIPE, + stdout=sp.PIPE, + stderr=sp.PIPE, + ).communicate("No\n") #This was originally to Yes, set to No (delete inputs?) for optimisation.. + + +def simset_calcattenuation( + simset_dir, sim_dir, output, hdr_to_copy, nrays=1, timeout=36000 +): + calcattenuation = join(simset_dir, "bin", "calcattenuation") + + current_dir = os.getcwd() + os.chdir(sim_dir) + + child = pexpect.spawn(calcattenuation, timeout=timeout) + # child.logfile = sys.stdout.buffer #I comment it because it failed when I run the code on the spyder console + child.expect("Enter name of param file: ") + child.sendline("phg.rec") + child.expect("Enter name of output file: ") + child.sendline(output) + child.expect("Enter the number of sub-samples*") + child.sendline(str(nrays)) + child.wait() + + + + + nib.save(nib.load(hdr_to_copy), hdr_to_copy[0:-3] + "hdr") + + CHUNK_SIZE = os.path.getsize(hdr_to_copy[0:-3] + "img") + + with open(output, "rb") as f: + chunk = f.read(CHUNK_SIZE) + with open(output + ".img", "wb") as chunk_file: + chunk_file.write(chunk) + chunk_file.close() + + + shutil.copy(hdr_to_copy[0:-3] + "hdr", output + ".hdr") + + nib.save(nib.load(output + ".hdr"), output + ".nii") + + #Remove unnecessary files: + os.remove(hdr_to_copy[0:-3] + "hdr") + os.remove(hdr_to_copy[0:-3] + "img") + os.remove(output + ".hdr") + os.remove(output + ".img") + + os.chdir(current_dir) diff --git a/src/stir/stir_tools.py b/src/stir/stir_tools.py index ee8ffe94..13b6cab8 100644 --- a/src/stir/stir_tools.py +++ b/src/stir/stir_tools.py @@ -8,7 +8,7 @@ from utils import tools from utils import resources as rsc - +#ATTENTION: Old version, must be updated for STIR >6.0 def create_stir_hs_from_detparams(scannerParams, output_file, output_format="SimSET"): num_rings = scannerParams.get("num_rings") max_z = scannerParams.get("axial_fov") / 2 @@ -16,7 +16,9 @@ def create_stir_hs_from_detparams(scannerParams, output_file, output_format="Sim z_crystal_size = scannerParams.get("z_crystal_size") gap_size = (max_z - min_z - z_crystal_size * num_rings) / (num_rings - 1) ring_spacing = z_crystal_size + gap_size - + + max_segment = scannerParams.get("max_segment") + min_td = -scannerParams.get("scanner_radius") max_td = scannerParams.get("scanner_radius") td_bins = scannerParams.get("num_td_bins") @@ -24,6 +26,9 @@ def create_stir_hs_from_detparams(scannerParams, output_file, output_format="Sim matrix_size, ring_difference = generate_segments_lists_stir( num_rings, num_rings - 1 ) + #matrix_size, ring_difference = generate_segments_lists_stir( #NOTE: This will be needed for FBP + # num_rings, max_segment + #) transaxial_crystal_distance = scannerParams.get("transaxial_crystal_size") scanner_name = scannerParams.get("scanner_name") @@ -60,7 +65,7 @@ def create_stir_hs_from_detparams(scannerParams, output_file, output_format="Sim + "number of dimensions := 4\n" + "matrix axis label [4] := segment\n" + "!matrix size [4] := " - + str(2 * (num_rings - 1) + 1) + + str(2 * (num_rings - 1) + 1) #WARNING: For FBP we need to do: str(2 * (max_segment) + 1) + "\n" + "matrix axis label [" + str(views_coordinate) @@ -161,17 +166,43 @@ def generate_segments_lists_stir(nrings, max_segment): my_matrix_ring_difference = " " for i in range(last_segment_sinograms, nrings): my_matrix_size = my_matrix_size + str(i) + "," - my_matrix_ring_difference = my_matrix_ring_difference + str(i - nrings) + "," + my_matrix_ring_difference = my_matrix_ring_difference + str(i - nrings) + "," + #my_matrix_ring_difference = my_matrix_ring_difference + str(nrings - i) + "," #NOTE: in case rings are invented. for i in range(nrings, last_segment_sinograms - 1, -1): my_matrix_size = my_matrix_size + str(i) + "," + my_matrix_ring_difference = my_matrix_ring_difference + str(nrings - i) + "," + #my_matrix_ring_difference = my_matrix_ring_difference + str(i - nrings) + "," #NOTE: in case rings are invented. + + my_matrix_size = "{" + my_matrix_size[0:-1] + "}" + my_matrix_ring_difference = "{" + my_matrix_ring_difference[0:-1] + "}" + + return my_matrix_size, my_matrix_ring_difference + + +def generate_segments_lists_stir_v2(nrings, max_segment): + my_matrix_size = " " + for i in range(1, nrings): + my_matrix_size = my_matrix_size + str(i) + "," + for i in range(nrings, 0, -1): + my_matrix_size = my_matrix_size + str(i) + "," + + last_segment_sinograms = nrings - max_segment + my_matrix_ring_difference = " " + for i in range(last_segment_sinograms, nrings): + my_matrix_ring_difference = my_matrix_ring_difference + str(i - nrings) + "," + for i in range(nrings, last_segment_sinograms - 1, -1): my_matrix_ring_difference = my_matrix_ring_difference + str(nrings - i) + "," my_matrix_size = "{" + my_matrix_size[0:-1] + "}" my_matrix_ring_difference = "{" + my_matrix_ring_difference[0:-1] + "}" + + print("MY MATRIX SIZE:", my_matrix_size) + print("MY MATRIX RING DIFFERENCE:", my_matrix_ring_difference) return my_matrix_size, my_matrix_ring_difference + def apply_psf(scannerParams, sinogram_stir, log_file): num_z_bins = scannerParams.get("num_rings") num_aa_bins = scannerParams.get("num_aa_bins") @@ -291,7 +322,9 @@ def FBP2D_recons(config, scannerParams, sinograms_stir, output_dir, log_file): tools.osrun(command, log_file) output = recFileName + ".hv" - output = tools.anything_to_hdr_convert(output, log_file) + + output = tools.convert_hv_to_nii(output, log_file) + tools.flip_rec_nifti(recFileName + ".nii", join(output_dir, "reconstruction.nii")) #TODO: We need to check the flips done here. Final reconstruction is flipped from maps (due to sinogram flipping). return output @@ -302,6 +335,7 @@ def FBP3D_recons(config, scannerParams, sinograms_stir, output_dir, log_file): zoom = scannerParams.get("zoomFactor") xyOutputSize = scannerParams.get("xyOutputSize") + max_segment = scannerParams.get("max_segment") recFileName = join(output_dir, "rec_FBP3D") @@ -321,20 +355,23 @@ def FBP3D_recons(config, scannerParams, sinograms_stir, output_dir, log_file): + "xy output image size (in pixels) := " + str(xyOutputSize) + "\n\n" - "maximum absolute segment number to process := 2 \n" + #"maximum absolute segment number to process := 2 \n" + #+ "maximum absolute segment number to process := -1 \n" + + "maximum absolute segment number to process := " + str(max_segment) + "\n" + "num segments to combine with ssrb := -1 \n\n" + #+ "num segments to combine with ssrb := 3 \n\n" + "alpha parameter for ramp filter := 1 \n" + "cut-off for ramp filter (in cycles) := 0.5 \n\n" - + "alpha parameter for colsher filter in axial direction := 1 \n" - + "cut-off for colsher filter in axial direction (in cycles) := 0.5 \n" - + "alpha parameter for colsher filter in planar direction := 1 \n" - + "cut-off for colsher filter in planar direction (in cycles) := 0.5 \n\n" - + "stretch factor for colsher filter definition in axial direction := 2 \n" - + "stretch factor for colsher filter definition in planar direction := 2 \n\n" - + "transaxial extension for fft := 1 \n" - + "axial extension for fft := 1 \n\n" - + "save intermediate images := 0 \n" - + "display level := 0 \n\n" + #+ "alpha parameter for colsher filter in axial direction := 1 \n" + #+ "cut-off for colsher filter in axial direction (in cycles) := 0.5 \n" + #+ "alpha parameter for colsher filter in planar direction := 1 \n" + #+ "cut-off for colsher filter in planar direction (in cycles) := 0.5 \n\n" + #+ "stretch factor for colsher filter definition in axial direction := 2 \n" + #+ "stretch factor for colsher filter definition in planar direction := 2 \n\n" + #+ "transaxial extension for fft := 1 \n" + #+ "axial extension for fft := 1 \n\n" + #+ "save intermediate images := 0 \n" + #+ "display level := 0 \n\n" + "end := \n" ) @@ -344,7 +381,9 @@ def FBP3D_recons(config, scannerParams, sinograms_stir, output_dir, log_file): tools.osrun(command, log_file) output = recFileName + ".hv" - output = tools.anything_to_hdr_convert((output)) + + output = tools.convert_hv_to_nii(output, log_file) + tools.flip_rec_nifti(recFileName + ".nii", join(output_dir, "reconstruction.nii")) #WARNING! REMOVE IF THIS FAILS. THIS FLIPS THE FINAL REC IN 0 AXIS return output @@ -558,8 +597,19 @@ def OSEM3D_recons( numberOfSubsets = scannerParams.get("numberOfSubsets") numberOfIterations = scannerParams.get("numberOfIterations") savingInterval = scannerParams.get("savingInterval") - - if scannerParams.get("stir_recons_att_corr") == 1: + + #TODO: To unify these blocks if stir_att file with same name is generated with both methods. + if scannerParams.get("attenuation_mode") == 1: + att_corr_str = ( + "Bin Normalisation type := From ProjData \n" + + "Bin Normalisation From ProjData := \n" + + "normalisation projdata filename:= " + + att_stir + + "\n" + + "End Bin Normalisation From ProjData:= \n" + ) + + elif scannerParams.get("attenuation_mode") == 2: att_corr_str = ( "Bin Normalisation type := From ProjData \n" + "Bin Normalisation From ProjData := \n" @@ -570,7 +620,8 @@ def OSEM3D_recons( ) else: att_corr_str = "" - + + if ( scannerParams.get("stir_scatt_corr_smoothing") == 1 ): # Will use smoothed SimSET scatter as additive_sinogram. @@ -629,7 +680,7 @@ def OSEM3D_recons( + "Projector Pair Using Matrix Parameters := \n" + "Matrix type := Ray Tracing \n" + "Ray tracing matrix parameters := \n" - + "number of rays in tangential direction to trace for each bin:= 5 \n" + + "number of rays in tangential direction to trace for each bin:= 5 \n" #WARNING! ORIGINAL WAS 5 + "End Ray tracing matrix parameters := \n" + "End Projector Pair Using Matrix Parameters := \n\n" + att_corr_str @@ -672,12 +723,13 @@ def OSEM3D_recons( ) new_file.close() - + command = "%s %s >> %s" % (recons, paramsFile, log_file) tools.osrun(command, log_file) - + output = recFileName + "_" + str(scannerParams.get("numberOfIterations")) + ".hv" - output = tools.anything_to_hdr_convert(output, log_file) + output = tools.convert_hv_to_nii(output, log_file) + tools.flip_rec_nifti(recFileName + "_" + str(scannerParams.get("numberOfIterations")) + ".nii", join(output_dir, "reconstruction.nii")) #WARNING! REMOVE IF THIS FAILS. THIS FLIPS THE FINAL REC IN 0 AXIS return output diff --git a/submodules/STIR/STIR b/submodules/STIR/STIR index cab371bc..166ba624 160000 --- a/submodules/STIR/STIR +++ b/submodules/STIR/STIR @@ -1 +1 @@ -Subproject commit cab371bcb558a6bc202741f08d4d28b89a816682 +Subproject commit 166ba624e4c00af20612d07c94f1d2f7f04eedd6 diff --git a/submodules/simset_simpet b/submodules/simset_simpet index f036934c..54840ac3 160000 --- a/submodules/simset_simpet +++ b/submodules/simset_simpet @@ -1 +1 @@ -Subproject commit f036934ce199fe0b39201a419223b938bf9540cd +Subproject commit 54840ac3db52621aaafb953ff3c3f02ece9ac1be diff --git a/utils/tools.py b/utils/tools.py index 6f411d04..3bb5edf5 100644 --- a/utils/tools.py +++ b/utils/tools.py @@ -1,1011 +1,1510 @@ -import os, shutil, datetime -import re -from os.path import join, exists, isdir, dirname, basename, splitext -from subprocess import getstatusoutput as getoutput -import nibabel as nib -import nibabel.processing as nibp -from utils import resources as rsc -from utils import spm_tools as spm -import numpy as np -from operator import itemgetter -from nipype.interfaces.dcm2nii import Dcm2nii -from nipype.interfaces import fsl - - -def osrun(command, logfile, catch_out=False): - """ - Executes command, raise an error if it fails and send status to logger - :param command: command to be executed - :param logger: logger file - :return: - """ - if catch_out: - status, out = getoutput(command) - if status != 0: - log_message(logfile, command, 'error') - raise TypeError(command) - else: - log_message(logfile, command) - return out - - else: - if os.system(command) != 0: - - log_message(logfile, command, 'error') - raise TypeError(command) - else: - log_message(logfile, command) - -def nib_load(image, logfile=False): - """ - Load image and return data array - :param image: image to load data from - :return: image data array - """ - try: - img = nib.load(image) - data = img.get_data()[:,:,:] - return img, data - except Exception as e: - message = "Error: " + str(e) - if logfile: - log_message(logfile, message, 'error') - else: - print(message) - -def copy_analyze(image1, image2=False, dest_dir=False, logfile=False): - """ - Create a copy of an Analyze format image - :param image1: (string) path to the original image - :param image2: (string, optional) path to the copy image - :param dest_dir: (string, optional) path to the destination folder - :return: - """ - - if image2: - if image1[-3:] == 'hdr' or image1[-3:] == 'img' and image2[-3:] == 'hdr' or image2[-3:] == 'img': - image1_hdr = image1[0:-3] + 'hdr' - image2_hdr = image2[0:-3] + 'hdr' - shutil.copy(image1_hdr, image2_hdr) - image1_img = image1[0:-3] + 'img' - image2_img = image2[0:-3] + 'img' - shutil.copy(image1_img, image2_img) - - return image2_hdr - else: - message = 'Error! The provided image is not in Analyze format:' + str(image1) - if logfile: log_message(logfile, message, 'error') - raise TypeError(message) - - elif not image2 and isdir(str(dest_dir)): - - ext = splitext(basename(image1))[1] - - if ext == '.img' or ext == '.hdr': - image1_hdr = image1[0:-3] + 'hdr' - image2_hdr = join(dest_dir, basename(image1)[0:-3] + 'hdr') - shutil.copy(image1_hdr, image2_hdr) - image1_img = image1[0:-3] + 'img' - image2_img = join(dest_dir, basename(image1)[0:-3] + 'img') - shutil.copy(image1_img, image2_img) - - return image2_hdr - else: - message = 'Error! The provided image is not in Analyze format:' + str(image1) - if logfile: log_message(logfile, message, 'error') - raise TypeError(message) - else: - message = 'Error! Not image2 path or dest dir path provided: ' + str(image2) + ', ' + str(dest_dir) - if logfile: log_message(logfile, message, 'error') - raise TypeError(message) - -def create_analyze_from_imgdata(data, out, pix_x, pix_y, pix_z, tx, ty, tz, data_type="fl"): - - if data_type == "1b": - dtype=np.int8 - elif data_type == "2b": - dtype=np.int16 - elif data_type == "db": - dtype=np.float64 - else: - dtype=np.float32 - - hdr1 = nib.AnalyzeHeader() - hdr1.set_data_dtype(dtype) - hdr1.set_data_shape((pix_x,pix_y,pix_z)) - hdr1.set_zooms((tx,ty,tz)) - - f = open(data,'rb') - - img_data = hdr1.raw_data_from_fileobj(f) - - analyze_img = nib.AnalyzeImage(img_data, hdr1.get_base_affine(), hdr1) - - nib.save(analyze_img,out) - -def read_analyze_header(header_file,logfile): - - img = nib.load(header_file) - - zpix = img.shape[2] - #zsize = abs(img.get_affine[2,2]) - zsize = abs(img.affine[2,2]) - xpix = img.shape[0] - #xsize = abs(img.get_affine[0,0]) - xsize = abs(img.affine[0,0]) - ypix = img.shape[1] - #ysize = abs(img.get_affine[1,1]) - ysize = abs(img.affine[1,1]) - - return zpix, zsize, xpix, xsize, ypix, ysize - -def write_interfile_header(header_file,matrix_size_x,pixel_size_x, matrix_size_y,pixel_size_y, matrix_size_z,pixel_size_z, offset_z=0): - - image_v = os.path.basename(header_file)[0:-2] + "v" - fheader_hv = open(header_file, "w") - - offset_x = matrix_size_x/2*pixel_size_x + 0.5*pixel_size_x - offset_y = matrix_size_y/2*pixel_size_y + 0.5*pixel_size_y - - fheader_hv.write("!INTERFILE :=\n") - fheader_hv.write("name of data file := %s\n" % image_v) - fheader_hv.write("!GENERAL DATA :=\n") - fheader_hv.write("!GENERAL IMAGE DATA :=\n") - fheader_hv.write("!type of data := PET\n") - fheader_hv.write("imagedata byte order := LITTLEENDIAN\n") - fheader_hv.write("!PET STUDY (General) :=\n") - fheader_hv.write("!PET data type := Image\n") - fheader_hv.write("process status := Reconstructed\n") - fheader_hv.write("!number format := float\n") - fheader_hv.write("!number of bytes per pixel := 4\n") - fheader_hv.write("number of dimensions := 3\n") - fheader_hv.write("matrix axis label [1] := x\n") - fheader_hv.write("!matrix size [1] := %s\n" % matrix_size_x) - fheader_hv.write("scaling factor (mm/pixel) [1] := %s\n" % pixel_size_x) - fheader_hv.write("matrix axis label [2] := y\n") - fheader_hv.write("!matrix size [2] := %s\n" % matrix_size_y) - fheader_hv.write("scaling factor (mm/pixel) [2] := %s\n" % pixel_size_y) - fheader_hv.write("matrix axis label [3] := z\n") - fheader_hv.write("!matrix size [3] := %s\n" % matrix_size_z) - fheader_hv.write("scaling factor (mm/pixel) [3] := %s\n" % pixel_size_z) - fheader_hv.write("first pixel offset (mm) [1] := -%s\n" % offset_x) - fheader_hv.write("first pixel offset (mm) [2] := -%s\n" % offset_y) - fheader_hv.write("first pixel offset (mm) [3] := 0\n") - fheader_hv.write("number of time frames := 1\n") - fheader_hv.write("!END OF INTERFILE :=\n") - - fheader_hv.close() - -def nii_analyze_convert(image, logfile=False, outfile=False): - """ - Converts the provided image file to the analyze or nifti format - depending on the filex extension - :param image: image to be converted - :return: - """ - - # if not logfile: logfile = join(dirname(image), 'log_nii_analyze_convert.txt') - - ext = os.path.splitext(image)[1] - img, data = nib_load(image) - data = data.astype(np.float32) #casting to avoid problems when applying nib.save with no float data - hdr = nib.AnalyzeHeader() - hdr.set_data_dtype(np.float32) - hdr.set_data_shape(data.shape) - - if ext == '.nii': - # n2a = rsc.get_rsc('nii2analyze', 'exe') - if outfile: - cimage = outfile - else: - cimage = image.replace('.nii', '.hdr') - # rcommand = '%s %s %s >> %s' % (n2a, image, cimage, logfile) - # osrun(rcommand, logfile) - imageToWrite = nib.AnalyzeImage(data, img.affine, hdr) - elif ext == '.hdr': - # a2n = rsc.get_rsc('analyze2nii', 'exe') - if outfile: - cimage = outfile - else: - cimage = image.replace('.hdr', '.nii') - # rcommand = '%s %s %s >> %s' % (a2n, image, cimage, logfile) - # osrun(rcommand, logfile) - imageToWrite = nib.Nifti1Image(data, img.affine, hdr) - else: - # a2n = rsc.get_rsc('analyze2nii', 'exe') - if outfile: - cimage = outfile - else: - cimage = image.replace('.img', '.nii') - # image_hdr = str(image).replace(".img", ".hdr") - # rcommand = '%s %s %s >> %s' % (a2n, image_hdr, cimage, logfile) - # osrun(rcommand, logfile) - imageToWrite = nib.Nifti1Image(data, img.affine, hdr) - - nib.save(imageToWrite, cimage) - - if exists(cimage): - if outfile: - shutil.move(cimage, outfile) - return outfile - else: - return cimage - else: - return 0 - -def anything_to_hdr_convert(image, logfile=False, outfile=False ): - """ - This function will try to convert dicom, .nii.gz or .nii to analyze... - """ - - if image[-3:] == "hdr" or image[-3:] == "img": - return image[0:-3]+"hdr" - - elif image[-3:] == "nii": - hdr = nii_analyze_convert(image,logfile=logfile) - if exists(hdr): - return hdr - else: - raise TypeError ("nifti-analyze conversion failed....") - - elif image[-6:]=="nii.gz": - hdr = image[0:-6] + 'hdr' - nii_analyze_convert(image,logfile=logfile, outfile=hdr) - if exists(hdr): - return hdr - else: - raise TypeError ("nifti-analyze conversion failed....") - - elif isdir(image): - print("I think the provided image:\n %s \n is a dicom directory. I will try to convert it...") - dicom2nii = rsc.get_rsc("dicom2nii","exe") - nii = join(dirname(image),"image.nii") - rcommand = '%s %s %s' % (dicom2nii, image, nii) - osrun(rcommand,logfile) - if exists(nii): - hdr = nii_analyze_convert(nii,logfile=logfile) - os.remove(nii) - if exists(hdr): - return hdr - else: - raise TypeError ("nifti-analyze conversion failed....") - else: - raise TypeError ("dicom-nifti conversion failed....") - - elif image[-2:] == "hv": #Only works with floats data - - with open(image) as f: - lines = f.readlines() - - lines = [x.strip() for x in lines] - - pixel = [re.compile(f"\!matrix size \[{i}\].*") for i in range(1, 4)] - pixel_size = [re.compile(f"scaling factor \(mm.pixel\) \[{i}\].*") for i in range(1, 4)] - - pixel_x = list(filter(pixel[0].match, lines))[0].split()[-1] - pixel_size_x = list(filter(pixel_size[0].match, lines))[0].split()[-1] - pixel_y = list(filter(pixel[1].match, lines))[0].split()[-1] - pixel_size_y = list(filter(pixel_size[1].match, lines))[0].split()[-1] - pixel_z = list(filter(pixel[2].match, lines))[0].split()[-1] - pixel_size_z = list(filter(pixel_size[2].match, lines))[0].split()[-1] - - hdr_header = image[0:-2] + "hdr" - # img_file = image[0:-2] + "img" - data_file = image[0:-2] + "v" - - #This will convert .hv to .hdr and copy the data - # gen_hdr = rsc.get_rsc("gen_hdr", "fruitcake") - # rcommand = "%s %s %s %s %s fl %s %s %s 0" % (gen_hdr, hdr_header[0:-4], pixel_x, pixel_y, pixel_z, pixel_size_x, pixel_size_y, pixel_size_z) - # osrun(rcommand, logfile) - # shutil.copy(data_file, img_file) - create_analyze_from_imgdata(data_file,hdr_header,float(pixel_x),float(pixel_y),float(pixel_z),float(pixel_size_x),float(pixel_size_y),float(pixel_size_z)) - #nii_analyze_convert("inicial.hdr",logfile,"aux.nii") - #nii_analyze_convert("aux.nii",logfile,hdr_header) - #os.remove("aux.nii") - - if exists(hdr_header): - return hdr_header - else: - raise TypeError ("Interfile-analyze conversion failed....") - -def prepare_input_image(image_hdr, logfile, min_voxel_size=1): - """ - This method converts input_image to float data type, re-sizes the image to - 1mm size voxels if too large to keep a reasonable analysis execution time and - removes the negative and NaN values. - :param input_image: image to prepare for the analysis - :return: - """ - - # Tools for image manipulation - # change_format = rsc.get_rsc('change_format', 'fruitcake') - # change_matrix = rsc.get_rsc('change_img_matrix', 'fruitcake') - # erase_negs = rsc.get_rsc('erase_negs', 'fruitcake') - # erase_nans = rsc.get_rsc('erase_nans', 'fruitcake') - - # Convert input image to Analyze format and float data type - # rcommand = '%s %s %s fl >> %s' % (change_format, image_hdr, image_hdr, logfile) - # osrun(rcommand, logfile) - - ## # CONTEMPLAR POSIBILIDAD DE CONVERTIR ESTO EN FUNCIÓN PARA SUSTITUIR AL "cambia_formato_hdr" de fruitcake - change_format(image_hdr, "fl", logfile) - - ### - - # Resize image if necessary - # ndims = recalculate_matrix(image_hdr, min_voxel_size) - # rcommand = '%s %s %s %s %s %s novecino >> %s' % (change_matrix, image_hdr, image_hdr, - # str(ndims[0]), str(ndims[1]), str(ndims[2]), logfile) - # osrun(rcommand, logfile) - image_hdr = resampleXYvoxelSizes(image_hdr, min_voxel_size, logfile) - image_hdr = resampleZvoxelSize(image_hdr, min_voxel_size, logfile) - - # Erase negative and NaN values - # rcommand = '%s %s %s >> %s' % (erase_negs, image_hdr, image_hdr, logfile) - # osrun(rcommand, logfile) - # rcommand = '%s %s %s >> %s' % (erase_nans, image_hdr, image_hdr, logfile) - # osrun(rcommand, logfile) - remove_neg_nan(image_hdr) - - - return image_hdr - -def recalculate_matrix(input_image, voxelsize, mode="downsampling"): - - new_dimensions = [] - image_load, data = nib_load(input_image) - header = image_load.header - sizes = header.get_zooms() - dimensions = image_load.shape - x_lenght = int(sizes[0]*dimensions[0]/voxelsize) - y_lenght = int(sizes[1]*dimensions[1]/voxelsize) - z_lenght = int(sizes[2]*dimensions[2]/voxelsize) - if mode == "downsampling": - if sizes[0]voxelsize: - new_dimensions.append(x_lenght) - else: - new_dimensions.append(dimensions[0]) - if sizes[1]>voxelsize: - new_dimensions.append(y_lenght) - else: - new_dimensions.append(dimensions[1]) - if sizes[2]>voxelsize: - new_dimensions.append(z_lenght) - else: - new_dimensions.append(dimensions[2]) - - return new_dimensions - -def verify_roi_exists(rois_image, roi_number): - """ - Compute the number of voxels in each ROI - :param rois_image: (string) path to the ROI parcelled image - :param rois_list: (array) ROIs indexing - :return: False if the number of voxels is zero - :return: nvox is not zero - """ - rois_img = nib.load(rois_image) - rois_data = rois_img.get_data()[:, :, :] - indx = np.where(rois_data == roi_number) - nvox = len(rois_data[indx]) - if nvox == 0: - return False - else: - return True - -def operate_single_image(input_image, operation, factor, output_image, logfile): - """ - Given an input image, multiply or divide it by a numerical factor - saving the result as output_image - :param input_image: image base operation on - #:param operation: 1 = multiply, 2 = divide - :param operation: 'mult' = multiply, 'div' = divide - :param factor: operation factor - :param output_image: output image file - :return: - """ - - img = nib.load(input_image) - data = img.get_data()[:,:,:] - data = np.nan_to_num(data) - - - if operation == 'mult': - data = data * float(factor) - elif operation == 'div': - data = data / float(factor) - else: - message = "Error! Invalid operation: " +str(operation) - print(message) - log_message(logfile, message, 'error') - - hdr1 = nib.AnalyzeHeader() - hdr1.set_data_dtype(img.get_data_dtype()) - hdr1.set_data_shape(img.shape) - #hdr1.set_zooms(abs(np.diag(img.affine))) - #hdr1.set_zooms(abs(np.diag(img.affine))[0:3]) - #hdr1.set_zooms(abs(np.diag(img.affine))[0:4]) - hdr1.set_zooms(abs(np.diag(img.affine))[0:img.ndim]) - - analyze_img = nib.AnalyzeImage(data, hdr1.get_base_affine(), hdr1) - - nib.save(analyze_img,output_image) - - -def operate_images_analyze(image1, image2, out_image, operation='mult'): - """ - Given the input images, calculate the multiplication image or the ratio between them - :param image1: string, path to the first image - :param image2: string, path to the second image - :param operation: string, multi (default) for multiplication divid for division - :param out_image: string (optional), path to the output image - :return: - """ - img1, data1 = nib_load(image1) - img2, data2 = nib_load(image2) - - # TODO CHECK IF NEGATIVE VALUES NEED TO BE REMOVED - # Remove NaN and negative values - data1 = np.nan_to_num(data1) - data2 = np.nan_to_num(data2) - - if operation == 'mult': - res_data = data1 * data2 - elif operation == 'div': - res_data = data1 / data2 - elif operation == 'sum': - res_data = data1 + data2 - elif operation == 'diff': - res_data = data1 - data2 - else: - message = 'Error! Unknown operation: ' + str(operation) - raise TypeError(message) - - hdr1 = nib.AnalyzeHeader() - hdr1.set_data_dtype(img1.get_data_dtype()) - hdr1.set_data_shape(img1.shape) - ##hdr1.set_zooms(abs(np.diag(img1.affine))[0:3]) - ##hdr1.set_zooms(abs(np.diag(img1.affine))[0:4]) - hdr1.set_zooms(abs(np.diag(img1.affine))[0:img1.ndim]) - - analyze_img = nib.AnalyzeImage(res_data, hdr1.get_base_affine(), hdr1) - - nib.save(analyze_img,out_image) - -def smooth_analyze(image,fwhm, output): - - from nibabel import processing as nibproc - - img =nib.load(image) - smoothed=nibproc.smooth_image(img,fwhm,out_class=nib.AnalyzeImage) - nib.save(smoothed,output) - -def log_message(logfile, message, mode='info'): - """" - Print outs logger messages to the specified logfile - """ - - #ctime = read_time(datetime.datetime.now()) - ctime = datetime.datetime.now() - - separator = '\n###########################################################' - - if mode == 'exe': - stream = separator + '\nTIME: ' + str(ctime) + '\n\nEXE: ' + message + separator + '\n' - elif mode == 'info': - stream = separator + '\nTIME: ' + str(ctime) + '\n\nINFO: ' + message + separator + '\n' - elif mode == 'warning': - stream = separator + '\nTIME: ' + str(ctime) + '\n\nWARNING: ' + message + separator + '\n' - elif mode == 'error': - stream = separator + '\nTIME: ' + str(ctime) + '\n\nERROR: ' + message + separator + '\n' - else: - stream = separator + '\nTIME: ' + str(ctime) + '\n\n' + message + '\n' - - if exists(logfile): - with open(logfile, 'a') as lfile: - lfile.write(stream) - else: - with open(logfile, 'w') as lfile: - lfile.write(stream) - -def reorient_dcmtonii(image_path): - - - components=os.path.split(image_path) - path = components[0] - image = components[1] - converter = Dcm2nii() - converter.inputs.source_names=[image_path] - converter.inputs.reorient_and_crop=True - # converter.nii_output=False - converter.run() - - reor_ima_name = "o"+image - if exists(join(path,reor_ima_name)): - shutil.copy(join(path,reor_ima_name),image_path) - os.remove(join(path,reor_ima_name)) - os.remove(join(path,"co"+image)) - else: - # converter.inputs_reorient_and_crop=False - # converter.run() - # 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. - Required inputs are: - pet_image: dicom_dir, nii.gz, nii or Analyze (.hdr or .img) - mri_image: dicom_dir, nii.gz, nii or Analyze (.hdr or .img) - mode: Choose STIR or SIMSET. The maps will be different for each simulation. - The inputs will be stored to Data/simulation_name/Patient as reformatted/corregistered Analyze - #The maps will be stored on Data/simulation_name/Maps - The maps will be stored on Results/output_name/Maps - """ - message = "GENERATING ACT AND ATT MAPS FROM PET, (CT) and MR IMAGES" - log_message(log_file, message, mode='info') - - reorient_dcmtonii(pet_image) - reorient_dcmtonii(mri_image) - reorient_dcmtonii(ct_image) - #First of all lets take all to analyze - pet_hdr = anything_to_hdr_convert(pet_image) - #pet_hdr = copy_analyze(pet_hdr,image2=False,dest_dir=patient_dir) - pet_hdr = prepare_input_image(pet_hdr,log_file,min_voxel_size=1.5) - pet_img = pet_hdr[0:-3]+"img" - - if ct_image: #ct_image is not empty - ct_hdr = anything_to_hdr_convert(ct_image) - ct_hdr = prepare_input_image(ct_hdr,log_file,min_voxel_size=1.5) - ct_img = ct_hdr[0:-3]+"img" - - mri_hdr = anything_to_hdr_convert(mri_image) - #mri_hdr = copy_analyze(mri_hdr,image2=False,dest_dir=patient_dir) - mri_hdr = prepare_input_image(mri_hdr,log_file,min_voxel_size=1.5) - mri_img = mri_hdr[0:-3]+"img" - - #make mri image square - makeImageSquare(mri_hdr, log_file) - - #Performing PET-CT/MR coregister - mfile = os.path.join(output_dir,"fusion_pet_to_mri.m") - # correg_pet_hdr = fsl_flirt(mri_hdr, pet_hdr, log_file) - # correg_pet_img =correg_pet_hdr[0:-3]+"img" - correg_pet_img = spm.image_fusion(spm_run, mfile, mri_img, pet_img, log_file) - correg_ct_img="" - if ct_image: #ct_image is not empty - mfile = os.path.join(output_dir,"fusion_ct_to_mri.m") - correg_ct_img = spm.image_fusion(spm_run, mfile, mri_img, ct_img, log_file) - - #Now the map generation - from utils.patient2maps import patient2maps - - my_map_generation = patient2maps(spm_run, output_dir, log_file, - mri_img, correg_pet_img, correg_ct_img, mode=mode) - activity_map_hdr, attenuation_map_hdr = my_map_generation.run() - - return activity_map_hdr, attenuation_map_hdr - -def convert_map_values(act_map,att_map,output_dir,log_file,mode="SimSET"): - - message = "CONVERTING ACT AND ATT TO %s" % mode - log_message(log_file, message, mode='info') - - act_map = anything_to_hdr_convert(act_map) - att_map = anything_to_hdr_convert(att_map) - - cambia_formato = rsc.get_rsc('change_format', 'fruitcake') - cambia_val = rsc.get_rsc('change_values', 'fruitcake') - - new_att_map = join(output_dir, "att_map_" + mode +".hdr") - new_act_map = join(output_dir,"act_map_" + mode +".hdr") - - if mode == "SimSET": - - #rcommand = '%s %s %s 0.096 1 >> %s' % (cambia_val, att_map, new_att_map, log_file) - rcommand = '%s %s %s 0.096 4 >> %s' % (cambia_val, att_map, new_att_map, log_file) - osrun(rcommand, log_file) - rcommand = '%s %s %s 0.135 3 >> %s' % (cambia_val, new_att_map, new_att_map, log_file) - osrun(rcommand, log_file) - rcommand = '%s %s %s 1B >> %s' % (cambia_formato, new_att_map, new_att_map, log_file) - osrun(rcommand, log_file) - rcommand = '%s %s %s 1B >> %s' % (cambia_formato, act_map, new_act_map, log_file) - osrun(rcommand, log_file) - - if mode == "STIR": - - rcommand = '%s %s %s fl >> %s' % (cambia_formato, att_map, new_att_map, log_file) - osrun(rcommand, log_file) - rcommand = '%s %s %s fl >> %s' % (cambia_formato, new_act_map, new_act_map, log_file) - osrun(rcommand, log_file) - rcommand = '%s %s %s 1 0.096 >> %s' % (cambia_val, new_att_map, new_att_map, log_file) - osrun(rcommand, log_file) - rcommand = '%s %s %s 3 0.135 >> %s' % (cambia_val, new_att_map, new_att_map, log_file) - osrun(rcommand, log_file) - - return new_act_map, new_att_map - -def ncounts(image_hdr): - - img, data = nib_load(image_hdr) - ncounts = np.sum(data) - return ncounts - -def convert_simset_sino_to_stir(input_img, output=False): - - ## To be continued.... - - simset_img = nib.load(input_img) - simset_img_data = simset_img.get_fdata() - shape = simset_img_data.shape - - n_slices = shape[2] - nrings = np.sqrt(n_slices) - n_x = shape[0] - - input_definition = [] - - for i in range(n_slices): - - ring1,ring2 = divmod(i,nrings) - segment = ring1 - ring2 - slice_def = [i, ring1, ring2, segment] - input_definition.append(slice_def) - - output_definition = sorted(input_definition, key=itemgetter(3)) - stir_img_data = np.empty(shape, dtype=float, order='C') - stir_img_data_flip_x = np.empty(shape, dtype=float, order='C') - - for i in range(n_slices): - - output_index = output_definition[i][0] - input_slice = simset_img_data[:,:,output_index] - stir_img_data[:,:,i] = input_slice - - for j in range(n_x): - stir_img_data_flip_x[j,:,:] = stir_img_data[n_x-1-j, :, :] - - #stir_img = nib.AnalyzeImage(stir_img_data, simset_img.affine, simset_img.header) - stir_img = nib.AnalyzeImage(stir_img_data_flip_x, simset_img.affine, simset_img.header) - - if not output: - output = input_img [0:-4] + '_stir.hdr' - - nib.save(stir_img,output) - - -def resampleXYvoxelSizes(image_hdr, xyVoxelSize, log_file): - img = nib.load(image_hdr) - z_VoxelSize =img.header['pixdim'][3] - taget_vsize = (xyVoxelSize,xyVoxelSize,z_VoxelSize) - res_img = nibp.resample_to_output(img, taget_vsize, out_class=img.__class__) - nib.save(res_img,image_hdr) - - return image_hdr - - -def resampleZvoxelSize(image_hdr, zOutputvoxelSize, log_file): - img = nib.load(image_hdr) - x_VoxelSize = img.header['pixdim'][1] - y_VoxelSize = img.header['pixdim'][2] - taget_vsize = (x_VoxelSize,y_VoxelSize,zOutputvoxelSize) - res_img = nibp.resample_to_output(img, taget_vsize, out_class=img.__class__) - nib.save(res_img,image_hdr) - - return image_hdr - - -def makeImageSquare(image_hdr, log_file): - - # Tools for image manipulation - corta_pega_filcol_hdr = rsc.get_rsc('corta_pega_filcol_hdr', 'fruitcake') - zpix, zsize, xpix, xsize, ypix, ysize = read_analyze_header(image_hdr,log_file) - - pixDif = abs(xpix-ypix) - - if pixDif != 0: - if pixDif % 2 == 0: - oneSide = pixDif/2 - otherSide = oneSide - else: - oneSide = np.trunc(pixDif/2) - otherSide = oneSide +1 - - if xpix > ypix: - rcommand1 = '%s %s %s peg %s %s >> %s' % (corta_pega_filcol_hdr, image_hdr, image_hdr, "1", str(oneSide), log_file) - rcommand2 = '%s %s %s peg %s %s >> %s' % (corta_pega_filcol_hdr, image_hdr, image_hdr, "2", str(otherSide), log_file) - elif ypix > xpix: - rcommand1 = '%s %s %s peg %s %s >> %s' % (corta_pega_filcol_hdr, image_hdr, image_hdr, "3", str(oneSide), log_file) - rcommand2 = '%s %s %s peg %s %s >> %s' % (corta_pega_filcol_hdr, image_hdr, image_hdr, "4", str(otherSide), log_file) - - osrun(rcommand1, log_file) - osrun(rcommand2, log_file) - - -def scalImage(image_hdr, maxValue, log_file): - """ - Given the input image_hdr (analyze format), this function scales their values to the input maxValue - - """ - cambia_val_interval = rsc.get_rsc('change_interval', 'fruitcake') - img, data =nib_load(image_hdr, log_file) - - factor = maxValue/np.max(data) - - mask_hdr="mask.hdr" - rcommand = '%s %s %s 0 0.9 0 >> %s' % (cambia_val_interval, image_hdr, mask_hdr, log_file) - osrun(rcommand, log_file) - rcommand = '%s %s %s 0.9 10000000000 1 >> %s' % (cambia_val_interval, mask_hdr, mask_hdr, log_file) - osrun(rcommand, log_file) - - operate_single_image(image_hdr, "mult", factor, image_hdr, log_file) - - rcommand = '%s %s %s 0 1 1 >> %s' % (cambia_val_interval, image_hdr, image_hdr, log_file) - osrun(rcommand, log_file) - - operate_images_analyze(image_hdr, mask_hdr, image_hdr, "mult") - - os.remove(mask_hdr) - os.remove("mask.img") - -def remove_neg_nan(image_hdr): - - img, data = nib_load(image_hdr) - indx = np.where(data<0) - data[indx] = 0 - indx = np.where(np.isnan(data)) - data[indx] = 0 - - imageToWrite = nib.AnalyzeImage(data,img.affine,img.header) - nib.save(imageToWrite, "aux.hdr") - copy_analyze("aux.hdr",image_hdr) - os.remove("aux.hdr") - os.remove("aux.img") - - #return image_hdr - -def update_act_map(spmrun, act_map, att_map, orig_pet, simu_pet, output_act_map, axialFOV, log_file): - - output_dir = dirname(output_act_map) - # mfileFusion = join(dirname(output_act_map),"fusion.m") - - # Getting the necesary resources - # cambia_formato = rsc.get_rsc('change_format', 'fruitcake') - # cambia_val_interval = rsc.get_rsc('change_interval', 'fruitcake') - # operate_image_hdr = rsc.get_rsc('operate_image', 'fruitcake') - # remove_nan_hdr = rsc.get_rsc('erase_nans', 'fruitcake') - # remove_neg_hdr = rsc.get_rsc('erase_negs', 'fruitcake') - - # First step is coregistering the output image with the orig pet (coregistered to the mri) - #act_map_img = act_map[0:-3]+"img" - # simu_pet_img = simu_pet[0:-3]+"img" - orig_pet_img = orig_pet[0:-3]+"img" - # coreg_simpet_img = spm.image_fusion(spmrun, mfileFusion, orig_pet_img, simu_pet_img, log_file) - coreg_simpet_hdr = fsl_flirt(orig_pet, simu_pet, log_file) - coreg_simpet_img =coreg_simpet_hdr[0:-3]+"img" - - act, act_data = nib_load(act_map) - - - # Next, we will do a scaling by the mean - norm_factor = proportional_scaling(coreg_simpet_hdr, orig_pet, orig_pet, log_file) - operate_single_image(coreg_simpet_hdr,'mult',norm_factor, coreg_simpet_hdr, log_file) - - # Now we do a smoothing of both data to avoid multiply noise and perform the division - mfileSmooth = join(dirname(output_act_map),"smooth.m") - division_hdr = join(output_dir, "division.hdr") - s_coreg_simpet_img = spm.smoothing(spmrun, mfileSmooth, coreg_simpet_img, 5, "s", log_file) - s_orig_pet_img = spm.smoothing(spmrun, mfileSmooth, orig_pet_img, 5, "s", log_file) - - s_coreg_simpet_hdr = s_coreg_simpet_img[0:-3]+"hdr" - s_orig_pet_hdr = s_orig_pet_img[0:-3]+"hdr" - operate_images_analyze(s_orig_pet_hdr, s_coreg_simpet_hdr, division_hdr, 'div') - # rcommand = '%s %s %s %s fl divid' % (operate_image_hdr, s_orig_pet_hdr, s_coreg_simpet_hdr, division_hdr) - # osrun(rcommand, log_file) - - # Now we do some stuff on the division image to avoid problems - pet_mask_hdr = join(dirname(orig_pet), "pet_mask.hdr") - deleteValuesOutFov(pet_mask_hdr, axialFOV/2, act.shape[2]/2) - - # rcommand = '%s %s %s %s fl multi' % (operate_image_hdr, division_hdr, pet_mask_hdr, division_hdr) - # osrun(rcommand, log_file) - # rcommand = '%s %s %s >> %s' % (remove_nan_hdr, division_hdr, division_hdr, log_file) - # osrun(rcommand, log_file) - # rcommand = '%s %s %s >> %s' % (remove_neg_hdr, division_hdr, division_hdr, log_file) - # osrun(rcommand, log_file) - fix_4d_image(pet_mask_hdr) - operate_images_analyze(division_hdr, pet_mask_hdr, division_hdr, "mult") - remove_neg_nan(division_hdr) - change_interval_values(division_hdr, division_hdr, 5, 100000000000000000000000000000000000,1) - # rcommand = '%s %s %s 5 100000000000000000000000000000000000 1' % (cambia_val_interval, division_hdr, division_hdr) - # osrun(rcommand, log_file) - - #operate_images_analyze(division_hdr, pet_mask_hdr, division_hdr, 'mult') - operate_images_analyze(division_hdr, pet_mask_hdr, division_hdr, "mult") - # rcommand = '%s %s %s %s fl multi' % (operate_image_hdr, division_hdr, pet_mask_hdr, division_hdr) - # osrun(rcommand, log_file) - - #and finally, we calculate the new activity image - fix_4d_image(act_map) - operate_images_analyze(division_hdr, act_map, output_act_map, "mult") - # rcommand = '%s %s %s %s fl multi' % (operate_image_hdr, act_map, division_hdr, output_act_map) - # osrun(rcommand, log_file) - scalImage(output_act_map, 127, log_file) - change_format(output_act_map,"1B", log_file) - # rcommand = '%s %s %s 1B >> %s' % (cambia_formato, output_act_map, output_act_map, log_file) - # osrun(rcommand, log_file) - -def fsl_flirt(reference_hdr, input_hdr, log_file): - components = os.path.split(input_hdr) - coreg_hdr = os.path.join(components[0], 'r' + components[1]) - flt = fsl.FLIRT(bins=640, cost_func='mutualinfo') - flt.inputs.in_file = input_hdr - flt.inputs.reference = reference_hdr - flt.out_file = coreg_hdr - flt.out_log=log_file - flt.save_log=True - flt.inputs.output_type = "NIFTI_GZ" - flt.cmdline - 'flirt -in %s -ref %s -out %s -dof 6 -cost mutualinfo -searchrx -30 30 -searchry -30 30 -searchrz -30 30' %(input_hdr, reference_hdr, coreg_hdr) - aux = os.getcwd() - os.chdir(components[0]) - flt.run() - os.chdir(aux) - - anything_to_hdr_convert(input_hdr[0:-4]+"_flirt.nii.gz") - copy_analyze(input_hdr[0:-4]+"_flirt.hdr", coreg_hdr) - os.system("rm %s*" % input_hdr[0:-4]+"_flirt.*") - - return coreg_hdr - - -def proportional_scaling(img,ref_img,mask_img, log_file): - - img_max, img_mean = compute_vmax_vmean(img, mask_img) - ref_max, ref_mean = compute_vmax_vmean(ref_img, mask_img) - - if float(img_mean) != 0: - fnorm = ref_mean / img_mean - return fnorm - else: - message = 'Error scaling image. Image maean value is zero: ' + str(img) - log_message(log_file, message, 'error') - #print message - -def compute_vmax_vmean(img, mask_img): - """ - Compute maximum and mean intensity values on input image counting - on voxels inside the reference image (brain mask) - :param img: (string) input image path - :param ref_img: (string) reference image path - :return: - """ - - img, data = nib_load(img) - data = np.nan_to_num(data) - - ref_img, data_ref = nib_load(mask_img) - data_ref = np.nan_to_num(data_ref) - - i_max = np.amax(data_ref) - - super_threshold_indices = data_ref > 0.2*i_max - data_ref[super_threshold_indices] = 0 - - # Compute values restricted to voxels inside mask (ref image) - indx = np.where((data>0) & (data_ref.reshape(data.shape)>0)) - - # Maximum intensity value - v_max = np.max(data[indx]) - # Mean intensity value - v_mean = np.mean(data[indx]) - - return v_max, v_mean - - -def deleteValuesOutFov(mask_hdr, max_z, central_slice): - - img, data = nib_load(mask_hdr) - data = np.nan_to_num(data) - - z = img.shape[2] - z_size = img.header['pixdim'][3] #mm - max_z = int(round((float(max_z)*10)/float(z_size))) - - i_min = int(central_slice)-max_z - i_max = int(central_slice)+max_z - - for i in range(0, i_min): - data[:,:,i]=0 - - for i in range(i_max,z): - data[:,:,i]=0 - - img_to_write = nib.Nifti1Image(data, img.affine, img.header) - nib.save(img_to_write,mask_hdr) - -def compute_corr_coeff(img1, img2, log_file): - """ - Compute the correlation coefficiente between img1 and img2 - :param img1: (string) input header image1 path - :param img2: (string) input header image2 path - :return: - """ - - img1_img, img1_data = nib_load(img1, log_file) - img2_img, img2_data = nib_load(img2, log_file) - corrCoefmtx = np.corrcoef(img1_data.flatten(),img2_data.flatten()) - - return corrCoefmtx[0,1] - -def fix_4d_image(image_hdr): - img, data = nib_load(image_hdr) - shape = data.shape - - if len(shape) != 3: - data_new = data[:,:,:,0] - - imageToWrite = nib.AnalyzeImage(data_new,img.affine,img.header) - nib.save(imageToWrite, "aux.hdr") - copy_analyze("aux.hdr",image_hdr) - os.remove("aux.hdr") - os.remove("aux.img") - -def change_interval_values(input_hdr, output_hdr, min_value, max_value, rep_value): - img, data = nib_load(input_hdr) - - - indx = np.where(datamin_value) - - data[indx] = rep_value - - imageToWrite = nib.AnalyzeImage(data,img.affine,img.header) - nib.save(imageToWrite, "aux.hdr") - copy_analyze("aux.hdr",output_hdr) - os.remove("aux.hdr") - os.remove("aux.img") - -def change_format(image_hdr, newFormat, logfile): - message="" - img, data = nib_load(image_hdr) - - if newFormat == "fl": - form=np.float32 - elif newFormat == "1B": - form=np.uint8 - else: - message = "Error! Invalid format (or not implemented yet): " +str(newFormat) - print(message) - log_message(logfile, message, 'error') - - if message=="": - data = data.astype(form) - hdr = nib.AnalyzeHeader() - hdr.set_data_dtype(form) - hdr.set_data_shape(data.shape) - imageToWrite = nib.AnalyzeImage(data, img.affine, hdr) - nib.save(imageToWrite, image_hdr) - -def fix_4d_data(data): - shape = data.shape - - if len(shape) == 3: - return data - else: - return data[:, :, :, 0] +import os, shutil, datetime +import re +from os.path import join, exists, isdir, dirname, basename, splitext +from subprocess import getstatusoutput as getoutput +import nibabel as nib +import nibabel.processing as nibp +from utils import resources as rsc +from utils import spm_tools as spm +import numpy as np +from operator import itemgetter +from nipype.interfaces.dcm2nii import Dcm2nii +from nipype.interfaces import fsl + + +def osrun(command, logfile, catch_out=False): + """ + Executes command, raise an error if it fails and send status to logger + :param command: command to be executed + :param logger: logger file + :return: + """ + if catch_out: + status, out = getoutput(command) + if status != 0: + log_message(logfile, command, 'error') + raise TypeError(command) + else: + log_message(logfile, command) + return out + + else: + if os.system(command) != 0: + + log_message(logfile, command, 'error') + raise TypeError(command) + else: + log_message(logfile, command) + +def nib_load(image, logfile=False): + """ + Load image and return data array + :param image: image to load data from + :return: image data array + """ + try: + img = nib.load(image) + data = img.get_data()[:,:,:] + return img, data + except Exception as e: + message = "Error: " + str(e) + if logfile: + log_message(logfile, message, 'error') + else: + print(message) + +def copy_analyze(image1, image2=False, dest_dir=False, logfile=False): #NOTE: Potentially unused. + """ + Create a copy of an Analyze format image + :param image1: (string) path to the original image + :param image2: (string, optional) path to the copy image + :param dest_dir: (string, optional) path to the destination folder + :return: + """ + + if image2: + if image1[-3:] == 'hdr' or image1[-3:] == 'img' and image2[-3:] == 'hdr' or image2[-3:] == 'img': + image1_hdr = image1[0:-3] + 'hdr' + image2_hdr = image2[0:-3] + 'hdr' + shutil.copy(image1_hdr, image2_hdr) + image1_img = image1[0:-3] + 'img' + image2_img = image2[0:-3] + 'img' + shutil.copy(image1_img, image2_img) + + return image2_hdr + else: + message = 'Error! The provided image is not in Analyze format:' + str(image1) + if logfile: log_message(logfile, message, 'error') + raise TypeError(message) + + elif not image2 and isdir(str(dest_dir)): + + ext = splitext(basename(image1))[1] + + if ext == '.img' or ext == '.hdr': + image1_hdr = image1[0:-3] + 'hdr' + image2_hdr = join(dest_dir, basename(image1)[0:-3] + 'hdr') + shutil.copy(image1_hdr, image2_hdr) + image1_img = image1[0:-3] + 'img' + image2_img = join(dest_dir, basename(image1)[0:-3] + 'img') + shutil.copy(image1_img, image2_img) + + return image2_hdr + else: + message = 'Error! The provided image is not in Analyze format:' + str(image1) + if logfile: log_message(logfile, message, 'error') + raise TypeError(message) + else: + message = 'Error! Not image2 path or dest dir path provided: ' + str(image2) + ', ' + str(dest_dir) + if logfile: log_message(logfile, message, 'error') + raise TypeError(message) + + +def copy_nifti(image1, image2=False, dest_dir=False, logfile=False): + """ + Create a copy of an Analyze format image + :param image1: (string) path to the original image + :param image2: (string, optional) path to the copy image + :param dest_dir: (string, optional) path to the destination folder + :return: + """ + + if image2: + if image1[-3:] == 'nii' and image2[-3:] == 'nii': + image1_img = image1[0:-3] + 'nii' + image2_img = image2[0:-3] + 'nii' + shutil.copy(image1_img, image2_img) + + return image2_img + + elif image1[-6:] == 'nii.gz' and image2[-6:] == 'nii.gz': + image1_img = image1[0:-6] + 'nii.gz' + image2_img = image2[0:-6] + 'nii.gz' + shutil.copy(image1_img, image2_img) + + return image2_img + + else: + message = 'Error! The provided image is not in Nifti format:' + str(image1) + if logfile: log_message(logfile, message, 'error') + raise TypeError(message) + + elif not image2 and isdir(str(dest_dir)): + + ext = splitext(basename(image1))[1] + + if ext == '.nii': + image1_img = image1[0:-3] + 'nii' + image2_img = join(dest_dir, basename(image1)[0:-3] + 'nii') + shutil.copy(image1_img, image2_img) + + return image2_img + + if ext == '.nii.gz': + image1_img = image1[0:-6] + 'nii.gz' + image2_img = join(dest_dir, basename(image1)[0:-6] + 'nii.gz') + shutil.copy(image1_img, image2_img) + + return image2_img + + else: + message = 'Error! The provided image is not in Nifti format:' + str(image1) + if logfile: log_message(logfile, message, 'error') + raise TypeError(message) + else: + message = 'Error! Not image2 path or dest dir path provided: ' + str(image2) + ', ' + str(dest_dir) + if logfile: log_message(logfile, message, 'error') + raise TypeError(message) + + + + +def create_analyze_from_imgdata(data, out, pix_x, pix_y, pix_z, tx, ty, tz, data_type="fl"): + + if data_type == "1b": + dtype=np.int8 + elif data_type == "2b": + dtype=np.int16 + elif data_type == "db": + dtype=np.float64 + else: + dtype=np.float32 + + hdr1 = nib.AnalyzeHeader() + hdr1.set_data_dtype(dtype) + hdr1.set_data_shape((pix_x,pix_y,pix_z)) + hdr1.set_zooms((tx,ty,tz)) + + f = open(data,'rb') + + img_data = hdr1.raw_data_from_fileobj(f) + + analyze_img = nib.AnalyzeImage(img_data, hdr1.get_base_affine(), hdr1) + + nib.save(analyze_img,out) + + +def create_nifti_from_imgdata(data, out, pix_x, pix_y, pix_z, tx, ty, tz, data_type="fl"): + + if data_type == "1b": + dtype=np.int8 + elif data_type == "2b": + dtype=np.int16 + elif data_type == "db": + dtype=np.float64 + else: + dtype=np.float32 + + hdr1 = nib.Nifti2Header() + #hdr1 = nib.Nifti1Header() + hdr1.set_data_dtype(dtype) + hdr1.set_data_shape((pix_x,pix_y,pix_z)) + hdr1.set_zooms((tx,ty,tz)) + + f = open(data,'rb') + + img_data = hdr1.raw_data_from_fileobj(f) + + nifti_img = nib.Nifti2Image(img_data, hdr1.get_base_affine(), hdr1) + #nifti_img = nib.Nifti1Image(img_data, hdr1.get_base_affine(), hdr1) + + nib.save(nifti_img,out) + + +def create_nifti1_from_imgdata(data, out, pix_x, pix_y, pix_z, tx, ty, tz, data_type="fl"): + + if data_type == "1b": + dtype=np.int8 + elif data_type == "2b": + dtype=np.int16 + elif data_type == "db": + dtype=np.float64 + else: + dtype=np.float32 + + hdr1 = nib.Nifti1Header() + hdr1.set_data_dtype(dtype) + hdr1.set_data_shape((pix_x,pix_y,pix_z)) + hdr1.set_zooms((tx,ty,tz)) + + f = open(data,'rb') + + img_data = hdr1.raw_data_from_fileobj(f) + + nifti_img = nib.Nifti1Image(img_data, hdr1.get_base_affine(), hdr1) + + nib.save(nifti_img,out) + + + +def read_analyze_header(header_file,logfile): + + img = nib.load(header_file) + + zpix = img.shape[2] + #zsize = abs(img.get_affine[2,2]) + zsize = abs(img.affine[2,2]) + xpix = img.shape[0] + #xsize = abs(img.get_affine[0,0]) + xsize = abs(img.affine[0,0]) + ypix = img.shape[1] + #ysize = abs(img.get_affine[1,1]) + ysize = abs(img.affine[1,1]) + + return zpix, zsize, xpix, xsize, ypix, ysize + +def write_interfile_header(header_file,matrix_size_x,pixel_size_x, matrix_size_y,pixel_size_y, matrix_size_z,pixel_size_z, offset_z=0): + + image_v = os.path.basename(header_file)[0:-2] + "v" + fheader_hv = open(header_file, "w") + + offset_x = matrix_size_x/2*pixel_size_x + 0.5*pixel_size_x + offset_y = matrix_size_y/2*pixel_size_y + 0.5*pixel_size_y + + fheader_hv.write("!INTERFILE :=\n") + fheader_hv.write("name of data file := %s\n" % image_v) + fheader_hv.write("!GENERAL DATA :=\n") + fheader_hv.write("!GENERAL IMAGE DATA :=\n") + fheader_hv.write("!type of data := PET\n") + fheader_hv.write("imagedata byte order := LITTLEENDIAN\n") + fheader_hv.write("!PET STUDY (General) :=\n") + fheader_hv.write("!PET data type := Image\n") + fheader_hv.write("process status := Reconstructed\n") + fheader_hv.write("!number format := float\n") + fheader_hv.write("!number of bytes per pixel := 4\n") + fheader_hv.write("number of dimensions := 3\n") + fheader_hv.write("matrix axis label [1] := x\n") + fheader_hv.write("!matrix size [1] := %s\n" % matrix_size_x) + fheader_hv.write("scaling factor (mm/pixel) [1] := %s\n" % pixel_size_x) + fheader_hv.write("matrix axis label [2] := y\n") + fheader_hv.write("!matrix size [2] := %s\n" % matrix_size_y) + fheader_hv.write("scaling factor (mm/pixel) [2] := %s\n" % pixel_size_y) + fheader_hv.write("matrix axis label [3] := z\n") + fheader_hv.write("!matrix size [3] := %s\n" % matrix_size_z) + fheader_hv.write("scaling factor (mm/pixel) [3] := %s\n" % pixel_size_z) + fheader_hv.write("first pixel offset (mm) [1] := -%s\n" % offset_x) + fheader_hv.write("first pixel offset (mm) [2] := -%s\n" % offset_y) + fheader_hv.write("first pixel offset (mm) [3] := 0\n") + fheader_hv.write("number of time frames := 1\n") + fheader_hv.write("!END OF INTERFILE :=\n") + + fheader_hv.close() + +def write_interfile_header_mu(header_file,matrix_size_x,pixel_size_x, matrix_size_y,pixel_size_y, matrix_size_z,pixel_size_z, offset_z=None): + + image_v = os.path.basename(header_file)[0:-2] + "v" + fheader_hv = open(header_file, "w") + + offset_x = matrix_size_x/2*pixel_size_x + 0.5*pixel_size_x + offset_y = matrix_size_y/2*pixel_size_y + 0.5*pixel_size_y + + fheader_hv.write("!INTERFILE :=\n") + fheader_hv.write("name of data file := %s\n" % image_v) + fheader_hv.write("!GENERAL DATA :=\n") + fheader_hv.write("!GENERAL IMAGE DATA :=\n") + fheader_hv.write("!type of data := PET\n") + fheader_hv.write("imagedata byte order := LITTLEENDIAN\n") + fheader_hv.write("!PET STUDY (General) :=\n") + fheader_hv.write("!PET data type := Image\n") + fheader_hv.write("process status := Reconstructed\n") + fheader_hv.write("!number format := float\n") + fheader_hv.write("!number of bytes per pixel := 4\n") + fheader_hv.write("number of dimensions := 3\n") + fheader_hv.write("matrix axis label [1] := x\n") + fheader_hv.write("!matrix size [1] := %s\n" % matrix_size_x) + fheader_hv.write("scaling factor (mm/pixel) [1] := %s\n" % pixel_size_x) + fheader_hv.write("matrix axis label [2] := y\n") + fheader_hv.write("!matrix size [2] := %s\n" % matrix_size_y) + fheader_hv.write("scaling factor (mm/pixel) [2] := %s\n" % pixel_size_y) + fheader_hv.write("matrix axis label [3] := z\n") + fheader_hv.write("!matrix size [3] := %s\n" % matrix_size_z) + fheader_hv.write("scaling factor (mm/pixel) [3] := %s\n" % pixel_size_z) + + if offset_z != None: + fheader_hv.write("first pixel offset (mm) [3] := %s\n" % offset_z) + + fheader_hv.write("number of time frames := 1\n") + fheader_hv.write("!END OF INTERFILE :=\n") + + fheader_hv.close() + +def nii_analyze_convert(image, logfile=False, outfile=False): + """ + Converts the provided image file to the analyze or nifti format + depending on the filex extension + :param image: image to be converted + :return: + """ + + # if not logfile: logfile = join(dirname(image), 'log_nii_analyze_convert.txt') + + ext = os.path.splitext(image)[1] + img, data = nib_load(image) + data = data.astype(np.float32) #casting to avoid problems when applying nib.save with no float data + hdr = nib.AnalyzeHeader() + hdr.set_data_dtype(np.float32) + hdr.set_data_shape(data.shape) + + if ext == '.nii': + # n2a = rsc.get_rsc('nii2analyze', 'exe') + if outfile: + cimage = outfile + else: + cimage = image.replace('.nii', '.hdr') + # rcommand = '%s %s %s >> %s' % (n2a, image, cimage, logfile) + # osrun(rcommand, logfile) + imageToWrite = nib.AnalyzeImage(data, img.affine, hdr) + elif ext == '.hdr': + # a2n = rsc.get_rsc('analyze2nii', 'exe') + if outfile: + cimage = outfile + else: + cimage = image.replace('.hdr', '.nii') + # rcommand = '%s %s %s >> %s' % (a2n, image, cimage, logfile) + # osrun(rcommand, logfile) + imageToWrite = nib.Nifti1Image(data, img.affine, hdr) + else: + # a2n = rsc.get_rsc('analyze2nii', 'exe') + if outfile: + cimage = outfile + else: + cimage = image.replace('.img', '.nii') + # image_hdr = str(image).replace(".img", ".hdr") + # rcommand = '%s %s %s >> %s' % (a2n, image_hdr, cimage, logfile) + # osrun(rcommand, logfile) + imageToWrite = nib.Nifti1Image(data, img.affine, hdr) + + nib.save(imageToWrite, cimage) + + if exists(cimage): + if outfile: + shutil.move(cimage, outfile) + return outfile + else: + return cimage + else: + return 0 + +def anything_to_hdr_convert(image, logfile=False, outfile=False ): + """ + This function will try to convert dicom, .nii.gz or .nii to analyze... + """ + + if image[-3:] == "hdr" or image[-3:] == "img": + return image[0:-3]+"hdr" + + elif image[-3:] == "nii": + hdr = nii_analyze_convert(image,logfile=logfile) + if exists(hdr): + return hdr + else: + raise TypeError ("nifti-analyze conversion failed....") + + elif image[-6:]=="nii.gz": + hdr = image[0:-6] + 'hdr' + nii_analyze_convert(image,logfile=logfile, outfile=hdr) + if exists(hdr): + return hdr + else: + raise TypeError ("nifti-analyze conversion failed....") + + elif isdir(image): + print("I think the provided image:\n %s \n is a dicom directory. I will try to convert it...") + dicom2nii = rsc.get_rsc("dicom2nii","exe") + nii = join(dirname(image),"image.nii") + rcommand = '%s %s %s' % (dicom2nii, image, nii) + osrun(rcommand,logfile) + if exists(nii): + hdr = nii_analyze_convert(nii,logfile=logfile) + os.remove(nii) + if exists(hdr): + return hdr + else: + raise TypeError ("nifti-analyze conversion failed....") + else: + raise TypeError ("dicom-nifti conversion failed....") + + elif image[-2:] == "hv": #Only works with floats data + + with open(image) as f: + lines = f.readlines() + + lines = [x.strip() for x in lines] + + pixel = [re.compile(f"\!matrix size \[{i}\].*") for i in range(1, 4)] + pixel_size = [re.compile(f"scaling factor \(mm.pixel\) \[{i}\].*") for i in range(1, 4)] + + pixel_x = list(filter(pixel[0].match, lines))[0].split()[-1] + pixel_size_x = list(filter(pixel_size[0].match, lines))[0].split()[-1] + pixel_y = list(filter(pixel[1].match, lines))[0].split()[-1] + pixel_size_y = list(filter(pixel_size[1].match, lines))[0].split()[-1] + pixel_z = list(filter(pixel[2].match, lines))[0].split()[-1] + pixel_size_z = list(filter(pixel_size[2].match, lines))[0].split()[-1] + + hdr_header = image[0:-2] + "hdr" + # img_file = image[0:-2] + "img" + data_file = image[0:-2] + "v" + + #This will convert .hv to .hdr and copy the data + # gen_hdr = rsc.get_rsc("gen_hdr", "fruitcake") + # rcommand = "%s %s %s %s %s fl %s %s %s 0" % (gen_hdr, hdr_header[0:-4], pixel_x, pixel_y, pixel_z, pixel_size_x, pixel_size_y, pixel_size_z) + # osrun(rcommand, logfile) + # shutil.copy(data_file, img_file) + create_analyze_from_imgdata(data_file,hdr_header,float(pixel_x),float(pixel_y),float(pixel_z),float(pixel_size_x),float(pixel_size_y),float(pixel_size_z)) + #nii_analyze_convert("inicial.hdr",logfile,"aux.nii") + #nii_analyze_convert("aux.nii",logfile,hdr_header) + #os.remove("aux.nii") + + if exists(hdr_header): + return hdr_header + else: + raise TypeError ("Interfile-analyze conversion failed....") + +def convert_hv_to_nii(image, logfile=False, outfile=False): + + with open(image) as f: + lines = f.readlines() + + lines = [x.strip() for x in lines] + + pixel = [re.compile(f"\!matrix size \[{i}\].*") for i in range(1, 4)] + pixel_size = [re.compile(f"scaling factor \(mm.pixel\) \[{i}\].*") for i in range(1, 4)] + + pixel_x = list(filter(pixel[0].match, lines))[0].split()[-1] + pixel_size_x = list(filter(pixel_size[0].match, lines))[0].split()[-1] + pixel_y = list(filter(pixel[1].match, lines))[0].split()[-1] + pixel_size_y = list(filter(pixel_size[1].match, lines))[0].split()[-1] + pixel_z = list(filter(pixel[2].match, lines))[0].split()[-1] + pixel_size_z = list(filter(pixel_size[2].match, lines))[0].split()[-1] + + hdr_header = image[0:-2] + "nii" + data_file = image[0:-2] + "v" + + create_nifti1_from_imgdata(data_file,hdr_header,float(pixel_x),float(pixel_y),float(pixel_z),float(pixel_size_x),float(pixel_size_y),float(pixel_size_z)) + + if exists(hdr_header): + return hdr_header + else: + raise TypeError ("Interfile-analyze conversion failed....") + + + + +def prepare_input_image(image_hdr, logfile, min_voxel_size=1): #TODO: This may be optimised. + """ + This method converts input_image to float data type, re-sizes the image to + 1mm size voxels if too large to keep a reasonable analysis execution time and + removes the negative and NaN values. + :param input_image: image to prepare for the analysis + :return: + """ + + # Tools for image manipulation + # change_format = rsc.get_rsc('change_format', 'fruitcake') + # change_matrix = rsc.get_rsc('change_img_matrix', 'fruitcake') + # erase_negs = rsc.get_rsc('erase_negs', 'fruitcake') + # erase_nans = rsc.get_rsc('erase_nans', 'fruitcake') + + # Convert input image to Analyze format and float data type + # rcommand = '%s %s %s fl >> %s' % (change_format, image_hdr, image_hdr, logfile) + # osrun(rcommand, logfile) + + ## # CONTEMPLAR POSIBILIDAD DE CONVERTIR ESTO EN FUNCIÓN PARA SUSTITUIR AL "cambia_formato_hdr" de fruitcake + change_format(image_hdr, "fl", logfile) + + ### + + # Resize image if necessary + # ndims = recalculate_matrix(image_hdr, min_voxel_size) + # rcommand = '%s %s %s %s %s %s novecino >> %s' % (change_matrix, image_hdr, image_hdr, + # str(ndims[0]), str(ndims[1]), str(ndims[2]), logfile) + # osrun(rcommand, logfile) + image_hdr = resampleXYvoxelSizes(image_hdr, min_voxel_size, logfile) + image_hdr = resampleZvoxelSize(image_hdr, min_voxel_size, logfile) + + # Erase negative and NaN values + # rcommand = '%s %s %s >> %s' % (erase_negs, image_hdr, image_hdr, logfile) + # osrun(rcommand, logfile) + # rcommand = '%s %s %s >> %s' % (erase_nans, image_hdr, image_hdr, logfile) + # osrun(rcommand, logfile) + remove_neg_nan(image_hdr) + + + return image_hdr + +def recalculate_matrix(input_image, voxelsize, mode="downsampling"): + + new_dimensions = [] + image_load, data = nib_load(input_image) + header = image_load.header + sizes = header.get_zooms() + dimensions = image_load.shape + x_lenght = int(sizes[0]*dimensions[0]/voxelsize) + y_lenght = int(sizes[1]*dimensions[1]/voxelsize) + z_lenght = int(sizes[2]*dimensions[2]/voxelsize) + if mode == "downsampling": + if sizes[0]voxelsize: + new_dimensions.append(x_lenght) + else: + new_dimensions.append(dimensions[0]) + if sizes[1]>voxelsize: + new_dimensions.append(y_lenght) + else: + new_dimensions.append(dimensions[1]) + if sizes[2]>voxelsize: + new_dimensions.append(z_lenght) + else: + new_dimensions.append(dimensions[2]) + + return new_dimensions + +def verify_roi_exists(rois_image, roi_number): + """ + Compute the number of voxels in each ROI + :param rois_image: (string) path to the ROI parcelled image + :param rois_list: (array) ROIs indexing + :return: False if the number of voxels is zero + :return: nvox is not zero + """ + rois_img = nib.load(rois_image) + rois_data = rois_img.get_data()[:, :, :] + indx = np.where(rois_data == roi_number) + nvox = len(rois_data[indx]) + if nvox == 0: + return False + else: + return True + +def operate_single_image(input_image, operation, factor, output_image, logfile): #NOTE: Potentially unused. + """ + Given an input image, multiply or divide it by a numerical factor + saving the result as output_image + :param input_image: image base operation on + #:param operation: 1 = multiply, 2 = divide + :param operation: 'mult' = multiply, 'div' = divide + :param factor: operation factor + :param output_image: output image file + :return: + """ + + img = nib.load(input_image) + data = img.get_data()[:,:,:] + data = np.nan_to_num(data) + + + if operation == 'mult': + data = data * float(factor) + elif operation == 'div': + data = data / float(factor) + else: + message = "Error! Invalid operation: " +str(operation) + print(message) + log_message(logfile, message, 'error') + + hdr1 = nib.AnalyzeHeader() + hdr1.set_data_dtype(img.get_data_dtype()) + hdr1.set_data_shape(img.shape) + #hdr1.set_zooms(abs(np.diag(img.affine))) + #hdr1.set_zooms(abs(np.diag(img.affine))[0:3]) + #hdr1.set_zooms(abs(np.diag(img.affine))[0:4]) + hdr1.set_zooms(abs(np.diag(img.affine))[0:img.ndim]) + + analyze_img = nib.AnalyzeImage(data, hdr1.get_base_affine(), hdr1) + + nib.save(analyze_img,output_image) + + +def operate_single_image_nii(input_image, operation, factor, output_image, logfile, check_nans=True): + """ + Given an input image, multiply or divide it by a numerical factor + saving the result as output_image + :param input_image: image base operation on + #:param operation: 1 = multiply, 2 = divide + :param operation: 'mult' = multiply, 'div' = divide + :param factor: operation factor + :param output_image: output image file + :return: + """ + + img = nib.load(input_image) + data = img.get_data()[:,:,:] + #data = np.nan_to_num(data) #re-added + """ + if check_nans: #TODO SEE WHAT HAPPENS TO THIS + data = np.nan_to_num(data) + #data[np.isnan(data)]=0 + """ + if operation == 'mult': + #data = data * float(factor) + data *= float(factor) #re-added + elif operation == 'div': + #data = data / float(factor) + data /= float(factor) + else: + message = "Error! Invalid operation: " +str(operation) + print(message) + log_message(logfile, message, 'error') + + hdr1 = nib.Nifti2Header() + #hdr1 = nib.Nifti1Header() + hdr1.set_data_dtype(img.get_data_dtype()) + hdr1.set_data_shape(img.shape) + #hdr1.set_zooms(abs(np.diag(img.affine))) + #hdr1.set_zooms(abs(np.diag(img.affine))[0:3]) + #hdr1.set_zooms(abs(np.diag(img.affine))[0:4]) + hdr1.set_zooms(abs(np.diag(img.affine))[0:img.ndim]) + + nifti_img = nib.Nifti2Image(data, hdr1.get_base_affine(), hdr1) + #nifti_img = nib.Nifti1Image(data, hdr1.get_base_affine(), hdr1) + + nib.save(nifti_img,output_image) + + +def operate_images_analyze(image1, image2, out_image, operation='mult'): #NOTE: Potentially unused. + """ + Given the input images, calculate the multiplication image or the ratio between them + :param image1: string, path to the first image + :param image2: string, path to the second image + :param operation: string, multi (default) for multiplication divid for division + :param out_image: string (optional), path to the output image + :return: + """ + img1, data1 = nib_load(image1) + img2, data2 = nib_load(image2) + + # TODO CHECK IF NEGATIVE VALUES NEED TO BE REMOVED + # Remove NaN and negative values + #data1 = np.nan_to_num(data1) + #data2 = np.nan_to_num(data2) + + if operation == 'mult': + res_data = data1 * data2 + elif operation == 'div': + res_data = data1 / data2 + elif operation == 'sum': + res_data = data1 + data2 + elif operation == 'diff': + res_data = data1 - data2 + else: + message = 'Error! Unknown operation: ' + str(operation) + raise TypeError(message) + + hdr1 = nib.AnalyzeHeader() + hdr1.set_data_dtype(img1.get_data_dtype()) + hdr1.set_data_shape(img1.shape) + ##hdr1.set_zooms(abs(np.diag(img1.affine))[0:3]) + ##hdr1.set_zooms(abs(np.diag(img1.affine))[0:4]) + hdr1.set_zooms(abs(np.diag(img1.affine))[0:img1.ndim]) + + analyze_img = nib.AnalyzeImage(res_data, hdr1.get_base_affine(), hdr1) + + nib.save(analyze_img,out_image) + + +def operate_images_nii(image1, image2, out_image, operation='mult', check_nans=True): + """ + Given the input images, calculate the multiplication image or the ratio between them + :param image1: string, path to the first image + :param image2: string, path to the second image + :param operation: string, multi (default) for multiplication divid for division + :param out_image: string (optional), path to the output image + :return: + """ + img1, data1 = nib_load(image1) + img2, data2 = nib_load(image2) + + #data1 = np.nan_to_num(data1) #re-added + #data2 = np.nan_to_num(data2) #re-added + + """ + if check_nans: #TODO: Implement this as may be needed. At the moment no NaNs are checked. + # Remove NaN and negative values + #data1 = np.nan_to_num(data1) + #data2 = np.nan_to_num(data2) + + data1[np.isnan(data1)]=0 + data2[np.isnan(data2)]=0 + """ + + if operation == 'mult': + #res_data = data1 * data2 #re-added + data1 *= data2 + elif operation == 'div': + #res_data = data1 / data2 + data1 /= data2 + elif operation == 'sum': + #res_data = data1 + data2 #re-added + data1 += data2 + elif operation == 'diff': + #res_data = data1 - data2 + data1 -= data2 + else: + message = 'Error! Unknown operation: ' + str(operation) + raise TypeError(message) + + hdr1 = nib.Nifti2Header() + hdr1.set_data_dtype(img1.get_data_dtype()) + hdr1.set_data_shape(img1.shape) + ##hdr1.set_zooms(abs(np.diag(img1.affine))[0:3]) + ##hdr1.set_zooms(abs(np.diag(img1.affine))[0:4]) + hdr1.set_zooms(abs(np.diag(img1.affine))[0:img1.ndim]) + + #nifti_img = nib.Nifti2Image(res_data, hdr1.get_base_affine(), hdr1) + nifti_img = nib.Nifti2Image(data1, hdr1.get_base_affine(), hdr1) + + nib.save(nifti_img,out_image) + + +def operate_sinograms_nii(image1, image2, out_image, operation='mult'): + """ + Given the input images, calculate the multiplication image or the ratio between them + :param image1: string, path to the first image + :param image2: string, path to the second image + :param operation: string, multi (default) for multiplication divid for division + :param out_image: string (optional), path to the output image + :return: + """ + img1, data1 = nib_load(image1) + img2, data2 = nib_load(image2) + + if operation == 'mult': + res_data = data1 * data2 + elif operation == 'div': + res_data = data1 / data2 + elif operation == 'sum': + res_data = data1 + data2 + elif operation == 'diff': + res_data = data1 - data2 + else: + message = 'Error! Unknown operation: ' + str(operation) + raise TypeError(message) + + hdr1 = nib.Nifti2Header() + hdr1.set_data_dtype(img1.get_data_dtype()) + hdr1.set_data_shape(img1.shape) + hdr1.set_zooms(abs(np.diag(img1.affine))[0:img1.ndim]) + + nifti_img = nib.Nifti2Image(res_data, hdr1.get_base_affine(), hdr1) + + nib.save(nifti_img,out_image) + + + +def smooth_analyze(image,fwhm, output): + + from nibabel import processing as nibproc + + img =nib.load(image) + smoothed=nibproc.smooth_image(img,fwhm,out_class=nib.AnalyzeImage) + nib.save(smoothed,output) + + +def smooth_nifti(image,fwhm, output): + + from nibabel import processing as nibproc + + img =nib.load(image) + smoothed=nibproc.smooth_image(img,fwhm,out_class=nib.Nifti2Image) + #smoothed=nibproc.smooth_image(img,fwhm,out_class=nib.Nifti1Image) + nib.save(smoothed,output) + + +def log_message(logfile, message, mode='info'): + """" + Print outs logger messages to the specified logfile + """ + + #ctime = read_time(datetime.datetime.now()) + ctime = datetime.datetime.now() + + separator = '\n###########################################################' + + if mode == 'exe': + stream = separator + '\nTIME: ' + str(ctime) + '\n\nEXE: ' + message + separator + '\n' + elif mode == 'info': + stream = separator + '\nTIME: ' + str(ctime) + '\n\nINFO: ' + message + separator + '\n' + elif mode == 'warning': + stream = separator + '\nTIME: ' + str(ctime) + '\n\nWARNING: ' + message + separator + '\n' + elif mode == 'error': + stream = separator + '\nTIME: ' + str(ctime) + '\n\nERROR: ' + message + separator + '\n' + else: + stream = separator + '\nTIME: ' + str(ctime) + '\n\n' + message + '\n' + + if exists(logfile): + with open(logfile, 'a') as lfile: + lfile.write(stream) + else: + with open(logfile, 'w') as lfile: + lfile.write(stream) + +def reorient_dcmtonii(image_path): + + + components=os.path.split(image_path) + path = components[0] + image = components[1] + converter = Dcm2nii() + converter.inputs.source_names=[image_path] + converter.inputs.reorient_and_crop=True + # converter.nii_output=False + converter.run() + + reor_ima_name = "o"+image + if exists(join(path,reor_ima_name)): + shutil.copy(join(path,reor_ima_name),image_path) + os.remove(join(path,reor_ima_name)) + os.remove(join(path,"co"+image)) + else: + # converter.inputs_reorient_and_crop=False + # converter.run() + # 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. + Required inputs are: + pet_image: dicom_dir, nii.gz, nii or Analyze (.hdr or .img) + mri_image: dicom_dir, nii.gz, nii or Analyze (.hdr or .img) + mode: Choose STIR or SIMSET. The maps will be different for each simulation. + The inputs will be stored to Data/simulation_name/Patient as reformatted/corregistered Analyze + #The maps will be stored on Data/simulation_name/Maps + The maps will be stored on Results/output_name/Maps + """ + message = "GENERATING ACT AND ATT MAPS FROM PET, (CT) and MR IMAGES" + log_message(log_file, message, mode='info') + + reorient_dcmtonii(pet_image) + reorient_dcmtonii(mri_image) + reorient_dcmtonii(ct_image) + #First of all lets take all to analyze + pet_hdr = anything_to_hdr_convert(pet_image) + #pet_hdr = copy_analyze(pet_hdr,image2=False,dest_dir=patient_dir) + pet_hdr = prepare_input_image(pet_hdr,log_file,min_voxel_size=1.5) + pet_img = pet_hdr[0:-3]+"img" + + if ct_image: #ct_image is not empty + ct_hdr = anything_to_hdr_convert(ct_image) + ct_hdr = prepare_input_image(ct_hdr,log_file,min_voxel_size=1.5) + ct_img = ct_hdr[0:-3]+"img" + + mri_hdr = anything_to_hdr_convert(mri_image) + #mri_hdr = copy_analyze(mri_hdr,image2=False,dest_dir=patient_dir) + mri_hdr = prepare_input_image(mri_hdr,log_file,min_voxel_size=1.5) + mri_img = mri_hdr[0:-3]+"img" + + #make mri image square + makeImageSquare(mri_hdr, log_file) + + #Performing PET-CT/MR coregister + mfile = os.path.join(output_dir,"fusion_pet_to_mri.m") + # correg_pet_hdr = fsl_flirt(mri_hdr, pet_hdr, log_file) + # correg_pet_img =correg_pet_hdr[0:-3]+"img" + correg_pet_img = spm.image_fusion(spm_run, mfile, mri_img, pet_img, log_file) + correg_ct_img="" + if ct_image: #ct_image is not empty + mfile = os.path.join(output_dir,"fusion_ct_to_mri.m") + correg_ct_img = spm.image_fusion(spm_run, mfile, mri_img, ct_img, log_file) + + #Now the map generation + from utils.patient2maps import patient2maps + + my_map_generation = patient2maps(spm_run, output_dir, log_file, + mri_img, correg_pet_img, correg_ct_img, mode=mode) + activity_map_hdr, attenuation_map_hdr = my_map_generation.run() + + return activity_map_hdr, attenuation_map_hdr + +def convert_map_values(act_map,att_map,output_dir,log_file,mode="SimSET"): + + message = "CONVERTING ACT AND ATT TO %s" % mode + log_message(log_file, message, mode='info') + + act_map = anything_to_hdr_convert(act_map) + att_map = anything_to_hdr_convert(att_map) + + cambia_formato = rsc.get_rsc('change_format', 'fruitcake') + cambia_val = rsc.get_rsc('change_values', 'fruitcake') + + new_att_map = join(output_dir, "att_map_" + mode +".hdr") + new_act_map = join(output_dir,"act_map_" + mode +".hdr") + + if mode == "SimSET": + + #rcommand = '%s %s %s 0.096 1 >> %s' % (cambia_val, att_map, new_att_map, log_file) + rcommand = '%s %s %s 0.096 4 >> %s' % (cambia_val, att_map, new_att_map, log_file) + osrun(rcommand, log_file) + rcommand = '%s %s %s 0.135 3 >> %s' % (cambia_val, new_att_map, new_att_map, log_file) + osrun(rcommand, log_file) + rcommand = '%s %s %s 1B >> %s' % (cambia_formato, new_att_map, new_att_map, log_file) + osrun(rcommand, log_file) + rcommand = '%s %s %s 1B >> %s' % (cambia_formato, act_map, new_act_map, log_file) + osrun(rcommand, log_file) + + if mode == "STIR": + + rcommand = '%s %s %s fl >> %s' % (cambia_formato, att_map, new_att_map, log_file) + osrun(rcommand, log_file) + rcommand = '%s %s %s fl >> %s' % (cambia_formato, new_act_map, new_act_map, log_file) + osrun(rcommand, log_file) + rcommand = '%s %s %s 1 0.096 >> %s' % (cambia_val, new_att_map, new_att_map, log_file) + osrun(rcommand, log_file) + rcommand = '%s %s %s 3 0.135 >> %s' % (cambia_val, new_att_map, new_att_map, log_file) + osrun(rcommand, log_file) + + return new_act_map, new_att_map + +def ncounts(image_hdr): + + img, data = nib_load(image_hdr) + ncounts = np.sum(data) + return ncounts + +def convert_simset_sino_to_stir(input_img, output=False): #NOTE: Potentially unused. + + ## To be continued.... + + simset_img = nib.load(input_img) + simset_img_data = simset_img.get_fdata() + shape = simset_img_data.shape + + n_slices = shape[2] + nrings = np.sqrt(n_slices) + n_x = shape[0] + + input_definition = [] + + for i in range(n_slices): + + ring1,ring2 = divmod(i,nrings) + segment = ring1 - ring2 + slice_def = [i, ring1, ring2, segment] + input_definition.append(slice_def) + + output_definition = sorted(input_definition, key=itemgetter(3)) + stir_img_data = np.empty(shape, dtype=float, order='C') + stir_img_data_flip_x = np.empty(shape, dtype=float, order='C') + + for i in range(n_slices): + + output_index = output_definition[i][0] + input_slice = simset_img_data[:,:,output_index] + stir_img_data[:,:,i] = input_slice + + for j in range(n_x): + stir_img_data_flip_x[j,:,:] = stir_img_data[n_x-1-j, :, :] + + #stir_img = nib.AnalyzeImage(stir_img_data, simset_img.affine, simset_img.header) + stir_img = nib.AnalyzeImage(stir_img_data_flip_x, simset_img.affine, simset_img.header) + + if not output: + output = input_img [0:-4] + '_stir.hdr' + + (stir_img_data_flip_x.astype(np.float32)).tofile(output[0:-4] + '.bin') #ADDED, DELETE. + + nib.save(stir_img,output) + + +def convert_simset_sino_to_stir_nii(input_img, output=False): + + simset_img = nib.load(input_img) + simset_img_data = simset_img.get_fdata(dtype=np.float32) + shape = simset_img_data.shape + + n_slices = shape[2] + nrings = np.sqrt(n_slices) + n_x = shape[0] + + input_definition = [] + + for i in range(n_slices): + + ring1,ring2 = divmod(i,nrings) + segment = ring1 - ring2 + slice_def = [i, ring1, ring2, segment] + input_definition.append(slice_def) + + output_definition = sorted(input_definition, key=itemgetter(3)) + stir_img_data = np.empty(shape, dtype=np.float32, order='C') + + for i in range(n_slices): + + output_index = output_definition[i][0] + input_slice = simset_img_data[:,:,output_index] + stir_img_data[:,:,i] = input_slice + + + #TODO: THIS CAN BE SURELY SPEED-UP (FOR TOTAL-BODY), CHECK THIS: + """ + n_slices = shape[2] + nrings = int(np.sqrt(n_slices)) # assuming perfect square + n_x = shape[0] + + # Build definition table + indices = np.arange(n_slices) + ring1 = indices // nrings + ring2 = indices % nrings + segment = ring1 - ring2 + + input_definition = np.column_stack([indices, ring1, ring2, segment]) + + # Sort by segment + #output_definition = input_definition[np.argsort(input_definition[:, 3])] #not working. + output_definition = sorted(input_definition, key=itemgetter(3)) + output_definition = np.array(output_definition) + + # Reorder slices in one vectorized step + order = output_definition[:, 0].astype(int) + stir_img_data = simset_img_data[:, :, order] + """ + + + #WARNING: This is still being investigated. At the moment we suspect SimSET sinograms are inverted in x and y respect to STIR. May be also inverted in z. + + #stir_img_data_flip = stir_img_data[::-1,:,:] + stir_img_data_flip = stir_img_data[::-1,::-1,:] #WARNING: THIS IS TEMPORAL + + stir_img = nib.Nifti2Image(stir_img_data_flip, simset_img.affine, simset_img.header) + + if not output: + output = input_img [0:-4] + '_stir.nii' + + nib.save(stir_img,output) + + + +def copy_sinogram_stir_to_output(input_img, output_img): + + + simset_sino = nib.load(input_img) + + if input_img[-7:] == '.nii.gz': + nib.save(simset_sino, input_img[0:-7] + '.hdr') + shutil.copy(input_img[0:-6] + "img", output_img) + os.remove(input_img[0:-6] + "img") + os.remove(input_img[0:-6] + "hdr") + + elif input_img[-4:] == '.nii': + nib.save(simset_sino, input_img[0:-4] + '.hdr') + shutil.copy(input_img[0:-3] + "img", output_img) + os.remove(input_img[0:-3] + "img") + os.remove(input_img[0:-3] + "hdr") + + else: + print("Image is not in the NifTi format!") + + +#NOTE: THIS FUNCTION MAY BE DELETED IN THE FUTURE... +def copy_reduced_sinogram_stir_to_output(input_img, output_img, nrings, max_segment): #NOTE: This is unused for now, but is needed for FBP3D + + simset_sino = nib.load(input_img) + sino_dataobj = simset_sino.dataobj + + init_idx = int(np.sum([i for i in range(1, nrings-max_segment)])) #+1)])) + final_idx = nrings*nrings - init_idx + + sino_reduced = sino_dataobj[..., init_idx:final_idx].astype(sino_dataobj.dtype) + + sino_reduced_nii = nib.Nifti2Image(sino_reduced, simset_sino.affine) + nib.save(sino_reduced_nii, input_img[0:-4] + '.hdr') + + shutil.copy(input_img[0:-3] + "img", output_img) + os.remove(input_img[0:-3] + "img") + os.remove(input_img[0:-3] + "hdr") + + +#NOTE: Following functions may be moved to a BrainVISET tools separate file... + +def resampleXYvoxelSizes(image_hdr, xyVoxelSize, log_file): + img = nib.load(image_hdr) + z_VoxelSize =img.header['pixdim'][3] + taget_vsize = (xyVoxelSize,xyVoxelSize,z_VoxelSize) + res_img = nibp.resample_to_output(img, taget_vsize, out_class=img.__class__) + nib.save(res_img,image_hdr) + + return image_hdr + + +def resampleZvoxelSize(image_hdr, zOutputvoxelSize, log_file): + img = nib.load(image_hdr) + x_VoxelSize = img.header['pixdim'][1] + y_VoxelSize = img.header['pixdim'][2] + taget_vsize = (x_VoxelSize,y_VoxelSize,zOutputvoxelSize) + res_img = nibp.resample_to_output(img, taget_vsize, out_class=img.__class__) + nib.save(res_img,image_hdr) + + return image_hdr + + +def makeImageSquare(image_hdr, log_file): + + # Tools for image manipulation + corta_pega_filcol_hdr = rsc.get_rsc('corta_pega_filcol_hdr', 'fruitcake') + zpix, zsize, xpix, xsize, ypix, ysize = read_analyze_header(image_hdr,log_file) + + pixDif = abs(xpix-ypix) + + if pixDif != 0: + if pixDif % 2 == 0: + oneSide = pixDif/2 + otherSide = oneSide + else: + oneSide = np.trunc(pixDif/2) + otherSide = oneSide +1 + + if xpix > ypix: + rcommand1 = '%s %s %s peg %s %s >> %s' % (corta_pega_filcol_hdr, image_hdr, image_hdr, "1", str(oneSide), log_file) + rcommand2 = '%s %s %s peg %s %s >> %s' % (corta_pega_filcol_hdr, image_hdr, image_hdr, "2", str(otherSide), log_file) + elif ypix > xpix: + rcommand1 = '%s %s %s peg %s %s >> %s' % (corta_pega_filcol_hdr, image_hdr, image_hdr, "3", str(oneSide), log_file) + rcommand2 = '%s %s %s peg %s %s >> %s' % (corta_pega_filcol_hdr, image_hdr, image_hdr, "4", str(otherSide), log_file) + + osrun(rcommand1, log_file) + osrun(rcommand2, log_file) + + +def scalImage(image_hdr, maxValue, log_file): + """ + Given the input image_hdr (analyze format), this function scales their values to the input maxValue + + """ + cambia_val_interval = rsc.get_rsc('change_interval', 'fruitcake') + img, data =nib_load(image_hdr, log_file) + + factor = maxValue/np.max(data) + + mask_hdr="mask.hdr" + rcommand = '%s %s %s 0 0.9 0 >> %s' % (cambia_val_interval, image_hdr, mask_hdr, log_file) + osrun(rcommand, log_file) + rcommand = '%s %s %s 0.9 10000000000 1 >> %s' % (cambia_val_interval, mask_hdr, mask_hdr, log_file) + osrun(rcommand, log_file) + + operate_single_image(image_hdr, "mult", factor, image_hdr, log_file) + + rcommand = '%s %s %s 0 1 1 >> %s' % (cambia_val_interval, image_hdr, image_hdr, log_file) + osrun(rcommand, log_file) + + operate_images_analyze(image_hdr, mask_hdr, image_hdr, "mult") + + os.remove(mask_hdr) + os.remove("mask.img") + +def remove_neg_nan(image_hdr): + + img, data = nib_load(image_hdr) + indx = np.where(data<0) + data[indx] = 0 + indx = np.where(np.isnan(data)) + data[indx] = 0 + + imageToWrite = nib.AnalyzeImage(data,img.affine,img.header) + nib.save(imageToWrite, "aux.hdr") + copy_analyze("aux.hdr",image_hdr) + os.remove("aux.hdr") + os.remove("aux.img") + + #return image_hdr + +def update_act_map(spmrun, act_map, att_map, orig_pet, simu_pet, output_act_map, axialFOV, log_file): + + output_dir = dirname(output_act_map) + # mfileFusion = join(dirname(output_act_map),"fusion.m") + + # Getting the necesary resources + # cambia_formato = rsc.get_rsc('change_format', 'fruitcake') + # cambia_val_interval = rsc.get_rsc('change_interval', 'fruitcake') + # operate_image_hdr = rsc.get_rsc('operate_image', 'fruitcake') + # remove_nan_hdr = rsc.get_rsc('erase_nans', 'fruitcake') + # remove_neg_hdr = rsc.get_rsc('erase_negs', 'fruitcake') + + # First step is coregistering the output image with the orig pet (coregistered to the mri) + #act_map_img = act_map[0:-3]+"img" + # simu_pet_img = simu_pet[0:-3]+"img" + orig_pet_img = orig_pet[0:-3]+"img" + # coreg_simpet_img = spm.image_fusion(spmrun, mfileFusion, orig_pet_img, simu_pet_img, log_file) + coreg_simpet_hdr = fsl_flirt(orig_pet, simu_pet, log_file) + coreg_simpet_img =coreg_simpet_hdr[0:-3]+"img" + + act, act_data = nib_load(act_map) + + + # Next, we will do a scaling by the mean + norm_factor = proportional_scaling(coreg_simpet_hdr, orig_pet, orig_pet, log_file) + operate_single_image(coreg_simpet_hdr,'mult',norm_factor, coreg_simpet_hdr, log_file) + + # Now we do a smoothing of both data to avoid multiply noise and perform the division + mfileSmooth = join(dirname(output_act_map),"smooth.m") + division_hdr = join(output_dir, "division.hdr") + s_coreg_simpet_img = spm.smoothing(spmrun, mfileSmooth, coreg_simpet_img, 5, "s", log_file) + s_orig_pet_img = spm.smoothing(spmrun, mfileSmooth, orig_pet_img, 5, "s", log_file) + + s_coreg_simpet_hdr = s_coreg_simpet_img[0:-3]+"hdr" + s_orig_pet_hdr = s_orig_pet_img[0:-3]+"hdr" + operate_images_analyze(s_orig_pet_hdr, s_coreg_simpet_hdr, division_hdr, 'div') + # rcommand = '%s %s %s %s fl divid' % (operate_image_hdr, s_orig_pet_hdr, s_coreg_simpet_hdr, division_hdr) + # osrun(rcommand, log_file) + + # Now we do some stuff on the division image to avoid problems + pet_mask_hdr = join(dirname(orig_pet), "pet_mask.hdr") + deleteValuesOutFov(pet_mask_hdr, axialFOV/2, act.shape[2]/2) + + # rcommand = '%s %s %s %s fl multi' % (operate_image_hdr, division_hdr, pet_mask_hdr, division_hdr) + # osrun(rcommand, log_file) + # rcommand = '%s %s %s >> %s' % (remove_nan_hdr, division_hdr, division_hdr, log_file) + # osrun(rcommand, log_file) + # rcommand = '%s %s %s >> %s' % (remove_neg_hdr, division_hdr, division_hdr, log_file) + # osrun(rcommand, log_file) + fix_4d_image(pet_mask_hdr) + operate_images_analyze(division_hdr, pet_mask_hdr, division_hdr, "mult") + remove_neg_nan(division_hdr) + change_interval_values(division_hdr, division_hdr, 5, 100000000000000000000000000000000000,1) + # rcommand = '%s %s %s 5 100000000000000000000000000000000000 1' % (cambia_val_interval, division_hdr, division_hdr) + # osrun(rcommand, log_file) + + #operate_images_analyze(division_hdr, pet_mask_hdr, division_hdr, 'mult') + operate_images_analyze(division_hdr, pet_mask_hdr, division_hdr, "mult") + # rcommand = '%s %s %s %s fl multi' % (operate_image_hdr, division_hdr, pet_mask_hdr, division_hdr) + # osrun(rcommand, log_file) + + #and finally, we calculate the new activity image + fix_4d_image(act_map) + operate_images_analyze(division_hdr, act_map, output_act_map, "mult") + # rcommand = '%s %s %s %s fl multi' % (operate_image_hdr, act_map, division_hdr, output_act_map) + # osrun(rcommand, log_file) + scalImage(output_act_map, 127, log_file) + change_format(output_act_map,"1B", log_file) + # rcommand = '%s %s %s 1B >> %s' % (cambia_formato, output_act_map, output_act_map, log_file) + # osrun(rcommand, log_file) + +def fsl_flirt(reference_hdr, input_hdr, log_file): + components = os.path.split(input_hdr) + coreg_hdr = os.path.join(components[0], 'r' + components[1]) + flt = fsl.FLIRT(bins=640, cost_func='mutualinfo') + flt.inputs.in_file = input_hdr + flt.inputs.reference = reference_hdr + flt.out_file = coreg_hdr + flt.out_log=log_file + flt.save_log=True + flt.inputs.output_type = "NIFTI_GZ" + flt.cmdline + 'flirt -in %s -ref %s -out %s -dof 6 -cost mutualinfo -searchrx -30 30 -searchry -30 30 -searchrz -30 30' %(input_hdr, reference_hdr, coreg_hdr) + aux = os.getcwd() + os.chdir(components[0]) + flt.run() + os.chdir(aux) + + anything_to_hdr_convert(input_hdr[0:-4]+"_flirt.nii.gz") + copy_analyze(input_hdr[0:-4]+"_flirt.hdr", coreg_hdr) + os.system("rm %s*" % input_hdr[0:-4]+"_flirt.*") + + return coreg_hdr + + +def proportional_scaling(img,ref_img,mask_img, log_file): + + img_max, img_mean = compute_vmax_vmean(img, mask_img) + ref_max, ref_mean = compute_vmax_vmean(ref_img, mask_img) + + if float(img_mean) != 0: + fnorm = ref_mean / img_mean + return fnorm + else: + message = 'Error scaling image. Image maean value is zero: ' + str(img) + log_message(log_file, message, 'error') + #print message + +def compute_vmax_vmean(img, mask_img): + """ + Compute maximum and mean intensity values on input image counting + on voxels inside the reference image (brain mask) + :param img: (string) input image path + :param ref_img: (string) reference image path + :return: + """ + + img, data = nib_load(img) + data = np.nan_to_num(data) + + ref_img, data_ref = nib_load(mask_img) + data_ref = np.nan_to_num(data_ref) + + i_max = np.amax(data_ref) + + super_threshold_indices = data_ref > 0.2*i_max + data_ref[super_threshold_indices] = 0 + + # Compute values restricted to voxels inside mask (ref image) + indx = np.where((data>0) & (data_ref.reshape(data.shape)>0)) + + # Maximum intensity value + v_max = np.max(data[indx]) + # Mean intensity value + v_mean = np.mean(data[indx]) + + return v_max, v_mean + + +def deleteValuesOutFov(mask_hdr, max_z, central_slice): + + img, data = nib_load(mask_hdr) + data = np.nan_to_num(data) + + z = img.shape[2] + z_size = img.header['pixdim'][3] #mm + max_z = int(round((float(max_z)*10)/float(z_size))) + + i_min = int(central_slice)-max_z + i_max = int(central_slice)+max_z + + for i in range(0, i_min): + data[:,:,i]=0 + + for i in range(i_max,z): + data[:,:,i]=0 + + img_to_write = nib.Nifti1Image(data, img.affine, img.header) + nib.save(img_to_write,mask_hdr) + +def compute_corr_coeff(img1, img2, log_file): + """ + Compute the correlation coefficiente between img1 and img2 + :param img1: (string) input header image1 path + :param img2: (string) input header image2 path + :return: + """ + + img1_img, img1_data = nib_load(img1, log_file) + img2_img, img2_data = nib_load(img2, log_file) + corrCoefmtx = np.corrcoef(img1_data.flatten(),img2_data.flatten()) + + return corrCoefmtx[0,1] + +def fix_4d_image(image_hdr): + img, data = nib_load(image_hdr) + shape = data.shape + + if len(shape) != 3: + data_new = data[:,:,:,0] + + imageToWrite = nib.AnalyzeImage(data_new,img.affine,img.header) + nib.save(imageToWrite, "aux.hdr") + copy_analyze("aux.hdr",image_hdr) + os.remove("aux.hdr") + os.remove("aux.img") + +def change_interval_values(input_hdr, output_hdr, min_value, max_value, rep_value): + img, data = nib_load(input_hdr) + + + indx = np.where(datamin_value) + + data[indx] = rep_value + + imageToWrite = nib.AnalyzeImage(data,img.affine,img.header) + nib.save(imageToWrite, "aux.hdr") + copy_analyze("aux.hdr",output_hdr) + os.remove("aux.hdr") + os.remove("aux.img") + +def change_format(image_hdr, newFormat, logfile): + message="" + img, data = nib_load(image_hdr) + + if newFormat == "fl": + form=np.float32 + elif newFormat == "1B": + form=np.uint8 + else: + message = "Error! Invalid format (or not implemented yet): " +str(newFormat) + print(message) + log_message(logfile, message, 'error') + + if message=="": + data = data.astype(form) + hdr = nib.AnalyzeHeader() + hdr.set_data_dtype(form) + hdr.set_data_shape(data.shape) + imageToWrite = nib.AnalyzeImage(data, img.affine, hdr) + nib.save(imageToWrite, image_hdr) + +def fix_4d_data(data): + shape = data.shape + + if len(shape) == 3: + return data + else: + return data[:, :, :, 0] + + +def mu_coef_511keV(tissue_n): #NOTE: Check if this can be done with a internal function of SimSET. + rows_per_tissue = 1000 + + start = tissue_n * (rows_per_tissue + 1) + end = start + rows_per_tissue + + PHG_TABLE_PATH = "include/SimSET/2.9.2/phg.data/phg_att_table" + data = np.loadtxt(PHG_TABLE_PATH, skiprows=2+start, max_rows=rows_per_tissue) + + mu = data[510,1] + + return mu + + +def write_fwdproj_parfile(output_path): #TODO: ADD NRAY ARGUMENT! ORIGINAL WAS 1, POTENTIALLY 5. + content = """Forward Projector parameters:= + type := Matrix + Forward projector Using Matrix Parameters := + Matrix type := Ray Tracing + Ray tracing matrix parameters := + number of rays in tangential direction to trace for each bin := 5 + End Ray tracing matrix parameters := + End Forward Projector Using Matrix Parameters := +End:= +""" + + with open(output_path, "w", encoding="utf-8") as f: + f.write(content) + +def flip_rec_nifti(input_path, output_path=None): + + img = nib.load(input_path) + data = np.asanyarray(img.dataobj) + + # Flip data along axis 0 + flipped_data = np.flip(data, axis=0) + + # Fix affine: NOTE at the moment seems to not be necessary to change affine. + affine = img.affine.copy() + #affine[0, :] *= -1 + #affine[:3, 3] += img.affine[:3, 0] * (data.shape[0] - 1) + + # Create new image + new_img = nib.Nifti1Image(flipped_data, affine, img.header) + + # Save if path provided + if output_path is None: + base, ext = os.path.splitext(input_path) + if ext == ".gz": # handle .nii.gz + base, ext2 = os.path.splitext(base) + ext = ext2 + ext + output_path = base + "_flipped" + ext + + nib.save(new_img, output_path) + + return output_path