Image Anal Stereol 2018;37:21-34 doi:10.5566/ias.1624 Original Research Paper ARFBF MORPHOLOGICAL IMAGE ANALYSIS – APPLICATION TO THE DISCRIMINATION OF CATALYST ACTIVE PHASES ZHANGYUN TAN1, MAXIME MOREAUDB,2, OLIVIER ALATA3 AND ABDOURRAHMANE M. ATTO1 1LISTIC, EA 3703, University Savoie Mont Blanc, France; 2IFP Energies Nouvelles, rond-point de l’echangeur de Solaize, BP 3, 69360 Solaize, France; 3Lab. Hubert Curien, UJM-Saint-Etienne, CNRS UMR 5516, IOGS, Université de Lyon, 42023 Saint-Etienne, France e-mail: zhangyun.tan@univ-paris13.fr, maxime.moreaud@ifpen.fr, olivier.alata@univ-st-etienne.fr, abdourrahmane.atto@univ-smb.fr (Received September 15, 2016; revised November 14, 2017; accepted December 18, 2017) ABSTRACT This paper addresses the characterization of spatial arrangements of fringes in catalysts imaged by High Resolution Transmission Electron Microscopy (HRTEM). It presents a statistical model-based approach for analyzing these fringes. The proposed approach involves Fractional Brownian Field (FBF) and 2-D Auto- Regressive (AR) modeling, as well as morphological analysis. The originality of the approach consists in identifying the image background as an FBF, subtracting this background, modeling the residual by 2-D AR so as to capture fringe information and, finally, discriminating catalysts from fringe characterizations obtained by morphological analysis. The overall analysis is called ARFBF (Auto-Regressive Fractional Brownian Field) based morphology characterization. Keywords: auto-regressive field, fractional Brownian field, HRTEM imaging, mathematical morphology, texture analysis. INTRODUCTION Texture analysis plays a key role in almost all imaging systems (radar, sonar, X-ray, microscopy. . . ) and human visual perception. Texture structure characterization has deserved a considerable amount of works in order to describe image patterns for feature recognition (see Qazi et al., 2011; Goncalves and Bruno, 2013, among other references). In recent literature, material texture characterization at nanoscale using image analysis has been considered either qualitative or quantitative methods. For instance, Da Costa et al. (2015) investigated carbon planes with orientation description, Pré et al. (2013) considered activated carbon with mathematical morphology approach, Moreaud et al. (2012) studied alumina platelets using morphological random model- based approach, Toth et al. (2013) investigated soot with multi-resolution approach, Zhang et al. (2011) considered active phases of unsupported hydrodesulphurization catalyst and Moreaud et al. (2008) investigated ceria active phase using local frequency segmentation approach and morphological analysis. In this paper, we are interested in the characterization of hydrotreating catalysts with sulphide phases supported on alumina (see Toulhoat and Raybaud 2013, Celce et al., 2008) using Fig. 1. Molybdenum sulfide sheets (black fringes) deposited on an alumina support (top image) and its sub-images (significant fringes selected, in the second row). The initial image consists of 1024×1024 pixels. The numerical resolution is 0.057 nm by pixel. 21 TAN Z ET AL: ARFBF morphological image analysis statistical model-based approach. Samples were observed by high resolution transmission electron microscopy, using a JEOL TEM 2100F operated at 200kV, whose point-to-point resolution is 0.23 nm. The so-called CoMoS1 active phase consists of nanolayers of MoS2, decorated at the edges with Co atoms. When their basal plane is oriented parallel to the electron beam axis, then they appear as black fringes. They can be either isolated of stacked with a stacking number up to 5 usually. If the catalyst presents a high loading of active phase, nanolayers can also form aggregates filling the support porosity. Length and stacking number of the nanolayers depend on several parameters: nature of the support, impregnation and sulfidation methods. They directly impact the activity and selectivity of the catalyst (see Toulhoat and Raybaud, 2013 and Nikulshin et al., 2014). Active phases shown in Fig. 1 can be seen as pseudo-periodic fields, a property which can be captured by using a 2-D Auto-Regressive model (2-D AR). Indeed, the 2-D AR random field model has shown relevancy for describing a wide class of pseudo- periodic textures (see Haralick, 1979; Souza, 1982; Mao and Jain, 1992; Alata and Ramananjarasoa, 2005 for the relevancy of this model in classification, segmentation and recognition of textural information). We thus consider the AR model for characterizing the active phases and propose using AR field power spectrum in order to derive their features. The spectrum of 2-D AR model can be calculated by the Harmonic Mean (HM) method (see Jackson and Chien, 1979; Alata et al., 1998). The nanolayers are deposited on a crystalline alumina support (Fig. 1). Observations at high resolution also produce fringes (corresponding to atomic planes) and Moire fringes. Nevertheless, their contrast and their periodicity differ from those of CoMoS nanolayers, preventing confusion. Moreover, post-treated areas correspond to the thinnest zones of the sample, where alumina contribution is minored. This background is considered here as a noisy field interacting with the intrinsic AR model and it has to be removed or attenuated significantly so as to make AR modeling efficient. For this purpose, we propose using a 2-D Fractional Brownian Field (FBF) for modeling this undesirable interaction, prior to AR modeling. FBF is the spatial extension of the fractional Brownian motion (fBm), a non-stationary process derived by Mandelbrot and Ness (1968). Several variants of FBF exist, among which the separable and isotropic extensions. We consider hereafter the isotropic FBF whose characterizations are given in Pesquet-Popescu and Larzabal (1997); Pesquet- Popescu and Véhel (2002); Huang and Li (2005); Richard and Bierme (2010); Atto et al. (2014). Isotropic FBF is a one-parameter model. Its parameter, called Hurst parameter, is representative of stochastic regularity (texture roughness, in practice). In this respect, it suits for the description of the heterogeneous HRTEM image background. In the literature, several methods exist for estimating the Hurst parameter of the 2-D FBF, for instance, maximum likelihood estimates (see Fieguth and Willsky, 1996; Tafti et al., 2009), box- counting approach (see Huang et al., 1994) and log- periodogram methods (see Geweke and Porter-Hudak, 1983; Robinson, 1995). Because FBF admits one pole located at zero frequency in the spectral domain (infinite spectral value at zero), these estimators lack robustness on small sample sizes. In order to obtain a robust Hurst parameter estimator, we consider the extension to FBF (see Tan et al., 2015) of the wavelet packet method of Atto et al. (2010) dedicated to fractional Brownian motion Hurst parameter estimation. This extension, based on 2-D wavelet packet spectrum, is called Log-Regression on Polar representation of Wavelet Packet spectrum (Log- RPWP). In this paper, we propose an ARFBF based morphology analysis so as to capture HRTEM fringe information. This method relies on: 1. the integration of 2-D FBF and 2-D AR by using a convolution operator in order to define HRTEM ARFBF (Auto-Regressive Fractional Brownian Field, see Tan et al., 2015) features and 2. morphological analysis on the basis of the HRTEM ARFBF features. Summarizing, the originality of the approach consists in identifying the image background as an FBF, subtracting this background, modeling the residual by a 2-D AR which is able to characterize the pseudo- periodic structure of the active phase deposited on the support of catalyst. The contributions presented in this work involve: – assessing, in comparison with the state of the art (see Van De Ville et al., 2005; Tafti et al., 2009), the performance of Log-RPWP Hurst parameter estimation method introduced in Tan et al. (2015) on different FBF generators, – deriving an analysis framework by integrating morphological analysis of the ARFBF description, 1CoMoS refers to catalysts consisting of nanolayers of MoS2, decorated at the edges with Co atoms. 22 Image Anal Stereol 2018;37:21-34 – addressing HRTEM micrograph morphological texture characterization by using this ARFBF morphological framework, – discriminating different catalysts from obtained characterizations. A comparison is proposed with a Fourier-based spectrum analysis in order to show the interest of using model-based HM method. In the next section, we present the 2-D ARFBF model after some recalls about 2-D AR model and FBF model, and after the new results about Hurst parameter estimation. In section “Results”, 2-D ARFBF morphological analysis of HRTEM is detailed and experimental results about catalyst discrimination are provided. In section “Conclusion”, the conclusion of this study is provided with some perspectives. MATERIALS AND METHODS ARFBF MODEL The 2-D ARFBF introduced in Tan et al. (2015) is defined as the convolution with respect to deterministic spatial variables, of AR field A and isotropic FBF BH . This 2-D ARFBF, hereafter denoted Z(x,y), is defined as follows, Z(x,y) = (A∗BH)(x,y) . (1) Particularly, Z behaves as an FBF when AR is a white noise; Z is an AR when FBF is a white noise (H = 0). When an image is associated to an ARFBF observation Z, this means that the image content can be seen as an AR based filtering of a fractional Brownian texture. The model applies on spatial texture level by describing the spatial covariance structure, in contrast with pixel level models such as marked point processes. The following Section “2-D AR model” presents the 2-D AR framework used in the paper. Then, Section “2-D FBF isotropic model and Hurst parameter estimation” presents the isotropic FBF, as well as an estimation procedure for the parameter of this model, in addition with a performance validation stage (main contribution of the Section) through a comparison with respect to the state of the art. 2-D AR MODEL Let us define a centered second-order stationary field as A = {A(x,y)} ,(x,y) ∈ Z2. A is a 2-D AR process if A(x,y) =− ∑ m∈D βmA(x−m1,y−m2)+E(x,y) , (2) where D is the prediction support, m = (m1,m2) ∈ D, E(x,y) is an independent and identically distributed process with zero mean and variance σ2e , βm are 2-D transverse coefficients. Different prediction supports can be defined (Alata and Olivier, 2003). In this paper, we will use two prediction supports with Quarter Plan (QP) forms respectively denoted QP̀ , ` = 1,2 (Fig. 2 shows an illustrative example of these QPs), with: DQP1 = {0≤m1 ≤M1,0≤m2 ≤M2}\{(0,0)} , (3) and DQP2 = {−M1 ≤ m1 ≤ 0,0≤ m2 ≤M2}\{(0,0)} . (4) (a) (b) Fig. 2. Quarter plan prediction supports denoted (a) DQP1 (Eq. 3) and (b) DQP2 (Eq. 4) with finite order (M1,M2) = (2,2) and (x,y) ∈ Z2. 23 TAN Z ET AL: ARFBF morphological image analysis (M1,M2) is the 2-D AR order. Let us denote the set of AR parameters associated with QP̀ , ` ∈ {1,2}, as θ QP̀M1,M2 = {σ 2 e,`,{βm,`,m ∈ DQP̀ }} . (5) At a fixed order, the parameters of θ QP̀M1,M2 , ` ∈ {1,2}, are estimated thanks to minimization in the least-squares sense of prediction errors, see Ranganath and Jain (1985); Alata and Cariou (2008a). This procedure involves Yule-Walker equations and it is equivalent to the maximum likelihood estimation, when the random variables are Gaussian. In practice, the selection of an accurate prediction support determines model performance. Although a number of works focus on the order (M1,M2) selection, see references Akaike (1974); Alata and Olivier (2003) among others, this paper considers, from preliminary experimental model validation, M1 = M2 = 10 as the order of the prediction support. The 2-D AR Power Spectral Density (PSD), SAR(u,v), is then derived as the harmonic mean of the two spectra associated with prediction supports QP̀ , `= 1,2 (Jackson and Chien, 1979): SAR(u,v) = 2SQP1(u,v)SQP2(u,v) SQP1(u,v)+SQP2(u,v) , (6) where SQP̀ (u,v) = σ2e,` |FQP̀ (u,v)|2 (7) and FQP̀ (u,v) = 1+ ∑ m=(m1,m2)∈DQP̀ βm,`e−i2πum1e−i2πvm2 . The 2-D spectrum estimated from this method is easy to compute and has good estimation properties with respect to the other existing methods, see Alata and Cariou (2008b) for details. 2-D FBF ISOTROPIC MODEL AND HURST PARAMETER ESTIMATION 2-D FBF model The 2-D isotropic FBF with Hurst parameter H, 0 < H < 1, denoted here as BH(x,y), is defined to be a non-stationary Gaussian zero-mean real-valued field with auto-correlation function defined as RBH (t,s) = σ2b 2 (‖t‖2H +‖s‖2H −‖t− s‖2H) , (8) where σ2b is a constant representing the variance of a white Gaussian noise. Although FBF is a non-stationary process, this random process has stationary-increments and stationary wavelet projections. The spectrum of FBF can then be defined by association with respect to the above stationary FBF instances and is given by (see Pesquet-Popescu and Véhel (2002); Huang and Li (2005); Richard and Bierme (2010); Atto et al. (2014)): SFBF(u,v) = ξ (H) 1 (u2 + v2)H+1 = ξ (H) 1 ‖(u,v)‖2H+2 , (9) where ‖(u,v)‖= √ u2 + v2, ξ (H) = 2−(2H+1)π2σ2b sin(πH)Γ2(1+H) and Γ is the standard gamma function. Generator #1 Generator #2 H = 0.2 H = 0.5 H = 0.8 Fig. 3. Examples of sample FBFs used for Monte- Carlo simulations in Table 1. 24 Image Anal Stereol 2018;37:21-34 Table 1. Mean values and standard deviations for estimated Hurst parameters from Monte-Carlo simulations with 10 FBF realizations. Sample FBF Log-RPWP Log-RPHW ML-PHW H Ĥ std Ĥ std Ĥ std 0.2 0.210 0.031 0.204 0.009 0.200 0.006 Generator #1 0.4 0.427 0.020 0.392 0.012 0.392 0.009 0.6 0.643 0.039 0.600 0.010 0.597 0.007 0.8 0.877 0.017 0.792 0.006 0.793 0.004 0.2 0.156 0.018 -0.029 0.004 -0.088 0.005 Generator #2 0.4 0.424 0.026 0.264 0.014 0.226 0.010 0.6 0.579 0.033 0.452 0.024 0.362 0.035 0.8 0.815 0.048 0.415 0.029 0.313 0.023 For Hurst parameter estimation, we propose the Log-RPWP method introduced in Tan et al. (2015). Log-RPWP Hurst parameter estimation method consists of three steps. – In the first step, the spectrum with polar coordinates Sp is computed as, Sp(r,θ) = T (ŜFBF(u,v)) , (10) where ŜFBF(u,v) is the wavelet packet spectrum (Atto et al. (2013)), in Cartesian coordinates (u,v), estimated from FBF samples and T is the Cartesian-to-polar transform. – In the second step, averages are done over the angles: Sp(ri) = 1 J J ∑ j=1 Sp(ri,θ j) , (11) with 1≤ i≤ N denoting the radial sampling index. – In the third step, H is estimated by: ĤRPWP = 1 2C ∑1≤i,k≤N i α1, 0, otherwise. (17) If a second lobe is located on S∗, it can be obtained using a similar procedure. The contribution of the first lobe is removed on S∗ using morphological dilation with a disc of radius R on Q1: S∗1(r,θ) = { S∗(r,θ), if δR(Q1) 6= 0, 0, otherwise. (18) From available HRTEM data (see samples given in Fig. 1 for illustration), R = 2 is efficient. Then a new threshold α2 is calculated and a binary image Q2 is obtained by: α2 = λ2S∗(P1) , (19) 3PSD estimated using periodogram which is based on the square modulus of discrete Fourier transform. The origin is at the center of the images which exhibit central symmetry. 27 TAN Z ET AL: ARFBF morphological image analysis with λ2 ∈ [0;λ1], P2 = argmax r,θ S∗1(r,θ) , (20) Q2 = γ S∗∗1 rec (P2) , (21) and S∗∗1 (r,θ) = { 1, if S∗1(r,θ)> α2, 0, otherwise. (22) In our application, we chose λ2 = 0.6 (selection from validation tests). In contrast with Q1, Q2 value is not necessarily non-zero (if there is no remaining significant second lobe). We limit the study to the detection of two lobes in this experimental setup, with the knowledge that observation of more than two layers of overlapping fringes is rare in our application. For each lobe on images Qi, i = 1 or 2, denoting a layer of fringes, an average spatial distance between atomic layers is estimated by the distance value G, G = Te r∗ , (23) associated to the maximum of the lobe Pi = [r∗,θ ∗], i = 1 or 2, Te is the sampling period (here, Te = 0.057). Lobe on binary image Qi, i = 1 or 2, can be considered as an ellipse embedded in a bounding box (rmin,rmax,θmin,θmax) (Fig. 5). The extension of the lobe allows us to describe the changes in the distance between atomic layers (regularity of spacing) and in the curvature of atomic layers (regularity of curvature). The regularity of spacing and regularity of curvature can be estimated by distance variation4G (Eq. 24) and tangential length Lθ (Eq. 27) respectively. 4G and Lθ are defined by: Fig. 5. The spectral lobe detected can be considered as an ellipse embedded in a bounding box (rmin,rmax,θmin,θmax). 4G = |Gmax−Gmin| , (24) with Gmax = Te rmin , (25) Gmin = Te rmax , (26) where Te = 0.057, and Lθ = |θmax−θmin|. (27) In the next section, we present our statistical analysis for the discrimination of catalyst active phases. STATISTICAL ANALYSIS FOR CATALYST DISCRIMINATION Based on the geometric features described previously, a comparison between two catalysts is presented. The point-to point resolution of the TEM is 0.23 nm. The pixel size of the images is 0.057 nm, at an image resolution of 1024× 1024. For catalyst 1 (Cat1), 78 sub-images containing active phases are taken from 21 HRTEM images. For catalyst 2 (Cat2), 88 sub-images containing active phases are taken from 19 HRTEM images. For each sample, we calculate geometric features G,4G and Lθ on detected lobes. Cat 1 is a CoMoS catalyst supported on alumina. It presents a high loading of active phase, corresponding to a Mo density of 7 at/nm2. Cat 2 is a CoMoS catalyst supported on silica. It presents a low loading of active phase, with a Mo density of 1 at/nm2. The mean length and stacking of slabs were measured by HRTEM on a minimum of 200 slabs, using an in-house image processing application (Celse et al., 2008): the image is pre-processed in order to stretch contrast, then a region growing algorithm and active contour techniques are used to obtain slab contours. Cat1 presents a mean slab length of 3.89 nm and a mean stacking number of 2.73; Cat 2 presents a mean slab length of 2.49 nm and a mean stacking number of 2.09. RESULTS STATISTICAL DISTRIBUTIONS OF G, 4G AND Lθ The distributions of the features are studied by using a Kernel smoothing function fw which is defined as follows (Rosenblatt, 1956; Parzen, 1962): fw(x) = 1 nw n ∑ i=1 K( x− xi w ) , (28) 28 Image Anal Stereol 2018;37:21-34 Table 2. Using AR-based method for PSD estimation in polar coordinates (Fig. 4 & 5), statistics of distance (Eq. 23), distance variation (Eq. 24) and tangential length (Eq. 27) features of the detected lobe of catalyst image databases (Cat1 and Cat2). Catalyst (n) Cat1 (93) Cat2 (109) Stat G 4G Lθ G 4G Lθ Mean 0.62 0.124 12.323 0.598 0.097 9.988 Variance 0.003 0.002 27.682 0.003 0.002 26.653 Skewness -0.991 1.222 1.116 -1.135 1.721 1.425 Kurtosis 7.21 4.784 4.572 4.434 9.681 5.913 Table 3. Using interpolated periodogram in polar coordinates (Fig. 7): statistics of distance (Eq. 23), distance variation (Eq. 24) and tangential length (Eq. 27) features of the detected lobe of catalyst image databases (Cat1 and Cat2). Catalyst (n) Cat1 (93) Cat2 (109) Stat G 4G Lθ G 4G Lθ Mean 0.63 0.067 13.296 0.61 0.056 12.776 Variance 0.002 0.0006 10.982 0.003 0.0005 13.879 Skewness -0.56 2.155 1.000 -1.254 1.745 0.634 Kurtosis 5.987 9.525 4.036 4.898 8.786 3.636 Fig. 6. Kernel distributions of distance (left), distance variation (center) and tangential length (right) of the detected lobe S∗LobeD of Cat1 (− in blue) and Cat2 (− in red). First row: using PSD estimated in polar coordinates with HM method based on the AR model. Second row: Using interpolated periodogram in polar coordinates. where n is the sample size, {xi}i=1,...,n are the set of considered samples, x ∈ R, K(·) is the kernel smoothing function and w is the bandwidth. Table 2 gives certain statistic information of distance G (measured in nm), distance variation 4G (measured in nm) and tangential lengths Lθ (measured in degree) in terms of means, variances, third and fourth standardized moments (known as skewness and kurtosis) of the detected lobes on Cat1 and Cat2 29 TAN Z ET AL: ARFBF morphological image analysis Table 4. Kolmogorov-Smirnov test for the discrimination of Cat1 and Cat2 based on features computed using PSD estimated with AR-based method in polar coordinates. Null hypothesis is rejected at the level α = 0.05 if K > Ks with Ks ≈ 1.36 √ (N1 +N2)/(N1×N2). With N1 = 78 (number of samples of Cat1) and N2 = 88 (number of samples of Cat2), Ks ≈ 0.212. G 4G Lθ Statistic Catk \Cat` Cat1 Cat2 Cat1 Cat2 Cat1 Cat2 H Cat1 0 0 0 1 0 1 Cat2 0 0 1 0 1 0 103K Cat1 0 191 0 316 0 259 Cat2 191 0 316 0 259 0 Table 5. Kolmogorov-Smirnov test for the discrimination of Cat1 and Cat2 based on features using interpolated periodogram in polar coordinates. Null hypothesis is rejected at the level α = 0.05 if K > Ks with Ks ≈ 1.36 √ (N1 +N2)/(N1×N2). With N1 = 78 (number of samples of Cat1) and N2 = 88 (number of samples of Cat2), Ks ≈ 0.212. G 4G Lθ Statistic Catk \Cat` Cat1 Cat2 Cat1 Cat2 Cat1 Cat2 H Cat1 0 0 0 1 0 0 Cat2 0 0 1 0 0 0 103K Cat1 0 204 0 295 0 127 Cat2 204 0 295 0 127 0 images. Catalyst (n) represents the number of samples in Cat1 and Cat2 respectively. For these catalysts, the inter-distance between two neighboring white/black fringes G is known to be close to 0.615 nm (see Geantet and Sorbier (2012)). In Table 2, G has almost same value for Cat1 and Cat2 (0.611 nm vs 0.593 nm), it confirms G as a physical characterization. For 4G and Lθ of Cat1 and Lθ of Cat2, the positive skewness values in Table 2 indicates that the tail on the right side is longer or fatter than that of the left side. When the skewness values are bigger, this phenomenon will be clearer (Fig. 6). In contrast, for the other cases, the tail on the left side is longer or fatter than that of the right side (Fig. 6) because of negative skewness values. The kurtosis measures provided ”tailedness” information of the distribution of variables G, 4G or Lθ . All kurtosis being larger than 3, one can conclude that G (for both of Cat1 and Cat2) and 4G (just for Cat1) distributions are far from being Gaussian. In addition, when kurtosis are large, this implies that several outliers can be present in the corresponding variables. These observations are confirmed by Fig. 6 which shows the statistical distributions of G,4G or Lθ for both catalysts. KOLMOGOROV-SMIRNOV TEST FOR CATALYST DISCRIMINATION In this section, we propose Kolmogorov-Smirnov (KS) test for comparing Cat1 and Cat2 ARFBF morphology characterizations. The KS measure is given by (see (Peacock, 1983) for details) K = maxx(| f̂1(x)− f̂2(x)|), (29) where f̂1, f̂2 are the empirical cumulative distribution functions of catalyst datasets indexed by ‘1’ and ‘2’. Function f̂ is hereafter the kernel distribution for one of the three features: G, ∆G and Lθ . A threshold, based on asymptotic of K , is used to derive a decision between alternative hypotheses: null hypothesis (H = 0) is that the two samples are drawn from the same catalyst or H = 1 means that the test rejects the null hypothesis at 5% significance level. First, Table 4 highlights that Kolmogorov-Smirnov test makes a clear discrimination between the two catalysts effective when considering parameters 4G and Lθ . Second, when performing Kolmogorov- Smirnov tests on sub-classes of catalysts Cat1 and Cat2 (2 sub-classes per catalyst due to the limited number of available samples), then the test still performs a relevant discrimination, as it can be seen in Table 6. 30 Image Anal Stereol 2018;37:21-34 Table 6. Kolmogorov-Smirnov test for the discrimination of sub-classes Catm1 and Cat n 2 of Cat1 and Cat2 respectively, m,n ∈ {1,2}, based on features computed using PSD estimated with AR-based method in polar coordinates. Null hypothesis is rejected at the level α = 0.05 if K >Ks with Ks ≈ 1.36 √ (N1 +N2)/(N1×N2): for Cat11 (N1 = 39) and Cat 2 1 (N2 = 39), Ks ≈ 0.308; for Cat11 (N1 = 39) and Cat12 (N2 = 44), Ks ≈ 0.299; for Cat11 (N1 = 39) and Cat 2 2 (N2 = 44), Ks ≈ 0.299; for Cat12 (N1 = 44) and Cat22 (N2 = 44), Ks ≈ 0.29. G 4G Lθ Stat Catmk \Catn`Cat11 Cat21 Cat12 Cat22 Cat11 Cat21 Cat12 Cat22 Cat11 Cat21 Cat12 Cat22 Cat11 0 0 0 0 0 0 1 1 0 0 1 1 H Cat21 0 0 1 0 0 0 1 1 0 0 0 0 Cat12 0 1 0 0 1 1 0 0 1 0 0 0 Cat22 0 0 0 0 1 1 0 0 1 0 0 0 Cat11 0 256 235 114 0 128 295 383 0 179 372 326 103K Cat21 256 0 366 265 128 0 306 374 179 0 206 195 Cat12 235 366 0 182 295 306 0 113 372 206 0 91 Cat22 113 265 182 0 383 374 113 0 326 195 91 0 Table 7. Kolmogorov-Smirnov test for the discrimination of sub-classes Catm1 and Cat n 2 of Cat1 and Cat2 respectively, m,n ∈ {1,2}, based on features computed using interpolated periodogram in polar coordinates. Null hypothesis is rejected at the level α = 0.05 if K > Ks with Ks ≈ 1.36 √ (N1 +N2)/(N1×N2): for Cat11 (N1 = 39) and Cat21 (N2 = 39), Ks ≈ 0.308; for Cat11 (N1 = 39) and Cat12 (N2 = 44), Ks ≈ 0.299; for Cat11 (N1 = 39) and Cat22 (N2 = 44), Ks ≈ 0.299; for Cat12 (N1 = 44) and Cat22 (N2 = 44), Ks ≈ 0.29. G 4G Lθ Stat Catmk \Catn`Cat11 Cat21 Cat12 Cat22 Cat11 Cat21 Cat12 Cat22 Cat11 Cat21 Cat12 Cat22 Cat11 0 0 0 0 0 0 0 1 0 1 0 0 H Cat21 0 0 0 1 0 0 1 1 1 0 0 1 Cat12 0 0 0 0 0 1 0 1 0 0 0 0 Cat22 0 1 0 0 1 1 1 0 0 1 0 0 Cat11 0 197 185 218 0 286 187 402 0 322 167 105 103K Cat21 197 0 199 300 286 0 314 458 322 0 224 342 Cat12 185 199 0 126 187 314 0 300 167 224 0 118 Cat22 218 300 126 0 402 458 300 0 105 342 118 0 The motivation of using the PSD estimated from an AR model (Eqs. 6, 13) is also evaluated by providing a comparison with a direct use of the PSD. Figure 7 provides the periodogram (proportional to the square of the modulus of the Discrete Fourier Transform, DFT) obtained by using zero-padding method for the same input image as in Fig. 4 – [third row]. This periodogram lacks smoothness and 2D Gaussian filters with standard deviations σ = 3,6 have been used for its regularization (Fig. 7, second and third rows). The three proposed informative parameters are then computed from this periodogram in polar coordinates via spline interpolation (Fig. 7). For the morphological analysis, different values of the threshold λ1 (Eq. 15) are used: 0.6, 0.7, 0.8 and 0.9. The best discrimination results have been obtained for σ = 3 and λ1 = 0.8. For these values, statistics of the informative parameters (Table 3), kernel distributions (Fig. 6 second row) and Kolmogorov-Smirnov tests (Tables 5, 7) are performed. The main conclusion remains the same (the two catalysts can be discriminated) but only with the distribution of the distance variation: the tangential length is no longer discriminant with the direct approach whereas the results obtained for some sub- populations show discriminations between catalysts for both the distance variation and the tangential length (Table 7). 31 TAN Z ET AL: ARFBF morphological image analysis Fig. 7. PSD estimated using the square of the modulus of the Discrete Fourier Transform. 1st row: left, periodogram with size 512× 512 obtained using zero- padding method (PSD1); right, interpolated PSD1 in polar coordinates. 2nd row: left, PSD1 smoothed by a Gaussian filter with σ = 3 (PSD2); right, interpolated PSD2 in polar coordinates. 3rd row: left, PSD1 smoothed by a Gaussian filter with σ = 6 (PSD3); right, interpolated PSD3 in polar coordinates. HRTEM observations of the two catalysts show important differences in the organization of the active phase on the support: Cat1, containing high amounts of active phase, presents as well stacks and aggregates of CoMoS slabs, in which slabs appear curved whereas Cat 2 presents isolated slabs of stacks of slabs, with a rather straight morphology. This difference is not clearly expressed by traditional measures of the length and stacking numbers. Lθ is higher in Cat 1 (the CoMoS/alumina catalyst with high Mo loading) as it can be observed in Table 2 and Fig. 6, first row. It is in coherence with the observation of aggregated slabs. Lθ translates the curvature of the slabs and can be a proper descriptor of the aggregation of the slabs. The discrimination of the catalysts using the distribution of Lθ is obtained using AR-based method for PSD estimation and not the Fourier-based method (Tables 4, 5). Thus, ARFBF morphology characterization proposed in the paper seems to be a promising feature- based description to separate the two sets of HRTEM images representing active phases of the two catalysts. CONCLUSION In this paper, we have proposed an ARFBF based morphological description of HRTEM image morphological textures corresponding to the observation of material micro-structures at nanometer scale. This morphological description is based on 2-D parametric spectrum estimation. We described the different components of the ARFBF model and we explained parameter estimation methods that lead to 2-D spectrum estimation. For the estimation of FBF Hurst parameter, we made a comparative study of different methods in order to select a robust one with respect to different sampling generators. The obtained Hurst and AR features provide a description of heterogeneous HRTEM image as a stochastic field, where the active phase fringes are associated with AR parameters and geometric fringe features have been derived by morphology analysis of the spectral fringe information. The analysis proposed allows to discriminate between two catalyst classes and opens some prospects on automatic monitoring of active phases associated with different materials, as well as material at nanoscale. It would be interesting to investigate further, the relationship between the catalytic performance and the value of Lθ parameter. ACKNOWLEDGMENT This research was supported by ARC6 grant from Rhône-Alpes French region. The authors would like to thank the reviewers (their feedback has allowed significant improvement of the paper quality) and also Anne Sophie Gay from IFPEN for her valuable help. LIST OF SYMBOLS – A: 2-D Auto-Regressive (AR) field. – QPl: prediction domains for A and associated with quarter plans, l = 1,2. – D⊂ Z2: prediction support for A. For instance, · DQP1 = {0≤ m1 ≤M1,0≤ m2 ≤M2}\(0,0); 32 Image Anal Stereol 2018;37:21-34 · DQP2 = {−M1 ≤ m1 ≤ 0,0≤ m2 ≤M2}\(0,0). – (M1,M2): order of the 2-D AR model (M1 = M2 = 10 are used for experimental results). – θ QPlM1,M2 set of parameters with the form {σ2e,l,{βm,l,m ∈ DQPl}} of the AR model. – BH : isotropic fractional Brownian field with Hurst parameter H, 0 < H < 1. – RH : auto-correlation function of BH . – Z: ARFBF, spatial deterministic convolution of field A and field BH . – S: power spectral density (PSD). – ŜB(u,v): wavelet packet spectrum in Cartesian coordinates (u,v). – Sp(r,θ): PSD in polar coordinate. – S∗(r,θ): PSD calculated from AR estimated parameters in polar coordinate. – S∗1(r,θ): PSD after removing the contribution of the first lobe on S∗. – P1, P2: frequencies providing maximum values of S∗, S∗1respectively. – α1, α2: thresholds for S∗, S∗1 respectively. – λ1 ∈ [0,1], λ2 ∈ [0,λ1]: α1 = λ1S∗(P1) and α2 = λ2S∗(P1). – S∗∗, S∗∗1 : binary images from S ∗, S∗1 when using thresholds α1, α2, respectively. – Q1, Q2: binary images obtained by using morphological reconstructions with markers P1, P2, respectively. – R: radius of disk in morphological dilation. – Te: sampling period. – G: distance value. – ∆G: distance variation (regularity of spacing). – Lθ : tangential length (regularity of curvature). – Γ: standard gamma special function. – T : Cartesian-to-polar transform. REFERENCES Akaike H (1974). A new look at the statistical model identification. IEEE T Automat Contr 19(6):716–23. Alata O, Cariou C (2008a). 2-D linear stochastic modeling. In: Garello R, Ed. Two-dimensional signal analysis, London: ISTE–Wiley. 65–114. Alata O, Cariou C (2008b). 2-D spectral analysis. In: Garello R, Ed. Two-dimensional signal analysis, London: ISTE–Wiley. 115–74. Alata O, Cariou C, Ramananjarasoa C, Najim M (1998). Classification of rotated and scaled textures using HMHV spectrum estimation and the Fourier-Mellin transform. In: Proc IEEE Int Conf Image Proc ICIP98, vol 1, pp. 53–6. Alata O, Olivier C (2003). Choice of a 2-D causal autoregressive texture model using information criteria. Pattern Recogn Lett 24:1191–201. Alata O, Ramananjarasoa C (2005). Unsupervised textured image segmentation using 2-D quarter plane autoregressive support with four prediction support. Pattern Recogn Lett 26:1069–81. Atto A, Berthoumieu Y, Bolon P (2013). 2-dimensional wavelet packet spectrum for texture analysis. IEEE T Image Process 22(6):2495–500. Atto A, Pastor D, Mercier G (2010). Wavelet packets of fractional Brownian motion: Asymptotic analysis and spectrum estimation. IEEE T Inform Theory 56:4741 – 53. Atto A, Tan Z, Alata O, Moreaud M (2014). Non-stationary texture synthesis from random field modeling. In: Proc IEEE Int Conf Image Process ICIP, pp. 4266–70. Celce B, Bres S, Moreau F, Gueroult P, Sorbier L (2008). Semi-automatic detection of sulfur slabs. In: Proc 8th Int Conf Stereol Image Anal Mater Sci STERMAT. Inż Mater 29(4):421–6. Da Costa JP, Weisbecker P, Farbos B, Leyssale J, Vignoles G, Germain C (2015). Investigating carbon materials nanostructure using image orientation statistics. Carbon 84:160–73. Fieguth PW, Willsky AS (1996). Fractal estimation using models on multiscale trees. IEEE T Signal Process 44:1297–300. Geantet C, Sorbier L (2012). Chapter 2.6 - Characterisation of catalysts. Geweke J, Porter-Hudak S (1983). The estimation and application of long memory time series models. J Time Ser Anal 4:221–37. Goncalves W, Bruno O (2013). Combining fractal and deterministic walkers for texture analysis and classification. Pattern Recogn 46:2953–68. Haralick RM (1979). Statistical and structural approaches to texture. Proc IEEE 67(5):786–804. Huang Q, Lorch JR, Dubes RC (1994). Can the fractal dimension of images be measured? Pattern Recogn 27(3):339–49. Huang Z, Li C (2005). Anisotropic fractional Brownian random fields as white noise functionals. Acta Math 33 TAN Z ET AL: ARFBF morphological image analysis Appl Sin 21(4):655–60. Jackson LB, Chien HC (1979). Frequency and bearing estimation by two-dimensional linear prediction. In: Proc IEEE Int Conf Acoust Speech Sig Proc ICASSP, pp. 665–8. Kroese DP, Botev ZI (2013). Spatial process generation. Statistics Computation. arXiv:1308.0399 Mandelbrot BB, Ness JWV (1968). Fractional Brownian motions, fractional noises and application. SIAM Review 10(4):422–37. Mao J, Jain A (1992). Texture classification and segmentation using multi-resolution simultaneous autoregressive (MR-SAR) models. Pattern Recogn 25(2):173–88. Moreaud M, Jeulin D, Morard V, Revel R (2012). TEM image analysis and modelling: application to boehmites nanoparticles. J Microsc 245(2):186–99. Moreaud M, Jeulin D, Thorel A, Chane-Ching JY (2008). A quantitative morphological analysis of nanostructured ceria–silica composite catalysts. J Microsc 232:293– 305. Nason G, Silverman B (1995). The stationary wavelet transform and some statistical applications. In: Antoniadis A, Oppenheim G, eds. Wavelets and Statistics. Lect Notes Stat 103:281–99. Nikulshin P, Ishutenko D, Mozhaev A, Maslakov K, Pimerzin A (2014). Effects of composition and morphology of active phase of CoMo/Al2O3 catalysts prepared using Co2Mo10-heteropolyacid and chelating agents on their catalytic properties in HDS and HYD reactions. J Catal 312:152–69. Parzen E (1962). On estimation of a probability density function and mode. Ann Math Statist 33:1065–76. Peacock J (1983). Two-dimensional goodness-of-fit testing in astronomy. Mon Not R Astron Soc 202(3):615–27. Pesquet-Popescu B, Larzabal P (1997). 2D self-similar processes with stationary fractional increments. In: Lévy Véhel J, Lutton E, Tricot C, eds. Fractals in Engineering. London: Springer, 138–51. Pesquet-Popescu B, Véhel JL (2002). Stochastic fractal models for image processing. IEEE Signal Proc Mag 19:48–62. Pré P, Huchet G, Jeulin D, Rouzaud J, Sennour M, Thorel A (2013). A new approach to characterize the nanostructure of activated carbons from mathematical morphology applied to high resolution transmission electron microscopy images. Carbon 52:239–58. Qazi IUH, Alata O, Burie JC, Moussa A, Fernandez- Maloigne C (2011). Choice of a pertinent color space for color texture characterization using parametric spectral analysis. Pattern Recogn Lett 44:16–31. Ranganath S, Jain AK (1985). Two-dimensional linear prediction models – part I: Spectral factorization and realization. IEEE T Acoust Speech 33(1):280–99. Richard F, Bierme H (2010). Statistical tests of anisotropy for fractional Brownian textures. Application to full- field digital mammography. J Math Imaging Vis 36(3):227 – 40. Robinson PM (1995). Log-periodogram regression of times series with long range dependence. Ann Stat 23(3):1048–72. Rosenblatt M (1956). Remarks on some nonparametric estimates of a density function. Ann Math Statist 27:832–7. Serra J (1988). Image analysis and mathematical morphology: Theoretical advances. London: Academic Press. Souza P (1982). Texture recognition via autoregression. Pattern Recogn 15(6):471–5. Tafti P, Van De Ville D, Unser M (2009). Invariances, Laplacian-like wavelet bases, and the whitening of fractal processes. IEEE T Image Process 18:689–702. Tan Z, Atto A, Alata O, Moreaud M (2015). ARFBF model for non-stationary random fields and application in HRTEM images. In: Proc IEEE Int Conf Image Process ICIP 2651–5. Toth P, Palotas A, Eddings E, Whitaker R, Lighty J (2013). A novel framework for the quantitative analysis of high resolution transmission electron micrographs of soot II. Robust multiscale nanostructure quantification. Combust Flame 160:920–32. Toulhoat H, Raybaud P (2013). Catalysis by transition metal sulphides – from molecular theory to industrial application. Paris: Technip. Van De Ville D, Blu T, Unser M (2005). Isotropic polyharmonic B-splines: Scaling functions and wavelets. IEEE T Image Process 14:1798–813. Yehliu K, Wal RV, Boehman A (2011). Development of an HRTEM image analysis method to quantify carbon nanostructure. Combust Flame 158:1837–51. Zhang BS, Yi YJ, Zhang W, Liang CH, Su DS (2011). Electron microscopy investigation of the microstructure of unsupported Ni-Mo-W sulfide. Mater Charact 62:684–90. 34