MC-GPU
Functions
MC-GPU_v1.3.cu File Reference

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 Main program of MC-GPU: initialize the simulation enviroment, launch the GPU kernels that perform the x ray transport and report the final results.
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.
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 (#).
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_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.
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).
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.
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.
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.
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 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.
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 IRND0 (float *W, float *F, short int *K, int N)
 Initialisation of Walker's aliasing algorithm for random sampling from discrete probability distributions.

Detailed Description

Author:
Andreu Badal (Andreu.Badal-Soler@fda.hhs.gov)
Date:
2012/12/12 -- MC-GPU v.1.3: 2012/12/12 -- MC-GPU v.1.2: 2011/10/25 -- MC-GPU v.1.1: 2010/06/25 -- MC-GPU v.1.0: 2009/03/17

Definition in file MC-GPU_v1.3.cu.

Function Documentation

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]numCharacters to read
[in]file_ptrPointer to the input file stream
[out]trimmed_lineTrimmed 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 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_espcFile containing the energy spectrum (lower energy value in each bin and its emission probability).
[in,out]source_energy_dataEnergy spectrum and other source data. The Walker alias and cutoffs are initialized in this function.
[out]mean_energy_spectrumMean 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 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_materialsArray with the names of the material files.
[in]density_maxmaximum density in the geometry (needed to set Woodcock trick)
[out]density_nominalArray with the nominal density of the materials read
[out]mfp_table_dataConstant values for the linear interpolation
[out]mfp_table_a_ptrFirst element for the linear interpolation.
[out]mfp_table_b_ptrSecond 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_voxelsName of the voxelized geometry file.
[out]density_maxArray with the maximum density for each material in the voxels.
[out]voxel_dataPointer to a structure containing the voxel number and size.
[out]voxel_mat_dens_ptrPointer to the vector with the voxel materials and densities.
[in]dose_ROI_x/y/z_maxSize 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 main ( int  argc,
char **  argv 
)

Main program of MC-GPU: initialize the simulation enviroment, launch the GPU kernels that perform the x ray transport and report the final results.

This function reads the description of the simulation from an external file given in the command line. This input file defines the number of particles to simulate, the characteristics of the x-ray source and the detector, the number and spacing of the projections (if simulating a CT), the location of the material files containing the interaction mean free paths, and the location of the voxelized geometry file.

Author:
Andreu Badal

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

References source_struct::direction, 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, linear_interp::e0, source_energy_struct::espc, linear_interp::ide, init_energy_spectrum(), voxel_struct::inv_voxel_size, load_material(), load_voxels(), MASTER_THREAD, MAX_MATERIALS, mfp_table_data_CONST, source_energy_struct::num_bins_espc, detector_struct::num_pixels, linear_interp::num_values, voxel_struct::num_voxels, source_struct::position, RAD2DEG, read_input(), report_image(), report_materials_dose(), report_voxels_dose(), detector_struct::sdd, set_CT_trajectory(), source_energy_data_CONST, track_particles(), update_seed_PRNG(), voxel_data_CONST, int2::x, int3::x, float3::x, ulonglong2::x, int2::y, int3::y, float3::y, ulonglong2::y, int3::z, and float3::z.

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]argcCommand line parameters
[in]argvCommand line parameters: name of input file
[out]total_historiesTotal number of particles to simulate
[out]seed_inputInput random number generator seed
[out]num_threads_per_blockNumber 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_outputFile where tallied image is reported
[in]detector_dataDetector description read from the input file (pointer to detector_struct)
[in]imageTallied image (in meV per pixel)
[in]time_elapsedTime elapsed during the main loop execution (in seconds)
[in]total_historiesTotal 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_projectionsNumber of projections simulated
[in]total_historiesTotal number of x-rays simulated per projection
[out]density_nominalArray with the nominal densities of materials given in the input file; -1 for materials not defined. Used to report only defined materials.
[in]materials_doseTallied 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_outputFile where tallied image is reported
[in]detector_dataDetector description read from the input file (pointer to detector_struct)
[in]imageTallied image (in meV per pixel)
[in]time_elapsedTime elapsed during the main loop execution (in seconds)
[in]total_historiesTotal number of x-rays simulated
[in]source_dataData 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().

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 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_lineInput sentence with blanks and a trailing comment
[out]file_nameTrimmed 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_numberElements to skip (eg, MPI thread_number).
[in]total_historiesHistories to skip.
[in,out]seedInitial 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().