Folding a 2-D powder diffraction image into a 1-D scan: a new procedure. 6 0 0 Antonio Cervellinoa,b, Cinzia Gianninib, Antonietta Guagliardib and Massimo Ladisab,⋆ 2 n a Paul Scherrer Institute CH-5232, Villigen PSI, Switzerland a J b Istituto di Cristallografia (IC-CNR), Via G. Amendola 122/O, I-70126, Bari, Italy 5 2 (February 6, 2008) ] i c s - l r Abstract t m . t a A new procedure aiming at folding a powder diffraction 2-D into a 1-D scan m - is presented. The technique consists of three steps: tracking the beam centre d n o by means of a Simulated Annealing (SA) of the diffraction rings along the c [ same axis, detector tilt and rotation determination by a Hankel Lanczos Sin- 1 v gularValueDecomposition(HLSVD)andintensityintegration byanadaptive 6 7 binning algorithm. The X-ray powder diffraction (XRPD) intensity profile of 5 1 0 the standard NIST Si 640c sample is used to test the performances. Results 6 0 show the robustnessof the method and its capability of efficiently tagging the / t a pixels ina2-D readoutsystem by matchingtheideal geometry ofthedetector m - d to the real beam-sample-detector frame. The whole technique turns out in a n o versatile and user-friendly tool for the 2ϑ scanning of 2-D XRPD profiles. c : v ⋆ Corresponding author, E mail: [email protected], Phone: +39 0805929166, Fax: i X +39 0805929170 r a 1 I. INTRODUCTION X-ray powder diffraction (XRPD) technique is nowadays a well known tool to study crystalline properties, which provide important information for applications in many fields such as crystal structure solution, microstructure characterization, phase quantification and, recently, also in nanotechnology. In these fields, charge-coupled device (CCD) and imaging plate (IP) detectors have been widely used for collecting X-ray diffraction images, expecially because oftheir fast readout system (Hanley et al. 2005; Hammersley et al. 1997, Muchmore 1999, Phillips et al. 2002, Suzuki et al. 1999). When the analysis of the 1-D powder pattern is rather sofisticated, as for the Rietveld techniques, a reliable pre-processing aimed at enhancing the quality of XRPD data is required. When working with 2-D detectors, like a CCD or an IP system, a crucial step is the folding of the 2-D powder diffraction collected image into a 1-D scan. In the present paper a new procedure aiming at folding a powder diffraction 2-D into a 1-D scan is presented. II. THE METHOD The ideal X-ray powder diffraction geometry is cylindrical along the beam-sample axis. The diffraction cone with an opening angle 2ϑ and the apex on the sample reaches the imaging plate (a CCD detector for instance). Ideally the imaging plate would be orthogonal to the beam-sample axis and it would collect a series of circles with radii corresponding to different momentum transfer: the larger the radius the higher the momentum transfer (i.e. the resolution). Unfortunately the imaging plate axis rarely fits the ideal geometry and the circles become ellipses with an eccentricity ǫ and a rotation φ parametrizing tilt and 0 rotation of the detector frame. Furthermore the beam rarely reaches the geometrical centre of the plate and a global shift of the rings is needed to perform a correct 2ϑ scanning. The problem reduces to the nonlinear fit of the parameters appearing in the cartesian equa- tion of an ellipse: x2/a2 +y2/b2 = 1, where a/b is the major/minor semiaxis. For a generic 2 beam-sample-plate geometry the following equations hold: x = x +∆x +a cos(φ φ )/ 1+ǫ sin2(φ φ ) , c c 0 0 − − q y = y +∆y +a sin(φ φ )/ 1+ǫ sin2(φ φ ) , (1) c c 0 0 − − q being ∆x ,∆y the discrepancies along x,y axes between the beam centre and the imaging c c plate one (x ,y ). ǫ is the ellipse eccentricity a2/b2 1. The cartesian equation for the ellipse c c − finally reads: ρ2 ∆x C ∆y S ∆x2 ∆y2 1+ǫS2 1 = 2ρ c + c c + c , (2) a2 h i− (cid:20) a2 b2 (cid:21)−" a2 b2 # where ρ = (x x )2 +(y y )2, S = sin(φ φ ), C = cos(φ φ ), and φ = φ + c c 0 0 0 − − − − q arctan((y y )/(x x )). c c − − ⌋ ⌋ In Eq. (2) the l.h.s. is (ǫ) while the r.h.s. is ·§ ·† and, since usually O O s ! ⊣ ⌊ ∆x ∆y c c ǫ , Eq. (2) finally decouples. It turns out that the problem of finding the ≪ s a b rings centroid can be solved independently from the ǫ/φ determination which is indeed a 0 fine tuning with respect to the former task. Needless to say that the attempt of finding the rings centroid cannot be accomplished with by using the eq. (2) itself where both ǫ and φ 0 are unknown at this stage. The decoupling of the two problems, i.e. centroid hunting and ellipses fitting, is the main advantage of our method in comparison to the ones nowadays available in the literature (see Hammersley et al., 1996) where the full set of parameters of an equation analogous to Eq. (2) is computed. A. Beam centroid hunting This task is carried out by using a Simulated Annealing (SA) algorithm (see Metropolis et al. 1958, Pincus 1970, Kirkpatrick et al. 1982, Kirkpatrick et al. 1983, Kirkpatrick 1984) to search for a minimum of the XRPD rings centroid position on the imaging plate. As known, the algorithm employs a random search which not only accepts changes that 3 decrease objective function f, but also some changes that increase it. Changes decreasing f are accepted with a probability p = exp( δf/T), where δf is the decrease in f and T is − a control parameter, which by analogy with the original application is known as the system temperature irrespective of the objective function involved. Here the objective function f is the sharpness (defined e.g. as the inverse normalized variance) of the radial distribution function of the highest intensity rings recorded on the imaging plate and the system tem- perature is the domain size around the guess beam centre during the random search. f is computed with respect to the guess rings centre. In Fig. 1 the sharpness of the radial distri- bution function is shown before (top right) and after (bottom left) the SA. The increasing of f sharpness is dramatic and it follows the correct positioning of the beam centre onto the imaging plate (see the two spots on the top left image of Fig. 1). B. Detector tilt and rotation determination Once the beam centre has been positioned, i.e. ∆x = ∆y = 0, the eccentricity ǫ c c and rotation φ computation is straightforward; infact the eq. (2) reads 1/ρ2 1/a2 0 ≃ × 1+ǫ sin2(φ φ ) . The detector tilt and rotation are determined by using a subspace- 0 − h i based parameter estimation method called Hankel Lanczos Singular Value Decomposition (HLSVD) technique (Laudadio et al., 2002). The HLSVD method works as follows. Let us model the 1/ρ2 by samples 1/ρ2 collected at angles φ , n = 0,...N 1 as the sum of K n n − exponentially damped complex sinusoids K 1/ρ2 a exp( d φ ) cos[2πf φ +ϕ ] , (3) n ≃ k − k n k n k k=1 X a is the amplitude, ϕ the phase, d the damping factor and f the frequency of the kth k k k k sinusoid, k = 1,...,K, with K the number of damped sinusoids. This multiparametric functions family has been successfully used (Laudadio et al., 2002) for modelling NMR signals, because these functions are intrinsically related to the physical description of a NMR signal. It has been shown, however, that this family has a great approximating 4 power for many kinds of signals, including some of the natural functions used to describe crystal diffraction patterns, namely the Dirichlet kernel, related to the the Laue interference function (see Beylkin & Monz´on 2005, Sec. 5.1). Practical tests (Ladisa et al. 2005a; Ladisa et al. 2005b; Cervellino et al. 2005a; Cervellino et al. 2005b) have also confirmed their effectiveness. To estimate the parameters involved in the approximation we proceed as follows. The N data points defined in (3) are arranged (see Beylkin & Monz´on 2005, Sec. 3) into a Hankel matrix H d=ef. Hm l = 1/ρ2 , m = 0,...,M 1, l = 0,...,L 1 L×M m+l − − , with L + M = N + 1 (M L N/2). The spectral decomposition of the Hankel ≃ ≃ matrix is performed by SVD. SVD gives this method its main advantage, that is flexibility, since the number of parameters is not fixed and it is determined only by the desidered agreement level between the model and the real peak profile shape. SVD of the Hankel matrix is the decomposition H = U Σ VH , where Σ = diag λ ,λ ,...,λ , L×M L×L L×M M×M { 1 2 r} λ λ ...λ , r = min(L,M), U and V are orthogonal matrices and the superscript H 1 2 r ≥ ≥ denotes the Hermitian conjugate. A fast SVD decomposition is achieved by the Lanczos bidiagonalizationalgorithmwithpartialreorthogonalization. Thisalgorithm,basedonFFT, computes the two matrix-vector products which are performed at each step of the Lanczos procedure in (( + )log ( + )) rather than in ( ). In order to obtain the O L M ∈ L M O LM signal subspace, the matrix H is truncated to a matrix H of rank K H = U Σ VH, K K K K K where U , V , and Σ are defined by taking the first K columns of U and V, and the K K K K K upper-left matrix of Σ, respectively. As subsequent step, the least-squares solution × of the following over-determined set of equations is computed V(top)EH V(bottom), where K ≃ K (bottom) (top) V and V are derived from V by deleting its first and last row, respectively. The K K K ˆ K eigenvalues zˆ of matrix E are used to estimate the frequencies f and damping factors k k ˆ ˆ ˆ d of the model damped sinusoids from the relationship zˆ = exp d +ı2πf ∆ϑ , with k k k k − h(cid:16) (cid:17) i k = 1,...,K. Values so obtained are inserted into the model equation (3) which yields the set of equations K 1/ρ2 a exp( dˆ φ ) cos[2πfˆφ +ϕ ] , (4) n ≃ k − k n k n k k=1 X 5 with n = 0,...,N 1. The least-squares solution of (4) provides the amplitude aˆ and k − phase ϕˆ estimates of the model sinusoids which are used in the next step. k The eq. (4) and the l.h.s. of the eq. (2) formally match for K=3 and the parameters ǫ and φ finally read: ǫ = 2 a /(a a ) and φ = (ϕ π)/2. In passing we also quote the major 0 1 0 1 0 1 − − semiaxis of the ellipse a = 1/√a a . Results of HLSVD are summarized in Tab. 1; for the 0 1 − sample test the average estimated values for ǫ and φ are respectively (0.93 0.10) 10−3 0 ± × and 82◦ 5◦. ± C. Tagging pixels and radial integration We shall focus on the pixel tagging problem. We do not apply data-reduction procedures (detector response, diffraction corrections, Bragg peaks removal, etc.) at this stage since we are not interested in them. Nonetheless to say that their inclusion will not affect our conclusions (see Hammersley et al. 1996 for instance). Tagging pixels onto the 2-D detector is a crucial step in the XRPD intensity profiling prob- lem. Infact the final XRPD intensity profile follows the radial integration over the ellipses domains. Pixels are actually labeled by an adaptive binning algorithm. Adaptive binning is an analysis technique that involves subdividing a signal into samples that are more homo- geneous than the whole signal (see Sanders et al. 2001 and references therein for instance). This technique reveals information about the structure of the signal. Here we have adopted the following scheme. The XRPD 2-D image is divided into a set of annuli with common centre and eccentricity/phase computed in the previous steps. Then each annulus is tested to see whether it meets some criterion of homogeneity (e.g., whether all the pixels in the annulus are within a specific dynamic range); here we have chosen a conservative criterion: the same statistical representation within each annulus, i.e. the same number of pixels per annulus. More complicated criteria can be devised to account for other features (detector response and diffraction correction for instance). If an annulus meets the criterion, it is not divided any further. If it does not meet the criterion, it is subdivided again and the 6 test criterion is applied to the new annuli. This process is repeated iteratively until each annulus meets the criterion. The result might have annuli of several different sizes. Finally the 2ϑ scanning is achieved by integrating the 2-D XRPD intensity profile over the annular domains tagged in the previous step. The resulting intensity profile corresponding to the bottom right image of Fig. 1 is shown in Fig. 2. III. RESULTS AND CONCLUSIONS We have considered a new approach tothe 2ϑscanning in powder diffraction experiments using 2-D detectors. Our method relies on the decoupling of two main problems: the beam centre positioning and the calculation of eccentricity and rotation of the set of ellipses recorded. The two tasks are respectively accomplished with by a Simulated Annealing and a Hankel Lanczos Singular Value Decomposition, two well known and feasible techniques widely applied in several optimization problems. Tagging pixels task is achieved by an adaptive binning algorithm by which the whole image is partitioned in a set of annuli with the same statistics. Integration over the annular domains finally provides the 2ϑ intensity profile. Thisapproachhasbeensuccesfully appliedonasetofXRPDCCDimagesofsamples of the standard NIST Si 640c collected with a KappaCCD Nonius diffractometer. Acknowledgements M.L. is indebted with A. Lamura, T. Laudadio and G. Nico for deep discussions on HLSVD method. 7 REFERENCES [1] Beylkin, G. & Monz´on, L. (2005), Appl. Comput. Harmon. Anal. 19, 17-48. [2] Cervellino, A., Giannini, C., Guagliardi, A. & Ladisa, M. (2005a), Physical Review B 72, 035412. [3] Cervellino, A., Giannini, C., Guagliardi, A. & Ladisa, M. (2005b), Journal of Applied Crystallography 38, 685-687. [4] Hammersley, A. P., Svensson, S. O., HanFland, M., Fitch, A. N., and H¨ausermann, D. (1996), High Pressure Research 14, 235-248. [5] Hammersley, A. P., Brown, K., Burmeister, W., Claustre, L., Gonzalez, A., Mc- Sweeney, S., Mitchell, E., Moy, J.-P., Svensson, S. O. and Thompson, A. W. (1997), Journal of Synchrotron Radiation 4, 67-77. [6] Hanley, Q. S. and Bonner Denton, M. (2005), Journal of Synchrotron Radiation 12, 618-625. [7] Kirkpatrick, S., Gerlatt, C. D. Jr., and Vecchi, M. P. (1982), IBM Research Report RC 9355. [8] Kirkpatrick, S., Gerlatt, C. D. Jr., and Vecchi, M. P. (1983), Science 220, 671-680. [9] Kirkpatrick, S. (1984), J. Stat. Phys. 34, 975-986. [10] Ladisa, M., Lamura, A., Laudadio, T. , Nico, G. (2005a) submitted to IEEE, [arXiv:math.SP/0502539]. [11] Ladisa, M., Lamura, A., Laudadio, T. , Nico, G. (2005b) to appear on Digital Signal Processing, [arXiv:math.SP/0504171]. [12] Laudadio, T., Mastronardi, N., Vanhamme, L., Van Hecke, P. and Van Huffel, S. (2002), Journal of Magnetic Resonance, 157, 292-297. 8 [13] Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H. and Teller, E. (1958), J. Chem. Phys. 21, 1087- 1092. [14] Muchmore, S. W. (1999), Acta Cryst. D55, 1669- 1671. [15] Phillips, W. C., Stewart, A., Stanton, M., Naday, I. and Ingersoli, C. (2002), Journal of Synchrotron Radiation 9, 36-43. [16] Pincus, M. (1970), Oper. Res. 18, 1225-1228. [17] Sanders, J. S., Fabian, A. C. (2001), Mon. Not. Roy. Astron. Soc. 325 178. [18] Suzuki, M., Yamamoto, M., Kumasaka, T., Sato, K., Toyokawa, H., Aries, I. F., Jerram, P. A., Gullick, D. and Ueki, T. (1999), Journal of Synchrotron Radiation 6, 6-18. 9 FIGURES FIG. 1. Top left. Original XRPD pattern: the beam centres before and after SA are spotted on the image. Top right. Radial distribution function before SA. Bottom left. Radial distribution function after SA. Bottom right. Original XRPD pattern: the beam centre after SA together with few ellipses tagged by the HLSVD procedure are enlightened. CCD plate ( 576 x 625 pixels ) Radial distribution function at (x ,y ) = (304,278) c c 1400 1200 1000 n bi 800 s/ nt ou 600 C 400 200 0 0 100 200 300 400 500 Radius ( pixels ) CCD plate ( 576 x 625 pixels ) Radial distribution function at (x ,y ) = (313,281) c c 1400 1200 1000 n bi 800 s/ nt ou 600 C 400 200 0 0 100 200 300 400 500 Radius ( pixels ) 10