TLBO-based Algorithms for Minimization of Multi-Ray Path Lengths in Volumetric Objects on the GPU

This master thesis presents an approach in the area of ray tracing to compute and optimize rays in heterogeneous media according to a specific definition of "optimal path", using an appropriate beam model, while also adhering to potential, additional constraints. This was done utilizing the GPU to parallelize the principle steps. The system was then applied to the application area of proton radiation therapy for beam angle optimization problems, which is an important building block for inverse treatment planning.


This master thesis presents an approach and implementation to solve a subproblem in the area of ray tracing in heterogeneous media: that of determining the "shortest path", where the definition of "shortest" depends on the application domain. For example, in nuclear physics, "shortest" usually means the least amount of energy required for a particle/wave to traverse a volume from point A to B. The application domain chosen for this thesis is radiation therapy with protons and how to determine the best beam angles depending on constraints and the media given.

First, we trace rays through a CT (using the parallelization capabilities of the GPU), or similar discretized spatial data in form of a voxel grid using the appropriate beam model to build a solution space called a "cost map". Then, from a voxel of interest in the interior (e.g., a tumor site), rays are traced to each and every surface voxel, accumulating (partial) densities according to the beam model along the ray paths. For memory-saving purposes, the finished cost map is then "hollowed out" by only considering the hull voxels (holding the accumulated values), essentially forming a direct sum of the cuboid face vectors holding the respective face voxels ("cuboid surface" data structure).

Second, the cuboid surface object (storing profit values), together with the domain-specific constraints, is then used to form a binary multidimensional knapsack problem (0/1-MKP). To solve it, a relatively new meta-heuristic called "teaching-learning-based optimization" (TLBO) is applied over a initially generated population of solution candidates. The population generation utilizes that an isomorphism between natural numbers and combinations (of subsets of a set of natural numbers) can be established, which in turn allows combinations to be lexicographically indexed. Using this together with an O(1)-method to generate n-random numbers from a range of numbers allows for rapid solution candidate generation.
A solution is encoded as a vector, of size akin to the number of elements the Cuboid Surface object holds, of binary elements. A 1-component at a specific index "selects" the corresponding profit value. The evaluation operation has been parallelized using the GPU, where for each constraint a grid with a kernel and a given constraint function is launched. The kernel itself consists of a modified inner-product kernel and a modified reduction kernel to allow for access to a function table according to the constraint queued.
Solution candidates are then combined and selected according to the phase and combination operation within TLBO. If a new solution candidate is deemed unfeasible, it gets repaired by the repair operator. This operator consists of two general sub-operations: Drop and add. The drop operation drops, from right to left, a 1-component until the solution is deemed feasible. The add operation then adds, from left to right, a 1-component until the maximum amount of items (number of beams) are reached or the solution is not feasible anymore. To be more effective in adding and dropping elements, the indices of the items in a solution are sorted by their utility ratio, which determines the overall "impact" of an element according to its weight in the calculation of the overall objective score.
All of this is attempted to be solved with performance in mind: Major components of every step ("cost map generation", "initial population generation", "evaluation operation", "repair operation (with drop and add)") are utilizing the GPU.


Component Timings

Cost Map: Average timings over 10 runs.
Initial Population Generation: Minimum, maximum and average timings are taken from 50 runs.
TLBO Main Loop: Minimum, maximum and average timings are taken from 50 iterations of the main loop.

Dose-Volume Histograms

For each scenario, three beams with angles as suggested by the system were applied using matRad.


Full version of the master thesis (English only)


This original work is copyright by University of Bremen.
Any software of this work is covered by the European Union Public Licence v1.2. To view a copy of this license, visit
Any other assets (3D models, movies, documents, etc.) are covered by the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of this license, visit
If you use any of the assets or software to produce a publication, then you must give credit and put a reference in your publication.
If you would like to use our software in proprietary software, you can obtain an exception from the above license (aka. dual licensing). Please contact zach at cs.uni-bremen dot de.