Computerized Medical Imaging and Graphics 42 (2015) 51–55
Contents lists available at ScienceDirect
Computerized Medical Imaging and Graphics journal homepage: www.elsevier.com/locate/compmedimag
Heterogeneity assessment of histological tissue sections in whole slide images夽 Philippe Belhomme a,∗ , Simon Toralba a,c , Benoît Plancoulaine a , Myriam Oger a,b , Metin N. Gurcan c , Catherine Bor-Angelier a,b a
Normandie Université; UNICAEN, CLCC F. Baclesse, PATHIMAGE BioTICLA EA 4656, Caen, France Pathology Department, CLCC F. Baclesse, Caen, France c CIALab, Department of Biomedical Informatics, OSU, Columbus, OH, USA b
a r t i c l e
i n f o
Article history: Received 2 April 2014 Received in revised form 10 October 2014 Accepted 10 November 2014 Keywords: Heterogeneity Whole slide image Breast cancer Dimensionality reduction Spectral graph theory
a b s t r a c t Computerized image analysis (IA) can provide quantitative and repeatable object measurements by means of methods such as segmentation, indexation, classification, etc. Embedded in reliable automated systems, IA could help pathologists in their daily work and thus contribute to more accurate determination of prognostic histological factors on whole slide images. One of the key concept pathologists want to dispose of now is a numerical estimation of heterogeneity. In this study, the objective is to propose a general framework based on the diffusion maps technique for measuring tissue heterogeneity in whole slide images and to apply this methodology on breast cancer histopathology digital images. © 2014 Elsevier Ltd. All rights reserved.
1. Introduction Research in signal and image analysis is going on for many decades now and is directly linked with the exceptional development of computer technologies. But after all these years, it must be admitted that there are not so many real working applications in practice, especially in the medicine area where the expert’s eye is still more accurate and faster than many automated systems dealing with large amounts of data. However, reliable automated systems could really help pathologists in their daily work as the number of pathological cases increases as far as the early screening campaigns do. We know that computerized image analysis (IA) can provide quantitative and repeatable object measurements by means of methods such as segmentation, indexation or classification. When IA is used to analyze images of histological sections in the medical research, the typical objects to be processed are nuclei, vessels, cell groups or tumors. IA allows these structures in histopathological images to be detected, to be analyzed automatically in terms of their size or shape, in order to assess their proportion or to compute their staining intensities. Indeed, IA can
夽 Special thanks to the Pathology Department, Caen, France & National Center of Pathology, Vilnius, Lithuania. ∗ Corresponding author. Tel.: +33 2 33 01 46 29. E-mail address: [email protected] (P. Belhomme). http://dx.doi.org/10.1016/j.compmedimag.2014.11.006 0895-6111/© 2014 Elsevier Ltd. All rights reserved.
contribute to more accurate determination of prognostic histological factors by pathologists [1] but its action is limited to provide quantitative information as no multi-purpose method exists for segmenting or classifying images with a wide range of color or shape configurations, with which pathologists are familiar. The advent of digital scanners leads to generate whole slide images (WSI) from histological sections acquired at a full resolution. A main challenge this technology presents is the size of data to be processed which are generally very large and requires a huge storage capacity and processing needs. Conversely, the main advantage is that the complete tissue structure is easily accessed. In parallel, the aggressiveness of a cancer could result in morphological and architectural changes that can be observed in the tissue structure, and so be characterized by the object distribution on the slide, by the cross relations between objects and by the texture. This kind of information could contribute to evaluate a well-known concept: heterogeneity. Frequently addressed in signal processing, especially in terms of “entropy”, but more rarely in the field of imaging [2], the objective here is to propose a framework for measuring tissue heterogeneity in WSI and apply this methodology on breast cancer histopathology digital images. The key idea in this work is to not rely on segmentation (e.g. [1,2]) of individual structures to characterize heterogeneity, but to make use of classification of squared sub-images later called ‘patches’. In some previous works dedicated to the development of a computer-aided diagnosis system (CADS) based on image retrieval and classification [3,5], we have used a
52
P. Belhomme et al. / Computerized Medical Imaging and Graphics 42 (2015) 51–55
method coming from spectral graph theory, the Diffusion Maps (DM) [6], to process WSI split in small squares called ‘patches.’ The DM algorithm, in which eigenvalues and eigenvectors of a Markov matrix defining a random walk on the data are computed, allows to both cluster non-linear input data thanks to its inner classification properties preserving local neighborhood relationships, but also to reduce the input data dimensionality in a space (usually a 2D or 3D space) where it is therefore possible to compute euclidean distances between the objects to be analyzed [6,8]. To briefly describe the complete CADS we are developing, a first step consists in building a knowledge database involving many features extracted from a set of well-known images; this is an ‘off-line’ procedure conducted once. These features are represented by vectors of non-linear data acting as a signature. In a second step, signatures are obtained from new unknown images and then compared with those already present in the database; this is an ‘on-line’ procedure that has to be conducted each time a new image is processed. In a last step, before rendering a qualitative measure of similarity/dissimilarity between the supervised images and new unknown images, a feedback procedure would have to be processed in order to eliminate most of the artifacts that usually corrupt initial knowledge databases. But this general approach, especially with the DM algorithm, can also be derived to analyze a set of image patches coming from any WSI, with just the goal to compare their feature vectors. In this paper, we focus on a way to characterize the tissue heterogeneity at a regional level by working on the projected coordinates of image patches in a reduced 3D space. 2. Materials and methods Images used to illustrate this study come either from histological sections or tissue microarrays (TMAs) of breast cancer stained according to the Ki67 protocol and the Hematoxylin–Eosin-Saffron protocol (HES). They are acquired at a 20× magnification on a microscopic scanner device (ScanScope CS; Aperio Technologies). The resolution of WSI is 0.5 microns/pixel and the typical image size of any histological section can reach 50,000 × 40,000 pixels, corresponding to an uncompressed file of size 5.6 Gb. They are stored in the TIFF 6.0 file format (compression 30%). In each WSI some regions of interest (ROI), chosen by the pathologist, are then split in connected blocks of size 50 × 50 pixels, called “patches”, and each patch is associated with a numerical signature expressed as a feature vector. Tools developed here are written in Python language with the help of specialized modules (PIL: Python Imaging Library, SciPy-Numpy, Matplotlib, Mahotas. . .). 2.1. Features extraction The feature vectors embed some statistical measures obtained from color components and some texture parameters. In this study, 61 numerical values per component were calculated from up to two color spaces (RGB for Ki67 images and RGB + H&E color deconvolution for HES images [11]) plus 18 values obtained from the excess-red component (2R-G-B). From any given component, the 9 statistical features are the mean, median, mode, Skewness, Kurtosis parameters with also the 20%–40%–60%–80% quantiles of its cumulated histogram (initially reduced to 64 values). The 52 texture features correspond to the classic 13 Haralick parameters obtained in four directions [12]. We also added to the texture features 18 new values coming from the intrinsic statistical parameters of regionalized variables in three directions (0◦ –45◦ –90◦ ). They are derived from the nugget, sill and range of a geostatistical method [13]. Finally, each feature vector thus contains either F = 201 (61 × 3 + 18) or F = 323 (61 × 5 + 18) numerical values, here with the predominance of texture features in order to be quite independent from the
color staining variations encountered between different laboratories (and often inside the same laboratory). In order to later compare feature vectors, and considering the sparse numerical range of their values, the symmetric Kullback–Leibler distance [3,4] has been retained for its ability to easily manage such a case, while remaining fast to implement. The distance between two vectors p1 ,p2 of length F is given by: 1 DKL (p1 , p2 ) = p1j · log 2 F
j=1
p1j p2j
+ p2j
p2j p1j
(1)
2.2. Dimensionality reduction In any classical CADS, one of the key components is a visualization tool showing relationships between supervised images, stored in a knowledge database, and new images that are presented to the system. Typically, these relationships may be expressed as a connected graph in a 2D or 3D space where one hopes to find distinctive clusters corresponding to histological types or sub-types. It is therefore mandatory to reduce dimensionality from F (201 or 323 in our application) to 2 or 3 dimensions. With feature vectors containing non linear data as we are faced with, authors in [7,8] have shown that it was not appropriate to perform a principal component analysis (PCA). Instead, methods based on Spectral Connectivity Analysis (SCA) such as the Diffusion Maps, involving eigenvalues and eigenvectors of a normalized graph Laplacian, are well suited for this task. Let X = {x1 ,x2 ,. . .,xn } be a set of n patches that we assimilate to a fully connected graph G, that means a distance function is computed for each pair {xi ,xj }. A n × n kernel P is obtained from a Gaussian function whose coefficients are given by: p(xi , xj ) =
w(xi , xj )
(2)
d(xi )
with d(xi ) =
w(xi , xk )
(3)
xk ∈ X
and
D
w(xi , xj ) = e−
KL (xi ,xj ) ∈
(4)
In fact, p(xi ,xj ) may be considered as the transition kernel of the Markov chain on G. In other words, p(xi ,xj ) defines the transition probability for going from xi to xj in one time step. The eigenvectors k of P, ordered by decreasing positive eigenvalues, give the practical observation space axes. It must be noticed that 0 is never used since linked to the eigenvalue = 1 (i.e. the data set mean or trivial solution). A 3D projection can be then obtained along axes (1 ,2 ,3 ). Choosing ε in w(xi ,xj ) is an empirical task which should permit a moderate decrease of the exponential in Eq. (4); some works [8] use the median value of all DKL (xi ,xj ) distances whereas other works [6] use the mean distance obtained in the k nearest neighbors from a subset of X, which finally yields to quite the median value in many cases. We have retained the first solution for its simplicity and to lessen the overall computational time. Fig. 1a–c show the resulting projections obtained from the DM algorithm for 300 patches of fibroadenoma with ε = 0.5 × medianValue, ε = medianValue and ε = 2 × medianValue. One can see that the general shapes of point clouds are not the same but, due to the inner classification properties of DM, we have shown that the local neighborhood relationships were preserved in terms of histological types/sub-types [9].
P. Belhomme et al. / Computerized Medical Imaging and Graphics 42 (2015) 51–55
53
2.3. Pseudo-color space At the end of the DM algorithm, the three main eigenvectors of matrix P provide an uncorrelated space in which all patches are represented by a point (this space is thereafter called “Diffusion Maps” cube or DM cube). The Diffusion Maps method, belonging to the field of unsupervised learning algorithms, has the interesting property to express the complex structure of a multidimensional data set (also called manifold) in a reduced space where the Euclidean distance can be used, while preserving initial local neighborhood relationships. Thus, patches that seem to be similar in the WSI will be projected in the same area within the DM cube. In order to ease displaying the projected patches, the DM cube is considered as a color space in which each patch is assigned to a pseudo-color corresponding to its physical location. To mimic human visual system, the first main eigenvector is taken for the green component, the next two eigenvectors respectively for the red and blue components. Therefore, estimating the heterogeneity of histological sections comes down to measure pseudo-color differences of their patches in the DM cube. 2.4. Heterogeneity measurement In this paper, three scale levels for rendering a heterogeneity measurement are explored. The first level is local and computes a color gradient image in 4-connexity, which brings to light spatial variations of pattern structures in tissue sections. The result image is a classic gradient image with low values where patch signatures are slightly different between neighbors, and higher values on pattern transitions. This allows separating areas such as stroma and epithelium. The heterogeneity measurement can then be obtained by computing the average of all the gradient values, associated with their standard deviation. In Fig. 4a the mean value is expressed. The second level proposed is regional as it provides a linear measurement of heterogeneity along Bresenham lines plotted in the WSI. For a line LBj whose extremities correspond to patches A and B and that includes a series of intermediate patches Pi , the algorithm [2] adds the color differences of pairwise connected patches Pi and Pi+1 along LBj , the sum being finally divided by LBj length. If changes in patch signatures along a Bresenham line are regular and proportional to the length, then patches should rather be aligned in the DM cube and H(LBj ) should be low. To ensure a relevant measure of heterogeneity, multiple Bresenham lines with variable lengths and directions are randomly drawn in the WSI. In Fig. 4b the mean value of all H(LBj ) is expressed. Finally, the third level proposed to compute the heterogeneity is global. It involves computing the Rao’s quadratic entropy (RQE) [10] from N patches according to the following formula: RQE =
N N−1
dij pi pj
(5)
i=1 j=i+1
where dij is the pseudo-color difference between two patches and pi ,pj are the probabilities for patches to be present, set up here such as pi = pj = 1/N. RQE usually quantifies the general disorder of a population as the average difference between individuals; in our case the population is represented by image patches. Fig. 1. (a) Projection of 300 patches in a 3D space with ε = 0.5 × medianValue. Colors have no sense here except to represent groups of 100 patches. (b) Projection of 300 patches in a 3D space with ε = medianValue. (c) Projection of 300 patches in a 3D space with ε = 2 × medianValue.
3. Results and discussion To simplify the displaying of our tissue heterogeneity quantification method in this abstract, tissue microarrays (TMAs) and sub-images extracted from WSIs of breast cancer are used. Fig. 2 shows four TMAs with their corresponding pseudo-color map in the DM cube. The homogeneous parts are visible and are separated
54
P. Belhomme et al. / Computerized Medical Imaging and Graphics 42 (2015) 51–55
Fig. 2. TMAs from virtual slides of breast cancer and their corresponding pseudo-color map.
Fig. 3. Four sub-images qualified as slightly heterogeneous (1, 2), greatly (3) and moderately (4).
by strong transitions in the tissue structure. On Fig. 3, four subimages of size 3000 × 2000 pixels were selected and qualitatively evaluated by a pathologist. Fig. 3(1) and Fig. 3(2), that have a similar aspect, are qualified as slightly heterogeneous, while Fig. 3(3) is said greatly and Fig. 3(4) moderately heterogeneous. Once the
100
10
1
0,8
25
Quadratic Rao's entropy
Regional Brensenham lines
Mean Color Gradient
1000
heterogeneity measure has been computed according to the three methods above mentioned (locally with the color gradient, regionally with Bresenham lines and globally with the Rao’s quadratic entropy), the returned relative values correlate with the qualitative pathologist evaluation (Fig. 4). In this study, we have exposed
20
15
10
5
0 Image 1
Image 2
Image 3
Image 4
0,7 0,6 0,5 0,4 0,3 0,2 0,1 0
Image 1
Image 2
Image 3
Image 4
Image 1
Image 2
Image 3
Image 4
Fig. 4. (a–c) Tissue heterogeneity measurement according to three scale levels. (a) Local with the mean value of a 4-connexity color gradient; (b) regional with the mean H(LBj ) values of randomly selected Bresenham lines; (c) global with the Rao’s quadratic entropy measure computed from the Diffusion Map cube.
P. Belhomme et al. / Computerized Medical Imaging and Graphics 42 (2015) 51–55
the ability to assess tissue heterogeneity in a WSI using a set of three different measurements. An interesting contribution of our approach is that we do not base the image processing task on object segmentation, which is a complex and very time-consuming procedure with such large images, but rather work on feature vectors that are extracted from patches. The tuning parameters then only correspond to the patch size and the features used in signatures which must be adapted to the input image. In our experience, texture features might be preferred to color features, especially for comparing images coming from different laboratories. Indeed, a wide range of color difference can be encountered, even with histological slides that have been automatically prepared according to the same staining protocol. That is why the number of texture parameters retained is higher than the number of statistical parameters (52 vs 9 per color component, with moreover 18 parameters per patch coming from the geostatistical method). 4. Conclusion Patch analysis, while it may appear as a basic approach, allows segmentation of areas with similar configurations, such as stroma or tumor epithelium, and thus to encode the spatial organization of the same structure within the entire image. In future works, with the goal to decrease the computational burden generated by the “connected patches” approach, we plan to work on sampled patches obtained from stereological regular grids.
55
References [1] Gurcan M, Boucheron L, Can A, Madabhushi A, Rajpoot N, Yener B. Histopathological image analysis: a review. IEEE Rev Biomed Eng 2009;2: 147–71. [2] Brooks FJ, Grigsby PW. Quantification of heterogeneity observed in medical images. BMC Med Imaging 2013;13:7. [3] Oger M, Belhomme P, Klossa J, Michels JJ, Elmoataz A. Automated region of interest retrieval and classification using spectral analysis. Diagn Pathol 2008;3(Suppl. 1):S17. [4] Kullback S, Leibler RA. On information and sufficiency. Ann Math Stat 1951;22(1):79–86. [5] Belhomme P, Oger M, Michels JJ, Plancoulaine B, Herlin P. Towards a computer aided diagnosis system dedicated to virtual microscopy based on stereology sampling and diffusion maps. Diagn Pathol 2011;6(Suppl. 1):S3. [6] Coifman RR, Lafon S. Diffusion maps. Appl Comput Harmonic Anal 2006;21(1):5–30. [7] Coifman RR, Lafon S, Lee AB, Maggioni M, Nadler B, Warner F, et al. Geometric diffusions as a tool for harmonics analysis and structure definition of data: diffusion maps. Proc Natl Acad Sci 2005;102(21):7426–31. [8] Belkin M, Niyogi P. Laplacian eigen maps for dimensionality reduction and data representation. Neural Comput 2003;15:1373–96. [9] Belhomme P, Oger M, Michels JJ, Plancoulaine B. Out-of-sample extension of diffusion maps in a computer aided diagnosis system. Application to breast cancer virtual slide images. Diagn Pathol 2013;8(Suppl. 1): S9. [10] Rao CR. Diversity and dissimilarity coefficients: a unified approach. Theor Pop Biol 1982;21:24–43. [11] Ruifrok AC, Johnston DA. Quantification of histochemical staining by color deconvolution. Anal Quant Cytol Histol 2001;23:291–9. [12] Haralick RM, Shanmugam K, Dinstein I. Textural features for image classification. IEEE Trans Syst Man Cybern 1973;SMC-3:610–21. [13] Cressie NAC. Statistics for Spatial Data. New York: John Wiley & Sons; 1993.
Contents lists available at ScienceDirect
Computerized Medical Imaging and Graphics journal homepage: www.elsevier.com/locate/compmedimag
Heterogeneity assessment of histological tissue sections in whole slide images夽 Philippe Belhomme a,∗ , Simon Toralba a,c , Benoît Plancoulaine a , Myriam Oger a,b , Metin N. Gurcan c , Catherine Bor-Angelier a,b a
Normandie Université; UNICAEN, CLCC F. Baclesse, PATHIMAGE BioTICLA EA 4656, Caen, France Pathology Department, CLCC F. Baclesse, Caen, France c CIALab, Department of Biomedical Informatics, OSU, Columbus, OH, USA b
a r t i c l e
i n f o
Article history: Received 2 April 2014 Received in revised form 10 October 2014 Accepted 10 November 2014 Keywords: Heterogeneity Whole slide image Breast cancer Dimensionality reduction Spectral graph theory
a b s t r a c t Computerized image analysis (IA) can provide quantitative and repeatable object measurements by means of methods such as segmentation, indexation, classification, etc. Embedded in reliable automated systems, IA could help pathologists in their daily work and thus contribute to more accurate determination of prognostic histological factors on whole slide images. One of the key concept pathologists want to dispose of now is a numerical estimation of heterogeneity. In this study, the objective is to propose a general framework based on the diffusion maps technique for measuring tissue heterogeneity in whole slide images and to apply this methodology on breast cancer histopathology digital images. © 2014 Elsevier Ltd. All rights reserved.
1. Introduction Research in signal and image analysis is going on for many decades now and is directly linked with the exceptional development of computer technologies. But after all these years, it must be admitted that there are not so many real working applications in practice, especially in the medicine area where the expert’s eye is still more accurate and faster than many automated systems dealing with large amounts of data. However, reliable automated systems could really help pathologists in their daily work as the number of pathological cases increases as far as the early screening campaigns do. We know that computerized image analysis (IA) can provide quantitative and repeatable object measurements by means of methods such as segmentation, indexation or classification. When IA is used to analyze images of histological sections in the medical research, the typical objects to be processed are nuclei, vessels, cell groups or tumors. IA allows these structures in histopathological images to be detected, to be analyzed automatically in terms of their size or shape, in order to assess their proportion or to compute their staining intensities. Indeed, IA can
夽 Special thanks to the Pathology Department, Caen, France & National Center of Pathology, Vilnius, Lithuania. ∗ Corresponding author. Tel.: +33 2 33 01 46 29. E-mail address: [email protected] (P. Belhomme). http://dx.doi.org/10.1016/j.compmedimag.2014.11.006 0895-6111/© 2014 Elsevier Ltd. All rights reserved.
contribute to more accurate determination of prognostic histological factors by pathologists [1] but its action is limited to provide quantitative information as no multi-purpose method exists for segmenting or classifying images with a wide range of color or shape configurations, with which pathologists are familiar. The advent of digital scanners leads to generate whole slide images (WSI) from histological sections acquired at a full resolution. A main challenge this technology presents is the size of data to be processed which are generally very large and requires a huge storage capacity and processing needs. Conversely, the main advantage is that the complete tissue structure is easily accessed. In parallel, the aggressiveness of a cancer could result in morphological and architectural changes that can be observed in the tissue structure, and so be characterized by the object distribution on the slide, by the cross relations between objects and by the texture. This kind of information could contribute to evaluate a well-known concept: heterogeneity. Frequently addressed in signal processing, especially in terms of “entropy”, but more rarely in the field of imaging [2], the objective here is to propose a framework for measuring tissue heterogeneity in WSI and apply this methodology on breast cancer histopathology digital images. The key idea in this work is to not rely on segmentation (e.g. [1,2]) of individual structures to characterize heterogeneity, but to make use of classification of squared sub-images later called ‘patches’. In some previous works dedicated to the development of a computer-aided diagnosis system (CADS) based on image retrieval and classification [3,5], we have used a
52
P. Belhomme et al. / Computerized Medical Imaging and Graphics 42 (2015) 51–55
method coming from spectral graph theory, the Diffusion Maps (DM) [6], to process WSI split in small squares called ‘patches.’ The DM algorithm, in which eigenvalues and eigenvectors of a Markov matrix defining a random walk on the data are computed, allows to both cluster non-linear input data thanks to its inner classification properties preserving local neighborhood relationships, but also to reduce the input data dimensionality in a space (usually a 2D or 3D space) where it is therefore possible to compute euclidean distances between the objects to be analyzed [6,8]. To briefly describe the complete CADS we are developing, a first step consists in building a knowledge database involving many features extracted from a set of well-known images; this is an ‘off-line’ procedure conducted once. These features are represented by vectors of non-linear data acting as a signature. In a second step, signatures are obtained from new unknown images and then compared with those already present in the database; this is an ‘on-line’ procedure that has to be conducted each time a new image is processed. In a last step, before rendering a qualitative measure of similarity/dissimilarity between the supervised images and new unknown images, a feedback procedure would have to be processed in order to eliminate most of the artifacts that usually corrupt initial knowledge databases. But this general approach, especially with the DM algorithm, can also be derived to analyze a set of image patches coming from any WSI, with just the goal to compare their feature vectors. In this paper, we focus on a way to characterize the tissue heterogeneity at a regional level by working on the projected coordinates of image patches in a reduced 3D space. 2. Materials and methods Images used to illustrate this study come either from histological sections or tissue microarrays (TMAs) of breast cancer stained according to the Ki67 protocol and the Hematoxylin–Eosin-Saffron protocol (HES). They are acquired at a 20× magnification on a microscopic scanner device (ScanScope CS; Aperio Technologies). The resolution of WSI is 0.5 microns/pixel and the typical image size of any histological section can reach 50,000 × 40,000 pixels, corresponding to an uncompressed file of size 5.6 Gb. They are stored in the TIFF 6.0 file format (compression 30%). In each WSI some regions of interest (ROI), chosen by the pathologist, are then split in connected blocks of size 50 × 50 pixels, called “patches”, and each patch is associated with a numerical signature expressed as a feature vector. Tools developed here are written in Python language with the help of specialized modules (PIL: Python Imaging Library, SciPy-Numpy, Matplotlib, Mahotas. . .). 2.1. Features extraction The feature vectors embed some statistical measures obtained from color components and some texture parameters. In this study, 61 numerical values per component were calculated from up to two color spaces (RGB for Ki67 images and RGB + H&E color deconvolution for HES images [11]) plus 18 values obtained from the excess-red component (2R-G-B). From any given component, the 9 statistical features are the mean, median, mode, Skewness, Kurtosis parameters with also the 20%–40%–60%–80% quantiles of its cumulated histogram (initially reduced to 64 values). The 52 texture features correspond to the classic 13 Haralick parameters obtained in four directions [12]. We also added to the texture features 18 new values coming from the intrinsic statistical parameters of regionalized variables in three directions (0◦ –45◦ –90◦ ). They are derived from the nugget, sill and range of a geostatistical method [13]. Finally, each feature vector thus contains either F = 201 (61 × 3 + 18) or F = 323 (61 × 5 + 18) numerical values, here with the predominance of texture features in order to be quite independent from the
color staining variations encountered between different laboratories (and often inside the same laboratory). In order to later compare feature vectors, and considering the sparse numerical range of their values, the symmetric Kullback–Leibler distance [3,4] has been retained for its ability to easily manage such a case, while remaining fast to implement. The distance between two vectors p1 ,p2 of length F is given by: 1 DKL (p1 , p2 ) = p1j · log 2 F
j=1
p1j p2j
+ p2j
p2j p1j
(1)
2.2. Dimensionality reduction In any classical CADS, one of the key components is a visualization tool showing relationships between supervised images, stored in a knowledge database, and new images that are presented to the system. Typically, these relationships may be expressed as a connected graph in a 2D or 3D space where one hopes to find distinctive clusters corresponding to histological types or sub-types. It is therefore mandatory to reduce dimensionality from F (201 or 323 in our application) to 2 or 3 dimensions. With feature vectors containing non linear data as we are faced with, authors in [7,8] have shown that it was not appropriate to perform a principal component analysis (PCA). Instead, methods based on Spectral Connectivity Analysis (SCA) such as the Diffusion Maps, involving eigenvalues and eigenvectors of a normalized graph Laplacian, are well suited for this task. Let X = {x1 ,x2 ,. . .,xn } be a set of n patches that we assimilate to a fully connected graph G, that means a distance function is computed for each pair {xi ,xj }. A n × n kernel P is obtained from a Gaussian function whose coefficients are given by: p(xi , xj ) =
w(xi , xj )
(2)
d(xi )
with d(xi ) =
w(xi , xk )
(3)
xk ∈ X
and
D
w(xi , xj ) = e−
KL (xi ,xj ) ∈
(4)
In fact, p(xi ,xj ) may be considered as the transition kernel of the Markov chain on G. In other words, p(xi ,xj ) defines the transition probability for going from xi to xj in one time step. The eigenvectors k of P, ordered by decreasing positive eigenvalues, give the practical observation space axes. It must be noticed that 0 is never used since linked to the eigenvalue = 1 (i.e. the data set mean or trivial solution). A 3D projection can be then obtained along axes (1 ,2 ,3 ). Choosing ε in w(xi ,xj ) is an empirical task which should permit a moderate decrease of the exponential in Eq. (4); some works [8] use the median value of all DKL (xi ,xj ) distances whereas other works [6] use the mean distance obtained in the k nearest neighbors from a subset of X, which finally yields to quite the median value in many cases. We have retained the first solution for its simplicity and to lessen the overall computational time. Fig. 1a–c show the resulting projections obtained from the DM algorithm for 300 patches of fibroadenoma with ε = 0.5 × medianValue, ε = medianValue and ε = 2 × medianValue. One can see that the general shapes of point clouds are not the same but, due to the inner classification properties of DM, we have shown that the local neighborhood relationships were preserved in terms of histological types/sub-types [9].
P. Belhomme et al. / Computerized Medical Imaging and Graphics 42 (2015) 51–55
53
2.3. Pseudo-color space At the end of the DM algorithm, the three main eigenvectors of matrix P provide an uncorrelated space in which all patches are represented by a point (this space is thereafter called “Diffusion Maps” cube or DM cube). The Diffusion Maps method, belonging to the field of unsupervised learning algorithms, has the interesting property to express the complex structure of a multidimensional data set (also called manifold) in a reduced space where the Euclidean distance can be used, while preserving initial local neighborhood relationships. Thus, patches that seem to be similar in the WSI will be projected in the same area within the DM cube. In order to ease displaying the projected patches, the DM cube is considered as a color space in which each patch is assigned to a pseudo-color corresponding to its physical location. To mimic human visual system, the first main eigenvector is taken for the green component, the next two eigenvectors respectively for the red and blue components. Therefore, estimating the heterogeneity of histological sections comes down to measure pseudo-color differences of their patches in the DM cube. 2.4. Heterogeneity measurement In this paper, three scale levels for rendering a heterogeneity measurement are explored. The first level is local and computes a color gradient image in 4-connexity, which brings to light spatial variations of pattern structures in tissue sections. The result image is a classic gradient image with low values where patch signatures are slightly different between neighbors, and higher values on pattern transitions. This allows separating areas such as stroma and epithelium. The heterogeneity measurement can then be obtained by computing the average of all the gradient values, associated with their standard deviation. In Fig. 4a the mean value is expressed. The second level proposed is regional as it provides a linear measurement of heterogeneity along Bresenham lines plotted in the WSI. For a line LBj whose extremities correspond to patches A and B and that includes a series of intermediate patches Pi , the algorithm [2] adds the color differences of pairwise connected patches Pi and Pi+1 along LBj , the sum being finally divided by LBj length. If changes in patch signatures along a Bresenham line are regular and proportional to the length, then patches should rather be aligned in the DM cube and H(LBj ) should be low. To ensure a relevant measure of heterogeneity, multiple Bresenham lines with variable lengths and directions are randomly drawn in the WSI. In Fig. 4b the mean value of all H(LBj ) is expressed. Finally, the third level proposed to compute the heterogeneity is global. It involves computing the Rao’s quadratic entropy (RQE) [10] from N patches according to the following formula: RQE =
N N−1
dij pi pj
(5)
i=1 j=i+1
where dij is the pseudo-color difference between two patches and pi ,pj are the probabilities for patches to be present, set up here such as pi = pj = 1/N. RQE usually quantifies the general disorder of a population as the average difference between individuals; in our case the population is represented by image patches. Fig. 1. (a) Projection of 300 patches in a 3D space with ε = 0.5 × medianValue. Colors have no sense here except to represent groups of 100 patches. (b) Projection of 300 patches in a 3D space with ε = medianValue. (c) Projection of 300 patches in a 3D space with ε = 2 × medianValue.
3. Results and discussion To simplify the displaying of our tissue heterogeneity quantification method in this abstract, tissue microarrays (TMAs) and sub-images extracted from WSIs of breast cancer are used. Fig. 2 shows four TMAs with their corresponding pseudo-color map in the DM cube. The homogeneous parts are visible and are separated
54
P. Belhomme et al. / Computerized Medical Imaging and Graphics 42 (2015) 51–55
Fig. 2. TMAs from virtual slides of breast cancer and their corresponding pseudo-color map.
Fig. 3. Four sub-images qualified as slightly heterogeneous (1, 2), greatly (3) and moderately (4).
by strong transitions in the tissue structure. On Fig. 3, four subimages of size 3000 × 2000 pixels were selected and qualitatively evaluated by a pathologist. Fig. 3(1) and Fig. 3(2), that have a similar aspect, are qualified as slightly heterogeneous, while Fig. 3(3) is said greatly and Fig. 3(4) moderately heterogeneous. Once the
100
10
1
0,8
25
Quadratic Rao's entropy
Regional Brensenham lines
Mean Color Gradient
1000
heterogeneity measure has been computed according to the three methods above mentioned (locally with the color gradient, regionally with Bresenham lines and globally with the Rao’s quadratic entropy), the returned relative values correlate with the qualitative pathologist evaluation (Fig. 4). In this study, we have exposed
20
15
10
5
0 Image 1
Image 2
Image 3
Image 4
0,7 0,6 0,5 0,4 0,3 0,2 0,1 0
Image 1
Image 2
Image 3
Image 4
Image 1
Image 2
Image 3
Image 4
Fig. 4. (a–c) Tissue heterogeneity measurement according to three scale levels. (a) Local with the mean value of a 4-connexity color gradient; (b) regional with the mean H(LBj ) values of randomly selected Bresenham lines; (c) global with the Rao’s quadratic entropy measure computed from the Diffusion Map cube.
P. Belhomme et al. / Computerized Medical Imaging and Graphics 42 (2015) 51–55
the ability to assess tissue heterogeneity in a WSI using a set of three different measurements. An interesting contribution of our approach is that we do not base the image processing task on object segmentation, which is a complex and very time-consuming procedure with such large images, but rather work on feature vectors that are extracted from patches. The tuning parameters then only correspond to the patch size and the features used in signatures which must be adapted to the input image. In our experience, texture features might be preferred to color features, especially for comparing images coming from different laboratories. Indeed, a wide range of color difference can be encountered, even with histological slides that have been automatically prepared according to the same staining protocol. That is why the number of texture parameters retained is higher than the number of statistical parameters (52 vs 9 per color component, with moreover 18 parameters per patch coming from the geostatistical method). 4. Conclusion Patch analysis, while it may appear as a basic approach, allows segmentation of areas with similar configurations, such as stroma or tumor epithelium, and thus to encode the spatial organization of the same structure within the entire image. In future works, with the goal to decrease the computational burden generated by the “connected patches” approach, we plan to work on sampled patches obtained from stereological regular grids.
55
References [1] Gurcan M, Boucheron L, Can A, Madabhushi A, Rajpoot N, Yener B. Histopathological image analysis: a review. IEEE Rev Biomed Eng 2009;2: 147–71. [2] Brooks FJ, Grigsby PW. Quantification of heterogeneity observed in medical images. BMC Med Imaging 2013;13:7. [3] Oger M, Belhomme P, Klossa J, Michels JJ, Elmoataz A. Automated region of interest retrieval and classification using spectral analysis. Diagn Pathol 2008;3(Suppl. 1):S17. [4] Kullback S, Leibler RA. On information and sufficiency. Ann Math Stat 1951;22(1):79–86. [5] Belhomme P, Oger M, Michels JJ, Plancoulaine B, Herlin P. Towards a computer aided diagnosis system dedicated to virtual microscopy based on stereology sampling and diffusion maps. Diagn Pathol 2011;6(Suppl. 1):S3. [6] Coifman RR, Lafon S. Diffusion maps. Appl Comput Harmonic Anal 2006;21(1):5–30. [7] Coifman RR, Lafon S, Lee AB, Maggioni M, Nadler B, Warner F, et al. Geometric diffusions as a tool for harmonics analysis and structure definition of data: diffusion maps. Proc Natl Acad Sci 2005;102(21):7426–31. [8] Belkin M, Niyogi P. Laplacian eigen maps for dimensionality reduction and data representation. Neural Comput 2003;15:1373–96. [9] Belhomme P, Oger M, Michels JJ, Plancoulaine B. Out-of-sample extension of diffusion maps in a computer aided diagnosis system. Application to breast cancer virtual slide images. Diagn Pathol 2013;8(Suppl. 1): S9. [10] Rao CR. Diversity and dissimilarity coefficients: a unified approach. Theor Pop Biol 1982;21:24–43. [11] Ruifrok AC, Johnston DA. Quantification of histochemical staining by color deconvolution. Anal Quant Cytol Histol 2001;23:291–9. [12] Haralick RM, Shanmugam K, Dinstein I. Textural features for image classification. IEEE Trans Syst Man Cybern 1973;SMC-3:610–21. [13] Cressie NAC. Statistics for Spatial Data. New York: John Wiley & Sons; 1993.
