Light transport in tissue by 3D Monte Carlo: Influence of boundary voxelization

Abstract

Monte Carlo (MC) based simulations of photon transport in living tissues have become the “gold standard” technique in biomedical optics. Three-dimensional (3D) voxel-based images are the natural way to represent human (and animal) tissues. It is generally believed that the combination of 3D images and MC based algorithms allows one to produce the most realistic models of photon propagation. In the present work, it is shown that this approach may lead to large errors in the MC data due to the “roughness” of the geometrical boundaries generated by the presence of the voxels. In particular, the computed intensity of the light detected on the tissue surface of a simple cubic tissue phantom may display errors from −80% to 120%. It is also shown that these errors depend in a complex manner on optical and geometrical parameters such as the interoptode distance, scattering coefficient, refractive index, etc. and on the degree of voxelization (“roughness”) of the boundaries. It is concluded that if one wants to perform reliable 3D Monte Carlo simulations on complex geometries, such as human brain, skin or trabecular bone, it is necessary to introduce boundary meshing techniques or other equivalent procedures in the MC code to eliminate the deleterious effect of voxelization.

Introduction

Simulation techniques to describe light transport in tissues are considered to be fundamental tools in biomedical optics. Among the different approaches available, the Monte Carlo (MC) based methods represent a “gold standard”. MC based methods are utilised for wavelengths from the visible to the near infrared (NIR) range, covering in this way a large domain of non-ionizing radiations and thus allowing the investigation of a large set of possible biomedical applications.

In the NIR domain, MC based methods have been successful because of their ability to test the validity of complex analytical algorithms obtained from the radiative transport equation. These analytical algorithms allow one to derive, from the experimental optical signals, quantities of fundamental medical interest such as tissue oxy- and deoxy-haemoglobin concentrations, oxygen saturation [1], tissue blood speed/flow [2] etc. and for this reason testing procedures have great importance. It must be noted that it is also possible to use the MC approach to test image reconstruction methods used in NIR topography and tomography [3], [4], [5], [6], [7]. On the other hand, MC simulations permit one to directly “generate” new numerical algorithms allowing one to solve the so-called “inverse problem”. In fact, intuitively, the synthetic data derived from MC simulations can be utilised to directly “fit” real experimental optical data and thus obtain the wanted parameters [8], [9], [10], [11], [12], [13], [14]. This latter approach is particularly interesting when analytical solutions are too complex or even impossible to obtain. Other areas where MC methods have found value include calculating power deposits in laser treatment of port wine stains [15] modelling photoacoustic imaging [16], laser-doppler perfusion and tissue flurorescence [17].

Finally, besides the former straightforward applications of MC, there is also a need to better understand the relationship between measured optical parameters and the underling tissue physiology. This can be done for example, by analysing the derived optical quantities against suitable theoretical models describing the physiology of the investigated tissue [18]. This approach was originally initiated by the international project known as the Physiome Project [19] and the biological mathematical models have already reached a high level of complexity, see Ref. [20]. For these models, hundreds of physiological parameters can be taken into account, the majority of which will probably not be (easily) measurable in vivo. By combining a MC based model, describing the transport of light in a tissue with complex structure, with a biological model [20] describing the physiology of the tissue, one can produce what we call here a Virtual Optical Human Phantom (VOHP). This includes also by definition the examples given in the previous paragraphs that may be seen as a virtual optical instrument allowing one to obtain optical data, generate or test algorithms or, “understand” the relationship existing between the underlying physiological variables and the experimental optical data through the physiological/biochemical aspects.

If these models are to generate reliable optical data, it is clear that the precision of the reproduction of the anatomical features (“geometry”) defining the VOHP becomes a key factor and as we will see in the present work, an approximated geometry might give significantly incorrect results and thus may not be regarded as a “gold standard”. In practice, the most common way to obtain larger scale anatomical information is to use standard medical images such as MRI, CT, PET, etc. For smaller features, as in the case of the description of human skin for studies of port-wine stains for NIR laser therapy, microscopic techniques such as optical coherence tomography (OCT) are of course more suitable. All of these techniques allow one to obtain personalised numerical 3D images related to a real and unique patient. A common feature of all these 3D images is that they are composed of a set of voxels which means that the boundaries of the anatomical structures (e.g. organs) are not smooth as they are in reality but are represented as irregular surfaces composed of the “corners” of the voxels and this might in principle influence the results of the MC simulations. It has already been shown that in the case of thin slabs representing the human skin, the presence of rough random surfaces at the interfaces (described by a stationary random Gaussian process) can lead to significant differences in the estimation of the tissue optical parameters [21]. Similar problems have also been found with periodic roughness of a sinusoidal type [22]. Thus, there is good evidence to believe that voxelization of human anatomy might also introduce distortion in MC results. Such potential errors may be of considerable importance especially when exact or comparable values are needed. This is the case for example for the theoretical evaluation of the laser power needed for port-wine stains treatment [23], [24]. Moreover, if the size of the errors varies for example with the spatial position of the light source/detectors, then it may become impossible to test/reproduce analytical models.

The problem of voxelization is usually neglected in biomedical optics or considered not important especially in the case where the different indexes of refraction of the tissues are set to the same value (i.e. no reflection/refraction is assumed or even made equal to the external medium) [15], [23], [25], [26], [27] Moreover, even if the study of the Fresnel reflectance of rough boundaries is a fruitful domain of research [28], to our knowledge, there are no studies explicitly describing the voxelization problem. In practice, the size of the voxels is usually large compared to the utilized wavelength (especially in the NIR domain) and thus the boundaries appear practically smooth to the incident light from the local point of view and this is considered not to be an interesting problem because this can be described with the usual reflection/refraction laws.

The aim of the present work was to investigate the influence of voxelization on the optical quantities obtained by the MC simulations and thus on the applications cited above. To be more precise this work does not seek to discuss the influence of the roughness of the surface of a “real” tissue/organ, but to consider the roughness generated by the 3D imaging systems with a finite precision and represent images as 3D voxels.

For simplicity, because of the almost infinite number of possible simulations, in this paper we only consider the case of CW radiation and a very simple geometrical form, where for comparison an exact solution can also be obtained (i.e. without the influence of the voxelization). This has the advantage of clearly identifying the different effects of the voxelization.