Header file containing the declarations for the MC-GPU code. More...

Classes
struct	int2
struct	int3
struct	float2
struct	float3
struct	double2
struct	double3
struct	short3
struct	ulonglong2
struct	source_struct
	Structure storing the data defining the source model (except for the energy spectrum). More...
struct	source_energy_struct
	Structure storing the source energy spectrum data to be sampled using the Walker aliasing algorithm. More...
struct	detector_struct
	Structure storing the data defining the x-ray detector. More...
struct	voxel_struct
	Structure defining a voxelized box with the back-lower corner at the coordinate origin. More...
struct	linear_interp
	Structure with the basic data required by the linear interpolation of the mean free paths: number of values and energy grid. More...
struct	compton_struct
	Structure storing the data of the Compton interaction sampling model (equivalent to PENELOPE's common block /CGCO/). More...
struct	rayleigh_struct
	Structure storing the data of the Rayleigh interaction sampling model (equivalent to PENELOPE's common block /CGRA/). More...
Defines
#define	MASTER_THREAD if(0==myID)
	MPI macro: mark commands to be executed only by the master thread (myID==0).
#define	MAX_NUM_PROJECTIONS 720
	Maximum number of projections allowed in the CT simulation (not limited by the constant memory because stored in global and shared memory):
#define	MAX_MATERIALS 15
	Constants values for the Compton and Rayleigh models:
#define	MAX_SHELLS 30
#define	NP_RAYLEIGH 128
#define	MAX_ENERGYBINS_RAYLEIGH 25005
#define	MAX_ENERGY_BINS 256
	Maximum number of energy bins in the input x-ray energy spectrum.
#define	PI 3.14159265358979323846
#define	RAD2DEG 180.0/PI
#define	DEG2RAD PI/180.0
#define	SCALE_eV 100.0f
	Value to scale the deposited energy in the pixels so that it can be stored as a long long integer instead of a double precision float.
#define	EPS_SOURCE 0.000015f
	Offset value to make sure the particles are completely inside, or outside, the voxel bounding box.
#define	NEG_EPS_SOURCE -EPS_SOURCE
#define	INF 500000.0f
#define	INF_minus1 499999.0f
#define	NEG_INF -500000.0f
#define	max_value(a, b) ( ((a) > (b)) ? (a) : (b) )
	Preprocessor macro to calculate maximum and minimum values:
#define	min_value(a, b) ( ((a) < (b)) ? (a) : (b) )
Typedefs
typedef struct int2	int2
typedef struct int3	int3
typedef struct float2	float2
typedef struct float3	float3
typedef struct double2	double2
typedef struct double3	double3
typedef struct short3	short3
typedef struct ulonglong2	ulonglong2
Functions
void	read_input (int argc, char *argv, int myID, unsigned long long int total_histories, int gpu_id, int seed_input, int num_threads_per_block, int histories_per_thread, struct detector_struct detector_data, unsigned long long int image_ptr, int image_bytes, struct source_struct source_data, struct source_energy_struct source_energy_data, char file_name_voxels, char file_name_materials[MAX_MATERIALS][250], char file_name_output, char file_name_espc, int num_projections, double D_angle, double angularROI_0, double angularROI_1, double initial_angle, ulonglong2 *voxels_Edep_ptr, int voxels_Edep_bytes, char file_dose_output, short int dose_ROI_x_min, short int dose_ROI_x_max, short int dose_ROI_y_min, short int dose_ROI_y_max, short int dose_ROI_z_min, short int dose_ROI_z_max, double SRotAxisD, double vertical_translation_per_projection, int flag_material_dose)
	Read the input file given in the command line and return the significant data.
void	load_voxels (int myID, char file_name_voxels, float density_max, struct voxel_struct voxel_data, float2 voxel_mat_dens_ptr, unsigned int voxel_mat_dens_bytes, short int dose_ROI_x_max, short int dose_ROI_y_max, short int *dose_ROI_z_max)
	Read the voxel data and allocate the material and density matrix.
void	load_material (int myID, char file_name_materials[MAX_MATERIALS][250], float density_max, float density_nominal, struct linear_interp mfp_table_data, float2 mfp_Woodcock_table, int mfp_Woodcock_table_bytes, float3 mfp_table_a_ptr, float3 mfp_table_b_ptr, int mfp_table_bytes, struct rayleigh_struct rayleigh_table_ptr, struct compton_struct *compton_table_ptr)
	Read the material input files and set the mean free paths and the "linear_interp" structures.
void	trim_name (char input_line, char file_name)
	Extract a file name from an input text line, trimming the initial blanks, trailing comment (#) and stopping at the first blank (the file name should not contain blanks).
char *	fgets_trimmed (char trimmed_line, int num, FILE file_ptr)
	Read a line of text and trim initial blancks and trailing comments (#).
int	report_image (char file_name_output, struct detector_struct detector_data, struct source_struct source_data, float mean_energy_spectrum, unsigned long long int image, double time_elapsed, unsigned long long int total_histories, int current_projection, int num_projections, double D_angle, double initial_angle, int myID, int numprocs)
	Report the tallied image in ASCII and binary form (32-bit floats).
void	set_CT_trajectory (int myID, int num_projections, double D_angle, double angularROI_0, double angularROI_1, double SRotAxisD, struct source_struct source_data, struct detector_struct detector_data, double vertical_translation_per_projection)
	Sets the CT trajectory: store in memory the source and detector rotations that are needed to calculate the multiple projections.
int	report_voxels_dose (char file_dose_output, int num_projections, struct voxel_struct voxel_data, float2 voxel_mat_dens, ulonglong2 voxels_Edep, double time_elapsed_total, unsigned long long int total_histories, short int dose_ROI_x_min, short int dose_ROI_x_max, short int dose_ROI_y_min, short int dose_ROI_y_max, short int dose_ROI_z_min, short int dose_ROI_z_max, struct source_struct *source_data)
	Report the total tallied 3D voxel dose deposition for all projections.
void	init_energy_spectrum (char file_name_espc, struct source_energy_struct source_energy_data, float *mean_energy_spectrum)
	Read the energy spectrum file and initialize the Walker aliasing sampling.
void	update_seed_PRNG (int batch_number, unsigned long long int total_histories, int *seed)
	Initialize the first seed of the pseudo-random number generator (PRNG) RANECU to a position far away from the previous history (leap frog technique).
void	IRND0 (float W, float F, short int *K, int N)
	Initialisation of Walker's aliasing algorithm for random sampling from discrete probability distributions.
int	report_materials_dose (int num_projections, unsigned long long int total_histories, float density_nominal, ulonglong2 materials_dose, double *mass_materials)
	Report the tallied dose to each material number, accounting for different densities in different regions with the same material number.
int	seeki_walker (float cutoff, short int alias, float randno, int n)
	Finds the interval (x(i),x(i+1)] containing the input value using Walker's aliasing method.
void	track_particles (int history_batch, int histories_per_thread, int num_p, int seed_input, unsigned long long int image, ulonglong2 voxels_Edep, float2 voxel_mat_dens, float2 mfp_Woodcock_table, float3 mfp_table_a, float3 mfp_table_b, struct rayleigh_struct rayleigh_table, struct compton_struct compton_table, struct detector_struct detector_data_array, struct source_struct source_data_array, ulonglong2 *materials_dose)
	Initialize the image array, ie, set all pixels to zero Essentially, this function has the same effect as the command: "cutilSafeCall(cudaMemcpy(image_device, image, image_bytes, cudaMemcpyHostToDevice))";.
void	source (float3 position, float3 direction, float energy, int2 seed, int absvox, struct source_struct source_data_SHARED, struct detector_struct *detector_data_SHARED)
	Source that creates primary x rays, according to the defined source model.
void	move_to_bbox (float3 position, float3 direction, float3 size_bbox, int *intersection_flag)
	Functions that moves a particle inside the voxelized geometry bounding box.
void	tally_image (float energy, float3 position, float3 direction, signed char scatter_state, unsigned long long int image, struct source_struct source_data_SHARED, struct detector_struct *detector_data_SHARED)
	Tally a radiographic projection image.
void	init_PRNG (int history_batch, int histories_per_thread, int seed_input, int2 *seed)
	Initialize the pseudo-random number generator (PRNG) RANECU to a position far away from the previous history (leap frog technique).
int	abMODm (int m, int a, int s)
	Calculate "(a1*a2) MOD m" with 32-bit integers and avoiding the possible overflow, using the Russian Peasant approach modulo m and the approximate factoring method, as described in: L'Ecuyer and Cote, ACM Trans.
float	ranecu (int2 *seed)
	Pseudo-random number generator (PRNG) RANECU returning a float value (single precision version).
double	ranecu_double (int2 *seed)
	Pseudo-random number generator (PRNG) RANECU returning a double value.
int	locate_voxel (float3 position, short3 voxel_coord)
	Find the voxel that contains the current position.
void	rotate_double (float3 *direction, double cos_theta, double phi)
	Rotates a vector; the rotation is specified by giving the polar and azimuthal angles in the "self-frame", as determined by the vector to be rotated.
void	GRAa (float energy, double costh_Rayleigh, int mat, float pmax_current, int2 seed, struct rayleigh_struct cgra)
	Sample a Rayleigh interaction using the sampling algorithm used in PENELOPE 2006.
void	GCOa (float energy, double costh_Compton, int mat, int2 seed, struct compton_struct *cgco_SHARED)
	Random sampling of incoherent (Compton) scattering of photons, using the sampling algorithm from PENELOPE 2006: Relativistic impulse approximation with analytical one-electron Compton profiles.
void	tally_voxel_energy_deposition (float Edep, short3 voxel_coord, ulonglong2 *dose)
	Tally the dose deposited in the voxels.
void	tally_materials_dose (float Edep, int material, ulonglong2 *materials_dose)
	Tally the depose deposited inside each material.
Variables
short int	dose_ROI_x_min_CONST
	Global variable to be stored in the GPU constant memory defining the coordinates of the dose deposition region of interest.
short int	dose_ROI_x_max_CONST
short int	dose_ROI_y_min_CONST
short int	dose_ROI_y_max_CONST
short int	dose_ROI_z_min_CONST
short int	dose_ROI_z_max_CONST
struct voxel_struct	voxel_data_CONST
	Global variable to be stored in the GPU constant memory defining the size of the voxel phantom.
struct source_struct	source_data_CONST
	Global variable to be stored in the GPU constant memory defining the x-ray source.
struct detector_struct	detector_data_CONST
	Global variable to be stored in the GPU constant memory defining the x-ray detector.
struct linear_interp	mfp_table_data_CONST
	Global variable to be stored in the GPU constant memory defining the linear interpolation data.
struct source_energy_struct	source_energy_data_CONST
	Global variable to be stored in the GPU constant memory defining the source energy spectrum.

Detailed Description

Header file containing the declarations for the MC-GPU code.

This file declares all the host and device functions and structures, the library files to include in the compilation, various constants parameters of the simulation, pre-processor macro functions, etc.

Author:: Andreu Badal (Andreu.Badal-Soler@fda.hhs.gov)

Date:: 2012/12/12

Definition in file MC-GPU_v1.3.h.

Define Documentation

#define DEG2RAD PI/180.0

Definition at line 73 of file MC-GPU_v1.3.h.

Referenced by read_input().

#define EPS_SOURCE 0.000015f

Offset value to make sure the particles are completely inside, or outside, the voxel bounding box.

For example, to start a particle track the source may have to translate a particle from the focal spot to the plane Y=0, but in reality the particle will be moved to Y=EPS_SOURCE to make sure it is unmistakenly located inside a voxel. If "EPS_SOURCE" is too small an artifact with weird concentric circular shadows may appear at the center of the simulated image.

Definition at line 85 of file MC-GPU_v1.3.h.

Referenced by locate_voxel(), and move_to_bbox().

#define INF 500000.0f

Definition at line 88 of file MC-GPU_v1.3.h.

#define INF_minus1 499999.0f

Definition at line 89 of file MC-GPU_v1.3.h.

#define MASTER_THREAD if(0==myID)

MPI macro: mark commands to be executed only by the master thread (myID==0).

Definition at line 56 of file MC-GPU_v1.3.h.

Referenced by load_material(), load_voxels(), main(), read_input(), and set_CT_trajectory().

#define MAX_ENERGY_BINS 256

Maximum number of energy bins in the input x-ray energy spectrum.

Definition at line 68 of file MC-GPU_v1.3.h.

Referenced by init_energy_spectrum().

#define MAX_ENERGYBINS_RAYLEIGH 25005

Definition at line 65 of file MC-GPU_v1.3.h.

Referenced by load_material().

#define MAX_MATERIALS 15

Constants values for the Compton and Rayleigh models:

Definition at line 62 of file MC-GPU_v1.3.h.

Referenced by GCOa(), load_material(), load_voxels(), main(), read_input(), report_materials_dose(), report_voxels_dose(), and track_particles().

#define MAX_NUM_PROJECTIONS 720

Maximum number of projections allowed in the CT simulation (not limited by the constant memory because stored in global and shared memory):

Definition at line 59 of file MC-GPU_v1.3.h.

Referenced by read_input(), and set_CT_trajectory().

#define MAX_SHELLS 30

Definition at line 63 of file MC-GPU_v1.3.h.

Referenced by GCOa(), and load_material().

#define max_value	(	a,
		b
	)	( ((a) > (b)) ? (a) : (b) )

Preprocessor macro to calculate maximum and minimum values:

Definition at line 93 of file MC-GPU_v1.3.h.

Referenced by GCOa(), and init_energy_spectrum().

#define min_value	(	a,
		b
	)	( ((a) < (b)) ? (a) : (b) )

Definition at line 94 of file MC-GPU_v1.3.h.

Referenced by GCOa(), GRAa(), and load_voxels().

#define NEG_EPS_SOURCE -EPS_SOURCE

Definition at line 87 of file MC-GPU_v1.3.h.

Referenced by move_to_bbox().

#define NEG_INF -500000.0f

Definition at line 90 of file MC-GPU_v1.3.h.

Referenced by move_to_bbox().

#define NP_RAYLEIGH 128

Definition at line 64 of file MC-GPU_v1.3.h.

Referenced by GRAa(), and load_material().

#define PI 3.14159265358979323846

Definition at line 71 of file MC-GPU_v1.3.h.

Referenced by read_input(), and set_CT_trajectory().

#define RAD2DEG 180.0/PI

Definition at line 72 of file MC-GPU_v1.3.h.

Referenced by main(), read_input(), report_image(), and set_CT_trajectory().

#define SCALE_eV 100.0f

Value to scale the deposited energy in the pixels so that it can be stored as a long long integer instead of a double precision float.

The integer values have to be used in order to use the atomicAdd function in CUDA.

Definition at line 79 of file MC-GPU_v1.3.h.

Referenced by report_image(), report_materials_dose(), report_voxels_dose(), tally_image(), tally_materials_dose(), and tally_voxel_energy_deposition().

Typedef Documentation

typedef struct double2 double2

Definition at line 134 of file MC-GPU_v1.3.h.

typedef struct double3 double3

Definition at line 135 of file MC-GPU_v1.3.h.

typedef struct float2 float2

Definition at line 132 of file MC-GPU_v1.3.h.

typedef struct float3 float3

Definition at line 133 of file MC-GPU_v1.3.h.

typedef struct int2 int2

Definition at line 130 of file MC-GPU_v1.3.h.

typedef struct int3 int3

Definition at line 131 of file MC-GPU_v1.3.h.

typedef struct short3 short3

Definition at line 136 of file MC-GPU_v1.3.h.

typedef struct ulonglong2 ulonglong2

Definition at line 137 of file MC-GPU_v1.3.h.

Function Documentation

int abMODm	(	int	m,
		int	a,
		int	s
	)		`[inline]`

Calculate "(a1*a2) MOD m" with 32-bit integers and avoiding the possible overflow, using the Russian Peasant approach modulo m and the approximate factoring method, as described in: L'Ecuyer and Cote, ACM Trans.

Math. Soft. 17 (1991).

This function has been adapted from "seedsMLCG.f", see: Badal and Sempau, Computer Physics Communications 175 (2006)

Parameters:

[in] m,a,s MLCG parameters

Returns:: (a1*a2) MOD m

Definition at line 919 of file MC-GPU_kernel_v1.3.cu.

Referenced by init_PRNG(), and update_seed_PRNG().

char* fgets_trimmed	(	char *	trimmed_line,
		int	num,
		FILE *	file_ptr
	)

Read a line of text and trim initial blancks and trailing comments (#).

Parameters:

[in]	num	Characters to read
[in]	file_ptr	Pointer to the input file stream
[out]	trimmed_line	Trimmed line from input file, skipping empty lines and comments

Definition at line 1868 of file MC-GPU_v1.3.cu.

Referenced by init_energy_spectrum(), and read_input().

void GCOa	(	float *	energy,
		double *	costh_Compton,
		int *	mat,
		int2 *	seed,
		struct compton_struct *	cgco_SHARED
	)		`[inline]`

Random sampling of incoherent (Compton) scattering of photons, using the sampling algorithm from PENELOPE 2006: Relativistic impulse approximation with analytical one-electron Compton profiles.

Parameters:

[in,out]	energy	incident and final photon energy (eV)
[out]	costh_Compton	cosine of the polar scattering angle
[in]	material	Current voxel material
[in]	seed	RANECU PRNG seed

Definition at line 1287 of file MC-GPU_kernel_v1.3.cu.

References compton_struct::fco, compton_struct::fj0, MAX_MATERIALS, MAX_SHELLS, max_value, min_value, compton_struct::noscco, ranecu(), and compton_struct::uico.

Referenced by track_particles().

void GRAa	(	float *	energy,
		double *	costh_Rayleigh,
		int *	mat,
		float *	pmax_current,
		int2 *	seed,
		struct rayleigh_struct *	cgra
	)		`[inline]`

Sample a Rayleigh interaction using the sampling algorithm used in PENELOPE 2006.

Parameters:

[in]	energy	Particle energy (not modified with Rayleigh)
[out]	costh_Rayleigh	Cosine of the angular deflection
[in]	material	Current voxel material

Definition at line 1181 of file MC-GPU_kernel_v1.3.cu.

References rayleigh_struct::aco, rayleigh_struct::bco, rayleigh_struct::itlco, rayleigh_struct::ituco, min_value, NP_RAYLEIGH, rayleigh_struct::pco, ranecu_double(), and rayleigh_struct::xco.

Referenced by track_particles().

void init_energy_spectrum	(	char *	file_name_espc,
		struct source_energy_struct *	source_energy_data,
		float *	mean_energy_spectrum
	)

Read the energy spectrum file and initialize the Walker aliasing sampling.

Parameters:

[in]	file_name_espc	File containing the energy spectrum (lower energy value in each bin and its emission probability).
[in,out]	source_energy_data	Energy spectrum and other source data. The Walker alias and cutoffs are initialized in this function.
[out]	mean_energy_spectrum	Mean energy in the input x-ray energy spectrum.

!Walker!! Calling PENELOPE's function to init the Walker method

Definition at line 3398 of file MC-GPU_v1.3.cu.

References source_energy_struct::espc, source_energy_struct::espc_alias, source_energy_struct::espc_cutoff, fgets_trimmed(), IRND0(), MAX_ENERGY_BINS, max_value, and source_energy_struct::num_bins_espc.

Referenced by main().

void init_PRNG	(	int	history_batch,
		int	histories_per_thread,
		int	seed_input,
		int2 *	seed
	)		`[inline]`

Initialize the pseudo-random number generator (PRNG) RANECU to a position far away from the previous history (leap frog technique).

Each calculated seed initiates a consecutive and disjoint sequence of pseudo-random numbers with length LEAP_DISTANCE, that can be used to in a parallel simulation (Sequence Splitting parallelization method). The basic equation behind the algorithm is: S(i+j) = (a**j * S(i)) MOD m = [(a**j MOD m)*S(i)] MOD m , which is described in: P L'Ecuyer, Commun. ACM 31 (1988) p.742

This function has been adapted from "seedsMLCG.f", see: A Badal and J Sempau, Computer Physics Communications 175 (2006) p. 440-450

Parameters:

[in]	history	Particle bach number.
[in]	seed_input	Initial PRNG seed input (used to initiate both MLCGs in RANECU).
[out]	seed	Initial PRNG seeds for the present history.

Definition at line 841 of file MC-GPU_kernel_v1.3.cu.

References a1_RANECU, a2_RANECU, abMODm(), LEAP_DISTANCE, m1_RANECU, m2_RANECU, int2::x, and int2::y.

Referenced by track_particles().

void IRND0	(	float *	W,
		float *	F,
		short int *	K,
		int	N
	)

Initialisation of Walker's aliasing algorithm for random sampling from discrete probability distributions.

Input arguments: N ........ number of different values of the random variable. W(1:N) ... corresponding point probabilities (not necessarily normalised to unity). Output arguments: F(1:N) ... cutoff values. K(1:N) ... alias values.

This subroutine is part of the PENELOPE 2006 code developed by Francesc Salvat at the University of Barcelona. For more info: www.oecd-nea.org/science/pubs/2009/nea6416-penelope.pdf

Definition at line 3575 of file MC-GPU_v1.3.cu.

Referenced by init_energy_spectrum().

void load_material	(	int	myID,
		char	file_name_materials[MAX_MATERIALS][250],
		float *	density_max,
		float *	density_nominal,
		struct linear_interp *	mfp_table_data,
		float2 **	mfp_Woodcock_table_ptr,
		int *	mfp_Woodcock_table_bytes,
		float3 **	mfp_table_a_ptr,
		float3 **	mfp_table_b_ptr,
		int *	mfp_table_bytes,
		struct rayleigh_struct *	rayleigh_table_ptr,
		struct compton_struct *	compton_table_ptr
	)

Read the material input files and set the mean free paths and the "linear_interp" structures.

Find the material nominal density. Set the Woodcock trick data.

Parameters:

[in]	file_name_materials	Array with the names of the material files.
[in]	density_max	maximum density in the geometry (needed to set Woodcock trick)
[out]	density_nominal	Array with the nominal density of the materials read
[out]	mfp_table_data	Constant values for the linear interpolation
[out]	mfp_table_a_ptr	First element for the linear interpolation.
[out]	mfp_table_b_ptr	Second element for the linear interpolation.

Definition at line 2110 of file MC-GPU_v1.3.cu.

References rayleigh_struct::aco, rayleigh_struct::bco, linear_interp::e0, compton_struct::fco, compton_struct::fj0, linear_interp::ide, rayleigh_struct::itlco, rayleigh_struct::ituco, MASTER_THREAD, MAX_ENERGYBINS_RAYLEIGH, MAX_MATERIALS, MAX_SHELLS, compton_struct::noscco, NP_RAYLEIGH, linear_interp::num_values, rayleigh_struct::pco, rayleigh_struct::pmax, compton_struct::uico, and rayleigh_struct::xco.

Referenced by main().

void load_voxels	(	int	myID,
		char *	file_name_voxels,
		float *	density_max,
		struct voxel_struct *	voxel_data,
		float2 **	voxel_mat_dens_ptr,
		unsigned int *	voxel_mat_dens_bytes,
		short int *	dose_ROI_x_max,
		short int *	dose_ROI_y_max,
		short int *	dose_ROI_z_max
	)

Read the voxel data and allocate the material and density matrix.

Also find and report the maximum density defined in the geometry.

Parameters:

[in]	file_name_voxels	Name of the voxelized geometry file.
[out]	density_max	Array with the maximum density for each material in the voxels.
[out]	voxel_data	Pointer to a structure containing the voxel number and size.
[out]	voxel_mat_dens_ptr	Pointer to the vector with the voxel materials and densities.
[in]	dose_ROI_x/y/z_max	Size of the dose ROI: can not be larger than the total number of voxels in the geometry.

Definition at line 1929 of file MC-GPU_v1.3.cu.

References voxel_struct::inv_voxel_size, MASTER_THREAD, MAX_MATERIALS, min_value, voxel_struct::num_voxels, voxel_struct::size_bbox, int3::x, float3::x, int3::y, float3::y, int3::z, and float3::z.

Referenced by main().

int locate_voxel	(	float3 *	position,
		short3 *	voxel_coord
	)		`[inline]`

Find the voxel that contains the current position.

Report the voxel absolute index and the x,y,z indices. The structure containing the voxel number and size is read from CONSTANT memory.

Parameters:

[in]	position	Particle position
[out]	voxel_coord	Pointer to three integer values (short3*) that will store the x,y and z voxel indices.

Returns:: Returns "absvox", the voxel number where the particle is located (negative if position outside the voxel bbox).

Definition at line 1033 of file MC-GPU_kernel_v1.3.cu.

References EPS_SOURCE, voxel_struct::inv_voxel_size, voxel_struct::num_voxels, voxel_struct::size_bbox, voxel_data_CONST, int3::x, float3::x, short3::x, int3::y, float3::y, short3::y, float3::z, and short3::z.

Referenced by track_particles().

void move_to_bbox	(	float3 *	position,
		float3 *	direction,
		float3	size_bbox,
		int *	intersection_flag
	)		`[inline]`

Functions that moves a particle inside the voxelized geometry bounding box.

An EPSILON distance is added to make sure the particles will be clearly inside the bbox, not exactly on the surface.

This algorithm makes the following assumtions:

The back lower vertex of the voxel bounding box is always located at the origin: (x0,y0,z0)=(0,0,0).
The initial value of "position" corresponds to the focal spot location.
When a ray is not pointing towards the bbox plane that it should cross according to the sign of the direction, I assign a distance to the intersection =0 instead of the real negative distance. The wall that will be crossed to enter the bbox is always the furthest and therefore a 0 distance will never be used except in the case of a ray starting inside the bbox or outside the bbox and not pointing to any of the 3 planes. In this situation the ray will be transported a 0 distance, meaning that it will stay at the focal spot.

(Interesting information on ray-box intersection: http://tog.acm.org/resources/GraphicsGems/gems/RayBox.c)

Parameters:

[in,out]	position	Particle position: initially set to the focal spot, returned transported inside the voxel bbox.
[out]	direction	Sampled particle direction (cosine vectors).
[out]	intersection_flag	Set to <0 if particle outside bbox and will not cross the voxels, not changed otherwise.
[out]	size_bbox	Size of the bounding box.

Definition at line 714 of file MC-GPU_kernel_v1.3.cu.

References EPS_SOURCE, NEG_EPS_SOURCE, NEG_INF, float3::x, float3::y, and float3::z.

Referenced by source().

float ranecu ( int2 * seed ) [inline]

Pseudo-random number generator (PRNG) RANECU returning a float value (single precision version).

Parameters:

[in,out] seed PRNG seed (seed kept in the calling function and updated here).

Returns:: PRN double value in the open interval (0,1)

Definition at line 965 of file MC-GPU_kernel_v1.3.cu.

References int2::x, and int2::y.

Referenced by GCOa(), source(), and track_particles().

double ranecu_double ( int2 * seed ) [inline]

Pseudo-random number generator (PRNG) RANECU returning a double value.

Definition at line 995 of file MC-GPU_kernel_v1.3.cu.

References int2::x, and int2::y.

Referenced by GRAa(), and track_particles().

void read_input	(	int	argc,
		char **	argv,
		int	myID,
		unsigned long long int *	total_histories,
		int *	seed_input,
		int *	gpu_id,
		int *	num_threads_per_block,
		int *	histories_per_thread,
		struct detector_struct *	detector_data,
		unsigned long long int **	image_ptr,
		int *	image_bytes,
		struct source_struct *	source_data,
		struct source_energy_struct *	source_energy_data,
		char *	file_name_voxels,
		char	file_name_materials[MAX_MATERIALS][250],
		char *	file_name_output,
		char *	file_name_espc,
		int *	num_projections,
		double *	D_angle,
		double *	angularROI_0,
		double *	angularROI_1,
		double *	initial_angle,
		ulonglong2 **	voxels_Edep_ptr,
		int *	voxels_Edep_bytes,
		char *	file_dose_output,
		short int *	dose_ROI_x_min,
		short int *	dose_ROI_x_max,
		short int *	dose_ROI_y_min,
		short int *	dose_ROI_y_max,
		short int *	dose_ROI_z_min,
		short int *	dose_ROI_z_max,
		double *	SRotAxisD,
		double *	vertical_translation_per_projection,
		int *	flag_material_dose
	)

Read the input file given in the command line and return the significant data.

Example input file:

1000000 [Total number of histories to simulate] geometry.vox [Voxelized geometry file name] material.mat [Material data file name]

Parameters:

[in]	argc	Command line parameters
[in]	argv	Command line parameters: name of input file
[out]	total_histories	Total number of particles to simulate
[out]	seed_input	Input random number generator seed
[out]	num_threads_per_block	Number of CUDA threads for each GPU block
[out]	detector_data
[out]	image
[out]	source_data
[out]	file_name_voxels
[out]	file_name_materials
[out]	file_name_output

Definition at line 1242 of file MC-GPU_v1.3.cu.

References detector_struct::center, detector_struct::corner_min_rotated_to_Y, source_struct::cos_theta_low, DEG2RAD, source_struct::direction, fgets_trimmed(), detector_struct::height_Z, detector_struct::inv_pixel_size_X, detector_struct::inv_pixel_size_Z, MASTER_THREAD, MAX_MATERIALS, MAX_NUM_PROJECTIONS, detector_struct::num_pixels, PI, source_struct::position, RAD2DEG, detector_struct::rot_inv, detector_struct::rotation_flag, detector_struct::sdd, detector_struct::total_num_pixels, trim_name(), detector_struct::width_X, int2::x, float3::x, int2::y, float3::y, and float3::z.

Referenced by main().

int report_image	(	char *	file_name_output,
		struct detector_struct *	detector_data,
		struct source_struct *	source_data,
		float	mean_energy_spectrum,
		unsigned long long int *	image,
		double	time_elapsed,
		unsigned long long int	total_histories,
		int	current_projection,
		int	num_projections,
		double	D_angle,
		double	initial_angle,
		int	myID,
		int	numprocs
	)

Report the tallied image in ASCII and binary form (32-bit floats).

Separate images for primary and scatter radiation are generated.

Parameters:

[in]	file_name_output	File where tallied image is reported
[in]	detector_data	Detector description read from the input file (pointer to detector_struct)
[in]	image	Tallied image (in meV per pixel)
[in]	time_elapsed	Time elapsed during the main loop execution (in seconds)
[in]	total_histories	Total number of x-rays simulated

Definition at line 2715 of file MC-GPU_v1.3.cu.

References source_struct::direction, detector_struct::inv_pixel_size_X, detector_struct::inv_pixel_size_Z, detector_struct::num_pixels, source_struct::position, RAD2DEG, SCALE_eV, int2::x, float3::x, int2::y, float3::y, and float3::z.

Referenced by main().

int report_materials_dose	(	int	num_projections,
		unsigned long long int	total_histories,
		float *	density_nominal,
		ulonglong2 *	materials_dose,
		double *	mass_materials
	)

Report the tallied dose to each material number, accounting for different densities in different regions with the same material number.

Parameters:

[in]	num_projections	Number of projections simulated
[in]	total_histories	Total number of x-rays simulated per projection
[out]	density_nominal	Array with the nominal densities of materials given in the input file; -1 for materials not defined. Used to report only defined materials.
[in]	materials_dose	Tallied dose and dose^2 arrays

Definition at line 3128 of file MC-GPU_v1.3.cu.

References MAX_MATERIALS, and SCALE_eV.

Referenced by main().

int report_voxels_dose	(	char *	file_dose_output,
		int	num_projections,
		struct voxel_struct *	voxel_data,
		float2 *	voxel_mat_dens,
		ulonglong2 *	voxels_Edep,
		double	time_total_MC_init_report,
		unsigned long long int	total_histories,
		short int	dose_ROI_x_min,
		short int	dose_ROI_x_max,
		short int	dose_ROI_y_min,
		short int	dose_ROI_y_max,
		short int	dose_ROI_z_min,
		short int	dose_ROI_z_max,
		struct source_struct *	source_data
	)

Report the total tallied 3D voxel dose deposition for all projections.

The voxel doses in the input ROI and their respective uncertainties are reported in binary form (32-bit floats) in two separate .raw files. The dose in a single plane at the level of the focal spot is also reported in ASCII format for simple visualization with GNUPLOT. The total dose deposited in each different material is reported to the standard output. The material dose is calculated adding the energy deposited in the individual voxels within the dose ROI, and dividing by the total mass of the material in the ROI.

Parameters:

[in]	file_dose_output	File where tallied image is reported
[in]	detector_data	Detector description read from the input file (pointer to detector_struct)
[in]	image	Tallied image (in meV per pixel)
[in]	time_elapsed	Time elapsed during the main loop execution (in seconds)
[in]	total_histories	Total number of x-rays simulated
[in]	source_data	Data required to compute the voxel plane to report in ASCII format: Z at the level of the source, 1st projection

Definition at line 2890 of file MC-GPU_v1.3.cu.

References voxel_struct::inv_voxel_size, MAX_MATERIALS, voxel_struct::num_voxels, SCALE_eV, int3::x, float3::x, ulonglong2::x, int3::y, float2::y, float3::y, and float3::z.

Referenced by main().

void rotate_double	(	float3 *	direction,
		double	costh,
		double	phi
	)		`[inline]`

Rotates a vector; the rotation is specified by giving the polar and azimuthal angles in the "self-frame", as determined by the vector to be rotated.

This function is a literal translation from Fortran to C of PENELOPE (v. 2006) subroutine "DIRECT".

Parameters:

[in,out]	(u,v,w)	input vector (=d) in the lab. frame; returns the rotated vector components in the lab. frame
[in]	costh	cos(theta), angle between d before and after turn
[in]	phi	azimuthal angle (rad) turned by d in its self-frame

Definition at line 1103 of file MC-GPU_kernel_v1.3.cu.

References float3::x, float3::y, and float3::z.

Referenced by track_particles().

int seeki_walker	(	float *	cutoff,
		short int *	alias,
		float	randno,
		int	n
	)		`[inline]`

Finds the interval (x(i),x(i+1)] containing the input value using Walker's aliasing method.

Input: cutoff(1..n) -> interval cutoff values for the Walker method cutoff(1..n) -> alias for the upper part of each interval randno -> point to be located n -> no. of data points Output: index i of the semiopen interval where randno lies Comments: -> The cutoff and alias values have to be previously initialised calling the penelope subroutine IRND0.

Algorithm implementation based on the PENELOPE code developed by Francesc Salvat at the University of Barcelona. For more info: www.oecd-nea.org/science/pubs/2009/nea6416-penelope.pdf

Definition at line 3526 of file MC-GPU_v1.3.cu.

void set_CT_trajectory	(	int	myID,
		int	num_projections,
		double	D_angle,
		double	angularROI_0,
		double	angularROI_1,
		double	SRotAxisD,
		struct source_struct *	source_data,
		struct detector_struct *	detector_data,
		double	vertical_translation_per_projection
	)

Sets the CT trajectory: store in memory the source and detector rotations that are needed to calculate the multiple projections.

The first projection (0) was previously initialized in function "read_input".

ASSUMPTIONS: the CT scan plane must be perpendicular to the Z axis, ie, the initial direction of the particles must have w=0!

Definition at line 3194 of file MC-GPU_v1.3.cu.

References detector_struct::center, detector_struct::corner_min_rotated_to_Y, source_struct::direction, detector_struct::height_Z, detector_struct::inv_pixel_size_X, detector_struct::inv_pixel_size_Z, MASTER_THREAD, MAX_NUM_PROJECTIONS, detector_struct::num_pixels, PI, source_struct::position, RAD2DEG, detector_struct::rot_inv, detector_struct::rotation_flag, detector_struct::sdd, detector_struct::total_num_pixels, detector_struct::width_X, float3::x, float3::y, and float3::z.

Referenced by main().

void source	(	float3 *	position,
		float3 *	direction,
		float *	energy,
		int2 *	seed,
		int *	absvox,
		struct source_struct *	source_data_SHARED,
		struct detector_struct *	detector_data_SHARED
	)		`[inline]`

Source that creates primary x rays, according to the defined source model.

The particles are automatically moved to the surface of the voxel bounding box, to start the tracking inside a real material. If the sampled particle do not enter the voxels, it is init in the focal spot and the main program will check if it arrives at the detector or not.

Parameters:

[in]	source_data	Structure describing the source.
[in]	source_energy_data_CONST	Global variable in constant memory space describing the source energy spectrum.
[out]	position	Initial particle position (particle transported inside the voxel bbox).
[out]	direction	Sampled particle direction (cosine vectors).
[out]	energy	Sampled energy of the new x ray.
[in]	seed	Current seed of the random number generator, requiered to sample the movement direction.
[out]	absvox	Set to <0 if primary particle will not cross the voxels, not changed otherwise (>0).

Definition at line 626 of file MC-GPU_kernel_v1.3.cu.

References source_struct::cos_theta_low, source_struct::D_cos_theta, source_struct::D_phi, source_energy_struct::espc, source_energy_struct::espc_alias, source_energy_struct::espc_cutoff, source_struct::max_height_at_y1cm, move_to_bbox(), source_energy_struct::num_bins_espc, source_struct::phi_low, source_struct::position, ranecu(), source_struct::rot_fan, detector_struct::rotation_flag, voxel_struct::size_bbox, source_energy_data_CONST, voxel_data_CONST, float3::x, float3::y, and float3::z.

Referenced by track_particles().

void tally_image	(	float *	energy,
		float3 *	position,
		float3 *	direction,
		signed char *	scatter_state,
		unsigned long long int *	image,
		struct source_struct *	source_data_SHARED,
		struct detector_struct *	detector_data_SHARED
	)		`[inline]`

Tally a radiographic projection image.

This function is called whenever a particle escapes the voxelized volume. The code checks if the particle would arrive at the detector if it kept moving in a straight line after exiting the voxels (assuming vacuum enclosure). An ideal image formation model is implemented: each pixel counts the total energy of the x rays that enter the pixel (100% detection efficiency for any energy). The image due to primaries and different kinds of scatter is tallied separately.

In the GPU, and atomicAdd() function is used to make sure that multiple threads do not update the same pixel at the same time, which would result in a lose of information. Since the atomicAdd function is only available for 'unsigned long long int' data, the float pixel values are scaled by a factor "SCALE_eV" defined in the header file (eg, #define SCALE_eV 10000.0f) and stored as unsigned long long integers in main memory.

WARNING! If the total tallied signal (for all particles) is larger than "1.8e19/SCALE_eV", there will be a bit overflow and the value will be reset to 0 giving bogus results.

WARNING! The detector plane should be located outside the voxels bounding box. However, since the particles are moved outside the bbox in the last step, they could cross the detector plane anyway. If the particles are less than 2.0 cm behind the detector, they are moved back and detected. Therefore the detector can be a few cm inside the bbox and still work. If the Woodcock mean free path is larger than the distance from the bbox to the detector, we may lose some particles behind the detector!

Parameters:

[in]	energy	X-ray energy
[in]	position	Particle position
[in]	direction	Particle direction (cosine vectors)
[in]	scatter_state	Flag marking primaries, single Compton, single Rayleigh or multiple scattered radiation
[out]	image	Integer array containing the image, ie, the pixel values (in tenths of meV)

Definition at line 482 of file MC-GPU_kernel_v1.3.cu.

References detector_struct::center, detector_struct::corner_min_rotated_to_Y, source_struct::direction, detector_struct::inv_pixel_size_X, detector_struct::inv_pixel_size_Z, detector_struct::num_pixels, detector_struct::rot_inv, detector_struct::rotation_flag, SCALE_eV, detector_struct::total_num_pixels, int2::x, float3::x, int2::y, float3::y, and float3::z.

Referenced by track_particles().

void tally_materials_dose	(	float *	Edep,
		int *	material,
		ulonglong2 *	materials_dose
	)		`[inline]`

Tally the depose deposited inside each material.

This function is called whenever a particle suffers a Compton or photoelectric interaction. The energy released in each interaction is added and later in the report function the total deposited energy is divided by the total mass of the material in the voxelized object to get the dose. This naturally accounts for multiple densities for voxels with the same material (not all voxels have same mass). Electrons are not transported in MC-GPU and therefore we are approximating that the dose is equal to the KERMA (energy released by the photons alone). This approximation is acceptable when there is electronic equilibrium and when the range of the secondary electrons is shorter than the organ size.

The function uses atomic functions for a thread-safe access to the GPU memory. We can check if this tally was disabled in the input file checking if the array materials_dose was allocated in the GPU (disabled if pointer = NULL).

Parameters:

[in]	Edep	Energy deposited in the interaction
[in]	material	Current material id number
[out]	materials_dose	ulonglong2 array storing the mateials dose [in eV/g] and dose^2 (ie, uncertainty).

Definition at line 1547 of file MC-GPU_kernel_v1.3.cu.

References SCALE_eV, ulonglong2::x, and ulonglong2::y.

Referenced by track_particles().

void tally_voxel_energy_deposition	(	float *	Edep,
		short3 *	voxel_coord,
		ulonglong2 *	voxels_Edep
	)		`[inline]`

Tally the dose deposited in the voxels.

This function is called whenever a particle suffers a Compton or photoelectric interaction. It is not necessary to call this function if the dose tally was disabled in the input file (ie, dose_ROI_x_max_CONST < 0). Electrons are not transported in MC-GPU and therefore we are approximating that the dose is equal to the KERMA (energy released by the photons alone). This approximation is acceptable when there is electronic equilibrium and when the range of the secondary electrons is shorter than the voxel size. Usually the doses will be acceptable for photon energies below 1 MeV. The dose estimates may not be accurate at the interface of low density volumes.

We need to use atomicAdd() in the GPU to prevent that multiple threads update the same voxel at the same time, which would result in a lose of information. This is very improbable when using a large number of voxels but gives troubles with a simple geometries with few voxels (in this case the atomicAdd will slow down the code because threads will update the voxel dose secuentially).

Parameters:

[in]	Edep	Energy deposited in the interaction
[in]	voxel_coord	Voxel coordinates, needed to check if particle located inside the input region of interest (ROI)
[out]	voxels_Edep	ulonglong2 array containing the 3D voxel dose and dose^2 (ie, uncertainty) as unsigned integers scaled by SCALE_eV.

Definition at line 418 of file MC-GPU_kernel_v1.3.cu.

References dose_ROI_x_max_CONST, dose_ROI_x_min_CONST, dose_ROI_y_max_CONST, dose_ROI_y_min_CONST, dose_ROI_z_max_CONST, dose_ROI_z_min_CONST, SCALE_eV, short3::x, ulonglong2::x, short3::y, ulonglong2::y, and short3::z.

Referenced by track_particles().

void track_particles	(	int	history_batch,
		int	histories_per_thread,
		int	num_p,
		int	seed_input,
		unsigned long long int *	image,
		ulonglong2 *	voxels_Edep,
		float2 *	voxel_mat_dens,
		float2 *	mfp_Woodcock_table,
		float3 *	mfp_table_a,
		float3 *	mfp_table_b,
		struct rayleigh_struct *	rayleigh_table,
		struct compton_struct *	compton_table,
		struct detector_struct *	detector_data_array,
		struct source_struct *	source_data_array,
		ulonglong2 *	materials_dose
	)

Initialize the image array, ie, set all pixels to zero Essentially, this function has the same effect as the command: "cutilSafeCall(cudaMemcpy(image_device, image, image_bytes, cudaMemcpyHostToDevice))";.

CUDA performs some initialization work the first time a GPU kernel is called. Therefore, calling a short kernel before the real particle tracking is performed may improve the accuracy of the timing measurements in the relevant kernel.

Parameters:

[in,out]	image	Pointer to the image array.
[in]	pixels_per_image	Number of pixels in the image (ie, elements in the array). Main function to simulate x-ray tracks inside a voxelized geometry. Secondary electrons are not simulated (in photoelectric and Compton events the energy is locally deposited).

The following global variables, in the GPU __constant__ memory are used: voxel_data_CONST, source_energy_data_CONST, detector_data_CONST, mfp_table_data_CONST.

Parameters:

[in]	history_batch	Particle batch number (only used in the CPU version when CUDA is disabled!, the GPU uses the built-in variable threadIdx)
[in]	num_p	Projection number in the CT simulation. This variable defines a specific angle and the corresponding source and detector will be used.
[in]	histories_per_thread	Number of histories to simulate for each call to this function (ie, for GPU thread).
[in]	seed_input	Random number generator seed (the same seed is used to initialize the two MLCGs of RANECU).
[in]	voxel_mat_dens	Pointer to the voxel densities and material vector (the voxelized geometry), stored in GPU glbal memory.
[in]	mfp_Woodcock_table	Two parameter table for the linear interpolation of the Woodcock mean free path (MFP) (stored in GPU global memory).
[in]	mfp_table_a	First element for the linear interpolation of the interaction mean free paths (stored in GPU global memory).
[in]	mfp_table_b	Second element for the linear interpolation of the interaction mean free paths (stored in GPU global memory).
[in]	rayleigh_table	Pointer to the table with the data required by the Rayleigh interaction sampling, stored in GPU global memory.
[in]	compton_table	Pointer to the table with the data required by the Compton interaction sampling, stored in GPU global memory.
[in,out]	image	Pointer to the image vector in the GPU global memory.
[in,out]	dose	Pointer to the array containing the 3D voxel dose (and its uncertainty) in the GPU global memory.

Definition at line 135 of file MC-GPU_kernel_v1.3.cu.

References dose_ROI_x_max_CONST, linear_interp::e0, GCOa(), GRAa(), linear_interp::ide, init_PRNG(), locate_voxel(), MAX_MATERIALS, mfp_table_data_CONST, rayleigh_struct::pmax, ranecu(), ranecu_double(), rotate_double(), source(), tally_image(), tally_materials_dose(), tally_voxel_energy_deposition(), float2::x, float3::x, float2::y, float3::y, and float3::z.

Referenced by main().

void trim_name	(	char *	input_line,
		char *	file_name
	)

Extract a file name from an input text line, trimming the initial blanks, trailing comment (#) and stopping at the first blank (the file name should not contain blanks).

Parameters:

[in]	input_line	Input sentence with blanks and a trailing comment
[out]	file_name	Trimmed file name

Definition at line 1840 of file MC-GPU_v1.3.cu.

Referenced by read_input().

void update_seed_PRNG	(	int	batch_number,
		unsigned long long int	total_histories,
		int *	seed
	)		`[inline]`

Initialize the first seed of the pseudo-random number generator (PRNG) RANECU to a position far away from the previous history (leap frog technique).

This function is equivalent to "init_PRNG" but only updates one of the seeds.

Note that if we use the same seed number to initialize the 2 MLCGs of the PRNG we can only warranty that the first MLCG will be uncorrelated for each value generated by "update_seed_PRNG". There is a tiny chance that the final PRNs will be correlated because the leap frog on the first MLCG will probably go over the repetition cycle of the MLCG, which is much smaller than the full RANECU. But any correlataion is extremely unlikely. Function "init_PRNG" doesn't have this issue.

Parameters:

[in]	batch_number	Elements to skip (eg, MPI thread_number).
[in]	total_histories	Histories to skip.
[in,out]	seed	Initial PRNG seeds; returns the updated seed.

Definition at line 3356 of file MC-GPU_v1.3.cu.

References a1_RANECU, abMODm(), LEAP_DISTANCE, and m1_RANECU.

Referenced by main().