DXMClib¶
DXMClib (Diagnostic X-ray Monte Carlo) is a radiation dose scoring library for diagnostic photon energies in voxelized geometry written in C++. The main goal for this library is to provide an accurate enough physics model to describe and model x-ray sources and estimate radiation doses by the Monte-Carlo method.
If you are looking for an application with graphical user interface to perform simulations, OpenDXMC uses this library as Monte Carlo engine and also allows for import CT images and phantoms as scoring volumes.
Example of usage¶
A brief example on how to use DXMClib to simulate a pencilbeam of 60 keV onto a box of aluminum surrounded by water and air is shown below.
1#include "dxmc/material.hpp"
2#include "dxmc/source.hpp"
3#include "dxmc/transport.hpp"
4#include "dxmc/world.hpp"
5
6#include <array>
7#include <iostream>
8#include <memory>
9#include <numeric>
10
11using namespace dxmc;
12
13template <typename T>
14double calculate()
15{
16 // Lets create a world that describes our voxelized model to score dose in.
17 World<T> world;
18 // We need to specify dimensions and voxel spacing (in millimeters)
19 std::array<std::size_t, 3> dimensions = { 56, 56, 56 };
20 world.setDimensions(dimensions);
21 std::array<T, 3> spacing = { 1.0, 1.0, 1.0 };
22 world.setSpacing(spacing);
23
24 // Now we fill the world box with some materials
25 // We specify three materials
26 Material air("Air, Dry (near sea level)"); // Material air("N0.76O0.23Ar0.01") is equivalent
27 Material water("Water, Liquid"); // Material water("H2O") is equivalent
28 Material aluminium(13); // Material aluminum
29
30 // Lets add the materials to the world material map
31 world.addMaterialToMap(air);
32 world.addMaterialToMap(water);
33 world.addMaterialToMap(aluminium);
34
35 // now we create the density array and material index array
36 auto arraySize = world.size();
37 // the type of density array is double
38 auto densityArray = std::make_shared<std::vector<T>>(arraySize);
39 // the type of material index array is std::uint8_t or old school unsigned char (up to 255 materials are supported)
40 auto materialIndexArray = std::make_shared<std::vector<std::uint8_t>>(arraySize);
41 // Now fill the arrays
42 for (std::size_t i = 0; i < dimensions[0]; ++i)
43 for (std::size_t j = 0; j < dimensions[1]; ++j)
44 for (std::size_t k = 0; k < dimensions[2]; ++k) {
45 auto index = i + j * dimensions[0] + k * dimensions[0] * dimensions[1];
46 if (k < dimensions[2] / 3) { // The first 1/3 part of the box is air.
47 densityArray->data()[index] = air.standardDensity();
48 materialIndexArray->data()[index] = 0;
49 } else if (k < dimensions[2] * 2 / 3) { // The second 1/3 part of the box is water.
50 densityArray->data()[index] = water.standardDensity();
51 materialIndexArray->data()[index] = 1;
52 } else { // The third 1/3 part of the box is aluminum.
53 densityArray->data()[index] = aluminium.standardDensity();
54 materialIndexArray->data()[index] = 2;
55 }
56 }
57 // Add the density array and the material index array to world.
58 world.setDensityArray(densityArray);
59 world.setMaterialIndexArray(materialIndexArray);
60 // The world is now complete, to test if it is valid and all material indices are correct
61 // and arrays have correct dimensions we can call:
62 world.makeValid();
63
64 // We need a beam source, in this case a simple pencil beam
65 PencilSource<T> pen;
66
67 // PencilSource only support monochromatic intensity, in this case we use 60 keV initial photon energy
68 pen.setPhotonEnergy(60.0); // keV
69
70 // Position the source above the world box
71 std::array<T, 3> source_position = { 0, 0, -(dimensions[2] * spacing[2]) };
72 pen.setPosition(source_position);
73 // We want the source plane normal aka the beam direction to point towards the box
74 // To set the direction of the beam we must specify the beam direction cosines.
75 // This is basically two orthonormal vectors on the beam plane [x, y].
76 std::array<T, 6> beam_cosines = { 1, 0, 0, 0, 1, 0 };
77 pen.setDirectionCosines(beam_cosines);
78 // The beam direction is then x cross y
79 std::array<T, 3> beam_direction;
80 vectormath::cross(beam_cosines.data(), beam_direction.data());
81 // The beam direction should be [0, 0, 1]
82
83 // Specify number of exposures, each exposure is run on a single thread
84 // so make sure there are more exposures than cores on your CPU for optimal performance.
85 // For CT sources, number of exposures are determined by the geometry.
86 pen.setTotalExposures(20);
87 // Set total number of histories per exposure
88 pen.setHistoriesPerExposure(1000000);
89
90 // Create a transport object
91 Transport<T> transport;
92 // Do the simulation, this might take some time depending on total number of histories.
93 auto res = transport(world, &pen);
94
95 // Print the dose along z direction of box
96 auto materials = world.materialMap();
97 auto time = std::chrono::duration_cast<std::chrono::milliseconds>(res.simulationTime).count() / 1000.0;
98 std::cout << "Simulation time, " << time << ", seconds\n";
99 std::cout << "Depth [mm], Dose, nEvents, Material index, Material\n";
100 auto dose = res.dose;
101 for (std::size_t k = 0; k < dimensions[2]; ++k) {
102 std::cout << spacing[2] * k << ", ";
103 const auto offset = dimensions[0] * dimensions[1];
104 const auto k_offset = dimensions[0] * dimensions[1] * k;
105 auto sum_start = dose.begin() + k_offset;
106 auto sum_end = dose.begin() + k_offset + offset;
107 auto sum = std::accumulate(sum_start, sum_end, 0.0);
108 std::cout << sum / static_cast<T>(offset) << ", ";
109 auto n_events = std::accumulate(res.nEvents.cbegin() + k_offset, res.nEvents.cbegin() + k_offset + offset, 0);
110 std::cout << n_events / static_cast<T>(offset) << ", ";
111 auto material_index = static_cast<std::size_t>(materialIndexArray->data()[k_offset]);
112 std::cout << material_index << ", " << materials[material_index].prettyName() << std::endl;
113 }
114 return time;
115}
116
117int main(int argc, char* argv[])
118{
119 auto n_secs_d = calculate<double>();
120 auto n_secs_f = calculate<float>();
121 std::cout << "Simulation time for double precision: ";
122 std::cout << n_secs_d << " seconds\n";
123 std::cout << "Simulation time for single precision: ";
124 std::cout << n_secs_f << " seconds\n";
125 std::cout << "Single precision was " << std::ceil(100 - n_secs_f / n_secs_d * 100);
126 std::cout << " percent faster\n";
127 return 1;
128}
The examble can be built by CMake with the following CMakeLists.txt file.
1# testing to see if dxmclib target is defined. If it's not we will download
2# dxmclib at configure time.
3IF( NOT TARGET libdxmc)
4 include(FetchContent)
5 ## Adding DXMClib package
6 FetchContent_Declare(
7 libdxmc
8 GIT_REPOSITORY https://github.com/medicalphysics/DXMClib.git
9 GIT_TAG master
10 )
11 FetchContent_MakeAvailable(libdxmc)
12ENDIF()
13
14
15add_executable(pencilbeam pencilbeam.cpp)
16target_link_libraries(pencilbeam PRIVATE libdxmc)
17
18install(TARGETS pencilbeam
19 RUNTIME DESTINATION examples)