1.

Introduction

The requirement of noninvasive techniques to study tissues has led to the development of several methods. Optical techniques show lot of potential and hence many medical diagnostic tools have been developed with them. In the field of optical imaging, techniques like optical coherence tomography, birefringence imaging, and confocal microscopy are under development.^{1, 2, 3, 4, 5} All these techniques provide images of subsurface tissues but with poor contrast as depth increases. Polarization discrimination techniques have been proposed.^{6, 7, 8} The results show that the visibility depth could be extended, by observing in the orthogonal or opposite polarized state when the illumination is linearly or circularly polarized. Spectral polarization difference imaging⁹ was proposed as an imaging method to image subsurface tissues at larger depths but has its drawbacks since it requires ultrafast lasers and complex procedures. The effectiveness of the spectral and polarization filtering technique¹⁰ to image tissues at a depth on the order of 1 cm was demonstrated recently. This method proposed that if the illumination of different wavelengths was used the number of photons reaching a specified depth without scattering increased as wavelength increased. But still the images obtained lacked significant contrast. In this paper, imaging depths of 2.5 cm are demonstrated by introducing an improvement to the spectral and polarization filtering technique.

The polarization discrimination technique is based on the assumption that weakly scattered light retains its initial polarization whereas highly scattered light does not. Consider a tissue illuminated by a linearly polarized source. A linear polarizer positioned after the light source ensures linear polarized illumination. As polarized photons penetrate the tissue, they interact with various tissue structures and some of the injected photons emerge from the tissue in the backscattering direction. Light that reflects from the surface (known as a spectral reflection) retains its polarization but the light that backscatters from somewhere below the surface of the tissue is depolarized. If a second polarizer is positioned in front of the camera lens with its polarization orientation perpendicular to that of the first polarizer, it will record only the perpendicular or orthogonal polarization image components. Therefore, the light viewed in the orthogonal state contains a large proportion of the light that interacted with the object than the copolarized light.

However, even for this orthogonal state the loss of contrast is still the light backscattered by the medium. So contrast could be enhanced if this backscattered contribution from the surrounding medium is removed.¹¹ Hence if the copolarized image, i.e., the image taken by positioning the second polarizer in front of the camera lens in the same orientation as the polarizer placed after the light source, is also recorded and a suitable fraction is subtracted from the orthogonal polarized image, significant contrast could be achieved. Now if spectral filtering, or subtraction of the contrast enhanced images (longer wavelength image from the shorter wavelength one), is done, significant improvement in the contrast of objects at deeper depths is observed.

Figure 1 illustrates the schematic of the proposed imaging system. A white light source along with an optical filter is used to obtain illumination at 632 nm and 990 nm. This is linearly polarized and irradiates the tissue sample. A target object whose image is to be taken is placed below the surface of the tissue sample. By changing the position of the second polarizer fitted to a CCD camera, copolarized and orthogonal polarized images can be taken. These images are subjected to polarization subtraction and wavelength subtraction and the resultant image is obtained.

Fig. 1

Schematic of the imaging system.

2.

Monte Carlo Simulation

The essential characteristic of Monte Carlo code is to study the propagation of light in biological tissue. The basic idea is that photons are launched from an isotropic point source of unit power $P=1\mathrm{W}$ within an infinite homogeneous medium with no boundaries. The medium has optical properties of absorption, scattering, and anisotropy. $N$ photons are launched, each with a “photon weight” initially set to 1. The photon takes steps between interactions with the tissue. The steps are based on the probability of photon movement before interaction by absorption and scattering. During each step as the photon propagates, the photon deposits a fraction of its weight into the local bin at its position. The photon is assumed to be dead if the weight of the photon falls below a fixed threshold. Now a new alive photon is launched.

In this paper this code is slightly modified taking into account the wavelength of the light transmitted and also its polarization state. The model used for the simulation is shown in Fig. 2 .

Fig. 2

Simulation model.

Photons are launched as a beam from the origin that initially propagates parallel to the $z$ axis. The illumination is assumed to be linearly polarized. In polarimetry, the Stokes vector $I$ of a light beam is given by

Each element of the Stokes matrix can be written in terms of the light intensities as

The radiance components parallel to the reference plane, i.e., copolarized and perpendicular to the reference plane, i.e., cross polarized components, can be expressed in terms of the Stokes parameters as

Hence the output Stokes vector varies as the polarization state of the incident beam and the Mueller matrix of the sample.

The Stokes vector for the four incident polarization states is

The far-field solution of Maxwell’s equations is expressed in terms of the two scattering functions¹³

Also,

The Mie scattering matrix is given by

After every scattering event in the tissue, the polarization state and the direction of propagation of the ray change. The scattering angle is chosen randomly according to a probability law depending on the phase function.^{14, 15} The scattering angle $\theta$ is generated by inversion and rejection methods.¹⁶ At each scattering position, the code randomly chooses the azimuthal angle $\psi$ for the rotation of the plane $xy$ . Under the assumption that the number of effective scattering events were larger, $\psi$ was assumed to be equiprobable.¹⁵ Finally performing the local coordinate transformation using the rotation matrix $\mathbf{T}$ ,¹⁴ the modified Stokes vector after scattering by the particle is

Next it is assumed that an object is placed at the boundary of the tissue, i.e., at a distance $z1$ in the $z$ direction. The object is modeled as a depolarizing diffuse reflector. The depolarization caused by the surface is modeled by the Mueller matrix¹¹

The photons after scattering by the tissue medium hit the object and exit the tissue. For each photon backscattered and returning to the top surface, the coordinates $(x_d,y_d)$ at which it intersected the detector plane were computed from¹⁷

A photon is counted as detected if its position falls within the area of the imaging lens, which is assumed to be circular and is determined by

The first two rows of the Stokes vector of these photons impinging on the lens yield the copolarized and the cross-polarized components and these are transformed to the image plane. The image plane was modeled as a grid and yields an image of size $60\times 60\mathrm{pixels}$ .

3.

Results

Simulations were run to investigate the effects of wavelength on the scattering of a photon in a tissue medium. Figure 3 illustrates the effect on the Mie scattering function by the wavelength of the photons for the wavelengths 632 and 990 nm. The scattering function was computed as a function of the scattering angle $\left(\theta \right)$ for the wavelengths and the results are shown. Comparison of the scattering function shows that longer wavelength photons are scattered less than shorter wavelength photons for most scattering angles. This also means that a photon at a longer wavelength, being less scattered, will reach a greater depth than a photon launched from a shorter wavelength source.

Fig. 3

Scattering function versus scattering angle for different wavelengths.

To demonstrate that scattering is less for longer wavelength photons, an experiment was conducted. Polarized light at the wavelengths 632 and 990 nm was sent to varying thicknesses of tissue phantom. A white light fitted with 632 nm and 990 nm bandpass filters and a linear polarizer acted as the source. Linearly polarized light at the wavelengths 632 nm and 990 nm was sent to a tissue phantom whose thickness varied from 0.5 to 2.5 cm. The phantom used for the experiment was $2.02\mu \mathrm{m}$ polystyrene spheres (Polysciences Inc. 19814) suspended in water. The phantom was created by dilution with distilled water to 0.001 spheres per cubic micron. The power of the emanating light was measured using an optical detector and power meter (Advantest, TO8215). Figure 4 illustrates the depth-dependent intensity decay for the wavelengths under investigation. The curves indicate that as thickness of tissue was increased, at 632 nm the power drops considerably, indicating more scatter. But at 990 nm, even for phantom thickness of 2.5 cm the power of the received light is considerable, indicating less scatter and therefore deeper penetration at 990 nm compared to 632 nm. The experiment therefore follows the theoretical conclusion.

Fig. 4

Normalized output power after scattering through varying thickness of phantom.

The simulation for the model illustrated in Fig. 2 was run. The photons were assumed to be launched as a point source at the origin. They were assumed to be horizontally linear polarized. Random values were taken for the scattering distance, scattering, and azimuth angles as described in the previous section. The refractive index for the outside medium (air) was assumed to be 1.0. The scattering medium was considered to be a homogenous suspension of $2.02\mu \mathrm{m}$ diameter polystyrene spheres in water. The density of scatterers is assumed to be 0.001 per cubic micron. The complex index of refraction of the spheres was 1.558 $-i$ 0.0006 and the refractive index of the medium (water) was assumed to be 1.33. The sample had a scattering coefficient $\left({\mu }_s\right)$ of $93.2{\mathrm{cm}}^{-1}$ , a reduced scattering coefficient $\left({{\mu }_s}^{\prime }\right)$ of $5.68{\mathrm{cm}}^{-1}$ , an anisotropy factor $\left(g\right)$ of 0.89, and an absorption coefficient $\left({\mu }_a\right)$ of $0.43{\mathrm{cm}}^{-1}$ at a wavelength of 990 nm. At 632 nm, for the same concentration of scatterers, the complex index of refraction of the spheres was 1.5892– $i$ 0.00113 and the sample’s scattering coefficient, reduced scattering coefficient, anisotropy, and absorption coefficient were $100{\mathrm{cm}}^{-1},8.9{\mathrm{cm}}^{-1}$ , 0.911, and $1.33{\mathrm{cm}}^{-1}$ , respectively. These values have been chosen so as to set to values obtained from real muscle tissue.¹⁸

The lens diameter, focal length, and height of the imaging lens were all taken to be 0.1 cm. The depth of the object was incremented in steps of 0.5 cm, and copolarized and cross-polarized images were recorded at wavelengths of 632 nm and 990 nm. The simulation was run for a launch of $1\times {10}^7$ photons.

Figure 5 shows the images from the simulation with a linearly polarized source and depolarizing diffuse reflector object. In each subfigure, each row represents the images obtained at a depth of 1, 1.5, 2, 2.5, 3, 3.5, 4, and 4.5 cm. In each row, from the left to the right are (a) orthogonal state at 632 nm, (b) copolarized state at 632 nm, (c) orthogonal state at 990 nm, and (d) copolarized state at 990 nm.

Fig. 5

Images from the simulation of diffuse reflector object illuminated with a linear polarized source: (a) orthogonal state at 632 nm, (b) copolarized state at 632 nm, (c) orthogonal state at 990 nm, and (d) copolarized state at 990 nm.

Figure 6 shows the images after being processed: (a) image after subtraction^{19, 20} of orthogonal state and fraction of copolarized state at 632 nm, (b) image after subtraction of orthogonal state and fraction of copolarized state at 990 nm, and (c) the image obtained after subtraction of (a) and (b), i.e., after the proposed method. Images (d) are the images obtained if only subtraction of the cross-polarized images at 990 and 632 nm (wavelength subtraction alone) is performed.

Fig. 6

Images after processing: (a) polarization filtering at 632 nm, (b) polarization filtering at 990 nm, (c) after the proposed method, and (d) image after wavelength filtering alone.

The improvement in contrast of the images in Fig. 6(c) over Fig. 6(d) is clearly visible, proving the merit of this model. The optimum fraction chosen for obtaining the images in Figs. 6(a), 6(b) was found to be 0.5. As can be seen, as the depth of the object increases, there is an increase in turbidity. The images formed after subtraction of the copolarized image from the opposite polarized state show better visibility of the object compared to the cross-polarized images. Once again if the wavelengths are compared, the images obtained at 990 nm show better contrast compared to the 632 nm images. This is because at 990 nm more photons have reached this depth without scattering than at 632 nm. Hence, choosing longer wavelengths for imaging gives better contrast images than shorter wavelengths. Still, the cross-polarized images at longer wavelengths have scattered components because the light backscattered from the medium reduces the contrast. Hence, subtraction of a fraction of the copolarized image from the orthogonal state removes this contribution.

Comparing the images Figs. 6(a), 6(b), as depth of the object increases, the visibility of the object decreases. The 632 nm images in Fig. 6(a) are due to photons backscattered from layers above the object layer and carry less of the object information. Hence, if this segment of image information due to photons that have not reached the object is removed, then the object information can be extracted. This is done in the images in Fig. 6(c) where a fraction of the 632 nm images in Fig. 6(a) is subtracted from the 990 nm images in Fig. 6(b). This final image subtraction process enhances the relative intensity of the object located at a deeper depth in the tissue.

This is best demonstrated by the intensity profile plots in Fig. 7 along the line containing the object for the cross-polarized images at 632 and 990 nm, the 990–632 nm image, and the image after the proposed method. If the wavelengths are compared, the images obtained at 990 nm show better contrast compared to the 632 nm images. Next, in comparison to the cross-polarized images, the images formed after wavelength subtraction (cross-polarized image at $990\mathrm{nm}-$ cross polarized image at 632 nm) show better visibility of the object. But still these images have scattered components from the medium. Hence, subtraction of a fraction of the copolarized image from the orthogonal state removes this, and the intensity profile of the image obtained through the proposed method shows only the object.

Fig. 7

Intensity profile along the line containing the object at 2.5 cm for the orthogonal polarized images at 632 and 990 nm; image obtained by wavelength subtraction and by the proposed method.

Comparing the images obtained at different depths, as depth of the object increases, the visibility of the object decreases. Figure 8 shows the digitized intensity profiles of the images obtained at different depths by the proposed method. The intensity profiles demonstrate the ability of the proposed method to image deep subsurface images in tissues. At lower depths of the target the profile is narrow and as the depth of the target is increased the intensity profile widens. This happens because at lower depths, a majority of the photons launched reach the target and the image is clearer. At higher values of the depths, the number of photons reaching the target is considerably reduced and hence the object is invisible in the light backscattered from tissue at layers above the target.

Fig. 8

Intensity profile along the line containing the object for the image obtained by the proposed method at depths of 1, 2.5, and 3.5 cm.

To demonstrate the effectiveness of the proposed method, during each simulation run, the number of photons actually hitting the target was counted. Figure 9 shows the plot between the percentages of photons reaching the object at different depths at the two wavelengths at which the images were recorded. It can be seen that for a launch of 100,000 photons the number hitting the object at each depth is higher at 990 nm compared to 632 nm. The object is more visible to photons at longer wavelength than photons at shorter wavelength. After a depth of 2.5 cm the number reaching the target was much less and not sufficient to provide the image of the object. Hence, the proposed method works effectively up to 2.5 cm in the tissue. And for depths greater than this, the object is no longer visible even by this method.

Fig. 9

Percentage of photons reaching a particular depth at 632 and 990 nm.

To validate the model of the object assumed, the DOLP was calculated from the Stokes vector of the photon after backscattering from the target. The DOLP was found to be approximately 0.8, a value less than unity meaning that the backscattered light is depolarized, justifying the model. It was also observed that the DOLP was the same throughout the range of optical depth since the target has been assumed to be a solid sample.

To further validate the combined polarization and wavelength subtraction technique, image contrast as a function of depth for the images obtained was computed. The contrast is defined as $(A_I-A_B)∕(A_I+A_B)$ where $A_I$ is the average intensity from light returning from the object and $A_B$ is that from light returning from the background. The computation for $A_I$ was done in the central nine $(3\times 3)$ pixels and for $A_B$ from the central $21\times 21\mathrm{pixels}$ . The computation was done on images at different depths and Fig. 10 illustrates how the contrast decays as target depth is increased.

Fig. 10

Graph showing how contrast decreases as depth increases for images obtained by the proposed method.

The image contrasts are listed in Table 1 for three different depths. On comparison of columns 2 and 3, the orthogonal images at 990 nm show better contrast values over their counterparts at 632 nm. Further examination reveals that between columns 4 and 5, polarization subtraction at 990 nm is better than at 632 nm. Values in columns 4 and 5 are larger than in columns 2 and 3, respectively, indicating improvement due to polarization subtraction. Column 6 indicates that if only wavelength subtraction was performed, the contrast is improved compared to the other images in columns 2, 3, and 5. The last column clearly demonstrates the superiority of the proposed method where the contrast is the highest among all other methods.

Table 1

Contrast for the images of the object obtained at depths of 1.5, 2.5, and 3.5 cm in the tissue.

4.

Conclusion

A method to enhance the contrast of tissue images at deeper depths has been proposed. Computer simulations based on a Monte Carlo model have been demonstrated for polarized photons launched from a point source at two different wavelengths. The resulting images indicate the effectiveness of combining the polarization subtraction along with wavelength subtraction. The employment of the different wavelengths as well as subtraction of a fraction of the copolarized image component allows for imaging at deeper depths of the tissue. This therefore promises to be a low-cost and effective method to diagnose abnormalities below the surface of the skin.