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)

Check also out

Docs

Indices and tables