Merge pull request #1682 from NicoleJurjew/Vision_script_Updates

Vision script bugfix: rotational shift between STIR and Siemens recons

Author: KrisThielemans

examples/python/Vision_files_preprocess.py

@@ -64,6 +64,7 @@
 
 
 def process_files(prompts_header_filename,
+    file_prefix = '',
     mu_map_header = 'umap_00.h33', # the *.h33 header
     randoms_data_filename = 'smoothed_rand_00.s',
     scatter_2D_header_filename = 'scat_00_00.s.hdr',
@@ -98,31 +99,31 @@ def process_files(prompts_header_filename,
     # header-name for prompts as we don't want to overwrite the Siemens header
     prompts_header_to_read_withSTIR = prompts_header_filename[:-6] + '_readwSTIR.s.hdr'
     # STIR writes the (DOI-adapted) prompts out to a STIR-file:
-    prompts_filename_STIR_corr_DOI = 'prompts.hs'
+    prompts_filename_STIR_corr_DOI = file_prefix + 'prompts.hs'
     # STIR writes a non-TOF sinogram, name:
     nonTOF_template_sinogram_name = 'template_nonTOF.hs'
     # header-name for attenuation correction factors as we don't want to overwrite the Siemens header
     attenuation_corr_factor_header_to_read_withSTIR = attenuation_corr_factor_data_filename[:-2] + '_readwSTIR.s.hdr'
     # header-name for randoms as we don't want to overwrite the Siemens header
     norm_sino_to_read_withSTIR = norm_sino_data_filename[:-2] + '_readwSTIR.s.hdr'
     # STIR writes the DOI-adapted, negative corrected norm-sino to
-    norm_filename_fSTIR = 'norm_sino_fSTIR.hs'
+    norm_filename_fSTIR = file_prefix + 'norm_sino_fSTIR.hs'
     # STIR writes the DOI-adapted detection efficiencies to
-    det_effs_filename_fSIRF = 'detection_efficiencies_forSIRF.hs'
+    det_effs_filename_fSIRF = file_prefix + 'detection_efficiencies_forSIRF.hs'
     # header-name for randoms as we don't want to overwrite the Siemens header
     randoms_header_to_read_withSTIR = randoms_data_filename[:-2] + '_readwSTIR.s.hdr'
     # header-name for randoms as we don't want to overwrite the Siemens header
     scatter_2D_header_to_read_withSTIR = scatter_2D_header_filename[:-6] + '_readwSTIR.s.hdr'
     # STIR writes the (DOI-adapted) randoms to
-    randoms_adapted_DOI_filename = randoms_data_filename[:-2] + '_fSTIR.hs'
+    randoms_adapted_DOI_filename = file_prefix + randoms_data_filename[:-2] + '_fSTIR.hs'
     # STIR writes the (DOI-adapted), iSSRBd, unnormalized scatter to
-    scatter_3D_unnorm_filename = 'scatter_3D_unnormalized.hs'
+    scatter_3D_unnorm_filename = file_prefix + 'scatter_3D_unnormalized.hs'
     # STIR writes the additive term (that's normalized scatter + normalized randoms, attenuation corrected) to:
-    additive_term_filename_fSTIR = 'additive_term.hs'
+    additive_term_filename_fSTIR = file_prefix + 'additive_term.hs'
     # STIR writes the multiplicative term (that's norm_sino * attenuation_CORRECTION_factors) to:
-    multi_term_filename_fSTIR = 'mult_factors_forSTIR.hs'
+    multi_term_filename_fSTIR = file_prefix + 'mult_factors_forSTIR.hs'
     # STIR writes the multiplicative term for SIRF (that's detection_efficiency_sino * attenuation_factors) to:
-    multi_term_filename_fSIRF = 'mult_factors_forSIRF.hs'
+    multi_term_filename_fSIRF = file_prefix + 'mult_factors_forSIRF.hs'
 
     #%%
     try:
@@ -148,6 +149,10 @@ def process_files(prompts_header_filename,
     ## the crystal depth of interaction (DOI) from 7mm to 10mm to minimize the differences.
     if apply_DOI_adaption: DOI_adaption(prompts_from_e7, 10)
 
+    ###################### View Offset ADAPTION ############################
+
+    view_offset_adaption(prompts_from_e7, 0.0324)
+
     ## write to file so you can use it later
     prompts_from_e7.write_to_file(os.path.join(STIR_output_folder, prompts_filename_STIR_corr_DOI))
 
@@ -182,13 +187,24 @@ def process_files(prompts_header_filename,
 
     #### now let's read it from file and plot to see if it worked
     mu_map = stir.FloatVoxelsOnCartesianGrid.read_from_file(mu_map_header[:-3]+'hv')
+    mu_map.write_to_file(os.path.join(STIR_output_folder, file_prefix + 'mu_map.hv'))
     mu_map_arr = stirextra.to_numpy(mu_map)
 
     plt.figure()
-    plot_2d_image([1,1,1],mu_map_arr[mu_map_arr.shape[0]//2,:,:],'mu-map')
-    plt.savefig(os.path.join(STIR_output_folder,'mu_map.png'), transparent=False, facecolor='w')
+    plot_2d_image([1,2,1],mu_map_arr[mu_map_arr.shape[0]//2,:,:],'mu-map', cmap='Greys_r')
+    plot_2d_image([1,2,2],mu_map_arr[:,mu_map_arr.shape[1]//2,:],'mu-map', cmap='Greys_r')
+    plt.savefig(os.path.join(STIR_output_folder, file_prefix + 'mu_map.png'), transparent=False, facecolor='w')
     plt.close()
 
+#%%
+    ###################### UNI 1 image ############################
+    ## to get a uniform 1 image for initialising, use mu-map
+    ## This is important to match Siemens Image dimensions!!
+
+    uni1 = mu_map.get_empty_copy()
+    uni1.fill(1)
+    uni1.write_to_file(os.path.join(STIR_output_folder, file_prefix + 'uni1_img.hv'))
+
 
     # %%
     ###################### NORMALIZATION ############################
@@ -472,6 +488,15 @@ def DOI_adaption(projdata, DOI_new):
     DOI = proj_info.get_scanner().get_average_depth_of_interaction()
     print('New Depth of interaction:', DOI)
 
+def view_offset_adaption(projdata, view_offset):
+    proj_info = projdata.get_proj_data_info()
+
+    VO = proj_info.get_scanner().get_intrinsic_azimuthal_tilt()
+    print('Current view offset (rad):', VO)
+    proj_info.get_scanner().set_intrinsic_azimuthal_tilt(view_offset)
+    VO = proj_info.get_scanner().get_average_depth_of_interaction()
+    print('New Depth of interaction(rad):', VO)
+
 def check_if_compressed(header_filename):
     with open(header_filename) as f:
         data = f.read()
@@ -609,6 +634,7 @@ def remove_tof_dimension(header_name_new, header_name):
                     prog='Vision_files_preprocess.py',
                     description='Converts e7tools sinogram files for the Vision into sinogram files that can be read by STIR')
     parser.add_argument('--prompts_filename_inclPath', required=True, help="The filename of the prompts file, including path")
+    parser.add_argument('--file_prefix', default='', help="A prefix added to all output files")
     parser.add_argument('--mu_map_header', default='umap_00.h33', help="The filename of the mu-map header")
     parser.add_argument('--randoms_data_filename', default='smoothed_rand_00.s', help="The filename of the randoms data")
     parser.add_argument('--scatter_2D_header_filename', default='scat_00_00.s.hdr', help="The filename of the 2D scatter header")
@@ -619,9 +645,10 @@ def remove_tof_dimension(header_name_new, header_name):
     args = parser.parse_args()
 
     process_files(args.prompts_filename_inclPath,
+          file_prefix=args.file_prefix,
           mu_map_header=args.mu_map_header,
           randoms_data_filename=args.randoms_data_filename,
           scatter_2D_header_filename=args.scatter_2D_header_filename,
           norm_sino_data_filename=args.norm_sino_data_filename,
           attenuation_corr_factor_data_filename=args.attenuation_corr_factor_data_filename,
-          STIR_output_folder=args.STIR_output_folder)
+          STIR_output_folder=args.STIR_output_folder)