PET/MR Simulation

.gitattributes

@@ -0,0 +1,4 @@
+src/xDynamicSimulation/cDynamicSimulation/Testing/tests_dynamicsimulation.h merge=ours
+src/xDynamicSimulation/cDynamicSimulation/Testing/tests_dynamicsimulation.cpp merge=ours
+src/xGadgetron/cGadgetron/gadgetron_data_containers.h merge=ours
+src/xGadgetron/cGadgetron/gadgetron_data_containers.cpp merge=ours

.gitignore

@@ -1,6 +1,11 @@
 ## Ignore Visual Studio temporary files, build results, and
 ## files generated by popular Visual Studio add-ons.
 
+# CMake
+*.cmake
+*.make
+*.include_cache
+*.internal
 # User-specific files
 *.suo
 *.user
@@ -20,6 +25,14 @@ bld/
 [Bb]in/
 [Oo]bj/
 
+# CMake
+CMakeCache.txt
+CMakeFiles*
+
+# Makefiles
+**/Makefile*
+
+
 # Roslyn cache directories
 *.ide/
 

CMakeLists.txt

@@ -65,7 +65,7 @@ else()
 endif()
 option(SIRF_INSTALL_DEPENDENCIES "Install dlls etc" ${DEFAULT_SIRF_INSTALL_DEPENDENCIES})
 ####### CMake path
-set (CMAKE_MODULE_PATH "${CMAKE_CURRENT_LIST_DIR}/cmake")
+set (CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH};${CMAKE_CURRENT_LIST_DIR}/cmake")
 
 # If OSX give the advanced option to use absolute paths for shared libraries
 if (APPLE)

README.md

@@ -43,4 +43,4 @@ what was found/built:
 [style-badge]: https://api.codacy.com/project/badge/Grade/392861b4085f4f438d12c41029f86b47
 [style-link]: https://www.codacy.com/gh/SyneRBI/SIRF?utm_source=github.com&utm_medium=referral&utm_content=SyneRBI/SIRF&utm_campaign=Badge_Grade
 [zenodo-badge]: https://zenodo.org/badge/DOI/10.5281/zenodo.2707911.svg
-[zenodo-link]: https://doi.org/10.5281/zenodo.2707911
+[zenodo-link]: https://doi.org/10.5281/zenodo.2707911

cmake/setup.py.cmake

@@ -22,6 +22,7 @@ if(BUILD_PYTHON)
   python_pkg_alias(pyiutilities "sirf.pyiutilities")
   python_pkg_alias(pReg "sirf.Reg")
   python_pkg_alias(pyreg "sirf.pyreg")
+  python_pkg_alias(pysimulation "sirf.pysimulation")
   # convert to python CSV tuple for setup.py configure_file
   string(REPLACE ";" "', '" PYTHON_SETUP_PKGS_CSV "${PYTHON_SETUP_PKGS}")
   set(PYTHON_SETUP_PKGS_CSV "'${PYTHON_SETUP_PKGS_CSV}'")
@@ -42,8 +43,13 @@ if(BUILD_PYTHON)
   if(PYTHONINTERP_FOUND)
     # python setup.py install
     if("${PYTHON_STRATEGY}" STREQUAL "SETUP_PY")
+      message("intalling setup.py")
       install(CODE "execute_process(COMMAND\n\
+<<<<<<< HEAD
+        \"${PYTHON_EXECUTABLE}\" setup.py install --prefix ${PYTHON_DEST}\n\
+=======
         \"${Python_EXECUTABLE}\" setup.py build install\n\
+>>>>>>> rpe-encoding
         WORKING_DIRECTORY \"${PYTHON_DEST}\")")
     endif()
   endif(PYTHONINTERP_FOUND)

examples/Python/Simulation/prepare_EPG_input.py

@@ -0,0 +1,26 @@
+import sys 
+sys.path.append("/home/sirfuser/devel/buildVM/sources/SIRF/src/xDynamicSimulation/pDynamicSimulation/")
+
+import numpy as np
+
+import TissueParameterList as TPL
+
+root_path = '/media/sf_CCPPETMR/TestData/Input/xDynamicSimulation/pDynamicSimulation/' 
+xml_path = root_path + 'Cube128/XCAT_TissueParameters_XML.xml'
+fpath_output = root_path + 'Fingerprints/XCAT_tissue_parameter_list.npz'
+
+
+# parse file
+tpl = TPL.TissueParameterList()
+tpl.parse_xml(xml_path)
+tpl.print_contents()
+
+# prepare EPG input array
+all_mr_params = tpl.mr_as_array()
+mr_params = all_mr_params[:,1:]
+mr_params[:,0] /= 100
+mr_params[:,-1] = 0
+
+mr_params_unique,idx_inverse = np.unique(mr_params, axis=0, return_inverse=True)
+
+np.savez(fpath_output, mr_params_full=all_mr_params, mr_parameters=mr_params_unique, unique_idx_inverse=idx_inverse)

examples/Python/Simulation/preprocess_speedup_testdata.py

@@ -0,0 +1,51 @@
+
+import numpy as np 
+from pathlib import Path
+import nibabel as nib
+
+
+
+num_slices = 10
+Nz = 128
+slab_start = Nz//2 -num_slices//2
+slab_end = Nz//2 + num_slices//2
+
+
+def preprocess_mvfs(rootpath, folder_pattern, prefix_output):
+
+    fpath_input = Path(rootpath + folder_pattern)
+    print("looking in {} ".format(fpath_input))
+    list_files = sorted(fpath_input.glob('mvf_*'))
+
+    for f in list_files:
+        print("loading {}".format(f))
+
+        mvf = nib.load(str(f))
+        print("The input motoinfield has size {}".format(mvf.shape))
+
+        mvf = mvf.slicer[:,:,slab_start:slab_end,:]
+        print("The output motoinfield has size {}".format(mvf.shape))
+
+        fname_output = str(prefix_output + folder_pattern + f.name)
+        print("Storing to {}".format(fname_output))
+        nib.save(mvf, fname_output)
+
+root_path = '/media/sf_CCPPETMR/TestData/Input/xDynamicSimulation/pDynamicSimulation/'
+
+fpath_input = root_path + 'Cube128/'
+fpath_output = root_path + 'Slab128/'
+
+foldername_resp = 'mvf_resp/'
+foldername_card = 'mvf_card/'
+
+preprocess_mvfs(fpath_input, foldername_resp, fpath_output)
+preprocess_mvfs(fpath_input, foldername_card, fpath_output)
+
+
+fname_segmentation = fpath_input + 'label_volume.nii'
+seg = nib.load(fname_segmentation)
+
+print("The input segmentation has size {}".format(seg.shape))
+
+seg = seg.slicer[:,:,slab_start:slab_end]
+nib.save(seg, fpath_output + 'label_volume.nii')

examples/Python/Simulation/simulate_3D_RPE.py

@@ -0,0 +1,185 @@
+'''
+bla bla
+
+Usage:
+cartesian_3D_simulation.py [--help | options]
+
+Options:
+--non-interactive           do not show plots
+'''
+
+## SyneRBI Synergistic Image Reconstruction Framework (SIRF).
+## Copyright 2015 - 2020 Rutherford Appleton Laboratory STFC.
+## Copyright 2015 - 2017 University College London.
+## Copyright 2015 - 2017 Physikalisch-Technische Bundesanstalt.
+##
+## This is software developed for the Collaborative Computational
+## Project in Synergistic Reconstruction for Biomedical Imaging (formerly CCP PETMR)
+## (http://www.ccpsynerbi.ac.uk/).
+##
+## Licensed under the Apache License, Version 2.0 (the "License");
+##   you may not use this file except in compliance with the License.
+##   You may obtain a copy of the License at
+##       http://www.apache.org/licenses/LICENSE-2.0
+##   Unless required by applicable law or agreed to in writing, software
+##   distributed under the License is distributed on an "AS IS" BASIS,
+##   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+##   See the License for the specific language governing permissions and
+##   limitations under the License.
+
+__version__ = '0.1.0'
+from docopt import docopt
+
+args = docopt(__doc__, version=__version__)
+
+from pUtilities import *
+import sirf.Reg as pReg
+import sirf.DynamicSimulation as pDS
+import sirf.Gadgetron as pMR
+
+# import engine module
+
+# process command-line options
+show_plot = not args['--non-interactive']
+
+from pathlib import Path
+import numpy as np
+import time
+
+
+def read_motionfields(fpath_prefix):
+	p = sorted( Path(fpath_prefix).glob('*.nii') )
+	files = [x for x in p if x.is_file()]
+
+	temp = []
+	for f in files:
+		print("Reading from {} ... ".format(f))
+		img = pReg.NiftiImageData3DDisplacement(str(f))
+		temp.append(img)
+
+	data = np.array(temp, dtype=object)
+	return data
+
+def set_motionfields_from_path(modyn, fpath_prefix):
+
+	assert_validity(modyn, pDS.MRMotionDynamic)
+	mvfs = read_motionfields(fpath_prefix)
+
+	for m in mvfs:
+		modyn.add_displacement_field(m)
+
+
+def get_normed_surrogate_signal(t0_s, tmax_s, Nt, f_Hz):
+
+	t_s = np.linspace(t0_s, tmax_s, Nt)
+	sig = 0.5 * (1 + np.sin( 2*np.pi*f_Hz*t_s))
+	return t_s, sig
+
+
+def main():
+
+    fpath_testdata_prefix = '/media/sf_CCPPETMR/TestData/'
+    input_fpath_prefix = fpath_testdata_prefix + 'Input/xDynamicSimulation/pDynamicSimulation/'
+    output_fpath_prefix = fpath_testdata_prefix + 'Output/xDynamicSimulation/pDynamicSimulation/'
+
+    fpath_xml = input_fpath_prefix + 'Cube128/XCAT_TissueParameters_XML.xml'
+
+    fpath_template_contrast_rawdata = input_fpath_prefix + 'Cube128/CV_nav_cart_128Cube_FLASH_T1_defaultorient.h5'
+    fpath_template_acquisition_rawdata = input_fpath_prefix + 'General/meas_MID00241_FID69145_Tho_T1_fast_ismrmrd_defaultorient.h5'
+
+    acquisition_ad = pMR.AcquisitionData(fpath_template_acquisition_rawdata)
+    acquisition_ad = pMR.set_grpe_trajectory(acquisition_ad)
+
+    # configure the simulation
+    contrast_ad = pMR.AcquisitionData(fpath_template_contrast_rawdata)
+    contrast_ad = pMR.preprocess_acquisition_data(contrast_ad)
+
+    labels = pReg.NiftiImageData3D( input_fpath_prefix + "Cube128/label_volume_rai.nii"	)
+    mrsim = pDS.MRDynamicSimulation(labels, fpath_xml)
+
+    mrsim.set_contrast_template_data(contrast_ad)
+    mrsim.set_acquisition_template_data(acquisition_ad)
+
+    offset_z_mm = 0
+    translation = np.array([0, 0, offset_z_mm])
+    euler_angles_deg = np.array([0,0,0])
+
+    offset_trafo = pReg.AffineTransformation(translation, euler_angles_deg)
+    mrsim.set_offset_trafo(offset_trafo)
+
+    # take CSM from the rawdata itself
+    # could be replaced if independent way of computing CSM is available
+    csm = pMR.CoilSensitivityData()
+    csm.calculate(acquisition_ad)
+    mrsim.set_csm(csm)
+
+    # set which tissue defines SNR
+    SNR = 10
+    SNR_label = 13
+
+    mrsim.set_snr(SNR)
+    mrsim.set_snr_label(SNR_label)
+
+    # configure the surrogates
+    Nt = 10000
+    t0_s = 0
+    tmax_s = 60* 5 
+
+
+    f_Hz_resp = 0.2
+
+    t_resp, sig_resp = get_normed_surrogate_signal(t0_s, tmax_s, Nt, f_Hz_resp)
+
+
+    # configure the motion
+    num_motion_states = 4
+    # RESP
+    num_sim_resp_states = num_motion_states
+    # resp_motion = pDS.MRMotionDynamic( num_sim_resp_states )
+    # resp_motion.set_dynamic_signal(t_resp, sig_resp)
+    # resp_motion.set_cyclicality(False)
+    # resp_motion.set_groundtruth_folder_prefix(output_fpath_prefix + "output_simulation_3D_motiondata_r_{}_gt_resp".format(num_sim_resp_states))		
+    # set_motionfields_from_path(resp_motion, input_fpath_prefix + 'Cube128/mvf_resp/')
+    # mrsim.add_motion_dynamic(resp_motion)
+
+    #
+    fname_simulation_output = "output_simulation_3DRPE_motiondata_traj_r{}".format(num_sim_resp_states)
+
+    fname_output = output_fpath_prefix + fname_simulation_output + ".h5"
+    simulated_file = Path(fname_output)
+    if not simulated_file.is_file():
+
+        tstart = time.time()
+        mrsim.simulate_data()
+        print("--- Required {} minutes for the simulation.".format( (time.time()-tstart)/60))
+        mrsim.write_simulation_results(str(simulated_file))
+    else:
+        print("Skipping simulation since output file already exists.")
+
+    # mrsim.save_motion_ground_truth()
+
+    simulated_data = pMR.AcquisitionData(str(simulated_file))
+
+    csm.calculate(simulated_data)
+
+    AM = pMR.AcquisitionModel()
+    AM.set_coil_sensitivity_maps(csm)
+    AM.set_up(simulated_data, csm)
+
+    recon_img = AM.inverse(simulated_data)
+    recon_nii = pReg.NiftiImageData3D(recon_img)
+    recon_nii = recon_nii.abs()
+
+    fname_output = output_fpath_prefix + "recon_" + fname_simulation_output + ".nii"
+    recon_nii.write(fname_output)
+
+    return 1
+
+try:
+    main()
+    print('\n=== done with %s' % __file__)
+
+except error as err:
+    # display error information
+    print('??? %s' % err.value)
+    exit(1)