Image Anal Stereol 2015;34:1-25 doi:10.5566/ias.1109 Original Research Paper MULTIVARIATE MATHEMATICAL MORPHOLOGY FOR DCE-MRI IMAGE ANALYSIS IN ANGIOGENESIS STUDIES GUILLAUME NOYELC,1, JESUS ANGULO1, DOMINIQUE JEULIN1, DANIEL BALVAY2 AND CHARLES ANDRÉ CUENOD2 1MINES ParisTech, Centre de Morphologie Mathématique, Mathématiques et Systemes – 77305 Fontainebleau, France; 2LRI-PARCC U970 Paris – Descartes University – APHP, HEGP, Service de Radiologie – Paris, France e-mail: guillaume.noyel@mines-paris.org, jesus.angulo@mines-paristech.fr, dominique.jeulin@mines-paristech.fr, daniel.balvay@inserm.fr, charles-andre.cuenod@egp.aphp.fr (Received December 30, 2013; revised April 26, 2014; accepted May 30, 2014) ABSTRACT We propose a new computer aided detection framework for tumours acquired on DCE-MRI (Dynamic Contrast Enhanced Magnetic Resonance Imaging) series on small animals. In this approach we consider DCE-MRI series as multivariate images. A full multivariate segmentation method based on dimensionality reduction, noise .ltering, supervised classi.cation and stochastic watershed is explained and tested on several data sets. The two main key-points introduced in this paper are noise reduction preserving contours and spatio temporal segmentation by stochastic watershed. Noise reduction is performed in a special way that selects factorial axes of Factor Correspondence Analysis in order to preserves contours. Then a spatio-temporal approach based on stochastic watershed is used to segment tumours. The results obtained are in accordance with the diagnosis of the medical doctors. Keywords: classi.cation, DCE-MRI series, multivariate mathematical morphology, segmentation, stochastic watershed, tumours. INTRODUCTION DCE-MRI (Dynamic Contrast Enhanced MRI) time series is a medical imaging modality useful to characterise the process of tissue vascularisation. As tumours correspond to zones of angiogenesis, where the vascularisation is increased, DCE-MRI series are a convenient way of identifying or characterising potential tumours. Hence, DCE-MRI provides an additive and functional information to more current morphological images. Due to the increasing amounts of images, the creation of tools to assist medical doctors is of great interest for the analysis of these images and the detection of candidate tumour regions. In this paper, our goal is to supply an automatic tool to help users to localise tumours, based on this vascular information. An automatic process is required to limit operator variability during tumour delineation. In current morphological images such question is very challenging because, currently, gray level differences between tumour and adjacent tissues are not suf.cient for a con.dent and stable separation between tissues. Conversely DCE imaging provides potential richer information to determine this difference. However this information is distributed among all the sequence of images, requiring mathematical development. Here we show the results on DCE­MRI images of recent developments in mathematical morphology segmentation for multivariate images. Our developments are particularly useful for a visual evaluation of the tumour extension, but also as a pre-processing step before pharmacokinetic modeling, evaluated then on images with less noise. The corresponding evaluation of functional parameters such as blood .ow, blood volume, permeably-surface of capillaries, etc. (Sourbron and Buckley, 2012) should be improved accordingly. The considered images are DCE-MRI series of L = 512 channels of size N × N, N = 128, pixels acquired at a regular step of t = 1 second, in time, on mice presenting tumours (Balvay et al., 2005). In DCE­ MRI imagery, tumours are regions corresponding to the accumulation of the contrast product. This accumulation is characterised by an increasing kinetics of the temporal signal for each pixel of the tumour. Our aim is to show the potential of a new method for segmentation purpose in medical imagery. The tests, made on 25 different series, were used to develop the new method based on hyperspectral mathematical morphology. The objective is to achieve computed aided diagnosis of tumours. From an image processing point of view, DCE­MRI images are time series which ful.l a fundamental temporal coherence hypothesis: at each time step, the MRI image is acquired on the same object of interest and the images are registered, i.e., any pixel has the same spatial position for every time step (i.e., for every image channel). In our case, in each experiment the channels are registered. Due to these assumptions, the sequences of DCE­MRI images may be interpreted as multivariate, i.e., hyperspectral images, of the time evolution of the scene observed. Some images of the sequence, interpreted as channels of a series, are shown in Fig. 1. The tumours (PC3 – human prostatic) are characterised by an hyper-vascularised ring and an hypoxic centre or even a necrotic centre which is due to the distance of the afferent vessels from the centre. This vessels are repulsed from the centre by tumoural proliferation. More details about acquisition conditions are given on the appendix section. Hyperspectral images are multivariate discrete functions with several tens or even hundreds of spectral bands. In a formal way, for each pixel of a 2D (or a 3D), hyperspectral image is viewed as a vector, with values associated with a given wavelength, time or any index j. Each wavelength, time or index has a corresponding image, named channel, in two (or three) dimensions. In the text below, we use the term of spectrum and spectral channel to describe temporal phenomena. The segmentation of hyperspectral images by hyperspectral mathematical morphology has been used for remote sensing images (Benediktsson et al., 2005; Noyel et al., 2011). However, we show here the generality of our approach since the methods are also useful for DCE-MRI series analysis. STATE-OF-THE-ART ON DCE-MRI SERIES ANALYSIS Dynamic Contrast Enhanced Imaging (DCE-Imaging) is an increasingly used non-invasive imaging strategy to analyse tissue micro-vascularisation and perfusion. This technique is based on 1) an injection of a bolus of contrast agent, 2) an imaging modality (CT, MRI or Ultrasound imaging) in a sequential mode (Ivancevic et al., 2001), and 3) a method to analyse the kinetics of tissue enhancement over time (Brix et al., 2012; Sourbron and Buckley, 2012). DCE-Imaging has been widely demonstrated to be useful in detection and characterisation of lesions such as ischemia and tumours, in a variety of tissues such as brain, heart, breast and liver (Van Dijke et al., 1996), as well as in prediction and evaluation of the effects of therapies (Zahra et al., 2007; O’Connor et al., 2008; Leach et al., 2012). DCE-Imaging is especially useful in analysing tumour angiogenesis, i.e., the fast growth of a new and chaotic capillary network in tumours induced by growing factors secreted by the fast multiplying tumour cells (Brasch et al., 2000). This functional imaging .eld is therefore growing rapidly an constitutes a valuable addition to traditional morphological imaging, signal intensity, shape and size of the lesions. However, several meta-analyses (Zahra et al., 2007; O’Connor et al., 2008) have underlined the dif.culties in comparing DCE imaging results between different centres due to differences in data acquisition and analysis of this technique, as well as to signal to noise limits. Signal to noise limits are due to fast dynamic acquisitions and due to motion artefacts induced by breathing, heart beating, bowel movements as well as involuntary patient movements. To improve reproducibility, the challenge consists in minimising the in.uence of the local experimental conditions on data. 2 f.1 f.12 f.13 0s 11s 12s f.256 f.512 255s 511s Fig. 1. Five channels of hyperspectral image of a mouse “serim447” which is a temporal series (128 × 128 × 512 pixels) with 512 channels acquired every second. In order to get a better understanding of the images, several portions are labeled in Fig. 2 : i) a portion of the tumour is in green, ii) a portion of the heart cavities is in blue, iii) a portion of the background is in red and iv) a portion of the lungs is in white. Fig. 2. Labeled portions of DCE-MRI series. Image Anal Stereol 2015;34:1-25 In Ding et al. (2009), a method based on Karhunen Loeve Transform has been presented in order to reduce noise on cardiac cine MRI. The factorial axes are selected according to the autocorrelation function of each eigenimage. The axes are retained in a contiguous way by an automatic criterion based on half-maximum height of the autocorrelation peak. In Ding et al. (2010), an extension of this method for spatially variant noise has been presented. It is based on the eigenvalue distribution of random matrices. In Balvay et al. (2011), an improvement of Signal to Noise Ratio for Dynamic Contrast-Enhanced Computed Tomography and Magnetic Resonance Imaging with PCA is explained. A new criterion, the fraction of residual information, is proposed to automatically select the factor axes. It takes into account the temporal order of the images in the series. In this paper we present the most interesting results on DCE-MRI series of a method that reduces the noise by Factor Correspondence Analysis, followed by a classi.cation and a segmentation stage. Readers interested in more details and proofs on our method are invited to consult (Noyel, 2008c). A short overview of some of the results of this study has been presented in (Noyel et al., 2008b). PAPER ORGANISATION The .rst part of this paper introduces some pre-requisites on hyper-spectral image processing. The second part is focused on .ltering and data reduction of hyperspectral images. In the third part the classi.cation step is explained and in the fourth part the segmentation by standard and stochastic watershed of DCE-MRI series is presented. The full image analysis process is illustrated and validated by an application to the automatic detection of tumours on animals. PRE-REQUISITES NOTATIONS In order to analyse DCE-MRI series with methods based on multivariate image processing, we use a speci.c notation for hyperspectral images. By using this generic notation, we consider the temporal dimension of the image as the “spectral dimension”. Let : E › T L f. : x › f. (x)=f.1 (x), f.2 (x),..., f.L (x) (1) be an hyperspectral image, where: – E . R2, T . Rand T L = T × T × ... × T – x = xi \ i .{1,2,...,P} is the spatial coordinates of a vector pixel f. (xi) (P is the number of pixels in E) – f. j \ j .{1,2,...,L} is a channel (L is the number of channels) – f. j (xi) is the value of vector pixel f. (xi) on channel f.j . Several approaches exist to analyse DCE-MRI: 1. spatial analysis: channel f.j by channel f.ij 2. spectral analysis: vector f. (xi) by vector f. (xii ) 3. spatio-spectral analysis: simultaneous use of both approaches f. . The methods are illustrated in Fig. 3. Generally, the spatio-spectral approach gives the best results. (a) Spatial (b) Spectral (c) Spatio-spectral Fig. 3. Three ways to analyse DCE-MRI sequences: (a) spatial, (b) spectral and (c) spatio-spectral approach. NECESSITY OF DATA REDUCTION DCE-MRI images are time series composed of several hundreds of time channels. As channels are not statistically independent, the reduction of spectral dimension is necessary: 1. to reduce “Hughes phenomenon” (Hughes, 1968) also called the “curse of dimension” 2. to reduce the amount of data and therefore to reduce the computational time. Hughes phenomenon has been studied in the case of hyperspectral images, among others in Landgrebe (2002) and Lennon (2002). To tackle this problem several data reduction methods exist, e.g., Correspondence Analysis, Principal Component Analysis, Independent Component Analysis, etc. Modeling the spectrum is also useful for reducing the spectral dimension as the example shown in this paper. METHODOLOGY TO SEGMENT DCE-MRI TIME SERIES The overall proposed methodology to segment DCE-MRI time series is composed of three steps (Fig. 4): 1. a .ltering of the images and a reduction of their spectral dimension 2. a classi.cation of the vector pixels to a given number of classes 3. a segmentation step based on a function to .ood by a watershed transform. This function combines the spatial and spectral information. It consists of a probability density function of contours The .rst two steps are included in the pre­processing step, which in this paper is based exclusively on the spectral information. respect to the spectral dimension of the original image space. This reduced space is composed of factor pixels cf . of the hyperspectral image f. . Factor Correspondence Analysis is useful to reduce spectrum dimension (Noyel et al., 2007c). Similar results can be obtained with other methods such as Principal Component Analysis (PCA) or Independent Component Analysis (ICA), etc. We have given priority to FCA because it is ef.cient in segmenting multivariate images with positive pixels values. The original image f. can be reconstructed from a limited number of factors leading to a good approximation ff. of this original image f. . The reconstructed image contains, under certain conditions, a spectral noise that is smaller than the noise on the original image. Therefore, a method which reduces the noise on the factor pixels cf will . be explained in this section. Additionally, the reduction of the spectral dimension of the image f. by .tting a spectrum model on the .ltered image f . will be presented. This is to take advantage of the prior knowledge of the spectrum. By modeling the spectra, some maps of the parameters of the model are obtained. These maps constitute a reduced space which is useful for further classi.cation and segmentation. DENOISING AND DIMENSIONALITY REDUCTION BY DATA ANALYSIS Introduction to FCA Data analysis is a transformation . of the space of the original image f. , with a dimension L, into a space of another hyperspectral image cf . , of reduced dimension K < L, and a set of parameters: Fig. 4. Framework to segment DCE-MRI time series as . . . . .. T L › T K such that K < L hyperspectral images. . . = cf . (x)= ff ....... (x) .K › (2) with: FILTERING AND DATA REDUCTION – c. f the factor pixels of the hyperspectral image. It is the coordinates of the vector pixels on the factorialMultivariate image denoising and dimensionality axes. reduction is addressed in this study with Factor Correspondence Analysis (FCA; Benzécri, 1973). A – df the factors of the channels. They are the .. factorial space of reduced dimension is obtained with coordinates of the channels on the factorial axes. 4 (x),...,c c .1 ....... df = d.f 1. ,...,df .. .K . The spectral approach is actually a global analysis on the image, because we compare all the vector pixels . : {µ. }.{.i. f. (x) =1...K .. . . . between them. On the other hand spatial approach is more local, because a given vector pixel is mainly }i=1...P {.. j}j=1...L f = .i .j fij compared to its nearest neighbours. Image Anal Stereol 2015;34:1-25 µ. is the inertia of factorial axis .. .i. the marginal frequency of the vector pixel f. j (xi)f. (xi): .i. = .Lj=1 .Lj=1 .P (xi) . i=1 f. j .. j the marginal frequency of the channel f.j : f. j (xi) .. j = .Pi=1. .Lj=1 .Pi=1 f. j (xi) f = .i .j fij = .i .jf. j (xi) the sum of all the values f.j (xi) of the image f. . –––– Fig. 5. Part of inertia of the thirty .rst factorial axes. For data analysis, a limited number K of factors is usually selected. Therefore, data analysis is a projection of the pixels of the original image f. into a space of smaller dimension K < L and often K « L. . -1 pseudo-inverse transform. It is an exact transform if all f of the image ff. The reconstruction is a c f .1 14.63% c f .2 5.73% c f .3 3.91% the axes are kept (K = L - 1). It consists in partially reconstructing the image f. from the pixels factors c f . and some other parameters. The reconstructed image ff. , with a limited number of factors, is an approximation of the original image: c f .4 3.33% c f .5 2.27% c f .6 2.12% . T K c . f › T L / K < L f. (x) . . . . . . (x) f d f .. ...... ...... . -1: › f with ff. (x)= ff.1 (x),..., ff.L (x) (3) {µ. }.=1...K {.i. . . . . }i=1...P {.. j}j=1...L f = .i .j fij c f .7 1.78% c f .8 1.71% c f .9 1.19% Fig. 6. The factor pixels on the 9 .rsts axes and their inertias. The factors pixels typeset upright are kept; = . those italicized are rejected due to a signal to noise ratio which is below a given threshold. Selection of the factor axes In Fig. 5 one can notice that the .rst 5 axes contains The number of factorial axes to be kept needs to be the main part of the total variance or inertia (about f f chosen. It depends of: 30%). However, by observing (Fig. 6) the images of the factor axes c we notice that the factor pixels c .k .7 f .8 and c contain mainly noise while the factor pixel the part of inertia (or variance) that they explained – f .9 contains mainly signal. In order to quantify the in the data cloud. Several tests, based on inertia, c amount of information contained in the factor pixels, exist that allow one to choose the correct number of their signal to noise ratio is estimated by the method axes such as the “Kaiser criterion” (Kaiser, 1960) proposed in Noyel et al. (2008a). or the “scree test” (Cattell, 1966). The channel c f .k (or f. j ) is considered as a realisation of a random function. For each factorial – the amount of information contained in the factor pixels. We have introduced a new criterion based f axis, the centred spatial covariance is estimated by on the signal to noise ratio of the factor pixels. assuming that c (x) is a stationary function: .k g.k (h)= E[c f .k (x)c f .k (x + h)] , (4) with c f ff .k the centred channel .k: c .k (x)= c .k (x) - E[c f .k (x)] and E[Y ] the expectation of the random f .k (x)] corresponds to the mean of the variable Y . E[c channel c f .k (x). It is known that the covariance of a noisy random function shows a discontinuity for h › 0. This discontinuity, equal to the variance of the noise, is called the “nugget” effect (Matheron, 1970; 1975). f This is illustrated in Fig. 7, where the covariance is Fig. 8. Signal to noise ratio of the factor pixels c. of plotted for a channel without much noise and for a the image and a threshold for a SNR of 0.3. noisy channel. In the latter case, a peak at the origin of the covariance can be noticed. The variance of the signal can be estimated by the difference between Table 1. The SNR of the factor pixels (greater than 0.3 g.k (0) and the nugget effect. It can be computed by an automatic extraction of the peak at the origin of the covariance image after a morphological opening .. The structuring element is chosen as small as possible. It is a square of size 3 × 3 pixels. From this extraction, a signal to noise ratio is estimated, according to the following de.nition: Var(signal) .g.k (0) SNR.k == , (5) Var(noise) g.k (0) - .g.k (0) with g the centred covariance and . the morphological opening. By observing the signal to noise ratio (SNR) of the different factors in Fig. 7, it appears that the channels for those kept) and their inertia part. The 16 factorial axes which are kept are typeset upright and the 3 rejected are italicized. Axes SNR Inertia (%) Axes SNR Inertia (%) 1 3.91 14.63 6 0.58 2.12 2 1.02 5.73 7 0.15 1.78 3 0.99 3.91 8 0.08 1.71 4 1.07 3.33 9 0.52 1.19 5 0.79 2.27 10 0.51 1.15 Axes SNR Inertia (%) Axes SNR Inertia (%) 11 0.52 1.08 16 0.38 0.75 12 0.51 1.01 17 0.32 0.67 13 0.39 0.92 18 0.35 0.57 14 0.29 0.82 19 0.33 0.54 15 0.33 0.78 c f .7 and c f .8 are noisier than others. In what follows, we propose channels with a SNR greater than 0.3 are The pixels factors typeset upright in Table 1 and in Fig. 6 are kept for data analysis, while those italicized are rejected. We notice that the selected factorial retained for reconstruction. axes, are not contiguous in terms of inertia, as only components with a SNR > 0.3 are kept. Reconstruction In Fig. 9, some channels of the image f. are displayed before and after reconstruction. On channel , we notice that one ventricle of the heart appears .g.1 f.12 c f .1 g.1 dark while the other appears bright. This opposition is preserved after the reconstruction of the channel using the selection of axes by SNR. One can f.12 also notice, that the image ff. is a good reconstruction of the image f. and that a part of the noise, in the f .100 g.100 .g.100 original image f. , is removed by FCA reconstruction. Fig. 7. Covariance before (g) and after a Therefore, the denoising has been made by a spectral f morphological opening (.g) channels c the factor pixels .ltering of FCA, which preserves the spatial structures on f .100 (without much noise) and c (noisy of the image. This is a crucial improvement before .1 Image Anal Stereol 2015;34:1-25 f.1 f.12 f.13 ff.1 ff.12 ff.13 f.256 f.512 ff.256 ff.512 Fig. 9. Five channels of the hyperspectral image (i.e., the DCE-MRI series) before ( f.i) and after reconstruction ( ff.i) from 16 axes of a FCA. Noise reduction Like Principal Component Analysis, FCA is useful to perform a spectral .ltering of the data. In Benzécri (1964; 1973) and in Orfeuil (1973), FCA is used to .lter arrays of “Euclidean data marred with errors”. The use of PCA to .lter hyperspectral data was primarily illustrated in Berman (1985) and in Green et al. (1988). A study about noise reduction by FCA on hyperspectral images is available in Noyel (2008c). In the following text, only the interesting results obtained in this previous study are used. In Noyel (2008c), we have shown two sequences of FCA and reconstruction steps reduce the noise more than just a single sequence. with t. a spectral translation of the data. Actually we add the constant . to the value of the image f. in order to recover positive values, as required for FCA: . T L › T L . t. : . f. › t. (f. ) such as . .i = 1...P , . j = 1...L . t. (f.j )(xi)= f.j (xi)+ . , . . R. (8) Therefore the second sequence of FCA-reconstruction is not identical to the .rst one because the data have been modi.ed by the translation t. . For . we propose to use the minimum of all the f(1) data in the .rst reconstructed image f. : ff(1) . = min (xi) , i = 1 ...P , j = 1 ...L . (9) .j i, j More details about the two sequences of FCA-reconstruction are given in the appendix. Notice that the inertia is evaluated on the axes of the original, i.e., .rst, FCA. In order to compare the original DCE-MRI series f. with the series after 2 sequences of FCA-reconstruction with 16 axes ff. = ff(2) , the residues are . computed channel by channel: .i . [1...P] , . j . [1...L] : r. j (xi)= | f.j (xi) - ff.j (xi)| . (10) An hyperspectral image of residues r. may be de.ned: r. = |f. -ff. | , (11) of which the channels r. j are equal to: . j . [1...L] : r.j = | f.j - ff.j | . (12) The centered spatial covariance is also gr. j measured on the residues of the channels: . j . [1...L] : = E[r.jr.j ] , (13) gr. j 1 with r.j = P .Pi=1 r.j (xi). f(1) A sequence of one FCA and a reconstructionf. is equal to the reconstructed image. It is written: ff(1)= .f-1 . . (f. ) . (6) . Two sequences of FCA-reconstruction are de.ned as: ff(2)= t. . .f-1 . . .t-. . .f-1 . .(f. ) , (7) . In Fig. 10, some channels after 2 FCA reconstructions and their residues are presented. Notice that some noise has been removed. This noise is due to the acquisition process and is mainly located at the centre of the image. The centered covariances contains a peak at their origin, which means that gr. j the residues correspond to noise. f(2) f.i, and after 2 FCA-reconstructions with 16 axes f, . their residues r. and the covariances on the residues gr. . The histogram of the residues has been normalised for visualisation. Some spectra of the reconstructed image after 2 FCA are plotted in Fig. 11. The .ltered spectra by 2 sequences of FCA-reconstructionff. (xi) have a smaller variability than the original spectra. Moreover after the .ltering stage the general trend of the spectrum is enhanced. For example, the signal of the spectra corresponding to the tumour signal is increasing. This corresponds to the accumulation of the product of contrast inside the tumour. In conclusion we propose a noise reduction method which consists in applying two series of FCA-reconstructions. This approach reduces the temporal (spectral) noise while preserving the contours. This is an important point for a further segmentation. STRONG DIMENSIONALITY REDUCTION BY SPECTRUM MODELING Following the noise reduction which preserves the spatial information by data analysis, an additional dimensionality reduction by spectrum modeling is proposed. In this approach, a parametric model is .tted to each spectrum and consequently, the images of parameters can be seen as maps. These maps are useful for classi.cation and morphological segmentation. We notice that the set of these maps of parameters (p1(x),..., pM(x)) constitutes a multivariate image with a reduced dimension: p(x)=(p1(x),..., pM(x)) . (14) The .tted model may take into account the physical phenomena, such as the physiological pharmacokinetic model used in Brochot et al. (2006). Their model is based on differential equations on six compartments: arterial and venous plasma, tumour (split into capillaries and interstitium), and the rest of the body (also split into capillaries and interstitium). However, in the current study we adopt a simpler linear model. More precisely, for each spectrum (time series) of the .ltered image ff. , a line model is .tted after removing the .rst 20 values which correspond to a transitory phenomenon (Fig. 12): Image Anal Stereol 2015;34:1-25 ff. (xi) ~ a(xi). + b(xi) , .i = 1...P ,with . = .j1 ....512 , (15) with j1 the .rst value after the peak ( j1 = 21 for our images). This model is made with two parameters: the slope p1 = a and the intercept p2 = b. In order to take into account some information contained in the transitory part of the time series, a third parameter is used. It corresponds to the amplitude of the signal over the twenty .rst values of the spectrum. This transitory part is characteristic of the injection of the contrast agent used in this imaging modality. This parameter is called the rise p3 = m and it is de.ned by: m(xi)= max ( ff.j (xi)) - min ( ff.j (xi)) , j.[1; j1-1] j.[1; j1-1] .i = 1 ...P . (16) Hence, by .tting the model for each spectrum f. (xi), three maps of parameters are obtained (Fig. 13). We remark that the dimensionality reduction is very important because the original multivariate image with 512 channels is transformed into another multivariate image with only 3 channels. Even with this reduced amount of information, the main morphological structures of the mouse (the heart cavities, the tumour and the lungs) are clearly apparent. To conclude, spectrum modeling is very ef.cient in reducing the number of channels for tumour segmentation. It could be interesting in further studies to .t a a model taking into account physical properties of the injection of the contrast agent. slope a intercept b rise m Fig. 13. Parameters maps of the linear model .t for each pixel. TEMPORAL CLASSIFICATION OF DCE-MRI TIME SERIES After performing a denoising step and a dimensionality reduction, the classi.cation of pixels is made in the temporal dimension. Supervised and unsupervised methods are considered. UNSUPERVISED APPROACH: REGIONAL IMPROVED K -MEANS For the unsupervised approach, a k-means classi.cation and a model classi.cation are compared. A k-means classi.cation of pixels is based on the Euclidean metric (Diday, 1971; Hartigan and Wong, 1979). This classi.er must be used in a space in which the associated metric is Euclidean. It is the case for the factor space of FCA or of PCA. However, using k-means in the image space would overweight channels with large dynamics. So classifying is done in the factor space instead. For DCE-MRI series, the k-means classi.cation is performed with 5 classes in the factor image space cf . of the second sequence of FCA-reconstruction. In Fig. 14 some anatomical parts of the mouse have been classi.ed into 5 classes. The number of classes has been empirically de.ned. The classes are as follows: (1) the green class corresponds largely to the tumour -see top right corner of the image, (2) the red class corresponds mainly to the background, (3) the blue class corresponds to the heart cavities, (4) the black class corresponds to the lung and (5) the cyan class is an intermediate class. reference k-means model Fig. 14. (a) Reference re f , (b) k-means classi.cation in 5 classes .kmeans,5, (c) model classi.cation .mod,5 cff . f. in 5 classes on the .ltered image ff. . The results have not been completely satisfactory, so we have introduced an alternative approach of classi.cation, which clusters each spectrum in comparison to some reference spectra obtained on the k-means classi.cation. For each class obtained by k-means, a mean .ltered kmeans spectrum is computed sp. On these mean spectra, ff. kmeanskmeans sp spf a line model is .tted .f . = {. f. }5 For kk=1. each spectrum ff. (xi) of the .ltered image by two sequences of FCA-reconstruction the line model is ff. .tted .ff. = {.}iP =1. Then, each point xi is classi.ed xi by minimisation of the L1 distance between the model kmeans sp of the mean spectra . ff. and the model for each k ff. pixel .. The class of the point xi is written C(xi) xi (Eq. 17 and Fig. 16): kmeans spf . ff. C(xi)= argmind1(. ,. k xi ) k kmeans L sp ff. ff. = argmin . |.k (.j) - .xi (.j)| . (17) kj=1 kmeans sp ff. model of the mean spectra .k and the model of ff. the spectra at each pixel .. xi . cf . As one can notice in Figs. 14 and 24, the classi.cation .mod,5 by the model approach seems f . to be more robust than the k-means classi.cation because a statistical noise .ltering has been made when .tting a model with 3 parameters. Then, the k-means classi.cation has been computed on 19 channels while the model classi.cation has been computed on 3 channels with few noise. Therefore, the model decreases the entropy of the image by introducing a prior information present in the shape of the spectra. Nevertheless, the k-means classi.cation is necessary as a .rst step, in order to estimate the mean kmeans sp ff. spectra . . k SUPERVISED APPROACH: LDA WITH HISTOGRAM NORMALISATION A semi-supervised classi.cation by Linear Discriminant Analysis (LDA) is also considered for the DCE-MRI series. In particular, LDA is based on a train set composed of four distinct parts of the anatomy of the mouse: – the tumour in green train(green)= t1 – the heart cavities in blue train(blue)= t2 – the background in red train(red)= t3 – the lungs in black or white train(black)= t4. Each class of the training set, T = train = (t1,t2,t3,t4), is made of 80 vector-pixels , f. (xi) with 512 components, selected by an operator. By measuring the mean spectra of the .ltered image ff. train on these classes sp, we notice that the kinetics of f . train classes are different (Fig. 18). Image Anal Stereol 2015;34:1-25 train sp(green) f . train train train sp(blue) sp(red) sp(black) ff. f . ff. Fig. 18. Classes of the training set, train, and mean spectra of these classes sptrain on the .ltered imageff. . ff. The LDA is performed into three different spaces: – the .ltered image space: {ff. (x)|x . T } – the PCA factor space of the spectra of the training set: .({ff. (x)|x . T }) – the three parameters space: {p(x)|x . T } In Fig. 19, the classi.cations to three different spaces are very similar. The LDA in the .ltered image space or the LDA in the PCA factor space of the training set are a bit better than the LDA in the parameters space. The training classi.cation error and the test errors are computed on 80 pixels vector of the training set by a 5-fold cross validation (Hastie et al., 2003). Both classi.cation errors are equal to zero. reference LDA on ff. LDA on .train(f . (x)) LDA on p Fig. 19. Semi-supervised classi.cation LDA in 4 classes in 3 different spaces: the .ltered image space, the PCA factor space of the spectra of the training set, the parameters space. CLASSIFICATION ON SIMILAR DCE-MRI SERIES As we want to classify several DCE-MRI series of large image databases, it is necessary to develop classi.cation methods which are robust in all series, without the needs for a training set for each image series. By testing our methods on another series called “serim460”, we can notice in Fig. 20 that the unsupervised classi.cations (k-means and model) are correct compared to the given reference. However, the LDA classi.cation on the new series is not correct whatever the image space. In this case the training set is from another series (“serim447”). In the series “serim447” and “serim460”, some mean spectra are measured into similar zones (tumour, heart cavities, background and lung). In Fig. 21 we notice that the range of the spectra are not the same for both images. This is the origin of the problem of the classi.cation for a supervised method such as LDA. .kmeans,5 .mod,5 re f cf . f . .LDA,4 .LDA,4 .(ff. )|Ts447 p|Ts447 cf . f . .ltered image .LDA,4 and into the parameter space .(ff. )|Ts447 .LDA,4 . The training set of the LDA is from the series p|Ts447 “serim447”. “serim447” “serim460” (a) (b) In order to use LDA on other image series, using the initial training pixels selected on the .rst image, the range of the grey levels of the images must be similar. Otherwise the projection of other series in the classi.cation space of the pixels would be incoherent. Consequently, we introduced a range normalisation method based on histogram anamorphosis. To get more robust results, the multivariate image of parameters p is used. For each parameter, the cumulative distribution function (cdf) of the values is estimated. The cdf is the primitive of the density function estimated by an histogram. It is composed of 255 classes de.ned on each map of parameters of the initial series. The cdf of each image is transformed by a numerical anamorphosis on the grey-tone values in order to be similar to the reference cdf of the initial series “serim447” (Fig. 22). .LDA,4 The LDA classi.cation on series p “serim460” after cdf normalisation (Fig. 23) gives cdf before normalisation cdf after normalisation Fig. 22. Cumulative distribution functions computed on the parameters a, b and m of each series before normalisation (top) and after normalisation (bottom). The references are the cdf of the parameters of the initial series “serim447”. To each colour is associated the cdf of a series. .LDA,4 re f p |Ts447 Fig. 23. LDA classi.cation in 4 classes .pLDA,4 on the p maps of parameters of the series “serim460” after much better results than the same classi.cation without cdf normalisation. The training has been made on the cdf normalisation (Fig. 20). series “serim447”. Image Anal Stereol 2015;34:1-25 “serim406” “serim415” “serim450” “serim457” “serim461” “serim1441” Fig. 24. Classi.cations kmeans .kmeans,5 , by model approach .mod,5 and LDA .LDA,4 on various series cffp . f. p|Ts447 “serimxxx”. For validation purposes, we compare the In the case of tumours starting to die in their classi.cations by k-means, model and LDA on other centre, which would make them potentially smaller series (Fig. 24). The classi.cations based on the model than their real size, our method identi.es viable approach are more robust than those obtained with k-tissues (i.e., functional tissues). The central zonesmeans. LDA classi.cations also give good results. The which are not included inside the tumours are zones heart cavities are correctly classi.ed and the tumours are characterised by extended classes in green. of severe ischemia or necrosis. Currently, more and more DCE maps are analysed in the following way: i) identi.cation of the “necrosis” ratio (volume of the necrosis divided by total volume); ii) characterisation of viable tissues of the tumour. Without both analysis, all the circulatory parameters are underestimated during the growth of the tumour which leaves in its centre more and more necrosis (or .brosis). SPATIO-TEMPORAL SEGMENTATION BY PROBABILISTIC WATERSHED After having .ltered the temporal noise and having reduced the temporal dimension, the classi.cation produced a partition of the image into non connected classes only based on temporal information. Our aim is now to segment the series with smoother contours and regular classes by combining the spatial and the spectral dimension in the segmentation process. In order to segment hyperspectral images, we introduced a general method based on deterministic watershed (WS) in Noyel et al. (2007c) and another one based on stochastic WS in Noyel et al. (2007a; 2011). Based on previous work we want to present the potential of application of our methods on DCE­MRI series in order to segment regions to be tumour candidates. PRINCIPLE OF THE SEGMENTATION OF MULTIVARIATE IMAGES BY WS The watershed transformation (WS) is one of the most powerful tools for segmenting images and was introduced in Beucher and Lantuéjoul (1979). According to the .ooding paradigm, the watershed lines associate a catchment basin to each minimum of the landscape to .ood (i.e., a scalar or greyscale image; Beucher and Meyer, 1992). Typically, the landscape to .ood is a gradient function which de.nes the transitions between the regions. Using the watershed on a scalar image without any preparation leads to a strong over-segmentation (due to a large number of minima). There are two alternatives in order to get rid of the over-segmentation. The .rst one consists in .rst determining markers for each region of interest. Then, using the homotopy modi.cation, only the local minima of the gradient function are imposed by the markers of the regions. The extraction of the markers, especially for generic images, is a dif.cult task. The second alternative involves hierarchical approaches either based on non-parametric merging of catchment basins (waterfall algorithm) or based on the selection of the most signi.cant minima. These minima are selected according to different criteria such as dynamics, area or volume extinction values. Extinction functions (Meyer, 2001) are used to remove the non selected minima. The general paradigm of WS-based segmentation of multivariate images (Fig. 25) requires two different inputs: (1) some markers for the regions of interest mrk and (2) a landscape to .ood g which describes the “likelihood” of the frontiers between the regions. The markers can be chosen interactively by a user, or automatically by means of a morphological criterion .N (Meyer, 2001), or with the classes of a previous spectral classi.cation. The landscape to .ood is a scalar function (i.e., a greyscale image). For the deterministic WS, it is usually a gradient (actually its norm), or a distance function. For the stochastic WS, the function to .ood is a probability density function (pdf) of the contours appearing in the image. The extracted markers are imposed as sources of the landscape to .ood and the WS is computed. The results denoted WS(g, mrk) or WS(g,.N ). SPECTRAL DISTANCES AND GRADIENT ON MULTIVARIATE IMAGES A gradient image, actually its norm, is usually chosen as a function to .ood. After normalisation, the norm of a gradient image is a scalar function with values in the reduced interval [0,1], i.e., .(x) : E › [0,1]. In order to de.ne a gradient, two approaches are considered: the standard symmetric morphological gradient on each marginal channel and a metric-based vectorial gradient on all channels (Noyel et al., 2007c). The morphological gradient is de.ned for scalar images f as the difference between a dilation and an erosion by a unit structuring element B, i.e., .( f. j (x)) = .B( f.j (x)) - .B( f.j (x)) = .[ f.j (y),y . B(x)] -.[ f. j (y),y . B(x)] . (18) Image Anal Stereol 2015;34:1-25 The morphological gradient can be generalised to multivariate functions (Hanbury and Serra, 2001) with the following metric-based gradient: .df. (x)= .[d(f. (x),f. (y)) / y . B(x),y = x] -.[d(f. (x),f. (y)) / y . B(x),y = x] . (19) Various metric distances d(f. (x), f. (y)) between two vector pixels, useful for multispectral images, are available for this gradient such as: – the Euclidean distance: INTRODUCTION TO STOCHASTIC WS In a classical watershed, small regions strongly depend on the position of the markers, or on the volume (i.e., the integral of the grey levels) of the catchment basins, associated with their minima. In order to improve segmentation results, stochastic watershed aims at enhancing the contours of signi.cant regions which are relatively independent of the position of the markers. Stochastic WS method is described in Angulo and Jeulin (2007) and Noyel et al. (2007a; 2011). Starting from a series of M realisations of N uniform or regionalized random germs (or markers) {mrki(x)}M mrk series of watershed segmentation {sg(x)}M are ii=1 L made on a landscape to .ood (for example a gradient). .( f.j (x) - f.j (y))2 i=1, dE (f. (x),f. (y)) = With these M segmentations, the probability density j=1 function of contours pd f (x) is estimated by the Parzen window method (Duda and Hart, 1973) with a gaussian kernel (typically with a 3 pixels standard deviation working on contours of one pixel width). Due to – the Chi-squared distance: d.2 (f. (xi),f. (xii )) = the smoothing effect of the method, the WS lines 2 with a very low probability, which correspond to non signi.cant boundaries, are removed. L Sf.j (xi) f.j (xii ) . - , f fxi. ..j fxii . To obtain closed contours, the pdf image is j=1 segmented by a watershed segmentation into R with f ..j = .iP =1 f. j (xi), fxi. = .Lj=1 f.j (xi) and regions. The stochastic WS needs two parameters: S = .Lj=1 .P (xi). i=1 f. j 1. M realisations of germs. The method is almost independent on M if it is large enough (between – the Mahalanobis distance: 20–50). dM(f. (x),f. (y)) = 2. N germs (or markers): if N is small, a segmentation in large regions is privileged; if N is too large, the mrk (f. (x) - f. (y))t .-1(f. (x) - f. (y)) , over-segmentation of sgleads to a very smooth i pdf, which looses its properties to select the R where . is the covariance matrix between variables (channels) of f. . If channels are uncorrelated, the covariance matrix is diagonal. The diagonal values are equal to the channels variance .2 \ j . .j {1,2,...,L}. Therefore, the Mahalanobis distance becomes the distance inverse of variances: regions. As shown in Angulo and Jeulin (2007), it is straightforward to use N > R. The originality of our approach of stochastic WS for DCE-MRI series is to condition the germs used to build the pdf by a previous classi.cation. 2 SEGMENTATION BY STOCHASTIC WS . L f.j (x) - f.j (y) d1/.2 (f. (x),f. (y)) = ..j j=1 An important point is to choose an appropriate distance depending on the space used for image representation: Chi-squared distance d.2 and distance of inverse variances d1/.2 are adapted to the image space and Euclidean distance to factorial space. More details on multivariate gradients are given in (Noyel et al., 2007c; 2011). Another example of a multivariate gradient is given in Scheunders (2002). Pre-processing of the temporal classi.cation As the classi.cation is based on temporal information we want to introduce it by conditioning the random markers used to generate the probability density function of contours. In order to do this, a pre-processing stage is necessary for two reasons: – to reduce the “segmentation noise” appearing as the smallest connected components of the classi.cation .. – to introduce the necessary degrees of freedom to perform each WS used to build the pdf. Therefore an anti-extensive transform is performed on each class of . by a morphological erosion with a structuring element of size 3 × 3 pixels. An alternative is to make an area opening (Soille, 1999). Then an extensive transformation such as a closing by reconstruction (Soille, 1999) is performed to .ll the holes inside the largest connected components. As the classes of the transformed classi.cation are not anymore a partition of the image, a “void” class is introduced. It corresponds to the class appearing in place of the transformed connected components. For the DCE-MRI sequences the complete transform . is processed on the LDA classi.cation on the parameters space .p LDA,4 (Fig. 26). .LDA,4 .(.LDA,4 p p ) Fig. 26. Classi.cation by LDA on the parameters space .LDA,4 p and pre-processing transform .. Extension of stochastic WS to multivariate images The extension of stochastic WS to multivariate images was introduced in Noyel et al. (2007a) and detailed in Noyel et al. (2011). Segmentation of DCE-MRI series by a standard WS on a distance-based gradient. Before introducing the way to extend stochastic WS to multivariate images, let us show that the segmentation by a standard WS on a distance based gradient has limited ef.ciency. For the gradient based segmentations, a distance adapted to the image space is used: – the Euclidean distance for the image of the factor pixels cf . – the distance of inverse variances for the image of parameters p – the Euclidean distance for the image of the PCA factor space obtained from the spectra of the training set ctrain . ß Comments: The distance-based gradient is computed after a morphological leveling on each channel of the considered image in order to get a smoother gradient. A morphological leveling is a morphological transformation that reduces the positive and negative “peaks" according to a reference while preserving the transitions of the objects. See Meyer (2004) for more details. The reference for the leveling is obtained by a gaussian .lter of size 11 × 11 pixels. A WS segmentation with a volume criterion is performed in Fig. 27. A marker-controlled WS is performed in Fig. 28. The markers are made by a morphological opening on each connected component of the classi.cation with an hexagonal structuring element of size 5. We notice that the segmentation is not perfect, especially for the tumour. However, the marker based segmentation in the PCA factor space of the spectra of the training set seems to be the best segmentation. f .E (ctrain .E (c. ) .1/.2 (p) ) ß R-vol (.E (cf R-vol (.1/.2 (p)) R-vol (.E (ctrain sg. )) sg sg)) ß (a) (b) (c) Fig. 27. Gradient-based distances and WS-segmentations with a volume criterion in R = 20 regions : (a) in the factorial space c. f , (b) in the parameters space p and (c) in the PCA space of the training set on the parameters ctrain ). ß Image Anal Stereol 2015;34:1-25 mrk. mrk. mrk. f sg(.E (ctrain sg(.E (c. ),mrk. ) sg(.1/.2 (p), mrk. ) ß ), mrk. ) (a) (b) (c) Fig. 28. Segmentations by marker-controlled WS into several spaces: (a) in the factorial space cf . , (b) in the parameters space p and (c) in the PCA space of the training set on the parameters ctrain . The gradient­ ß based distances are the same as in Fig. 27. The markers come from the LDA classi.cation for the different image spaces. Probability density function for multivariate images In Noyel et al. (2007a), we studied two ways to extend the probability density function of contours to multispectral images: 1. the .rst one is a marginal approach (i.e., channel by channel) called marginal pdf mpdf (algorithm in Table 2) 2. the second one is a vectorial approach (i.e., vector pixel by vector pixel) called vectorial pdf vpdf (algorithm in Table 3). Table 2. Algorithm: mpdf 1: For the morphological gradient of each channel .( f. j ), j . [1,...,L], throw M realisations of N uniform random germs, i.e., the markers j=1...L {mrk ij}i=1...M, generating M × L realisations. Get jj=1...L the series of segmentations, {sgi (x)}i=1...M, by watershed associated to morphological gradients of each channel .( f.j ). 2: Get the marginal pdfs on each channel by Parzen method: pd f j(x)= M 1 .Mi=1 sgij(x) * G. . 3: Obtain the weighted marginal pdf: L mpd f (x)= . wjpd f j(x) , (20) j=1 with wj = 1/L, j . [1,...,L] in the image space and wj equal to the inertia axes in the factorial space. Table 3. Algorithm: vpdf 1: For the vectorial gradient .d(f. ), throw M × L realisations of N uniform random germs, i.e., the markers {mrki}i=1...M×L, with L the channels number. Get the segmentation, {sgi(x)}i=1...M×L, by watershed associated to the vectorial gradient .d (f. ), with d = d.2 in the image space or d = dE in the factorial space. 2: Obtain the probability density function: M×L vpd f (x)= 1 . sgi(x) * G. . (21) M × L i=1 A probabilistic gradient was also de.ned in Angulo and Jeulin (2007) to ponder the enhancement of the largest regions by the introduction of smallest regions. It is de.ned as .prob = mpdf + .d : after normalisation in [0,1] of the weighted marginal pdf mpdf and the metric-based gradient .d . In order to obtain a partition from the mpdf, the vpdf or the gradient .prob, these probabilistic functions are segmented, for instance by a hierarchical WS with a volume criterion, as studied in Noyel et al. (2007a). In such a case, the goal is not to .nd all the regions. The stochastic WS addresses the problem of image segmentation in few pertinent regions according to a combined criterion of contrast and size. In the present study, as we discuss below, the segmentation of the pdf is obtained from the WS with a volume criterion. CONDITIONING OF THE GERMS OF THE PDF BY A PREVIOUS CLASSIFICATION The pdf of contours with uniform random germs contains only spatial information. By conditioning the random germs by the spectral classi.cation, we introduce a spatio-spectral pdf. These germs are going to be regionalized by a pre-segmentation obtained by a pre-processing of the spectral classi.cation. Several kinds of germs have been tested: 1. uniform random point germs mrki(x) 2. random germs regionalized by a pre-segmentation: a) as point-germs mrk.-pt (x) i b) as ball-germs where: – each connected class may be hit one time mrk.-b(x) i – each connected class may be hit several times and the union of balls is made in each connected class of the pre-segmentation mrki .-.b(x) – each connected class may be hit several times and the union of connected balls is made in each connected class of the pre-segmentation mrk.-.b-conn(x). i The pdf of contours with uniform random germs is constructed without any prior information about the spatial/spectral distribution of the image. Spectral information is introduced in the pdf by conditioning the germs by the previous transformed classi.cation .f. To do this, it is possible to use point germs or random ball germs whose location is conditioned by the classi.cation. An exhaustive study of the germs is presented in Noyel (2008c) and Noyel et al. (2010). Below, we present random ball germs regionalized by a classi.cation where each connected class may be hit one time, mrk.-b(x). For the detection of tumours i in DCE-MRI series we prefer to use the last kind of regionalized random balls germs mrk.-.b-conn(x). i The procedure is as follows: the transformed classi.cation .fis composed of connected classes, .f= .kCk with Ck .Cki = /0, for k = ki. The new void class, which appears after the transformation of ., is written C0. Then random germs are drawn conditionally to the connected components Ck of the .ltered classi.cation .f. To do this, the following rejection method is used: random point germs are uniformly distributed. If a point germ m falls inside a connected component Ck of minimal area S and not yet marked, then it is kept, otherwise it is rejected. Therefore not all the germs are kept. These point germs are called random point germs regionalized by the classi.cation .. However, these regionalized point germs are sampling all the classes, independently of their prior estimate of class size/shape. In order to address this limitation, we propose to use random balls as germs. The centres of the balls are the random point germs and the radii r are uniformly distributed between 0 and a maximum radius Rmax: U [1,Rmax]. At each step, only the intersection, B(m,r) . Ck, between the ball B(m,r) and the connected component Ck is kept as a germ. Then the union is made with the previous germs. At the end of the “fall” of the random germs, the connected classes are considered as markers for the watershed used to build the pdf of contours. These balls are called random balls germs regionalized by the classi.cation . and noted mrk.-.b-conn (x). The algorithm in table 4 sketches the process. Note that if N is the number of random germs to be generated, the effective number of implanted germs is lower than N. .LDA,4 .LDA,4 ctrain fctrain mpd f (p, mrk. ) ß ß sgR-vol (mpd f ) sgR-vol (mpd f ) with R = 30 with R = 20 Fig. 29. Segmentations by stochastic pdf with a volume criterion in R = 30 (or 20) regions. The pdf mpdf is .LDA,4 conditioned by the transformed classi.cation f train c ß on the parameters space p. The parameters used to build the mpdf with regionalized random balls-germs mrk.-.b-conn are N = 100 points, M = 100 i realisations, area S = 2 pixels, Rmax = 30 pixels. Table 4. Regionalized random balls-germs: each connected class may be hit several times and the union of connected balls is made in each connected class of the pre-segmentation mrk.-.b-conn(x) i 1: Given N the number of drawn germs m, {Ck}the set of all the connected components of the transform classi.cation .f, S the minimal area of a connected component Ck, and a boolean array of size equals to the number of connected component Ck (the array values are equal to marked or not marked) 2: Set the image of germs mrk.-.b-conn(x) equals to i zero 3: Set the background class and the void class C0 to marked 4: Set the class Ck of which the area is lower than S to marked 5: for all drawn germs m from 1 to N do 6: if Ck, such as m . Ck, is not marked then 7: r = U [1,Rmax] mrk.-.b-conn 8: (x)=(B(m,r) . Ck) . i mrk.-.b-conn (x) i 9: end if 10: end for 11: Label each connected regions in the image of markers Image Anal Stereol 2015;34:1-25 After computing the marginal pdf of contours mpdf with random balls germs mrk.-.b-conn(x), from the i representation of pixels in the parameters space, the pdf are segmented by a hierarchical WS with a volume criterion (Fig. 29). The classi.cation used to build the pdf is the LDA classi.cation in the PCA space of the training set of parameters ctrain . ß In Noyel et al. (2010), the pdf built with these germs seemed to give better results than others. Therefore, the results shown use these germs. In order to better understand the process to build the regionalized random balls-germs mrk.-.b-conn , i in Fig. 30 some realisations of germs and their associated contours are presented. For the watershed segmentation, a morphological gradient is used on each channel of the PCA space of the training set of parameters. The random balls-germs are conditioned by the LDA classi.cation in this space. i = 1 i = 5 i = 9 .LDA,4 f ctrain ß i .LDA,4 by the transformed classi.cation LDA f(with N train c ß = 100 points, M = 100 realisations and Rmax = 30 pixels). VALIDATION OF THE METHOD: APPLICATION TO COMPUTER AIDED DETECTION OF TUMOURS After presenting the way to compute the stochastic WS with regionalized random balls-germs, we are going to apply it to computer aided detection of tumours on several DCE-MRI series. In order to detect potentially tumourous areas, the DCE-MRI series are .rst segmented by stochastic WS. Then the regions of the segmentation are classi.ed in potentially tumourous (or not tumourous) areas. The whole analysis .owchart is presented in Fig. 31. It combines the different parts introduced earlier in this paper: – pre-processing stage: a noise reduction by FCA and model .tting – training stage of the classi.er: the LDA classi.er is trained on some reference pixels selected on a reference image of the parameters. – classi.cation stage: a LDA after normalising the histograms of the maps of the parameters which model the spectra. The histograms are normalised in order to match the histogram of the reference image. – segmentation stage: a stochastic WS with regionalized random-balls germs mrk.-.b-conn(x) i conditioned by the classi.cation. Fig. 31. Flowchart of tumour detection. Starting from the segmentation by stochastic watershed, the detection of potential tumours depends on two criteria which have been empirically determined according to a prior knowledge on DCE­MRI series: 1. a positive mean slope parameter a because the contrast agent tends to accumulate in these areas during the acquisition 2. a mean intercept b higher than given a threshold (800) after histogram normalisation. With this parameter, the areas of the background with a small positive slope are removed from the detection. So that medical doctors can evaluate the pertinence of the detected zones, con.dence maps on the parameters were built. For each zone, some coef.cients of variation were computed. These coef.cients ß are de.ned as the ratio between the standard deviation, ., and the mean, mean, of the parameters for the considered region: ßa = .a E[a] and ßb = .b E[b] . (22) Then, the con.dence maps were thresholded: for ßa at 5 and for ßb at 1. For coef.cients close to zero, the considered region is more likely classi.ed as cancerous. A look up table is applied on the con.dence maps: in blue is the highest risk (ß = 0) and in red is the lowest risk (ßa . 5 or ßb . 1). The detection det(x) and the con.dence maps for the series “serim447” are in Fig. 32. We notice that for .LDA,4 .LDA,4 p fp = mrk. mpd f (p,mrk. ) sgR-vol (mpd f ) ref det(x) ßa ßb Fig. 32. Detection of potentially cancerous areas. .LDA,4 LDA classi.cation in the parameter space fp , morphological transform of the classi.cation used to condition the germs of the pdf mrk. , mpdf with random balls germs regionalized by the classi.cation mpd f (p,mrk. ) (N = 100 points, M = 100 realisations, area S = 10 pixels, Rmax = 30 pixels), segmentation by volumic WS in R = 20 regions sgR-vol (mpd f ), reference re f , detection det(x) of potentially cancerous zones (mean(a) > 0 and mean(b) > 800), con.dence maps for the slope ßa and for the intercept ßb. the largest potentially cancerous zone, the risk is high (in blue). This corresponds to the tumour speci.ed by the medical doctors. On the other hand, the smallest zone in red, for which the risk is low, is not a tumour. Therefore, our detection method works. The same approach has been applied to 25 series of images. In the Fig. 33 some results are presented for 6 series. The potentially cancerous detected zones correspond to the references given by the doctors. In the series “serim450” and “serim457” some potentially cancerous zones with a higher risk are even detected while they were not selected in the reference. In the others series, similar results have been obtained. All the tumours marked by clinicians have been detected. CONCLUSION In this paper, an automatic method of detection of potentially cancerous zones on DCE-MRI series is presented. The results have been tested on a limited number of images. They are very promising and in agreement with the references given by medical doctors. Our method is composed of four stages. In the pre-processing stage, a dimensionality reduction and a noise reduction are .rst performed with the Factor Correspondence Analysis and the pixel-based spectrum modeling. These operations preserve the spatial contours of the image. Then in the classi.cation stage Linear Discriminant Analysis is performed on a subset of training pixels. In the third stage, the image is segmented by stochastic watershed with random-balls markers regionalized by the previous classi.cation. The originality of this approach is the combination of the spatial and temporal information to produce a “multivariate gradient” representing the probability density function of contours. These probability maps are segmented by stochastic watershed which is very useful when segmenting the low contrasted regions corresponding to tumours, since it regularises the contours. The last stage is a detection of potentially tumourous zones by statistical criteria. A con.dence maps is associated to the selected zones. These maps show the risk of the regions to be cancerous. A method of computer aided detection of potentially cancerous zones on DCE MRI sequences of small animal has been detailed. It seems to be very promising for low contrasted data sets. Further systematic tests should be performed, in order to validate the method on a larger data sets. In the future, some physical models could be .tted on the temporal series in place of the actual model. Image Anal Stereol 2015;34:1-25 re f (x) mpd f (x) det(x) ßa(x) ßb(x) “serim1441” Fig. 33. Detection results: Reference ref, marginal pdf mpdf (p,mrk.LDA )(x), detection det(x) of potentially cancerous zones in the parameters space after histogram normalisation (mean(a) > 0 and mean(b) > 800), con.dence maps on the slope ßa and on the intercept ßb. APPENDIX EXPERIMENTAL CONDITIONS As explained in Brochot et al. (2006), here are the details about experimental conditions: animals used and MRI examination. Animals. Experiments were performed on nude nu/nu male mice (Laboratoire Iffa Credo, L’Arbresle, France), in full compliance with the National Institutes of Health recommendations for animal care. Approximately 1.5 × 106 PC-3 human tumour cells were implanted subcutaneously into the .ank of each mouse, as described in Pradel et al. (2003). For imaging, the animals were anesthetised by a peritoneal injection of ketamine (Rompun, Bayer, Leverkusen, Germany) and xylazine (Imalgene, Mérial, Lyon, France). MRI examination. MRI examination was performed using a 1.5-T system (Sigma, General Electrics, Milwaukee, WI, USA) and a custom small-animal dedicated coil. A sagittal 2D T1-weighted spin echo sequence (TE 11 ms, TR 400 ms, FOV 8×8 cm, 256 ×128 matrix, 1 NEX) was used to check adequate positioning of the animal and to select the axial plane level containing the left ventricle cavity and one or two .ank tumours. The dynamic acquisition was performed using a single-slice T1-weighted 2D fast spoiled gradient recalled (FSPGR) sequence: TR 15 ms, TE 2.2 ms, .ip angle of 608, bandwidth 31.25 kHz, 256×76 matrix for the asymmetric FOV of 7 × 3 cm, 5 mm slice. The single slice was positioned at the level selected by the previous sagittal sequence, and dynamic acquisition was performed with 10 baseline images and after a caudal vein bolus injection of 0.045 mmol Gd/kg of a macromolecular contrast agent (Vistarem, Guerbet, Aulnay-Sous-Bois, France). NOISE REDUCTION In this section, more details are given about our method of noise reduction based on two sequences of FCA and reconstruction (see section “Noise reduction” on page 7). This approach needs the subtraction of a constant . as shown in Eqs. 7, 8 and 9. The idea of subtracting a constant from the data is based on the fact that after one sequence of FCA-reconstruction some values of the reconstructed image ff(1) turn out negative. This is due to the fact that only . a limited number of the factorial axes are kept for the reconstruction. These negative values, however, have no physical meaning. They are corrected in the .rst reconstructed image. Then a second sequence of FCA-reconstruction is applied because the data set has been modi.ed. In order to be consistent, the constant is added after the second sequence of FCA-reconstruction also. The data are classi.ed by k-means on the factor pixels of the .rst FCA cf(1) and on the factor pixels of . the second FCA cf(2) . One can notice, in Fig. 34, that . the classi.cation is better with two FCA than with one FCA. (a) (b) Fig. 34. Classi.cations kmeans into 5 classes in the space of the factor pixels (a) of FCA 1, cf(1) , and (b) of . FCA 2, cf(2) . . Starting from the sixteen retained factorial axes of the .rst FCA (see section 3.1.2), the image is partially reconstructed and a second sequence of FCA-reconstruction is applied with a spectral translation (equation 7). In order to verify the importance of two sequences of FCA-reconstruction, the SNR are estimated: on the channels of the original image f. , on the channels of f(1) the .rst reconstructed imagef. and on the channels of f(2) the second reconstructed image f. . The SNR are also estimated on the factor pixels of the .rst FCA cf(1) and . of the second FCA cf. (2) (.g. 35). After performing the .rst sequence of FCA-reconstruction, an improvement of the SNR is noticed in the image space. However, the second sequence of FCA does not improve the SNR in the image space. In the factor space, the SNR is improved after the second FCA in comparison with the SNR, in the factor space, after the .rst FCA. Why is it necessary to apply two FCA to .lter the noise in the factor space, while only one FCA is necessary in the image space? During the reconstruction stage the factor pixels are weighted by their inertia. Therefore the weight of the noise is reduced because it appears on the factor pixels with a small inertia. Moreover, the image reconstruction is based on the product of the marginal frequencies, .i... j, which corresponds to the reconstruction of the barycentre. This barycentre gives the general appearance of the image. However, in order to remove the noise on the factor pixels two FCA are necessary. Image Anal Stereol 2015;34:1-25 (b) Fig. 35. (a) SNR on the channels of the image. The red curve and the blue curve are superimposed. (b) SNR on the factor pixels for different sequences of FCA-reconstruction. A hyperspectral signal to noise ratio between the original image (with noise) and the reconstructed image (.ltered) may be de.ned as the ratio between the sum of the variance of the signal for each channel and the sum of the variance of the noise for each channel: .Lj=1 var( ff. j )SNRhyper(f. )= . (23).Lj=1 var( f. j - ff. j ) REFERENCES Angulo J, Jeulin D (2007). Stochastic watershed segmentation. In: Mathematical morphology and its applications to signal and image processing. Proc 8th Int Symp Math Morpho. Instituto Nacional de Pesquisas Espaciais (INPE), 265–76. Balvay D, Frouin F, Calmon G, Bessoud B, Kahn E, Siauve N, Clément O, Cuenod CA (2005). New Criteria for assessing .t quality in dynamic contrast-enhanced T1­weighted MRI for perfusion and permeability imaging. Magnet Reson Med 54:868–77. Balvay D, Kachenoura N, Espinoza S, Thomassin-Naggara I, Fournier LS, Clément O, Cuenod CA (2011). Signal-to-noise ratio improvement in dynamic contrast-enhanced CT and MR imaging with automated principal component analysis .ltering. Radiology 258(2):435– 45. Benediktsson JA, Palmason JA, Sveinsson JR (2005). Classi.cation of hyperspectral data from urban areas based on extended morphological pro.les. IEEE T Geosci Remote 42:480–91. Benzécri JP (1964). Sur l’Analyse Proximités. Publications Institut Université de Paris 13:235–82. Facde torielle Statistique­ des Benzécri JP (1973). L’Analyse Des Données, L’Analyse des Correspondances. Dunod-Paris 2. Berman M (1985). The statistical properties of three noise removal procedures for multichannel remotely sensed datas. CSIRO Div. of Mathematics and Statistics-Australia. Consulting Report NSW/85/31/MB9. Beucher S, Lantuéjoul C (1979). Use of watersheds in contour detection. In: Proc. Int. Workshop on image processing, real-time edge and motion detection-estimation 17–21. Beucher S, Meyer F (1992). The morphological approach to segmentation: The watershed transformation. In: Mathematical morphology in image processing. Marcel-Dekker, New York 433–81. Brasch RC, Li KCP, Husband JE, Keogan MT, Neeman M, Padhani AR et al.(2000). In vivo monitoring of tumor angiogenesis with MR imaging. Acad Radiol 7(10):812–23. Brix G, Ravesh MS, Zwick S, Griebel J, Delorme S (2012). On impulse response functions computed from dynamic contrast-enhanced image data by algebraic deconvolution and compartmental modeling. Phys Medica 28(2):119–28. Brochot C, Bessoud B, Balvay D, Cuénod CA, Siauve N, Bois FY. Evaluation of antiangiogenic treatment effects on tumors’ microcirculation by Bayesian physiological pharmacokinetic modeling and magnetic resonance imaging (2006). Magn Reson Imaging 24:1059—67. Cattell RB (1966). The scree test for the number of factors. Multivar Behav Res 245–76. Diday E (1971). Une nouvelle méthode en classi.cation automatique et reconnaissance des formes la méthode des nuées dynamique. Rev Stat Appl 19(2):19-33. Ding Y, Chung YC, Raman SV, Simonetti OP (2009). Application of the Karhunen-Loeve transform temporal image .lter to reduce noise in real-time cardiac cine MRI. Phys Med Biol 54(12):3909–22. Ding Y, Chung YC,Simonetti OP (2010). A method to assess spatially variant noise in dynamic MR image series. Magnet Reson Med 63(3):782–89. Duda RO, Hart PE (1973). Pattern classi.cation and scene analysis. Wiley-New York. Green A, Berman M, Switzer P, Craig MD (1988). A transformation for ordering multispectral data in terms of image quality with implications for noise removal. IEEE T Geosci Remote 26(1):65—74. Hanbury A, Serra J (2001). Morphological operators on the unit circle. IEEE T Image Process 10(12):1842–50. Hartigan JA, Wong MA (1979). A k-means clustering algorithm. Appl Stat 28:100–8. Hastie T, Tibshirani R, Friedman J (2003). The elements of statistical learning: Data mining, inference, and prediction. New York: Springer. Hughes G (1968). On the mean accuracy of statistical pattern recognizers. IEEE T Inform Theory 14:55–63. Ivancevic MK, Zimine I, Lazeyras F, Foxall D, Vallée JP (2001). FAST sequences optimization for contrast media pharmacokinetic quanti.cation in tissue. J Magn Reson Im 14(6):771–8. Kaiser H (1960). The application of electronic computers to factor analysis. Educ Psychol Meas 20:141–51. Landgrebe D (2002). Hyperspectral image data analysis. IEEE Signal Proc Mag 19(1):17–28. Leach MO, Morgan B, Tofts PS, Buckley DL, Huang W, Hors.eld MA et al.(2012). Imaging vascular function for early stage clinical trials using dynamic contrast-enhanced magnetic resonance imaging. Eur Radiol 22(7):1451-64. Lennon M (2002). Méthodes d’analyse d’images hyperspectrales. Exploitation du capteur aéroporté CASI pour des applications de cartographie agro­environnementale en Bretagne. Ph.d. Thesis-Université de Rennes, France. Matheron G (1970). La théorie des variables régionalisées et ses applications. Les cahiers du Centre de Morphologie Mathématique de Fontainebleau. Mines-ParisTech, France. http://www.cg.ensmp.fr/bibliotheque/index_ byyear.html Matheron G (1975). Random Sets and Integral Integral Geometry. Wiley, New York. Meyer F (2001). An overview of morphological segmentation. Int J Pattern Recogn 15(7):1089–118. Meyer F (2004). Levelings, image simpli.cation .lters for segmentation. J Math Imaging Vis 20:59–72. Noyel G, Angulo J, Jeulin D (2007). Random germs and stochastic watershed for unsupervised multispectral image segmentation. In: Apolloni B, Howlett RJ, Jain L, eds. Proc 11th Int Conf KES 2007/ WIRN 2007. Lect Not Comput Sci 4694:17–24. Noyel G, Angulo J, Jeulin D (2007). Morphological segmentation of hyperspectral images. In: Proc 12th Int Congr Stereol (ICS XII). Saint Etienne, France. Noyel G, Angulo J, Jeulin D (2007). Morphological segmentation of hyperspectral images. Image Anal Stereol 26:101–9. Noyel G, Angulo J, Jeulin D (2007). On distances, paths and connections for hyperspectral image segmentation. In: Mathematical morphology and its applications to signal and image processing, Proc 8th Int Symp Math Morpho. Instituto Nacional de Pesquisas Espaciais (INPE) 1:399–10. Noyel G, Angulo J, Jeulin D (2008). Classi.cation-driven stochastic watershed. Application to multispectral segmentation. In: Proc 4th Eur Conf Color Graph Imaging Vision (CGIV 2008) 471–6. Noyel G, Angulo J, Jeulin D (2008). Filtering, segmentation and region classi.cation by hyperspectral mathematical morphology of DCE-MRI series for angiogenesis imaging. In: Proc 5th IEEE Int Symp Biomed Imaging (ISBI 2008) 1517–20. Noyel G (2008). Filtrage, réduction de dimension, classi.cation et segmentation morphologique hyperspectrale. Mines-ParisTech, France. Noyel G, Angulo J, Jeulin D (2010). Regionalized random germs by a classi.cation for probabilistic watershed. Application: angiogenesis imaging segmentation. In: Progress in Industrial Mathematics at ECMI 2008. Math Indust 15:211–16. Noyel G, Angulo J, Jeulin D (2011). A new spatio-spectral morphological segmentation for multi-spectral remote-sensing images. Int J Remote Sens 31(22):5895–920. O’Connor JPB, Jackson A, Asselin MC, Buckley DL, Parker GJM, Jayson GC (2008). Quantitative imaging biomarkers in the clinical development of targeted therapeutics: current and future perspectives. Lancet Oncol 9(8):766–76. Orfeuil JP (1973). Une méthode de .ltrage des données. CG, Mines-ParisTech. Report N-314. Pradel C, Siauve N, Bruneteau G, et al (2003). Reduced capillary perfusion and permeability in human tumour xenografts treated with the VEGF signalling inhibitor ZD4190: an in vivo assessment using dynamic MR imaging and macromolecular contrast media. Magn Reson Imaging 21(8):845—51. Scheunders P (2002). A multivalued image wavelet representation based on multiscale fundamental forms. IEEE T Image Process 11(5):568–75. Soille P (1999). Morphological image analysis. Berlin: Springer. Sourbron SP, Buckley DL (2012). Tracer kinetic modelling in MRI: estimating perfusion and capillary permeability. Phys Med Biol 57(2):R1–33. Image Anal Stereol 2015;34:1-25 Zvan Dijke CF, Brasch RC, Roberts TP, Weidner N, Mathur A, Shames DM et al.(1996). Mammary carcinoma model: correlation of macromolecular contrast-enhanced MR imaging characterizations of tumor microvasculature and histologic capillary density. Radiology 198(3):813–8. Zahra MA, Hollingsworth KG, Sala E, Lomas DJ, Tan LT (2007). Dynamic contrast-enhanced MRI as a predictor of tumour response to radiotherapy. The Lancet Oncology 8(1):63–74.