An unsupervised bayesian approach for the joint reconstruction and classification of cutaneous reflectance confocal microscopy images

03/04/2017 ∙ by Abdelghafour Halimi, et al. ∙ 0

This paper studies a new Bayesian algorithm for the joint reconstruction and classification of reflectance confocal microscopy (RCM) images, with application to the identification of human skin lentigo. The proposed Bayesian approach takes advantage of the distribution of the multiplicative speckle noise affecting the true reflectivity of these images and of appropriate priors for the unknown model parameters. A Markov chain Monte Carlo (MCMC) algorithm is proposed to jointly estimate the model parameters and the image of true reflectivity while classifying images according to the distribution of their reflectivity. Precisely, a Metropolis-whitin-Gibbs sampler is investigated to sample the posterior distribution of the Bayesian model associated with RCM images and to build estimators of its parameters, including labels indicating the class of each RCM image. The resulting algorithm is applied to synthetic data and to real images from a clinical study containing healthy and lentigo patients.



There are no comments yet.


page 3

page 19

page 20

page 21

page 23

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

The lentigo is a hyperplasia that affects the skin. It comes from the proliferation of melanocyte cells at the dermo-epidermic 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 follow-up [5, 6, 7] and guidance [8] have also been published.

Fig. 1: Images (at the depth 49.5 ) from healthy (patient , , , , , ) and lentigo (patient , , , , , ) patients at the DEJ depth. One can observe more textured images in the presence of lentigo.

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 wavelet-based 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 V-A validates the proposed technique using simulated data with different noise levels. Section V-B shows results obtained using real data obtained from a clinical study. Conclusions and future work are finally reported in Section VI.

Ii Problem formulation

Ii-a 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


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


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.

Iii-a Likelihood

The multiplicative speckle noise is known to have a gamma distribution. Thus, the observation model (1) leads to


where means "is distributed according to",

is the gamma distribution whose probability density function (pdf) is


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

Iii-B Priors for the signal of interest

To ensure the positivity of the noiseless images, a truncated Gaussian distribution is assigned to for



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.

Iii-C Prior for the noise variances

A non-informative conjugate inverse gamma prior (denoted as ) is classically selected for the scale parameter [18]


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 .

Iii-D 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


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).

Iii-E 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 non-informative inverse gamma prior for the variance (see [19, 18] for motivations)


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.

Iii-F 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


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 a-posteriori (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.

Fig. 2: DAG for the parameter and hyperparameter priors. The user fixed hyperparameters appear in boxes (continuous line).

Iv Metropolis-within-Gibbs algorithm

This section studies a hybrid-Gibbs-sampler, 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 random-walk Metropolis-Hastings (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:

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

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

  3. 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].

1:  Input:
2:  Initialization
3:  Initialize
4:  for =1 to  do
5:     Parameter update
6:     Sample according to (12) using an RWMH with a truncated Gaussian proposal
7:     Sample according to (13) using an RWMH with a truncated Gaussian proposal
8:     Sample according to (14) using an RWMH with a Gaussian proposal
9:     Sample according to (15) using an RWMH with a truncated Gaussian proposal
10:     Sample from the pdf (16)
11:  end for
12:  Result: for .
Algorithm 1 Metropolis-within-Gibbs algorithm

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.

Iv-a Sampling the parameter :


The proposal law is :

We generate , such that :

The acceptance-rejection rule is :


Iv-B Sampling the parameter :


The proposal law is :

where is the cumulative normal distribution function

We generate , such that :

The acceptance-rejection rule is :


Iv-C Sampling the parameter :


The proposal law is :

We generate , such that :

The acceptance-rejection rule is :


Iv-D Sampling the parameter :


The proposal law is :

We generate , such that :

The acceptance-rejection rule is :


Iv-E Sampling the parameter :


Iv-F Bayesian inference and parameter estimation

The main steps of the proposed Metropolis-within-Gibbs 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


where is the number of burn-in iterations and is the total number of Monte Carlo iterations. Finally, the following maximum a-posteriori (MAP) estimator is considered for the label


where and denote the numbers of samples satisfying the conditions and in the interval .

Iv-G Convergence:

Running multiple chains with different initializations allows to define various convergence measures for MCMC methods [23]. The popular between-within 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 between-sequence variance and within-sequence variance for the Markov chains are defined by




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 so-called potential scale reduction factor defined as [27, p. 332]


V Simulation results

V-a 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 (IV-G) 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 Metropolis-within-Gibbs 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 burn-in 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


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
TABLE I: Performance of the proposed algorithm for denoising and classification of synthetic data for three corrupted data

V-B 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 dermis-epidermis 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 mis-classified 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
TABLE II: Classification performance on real data (45 patients) using the proposed method.
Confusion matrix
Lentigo 24 3 88.8
Healthy 5 13 72.2
Precision 82.7 81.2
Accuracy 82.2
TABLE III: Classification performance on real data (45 patients) using the CART method.
Fig. 3: Evolution of the convergence of the Markov chains for the different parameters estimated for the synthetic data with .
Fig. 4: Evolution of the convergence of the Markov chains for the different parameters estimated for the synthetic data with .
Fig. 5: Evolution of the convergence of the Markov chains for the different parameters estimated for the synthetic data with .
Fig. 6: Evolution of the convergence of the Markov chains for the different parameters estimated for the real RCM images.
Fig. 7: Images (at the depth 49.5 ) from the patient who is badly classified compared to a healthy and lentigo patient (well classified). One can observe more similarity between this patient and the healthy one then with the lentigo.
Fig. 8: Examples of noisy images (at the depth 49.5 ) and their estimated true reflectivities.

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.


  • [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. Hofmann-Wellenhof, E. M. Wurm, V. Ahlgrimm-Siess, S. K. E. Richtig, J. Smolle, and A. Gerger, “Reflectance confocal microscopy-state-of-art 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 real-time 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. Richards-Kortum, “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, “Semi-automated 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. Alessi-Fox, 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 in-vivo,” in

    SPIE BiOS.   International Society for Optics and Photonics, 2016, pp. 968 908–968 908.
  • [15] S. Koller, M. Wiltgen, V. Ahlgrimm-Siess, W. Weger, R. Hofmann-Wellenhof, 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 Metropolis-Hastings 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, “Semi-supervised 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 Metropolis-Hastings algorithms,” Statistical science, vol. 16, no. 4, pp. 351–367, 2001.
  • [22] C. P. Robert and G. Casella, Monte Carlo Statistical Methods.   New York: Springer-Verlag, 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.