1.

Introduction

Mapping the degree of polarization (DP) of backscattered light with respect to its initial polarization can provide valuable diagnostic information about the superficial and subsurface structures of the skin. It can be used to extract tissue characteristics that are impossible to discern from conventional visual or photographic methods. One of the major strengths of polarization measurements is the ability to separate the light backscattered from superficial tissue layers and the multiple scattered lights coming from much deeper tissue layer.¹ This distinction is based on conservation of initial polarization state in the superficial tissue layer while in the deeper tissue the light becomes depolarized. For this reason, polarization analysis can be used to extract tissue characteristics that conventional visual or photographic methods are not capable of discerning.^{2, 3, 4, 5} Recent technological advances have resulted in the design of a handheld polarized light camera that may further expand medical applications of the polarization degree mapping.⁶ Several optical characteristics of the tissue determine changes in polarization state of the backscattered light relative to the incident light. First, the polarization state of the incident beam can be modified by the specular reflection from the surface heterogeneities. Second, the spatial variations of the light absorption coefficient influence the mean path of the detected photons inside the tissue. Higher absorption coefficients result in smaller average pathlengths of detected photons and, therefore, a higher polarization degree of light, backscattered from the regions with increased absorption⁷ (to make this effect visible, one should exclude multiply scattered photons from detection, because scattering quickly randomizes the polarization state of photons). And third, the tissue may have an anisotropic birefringence effect due to structural anisotropy. For example, rods or fibers could change the degree of polarization according to the mutual orientation of their axes and polarization vectors.⁸

Analysis of the 2-D distributions of DP can help to reveal such anisotropic structural features. For example, to evaluate the presence and extent of skin fibrosis developing after an injury (initiated either by burns, trauma, or ionizing radiation used to treat cancer), polarization images can be investigated. The skin dermis can be characterized by the density and preferential orientation of collagen fibers. As a result of an incurred injury, collagen fibers within the dermis may organize into bundles that are not observed as structural components of the normal skin.^{9, 10} The ability to determine the structural components of skin with incipient fibrosis may enable physicians to avoid keloid formation and radiation treatment overdoses, or intervene with an antifibrotic treatment at an earlier stage.^{11, 12, 13, 14} Unfortunately, due to the highly scattering nature of tissue, in many cases, the relationship between the observed polarization of detected light and internal skin structure is distorted by depolarization, resulting from multiple scattering of photons. Surface roughness as well as speckle may also blur the image. For example, in the case of optical coherent tomography Rogowska and Brezinski observed image distortion due to speckles.¹⁵ Additional distortions of noisy polarization degree images can be easily introduced by data processing. Therefore, it is important to develop special procedures to enhance the internal skin structures, such as fibrosis, from noisy digital images. This paper demonstrates the effectiveness of polarized digital photography for visualization of skin fibrosis structures. The algorithm, based on Pearson correlation analysis, was developed for the polarization imaging of the tissue surface, and characterization of the size and orientation of heterogeneities within the skin fibrosis zone.

2.

Materials and Methods

The protocol used in this study was institutionally approved in accordance with the guidelines of the Institute of Laboratory Animal Resources, National Research Council. For this study, athymic nude mice were immobilized and a selected area of skin on the right flank was irradiated with an x-ray irradiation dose of $30\mathrm{Gy}$ . The area of treatment was approximately $20\times 20{\mathrm{mm}}^2$ . Such treatment typically causes earlier fibrosis of the dermal layer of the skin,^{10, 16} which developed during 2 to 3 weeks after irradiation. The measurements of polarization degree were performed 20 days after x-ray exposure.

The mouse skin was illuminated by an expanded linearly polarized probe light from a diode laser of $650\mathrm{nm}$ wavelength. A digital camera (EOS D60, Canon Inc. Japan) with an objective of $0.31\mathrm{m}$ focal length was used to produce images of the mouse skin. The images were recorded into files that showed the intensity distribution of the detected red light. Each image was $3072\times 2048$ pixels and the output of each pixel consisted of a 16-bit digital signal. The distance between the neighboring pixels corresponded to a physical distance of approximately $0.0075\mathrm{mm}$ on the surface of the object. Calibration of the camera response for wide beams was performed using a set of neutral density filters. The camera response was linear in the range of light intensities used in the experiments. The probe beam was aimed at an angle of approximately $15\mathrm{deg}$ relative to the direction of normal to the mouse skin surface, while the camera objective axis was parallel to this normal. Such orientation prevented a specularly reflected probe beam from reaching the objective. A polarization analyzer was positioned in front of the camera objective. The transmission axis was either parallel or perpendicular to the polarization vector of the incident beam, providing copolarized or cross-polarized images, respectively. The light of the diode laser backscattered from mouse skin created speckles in the front of the camera matrix. The speckle was estimated to be 1 to 2 pixels, while the appropriate intensity variance was 5 to 10%. This enabled us to ignore the influence of the speckles on the recorded image.

3.

Fast Fourier Transform Filtering

A conventional visual image of the x-ray-irradiated zone 20 days after radiation is shown in Fig. 1(a) . No feature changes of the skin structure can be observed.

Fig. 1

(a) Gray-scale visual image of the nude mouse skin 20 days after x-ray irradiation and (b) the image of the DP in 8-bit gray-scale format that corresponds to the same area as in (a), where points $A$ and $B$ represent x-ray irradiated and nonirradiated areas, respectively.

The matrix of the degree of polarization of the backscattered light with respect to its initial polarization can be shown as¹

The image of matrix $\mathsf{G}$ , corresponding to the conventional visual image of Fig. 1(a), is presented in Fig. 1(b), where the irradiated zone in the central part of the field marked $A$ can be distinguished from the background zones of the periphery marked $B$ . Similar patterns were also observed in the x-ray-irradiated area of additional mice 20 days after irradiation. A typical image of nonirradiated mouse skin [see Fig. 2(a) ] revealed no actual difference throughout the pattern of the degree of polarization as shown in Fig. 2(b). The observed difference in Fig. 1(b) may be attributed to the presence of collagen bundles observed in histological patterns.⁹ These bundles are formed in the dermis as a result of x-ray irradiation. Despite visualizing some features, the polarization degree pattern is significantly distorted by high-spatial-frequency noise. The $\mathsf{G}$ matrix is formed from the raw data in the space domain. It has been transformed into a spatial frequency domain matrix ${\mathsf{G}}^{\prime }$ by discrete fast Fourier transform (FFT). This transform has been realized by the cfft function of Mathcad 12, applicable to matrices of an arbitrary size. Particular spatial frequencies were filtered out from the signal by replacing corresponding elements of matrix ${\mathsf{G}}^{\prime }$ with zeroes:

To illustrate the results of suggested data processing, i.e., elimination of the high spatial frequencies and zero-frequency components present in the raw signal, a modified image of the same region of interest, as in Fig. 1(b), is presented in Fig. 3 . Chosen values of $k=1$ and $K=100$ [see Eq. 2] correspond to spatial frequencies of $\sim 20$ and $\sim 0.2{\mathrm{mm}}^{-1}$ , respectively [estimated, as $1∕2\pi k\left(0.0075\right)$ ]. Thus, the data processing removes image irregularities less than 0.05 and greater than $5\mathrm{mm}$ . To make a reasonable choice of $k$ and $K$ indices in this case, a subjective perception of the image structure of Fig. 3 was accepted. Despite the loss of some information, this procedure made visible structural characteristics that were previously masked by the high-spatial-frequency noise in the original image [compare Figs. 1(b) and 3]. Intensity noise can originate either from speckles or skin surface roughness, which may give some specular reflection in the direction of the camera lens. Additional noise in the polarization degree image may be caused by a reorientation of polarized light due to light scattering or effects of randomly oriented birefringent tissue fibers. Besides, amplitude of noise, presented in matrix $\mathsf{G}$ , should be at least doubled, in accordance with Eq. 1.

Fig. 2

(a) Conventional visual image of untreated area of a mouse and (b) the corresponding image of the DP.

Fig. 3

Image of DP of the region presented in Fig. 2 after data processing involving Fourier transformations and filtering, as described in the text.

4.

Pearson Correlation Analysis

It is worth estimating the size and orientation of heterogeneities that stand out from the noise in Fig. 1(b), but still are noticeably distorted by it. To do this, we developed a fractional statistical analysis of the digital pattern, based on calculating the Pearson correlation coefficient. To evaluate spatial correlation in the raw data, we introduce two similar submatrices $\mathsf{A}$ and $\mathsf{B}$ of size $(2m+1)\times (2n+1)$ that are formed from the original matrix $\mathsf{G}$ . Matrix $\mathsf{A}$ is centered at a given point with coordinates (indices) $(i_0,j_0)$ , while matrix $\mathsf{B}$ is shifted relative to $\mathsf{A}$ by some vector $\mathbf{r}$ . The Pearson correlation coefficient between two matrices $\mathsf{A}$ and $\mathsf{B}$ is defined by the expression

Evidently, the easiest analysis may be performed by scanning matrix $\mathsf{B}$ throughout pattern over all integer values along $x$ and $y$ axes. To estimate the dependence of the Pearson correlation coefficient on the arbitrary scale and azimuthal angle, we superimpose an auxiliary submatrix $\mathsf{SM}$ rotated at this angle, which includes the region of interest, onto the original matrix $\mathsf{G}$ . Values of the nearest pixels of the original matrix $\mathsf{G}$ are attributed to the corresponding pixel values of submatrix $\mathsf{SM}$ . Figure 4 geometrically demonstrates this procedure. To create submatrix $\mathsf{SM}$ of size $(2m+l)\times (2n+1)$ , centered at a point with coordinates (indices) $(i_0,j_0)$ and rotated at an angle $\alpha$ relative to the direction of the rows of the matrix $\mathsf{G}$ , we used the following shift and rotation transformations:

Fig. 4

Illustration of the procedure to design the submatrix of size $(2n+1)\times (2m+1)$ , centered at an arbitrary point $(j_0,i_0)$ of original matrix and tilted at an angle $\alpha$ with respect to the rows of original matrix ( $x$ axis).

Figure 5 shows the behavior of the Pearson correlation coefficient in the vicinity of point A [shown in Fig. 1(b)], close to the center of the irradiated skin region, for two perpendicular directions that were tilted relative to the horizontal axis at angles of 10 and $100\mathrm{deg}$ , respectively. The presented data were calculated for submatrices of $29\times 31$ pixels that correspond to a size of $0.22\times 0.23\mathrm{mm}$ as the region of interest. Figure 5 shows that for $10\mathrm{deg}$ , the correlation coefficient drops to a level of about 0.6 and then remains almost constant up to a distance of about $0.5\mathrm{mm}$ ; for $100\mathrm{degs}$ there is an obvious periodicity with spatial period of approximately $0.26\mathrm{mm}$ . This periodicity may be explained by the existence of some elongated structure in the vicinity of this point. The azimuthal dependence of the Pearson correlation coefficient presented in Fig. 6 not only supports this conclusion, but also helps to determine the directionality of this structure. In fact, the azimuth variation of the Pearson coefficient may be realized by scanning over a circumference of a given radius. It may provide additional information regarding the structural characteristics in the vicinity of the point of interest.

Fig. 6

Pearson correlation coefficient in the vicinity of point $A$ of Fig. 1(b), as a function of azimuthal angles for two distances 0.13 and $0.26\mathrm{mm}$ .

Fig. 5

Pearson correlation coefficient in the vicinity of point $A$ of Fig. 1(b), as a function of distance for azimuthal angles 10 and $100\mathrm{deg}$ .

Our quantitative analysis shows that the spatial correlation coefficient behaves differently, as a function of scale, far from the x-ray-irradiation zone [undamaged skin, e.g., close to point $B$ in Fig. 1(b)]. As presented in Fig. 7 , the correlation coefficient decreases with increasing distance, reaching almost zero at distances greater than $0.05\mathrm{mm}$ . This behavior demonstrates a low (practically nonexistent) level of any regular structure at the scales shown for both azimuths. Angular dependence of the correlation function, centered at point $B$ , can substantiate this conclusion.

Fig. 7

Pearson correlation coefficient in the vicinity of point $B$ of Fig. 1(b), as a function of distance for azimuthal angles $-10$ and $80\mathrm{deg}$ .

The described imaging technique has a potential use for estimating characteristic size and directionality of superficial hidden structures, for example, in dermatology or stomatology. How deeply a hidden structure may be traced by this analysis strongly depends on the scattering coefficient of the tissue and on the birefringence of superficial tissue layers. Polarized photons lose their initial polarization state in about 5 to 10 scattering events. Thus, residual degree of polarization remains detectable for photon path lengths of a few millimeters in tissues.²² The scattering coefficient of the human derma is about $200{\mathrm{cm}}^{-1}$ in the visible range.²³ Correspondingly, tissue structures at depths of 250 to $500\mathrm{\mu }\mathrm{m}$ below the skin surface are expected to be accessible for polarization analysis. When fibrosis has developed, the derma becomes thinner. This fact opens a promising opportunity to discover the early stage of fibrosis. Note that the probe light with wavelengths longer than the one used here $\left(630\mathrm{nm}\right)$ can penetrate deeper into the tissue, since scattering coefficient decreases with wavelength.

5.

Conclusions

Measurements of the DP enable digital polarized photography to reveal hidden structures beneath the skin surface, when the skin is characterized by special conditions such as fibrosis from x-ray irradiation. The Fourier transformation of DP data with subsequent elimination of high spatial frequencies, followed by inverse Fourier transformation, substantially enhances the informational content of the obtained images. Such a methodology makes it possible to estimate the characteristic sizes and directionality of possible regular structures. In this way, Pearson correlation analysis enables us to better investigate internal skin structures.