Fast and precise image generation of blood vessels embedded in skin

Abstract. A software for fast rendering the visual appearance of a blood vessel located in human skin was developed based on a numerical solution of the radiative transfer equation. The user can specify geometrical properties, such as the depth and the diameter of the vessel, and physiological properties, such as the oxygen saturation of the vessel or the blood concentration in the skin. From these data, the spatially and spectrally resolved reflectance from the skin containing the blood vessel is calculated via Monte Carlo simulations, by which a two-dimensional image is generated. The short calculation time of about a second is achieved by precalculating and storing the spatially resolved reflectance for a variety of combinations of the optical and geometrical properties. This concept gives the user the opportunity to rapidly explore the influence of the physiological and geometrical properties of the investigated blood vessel on its visual appearance. The correctness of the lookup table was validated by comparison with independent Monte Carlo simulations. Rendering examples of different blood vessels in human skins are given. The current version of the software can be downloaded at https://www.ilm-ulm.de/software.


Introduction
For a long time, color phenomena in nature, such as the changing colors of the sky from sunset to twilight, have been fascinating and inspiring people. Knowing the basic physics of the involved effects, such as the single or multiple scattering of light by particles or inhomogeneities, allows us to understand these phenomena. 1 About 20 years ago, the physical explanation for the bluish appearance of veins, besides the fact that blood itself is red, was quantitatively given. 2 It was found that the spectral characteristics of light propagation in combination with Land's retinex theory can be used to explain the color perception of blood vessels. Interestingly, the reason for the bluish color of veins is not a larger reflection of blue light compared to red light from the skin above the vein. It rather is the larger decrease of the reflection for the red light than for the blue light from the skin above the vessel compared to the surrounding skin. Red light has a larger penetration depth compared to blue light in the involved tissues and is, thus, affected to a larger extent by a sufficiently deep vessel. Furthermore, it was calculated that not only veins but also arteries can have a bluish appearance for certain conditions. 2 In general, however, arteries look less bluish compared to veins having the same radius and depth due to the decreased absorption of oxygenated hemoglobin compared to deoxygenated hemoglobin in the red wavelength range. Recently, it was postulated that Rayleigh scattering in the skin has a major impact on the bluish appearance of blood vessels. 3 Studies of other authors focused on related topics, for example, the multiple influences on the color appearance of skin, such as melanin concentration as well as different melanin types, the oiliness of skin or the blood concentration in skin. 4,5 Besides the interests to render skin correctly, diagnostics methods, such as the evaluation of erythema for different skin types, efficient port wine stain treatment, or forensic applications, where the colors of bruises are studied, also benefit from the understanding of light propagation in human skin. [6][7][8][9] Researchers have applied different techniques to calculate the light propagation in biological tissue. One of the most commonly used methods, especially in the field of computer graphics, is the so called diffusion approximation. [10][11][12] The diffusion approximation has its benefits in the easy implementation and high computational speed, but it lacks in accuracy for certain circumstances, for example, if low scattering and high absorption is present. If a higher precision is needed, the radiative transfer equation (RTE) has to be solved exactly. Analytical solutions have been found for certain geometries, such as layered media with and without fluorescence. 13 Unfortunately, the RTE cannot be solved analytically yet for more complex geometries, such as a vessel embedded in a layered skin model. The Monte Carlo (MC) simulation as an, in principle, exact numerical solution of the RTE has been established as the gold standard to calculate the light propagation for such problems. 14,15 A big disadvantage of the MC method is the potentially enormous computational time, due to its statistical nature. To reduce this handicap, massive parallelization using graphic cards (GPU), as well as multiple variance reduction methods, have been developed. [16][17][18][19] In this study, a software is presented to render the visual appearance of a blood vessel embedded in skin. The skin is modeled by two layers representing the epidermis and the dermis. A MATLAB (Mathworks Inc., Natick, Massachusetts) interface is provided, where multiple physiological and geometrical parameters, such as the oxygen saturation of the vessel or the radius of the vessel, can be specified by the user. The image generation is based on the spectrally and spatially resolved reflectance above the vessel, within the visible range of light, which is converted to an RGB-image using the CIE-XYZ color space. 20 The presented software can calculate the spatially and spectrally resolved reflectance, and hence, RGB-image in roughly one second using a lookup table (LUT). To accelerate the generation of the LUT, a variance reduction method was implemented. In the second part of this study, exemplary images are generated by the proposed method. Hereby, the results found by Kienle et al. 2 could be confirmed, such as the fact that a small vessel close to the surface will appear red but can look bluish by increasing its diameter, or the finding that veins are more likely to appear bluish compared to arteries. Beyond that, further interesting cases are investigated. First, the influence of the radius and depth for veins and arteries is systematically investigated. Second, the appearance of superficial vessels depending on the radius and the oxygen saturation of the blood is discussed. Finally, the influence of different spectral dependencies of the reduced scattering coefficient are presented.

Method
The visual appearance of a blood vessel depends on multiple properties, such as its depth, its radius, or the absorption and scattering coefficients of the tissue close to the vessel and of the vessel itself. To predict the appearance of a blood vessel and the surrounding skin, the spatially and spectrally resolved reflectance for the skin surface near the vessel is needed in the visible wavelength range between 400 and 800 nm. In this section, the generation of an efficient and accurate LUT, delivering the spectrally and spatially resolved reflectance I srr ðλ; iÞ detected at the skin surface along a vessel, is described. A variance reduction method was applied to accelerate the generation of the LUT. Curve fitting was performed to reduce the amount of data stored in the LUT. A MATLAB user interface is provided, where the operator can set multiple geometrical and physiological properties. These properties are converted to an image in sRGB color space using the generated LUT. Thus, the appearance of the blood vessel can be rendered in a two-dimensional image.

LUT Generation
A graphic card-based Monte Carlo (MC) simulation, as presented previously, 21 was used to generate an LUT, storing the spatially resolved reflectance I srr ðiÞ for different scenarios of a three-dimensional model as shown in Fig. 1. Note that the model is symmetric along its y axis. The epidermal layer was assumed to be a slab, the dermis to be semi-infinitely thick, and the blood vessel to be extended to infinity along the y axis. The blood vessel's x-position was set to be x ¼ 0. Photons were started homogeneously on a line from x ¼ −0.5l s to x ¼ 0.5l s at y ¼ 0 in z direction where l s is the light source size. Reflected photons were recorded with respect to their x-position within a certain range. The light source size l s was set three times as long as the detection range to approach an infinitely extended light source, such that problems at the detector-boundaries were avoided. For the following simulations, l s was fixed at l s ¼ 45 mm and the detection bin size was kept constant at jx iþ1 − x i j ¼ 0.1 mm. The spatially resolved reflectance E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 6 3 ; 1 0 5 I srr ðiÞ ¼ can be derived, where PðiÞ is the number of photons received between x i and x iþ1 , P is the total number of photons started, and N is the number of spatial bins. The factor 3 is due to the three times larger light source size compared with the detection range.
Even though a variance reduction method was used, as described in Sec. 2.2, only a limited number of dimensions and corresponding grid points could be chosen for the LUT, to keep the computational time as well as the memory demand of the final LUT in an acceptable range. The LUT has six-dimensions representing properties that have a major influence on the color impression: 1. Absorption coefficient of the epidermis μ a;1 in mm −1 .
3. Absorption coefficient of the blood vessel μ a;3 in mm −1 .
4. Reduced scattering coefficient of the dermis μ 0 s in mm −1 . The same reduced scattering coefficient was used for the blood vessel. A two times higher reduced scattering coefficient was applied for the epidermis. 4,22 5. Depth of the blood vessel z in mm. The depth is defined as the shortest distance of the blood vessel (top of the vessel) to the surface. The minimal depth is limited by the thickness of the epidermis since the vessel is always placed in the dermis.
6. Radius of the blood vessel r in mm.
The used ranges of the dimensions as well as the corresponding number of grid points are listed in Table 1. The grid points for z and r were spaced linearly across the corresponding ranges. The grid points for μ 0 s were also spaced linearly, but a finer linear spacing for μ 0 s < 1 mm −1 than μ 0 s > 1 mm −1 was chosen. Values for the absorption coefficients μ a;m , m ¼ 1;2; 3 were spaced logarithmically such that the difference between two grid points is smaller for small μ a;m and is increasing for larger μ a;m .
The values for the fixed parameters can be found in Table 2.

Scaled Monte Carlo
A variance reduction method, the so called scaled MC simulation, was used to reduce the computational time when the LUT is created. 16,19 Implementing this method allows us to generate I srr ðiÞ for arbitrary combinations of μ a;m with a single MC simulation for given z, r, and μ 0 s . At first, an MC simulation with rather small absorptions for each tissue type, namelyμ a;1 ¼ 10 −4 mm −1 ,μ a;2 ¼ 10 −4 mm −1 , andμ a;3 ¼ 10 −2 mm −1 , was performed. During the simulation the path lengths s m a photon spent in the corresponding tissues m ¼ 1;2; 3 were saved for each photon. Then

Curve Fitting
Saving the whole spatially resolved reflectance with the reported resolution would result in roughly 10 gigabytes of memory for the LUT in single precision. To reduce the memory demand as well as to increase the interpolation speed, the function E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 6 3 ; 1 5 7 fðiÞ ¼ a þ was fitted on the spatially resolved reflectance I srr ðiÞ, generated by the scaled MC simulation [see Eq. (2)]. The constraint 5 can be applied since a dip of I srr ðiÞ at x ¼ 0 mm is expected, if the absorption of the blood vessel μ a;3 is larger than the absorption of the dermis μ a;2 . In total 2K þ 1 individual parameters were stored within the LUT for one spatially resolved reflectance. Curve fitting was performed in C++ using the library ceres. 24 It was found that K ¼ 3 results in a good trade-off between accuracy and memory demand. In bðlÞ þ bðkÞ ≈ const: (9) lead to the same, correct, fit result. To improve the robustness of the interpolation, a postprocessing step was performed if cðlÞ ≈ cðkÞ where E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 3 2 6 ; 2 3 1 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 3 2 6 ; 1 9 9 bðlÞ ¼ bðkÞ ¼ 0.5ðbðlÞ þ bðkÞÞ:

RGB-Image Generation
Based on the wavelength-dependent optical properties μ 0 s ðλÞ and μ a;m ðλÞ, the visual appearance, in terms of RGB-values, of a blood vessel and the surrounding tissue can be displayed on a screen. At first a spectrally and spatially resolved reflectance I srr ðλ; iÞ is delivered by the LUT. Hereby linear, n-dimensional interpolation is applied separately for all 2K þ 1 parameters for each wavelength by the MATLAB function interpn (Mathworks Table 1 Value ranges and the corresponding number of grid points for each dimension of the LUT.

Variable
Range # grid points where SðλÞ is the spectrum of the incident light and xðλÞ, yðλÞ, and zðλÞ are the CIE standard observer color matching functions. 20,25 Throughout this study 5-nm steps were used for the calculations. RGB-values

User Interface
Typically, the six-dimensions of the LUT are rather unintuitive to discuss the visual appearance of a blood vessel. It is more convenient to use properties that are directly linked to the physiology of skin, such as melanin concentration of the epidermis or oxygen saturation of the blood vessel. Therefore, a user interface was designed, where such physiological properties can be specified. Similar to the existing literature, some assumptions about the chromophores of the skin model were made. 4,22,[28][29][30] It was assumed that the main absorber present in the epidermis is melanin, the main absorber inside the dermis is hemoglobin, which also is the only absorber for the blood vessel. The absorptions of the dermis and epidermis were extended by a baselineabsorption.
The following eight parameters p 1 − p 8 can be entered by the user in the current version of the MATLAB interface. The reduced scattering coefficient of the dermis is set to be always half of μ 0 s of the epidermis.
• p 7 specifies the depth of the blood vessel: z ¼ p 7 mm.
• p 8 specifies the radius of the blood vessel: r ¼ p 8 mm.
Note that the proposed equations can easily be changed within the provided MATLAB script as desired since the LUT is generated completely independent of the parameters p 1 − p 6 . An example of the user interface can be seen in Fig. 3 with fp 1 ; : : : ; p 8 g ¼ f0.66; 0.9; 0.005; 0.006; 0.80; 0.60; 0.64; 0.55g.

LUT Quality Test
To verify the LUT and especially to check the quality of the n-dimensional interpolation, 450 sample sets fz; r; μ 0 s ; μ a;m g were randomly drawn among the allowed six-dimensional parameter space as listed in Table 1, such that all combinations in between the grid points could be chosen. For each set a regular, independent MC simulation was performed with the exact six parameters and the spatially resolved reflectance I MC srr was recorded using Eq. (1). Next, the LUT was used to generate the 450-corresponding data I LUT srr . The mean relative difference E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 4 ; 3 2 6 ; 1 4 9 jI MC srr ðiÞ − I LUT srr ðiÞj I MC srr ðiÞ was computed for all sample sets. A histogram of the δ I distribution, considering all 450 sample sets, can be found in Fig. 4(a). In Fig. 4(b), the δ I distribution for all sample sets with μ 0 s > 0.5 mm −1 is shown. For small scattering coefficients, the interpolation may occasionally lead to relative errors up to 9%, but only 4 out of 450 sample sets had a relative error of δ I > 2%. The four outliers were the results for combinations of high absorption coefficients with small reduced scattering coefficients of the dermis, causing I srr to be smaller than 0.05, which is an irrelevant range of the reflected intensity for this study. Overall, a very good agreement of the spatially resolved reflectance calculated with a separate MC simulation and the corresponding reflectance obtained by the LUT was found for the vast majority of the parameter range.

Sample Images
In this section, a variety of influences on the visual appearance of a blood vessel are presented using the method introduced before. Further investigations of interest can easily be performed by each user individually after downloading the LUT and the user interface.

Depth and Radius of Blood Vessels
A major influence on the appearance of a blood vessel have its geometrical properties, namely its depth z and radius r, especially to elaborate on the question, why veins appear blue although blood itself is red. First, nine images were generated, where all combinations of three distinct depths and three radii were varied. The results for veins with an oxygen saturation of 0.6 for the vessel are shown in Fig. 5 and the images for arteries with an oxygen saturation of 1.0 for the vessel can be found in Fig. 6. Figure 5 shows that for the bluish appearance of a vein a sufficient depth and radius are necessary, compare Figs. 5(e), 5(f), 5(h), and 5(i). This behavior can be explained by the fact that the reduced scattering coefficient and, more pronounced, the absorption coefficient is larger for the blue part of the visible light compared to the red part. First of all, the increasing fraction μ 0 s ðλÞ μ a ðλÞ for larger wavelengths causes the typical light reddish color of the skin, for areas without a vessel, where more red light is reflected compared to blue light. To understand the bluish appearance of veins, it is important to notice that the spectral course of the optical properties lead to a larger penetration depth of red light and thereby result in a higher probability to reach the vessel compared with blue light. The vein has a high absorption, even for red light, compared to skin and thus absorbs a large part of the light incident onto the vein. Combining these observations induces that the reflection of red light is affected to a larger extend by the vessel compared to blue light in a way that a dip of the fraction RðiÞ BðiÞ is formed for positions i above the vessel. Although the reflection in the red RðiÞ above the vessel might be larger then the reflection in the blue BðiÞ above the blood vessel, a vein typically appears bluish according to Land's retinex theory due to the dip of the fraction RðiÞ BðiÞ , see for details Refs. 2 and 32. Blood vessels that are very close to the surface and have a small radius, such as Figs. 5(a) and 6(a), appear reddish. An explanation is given in Sec. 4.2, where superficial vessels are investigated in detail. For arteries, it can be seen that the bluish appearance is not that pronounced but still present in Figs. 6(h) and 6(i). In general, the arteries appear less bluish compared to veins because the absorption in the vessel is decreasing in the red wavelength range when the oxygen saturation is increased. On top of the reddish appearance, the decreasing absorption for red light is also a reason for the overall worse visibility of arteries compared to veins. For arteries as well as veins, vessels with a rather small radius become almost invisible for larger depths, see Figs. 5(b) and 5(c) as well as Figs. 6(b) and 6(c).

Oxygen Saturation and Radius of Superficial Vessels
Next, superficial blood vessels are investigated with a fixed depth at z ¼ 0.1 mm. The oxygen saturation as well as the radius of the vessels were varied in three steps. The resulting images of the nine combinations can be found in Fig. 7. As described previously, the color of a vessel can turn from a bluish to a reddish appearance, when the oxygen saturation is increased. Superficial vessels appear darker by increasing the radius because more light is absorbed in total. In addition, the amount of red light that is absorbed increases at a faster rate compared to blue light, when the radius is increased, due to the larger penetration depth of red light and the color is eventually becoming more and more purple for arteries and bluish for veins. For superficial vessels with a small diameter, as shown in Figs. 7(a)-7(c), the vessels always appear reddish, independent of the oxygen saturation because major parts of the blue light penetrate to the vessel, where they are strongly absorbed. Contrarily, red light is not as strongly absorbed and has a higher chance to transmit through the vessel and ultimately reflect from the tissue.

Spectral Dependency of the Reduced Scattering Coefficient
To evaluate the influence of the spectral dependency of the reduced scattering coefficient on the appearance of the blood vessel, different equations found in literature were applied for μ 0 s :  • Verdel18: 22 μ 0 s ¼ 5.5ð λ 500 nm Þ −1.5 mm −1 • Leeuwen18: 3 μ 0 s ¼ 5.5ð λ 500 nm Þ −4 mm −1 which are plotted in Fig. 8. "Donner06" was used throughout this study. The equation "Leeuwen18" is not explicitly given in their publication. Here, the spectral dependency was kept at λ −4 and the magnitude was adjusted. The corresponding images can be found in Fig. 9. The spectral dependency of μ 0 s ðλÞ found by Verdel et al. has only a λ −1.5 dependency, contrarily to the first two, which are given by a combination of a λ −4 dependency with a λ −0.22 dependency and a λ −1.5 dependency, respectively. Therefore, the reduced scattering coefficient "Verdel18" has the smallest difference between the red and blue wavelength range, resulting in an overall lighter reddish appearance of the skin compared to the other images. In addition, the blood vessel appears more blurred since the reduced scattering coefficient "Verdel18" is the largest for the red part of the spectrum. The reduced scattering coefficients "Jacques98" and "Leeuwen18" show the biggest difference between the red and blue wavelength ranges, resulting in an overall bluish touch of the skin. Leeuwen et al. postulated that the Rayleigh scattering of skin has a major influence on the typical bluish color of veins. 3 By comparing Figs. 9(c) and 9(d), one can clearly see that even a λ −1.5 dependency of μ 0 s ðλÞ leads to a bluish appearance of the vessel. Although, it was stated that the λ −4 dependency of μ 0 s ðλÞ is the decisive reason that leads to a significantly higher reflection for the blue channel B over the red channel R above the vessel. However, as Figs. 9(a)-9(d) show, this is not important, because the vessel appears bluish in all images.

Model Selection
In principle, skin consists of multiple layers, such as the stratum corneum, the stratum granulosum, dermal papillae, followed by different subcutaneous tissue layers. 23 In addition, a blood vessel is not only a tube but also is a multilayer system considering the vein wall. On the one hand, this study focuses on the major parameters that influence the visual appearance, such as oxygen saturation or depth and radius of the vessel. On the other hand, an efficient way to obtain a spatially resolved spectrum within a second should be achieved. The model was chosen to optimize these two opposing goals.
One of the most debatable restrictions of the used model is the assumption of an infinitely thick dermis layer. In order to investigate the influence of a subcutaneous layer, two MC simulations were performed. The first simulation calculated an image applying the three different tissue types, epidermis, dermis, and blood vessel, as it was used throughout this study. For the second simulation, a subcutaneous layer was added such that the dermis becomes a slab with a thickness of 1.0 mm and the subcutaneous layer was assumed to be semiinfinitely thick. The optical properties for the subcutaneous layer were set to μ a ðλÞ ¼ 0.01 mm −1 and μ 0 s ðλÞ ¼ 1 mm −1 for all wavelengths. 34 For both simulations, the same parameters fp 1 ; : : : ; p 8 g ¼ f1; 0.9; 0.005; 0.006; 0.8; 0.6; 0.1; 0.35g were used and the results are shown in Fig. 10. No major difference of the visual appearance can be seen in Fig. 10(a) when a subcutaneous layer is included (bottom) compared to the original   Fig. 10(b) and the corresponding relative difference RGB−RGB sub RGB in Fig. 10(c), only small alterations were found. Therefore, it can be concluded that a subcutaneous layer has only a minor influence on the visual appearance of the model.

Computational Time
The time needed in MATLAB to generate one image is roughly one second on a single core of a standard computer. A future goal is the usage of the presented LUT to fit a measurement hyperspectrally along a blood vessel. Therefore, an implementation, e.g., in c++ of the multidimensional interpolation might be needed to reduce the fitting time. Another possibility would be the use of neural networks. Previous studies successfully showed the determination of optical properties of tissue using neural networks. 35,36

Conclusion
In this study, a software is presented to render the visual appearance of a blood vessel in skin for a variety of relevant parameters of the vessel and the surrounding tissue. To achieve fast image generation, an LUT was generated that stores the spatially resolved reflectance I srr ðiÞ from a two-layered skin model containing a blood vessel. Six parameters can be chosen arbitrarily, namely the depth and radius of the vessel, the scattering coefficient of the epidermis/dermis, and the absorption coefficients for each of the three materials. To significantly reduce the computational time when generating the LUT, a variance reduction method was applied, where the results for all absorptions are generated from a single MC simulation. To reduce the memory demand as well as to increase the interpolation speed, a function, which is the sum of three Gaussian, was fitted to each I srr ðiÞ resulting in an almost exact representation with seven parameters. By obtaining I srr ðλ; iÞ for λ ¼ 400 nm to 800 nm the visual appearance was calculated using the CIE-XYZ color space.
A MATLAB user interface is provided, which allows the user to enter physiological and geometrical properties of the model, such as oxygen saturation of the vessel or its radius, which are converted and displayed as an sRGB-image. The correctness of the LUT, especially of the n-dimensional interpolation, was proven. The LUT was used to demonstrate a variety of influences on the visual appearance of veins and arteries.
A main goal of this study was to introduce the new software and demonstrate some applications to give an idea of the opportunities. The combination of the fast evaluation as well as the physical correctness granted by more than 16 millions MC simulations makes this framework applicable for future applications, i.e., for diagnostics based on hyperspectral imaging or educational purposes. Since the LUT depends only on basic physical properties, four optical properties, and two geometrical properties, the users can adjust and implement their own equations linking physiological properties to the physical ones. The short calculation time makes the software usable in real-time applications.

Disclosures
The authors declare that there are no conflicts of interest related to this article.