2.1.

Simplified Diffusion Approximate Model

The optical property of biological tissue is usually described by its absorption coefficient ${\mu }_a$ and reduced scattering coefficient ${\mu }_s^{\prime }$. The diffuse reflectance can be calculated with the MC simulation or the diffusion equation. Equation (1) shows the diffusion approximation for the diffuse reflectance, $R_d$, at a distance $r$ from the light source:²¹

diffusion coefficient $D={[3({\mu }_a+{\mu }_s^{\prime })]}^{-1}$,

effective attenuation ${\mu }_{\mathrm{eff}}=\sqrt{{\mu }_a/D}$,

transport mean free path $z^{\prime }={({\mu }_a+{\mu }_s^{\prime })}^{-1}$,

${\rho }_1=\sqrt{r^2+z^{\prime 2}}$ and

${\rho }_2=\sqrt{r^2+{(z^{\prime }+4D)}^2}$.

With the assumption of the diffusion equation being isotropic scattering and large in transport albedo, the result from Eq. (1) differs from the experimental result and the MC simulation at small $r$.⁷ Previous simulation studies have shown that the incident angle²² and width²³ of the light beam can alter the diffuse reflectance. Nonetheless, Eq. (1) gives a close approximation of $R_d$ in homogeneous medium in cases when $a^{\prime }$ is close to 1. The solving of ${\mu }_{\mathrm{eff}}$ from $R_d$ is commonly known as the inverse problem and the solution to Eq. (1) is discussed in Ref. 24. Since the solution has a high computational complexity, our goal is to simplify the solution to the inverse problem by simplifying Eq. (1). MC simulations show that the slope of log ($R_d$) is steeper as ${\mu }_{\mathrm{eff}}$ increases. It is proposed in Ref. 25, and has been verified by our MC simulations that the time-independent spatial diffuse reflectance of a homogeneous medium can be approximated in the form of:

2.2.

Homogeneous Medium

Figure 1 shows an example of the curve fitting with Eq. (2) for the diffuse reflectance from an MC simulation on a homogeneous medium with ${\mu }_a=0.1{\mathrm{cm}}^{-1}$ and ${\mu }_s^{\prime }=10{\mathrm{cm}}^{-1}$. Equation (2) with the derived constants was found to be:

Fig. 1

Fitting diffuse reflectance using simplified diffusion approximate model.

Figure 1 depicts that Eq. 2(a) fits the diffuse reflectance found by the MC simulation in most of the region, except where $r$ is small, $r<0.3\mathrm{cm}$. We repeated this curve fitting with different pairs of absorption and reduced scattering coefficient values and verified that the simulated diffuse reflectance can be fitted with Eq. (2). Here, we will show how the constants in the diffuse reflectance approximation model for homogeneous medium are derived from the first derivative of log $R_d(r)$. From Eq. (2), taking logarithm and differentiation on both sides, we have:

By letting $y(r)=r\frac{d\ln R_d^{\prime }(r)}{dr}$, one can obtain $k_2$ and $k_3$ via linear regression on $y$ against $r$. The values of $\frac{d\ln R_d^{\prime }(r)}{dr}$ can be found from $\ln R_d(r)$ and a number of adjacent points around $r$. $k_1$ is obtained from Eq. (3) after finding $k_2$ and $k_3$. After repeating the curve fitting with a medium of different absorption and scattering coefficients, we plot the found $k_3$ constant against the effective attenuation coefficient value ${\mu }_{\mathrm{eff}}$ and the found $k_2$ constant against the transport albedo value $a^{\prime }$ in Fig. 2.

Fig. 2

(a) Plot of $k_3$ against ${\mu }_{\mathrm{eff}}$; (b) Plot of $k_2$ against $a^{\prime }$.

In Fig. 2(a), MC simulations were implemented with the value of ${\mu }_s^{\prime }$ from 1 to $20{\mathrm{cm}}^{-1}$, while ${\mu }_a$ varies from 0.1 to $1.0{\mathrm{cm}}^{-1}$. For each pair of ${\mu }_a$ and ${\mu }_s^{\prime }$, $k_2$ and $k_3$ were found by the aforementioned linear regression method with $r>1.5\mathrm{cm}$. The plot of $k_3$ versus ${\mu }_{\mathrm{eff}}$ in Fig. 2(a) reveals that the value of $k_3$ is almost equal to ${\mu }_{\mathrm{eff}}$ with an average difference of 1.3%. Therefore, $k_3$ can be used to estimate ${\mu }_{\mathrm{eff}}$. The value of $k_2$ is found to be within the range of 2.0 and 2.2 when the transport albedo value $a^{\prime }>0.9$. The approximate values of $k_2$ and $k_3$ are comparable with those in Ref. 7, where $k_3={\mu }_{\mathrm{eff}}$ and $k_2=2$ for large $r$.

2.3.

Non-Homogeneous Medium

In reality, biological tissue does not have a constant ${\mu }_{\mathrm{eff}}$ along $r$. Instead of applying Eq. (1) for a non-homogeneous medium, we will exploit Eq. (2) for this purpose because of its simplicity. Here, we will consider a medium which is non-homogeneous along the $x\text{-}y$ plane and homogeneous along the depth ($z$) axis. A point light source is used on the origin and the simulated profile has cylindrical symmetry about the $z$-axis. We have shown, in Sec. 2.2, that the ${\mu }_{\mathrm{eff}}$ can be resolved from the first derivative of logarithm of diffuse reflectance for homogeneous medium. For simplicity’s sake, we will first assume that the transport albedo value $a^{\prime }$ is close to 1 in the non-homogeneous medium and, therefore, $k_2$ is relatively constant along $r$. Secondly, we will adopt the assumption of $k_3(r)={\mu }_{\mathrm{eff}}(r)$. With $k_3$ rewritten as an effective attenuation profile and $k_2$ being a constant, Eq. (4) becomes:

Equation (5) is the simplified diffusion approximation for non-homogeneous medium. It relates the first derivative of diffuse reflectance at a distance $r_0$ from the light source to the effective attenuation ${\mu }_{\mathrm{eff}}$ at $r_0$. The simplified diffusion approximation for diffuse reflectance in Eq. (2) becomes:

In order to validate the assumption of constant $k_2$ and examine how well Eq. (6) will fit the diffuse reflectance of a non-homogeneous medium, several MC simulations on media with different non-uniform effective attenuation profiles were implemented. These simulated profiles were selected to investigate the diffuse reflectance modeling of Eq. (6) for various simulated ${\mu }_{\mathrm{eff}}$ profiles with steep and flat slopes. The simulated spatial diffuse reflectance is then fitted with Eq. (6) and the predefined ${\mu }_{\mathrm{eff}}(r)$ from the simulated profile. The fitting results are shown in Fig. 3.

Fig. 3

Fitting diffuse reflectance using simplified diffusion approximate model for non-homogeneous media.

Figure 3 shows that Eq. (6) fits the simulated diffuse reflectance in four cases, but not in the regions with small $r$ or sharp changes of ${\mu }_{\mathrm{eff}}$. The diffusion process acts like a low-pass filter on the effective attenuation profile and, therefore, fast change of ${\mu }_{\mathrm{eff}}(r)$ cannot be observed in the simulated diffuse reflectance. Nonetheless, it demonstrates that the simulated diffuse reflectance of a non-homogeneous medium can be fitted, showing the possibility of approximating ${\mu }_{\mathrm{eff}}(r)$ from $R_d(r)$ based on Eq. (6). Derived from Eq. (5), the effective attenuation profile can be approximated from the diffuse reflectance curve for large $r(>10z^{\prime })$, as shown in Eq. (7). The transport mean free path $z^{\prime }$ is about 0.1 cm in Fig. 3.

Figure 4 shows the resolved profile ${\mu }_{\mathrm{eff}}^{\prime }(r)$ from the simulated diffuse reflectance $R_d(r)$ as in Fig. 3. The value of $k_2$ is fixed at 2.15, as found in Fig. 2(b), and is assumed to be a constant. On the same graph, the respective simulated profile ${\mu }_{\mathrm{eff}}(r)$ is plotted for comparison.

Fig. 4

Resolved effective attenuation ${\mu }_{\mathrm{eff}}(r)$ profiles from diffuse reflectance $R_d(r)$.

The reconstructed effective attenuation profiles in all four cases are low-pass filtered. The average differences between the reconstructed and input profiles are $-1.60$%, 2.11%, 0.257%, and 5.07% for Fig. 4(a)–4(d), respectively, where $r$ is between 1.0 and 4.0 cm. The small average differences show the feasibility of the reconstruction algorithm. However, the resolved profile is not perfectly fitted because of the low-pass filtering of the diffusion process and $a^{\prime }$ is assumed to be a constant along the radial distance $r$. Apart from the region with small $r$, the biggest difference is found where there is a sharp change in the input profile. The simulation result shows that the value of $k_2$ is dependent on both $a^{\prime }$ and $r$. However, the assumption of constant $k_2$ for cases with limited $r$ ($<4\mathrm{cm}$) and small variation of $a^{\prime }$ (0.95 to 1.00) provides simplicity for the reconstruction algorithm and we can still achieve reconstruction that is close to the actual value. In the following section, we extend the algorithm to a 2-D case and enhance the accuracy of the algorithm by averaging the reconstructed profile from multiple light sources.

3.1.

Algorithm and Simulation

We have built a model to simulate the 2-D diffuse reflectance on a semi-infinite medium that is non-homogeneous on the $x\text{-}y$ plane and homogeneous along the $z$-axis. An example of simulated tissue is shown in Fig. 5(a) with a high absorption crucifixion cross. The simulation and reconstruction procedures are illustrated in Fig. 5. The diffuse reflectance images from eight light sources, which were placed on the boundary of a squared area of $4\mathrm{cm}\times 4\mathrm{cm}$ in size, were simulated by the MC simulation. The resultant 2-D diffuse reflectance images in Fig. 5(b) were then input into our 2-D reconstruction program which treated the image as many 1-D radial lines originating from the point of light source. The outcome was the resolved 2-D effective attenuation profiles as in Fig. 5(c). As shown in Fig. 4, the reconstructed value at a distance smaller than 0.5 cm is invalid and cannot be used. At large $r$, the simulated light intensity is too small, leading to low signal-to-noise ratio. The diffuse reflectance is a decreasing function and, therefore, if a non-decreasing profile along $r$ is found, it will be treated as an invalid input value. These invalid areas are indicated in black in Fig. 5(c). By combining the resolved profiles of the 8 light-sources by weighted averaging, we obtained the resolved effective attenuation of the whole $4\mathrm{cm}\times 4\mathrm{cm}$ area, as illustrated in Fig. 5(d), which resembles the profile in Fig. 5(a).

Fig. 5

2-D diffuse reflectance simulation and effective attenuation profile reconstruction.

The 2-D reconstruction algorithm was repeated on four simulated non-homogeneous semi-infinite tissues of $4\mathrm{cm}\times 4\mathrm{cm}$ in size. These cases were created to see the effects of sharp attenuation change and small object size on the 2-D reconstruction. Figure 6 shows, from top to bottom, the simulation inputs, the resolved effective attenuation profiles and their differences in these four cases. In the input and resolved profiles, the red color on the image represents the area of high effective attenuation, while the blue color represents the area of low effective attenuation. In these simulations, the effective attenuation profile on the $x\text{-}y$ plane are shown in Fig. 6(a)–6(d) while it is homogeneous along the depth ($z$) axis.

Fig. 6

Simulations on 2-D non-homogeneous media. Top (a), (b), (c), (d) simulated profiles on $x\text{-}y$ plane; middle (e), (f), (g), (h) reconstructed effective attenuation profiles; bottom (i), (j), (k), (l) differences between the simulated and reconstructed profile.

In all four cases, the reconstructed 2-D profiles, as illustrated in Fig. 6(e)–6(h), agree well with the corresponding input profiles. The differences between the inputs and the reconstructed profiles are shown in Fig. 6(i)–6(l). Larger differences are found in the regions with abrupt changes in effective attenuation coefficients. It has been observed that the reconstructed profile is low-pass filtered by the diffusion process and therefore, these sharp changes cannot be reconstructed precisely. In Fig. 6(i)–6(l), the error percentage between the input profile and reconstructed profile are shown. In the case of Fig. 6(l), where there is no sharp change in the input profile, the average difference is $-1.77$% with a standard deviation of 10.8%. With rapid changes present in the profile, the average difference rises to around 22% in both Fig. 6(i) and 6(j). Similarly, a large difference is found in Fig. 4 where there is a sharp edge. The noise in the reconstructed profile is mostly attributable to the limited simulation time and can be reduced by increasing the number of light packets being simulated. The largest difference is exhibited where there is a sharp change in effective attenuation. It has also been observed, that the resolved value of effective attenuation is smaller than the simulated input value when the transport albedo value $a^{\prime }$ is small, as discussed in Sec. 2.3.

To examine the effect when tissue is non-homogeneous along the depth axis, a series of simulations is carried out. Figure 7 shows the resolved images when the medium inFig. 6(c) is covered by a thin layer, 0 to 5 mm, of homogeneous tissue with ${\mu }_a=0.1{\mathrm{cm}}^{-1}$ and ${\mu }_s^{\prime }=10{\mathrm{cm}}^{-1}$. The resolved value of effective attenuation becomes closer to that of the top homogeneous layer as the thickness of that layer increases. Nevertheless, it demonstrates that even when a thin homogeneous tissue is placed on top of the non-homogeneous medium, the 2-D reconstruction algorithm can still depict the actual pattern based on the diffuse reflectance distribution. As the thickness of the covered layer increases, the resultant pattern fades out and is similar to other diffuse reflectance resolving methods.

Fig. 7

Resolved effective attenuation profiles with a thin homogeneous layer on top. Thickness of the homogeneous layer: (a) 0 mm, (b) 2 mm, (c) 3.5 mm, and (d) 5 mm.

3.2.

Imaging Instrument

Experiments were carried out to test the reconstruction algorithm for the effective attenuation profile on real biological objects. A commercially available 8-bit monochromic CCD camera was used to capture the 2-D diffuse reflectance images on a flat biological tissue of $4\mathrm{cm}\times 4\mathrm{cm}$ in size. The design of the experiment is illustrated in Fig. 8(a) and will be discussed below.

Fig. 8

(a) Imaging Instrument; (b) LED Panel.

The LED panel, which was composed of 8 pairs of surface-mount LED light sources, was in direct contact with the tissue under test. An 8-bit monochromic CCD camera, Mightex MCE-B013-U, which had an output image resolution of $1280\times 1024\text{pixels}$, was attached to the top. The lens, Computar H0514-MP, attached to the camera had a focal length of 5 mm. The camera was controlled by a personal computer, PC, via USB 2.0 interface. Each pair of LED light sources had two LEDs with central wavelength at 660 and 880 nm respectively, and 3-dB spectral width of approximately 18 nm. The LED panel was controlled by a microcontroller unit through a constant current LED sink driver IC with gain control, Macroblock MBI5170. The PC communicated with the microcontroller via RS232/UART interface to instruct the microcontroller to turn the LEDs on one by one. Since the CCD camera had a relatively low 8-bit gray level, we overcame this limitation by bracketing the exposure time. Five images were captured for each LED, with different exposure values, and they were combined in the PC. The captured images were analyzed by the PC in order for the 2-D effective attenuation profile to be calculated. The design of the LED panel is shown in Fig. 8(b). Eight pairs of 660 and 880-nm surface mount typed LEDs were evenly distributed around the perimeter of the square board. With the effective distance of illumination of each LED being approximately 2 to 3 cm, the LEDs were able to cover the whole area bounded by the square board of the LED panel. The actual number of pixels that covers the $4\mathrm{cm}\times 4\mathrm{cm}$ area was around $760\times 760$, which translates to around 19 pixels for each millimeter.

3.3.

System Calibration Using Phantom

A phantom was constructed as in Fig. 9(a) to validate the proposed reconstruction of the effective attenuation profile using the imaging instrument. A piece of solid soap, measuring $9\mathrm{cm}\times 9\mathrm{cm}\times 4\mathrm{cm}$, was used. By using the spatially resolved method in Ref. 7, the absorption and reduced scattering coefficient at 660 nm for the soap were found to be $0.006$ and $0.80{\mathrm{cm}}^{-1}$, respectively. As soap is made of fatty acid and alkaline, its optical properties are comparable to that of fat tissue of which the absorption and reduced scattering coefficients were found to be around $0.02$ and $1.0{\mathrm{cm}}^{-1}$ at 800 nm.²⁶ A triangular prism was cut off on one side of the soap and the gap was filled with light absorbing mud which has a high absorption coefficient of over $2{\mathrm{cm}}^{-1}$. The aforementioned imaging instrument was placed on the phantom. The captured diffuse reflectance images were analyzed using the proposed algorithm and the reconstructed effective attenuation profile is shown in Fig. 9.

Fig. 9

(a) Phantom structure; (b) resolved effective attenuation profile.

The reconstructed profile in Fig. 9(b) successfully reveals the high absorption area gap. The calculated effective attenuation coefficient value found for the soap was around $0.122{\mathrm{cm}}^{-1}$, which is close to the expected value of $0.120{\mathrm{cm}}^{-1}$. A larger error was found near the tip of the triangular shape due to the low-pass filtering effect of the diffuse process. The result agrees well with the effective attenuation profile of the phantom and the simulation result in Sec. 3.1, thus, demonstrating the feasibility of reconstructing the effective attenuation profile using our imaging instrument. The experiment was repeated and was performed in a darkroom with the room temperature maintained at 25°C throughout the experiment.

3.4.

Preliminary Study on Biological Tissue

Using the effective attenuation imaging instrument as described, we performed a preliminary study on living tissue. Prior to the experiment on non-homogeneous tissue, the reduced scattering coefficient and absorption coefficient of the skin of our volunteer was found to be $1.33$ and $0.09{\mathrm{cm}}^{-1}$, respectively, at 660 nm by using the spatially resolved method in Ref. 7. These numbers translate to an effective attenuation coefficient, ${\mu }_{\mathrm{eff}}$, of $0.62{\mathrm{cm}}^{-1}$. Diffuse reflectance images were captured from the lower leg, using 660-nm LEDs, and were analyzed using the above algorithm. The resultant effective attenuation profile is shown in Fig. 10. The reconstructed effective attenuation profile shows a branch-shaped area with high attenuation corresponding to the veins. The veins are usually near the skin surface. Venous blood contains a high concentration of deoxygenated hemoglobin which leads to a high absorption of red light. The value of ${\mu }_{\mathrm{eff}}$ computed for skin tissue, reperesented by a blue area in Fig. 10, is around $0.64{\mathrm{cm}}^{-1}$, whereas, it is about $2.0{\mathrm{cm}}^{-1}$ for the vein. The ${\mu }_{\mathrm{eff}}$ is close to the aforementioned result obtained from the spatially resolved method. The black color in Fig. 10 indicates that the pixel reading of that particular point is not valid under the reconstruction algorithm, as described in Sec. 3.1. The area in black can be reduced by increasing the number of light sources, optimizing the placement of the light source, or increasing the exposure bracketing per light source. There are some horizontal artifacts on the resultant image caused by the CCD noise that occurs when capturing with high gain factor. Compared to the previous spatially resolved studies, our resolved algorithm is able to provide a highly detailed 2-D effective attenuation profile.¹¹ The reconstruction algorithm takes about 200 sec to resolve the effective attenuation profile on a 2.4-GHz single-core PC, and the resolution for the image is 20 pixels per millimeter. The processing time can be drastically reduced to 3 or 4 sec when the resolution of the image is reduced to 2 pixels per millimeter, a higher resolution compared to that in Ref. 11.

Fig. 10

Resolved effective attenuation profile of skin tissue.