1. Introduction

There is a higher incidence of skin cancer than the combined occurrence of breast, prostate, lung and colon cancers [1]. Melanoma, which accounts for an estimated 4% of skin cancer cases, is responsible for approximately 75% of all deaths from skin cancer. The total deaths in the United States due to melanomas and other types of skin cancer are estimated to be more than 12,000 for 2014 [2]. Currently, the detection of melanoma relies on a subjective ABCDE (Asymmetry, Border, Color, Diameter and Evolution) test performed visually by dermatologists, general practitioners (GP) or primary care physicians (PCP) [3]. However, the ABCDE test provides a qualitative guideline and it requires a trained specialist to actually distinguish malignant lesions from benign nevi. Moreover, the ABCDE approach has a relatively high false-alarm probability (0.35–0.44, i.e., a specificity in the range 56% to 65%) and moderate detection probability (0.47–0.89) [3–5]. Since a false negative (i.e., a patient with malignant condition that is declared to have benign condition) could lead to metastasis (spreading to other parts of the body) and death, excisional biopsies are routinely performed even on lesions that are non-cancerous. It was estimated that the number of biopsies undertaken in nine geographical areas of the US between 1986 and 2001 was close to 60 for every melanoma detected [6]. One of the critical barriers in early skin-cancer detection is the lack of reliable non-invasive techniques [7] that can detect the cancer at an early stage with high detection probability (i.e., the probability of correctly detecting a malignant lesion) and low false-alarm probability (i.e., the probability of declaring a benign lesion as malignant).

Non-invasive techniques for skin-cancer detection include multispectral (MS) imaging [8–10], digital dermatoscopy and videodermatoscopy (sequential digital dermatoscopy) [11, 12], reflectance-mode confocal microscopy [13], ultrasound [14, 15], laser Doppler perfusion imaging [16], and optical coherence tomography (OCT) [17, 18]. Each technology presents some restrictions and limitations to non-invasively detect skin cancer with high detection probability and low false-alarm probability. For example, on one hand MelaFind, which is a device approved by the Food and Drug Administration (FDA), presents a high-level of detection probability (> 95%) [19], but high false-alarm probability (> 90%) [20]. On the other hand, Vivosight Multi – Beam System, another FDA-approved device, achieves a detection probability between 79–94% and a false-alarm probability between 4–15% for non-melanoma skin cancer lesions [21], but the challenge is that the suspicious lesion must be probed several times before such an accuracy is achieved, which makes the acquisition time prohibitively high.

Here we report a method for statistical inference, which uses the technique of dynamic thermal imaging (DTI) and it demonstrates a rapid, accurate and non-invasive imaging system for detection of skin cancer. DTI has already demonstrated to have high potential for the skin cancer diagnosis [22, 23] and it is a technique in which a thermal stimulus is applied to the suspected lesion and the thermal recovery is captured as function of time using an infrared (IR) camera [24]. Even though several groups have reported that the thermal recovery curves (TRCs) of a skin-cancer lesion and the surrounding healthy skin is different [25–30], these methods only partially extracted the information present in the temporal evolution of the recovery process. More specifically, the existing DTI techniques have neglected the temporal statistical features inherent in the thermal recovery process. To fully extract the vital information present in the recovery process, which will enable us to make a reliable inference on the malignancy of lesion, two problems must be solved. First, the recovery process must be viewed as a random function of time and its temporal statistical properties, such as its temporal correlations, must be mathematically characterized. Second, such complete statistical understanding of the thermal recovery process must, in turn, be utilized in a statistical-inference framework that yields the optimal decision rule for classifying a lesion as malignant or benign. Both of these problems are formulated and solved in the present work. The mathematical method reported here optimally extracts all the temporal information present in the DTI time series, and subjects the extracted vital information to optimal statistical decision theory. A pilot study is also undertaken on 140 human subjects at the University of New Mexico (UNM) Dermatology clinic to demonstrate the effectiveness of the method. We have demonstrated > 99% sensitivity and > 90% specificity while showing excellent robustness to statistical variation in algorithm-training and patient data collection. To the best of our knowledge, this is by far the highest reported accuracy and robustness for any non-invasive method for detection of skin cancer.

This paper is organized as follows: In Section 2 we utilize the knowledge of the TRCs to postulate the classification of a lesion as benign or malignant as a continuous-time detection problem. We explain the parametric stochastic model for the TRCs and how we obtained an analytical form of the auto-covariance functions (ACF), which we use to solve the continuous-time detection problem. In Section 3 we study the inclusion of a reference signal to the detection problem, obtained locally from the very same patient under study. More precisely, for each hypothesis, we define a self-reference signal from the tissue that surrounds the suspicious lesions. In Section 4 we present the results over a sample of 140 patients and a robustness analysis of the methodology. In Section 5 we present the discussion and summarize our main conclusions.

2. Postulation and solution of the detection problem

2.1. Use of dynamic thermal imaging for skin cancer detection

Skin cancers, like all solid malignant tumours, require a blood supply in order to grow larger than a few millimeters in diameter [31, 32]. Tumours induce the growth of new capillary blood vessels by producing specific angiogenesis-promoting growth factors. New blood-vessel growth continues through the progression from precancerous skin lesions to full-blown skin cancer as depicted in Fig. 1. The presence of new blood vessels and the increased blood supply somehow change the thermal response of the tumor cells when a stimulus is applied.

 

Fig. 1 Tumor angiogenesis in cancer at different stages: (a) The tumor release growth factors that activate the growing cells generating blood vessel sprouts. (b) The blood vessels feed the tumor that growths thanks to cell proliferation. (c) The tumor becomes vascularized and it starts to metastasize through the blood stream (from webpage [33]).

Download Full Size | PDF

Under this scenario, we assume that the patient condition is hidden within TRCs of suspicious lesions. Moreover, we assume that the malignancy of a lesion can be inferred only by monitoring the tissue of the mole. Later in Section 3, the proposed method is further generalized to include TRCs of the tissue that surrounds the suspicious lesion as a local reference. Such a local reference permits the compensation of any possible anomalous behaviour in the lesion thermal recovery, which, in turn, improves both the theoretical and empirical performance of the method.

In the next section we propose a physics-based stochastic model for the skin thermal recovery, which is later utilized in the context of continuous-time detection theory in order to define an optimal decision rule for classifying a lesion as malignant or benign. To the best of our knowledge, the TRCs have never been modeled as continuous-time random processes, and, as such, there is no known relationship between the tumor morphology or blood vessel development and the statistical properties of the TRCs. As a consequence, the model we propose in the next section can pave the way to future research in such a relationship in order to further improve the performance of the proposed method.

2.2. Physics-based stochastic model for the skin thermal recovery

Let us assume that the lesion boundary is defined by visual inspection over the mole region in a visible image. In Section 4.2 we explain how this is achieved. Since the area that defines the mole may contain malignant and benign tissue we need to look at the aggregated effect of the lesion TRC. The most natural and simple way to do such an aggregation is by computing the average TRC over all the pixels within the region that defines the mole.

In the transfer of heat in biological environments, Pennes bioheat equation [34] is generally considered as the most suitable method in all the bioheat models so far [35]. The Pennes bioheat equation is normally simplified as $k{\nabla }^{2}s+q_v=\rho C\frac{\partial s}{\partial t}$, where k, C and ρ are the coefficient of heat conduction of the skin, the thermal capacity of the skin and the tissue density, respectively. Here, s is the distribution function of internal temperature (refered in this work as the TRC) and q_v represents the internal heat source of the body. Even though this model is the most appropriated for the purpose of this work, we want the mathematical model of the TRCs to be simple enough to ensure a feasible solution with the available information. As such, we assume that the physics of the problem is simply governed by a heat equation, where an effective difussion constant, $\mathcal{D}$, captures the cumulative effect of all the subcutaneous thermal processes originated by k, C, ρ and q_v. We further assume that the model only will be affected by the variations on the depth of the lesion, x, due to the agregation of TRCs previously discussed. Therefore, the temperature of the skin sample is assumed to be modeled by the one-dimensional heat equation,

From the available subject data we have used in this work, the realizations of the random variables θ_j,₁, θ_j,₂ and θ_j,₄, with j = 0, 1, follow a Gaussian distribution, which is an assumption we use later. The realizations of the random variables θ_j,₃ and θ_j,₅ also follow a Gaussian distribution, but their variances are at least two orders of magnitud lower than their mean value. In this work we have assumed that θ_j,₃ and θ_j,₅ area random variables, keeping in mind that their mean values are positive (0.01 and 0.12, respectively) and that their variances are very small (2.2 × 10⁻⁵, and 1.5 × 10⁻³, respectively). Comparing θ_j,₃ and θ_j,₅ variances with respect to the variances of the other three parameters, we notice that θ_j,₁, θ_j,₂ and θ_j,₄ always dominate the behaviour of Eq. (3). An unexplored reasonable option is to also assume θ_j,₃ and θ_j,₅ as deterministic, which is one of the future lines of research we are exploring now.

The binary hypothesis-testing problem (also termed detection problem) of determining whether a measured TRC, say Y (t), is either from benign tissue (null hypothesis, H₀) or malignant tissue (alternative hypothesis, H₁) can be now formulated by the continuous-time binary detection problem

The distributions and the correlations between the random vectors Θ_j = [θ_j,₁…θ_j,₅], j = 0, 1 must be determined, in principle, from patient data with known condition. It is possible, in principle, to utilize the vector of random parameters Θ_j to perform a statistical decision regarding the malignancy of the suspected lesions, but it will not be necessarily optimal. It can be optimal if exact extraction from the TRCs of these parameters is possible, but this is not the case here due to the noise inherent in infrared imaging. Thus, we solve the continuous-time detection problem in Eq. (4) by first constructing an analytical auto-covariance function of the TRC under each hypothesis (and hence account for the majority of the statistical correlation in the TRC), and then following Grenander’s approach [37] to perform the statistical decision.

2.3. Analytical form of the auto-covariance function

The auto-covariance function (ACF) of the stochastic processes S(t; Θ_j) is defined by C_j (t, u) = E [S(t; Θ_j)S(u; Θj)]−E [S(t; Θ_j)] E [S(u; Θ_j), where E [·] denotes the (ensemble) expectation operation. By expanding the first expression, the jth (j = 0, 1) ACF can be recast as

The first term in Eq. (5) is E[θ_j,₁²], which can be easily estimated from the patients with known diagnosis. Following Bohrnstedt and Goldberger ideas [38] we have found that the terms with three random variables can be expressed by

As mentioned earlier, it was observed from the available subject data that these random parameters follow a joint Gaussian distribution, i.e. ${\theta }_{j,n}~\mathcal{N}\left({\mu }_{j,n},{\sigma }_{j,n}^2\right)$, then E[θ_j,_n] = μ_j,n, for j = 0, 1 and n = 1, 2, 3, 4, 5. Now, for a fixed t ∊ ℝ, E [exp (−θ_j,ℓt)] = exp $\left(-t{\mu }_{j,\ell }+{\sigma }_{j,\ell }^2{t}^2/2\right)$ also for j = 0, 1 but for ℓ = 3, 5. Moreover, the random variables given by the products of the form θ_j,_n exp (−θ_j,ℓt) can be approximated to as a random variable that follows a Normal–log-Normal distribution [39–41]. In particular, from Chen’s work we know that the covariance of these products can be approximated by [39]

Assembling these results into Eq. (5) we observe that one can compute the ACF analytically by simply estimating pairwise moments of the random parameters of the TRCs. These moments are estimated from a collection of TRCs by assuming that each measured TRC is a realization of the random process of interest. The resulting ACF from the parameters of patients with known benign and malignant conditions are depicted in Fig. 2(a) and Fig. 2(b), respectively. More information on the patient data is detailed later in Section 4. In both cases, the acquisition time for each TRC was 100 seconds. By looking the ACFs from the left plane, as shown in Fig. 2(c) and Fig. 2(d), one can clearly note some of the differences between them. Thus, each ACF encapsules the different statistical correlations that describe the TRCs under each hypothesis and, therefore both can be used to statistically describe the detection problem under Grenander’s approach.

 

Fig. 2 (a) Auto-covariance function for the null-hypothesis (H₀) estimated from patient data with known benign condition. (b) Auto-covariance function for the alternative-hypothesis (H₁) estimated from patient data with known malignant condition. In order to highlight their differences, (c) and (d) show the projection onto one of the left plane of the Auto-covariance function for the null-hypothesis and the alternative-hypothesis, respectively.

Download Full Size | PDF

2.4. Mercer’s theorem and the Karhunen-Loève expansion of the thermal recovery curves

According to Mercer’s theorem [42], the jth ACF (j = 0, 1) can then be expanded by the absolutely-convergent series

The two sets of eigenfunctions, namely ${\left\{{\varphi }_{0,k}\right\}}_{k=1}^{\infty }$ and ${\left\{{\varphi }_{1,k}\right\}}_{k=1}^{\infty }$, are two complete sets because the corresponding ACFs, C₀ and C₁, are symmetric and positive definite [42, 44]. The completeness of these two sets allow us to represent (in the mean-square sense) any process with either of these sets. We choose to represent each signal with its own set; namely, for the jth hypothesis

The main challenge in solving Eq. (11) by optimal-inference theory is the need to know the probability density function of the KL coefficients under each hypothesis and for each k = 1, 2, … (in order to be able to describe the probability density function of the likelihood-ratio function, to be introduced later); this is extremely difficult in general, if not impossible [42]. An important exception to this difficulty is when the random processes S(t; Θ_j) are Gaussian, which is an assumption we adopt here given the dominance of the Gaussian parameters θ_j,_n with j = 0, 1 and n = 1, 2, 4. Under this assumption and since the KL coefficients are always uncorrelated, they are also independent random variables, with $S_{j,k}~\mathcal{N}\left(0,{\lambda }_{j,k}\right)$; thus, the original continuous-time detection problem becomes the discrete (but infinite) detection problem between two Gaussian distributions with different (diagonal) covariance matrices. Hence, for the observation process Y(t), we have the discrete detection problem

The discrete-time detection problem cast in Eq. (12) has a likelihood-ratio function defined by

The test-statistic associated with (13) is obtained as usual by separating the terms that depend on the KL coefficients of the observation process and letting the remaining terms be absorbed by the threshold. After some algebra, we find that the test-statistic, Z is defined by

It is worth mentioning that there is a practical limit in the number of eigenvalue-eigenfunction pairs one can reliably extract from the estimated ACFs. Nevertheless, the KL expansion offers the optimality-under-truncation property, i.e., the mean-square error resulting from a finite representation of the process is minimized [44]. Such a property allows us to still optimally represent our processes and our resulting test-statistic, when the most important eigenvalue-eigenfunction pairs are used.

2.5. Neyman-Pearson decision rule

In this section, we describe how to optimally define the threshold, η such that detection probability is maximized for a fixed, prescribed false-alarm probability. By detection probability we mean the probability of announcing H₁ when H₁ is true; in symbols P_D = Pr(H₁|H₁) = Pr(L(Y) > η|H₁). By false-alarm probability we mean the probability of announcing H₁ when H₀ is true; in symbols P_F = Pr(H₁|H₀) = Pr(L(Y) > η|H₀). It is clear that the definitions of P_F and P_D require us to know the distribution of either the likelihood ratio or the test-statistic, which is what we address next.

The KL coefficients under H₀ are Y_0,k = Y_1,k = S_0,k, k = 1, 2, …. As such, the test-statistic under H₀, denoted here by Z₀, is given by

In summary, the test-statistic under each hypothesis is a linear combination of χ²-distributed random variables, with a different set of positive coefficients for each hypothesis. For a finite number of KL coefficients (16) and (17) are quadratic forms of the Gaussian random variables, S_0,k and S_1,k, k = 1, 2, …, K. The distribution of quadratic forms of independent and identically-distributed Gaussian random variables with positive coefficients was studied by Pachares [45]. Pachares’ main result was that the cumulative distribution function (CDF) of any finite linear combination of independent χ²-distributed random variables with one degree of freedom, X_k, with positive coefficients c = [c₁ c₂ ⋯ c_K ], i.e., the CDF of $Z_j=\frac{1}{2}(c_1{X}_1+c_2{X}_2+\cdots +c_K{X}_K)$, is given by

 

Fig. 3 False-alarm and detection probabilities parameterized by the threshold value, η, for different number of eigenfunctions used in the construction of the test-statistic (14)

Download Full Size | PDF

Once the desired false-alarm probability is specified, say P_F = α, the parameterized optimum threshold, η₀ can be obtained from Fig. 3 and Eq. (20). Such an optimum threshold is later used to classify patient data by comparing the test-statistic of a patient with unknown diagnosis: if the test-statistic exceeds the optimum threshold then is classified as malignant. More details are given in the following section. With the parameterized false-alarm and detection probabilities, one can construct the so-called receiver-operating characteristic (ROC) curve, a standard measure of the decision-rule performance that depicts the direct relationship between the theoretical false-alarm probability and the corresponding (theoretical) detection probability. The ROC curve of a perfect classifier correspond to a line that achieves 100% detection (sensitivity) for any value of the false-alarm probability (i.e., P_D = 1.0 for P_F ∈ (0, 1)). We show the ROC curve corresponding to the parameterized probabilities in Fig. 4, where it can be noted that, as expected, as we include more eigenvalue-eigenfunction pairs, more features are extracted from the TRC and the theoretical performance is improved, plateauing at a level, around the twelfth pair.

 

Fig. 4 The theoretical receiver-operating characteristic (ROC) curve graphically shows the expected performance of the detector as we increase the number of eigenvalue-eigenfunction pairs. The larger the number of the pairs utilized to construct the test-statistic, the more statistical features utilized and the better the performance of the algorithm

Download Full Size | PDF

2.6. Eigenvalue-eigenfunction pairs selection

Here we investigate the maximum number of reliable eigenvalue-eigenfunction pairs that can be extracted from estimated ACFs. In our study, we sorted the largest sixteen eigenvalues and examined the corresponding eigenfunctions for both hypotheses. It was observed that both eigenfunction sets are essentially the same under each hypothesis. However, as we require more eigenfunctions from the ACFs a slight phase difference begins to appear. Once the phase difference becomes notorious (e.g., for K > 10) the eigenfunctions lose their similarity, and as a consequence, the cross-orthogonality (among eigenfuctions corresponding to different hypothesis) is also lost. Our results indicated that after the twelfth eigenvalue, the eigenfunction start to become noisy. The source of this problem is estimation error in associated with the ACFs. It is expected that after having more patient data available with sufficient variability, higher K value can be considered.

To produce a systematic method for finding the highest reliable value for K, we first computed the Discrete Fourier transform (DFT) of each eigenfunction and for each hypothesis. It is observed that as the eigenvalue number is increased, the peak spectrum of each eigenfunction is slightly shifted and reduced in amplitude. The maximum K value can be determined by looking for abrupt changes in the peak of the spectrum of the set of eigenfunctions. In our example such abrupt drop occurs at K = 12.

2.7. Summary of the procedure for detecting malignancy of a suspicious lesion

To this end, we assume that the auto-covariance functions were correctly estimated from patient data with known condition (i.e., training patient data), and that the eigenvalues-eigenfunction pairs were also obtained and sorted based on the value of the corresponding eigenvalues. We denote the number of the stored (and sorted) pairs under each hypothesis by K. We also assume that the optimum threshold, η₀ was already defined by means of the NP decision rule as described in the previous section.

Let the aggregated TRC of the patient with unknown condition be denoted by Y (t). We first utilize the eigenfunctions to compute the KL coefficients under each hypothesis, i.e., $Y_{j,k=}{\displaystyle {\int }_0^T{\varphi }_{j,k}(t)}Y(t)$, for j = 0, 1 and k = 1, 2, …, K. This computation is depicted graphically in Fig. 5 (left dashed block). We next utilize the KL coefficients and the corresponding eigenvalues to construct the test-statistic as defined in Eq. (14). Such a calculation is depicted in Fig. 5 (middle dashed block). The test-statistic is later compared with the optimum threshold, η₀. This last step is termed the decision stage and it is depicted in Fig. 5 (right dashed block). By the extraction of the KL coefficients and the thresholding of the test-statistic by the optimally defined threshold, the detection problem is guaranteed to ensure the maximum achievable detection probability, and, as such, it presents an upper bound in the expected performance of the detector.

 

Fig. 5 Block diagram of the detection stage of the proposed algorithm. The KL coefficients are computed by using the eigenfunctions of each hypothesis. These coefficients and the eigenvalues are used to compute the patient’s test-statistic, which is later compared with the optimum threshold to declare the malignancy

Download Full Size | PDF

Next, we generalize the proposed approach to include a self-reference TRC of the very same patient under analysis. The idea is that the additional information introduced by this new TRC can compensate out anomalous behaviour of the lesion TRCs either improving the overall performance of the algorithm or achieving the same performance as with one TRC but requiring less eigenvalue-eigenfuncion pairs.

3. Generalization to a self-referenced approach

3.1. Theoretical rationale and generalized KL expansion for vector random processes

In this section, we study the inclusion of a reference signal to the detection problem, obtained locally from the very same patient under study. More precisely, for each hypothesis, we define a self-reference signal from the tissue that surrounds the suspicious lesions. The hope is that the by self-referencing the patient’s TRC, abnormal features from the lesion TRCs can be compensated out and, as a consequence, a reduced number of KL coefficients is required to correctly detect the malignant lesions.

For simplicity, we utilize one aggregated TRC of the lesion and one of the surrounding skin. As before, we define the lesion TRC as the average TRC over all the pixels within the lesion boundary. Similarly, we define the reference TRC as the average TRC over all the pixels outside the lesion boundary. The lesion boundary can be selected manually only by visual inspection by the practitioner that conducts the analysis or by a color-based image classifier.

Let us define then the new random process that describe the jth hypothesis as the vectorial random process

 

Fig. 6 ACFs for the case of vectorial random processes: (a) Autocorrelation function for the null-hypothesis (H₀) estimated from patient data with known benign condition. (b) Autocorrelation function for the alternative-hypothesis (H₁) estimated from patient data with known malignant condition.

Download Full Size | PDF

3.2. Karhunen-Loève expansion for vector processes

Once again, we wish to represent this matrix ACF as a linear combination of eigenvalues and eigenvectors as we did for the single-signal approach. Van Tress [44] and Oya et al. [46] stated that the optimal way to address this problem is by defining vector eigenfunctions, since it was proved by Kelly and Root [47] that the vector eigenfunction with scalar eigenvalue representation is a generalization of the KL expansion, and, as such, is optimal. Hence, let us assume that we have a complete set of vector eigenfunctions ${\left\{{\mathbf{\Phi }}_{j,k}\right\}}_{k=1}^{\infty }$, j = 0, 1 to represent the vector random processes of the jth hypothesis, where each vector eigenfunction is defined by ${\mathbf{\Phi }}_{j,k}(t)={\left[{\varphi }_{j,k}^{(1)}(t){\varphi }_{j,k}^{(2)}(t)\right]}^T$, and the corresponding vector random process of the jth hypothesis X(t; Θ_j) can be expanded over this vector eigenfunctions by

3.3. Solution of the dual-signal detection problem

The resulting vector eigenfunctions allow us to represent the dual TRC under each hypothesis with scalar expansion coefficients by

In the following sections, we evaluate the performance for the single- and the dual-TRC approach.

4. Results

In this section we demonstrate the efficacy and delimit the scope and robustness of the skin cancer detection algorithm by measuring its empirical performance. The metrics we utilize to evaluate the empirical performance of the algorithm are the empirical false-alarm and detection probabilities. The empirical false-alarm probability, P_F,e, is defined as the ratio between the true-negatives (patients with benign condition as dictated by the biopsy that are declared as benign by the algorithm) and the total number of benign patients. The empirical detection probability, P_D,e, is defined as the ratio between the true-positives (patients with malignant condition as dictated by the biopsy that are declared as malignant by the algorithm) and the total number of malignant patients.

A cohort study with 140 subjects was performed to investigate the proposed approaches. Fifty eight percent of the subjects were male and, from the biopsy result, out of the 140 subjects 82 had benign condition and 58 had malignant condition. Out of those 58 subjects with malignant condition, 6 were diagnosed with malignant-melanomas (MM), 42 with basal-cell carcinoma (BCC) and 10 with squamous-cell carcinoma (SCC). The subjects were diagnosed by means of excisional biopsies performed at the University of New Mexico (UNM) Dermatology Clinic, New Mexico, USA, and all patient data was acquired at the same center. The acquisition harware used in this setting is described next as well as the data pre-processing stages required to apply the proposed techniques.

4.1. Image acquisition hardware

We performed DTI with three components. The first component is a cooling unit that is used to impart the temperature stimulus to the lesion and the surrounding skin tissue. Two different cooling unit were used in our study. The first one was a a Ranque-Hilsch vortex tube that generates an oil-free, moisture-free, ultra-quiet air flow. It was later replaced by a commercially available air-conditioning (AC) unit due to its portability. It was observed that by properly modifying the time the cool air was applied to the skin, the imparted temperature was almost the same for both cooling units.

The second component is an infrared marker, which is used for correction of involuntary movement of the subject (i.e., image registration); the IR evolved from a canvas paper marker to a square piece of plastic with a square opening in the middle. Since the only purpose of the marker is to aid in the registration of the infrared sequence of frames, changing the material of the marker did not change the acquisition protocol.

The third component includes the imagers. The first imager is a commercial visible still camera that is used to capture a reference image before the DTI acquisition commences. The second and most important imager is a longwave infrared (LWIR) camera that is used to capture a sequence of frames of the thermal recovery of the skin after the cool stimulus is applied. The LWIR camera consists of a 320×256 focal-plane array (FPA) of quantum-well infrared photodetectors (QWIP) operating at 60K. The noise equivalent temperature difference (NEDT) of the FPA is 20mK and the QWIP camera is fitted with a 50mm, f/2 LWIR lens, yielding an approximate spatial resolution of 300 microns per pixel. The QWIP camera was chosen for our study because it has higher array uniformity, lower NETD and high spatial resolution as compared with other IR camera technologies [22, 48].

All the components of the acquisition hardware are pictorially represented in Fig. 7(a) whereas Fig. 7(b) shows the infrared imager as well as the acquisition software. Next we explain the imaging procedure.

 

Fig. 7 Acquisition hardware utilized to acquire the patient datasets. (a) Prototype and (b) Infrared imager and aquisition software

Download Full Size | PDF

4.2. Imaging procedure

After informed consent, each subject was escorted to a designated room in the UNM Dermatology Clinic to perform the imaging procedure. The temperature of the room was controlled to be between 20°C to 22°C to make sure that all the patients were exposed to the same temperature before applying the cooling stimulus to the area of interest. At the beginning of the procedure, the square registration marker was placed around the lesion with the lesion centered in the opening, as shown in Fig. 8(a). A visible image of the lesion was then taken with the digital camera for reference. After collection of the visible image, a 15 second infrared image sequence of the marked area was collected to serve as a baseline. Later, the subjectâĂŹs skin within the marker opening was cooled for 15 or 110 seconds, depending of the cooling unit used. After cooling, the exposed area was allowed to warm up naturally to ambient temperature. During the warm-up phase, thermal images of the skin were captured for a total of 2 minutes at a rate of 60 frames per second with the QWIP camera. All the thermal images were recorded using an uncompressed 14-bit format. The total time required to complete the entire imaging procedure was less than five minutes.

 

Fig. 8 Example of a patient dataset: (a) example of one square plastic marker used in the data acquisition step; (b) first frame of the infrared sequence, note that the visible and this frame are spatially aligned; and (c) the thermal recovery curves (TRCs) for the labeled pixels in (b).

Download Full Size | PDF

If the subject was scheduled for a biopsy, the biopsy was performed following the data collection by the attending dermatologist and sent to pathology for diagnosis. The biopsy results were delivered to us within the two weeks following the imaging procedure. Some patients were clinically diagnosed with a benign condition by the staff, and, therefore, no biopsy was performed. These patients are considered as control patients and are included in the set of benign patients.

4.2.1. Image registration

Since involuntary movements of the patients cannot be avoided, image registration must be performed over the infrared sequence of images. Moreover, to correctly reference the lesion location within the IR sequence (i.e., the mole, which not necessarily can be spotted in the IR sequence), the visible picture must also be spatially aligned with the IR sequence. Therefore, the visible image is considered as an additional frame for the purposes of the registration process.

The registration is conducted as follows. First, we use the Harris corner-detector algorithm [49] to automatically detect the four corners of the plastic marker amongst the entire sequence of frames. Second, by assuming rigid movement of the scene, we estimate an affine transformation matrix that maps such a movement between the corners of consecutive frames (one matrix is estimated for each pair of consecutive frames) [50]. Third, we utilize the inverse of each transformation matrix to align each frame with respect to the first frame of the sequence [51]. After registration, both the visible image and the entire IR sequence are spatially aligned, generating a three-dimensional (3D) array real numbers that we term the patient dataset. Figure 8(b) depicts the first IR frame after the cooling was removed of the same example case presented in Fig. 8(a); note that both the visible and the first IR frame are spatially aligned. The thermal recovery curves (TRCs) of the labeled pixels are shown in Fig. 8(c), where it can be noted that there is some non-uniformity in the cooling process that make these TRCs to start at different initial temperature.

4.2.2. Camera calibration

In order to have a temperature measurement of the skin surface as accurate as possible, the QWIP camera must be radiometrically calibrated. As in any FPA, the camera suffers of the nonuniform response of its detectors (a problem known in the literature as non-uniformity) and it is compensated by means of non-uniformity correction (NUC) tables performed and stored during the factory calibration process.

The radiometric calibration is achieved by means of the two-point calibration technique [52]. This calibration is performed by placing, in the field-of-view (FOV) of the camera, a uniform-intensity calibration device such as a black-body source at two distinct and known temperatures [53]. The gain and the bias of each detector are then calibrated across the array so that all detectors produce a radiometrically accurate and uniform readout at the two reference temperatures. The reference temperatures where chosen to be within the normal temperature of the skin, i.e., 25°C and 40°C. Examples of thermal recovery curves after the temperature calibration was performed were already shown in Fig. 8(c).

4.3. Clinical application, analysis and discussion

Recall we have 140 subjects, out of which 58 had a malignant condition (including melanoma, basal-cell and squamous-cell carcinoma) and 82 subjects had benign conditions. To reiterate, the subjects were diagnosed by means of excisional biopsies performed at the UNM Dermatology Clinic, New Mexico, USA.

By performing the training over different sizes and permutations of datasets from the 140 patients while testing the method on the remaining patients that where not used in the training, we found that at least 110 patients with known conditions and ten KL coefficients may be used to train the method in order to perfectly separate both benign and malignant conditions for the patient datasets used in the training stage. In order to assess the robustness of using a size of 110 patients in the training dataset and 10 KL coefficients, we repeated the training 200 times, each time with a distinct (but randomly selected) permutation of 110 patients from the totality of 140 patients. Each selection of the 200 training dataset permutations yielded a lesion classifier. We then tested the method by studying the performance of each of the 200 classifiers, namely, we used each classifier to determine the condition of each of the 30 remaining patients (outside of its training dataset). For the purpose of this study, we selected the theoretical false-alarm probability to be P_F = 0.01. Since the eigenvalues were computed for each one of the trained classifiers, the corresponding optimal decision threshold varied for each classifier.

The empirical results demonstrate that the method achieves 100% accuracy (i.e., 100% sensitivity and 100% specificity) for 46% of the 200 training selections, which we term the highly reliable training data-sets, when we specify a theoretical false-alarm rate of 0.01. Moreover, the tested methodology achieves 100% sensitivity and 95% specificity (i.e., 5% false-alarm rate) for 76% of the 200 training datasets, and 100% sensitivity and 90% specificity for 93% of the 200 training datasets. These results demonstrate that by using highly reliable training datasets, the proposed technique is capable of correctly classifying both benign and malignant skin-cancer conditions with unprecedented accuracy. In addition, we have observed variability in the sensitivity (below 100%) within only 3% across all the 200 training datasets. Moreover, for different possible lesion boundaries (required to define the pixel-to-pixel averaged TRCs), we observed that the method presents a variability of 0.1% due to a change of 30% in the radius of the selected lesion, proving the robustness of the proposed method to variability in the selection of the lesion boundary by the operator or the medical practitioner performing the test.

After completing the above study, 11 new patient TRCs with known biopsy results were acquired: by means of biopsies 5 were diagnosed as malignant and 6 as benign. By testing the method on this new set of patients while using the same classifiers that were generated by the 200 training datasets from the previous study, we have correctly identified all malignant lesions and misclassified only one benign lesion for 76% of the original 200 classifiers. Moreover, we were able to correctly classify all new 11 lesions for 36% (71 classifiers) of the original 200 classifiers. Out of these 71 classifiers, 58 were from the “highly reliable dataset” identified in the earlier study (original set of 140 patients) as described in the previous paragraph, and 13 of the 71 were from the data set permutations that misclassified only one benign in the original set of patients. This observation shows consistency in the performance of the classifiers that were trained by the highly reliable datasets.

When the dual-TRC approach was utilized to classify the lesions we have observed that the mean theoretical performance over the same 200 permutations is improved with the help of the self-reference signal introduced as shown in Fig. 9. Here we utilize the area-under-the-ROC-curve (AUC) metric to compare the performance and it is defined as the computed area below the receiver-operating-characteristic (ROC) curve of the detector. (A perfect classifier, i.e., with 100% accuracy, will have an AUC of 1.) The most important characteristic of the proposed alternatives is that the detection is still performed over an scalar expansion coefficient, which has equivalent statistical features of the dual TRCs in a similar fashion as the KL coefficients characterized the single TRCs. The real difference in performance is observed in the empirical performance when 110 patients are utilized to train the algorithm. In this case, we observed that perfect classification of all malignant cases is achieved in all the permutations with only 5 eigenvalue-eigenfunction pairs per hypothesis when the highly reliable training datasets are utilized. (Again, this accuracy is achieved when the theoretical false-alarm rate is specified to be 0.01.) The same performance was achieved by the single-TRC approach with 110 training patients with exactly the double of the eigenvalue-eigenfunction pairs. It seems, therefore, that the inclusion of the reference signal can actually compensate some anomalous behaviour of the the lesion TRCs. Nevertheless, there some numerical issues when computing the eigenvalue-eigenfunction pairs that must be compensated in order to further validate this extension of the proposed methodology.

 

Fig. 9 Comparison of the mean theoretical ROC curves over 200 permutations when 110 training are used to train the single-TRC algorithm (blue) and the dual-TRC algorithm (red). Comparison is made by using the mean AUC for different number of used eigenvalue-eigenfunction pairs, using 110 patients to train the algorithm.

Download Full Size | PDF

In Table 1 we compare the results between the proposed algorithm (last row) and other methodologies previously discussed. The ABCDE test is only applied to melanomas, but it represents a classical approach to non-invasively detect skin cancer. Melafind achieves good performance in terms of sensitivity (high detection probability) but poor specificity (high false-alarm probability). Vivosight, on the other hand, achieves good sensitivity and specificity (i.e., high detection probability and low false-alarm probability). As it was aforementioned, the suspicious lesion must be probed several times before such an accuracy is achieved, which makes the acquisition time prohibitively high for clinical applications. The only approach found in the literature that utilized DTI with reported sensitivity and specificity is our previous work [30]. In that work, we performed classification by a distance-based classifier. To do so, we computed the aggregated TRCs inside and outside the suspicious lesions and then compared the normalized Euclidean norm of their distance. Malignant lesions were expected to have a bigger difference between these two curves, and, as such, have bigger values of their Euclidean distance. Even though this approach obtained good results, it did not extracted all the statistical information from the TRCs and neglected the order in the data, i.e., the temporal evolution of the TRCs; this explains why the method proposed here achieves better results than all other non-invasive techniques. Moreover, the rapid acquisition time and relatively fast processing are two other clear advantages of the proposed method as compared with other techniques for skin cancer detection.

Table 1. Comparison between the proposed methodology and other non-invasive techniques

View Table

5. Conclusions

In order to construct a suitable stochastic model for the TRCs we have simplified the bioheat equation and solved it for properly established boundary conditions. We used this model to derive analytical auto-covariance functions for the proposed hypotheses whose parameters are estimated from patients with known conditions. An optimal statistical decision rule was then derived for which the sensitivity is guaranteed to be at a maximum for every prescribed false-alarm probability. The algorithm was trained using TRCs from 110 human subjects and tested on 30 human subjects with unknown diagnosis. All the malignant subjects were correctly identified (100% sensitivity) for a false-alarm probability below 1% (specificity > 99%). Robustness analysis was undertaken with various permutations of testing and training datasets, and under different selections of the lesion boundaries. The maximum variability in the sensitivity of the method was 3% and 0.1%, respectively, for 200 different random permutations in selecting the training set and different radii of the region that determine the lesion boundary.

The detection method exploits, in a probabilistic sense, the temporal evolution of the warming-up process, as obtained from the thermal-image sequence, in order to extract statistical features that can be utilized to discriminate a case of a benign lesion from a malignant lesion. The term evolution is emphasized to highlight the importance of the temporal order of the warming-up data obtained from the image sequence. This is why the performance is far better than that obtained when the Euclidean norm of the thermal-image sequence is used because the latter, unlike the proposed method, does not depend on the ordering of the images in the sequence [30].

The success of the proposed methodology relies on knowledge of the auto-covariance function of the warming-up process under each hypothesis, which has not been known (to the best of our knowledge) prior to this work. We have presented the construction of such auto-covariance function by parameterizing the warming-up process using a few random variables and by estimating (from patient-data with known conditions) the probability distribution for these random variables and their relevant correlation functions.

We have also introduced and discussed an approach to represent vector random processes that in the literature is treated as the vector generalization of the KL expansion. Using this representation we have solved the detection problem in which each hypothesis is characterized by a dual TRC. The approach followed has the property of using scalar expansion coefficients, equivalent to the KL coefficients for the scalar TRC approach. Hence, after the expansion is applied to the vector random processes, the dual-TRC detection problem becomes equivalent to the single-TRC detection problem. As a consequence, once the expansion coefficients are characterized, they contain most of the statistical information of both of the signals of each hypothesis, and the solution is trivial because all the properties and limitations discussed for the single TRC approach hold for this new alternative. We explored both the theoretical and empirical performance of the proposed alternative using the same training settings utilized to test the robustness of the single-TRC approach. By comparing the theoretical performance for different number of eigenvalue-eigenfunction pairs, the dual-TRC approach achieves the same performance as the single-TRC approach by requiring less eigenvalue-eigenfunction pairs. The empirical performance follows the same trend, allowing, for example, 100% detection of malignant patients using one half of the eigenvalue-eigenfunction pairs that the single-TRC approach required for the same training/testing setting.

To the best of our knowledge, this work reports the highest accuracy and robustness for any non-invasive method for detection of skin cancer. As a consequence, the method will be valuable in clinical applications where the accessibility to trained dermatologists is difficult. Moreover, the methodology is applicably to any continuous-time hypothesis testing problem for which an analytical parametric form of the auto-covariance function is obtainable.