###############################################################################

## Vision_files_preprocess.py

##

## This script converts

## - (binary) sinogram files, generated by the Siemens research SW e7tools,

## measured on the Siemens Biograph Vision

## - into (binary) sinogram files that can be read by STIR.

###############################################################################

## First, you need to pre-process your data with JSRecon. This will generate

## some batch files for histogramming and reconstruction. Run the ones for

## histogramming.

##

## Then, in any batch-file for reconstruction, add the following 2 lines, then

## execute

## set cmd= %cmd% -d ./Debug

## set cmd= %cmd% --os scatter_520_2D.mhdr

##

## This way, it will write smoothed randoms, a mu-map and a sinogram for

## normalization to file.

##

## The e7tools provide the prompts sinogram in a compressed file-format.

## STIR can't read that, so you'll have to uncompress it first:

## C:\Siemens\PET\bin.win64-VG80\intfcompr.exe ^

## -e path\to\compressed\sinogram\filename.mhdr ^

## --oe path\to\UNcompressed\sinogram\NEWfilename.mhdr

##

## Here's a list of files you need (from the e7tools). Make sure they are in

## the same folder.

## - prompts-sino_uncompr_00.s.hdr & .s

## - smoothed_randoms_00.h33 & .s (in the "Debug" folder)

## - umap_00.h33 & .i (in the "Debug" folder)

## note: we highly recommend using this mu-map, as it is "cut to the FOV"

## and is positioned at scanner center.

## - norm3d_00.h33 & .a (in the "Debug" folder)

## - scatter_520_2D.s.hdr & .s (wherever you set the path to in the batch file)

## - acf_00.h33 & .a (in the "Debug" folder); these are attenuation correction

## factors

##

## To be able to read Siemens IMAGE data with STIR, you might have to adapt the

## Interfile header. In the header:

##

## - change the word "image data" to "imagedata"

## - remove line "!image relative start time (sec)...""

## - remove line "!image duration (sec):=1200..."

#%%

#

#

# Author Nicole Jurjew

# Copyright (C) 2025 University College London

# This file is part of STIR.

#

# SPDX-License-Identifier: Apache-2.0

#

# See STIR/LICENSE.txt for details

#

#%%

import os, re, argparse

import numpy as np

import matplotlib.pyplot as plt

import stir

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',

norm_sino_data_filename= 'norm3d_00.a',

attenuation_corr_factor_data_filename = 'acf_00.a',

STIR_output_folder = 'processing'):

"""This script converts

- (binary) sinogram files, generated by the Siemens research SW e7tools,

measured on the Siemens Biograph Vision

- into (binary) sinogram files that can be read by STIR.

Args:

prompts_header_filename (str): filename of the prompts sinogram header file, including path

mu_map_header (str, optional): filename of mu-map header file, must be in the same folder as prompts file. Defaults to 'umap_00.h33'.

randoms_data_filename (str, optional): filename of randoms sinogram header file, must be in the same

folder as prompts file. Defaults to 'smoothed_rand_00.s'.

scatter_2D_header_filename (str, optional): filename of 2D scatter header file, must be in the same

folder as prompts file. Defaults to 'scat_00_00.s.hdr'.

norm_sino_data_filename (str, optional): filename of norm sinogram header file, must be in the same

folder as prompts file. Defaults to 'norm3d_00.a'.

attenuation_corr_factor_data_filename (str, optional): filename of attenuation correction factor header

file, must be in the same folder as prompts file. Defaults to 'acf_00.a'.

STIR_output_folder (str, optional): folder where files are written to. Defaults to 'processing', within prompts folder.

"""

data_folder_PATH = os.path.dirname(prompts_header_filename)

os.chdir(data_folder_PATH)

prompts_header_filename = os.path.basename(prompts_header_filename)

apply_DOI_adaption = False

####### file-names, adapt as needed #######

# 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 = 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 = file_prefix + 'norm_sino_fSTIR.hs'

# STIR writes the DOI-adapted detection efficiencies to

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 = file_prefix + randoms_data_filename[:-2] + '_fSTIR.hs'

# STIR writes the (DOI-adapted), iSSRBd, unnormalized scatter to

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 = file_prefix + 'additive_term.hs'

# STIR writes the multiplicative term (that's norm_sino * attenuation_CORRECTION_factors) to:

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 = file_prefix + 'mult_factors_forSIRF.hs'

#%%

try:

os.mkdir(STIR_output_folder)

except FileExistsError:

print("STIR output folder exists, files are overwritten")

#%%

###################### PROMPTS FILE ############################

##### first, we check if the prompts file is compressed

check_if_compressed(prompts_header_filename)

##### as the e7tools run on Windows, the data file-name needs to be changed

change_datafilename_in_interfile_header(prompts_header_to_read_withSTIR, prompts_header_filename ,prompts_header_filename[:-4])

# ## we're ready to read the prompts with STIR now

prompts_from_e7 = stir.ProjData.read_from_file(prompts_header_to_read_withSTIR)

prompts_from_e7 = stir.ProjDataInMemory(prompts_from_e7)

#%%

###################### DOI ADAPTION ############################

## after comparing e7tools and STIR forward projections, we've found out we have to change

## the crystal depth of interaction (DOI) from 7mm to 10mm to minimize the differences.

if apply_DOI_adaption: DOI_adaption(prompts_from_e7, 10)

## 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))

#%%

###################### NON-TOF TEMPLATE ############################

## you might need a non-TOF sinogram with the same geometry (for attenuation factors e.g.)

## so we're creatiing one here

proj_info = prompts_from_e7.get_proj_data_info()

nonTOF_proj_info = proj_info.create_non_tof_clone()

stir.ProjDataInterfile(prompts_from_e7.get_exam_info(), nonTOF_proj_info, os.path.join(STIR_output_folder,nonTOF_template_sinogram_name))

#%%

######## select display variables

central_slice = proj_info.get_num_axial_poss(0)//2

TOF_bin = proj_info.get_num_tof_poss()//2

# to draw line-profiles, we'll average over a few slices, specify how many:

thickness_half = 5

# %%

###################### MU-MAP ############################

## to read in the mu-map, we have to convert the Siemens *.h33 header to STIR *.hv header via a STIR -script

# note: here, a file .h33.tmp is created, I don't know why, but can be deleted

cmd_string_conv_intf_to_stir = 'convertSiemensInterfileToSTIR.sh \

{} {}'.format(mu_map_header, mu_map_header[:-3]+'hv')

try:

exit_status = os.system(cmd_string_conv_intf_to_stir)

if exit_status != 0:

raise Exception('Command failed with status: {}'.format(exit_status))

except Exception as e:

print('An error occurred: {}'.format(e))

#### 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 = mu_map.as_array()

plt.figure()

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 ############################

### we're using the prompts-header as a template for the norm-header

change_datafilename_in_interfile_header(norm_sino_to_read_withSTIR, prompts_header_filename, norm_sino_data_filename)

change_datatype_in_interfile_header(norm_sino_to_read_withSTIR, 'float', 4)

#%%

#### ready to read in norm-sino with STIR

norm_sino = stir.ProjData.read_from_file(norm_sino_to_read_withSTIR)

norm_sino = stir.ProjDataInMemory(norm_sino)

###################### DOI ADAPTION ############################

## after comparing e7tools and STIR forward projections, we've found out we have to change

## the crystal depth of interaction (DOI) from 7mm to 10mm to minimize the differences.

if apply_DOI_adaption: DOI_adaption(norm_sino, 10)

norm_sino_arr = norm_sino.as_array()

##### in case there were bad miniblocks during your measurement, the norm-file

##### might contain negative values. We'll set them to a very high value here, such

##### that the detection efficiencies (1/norm-value) will be 0 (numerically)

norm_sino_arr[norm_sino_arr<=0.] = 10^37

norm_sino.fill(norm_sino_arr.flat)

#### this is the data STIR needs in a bin-normalisation from projdata, so we'll write it out

norm_sino.write_to_file(os.path.join(STIR_output_folder,norm_filename_fSTIR))

#%%

for i in range(33):

plt.figure()

plot_2d_image([1,1,1],norm_sino_arr[i, central_slice,:,:],'norm TOF bin {}'.format(i))

if i == TOF_bin: plt.savefig(os.path.join(STIR_output_folder,norm_filename_fSTIR[:-3] + '.png'), transparent=False, facecolor='w')

plt.close()

#%%

#### to compute scatter and for use in SIRF, we'll need the detection efficiencies

#### which are computed as 1/norm-sino

det_eff = stir.ProjDataInMemory(norm_sino)

det_eff_arr = np.nan_to_num(np.divide(1, norm_sino_arr))

det_eff.fill(det_eff_arr.flat)

det_eff.write_to_file(os.path.join(STIR_output_folder,det_effs_filename_fSIRF))

#%%

for i in range(33):

plt.figure()

plot_2d_image([1,1,1],det_eff_arr[i, central_slice,:,:],'det. eff., TOF bin {}'.format(i))

if i == TOF_bin: plt.savefig(os.path.join(STIR_output_folder,det_effs_filename_fSIRF[:-3]+'.png'), transparent=False, facecolor='w')

plt.close()

#%%

###################### RANDOMS ############################

### below, we will use the prompts-header as a template for a header to read in Siemens randoms

change_datafilename_in_interfile_header(randoms_header_to_read_withSTIR, prompts_header_filename, randoms_data_filename)

### need to change the file type, because prompts are in unsigned short, randoms are in float

change_datatype_in_interfile_header(randoms_header_to_read_withSTIR, 'float', 4)

### The data type descriptions are true for the prompts, but not

### for the randoms, so we need to remove them (and some other info) from the header.

remove_scan_data_lines_from_interfile_header(randoms_header_to_read_withSTIR, randoms_header_to_read_withSTIR)

remove_IMGDATADESC_lines_from_interfile_header(randoms_header_to_read_withSTIR, randoms_header_to_read_withSTIR)

### The first set of sinograms in the randoms file is the "delayeds", so we need

### to tell STIR to ignore that by setting a data offset. The other offset fields,

### we can just remove, in line with the "data types" we removed above.

remove_data_offset(randoms_header_to_read_withSTIR, randoms_header_to_read_withSTIR)

add_data_offset(randoms_header_to_read_withSTIR, randoms_header_to_read_withSTIR)

# %%

# #### read in again & plot to see if it worked

randoms = stir.ProjData.read_from_file(randoms_header_to_read_withSTIR)

if apply_DOI_adaption: DOI_adaption(randoms, 10)

randoms.write_to_file(os.path.join(STIR_output_folder,randoms_adapted_DOI_filename))

randoms_arr = randoms.as_array()

for i in range(33):

plt.figure()

plot_2d_image([1,1,1],randoms_arr[i, central_slice,:,:],'randoms, TOF bin {}'.format(i))

if i == TOF_bin: plt.savefig(os.path.join(STIR_output_folder,randoms_adapted_DOI_filename[:-3]+'.png'), transparent=False, facecolor='w')

plt.close()

#%%

###################### SCATTER ############################

# alter e7-header of scatter manually to read 2D scatter

seg0_max_rd = 9 # TODO: get this via proj-data-info; currently cannot “downcast” to ProjDataInfoCylindrical

remove_data_offset(scatter_2D_header_to_read_withSTIR, scatter_2D_header_filename)

remove_scan_data_lines_from_interfile_header(scatter_2D_header_to_read_withSTIR, scatter_2D_header_to_read_withSTIR)

replace_siemens_convention_in_interfile_header(scatter_2D_header_to_read_withSTIR, scatter_2D_header_to_read_withSTIR)

change_max_ring_distance(scatter_2D_header_to_read_withSTIR, scatter_2D_header_to_read_withSTIR, seg0_max_rd)

#%%

## read in 2D scatter with STIR

scatter_2D_normalized = stir.ProjData.read_from_file(scatter_2D_header_to_read_withSTIR)

if apply_DOI_adaption: DOI_adaption(scatter_2D_normalized, 10)

#%%

## we hand the object of the final 3D sinogram over to the inverse SSRB function.

scatter_3D_normalized = stir.ProjDataInMemory(prompts_from_e7)

stir.inverse_SSRB(scatter_3D_normalized, scatter_2D_normalized)

#%%

# plot to see if it worked

scatter_3D_norm_arr = scatter_3D_normalized.as_array()

for i in range(33):

plt.figure()

plot_2d_image([1,1,1],scatter_3D_norm_arr[i, central_slice,:,:],'scatter, normalized, TOF bin {}'.format(i))

plt.close()

#%%

## the Siemens output scatter is normalized, so we need to apply the detection efficiencies

scatter_3D_unnormalized = stir.ProjDataInMemory(scatter_3D_normalized)

de = stir.BinNormalisationFromProjData(det_eff)

de.set_up(scatter_3D_unnormalized.get_exam_info(), scatter_3D_unnormalized.get_proj_data_info()) #, randoms.get_proj_data_info())

de.apply(scatter_3D_unnormalized)

scatter_3D_unnormalized.write_to_file(os.path.join(STIR_output_folder,scatter_3D_unnorm_filename))

#%%

scatter_3D_unnormalized_arr = scatter_3D_unnormalized.as_array()

for i in range(33):

plt.figure()

plot_2d_image([1,1,1],scatter_3D_unnormalized_arr[i, central_slice,:,:],'scatter, unnormalized, TOF bin {}'.format(i))

plt.close()

#%%

###################### ADDITIVE TERM ############################

#### the additive term is what is added to the projected estimate, before any multiplicative factors

#### are applied. Therefore, we need to normalize randoms and add (normalized) scatter to it.

add_sino = stir.ProjDataInMemory(prompts_from_e7) #randoms)

add_sino.fill(randoms_arr.flat)

# normalise

norm_projdata = stir.ProjData.read_from_file(os.path.join(STIR_output_folder,norm_filename_fSTIR))

norm = stir.BinNormalisationFromProjData(norm_projdata)

norm.set_up(add_sino.get_exam_info(), add_sino.get_proj_data_info()) #, randoms.get_proj_data_info())

norm.apply(add_sino)

# normalised randoms + scatter

add_sino.sapyb(1, scatter_3D_normalized, 1)

#%%

###################### ATTENUATION CORRECTION FACTORS ############################

# to get the correct additive term, we need to apply attenuation correction.

# we'll use the ones from the e7tools here.

change_datafilename_in_interfile_header(attenuation_corr_factor_header_to_read_withSTIR, prompts_header_filename, attenuation_corr_factor_data_filename)

change_datatype_in_interfile_header(attenuation_corr_factor_header_to_read_withSTIR, 'float', 4)

remove_scan_data_lines_from_interfile_header(attenuation_corr_factor_header_to_read_withSTIR, attenuation_corr_factor_header_to_read_withSTIR)

remove_IMGDATADESC_lines_from_interfile_header(attenuation_corr_factor_header_to_read_withSTIR, attenuation_corr_factor_header_to_read_withSTIR)

remove_data_offset(attenuation_corr_factor_header_to_read_withSTIR, attenuation_corr_factor_header_to_read_withSTIR)

remove_tof_dimension(attenuation_corr_factor_header_to_read_withSTIR, attenuation_corr_factor_header_to_read_withSTIR)

#%%

acf_sino_nonTOF = stir.ProjData.read_from_file(attenuation_corr_factor_header_to_read_withSTIR)

acf_sino_nonTOF = stir.ProjDataInMemory(acf_sino_nonTOF)

if apply_DOI_adaption: DOI_adaption(acf_sino_nonTOF, 10)

#%%

#### expand to TOF as normalization data is TOF

acf_sino = stir.ProjDataInMemory(prompts_from_e7)

ai_arr = acf_sino_nonTOF.as_array()

expanded_arr = np.repeat(ai_arr, 33, axis=0)

acf_sino.fill(expanded_arr.flat)

#%%

###################### ATTENUATION FACTORS ############################

afs = np.divide(1,expanded_arr)

for i in range(33):

plt.figure()

plot_2d_image([1,1,1],afs[i, central_slice,:,:],'afs, TOF bin {}'.format(i))

plt.close()

#%%

AC_correction = stir.BinNormalisationFromProjData(acf_sino)

AC_correction.set_up(add_sino.get_exam_info(), add_sino.get_proj_data_info())

AC_correction.apply(add_sino)

add_sino.write_to_file(os.path.join(STIR_output_folder,additive_term_filename_fSTIR))

#%%

# let's see what it looks like

additive_term_arr = add_sino.as_array()

for i in range(33):

plt.figure()

plot_2d_image([1,1,1],additive_term_arr[i, central_slice,:,:],'additive term, TOF bin {}'.format(i))

if i == TOF_bin: plt.savefig(os.path.join(STIR_output_folder,additive_term_filename_fSTIR[:-3]+'.png'), transparent=False, facecolor='w')

plt.close()

#%%

######### now let's write out the multiplicative factors f. STIR

multi_factors_STIR = stir.ProjDataInMemory(norm_sino)

AC_correction = stir.BinNormalisationFromProjData(acf_sino)

AC_correction.set_up(multi_factors_STIR.get_exam_info(), multi_factors_STIR.get_proj_data_info())

AC_correction.apply(multi_factors_STIR)

multi_factors_STIR.write_to_file(os.path.join(STIR_output_folder,multi_term_filename_fSTIR))

#%%

######### now let's write out the multiplicative factors f. SIRF

multi_factors_SIRF = stir.ProjDataInMemory(det_eff)

AC_correction = stir.BinNormalisationFromProjData(acf_sino)

AC_correction.set_up(multi_factors_SIRF.get_exam_info(), multi_factors_SIRF.get_proj_data_info())

AC_correction.undo(multi_factors_SIRF)

multi_factors_SIRF.write_to_file(os.path.join(STIR_output_folder,multi_term_filename_fSIRF))

#%%

##### to compare the additive term with the acquisition data, we need to

##### pre-correct the prompts with the multiplicative factors

prompts_precorr_f_multi_fact = stir.ProjDataInMemory(prompts_from_e7)

AC_correction = stir.BinNormalisationFromProjData(acf_sino)

AC_correction.set_up(prompts_precorr_f_multi_fact.get_exam_info(), prompts_precorr_f_multi_fact.get_proj_data_info())

AC_correction.apply(prompts_precorr_f_multi_fact)

norm = stir.BinNormalisationFromProjData(norm_projdata)

norm.set_up(prompts_precorr_f_multi_fact.get_exam_info(), prompts_precorr_f_multi_fact.get_proj_data_info())

norm.apply(prompts_precorr_f_multi_fact)

#%%

#### PLOT ADDITIVE TERM

#### draw line-profiles to check if all's correct

prompts_precorr_arr = prompts_precorr_f_multi_fact.as_array()

additive_term_arr = add_sino.as_array()

_, ax = plt.subplots(figsize = (8,6))

prompts_mean = np.mean(prompts_precorr_arr[TOF_bin, central_slice-thickness_half:central_slice+thickness_half, 0, :], axis=(0))

ax.plot(prompts_mean, label="Prompts, pre-corrected f. multi. factors")

add_term_mean = np.mean(additive_term_arr[TOF_bin, central_slice-thickness_half:central_slice+thickness_half, 0, :], axis=(0))

ax.plot(add_term_mean, label="additive term")

ax.set_xlabel('Radial distance (bin)')

ax.set_ylabel('total counts')

ax.set_title('TOF bin:' + str(TOF_bin))

ax.legend()

plt.tight_layout()

plt.suptitle('Lineprofiles - Avg over 10 central slices')

plt.savefig(os.path.join(STIR_output_folder,'additive_term_lineprofiles.png'), transparent=False, facecolor='w')

plt.tight_layout()

# %%

#### PLOT BACKGROUND TERM

#### draw line-profiles to check if all's correct

prompts_arr = prompts_from_e7.as_array()

BG_arr = scatter_3D_unnormalized_arr + randoms_arr

#%%

_, ax = plt.subplots(figsize = (8,6))

ax.plot(np.mean(prompts_arr[TOF_bin, central_slice-thickness_half:central_slice+thickness_half, 0, :], axis=(0)), label='Prompts')

ax.plot(np.mean(BG_arr[TOF_bin, central_slice-thickness_half:central_slice+thickness_half, 0, :], axis=(0)), label='BG term')

ax.set_xlabel('Radial distance (bin)')

ax.set_ylabel('total counts')

ax.set_title('TOF bin:' + str(TOF_bin))

ax.legend()

plt.tight_layout()

plt.suptitle('Lineprofiles - Avg over 10 central slices')

plt.savefig(os.path.join(STIR_output_folder,'BG_term_lineprofiles.png'), transparent=False, facecolor='w')

plt.tight_layout()

# %%

def change_datafilename_in_interfile_header(new_header_filename, header_filename, data_filename):

with open(header_filename) as f:

data = f.read()

poss = re.search(r'name of data file\s*:=[^\n]*', data).span()

data = data.replace(data[poss[0]:poss[1]], \

'name of data file:={}'.format(data_filename))

with open(new_header_filename, 'w') as f2:

f2.write(data)

def DOI_adaption(projdata, DOI_new):

proj_info = projdata.get_proj_data_info()

DOI = proj_info.get_scanner().get_average_depth_of_interaction()

print('Current Depth of interaction:', DOI)

proj_info.get_scanner().set_average_depth_of_interaction(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_intrinsic_azimuthal_tilt()

print('New view offset (rad):', VO)

def check_if_compressed(header_filename):

with open(header_filename) as f:

data = f.read()

match = re.search(r'compression\s*:=\s*(\w+)', data)

if match:

if match.group(1) == 'off':

print('Compression is off, can proceed')

else:

raise ValueError('You are trying to read e7tools compressed data. Please uncompress first!')

else:

print('No compression info found in header!')

def plot_2d_image(idx,vol,title,clims=None,cmap="viridis"):

"""Customized version of subplot to plot 2D image"""

plt.subplot(*idx)

plt.imshow(vol,cmap=cmap)

if clims is not None:

plt.clim(clims)

plt.colorbar(shrink=.5, aspect=.9)

plt.title(title)

plt.axis("off")

def change_datatype_in_interfile_header(header_name, data_type, num_bytes_per_pixel):

with open(header_name) as f:

data = f.read()

poss = re.search(r'number format\s*:=[^\n]*', data).span()

data = data.replace(data[poss[0]:poss[1]], \

'!number format:={}'.format(data_type))

poss = re.search(r'!number of bytes per pixel\s*:=[^\n]*', data).span()

data = data.replace(data[poss[0]:poss[1]], \

'!number of bytes per pixel:={}'.format(num_bytes_per_pixel))

with open(header_name, 'w') as f2:

f2.write(data)

def remove_scan_data_lines_from_interfile_header(header_filename_new, header_filename_old):

with open(header_filename_old) as f:

data = f.read()

data_type_string = r'scan data type description[^\n]*\s*:=\s*[^\n]*\n'

data = re.sub(data_type_string, '', data)

num_data_types_string = r'number of scan data types[^\n]*\s*:=\s*[^\n]*\n'

data = re.sub(num_data_types_string, '', data)

with open(header_filename_new, 'w') as f2:

f2.write(data)

def remove_IMGDATADESC_lines_from_interfile_header(header_filename_new, header_filename_old):

with open(header_filename_old) as f:

data = f.read()

pattern = r'!IMAGE DATA DESCRIPTION:=.*'

data = re.sub(pattern, '', data, flags=re.DOTALL)

with open(header_filename_new, 'w') as f2:

f2.write(data)

def remove_data_offset(header_filename_new, header_filename_old):

with open(header_filename_old) as f:

data = f.read()

data_type_string = r'data offset in bytes[^\n]*\s*:=\s*[^\n]*\n'

data = re.sub(data_type_string, '', data)

with open(header_filename_new, 'w') as f2:

f2.write(data)

def add_data_offset(header_filename_new, header_filename_old):

with open(header_filename_old) as f:

data = f.read()

offset_string = r'\ndata offset in bytes[1]:= 84760000'

pattern = r'(%TOF mashing factor\s*:=[^\n]*)'

data = re.sub(pattern, r'\1' + offset_string, data)

with open(header_filename_new, 'w') as f2:

f2.write(data)

def replace_siemens_convention_in_interfile_header(header_name_new, header_name):

with open(header_name) as f:

data = f.read()

poss = re.search(r'matrix axis label\[2\]:=plane', data).span()

data = data.replace(data[poss[0]:poss[1]], \

'matrix axis label[2]:=sinogram views')

poss = re.search(r'matrix axis label\[3\]:=projection', data).span()

data = data.replace(data[poss[0]:poss[1]], \

'matrix axis label[3]:=number of sinograms')

with open(header_name_new, 'w') as f2:

f2.write(data)

def change_max_ring_distance(header_name_new, header_name, max_ring_diff):

with open(header_name) as f:

data = f.read()

poss = re.search(r'%maximum ring difference\s*:=[^\n]*', data).span()

data = data.replace(data[poss[0]:poss[1]], \

'%maximum ring difference:={}'.format(int(max_ring_diff)))

with open(header_name_new, 'w') as f2:

f2.write(data)

def remove_tof_dimension(header_name_new, header_name):

with open(header_name) as f:

data = f.read()

poss = re.search(r'number of dimensions\s*:=[^\n]*', data).span()

data = data.replace(data[poss[0]:poss[1]], \

'number of dimensions:={}'.format(3))

data_type_string = r'matrix size\[4\]\s*:=[^\n]*\n'

data = re.sub(data_type_string, '', data)

data_type_string = r'matrix axis label\[4\]*\s*:=TOF bin*\n'

data = re.sub(data_type_string, '', data)

data_type_string = r'scale factor \(ps\/bin\) (.*)\n'

data = re.sub(data_type_string, '', data)

data_type_string = r'%TOF mashing factor[^\n]*\s*:=\s*[^\n]*\n'

data = re.sub(data_type_string, '%TOF mashing factor :=0\n', data)

data_type_string = r'%number of TOF time bins[^\n]*\s*:=\s*[^\n]*\n'

data = re.sub(data_type_string, '', data)

with open(header_name_new, 'w') as f2:

f2.write(data)

if __name__ == '__main__':

parser = argparse.ArgumentParser(

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")

parser.add_argument('--norm_sino_data_filename', default='norm3d_00.a', help="The filename of the norm sinogram data")

parser.add_argument('--attenuation_corr_factor_data_filename', default='acf_00.a', help="The filename of the attenuation correction factor data")

parser.add_argument('--STIR_output_folder', default='processing', help="The folder where the output files are written to")

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)