I Introduction
The lentigo is a hyperplasia that affects the skin. It comes from the proliferation of melanocyte cells at the dermoepidermic junction. This leads to the disorganization of the regular cellular network [1]. Clinically, this disorder is assessed visually on the skin surface or through biopsy. Reflectance confocal microscopy (RCM) imaging is increasingly used to explore various skin lesions [2, 3], including lentigo. Figure 1 shows examples of images from patients with and without lentigo. Various studies have attested of the usefulness of RCM for cancer and other tumor diagnosis [4]. In [1], the authors reported good correlation between RCM and histology in the case of melanoma. Studies of RCM for treatment followup [5, 6, 7] and guidance [8] have also been published.
However, RCM images are up to now mainly analyzed visually. Image processing methods could be helpful to exploit their potential and provide aid for medical decision making. Few of such methods were reported in the literature. In [9]
, Luck et al. developed a nuclei segmentation method based on a Gaussian model for the nuclei reflectivity and a truncated Gaussian distribution for the intensity of the cytoplasm fibers. Their bayesian classification algorithm relies on a Gaussian Markov random field to ensure spatial correlation. Another application of RCM was developed and validated by Kurugol et al. to identify the dermoepidermal junction. This method is based on statistical classification of texture features
[10, 11]. Hames et al. [12, 13]proposed a skin layer segmentation method based on a logistic regression classifier. An SVM classification method was developed in
[14] based on speeded up robust features. This method was applied to identifying skin morphological patterns using RCM image texture. Finally, a waveletbased classification method was developed in [15]to distinguish benign and malignant melanocytic skin tumors. This method, with which we compare our’s, consists of classifying a vector of 39 features using a decision tree approach. In this paper we proposed a bayesian method to jointly reconstruct RCM true reflectivity images while classifying images as lentigo or healthy.
The first contribution of this paper is a hierarchical Bayesian model that allows a set of RCM images to be classified into healthy and lentigo classes. Each image is assumed to be corrupted by a multiplicative speckle noise with a gamma distribution. A truncated Gaussian distribution is then assigned to each image to classify, constraining these images to be positive. Prior distributions are finally assigned to the means and variances of these truncated Gaussian distributions, to the noise variances, and to the image labels. The joint posterior distribution of the proposed model is finally determined and will be used for image classification and parameter estimation. The second contribution of this paper is the derivation of an estimation algorithm associated with the proposed hierarchical Bayesian model. As the minimum mean square error (MMSE) and maximum a posteriori (MAP) estimators of this model cannot be easily computed from its joint posterior, we investigate a hybrid Gibbs sampler allowing the posterior of interest to be samples (see
[16, 17] for details). The proposed Bayesian model and estimation algorithm are validated using synthetic and real RCM images, resulting from a clinical study containing healthy and lentigo patients. The obtained results are very promising and show the potential of the proposed denoising and classification strategy.The paper is structured as follows. The classification problem studied in this work is introduced in Section II. The proposed hierarchical Bayesian model and its estimation algorithm are studied in Sections III and IV. Section VA validates the proposed technique using simulated data with different noise levels. Section VB shows results obtained using real data obtained from a clinical study. Conclusions and future work are finally reported in Section VI.
Ii Problem formulation
Iia Observation model
Consider noise free images, containing N pixels, gathered in the matrix , where , denotes the image associated with the th patient. Denote by the corresponding noisy images. Using these notations, the observation model is given by
(1) 
where and are () vectors representing the th observed and noiseless images, is a gamma noise () vector with a shape parameter and a scale parameter and denotes the termwise product. In order to ensure that the proposed model (1) is identifiable, the mean of the gamma noise is supposed to equal , leading to
(2) 
The problem addressed in this paper is to classify these images , into two classes representing healthy and lentigo patients. The next section introduces a hierarchical Bayesian model that is used for this classification.
Iii Hierarchical Bayesian model
This section introduces a hierarchical Bayesian model that can be used to estimate the unknown matrix of noiseless images , the vectors containing the class labels and the noise variances associated with the observed images from the matrix
. This model is defined by a likelihood, and by parameter and hyperparameter priors defined below.
Iiia Likelihood
The multiplicative speckle noise is known to have a gamma distribution. Thus, the observation model (1) leads to
(3) 
where means "is distributed according to",
is the gamma distribution whose probability density function (pdf) is
(4) 
with the indicator function on , means “proportional to” and denotes the gamma function. Assuming independence between the observed signals, the likelihood of the observed images can be written
IiiB Priors for the signal of interest
To ensure the positivity of the noiseless images, a truncated Gaussian distribution is assigned to for
(5) 
where
denotes the truncated normal distribution on
, takes the two values and depending on the patient class, and are the means and variance of the two truncated Gaussian distributions.IiiC Prior for the noise variances
A noninformative conjugate inverse gamma prior (denoted as ) is classically selected for the scale parameter [18]
(6) 
where and are fixed hyperparameters, that are adjusted to reflect the absence of prior knowledge on , i.e., the mean and variance of were fixed to and in order to obtain a flat prior. The joint prior for the vector of noise variances denoted as is finally obtained as the product of the marginal densities .
IiiD Prior for the label vector
The parameter vector is a label vector that associates each image to a healthy or lentigo skin. Because of the absence of prior knowledge about this parameter, it is assigned a uniform prior defined as
(7) 
The labels associated with the different patients are supposed to be a priori independent, i.e., the joint prior of denoted as f() is the product of the probabilities defined in (7).
IiiE Hyperparameter priors
In order to complete the description of the proposed hierarchical Bayesian model and to allow hyperparameters to be estimated directly from the data, we propose to assign priors for the different hyperparameters. A Gaussian prior has been selected for the mean and a noninformative inverse gamma prior for the variance (see [19, 18] for motivations)
(8) 
(9) 
where , , , are fixed in order to obtain flat priors, i.e., , whereas the mean and variance of were fixed to and . The joint pdfs and are finally obtained as the product of their marginal densities assuming prior independency between the components of these two vectors.
IiiF Joint posterior distribution
The proposed Bayesian model is illustrated by the directed acyclic graph (DAG) displayed in Fig. 2, which highlights the relationships between the observations , the parameters and the hyperparameters . Assuming prior independence between the different components of the parameter vector , the joint posterior distribution of this Bayesian model can be computed using the following hierarchical structure
(10) 
(11) 
The complexity of the proposed Bayesian model summarized in the DAG of Fig. 2 and its resulting posterior (10) prevent a simple computation of the maximum aposteriori (MAP) or minimum mean square (MSE) estimators of the unknown model parameters. The next section studies an MCMC method that is used to sample the posterior (10) and to build estimators of the parameters involved in the proposed Bayesian model using the generated samples.
Iv MetropoliswithinGibbs algorithm
This section studies a hybridGibbssampler, which is guaranteed to generate samples asymptotically distributed according to the target distribution (10). The Gibbs sampler described in Algo. 1, iteratively generates samples distributed according to the conditional distributions of (10). These conditional distributions are detailed in the rest of this section. Because of the complexity of the conditional distributions, we consider randomwalk MetropolisHastings (RWMH) [17, 16] moves within the Gibbs sampler, which consists in generating samples distributed according to the complex conditional distribution of each parameter of interest. This is achieved using the conditional distributions , for , and their associated proposal distributions . The first step is to initialize the sample value for each parameter , for . The main loop of the RWMH algorithm consists of three components:

Generate a candidate from the proposal distribution . The distribution is the truncated Gaussian distribution [20] for the parameters and the Gaussian distribution for .

Compute the acceptance probability using the acceptance function based upon the proposal distribution and the conditional density for each parameter

Accept the candidate with probability .
In order to maximize the efficiency of the algorithm, the variances of the proposal distributions have been adjusted such that the acceptance rate is between 0.3 and 0.6 as suggested in [16] and [21].
The conditional distributions of the parameters of interest are obtained by multiplying the likelihood with the different priors and by removing the multiplicative terms that do not depend on the variable of interest. The algorithm iteratively updates each parameter by using its conditional distribution detailed in the following paragraphs.
Iva Sampling the parameter :
(12) 
The proposal law is :
We generate , such that :
The acceptancerejection rule is :
with
IvB Sampling the parameter :
(13) 
The proposal law is :
where is the cumulative normal distribution function
We generate , such that :
The acceptancerejection rule is :
with
IvC Sampling the parameter :
(14) 
The proposal law is :
We generate , such that :
The acceptancerejection rule is :
with
IvD Sampling the parameter :
(15) 
The proposal law is :
We generate , such that :
The acceptancerejection rule is :
with
IvE Sampling the parameter :
(16) 
IvF Bayesian inference and parameter estimation
The main steps of the proposed MetropoliswithinGibbs sampler are summarized in Algo. 1. This algorithm provides a sequence of samples of the vector denoted as that are used to approximate the MMSE estimators by using Monte Carlo integration [22] as
(17) 
where is the number of burnin iterations and is the total number of Monte Carlo iterations. Finally, the following maximum aposteriori (MAP) estimator is considered for the label
(18) 
where and denote the numbers of samples satisfying the conditions and in the interval .
IvG Convergence:
Running multiple chains with different initializations allows to define various convergence measures for MCMC methods [23]. The popular betweenwithin variance criterion has shown interesting properties for diagnosing convergence of MCMC methods. This criterion was initially studied by Gelman and Rubin in [24] and has been used in many studies including [23, p. 33], [25] and [26]. The main idea is to run parallel chains of length + for each data set with different starting values and to evaluate the dispersion of the estimates obtained from the different chains. The betweensequence variance and withinsequence variance for the Markov chains are defined by
(19) 
(20) 
with
(21) 
where is the parameter of interest and is its estimate at the th run of th chain. The convergence of the chain can then be monitored by the socalled potential scale reduction factor defined as [27, p. 332]
(22) 
V Simulation results
Va Synthetic data
This section evaluates the performance of the proposed algorithm on synthetic data. Different experiments were conducted using three values of the signal to noise ratio , allowing the algorithm performance to be appreciated for different noise levels. This section considers synthetic images. Each image contains pixels and was generated according to (3). These images were separated into healthy and lentigo classes containing images. The noiseless images of the two classes were respectively generated according to the truncated Gaussian distributions and , with . The sampler convergence of the algorithm is monitored by computing the potential scale reduction factor introduced in (IVG) for an appropriate parameter of interest. Different choices for the parameter could be considered for the proposed method. This paper proposes to monitor the convergence of the MetropoliswithinGibbs sampler by checking the noise variance (see [19, 25] for a similar choice). The potential scale reduction factor for parameter computed for Markov chains is equal to 1.01. This value of confirms the good convergence of the sampler (a recommendation for convergence assessment is a value of [27, p. 332] ). Figs. 3, 4, 5 and 6 show the evolution of the Markov chains for the different parameters estimated for synthetic data with and the real data using the RCM images, respectively. Algo. 1 was run for iterations and the different model parameters were estimated using (17) and (18) using a burnin period of length . The performance of the algorithm was evaluated by computing the mean square errors (MSEs) of the different parameters and the signal to noise ratios (SNRs) defined as
(23) 
(24) 
Quantitative results are presented in Table I for the three experiments. This table shows good estimation results of the parameters when considering different noise levels. This table also shows excellent classification results for dB, and when considering the challenging case dB. These results highlight the potential of the proposed strategy in denoising and classifying the images obtained from model (3) and improving the estimation of the different parameters of this model.
SNR (dB)  SNR (dB)  SNR (dB)  
0.56  30.12  62.72  70.4  
0.95  21.42  73.24  67.79  
2.91  1.01  0.015  18.07  0.011  25.7  
7.14  2.57  4.58  5.42  0.006  22.07  
20.44  26.56  30.44  
5.48  16.53  2.88  20.81  0.7093  26.87  
Accuracy  91  100  100  
Accuracy (CART)  83  100  100 
VB Real data
This section is devoted to the validation of the proposed denoising and classification algorithm when applied to real RCM images. These RCM images were acquired with apparatus Vivascope and correspond to the stratum corneum, the epidermis layer, the dermisepidermis junction (DEJ) and the upper papillary dermis. Each RCM image shows a field of view with pixels. A set of women aged years and over were recruited. All the volunteers gave their informed consent for examination of skin by RCM. According to the clinical evaluation performed by a physician, volunteers were divided into two groups. The first group was formed by women with at least lentigines on the back of the hand whereas women without lentigo constituted the control group. Images were taken on lentigo lesions for volunteers of the first group and on healthy skin on the back of the hand for the control group. Consequently, our database contained = 45 patients. An examination of each acquisition was performed in order to locate the stratum corneum and the DEJ precisely in each image. Given the large size of the images, we preferred to select and apply our algorithm to patches of
pixels for each image to reduce the computational cost. The obtained results were then used to calculate the confusion matrix and four indicators (sensitivity, specificity, precision, accuracy) shown in Tables
II and III. These indicators are defined as Sensitivity = TP/(TP+FN), Specificity = TN/(FP+TN), Precision = TP/(TP+FP), Accuracy = (TP+TN)/(TP+FN+FP+TN), where TP, TN, FP and FN are the numbers of true positives, true negatives, false positives and false negatives. This table allows us to evaluate the classification performance of the proposed strategy. The accuracy of the proposed method equals , which corresponds to a single mistake for the lentigo patient . Fig. 7 shows that the texture of this misclassified image is not very destructed as for other lentigo patients, and is visually similar to the texture of healthy patients. Fig. 8 shows examples of noisy RCM images and their estimated true reflectivity, illustrating the denoising part of the proposed algorithm. We can observe that the estimated images have low intensities compared to the noisy images which is due to the fact that the noise is multiplicative. To assess the significance of our results, our algorithm was then compared to the method presented in [15]. This method consists in extracting from each RCM image a set of 39 analysis parameters (further technical details are available in [28]) and to apply to these features a classification procedure based on classification and regression trees (CART). Note that the CART algorithm was tested on the real RCM images using a leave one out procedure. As shown in Table III, the accuracy obtained with the CART algorithm is , i.e., it is slightly smaller that the one obtained with the proposed method. Moreover, the proposed Bayesian model can be used for the characterization of RCM images thanks to its estimated parameters.Confusion matrix 


Lentigo  26  1  96.2  
Healthy  0  18  100  
Precision  100  94.7  
Accuracy  97.7 
Confusion matrix 


Lentigo  24  3  88.8  
Healthy  5  13  72.2  
Precision  82.7  81.2  
Accuracy  82.2 
Vi Conclusions
This paper presented a new Bayesian strategy as well as an MCMC algorithm for classifying RCM images as healthy or lentigo images. A Bayesian model was introduced based on a gamma distribution for the multiplicative speckle noise and on various priors assigned to the unknown model parameters. A hybrid Gibbs sampler was then considered to sample the posterior of this Bayesian model and to build Bayesian estimators. Simulation results conducted on synthetic and real data allowed the good performance of the proposed classifier to be appreciated. Future work includes the introduction of spatial correlation on the estimated noiseless images to improve their quality.
References
 [1] T. D. Menge, B. P. Hibler, M. A. Cordova, K. S. Nehal, and A. M. Rossi, “Concordance of handheld reflectance confocal microscopy (rcm) with histopathology in the diagnosis of lentigo maligna (lm): A prospective study,” Journal of the American Academy of Dermatology, vol. 74, no. 6, pp. 1114–1120, 2016.
 [2] K. S. Nehal, D. Gareau, and M. Rajadhyaksha, “Skin imaging with reflectance confocal microscopy,” Seminars in Cutaneous Medecine and Surgery, Elsevier, vol. 27, pp. 37–43, 2008.
 [3] R. HofmannWellenhof, E. M. Wurm, V. AhlgrimmSiess, S. K. E. Richtig, J. Smolle, and A. Gerger, “Reflectance confocal microscopystateofart and research overview,” Seminars in Cutaneous Medecine and Surgery, Elsevier, vol. 28, pp. 172–179, 2009.
 [4] I. Alarcon, C. Carrera, J. Palou, L. Alos, J. Malvehy, and S. Puig, “Impact of in vivo reflectance confocal microscopy on the number needed to treat melanoma in doubtful lesions,” British journal of Dermatology., vol. 170, pp. 802–808, 2014.
 [5] I. Alarcon, C. Carrera, L. Alos, J. Palou, J. Malvehy, and S. Puig, “In vivo reflectance confocal microscopy to monitor the response of lentigo maligna to imiquimod,” J Am Acd Dermatol., vol. 71, pp. 49–55, 2014.
 [6] P. Guitera, L. Haydu, S. Menzies, and al, “Surveillance for treatment failure of lentigo maligna with dermoscopy and in vivo confocal microscopy: new descriptors.” British journal of Dermatology., vol. 170, pp. 1305–1312, 2014.
 [7] J. Champin, J.L. Perrot, E. Cinotti, B. Labeille, C. Douchet, G. Parrau, F. Cambazard, P. Seguin, and T. Alix, “In vivo reflectance confocal microscopy to optimize the spaghetti technique for defining surgical margins of lentigo maligna,” Dermatologic Surgery, vol. 40, no. 3, pp. 247–256, 2014.
 [8] B. Hibler, M. Cordova, R. Wong, and A. R. ., “Intraoperative realtime reflectance confocal microscopy for guiding surgical margins of lentigo maligna melanoma,” Dermatol Surg., vol. 41, pp. 980–983, 2015.
 [9] B. L. Luck, K. D. Carlson, A. C. Bovik, and R. R. RichardsKortum, “An image model and segmentation algorithm for reflectance confocal images of in vivo cervical tissue,” IEEE Trans. Image Processing, vol. 14, no. 9, pp. 1265–1276, 2005.
 [10] S. Kurugol, J. G. Dy, M. Rajadhyaksha, K. W. Gossage, J. Weissmann, and D. H. Brooks, “Semiautomated algorithm for localization of dermal/epidermal junction in reflectance confocal microscopy images of human skin,” in SPIE BiOS. International Society for Optics and Photonics, 2011, pp. 79 041A–79 041A.
 [11] S. Kurugol, M. Rajadhyaksha, J. G. Dy, and D. H. Brooks, “Validation study of automated dermal/epidermal junction localization algorithm in reflectance confocal microscopy images of skin,” in SPIE BiOS. International Society for Optics and Photonics, 2012, pp. 820 702–820 702.
 [12] S. C. Hames, M. Ardigo, H. P. Soyer, A. P. Bradley, and T. W. Prow, “Anatomical skin segmentation in reflectance confocal microscopy with weak labels,” in Digital image computing: techniques and applications (dICTA), 2015 international conference on. IEEE, 2015, pp. 1–8.
 [13] S. C. Hames, M. Ardigò, H. P. Soyer, A. P. Bradley, and T. W. Prow, “Automated segmentation of skin strata in reflectance confocal microscopy depth stacks,” PloS one, vol. 11, no. 4, p. e0153208, 2016.

[14]
K. Kose, C. AlessiFox, M. Gill, J. G. Dy, D. H. Brooks, and M. Rajadhyaksha, “A machine learning method for identifying morphological patterns in reflectance confocal microscopy mosaics of melanocytic skin lesions invivo,” in
SPIE BiOS. International Society for Optics and Photonics, 2016, pp. 968 908–968 908.  [15] S. Koller, M. Wiltgen, V. AhlgrimmSiess, W. Weger, R. HofmannWellenhof, E. Richtig, J. Smolle, and A. Gerger, “In vivo reflectance confocal microscopy: automated diagnostic image analysis of melanocytic skin tumours,” Journal of the European Academy of Dermatology and Venereology, vol. 25, no. 5, pp. 554–558, 2011.
 [16] S. Chib and E. Greenberg, “Understanding the MetropolisHastings algorithm,” The American Statistician, vol. 49, no. 4, pp. 327–335, 1995.
 [17] W. R. Gilks, S. Richardson, and D. Spiegelhalter, Markov chain Monte Carlo in practice. CRC press, 1995.

[18]
D. Fink, “A compendium of conjugate priors,”
See http://www. people. cornell. edu/pages/df36/CONJINTRnew% 20TEX. pdf, p. 46, 1997.  [19] N. Dobigeon, J.Y. Tourneret, and C.I Chang, “Semisupervised linear spectral unmixing using a hierarchical Bayesian model for hyperspectral imagery,” IEEE Trans. Signal Process., vol. 56, no. 7, pp. 2684–2695, July 2008.
 [20] C. P. Robert, “Simulation of truncated normal variables,” Statistics and computing, vol. 5, no. 2, pp. 121–125, 1995.
 [21] G. O. Roberts, J. S. Rosenthal, et al., “Optimal scaling for various MetropolisHastings algorithms,” Statistical science, vol. 16, no. 4, pp. 351–367, 2001.
 [22] C. P. Robert and G. Casella, Monte Carlo Statistical Methods. New York: SpringerVerlag, 1999.
 [23] C. P. Robert and D. Cellier, “Convergence control of MCMC algorithms,” in Discretization and MCMC Convergence Assessment, C. P. Robert, Ed. New York: Springer Verlag, 1998, pp. 27–46.
 [24] A. Gelman and D. Rubin, “Inference from iterative simulation using multiple sequences,” Statistical Science, vol. 7, no. 4, pp. 457–511, 1992.
 [25] S. J. Godsill and P. J. W. Rayner, “Statistical reconstruction and analysis of autoregressive signals in impulsive noise using the Gibbs sampler,” IEEE Trans. Speech, Audio Processing, vol. 6, no. 4, pp. 352–372, 1998.

[26]
P. M. Djurić and J.H. Chun, “An MCMC sampling approach to estimation of nonstationary hidden markov models,”
IEEE Trans. Signal Process., vol. 50, no. 5, pp. 1113–1123, 2002.  [27] A. Gelman, J. B. Carlin, H. S. Stern, and D. B. Rubin, Bayesian Data Analysis. London: Chapman & Hall, 1995.
 [28] M. Wiltgen, A. Gerger, C. Wagner, J. Smolle, et al., “Automatic identification of diagnostic significant regions in confocal laser scanning microscopy of melanocytic skin tumors,” Methods of Information in Medicine, vol. 47, no. 1, pp. 14–25, 2008.
Comments
There are no comments yet.