- Email: [email protected]
Circular Hough transform diffraction analysis: A software tool for automated measurement of selected area electron diffraction patterns within Digital Micrograph™
ARTICLE IN PRESS
Ultramicroscopy 108 (2008) 367–374 www.elsevier.com/locate/ultramic
Circular Hough transform diffraction analysis: A software tool for automated measurement of selected area electron diffraction patterns within Digital MicrographTM D.R.G. Mitchell Institute of Materials and Engineering Science, ANSTO, PMB 1, Menai, NSW 2234, Australia Received 2 March 2007; received in revised form 14 June 2007; accepted 20 June 2007
Abstract A software tool (script and plugin) for computing circular Hough transforms (CHT) in Digital MicrographTM has been developed, for the purpose of automated analysis of selected area electron diffraction patterns (SADPs) of polycrystalline materials. The CHT enables the diffraction pattern centre to be determined with sub-pixel accuracy, regardless of the exposure condition of the transmitted beam or if a beam stop is present. Radii of the diffraction rings can also be accurately measured with sub-pixel precision. If the pattern is calibrated against a known camera length, then d-spacings with an accuracy of better than 1% can be obtained. These measurements require no a priori knowledge of the pattern and very limited user interaction. The accuracy of the CHT is degraded by distortion introduced by the projector lens, and this should be minimised prior to pattern acquisition. A number of optimisations in the CHT software enable rapid processing of patterns; a typical analysis of a 1k 1k image taking just a few minutes. The CHT tool appears robust and is even able to accurately measure SADPs with very incomplete diffraction rings due to texture effects. This software tool is freely downloadable via the Internet. Crown Copyright r 2007 Published by Elsevier B.V. All rights reserved. PACS: 61.14.x; 07.05.Pj Keywords: Data processing; Image processing
1. Introduction The measurement of electron diffraction patterns poses a number of problems in terms of obtaining high accuracy. The centre of the pattern is not fixed but varies with projector lens deflector settings and pattern shifts also invariably occur with changes in camera length. The diffraction rings may be elliptical due to distortion introduced by the projector lens [1]. The centre is poorly defined by the transmitted spot, since this is typically much broader than the Bragg reflections and/or may be obscured by a beam stop. Finally, diffraction rings may be incomplete due to grain size and orientation effects, as well as the presence of a beam stop. Tel.: +61 2 9717 3456; fax: +61 2 9543 7179.
E-mail address: [email protected]
A number of approaches to dealing with these difficulties have been used. Gatan’s Digital MicrographTM (DM) software [2] is widely used in the electron microscopy community, and Gatan’s proprietary DIFPACKTM suite [3] enables automatic centre location. However, this uses a centre of mass approach and is therefore reliant upon the pattern being centro-symmetrical. Therefore, it may fail where intensities are unevenly distributed about the centre. The DM software has the advantage of a built-in scripting language, which permits the development of acquisition and analysis scripts by users, which can then be seamlessly integrated into the software [4]. A number of such scripts are freely available to users via the Internet [5]. One such script (Diffraction Ring Distortion-Analysis) by Hou [6] can iteratively locate the pattern centre with high accuracy by minimising the difference between a pattern and itself rotated about a guessed centre. Once the centre
0304-3991/$ - see front matter Crown Copyright r 2007 Published by Elsevier B.V. All rights reserved. doi:10.1016/j.ultramic.2007.06.003
ARTICLE IN PRESS 368
D.R.G. Mitchell / Ultramicroscopy 108 (2008) 367–374
has been located, the SADP Rotational Average script by Mitchell [7] can be used to perform a rotational average about this centre. This has the advantage of eliminating the influence of ‘diffraction ring spottiness’ due to crystal sampling effects, which would otherwise make the determination of the radial intensity profile very sensitive to the position of the profile extraction line. The resulting radial profile shows the position of the Bragg peaks. However, it is necessary for the user to measure the d-spacings manually. A fully automated approach to diffraction pattern analysis is therefore desirable. Such an approach should ideally be unaffected by typical diffraction pattern defects and characteristics mentioned earlier. It should also require minimal user interaction and no a priori knowledge of the pattern centre and ring positions. Variants of the Hough transform [8] are widely used in feature recognition applications, and have been applied to electron diffraction patterns previously [9,10]. The linear Hough transform maps features from real space (x and y Cartesian coordinates) into Hough space (polar coordinates—each point is represented as a radius and an angle). A coordinate in real space gives rise to a sinusoid in Hough space. Where points lie on a straight line in real space, they produce a series of sinusoids which overlap in Hough space. This is an accumulator space which adds votes cast by each real space point and local maxima arising from overlap of sinusoids identify straight-line features in the spatial domain [8,9,11]. Hough transforms can be extended to identify a range of geometrical shapes such as circles. However, the analysis is more complex, since in addition to the x and y coordinates of points in real space, a third term (radius) is introduced. This necessitates a three-dimensional Hough space, greatly increasing the computational demands compared to the linear transform [10,11]. At each x, y, z coordinate, circles of radius z, centred on x, y are drawn. Thus, a point in real space maps to a cone in circular Hough space. The apex of this cone occurs within the slice corresponding to a radius (z coordinate) of zero and at the x, y location corresponding to that of the pixel in real space. In practice, for each radius tested, circles of radius z are drawn centred on the coordinates corresponding to real space pixels. These circles are accumulated into a slice of circular Hough space corresponding to radius z. This is repeated for a range of radii, which should encompass the radii of the real space circular features being evaluated. When the radius being tested matches that of the circular feature being evaluated, then the circles drawn in circular Hough space will produce an intense overlap. The x and y coordinates of this overlap give the location of the centre of the circular feature being tested, and the z coordinate corresponds to the radius of the feature. Thus, it is possible to determine both the centre and radius of a circular feature with no a priori knowledge about it. However, it is useful to know (or guess) the range of possible radii present in real space, as this helps limit the number of radii which must be tested, and thus the computational overhead.
While standalone methods for carrying out circular Hough analysis have been written [9,10], such custom software may have limited user appeal due to specific hardware and software requirements and an unfamiliar user interface. The widespread use of DM amongst the TEM community provides a common software platform for microscopy users to share software tools. The scripting capability within DM provides a powerful and versatile software development environment [4]. Scripts are readily integrated into the existing DM menu functions, and users without scripting knowledge can use them in a very familiar user interface. The Digital Micrograph Script Database [5] provides a web-based focal point for users to share and exchange scripts and ideas. In this paper, the implementation of a Digital Micrograph Script to perform circular Hough transform (CHT) analysis of selected area electron diffraction patterns of polycrystalline materials is described. This script locates the pattern centre and the radii of Bragg diffraction rings with very limited user interaction required. The computationally intensive nature of the CHT has required the development of a compiled plugin containing the key CHT function and a number of software optimisations. All the software described here can be downloaded freely from the above database web site. 2. Software and hardware description Script development was carried out using DM version 3.11.0, and tested for compatibility on versions 3.9 upwards. A user-friendly dialog interface has been incorporated into the script to give the user access to all the key parameters and functions controlling the operation of the script. Compiled code was developed using the templates from Gatan’s Software Development Kit 3.6.3beta (available from the DM Script Database [5]) and Microsoft Visual Studio C++ version 6.0. The software developed here comprises a Digital Micrograph Script (CHT diffraction analysis)—referred to here as ‘CHT script’, a PlugIn (circularHoughtransform.dll), a comprehensive set of instructions and some demonstration diffraction patterns. These can be downloaded freely from the DM Script Database. Software performance was evaluated on an Intel 2 GHz Dual Core PC with 1 GB of RAM, running Windows XP. Electron diffraction patterns were recorded on a JEOL 2000FXII equipped with a Gatan DualVision camera using the maximum CCD resolution (1300 1030). 3. Software operation In principle, the CHT is capable of analysing a diffraction pattern with no user intervention at all. However, to do so, the routine would have to iterate through the complete range of circle radii which may be present in the image, i.e., from 0 up to the image diagonal. For a 1k 1k image, this would necessitate over 700 iterations, produce a data suite of several gigabytes and
ARTICLE IN PRESS D.R.G. Mitchell / Ultramicroscopy 108 (2008) 367–374
take tens of minutes to process. To avoid this, the user needs to define a range of possible radii for each diffraction ring. This is achieved in the first stage of the analysis— generating a radial intensity profile. The CHT diffraction analysis dialog is shown in Fig. 1. Radial intensity profile generation is achieved by the user selecting a SADP with ‘Get Image’ then guessing where the centre of the pattern is by positioning a central region of interest (ROI) on the diffraction pattern with the position buttons. The ROI should be large enough to guarantee that the true pattern centre is contained within it. Both the position and size of this central ROI can be changed interactively in the ‘SADP Centre’ section of the dialog (Fig. 1). This ROI is used subsequently to identity maxima in Hough space which lie within it—false maxima outside this region are rejected. The central ROI can be toggled on and off with the ‘Centre’ check box. The ‘Mask’ check box allows the user to superimpose a mask of concentric rings onto the SADP, the number of rings is set in the ‘Rings’ field. This mask moves with the central ROI and by aligning the mask rings with the diffraction rings, the user can be assured of the true pattern centre lying within the central ROI. A radial intensity profile is drawn from the centre of this ROI to the nearest edge of the pattern when the ‘Profile’ button is pressed (‘Profile’ section of the dialog (Fig. 1)). The user can set the integration width of this profile with the ‘Width’ parameter, to help reduce the impact of noise on subsequent top hat filtering (THF). Holding down the ALT key whilst pressing ‘Profile’ generates a profile which has been rotationally averaged about the guessed centre. This is useful where incomplete diffraction rings are present, and a single line profile might not intersect some of the diffraction features. Fig. 2 shows the resulting single
Fig. 1. Circular Hough transform diffraction analysis dialog. Functions are described in the text.
369
line radial intensity profile from a SADP of gold, integrated over a width of six pixels. Bragg reflections appear as intense peaks, although the indicated peak position may not represent the true position, since this is contingent upon the user guessing the centre correctly. However, the CHT script can correctly locate the centre regardless of the accuracy of the guess, provided that the true centre lies within the central ROI. The intensity profile’s purpose is to identify the range of candidate radii to test for each Bragg diffraction ring. In order to do so, peak ROIs are applied which straddle each of the peaks. These peak ROIs represent the range of radii through which the CHT script must step for each ring, to be assured of including the true radius of the ring. The peak ROI width can be set in the ‘ROI’ field in the CHT dialog (Fig. 1). Where the uncertainty in the pattern centre is high, then both the size of the central and peak ROIs should be increased to ensure that peak true positions are not missed. In a pattern with say 10 diffraction rings, and a peak ROI width of say 10 pixels, the CHT computation involves testing only 100 candidate radii which is much more manageable than processing the entire image (4700 radii). Location of the Bragg peaks has been automated by incorporating a THF into the script, which automatically processes the radial intensity profile. Top hat filtering is widely used for extracting peaks from backgrounds which may be sloped [11]. The THF parameters (‘Hat Width’, ‘Brim Width’ and ‘Threshold’) are accessible in the ‘Profile’ section of the CHT script’s dialog (Fig. 1), so that the user can tailor the THF’s response. The peak ROIs in Fig. 2 are those automatically generated by the THF. In general, well-resolved sharp Bragg peaks are readily identified. However, weak or unresolved peaks and false peaks (see later) require user intervention to add, delete, reposition and resize ROIs to deal with such features. The next stage is segmentation and thresholding of the diffraction pattern with the ‘Thresh’ button in the ‘Hough’ section of the dialog (Fig. 1). This extracts each diffraction ring for CHT analysis. The pattern is first filtered with an unsharp mask filter. The performance of this filter is controlled through the ‘USM’ parameter (kernel halfwidth), and a weighting factor (accessible within Digital Micrograph’s Global Info settings). Unsharp masking greatly enhances the visibility of the weaker diffraction rings [11]. Even after unsharp masking, thresholding the pattern as a whole would be challenging due to the large difference in intensities typically found between low and high index reflections. To overcome this, the pattern is segmented into individual diffraction rings. An annular region encompassing each diffraction ring, centred on the ‘guessed’ centre, and of the width and position defined by the corresponding peak ROI (Fig. 2) is extracted. The intensity mean and standard deviation of this extracted ring is then measured and the image greyscale thresholded so that all values below the threshold of the image’s mean intensity+standard deviation are set to zero. The ‘Thresh’ field in the dialog (Fig. 1) enables this value (the percentage
ARTICLE IN PRESS D.R.G. Mitchell / Ultramicroscopy 108 (2008) 367–374
370
of the image’s mean+standard deviation) to be adjusted to suit the pattern. During thresholding, the extracted rings are displayed sequentially. If the thresholding leads to missed rings, the user can iteratively reposition the ROIs in the intensity profile or adjust the ‘Thresh’ parameter, until all rings have been successfully thresholded. The ‘Est’ field displays an estimate of the CHT processing time on the basis of the number of pixels thresholded and the parameters selected. This allows the user to tailor the processing parameters to produce a reasonable processing time. The final step is to carry out the CHT. Initial experiments with a CHT function written in the DM Script language showed that the interpreted script language was too slow to carry out such a large computational task in a reasonable timescale. The function was therefore written in C++ and compiled as a plugin. As a result, an improvement in performance of a factor of 50 was observed, and 1k 1k patterns could be processed in several minutes. To further expedite the processing, a
resolution factor (‘Res’ in Fig. 1) was included in the function. This is a sampling frequency, where a value of 1 results in every pixel with a value 40 in the SADP being CHT processed, a value of 2 ¼ every second data point processed, 3 ¼ every third, etc. The effect of resolution on the accuracy of the CHT analysis is discussed later. Another refinement to the CHT function was the inclusion of a CHT weighting factor (‘Weighted’ in Fig. 1). Normally, the CHT function treats every pixel in the image as contributing equally to the Hough accumulator space. With the weighted option selected, the greyscale value of the pixel in real space is the value of the circle drawn in Hough space, allowing the user to weight the Hough space result according to real space intensity. CHT processing is initiated with the ‘Hough’ button. Pattern centres are determined with sub-pixel accuracy, by determining the maximum in Hough space, and then performing a centre of gravity analysis about that point over a ROI of size defined by the central ROI. Fig. 3 shows
35000 +
30000
+ +
25000
+ +
20000
+
15000
+ 10000
+
5000
+ 0
50
100
150
200
250
300
350
400
450
Fig. 2. Radial intensity profile from a SADP of gold. Bragg diffraction peaks have been automatically identified with a top hat filter, which brackets peaks with ROIs. Erroneous peaks due to noise/steeply sloping backgrounds occur (radius E20) and can be deleted. Unresolved doublets (e.g., radius E275) and weak peaks (e.g., radius E250) require manual identification.
Fig. 3. Central regions of Hough space calculated for a diffraction ring where (a) the guessed radius is 5diffraction ring radius—there is no central convergence in Hough space; (b) the guessed radius ¼ diffraction ring radius. The strong convergence identifies the guessed radius as correct and the convergence position as the centre of the diffraction ring.
ARTICLE IN PRESS D.R.G. Mitchell / Ultramicroscopy 108 (2008) 367–374
the central portion of the calculated Hough space of Synthetic Pattern 1 (see later), for two radii. Where the Hough analysis uses a radius smaller than that of the diffraction ring, there is no central convergence and a dark central region is present (Fig. 3a). Where the target radius matches that being analysed a strong central convergence occurs at the centre of the diffraction ring (Fig. 3b). The Hough space therefore furnishes both the centre of the circular feature being analysed and its radius. The raw data for each ring/radius assessed is displayed in the DM Results window. This is refined to exclude any Hough space maxima which lie outside the central ROI defined by the user. The true radius of a given diffraction ring is determined with sub-pixel precision by interpolation about the highest valid accumulated value in Hough space. The processed data are then displayed showing the radius and pattern centre for each diffraction ring. The mean values of the centre in x and y coordinates are also shown. Some plotting options are also provided in the script. The ‘Radial’ button (Fig. 1) creates a plot of Hough space intensity for all measured radii at the computed centre. The ‘Rotate’ function uses the computed pattern centre to rotationally average the diffraction pattern to produce a radial intensity profile, which provides a plot of diffracted intensity versus reciprocal spacing. 4. Assessment with simulated patterns The performance of the script was assessed by first testing it on synthetic diffraction patterns. Synthetic Pattern 1 consisted of a series of concentric rings one pixel wide (white), of various diameters, centred on x ¼ 474 and y ¼ 519, in an otherwise black image of 1024 1024 pixels. Synthetic Pattern 2 was similar to Pattern 1, except each ring was five pixels wide with a roughly gaussian intensity profile across each ring. This was a closer approximation of a real diffraction pattern. The results of CHT analysis of these patterns are shown in Table 1. The Hough resolution (‘Res’ in Fig. 1) was varied, to understand the influence of sampling frequency on the analysis. A Hough resolution of 64 results in only 1/64th of the image data being sampled. In neither pattern did resolution exert much influence on Table 1 CHT analysis of two synthetic SADPs comprising seven concentric rings of various diameters, centred on x ¼ 474, y ¼ 519 in a 1024 1024 image Hough resolution
Pattern 1 calculated centre (x:y)
Pattern 1 Pattern 2 time (s) calculated centre (x:y)
Pattern 2 time (s)
1 5 10 20 64
474.0:519.0 474.0:519.0 474.0:519.0 474.0:519.0 474.1:519.0
230 138 131 127 121
547 189 157 130 116
474.0:519.0 474.0:519.0 474.0:519.0 474.0:519.1 474.0:519.0
Ring widths: Pattern 1 ¼ 1 pixel; Pattern 2 ¼ 5 pixels with a E gaussian intensity profile.
371
the accuracy with which the centre was located. In all cases, the centre was located with a maximum error of 0.02%, and all ring radii were correctly determined to within70.1 pixels. The Hough resolution did influence the time taken to process the patterns. For Pattern 1, analysis took 230 s at Res ¼ 1, while at Res ¼ 64, which is the lowest sampling frequency permitted, the analysis time was only 121 s. Pattern 2 had longer processing times, since it contained more pixels in the wider rings. There was a diminishing advantage in using very low sampling frequencies to minimise processing time, since the script has fixed processing overheads which are independent of the pattern being analysed. Assessment of the influence of elliptical distortion of the concentric pattern on the success of the Hough transform was made (Table 2). Such distortion is introduced into experimental patterns by the projector lens [1], and Hou has written a script to assess this [6]. A simulated pattern similar to Pattern 2 was produced, but with elliptical distortion added. This distortion comprised a contraction on the y-axis and an expansion of the x-axis. The distortion indicated in Table 2 is the sum of these two distortions. Distortions of 0.8% and larger were visually apparent to the author when the Concentric Mask was overlaid, with the script’s ‘Mask’ function. The Hough space image was a more sensitive indicator of elliptical distortion. The central dark region of the Hough space image (e.g., Fig. 3a) for the largest diffraction rings appeared elliptical (to the author) rather than round when distortions as small as 0.4% were present. The centre determination remained reasonably accurate up to 0.6% distortion (Table 2). At 1% distortion, the maximum error in the centre determination reached 0.3%. Errors in the accuracy of individual ring radii determination became apparent at distortions of 0.4% and higher. At a distortion of 0.4%, the calculated radius of four of the seven rings was in error, albeit by a cumulative error of less than four pixels. The centre error—the standard deviation of the average centre position determined from all the rings, generally increased with distortion. It is clear that Table 2 Influence of elliptical distortion on CHT analysis of Synthetic Pattern 2 comprising seven concentric rings of various diameters, centred on x ¼ 474, y ¼ 519 in a 1024 1024 image Distortion Calculated (%) centre (x:y)
Maximum centre error (%)
No. of ring Cumulative ring error (pixels)b errorsa
0 0.1 0.2 0.4 0.6 0.8 1.0
0.00 0.00 0.00 0.07 0.08 0.52 0.34
0 0 0 4 5 7 7
a
474.0:519.0 474.0:519.0 474.0:519.0 474.2:518.6 474.4:519.2 476.5:519.1 475.6:519.4
0 0 0 3.5 9.7 13.5 11.3
No. of rings where the calculated radius is in error by 470.1 pixels. absolute value (calculated radiusactual radius).
bP
ARTICLE IN PRESS 372
D.R.G. Mitchell / Ultramicroscopy 108 (2008) 367–374
elliptical distortion has the potential to degrade the accuracy of the CHT analysis, and should therefore be minimised prior to pattern acquisition. 5. Assessment with experimental patterns Fig. 4a shows a SADP of a polycrystalline gold film. This represents the best case pattern for this type of analysis, since the rings are sharp and continuous and there is no beam stop present. The pattern was calibrated using the camera length previously determined for our instrument. It was then analysed using the CHT script. Fig. 2 shows the radial intensity profile from this pattern, and the ROIs fitted by the automatic THF routine. The THF may produce false peak identifications due to noise or steeply sloping backgrounds. The ‘peak’ at a radius of E20 pixels is an example of this, and this was deleted prior to CHT analysis. The THF may not locate all the peaks which are present, such as the weak peaks at radii of E250 pixels and E440 pixels and the unresolved doublets at radii of E210 and E275. Some refinement of the THF parameters (‘Brim’, ‘Hat’ and ‘Level’ (Fig. 1)) may help improve THF performance. The peaks missed by the THF were manually located. It is essential that the automatic thresholding routine invoked by the ‘Thresh’ button, can threshold the peak from the background. Very weak or noisy peaks which are difficult to threshold can be omitted from the CHT analysis. Following accurate determination of the pattern centre using the intense peaks, the ‘Rotate’ function can be used to produce the rotationally averaged radial intensity profile. The position of the omitted weak peaks can then be determined manually. A total of 14 Bragg peaks were measured for the gold pattern (Fig. 4) (Note: six additional peaks not located by the THF in Fig. 2 were manually located). Table 3 shows results of a series of independent measurements on the gold pattern at various Hough resolutions (‘Res’ in Fig. 1). The accuracy of the CHT analysis was judged by plotting the
measured d-spacings of the gold versus the actual d-spacings, and measuring the slope (‘Slope’ in Table 3). The quality of the fit was assessed by the R-value of the linear regression (R in Table 3). Clearly, in the first instance, the accuracy of any such d-spacing measurement on a TEM is dependent upon the accuracy of the original camera length calibration. However, any errors in that calibration are common to all the data determined here. At the highest resolution (Res ¼ 1), which took 500 s to complete, the slope of the plot of measured d-spacing versus actual d-spacing was 1.0079. At a resolution of 64 (the lowest resolution), the analysis time was 73 s, and the measured slope was 1.0062. R-values at all resolutions were also very close to 1 in all cases. The measured centres of the pattern at all resolutions were also very consistent, with the largest standard deviation of the average centre position corresponding to just 0.2%. Measurement errors affect the d-spacing values of the innermost diffraction rings (smallest radii) the most. At all resolutions, the error in the reported d-spacing for the innermost diffraction ring, was 0.7%. This error includes any errors arising from the camera length calibration, as well those from the CHT analysis. Fig. 4b shows a worst case experimental diffraction pattern for CHT analysis. It is from a heavily textured TiO2 Table 3 Influence of Hough resolution on the accuracy of d-spacing determination from CHT analysis of a SADP of polycrystalline gold Resolution
Time (s)
Centre (x:y)
Max centre error (%)
Slopea
Rb
1 5 10 32 64
500 180 100 82 73
719.2:480.7 719.2:481.0 719.2:480.9 719.4:480.7 719.2:480.9
0.08 0.07 0.08 0.08 0.14
1.0079 1.0079 1.0084 1.0054 1.0062
0.99999 0.99999 0.99999 0.99999 0.99998
a
Slope of plot of CHT measured d-spacings versus actual d-spacings. Linear regression R-value of plot in footnote a.
b
Fig. 4. The CHT analysis superimposes a coloured mask on the SADP, showing the position and size of the rings: (a) polycrystalline gold: the SADP (left) and the mask (right) are shown separately for clarity; (b) textured TiO2 (anatase) showing that despite incomplete rings the CHT analysis has located the majority of Bragg peaks.
ARTICLE IN PRESS D.R.G. Mitchell / Ultramicroscopy 108 (2008) 367–374
(anatase) film [12], with two preferred orientational variants. The Bragg rings are incomplete and a beam stop outline runs through the pattern. Due to the incomplete diffraction rings it was necessary to use the rotational averaging facility during profile generation, to locate diffraction features which would not otherwise be intersected by a simple linear profile. The CHT analysis correctly located most of the principal Bragg peaks. The standard deviation on the average centre position (0.35%) was somewhat higher than that of the gold pattern under similar conditions (0.1%). The error in the d-spacing determined for the innermost diffraction ring (most sensitive to measurement error) was 1.6%, which compared with 0.7% for the innermost diffraction ring in the gold pattern. 6. Discussion Accurate measurement of electron diffraction patterns poses a number of problems. The uncertainty in the location of the pattern centre precludes accurate measurement of ring radii [9]. Centre of mass approaches to determining the pattern centre [3] are critically dependent upon a centro-symmetrical distribution of intensities in the patterns, and this symmetry may be lost due to the presence of a beam stop or an off-axis condition. Even with a known centre, other pattern effects, such as ring spottiness due to grain size/sampling effects and preferred orientation makes any radial intensity profile very sensitive to the positioning of the extraction line. Rotationally averaging a diffraction pattern about a centre can overcome these two effects [7]. However, this approach is critically dependent upon accurately knowing the centre, and errors in the pattern centre will lead to broadening of peaks in the intensity profile. Pattern centres can be crudely estimated via superposition of a circular feature onto a diffraction ring or by taking the average intercept of g vectors. However, the likely uncertainty of these approaches is probably of the order 75 pixels for a typical pattern captured at 1k 1k resolution. For a high index diffraction ring of radius E100 pixels, this leads to an error of 5%, in the d-spacing determination, which for many measurements is unsatisfactory. Distortion induced by the projector lens [1] can also compound these errors, due to the direction dependence of the measured d-spacing. Pattern measurement can also be time consuming, and so a readily implemented method of performing these measurements accurately and with limited user involvement is desirable. The approach of Hou [6] goes some way towards addressing this need. The approach has been to determine the pattern centre by iterating from a guessed starting position, to minimise the difference between the original image, and one rotated about the centre. This provides not only the pattern centre with high accuracy and precision, but also a measure of the pattern ellipsicity induced by the projector lens. Comparison of the centre determination of
373
the gold pattern (Fig. 4) using the CHT approach and that of Hou’s showed consistent sub-pixel agreement. The centre of mass approach was wildly inaccurate on this pattern. While the gold pattern was centro-symmetric, the weak ‘comet tail’ caused by the electronic shuttering of the CCD camera deflecting the intense transmitted beam off the CCD, was sufficient to shift the centre of mass away from the true pattern centre. As well as providing the pattern centre, the CHT automatically locates the diffraction rings and measures them with sub-pixel precision. Provided the pattern is calibrated beforehand, this yields a listing of d-spacings. Comparison of the measured d-spacings for gold with those expected showed excellent agreement (Table 3). Errors in pattern measurement are most acute for low index diffraction rings of small radius. In the case of the gold measurements in Table 3, the difference between the measured d-spacing and the actual d-spacing of the innermost diffraction ring was just 0.7%. Some of this error is due to the error in the camera length calibration. The radius of this ring was 106 pixels, and referring to the earlier discussion on the likely accuracy of crude estimates of the centre position, an uncertainty of 75 pixels in the centre position, corresponds to a potential error in the d-spacing of up to 5%. The results of the CHT are clearly quite superior. The CHT is also robust when dealing with patterns containing incomplete rings. Russ et al. [10] have shown that textured or single crystal patterns from asbestos polymorphs can be analysed using this method. The highly textured anatase pattern (Fig. 4b) was successfully analysed using the CHT method. The uncertainty in the location of the pattern centre was higher (0.35%) than that of the gold (0.1%), and the error in the d-spacing of the innermost diffraction ring (1.6%) was also higher than that of gold (0.7%). A perfectly circular feature in real space, will produce an intensity distribution in circular Hough space which is radially symmetrical, and contains a sharply defined maximum—the position of which defines the circular feature’s centre. Fig. 3b shows a symmetrical Hough space derived from such a circular feature. Circular Hough space is an accumulator space derived from overlapping circles, and the more pixels that are present on any given radius, the greater the peak value in the corresponding Hough space. If one considers two touching circles, the centres of which are at the 12 and 6 o’clock positions on a circular feature, the position of the overlap of these two circles (the clock face centre) is very sharply defined along the vector joining 12–6 o’clock. However, it is poorly defined at right angles to that vector (9–3 o’clock). If additional circles are added at the 3 and 9 o’clock positions, the centre becomes much more sharply defined, and improves with additional circles at other clock positions. The textured anatase pattern (Fig. 4b) consists of limited numbers of spot pairs which results in a lower accuracy of centre and radius determination compared with the polycrystalline gold ring pattern (Fig. 4b) for the
ARTICLE IN PRESS 374
D.R.G. Mitchell / Ultramicroscopy 108 (2008) 367–374
above reasons. Nevertheless, the performance of the CHT on the textured anatase pattern was satisfactory. The CHT is a computationally intensive method, and its implementation as a useable script for DM has necessitated a number of steps to ensure processing time is not excessive. The creation of a compiled plugin to handle the Hough processing resulted in a 50-fold increase in processing speed. The use of an intensity profile in which the user selects the radii of interest, rather than sampling the entire range of radii which may be present, results in a further 5–10-fold decrease in processing time. The inclusion of a resolution (‘Res’) parameter to allow sub-sampling of real space allows the user to tailor the processing time to keep analysis times reasonable. Significant sub-sampling can be used without accuracy adversely suffering (Tables 1 and 3). The default value (Res ¼ 5) results in a E2–3-fold reduction in processing time, compared with the maximum resolution. The result of these optimisations is that a typical 1k 1k diffraction pattern can be analysed with a processing time of just a few minutes, which certainly compares very favourably with fully manual approaches to measurement. The compiled plugin makes the CHT function available to the DM scripting language, and examples of its use are included in the instructions which can be downloaded with the package. This permits users to incorporate this functionality into their own scripts. Potential applications for this include automatic location/measurement of any circular features, such a convergent beam electron diffraction pattern discs, circular holes/pores/particles, etc. 7. Conclusions The circular Hough transform (CHT) has been successfully incorporated into a user-friendly software tool for automatic analysis of selected area electron patterns within Digital MicrographTM. This software tool can be freely downloaded from the Internet. The computationally intensive CHT method has necessitated the development of a compiled plugin to handle the key Hough transform processing, and a number of software optimisations to permit analysis of diffraction patterns in just a few minutes. The accuracy and robustness of the software has been tested with simulated diffraction patterns. The influence of projector lens distortion on the CHT analysis has been assessed, and the need to minimise any such distortion highlighted. Assessment of best case (polycrystalline gold) and worst case (textured anatase) experimental selected
area diffraction patterns, showed that the CHT software tool was robust and accurate, determining the pattern centres with a maximum uncertainty of 0.1% and 0.35%, respectively. Errors in the measurement of the d-spacing of the innermost diffraction ring (most sensitive to measurement errors), were 0.7% and 1.6%, respectively, although this error did include the error from the camera length calibration of the microscope. The incorporation a CHT function into the Digital MicrographTM scripting language will enable users to develop their own scripts for location/ measurement of circular features. Acknowledgements Dr. Bernhard Schaffer (Graz University of Technology, Austria) is thanked for provision of a C++ function to sort images into a list. This function was incorporated into the CircularHoughTransform.dll described here. Christoph Koch and Douglas Hague are thanked for contributing useful demonstration code of C++ and its compilation to the Digital Micrograph Script Database [5]. Dr. Bill Bertram (ANSTO) is thanked for mathematical assistance with algorithm development. Darren Attard (ANSTO) is thanked for testing this script. References [1] G.C. Capitani, P. Oleynikov, S. Hovmo¨ller, M. Mellini, Ultramicroscopy 106 (2006) 66. [2] Digital MicrographTM, /http://www.gatan.com/imaging/downloads. phpS, Gatan Inc., Pleasanton, CA, 2007. [3] DIFPACKTM, Diffraction analysis package for Digital MicrographTM, Gatan Inc., Pleasanton, CA, 2007. [4] D.R.G. Mitchell, B. Schaffer, Ultramicroscopy 103 (2005) 319. [5] Digital MicrographTM Script Database: /http://www.felmi-zfe. tugraz.at/dm_scripts/dmscript1.htmlS, 2007. [6] V. Hou, Diffraction Rings Distortion-Analysis, Digital Micrograph Script in Digital MicrographTM Script Database: /http://www.felmi-zfe. tugraz.at/dm_scripts/dmscript1.htmlS, 2006. [7] D.R.G. Mitchell, SADP Rotational Average, Digital Micrograph Script in Digital MicrographTM Script Database: /http://www.felmi-zfe. tugraz.at/dm_scripts/dmscript1.htmlS, 2004. [8] P.V.C. Hough, Method and means for recognizing complex patterns, US Patent, 3069654, 1962. [9] J.C. Russ, D.S. Bright, J.C. Russ, M. Hare, J. Comput. Assist. Microsc. 1 (1989) 3. [10] J.C. Russ, T. Taguchi, P.M. Peters, E. Chatfield, J.C. Russ, W.D. Stewart, Adv. X-Ray Anal. 32 (1989) 593. [11] J.C. Russ, The Image Processing Handbook, CRC Press (1996). [12] D.R.G. Mitchell, D.J. Attard, G. Triani, J. Cryst. Growth 285 (2005) 208.
Ultramicroscopy 108 (2008) 367–374 www.elsevier.com/locate/ultramic
Circular Hough transform diffraction analysis: A software tool for automated measurement of selected area electron diffraction patterns within Digital MicrographTM D.R.G. Mitchell Institute of Materials and Engineering Science, ANSTO, PMB 1, Menai, NSW 2234, Australia Received 2 March 2007; received in revised form 14 June 2007; accepted 20 June 2007
Abstract A software tool (script and plugin) for computing circular Hough transforms (CHT) in Digital MicrographTM has been developed, for the purpose of automated analysis of selected area electron diffraction patterns (SADPs) of polycrystalline materials. The CHT enables the diffraction pattern centre to be determined with sub-pixel accuracy, regardless of the exposure condition of the transmitted beam or if a beam stop is present. Radii of the diffraction rings can also be accurately measured with sub-pixel precision. If the pattern is calibrated against a known camera length, then d-spacings with an accuracy of better than 1% can be obtained. These measurements require no a priori knowledge of the pattern and very limited user interaction. The accuracy of the CHT is degraded by distortion introduced by the projector lens, and this should be minimised prior to pattern acquisition. A number of optimisations in the CHT software enable rapid processing of patterns; a typical analysis of a 1k 1k image taking just a few minutes. The CHT tool appears robust and is even able to accurately measure SADPs with very incomplete diffraction rings due to texture effects. This software tool is freely downloadable via the Internet. Crown Copyright r 2007 Published by Elsevier B.V. All rights reserved. PACS: 61.14.x; 07.05.Pj Keywords: Data processing; Image processing
1. Introduction The measurement of electron diffraction patterns poses a number of problems in terms of obtaining high accuracy. The centre of the pattern is not fixed but varies with projector lens deflector settings and pattern shifts also invariably occur with changes in camera length. The diffraction rings may be elliptical due to distortion introduced by the projector lens [1]. The centre is poorly defined by the transmitted spot, since this is typically much broader than the Bragg reflections and/or may be obscured by a beam stop. Finally, diffraction rings may be incomplete due to grain size and orientation effects, as well as the presence of a beam stop. Tel.: +61 2 9717 3456; fax: +61 2 9543 7179.
E-mail address: [email protected]
A number of approaches to dealing with these difficulties have been used. Gatan’s Digital MicrographTM (DM) software [2] is widely used in the electron microscopy community, and Gatan’s proprietary DIFPACKTM suite [3] enables automatic centre location. However, this uses a centre of mass approach and is therefore reliant upon the pattern being centro-symmetrical. Therefore, it may fail where intensities are unevenly distributed about the centre. The DM software has the advantage of a built-in scripting language, which permits the development of acquisition and analysis scripts by users, which can then be seamlessly integrated into the software [4]. A number of such scripts are freely available to users via the Internet [5]. One such script (Diffraction Ring Distortion-Analysis) by Hou [6] can iteratively locate the pattern centre with high accuracy by minimising the difference between a pattern and itself rotated about a guessed centre. Once the centre
0304-3991/$ - see front matter Crown Copyright r 2007 Published by Elsevier B.V. All rights reserved. doi:10.1016/j.ultramic.2007.06.003
ARTICLE IN PRESS 368
D.R.G. Mitchell / Ultramicroscopy 108 (2008) 367–374
has been located, the SADP Rotational Average script by Mitchell [7] can be used to perform a rotational average about this centre. This has the advantage of eliminating the influence of ‘diffraction ring spottiness’ due to crystal sampling effects, which would otherwise make the determination of the radial intensity profile very sensitive to the position of the profile extraction line. The resulting radial profile shows the position of the Bragg peaks. However, it is necessary for the user to measure the d-spacings manually. A fully automated approach to diffraction pattern analysis is therefore desirable. Such an approach should ideally be unaffected by typical diffraction pattern defects and characteristics mentioned earlier. It should also require minimal user interaction and no a priori knowledge of the pattern centre and ring positions. Variants of the Hough transform [8] are widely used in feature recognition applications, and have been applied to electron diffraction patterns previously [9,10]. The linear Hough transform maps features from real space (x and y Cartesian coordinates) into Hough space (polar coordinates—each point is represented as a radius and an angle). A coordinate in real space gives rise to a sinusoid in Hough space. Where points lie on a straight line in real space, they produce a series of sinusoids which overlap in Hough space. This is an accumulator space which adds votes cast by each real space point and local maxima arising from overlap of sinusoids identify straight-line features in the spatial domain [8,9,11]. Hough transforms can be extended to identify a range of geometrical shapes such as circles. However, the analysis is more complex, since in addition to the x and y coordinates of points in real space, a third term (radius) is introduced. This necessitates a three-dimensional Hough space, greatly increasing the computational demands compared to the linear transform [10,11]. At each x, y, z coordinate, circles of radius z, centred on x, y are drawn. Thus, a point in real space maps to a cone in circular Hough space. The apex of this cone occurs within the slice corresponding to a radius (z coordinate) of zero and at the x, y location corresponding to that of the pixel in real space. In practice, for each radius tested, circles of radius z are drawn centred on the coordinates corresponding to real space pixels. These circles are accumulated into a slice of circular Hough space corresponding to radius z. This is repeated for a range of radii, which should encompass the radii of the real space circular features being evaluated. When the radius being tested matches that of the circular feature being evaluated, then the circles drawn in circular Hough space will produce an intense overlap. The x and y coordinates of this overlap give the location of the centre of the circular feature being tested, and the z coordinate corresponds to the radius of the feature. Thus, it is possible to determine both the centre and radius of a circular feature with no a priori knowledge about it. However, it is useful to know (or guess) the range of possible radii present in real space, as this helps limit the number of radii which must be tested, and thus the computational overhead.
While standalone methods for carrying out circular Hough analysis have been written [9,10], such custom software may have limited user appeal due to specific hardware and software requirements and an unfamiliar user interface. The widespread use of DM amongst the TEM community provides a common software platform for microscopy users to share software tools. The scripting capability within DM provides a powerful and versatile software development environment [4]. Scripts are readily integrated into the existing DM menu functions, and users without scripting knowledge can use them in a very familiar user interface. The Digital Micrograph Script Database [5] provides a web-based focal point for users to share and exchange scripts and ideas. In this paper, the implementation of a Digital Micrograph Script to perform circular Hough transform (CHT) analysis of selected area electron diffraction patterns of polycrystalline materials is described. This script locates the pattern centre and the radii of Bragg diffraction rings with very limited user interaction required. The computationally intensive nature of the CHT has required the development of a compiled plugin containing the key CHT function and a number of software optimisations. All the software described here can be downloaded freely from the above database web site. 2. Software and hardware description Script development was carried out using DM version 3.11.0, and tested for compatibility on versions 3.9 upwards. A user-friendly dialog interface has been incorporated into the script to give the user access to all the key parameters and functions controlling the operation of the script. Compiled code was developed using the templates from Gatan’s Software Development Kit 3.6.3beta (available from the DM Script Database [5]) and Microsoft Visual Studio C++ version 6.0. The software developed here comprises a Digital Micrograph Script (CHT diffraction analysis)—referred to here as ‘CHT script’, a PlugIn (circularHoughtransform.dll), a comprehensive set of instructions and some demonstration diffraction patterns. These can be downloaded freely from the DM Script Database. Software performance was evaluated on an Intel 2 GHz Dual Core PC with 1 GB of RAM, running Windows XP. Electron diffraction patterns were recorded on a JEOL 2000FXII equipped with a Gatan DualVision camera using the maximum CCD resolution (1300 1030). 3. Software operation In principle, the CHT is capable of analysing a diffraction pattern with no user intervention at all. However, to do so, the routine would have to iterate through the complete range of circle radii which may be present in the image, i.e., from 0 up to the image diagonal. For a 1k 1k image, this would necessitate over 700 iterations, produce a data suite of several gigabytes and
ARTICLE IN PRESS D.R.G. Mitchell / Ultramicroscopy 108 (2008) 367–374
take tens of minutes to process. To avoid this, the user needs to define a range of possible radii for each diffraction ring. This is achieved in the first stage of the analysis— generating a radial intensity profile. The CHT diffraction analysis dialog is shown in Fig. 1. Radial intensity profile generation is achieved by the user selecting a SADP with ‘Get Image’ then guessing where the centre of the pattern is by positioning a central region of interest (ROI) on the diffraction pattern with the position buttons. The ROI should be large enough to guarantee that the true pattern centre is contained within it. Both the position and size of this central ROI can be changed interactively in the ‘SADP Centre’ section of the dialog (Fig. 1). This ROI is used subsequently to identity maxima in Hough space which lie within it—false maxima outside this region are rejected. The central ROI can be toggled on and off with the ‘Centre’ check box. The ‘Mask’ check box allows the user to superimpose a mask of concentric rings onto the SADP, the number of rings is set in the ‘Rings’ field. This mask moves with the central ROI and by aligning the mask rings with the diffraction rings, the user can be assured of the true pattern centre lying within the central ROI. A radial intensity profile is drawn from the centre of this ROI to the nearest edge of the pattern when the ‘Profile’ button is pressed (‘Profile’ section of the dialog (Fig. 1)). The user can set the integration width of this profile with the ‘Width’ parameter, to help reduce the impact of noise on subsequent top hat filtering (THF). Holding down the ALT key whilst pressing ‘Profile’ generates a profile which has been rotationally averaged about the guessed centre. This is useful where incomplete diffraction rings are present, and a single line profile might not intersect some of the diffraction features. Fig. 2 shows the resulting single
Fig. 1. Circular Hough transform diffraction analysis dialog. Functions are described in the text.
369
line radial intensity profile from a SADP of gold, integrated over a width of six pixels. Bragg reflections appear as intense peaks, although the indicated peak position may not represent the true position, since this is contingent upon the user guessing the centre correctly. However, the CHT script can correctly locate the centre regardless of the accuracy of the guess, provided that the true centre lies within the central ROI. The intensity profile’s purpose is to identify the range of candidate radii to test for each Bragg diffraction ring. In order to do so, peak ROIs are applied which straddle each of the peaks. These peak ROIs represent the range of radii through which the CHT script must step for each ring, to be assured of including the true radius of the ring. The peak ROI width can be set in the ‘ROI’ field in the CHT dialog (Fig. 1). Where the uncertainty in the pattern centre is high, then both the size of the central and peak ROIs should be increased to ensure that peak true positions are not missed. In a pattern with say 10 diffraction rings, and a peak ROI width of say 10 pixels, the CHT computation involves testing only 100 candidate radii which is much more manageable than processing the entire image (4700 radii). Location of the Bragg peaks has been automated by incorporating a THF into the script, which automatically processes the radial intensity profile. Top hat filtering is widely used for extracting peaks from backgrounds which may be sloped [11]. The THF parameters (‘Hat Width’, ‘Brim Width’ and ‘Threshold’) are accessible in the ‘Profile’ section of the CHT script’s dialog (Fig. 1), so that the user can tailor the THF’s response. The peak ROIs in Fig. 2 are those automatically generated by the THF. In general, well-resolved sharp Bragg peaks are readily identified. However, weak or unresolved peaks and false peaks (see later) require user intervention to add, delete, reposition and resize ROIs to deal with such features. The next stage is segmentation and thresholding of the diffraction pattern with the ‘Thresh’ button in the ‘Hough’ section of the dialog (Fig. 1). This extracts each diffraction ring for CHT analysis. The pattern is first filtered with an unsharp mask filter. The performance of this filter is controlled through the ‘USM’ parameter (kernel halfwidth), and a weighting factor (accessible within Digital Micrograph’s Global Info settings). Unsharp masking greatly enhances the visibility of the weaker diffraction rings [11]. Even after unsharp masking, thresholding the pattern as a whole would be challenging due to the large difference in intensities typically found between low and high index reflections. To overcome this, the pattern is segmented into individual diffraction rings. An annular region encompassing each diffraction ring, centred on the ‘guessed’ centre, and of the width and position defined by the corresponding peak ROI (Fig. 2) is extracted. The intensity mean and standard deviation of this extracted ring is then measured and the image greyscale thresholded so that all values below the threshold of the image’s mean intensity+standard deviation are set to zero. The ‘Thresh’ field in the dialog (Fig. 1) enables this value (the percentage
ARTICLE IN PRESS D.R.G. Mitchell / Ultramicroscopy 108 (2008) 367–374
370
of the image’s mean+standard deviation) to be adjusted to suit the pattern. During thresholding, the extracted rings are displayed sequentially. If the thresholding leads to missed rings, the user can iteratively reposition the ROIs in the intensity profile or adjust the ‘Thresh’ parameter, until all rings have been successfully thresholded. The ‘Est’ field displays an estimate of the CHT processing time on the basis of the number of pixels thresholded and the parameters selected. This allows the user to tailor the processing parameters to produce a reasonable processing time. The final step is to carry out the CHT. Initial experiments with a CHT function written in the DM Script language showed that the interpreted script language was too slow to carry out such a large computational task in a reasonable timescale. The function was therefore written in C++ and compiled as a plugin. As a result, an improvement in performance of a factor of 50 was observed, and 1k 1k patterns could be processed in several minutes. To further expedite the processing, a
resolution factor (‘Res’ in Fig. 1) was included in the function. This is a sampling frequency, where a value of 1 results in every pixel with a value 40 in the SADP being CHT processed, a value of 2 ¼ every second data point processed, 3 ¼ every third, etc. The effect of resolution on the accuracy of the CHT analysis is discussed later. Another refinement to the CHT function was the inclusion of a CHT weighting factor (‘Weighted’ in Fig. 1). Normally, the CHT function treats every pixel in the image as contributing equally to the Hough accumulator space. With the weighted option selected, the greyscale value of the pixel in real space is the value of the circle drawn in Hough space, allowing the user to weight the Hough space result according to real space intensity. CHT processing is initiated with the ‘Hough’ button. Pattern centres are determined with sub-pixel accuracy, by determining the maximum in Hough space, and then performing a centre of gravity analysis about that point over a ROI of size defined by the central ROI. Fig. 3 shows
35000 +
30000
+ +
25000
+ +
20000
+
15000
+ 10000
+
5000
+ 0
50
100
150
200
250
300
350
400
450
Fig. 2. Radial intensity profile from a SADP of gold. Bragg diffraction peaks have been automatically identified with a top hat filter, which brackets peaks with ROIs. Erroneous peaks due to noise/steeply sloping backgrounds occur (radius E20) and can be deleted. Unresolved doublets (e.g., radius E275) and weak peaks (e.g., radius E250) require manual identification.
Fig. 3. Central regions of Hough space calculated for a diffraction ring where (a) the guessed radius is 5diffraction ring radius—there is no central convergence in Hough space; (b) the guessed radius ¼ diffraction ring radius. The strong convergence identifies the guessed radius as correct and the convergence position as the centre of the diffraction ring.
ARTICLE IN PRESS D.R.G. Mitchell / Ultramicroscopy 108 (2008) 367–374
the central portion of the calculated Hough space of Synthetic Pattern 1 (see later), for two radii. Where the Hough analysis uses a radius smaller than that of the diffraction ring, there is no central convergence and a dark central region is present (Fig. 3a). Where the target radius matches that being analysed a strong central convergence occurs at the centre of the diffraction ring (Fig. 3b). The Hough space therefore furnishes both the centre of the circular feature being analysed and its radius. The raw data for each ring/radius assessed is displayed in the DM Results window. This is refined to exclude any Hough space maxima which lie outside the central ROI defined by the user. The true radius of a given diffraction ring is determined with sub-pixel precision by interpolation about the highest valid accumulated value in Hough space. The processed data are then displayed showing the radius and pattern centre for each diffraction ring. The mean values of the centre in x and y coordinates are also shown. Some plotting options are also provided in the script. The ‘Radial’ button (Fig. 1) creates a plot of Hough space intensity for all measured radii at the computed centre. The ‘Rotate’ function uses the computed pattern centre to rotationally average the diffraction pattern to produce a radial intensity profile, which provides a plot of diffracted intensity versus reciprocal spacing. 4. Assessment with simulated patterns The performance of the script was assessed by first testing it on synthetic diffraction patterns. Synthetic Pattern 1 consisted of a series of concentric rings one pixel wide (white), of various diameters, centred on x ¼ 474 and y ¼ 519, in an otherwise black image of 1024 1024 pixels. Synthetic Pattern 2 was similar to Pattern 1, except each ring was five pixels wide with a roughly gaussian intensity profile across each ring. This was a closer approximation of a real diffraction pattern. The results of CHT analysis of these patterns are shown in Table 1. The Hough resolution (‘Res’ in Fig. 1) was varied, to understand the influence of sampling frequency on the analysis. A Hough resolution of 64 results in only 1/64th of the image data being sampled. In neither pattern did resolution exert much influence on Table 1 CHT analysis of two synthetic SADPs comprising seven concentric rings of various diameters, centred on x ¼ 474, y ¼ 519 in a 1024 1024 image Hough resolution
Pattern 1 calculated centre (x:y)
Pattern 1 Pattern 2 time (s) calculated centre (x:y)
Pattern 2 time (s)
1 5 10 20 64
474.0:519.0 474.0:519.0 474.0:519.0 474.0:519.0 474.1:519.0
230 138 131 127 121
547 189 157 130 116
474.0:519.0 474.0:519.0 474.0:519.0 474.0:519.1 474.0:519.0
Ring widths: Pattern 1 ¼ 1 pixel; Pattern 2 ¼ 5 pixels with a E gaussian intensity profile.
371
the accuracy with which the centre was located. In all cases, the centre was located with a maximum error of 0.02%, and all ring radii were correctly determined to within70.1 pixels. The Hough resolution did influence the time taken to process the patterns. For Pattern 1, analysis took 230 s at Res ¼ 1, while at Res ¼ 64, which is the lowest sampling frequency permitted, the analysis time was only 121 s. Pattern 2 had longer processing times, since it contained more pixels in the wider rings. There was a diminishing advantage in using very low sampling frequencies to minimise processing time, since the script has fixed processing overheads which are independent of the pattern being analysed. Assessment of the influence of elliptical distortion of the concentric pattern on the success of the Hough transform was made (Table 2). Such distortion is introduced into experimental patterns by the projector lens [1], and Hou has written a script to assess this [6]. A simulated pattern similar to Pattern 2 was produced, but with elliptical distortion added. This distortion comprised a contraction on the y-axis and an expansion of the x-axis. The distortion indicated in Table 2 is the sum of these two distortions. Distortions of 0.8% and larger were visually apparent to the author when the Concentric Mask was overlaid, with the script’s ‘Mask’ function. The Hough space image was a more sensitive indicator of elliptical distortion. The central dark region of the Hough space image (e.g., Fig. 3a) for the largest diffraction rings appeared elliptical (to the author) rather than round when distortions as small as 0.4% were present. The centre determination remained reasonably accurate up to 0.6% distortion (Table 2). At 1% distortion, the maximum error in the centre determination reached 0.3%. Errors in the accuracy of individual ring radii determination became apparent at distortions of 0.4% and higher. At a distortion of 0.4%, the calculated radius of four of the seven rings was in error, albeit by a cumulative error of less than four pixels. The centre error—the standard deviation of the average centre position determined from all the rings, generally increased with distortion. It is clear that Table 2 Influence of elliptical distortion on CHT analysis of Synthetic Pattern 2 comprising seven concentric rings of various diameters, centred on x ¼ 474, y ¼ 519 in a 1024 1024 image Distortion Calculated (%) centre (x:y)
Maximum centre error (%)
No. of ring Cumulative ring error (pixels)b errorsa
0 0.1 0.2 0.4 0.6 0.8 1.0
0.00 0.00 0.00 0.07 0.08 0.52 0.34
0 0 0 4 5 7 7
a
474.0:519.0 474.0:519.0 474.0:519.0 474.2:518.6 474.4:519.2 476.5:519.1 475.6:519.4
0 0 0 3.5 9.7 13.5 11.3
No. of rings where the calculated radius is in error by 470.1 pixels. absolute value (calculated radiusactual radius).
bP
ARTICLE IN PRESS 372
D.R.G. Mitchell / Ultramicroscopy 108 (2008) 367–374
elliptical distortion has the potential to degrade the accuracy of the CHT analysis, and should therefore be minimised prior to pattern acquisition. 5. Assessment with experimental patterns Fig. 4a shows a SADP of a polycrystalline gold film. This represents the best case pattern for this type of analysis, since the rings are sharp and continuous and there is no beam stop present. The pattern was calibrated using the camera length previously determined for our instrument. It was then analysed using the CHT script. Fig. 2 shows the radial intensity profile from this pattern, and the ROIs fitted by the automatic THF routine. The THF may produce false peak identifications due to noise or steeply sloping backgrounds. The ‘peak’ at a radius of E20 pixels is an example of this, and this was deleted prior to CHT analysis. The THF may not locate all the peaks which are present, such as the weak peaks at radii of E250 pixels and E440 pixels and the unresolved doublets at radii of E210 and E275. Some refinement of the THF parameters (‘Brim’, ‘Hat’ and ‘Level’ (Fig. 1)) may help improve THF performance. The peaks missed by the THF were manually located. It is essential that the automatic thresholding routine invoked by the ‘Thresh’ button, can threshold the peak from the background. Very weak or noisy peaks which are difficult to threshold can be omitted from the CHT analysis. Following accurate determination of the pattern centre using the intense peaks, the ‘Rotate’ function can be used to produce the rotationally averaged radial intensity profile. The position of the omitted weak peaks can then be determined manually. A total of 14 Bragg peaks were measured for the gold pattern (Fig. 4) (Note: six additional peaks not located by the THF in Fig. 2 were manually located). Table 3 shows results of a series of independent measurements on the gold pattern at various Hough resolutions (‘Res’ in Fig. 1). The accuracy of the CHT analysis was judged by plotting the
measured d-spacings of the gold versus the actual d-spacings, and measuring the slope (‘Slope’ in Table 3). The quality of the fit was assessed by the R-value of the linear regression (R in Table 3). Clearly, in the first instance, the accuracy of any such d-spacing measurement on a TEM is dependent upon the accuracy of the original camera length calibration. However, any errors in that calibration are common to all the data determined here. At the highest resolution (Res ¼ 1), which took 500 s to complete, the slope of the plot of measured d-spacing versus actual d-spacing was 1.0079. At a resolution of 64 (the lowest resolution), the analysis time was 73 s, and the measured slope was 1.0062. R-values at all resolutions were also very close to 1 in all cases. The measured centres of the pattern at all resolutions were also very consistent, with the largest standard deviation of the average centre position corresponding to just 0.2%. Measurement errors affect the d-spacing values of the innermost diffraction rings (smallest radii) the most. At all resolutions, the error in the reported d-spacing for the innermost diffraction ring, was 0.7%. This error includes any errors arising from the camera length calibration, as well those from the CHT analysis. Fig. 4b shows a worst case experimental diffraction pattern for CHT analysis. It is from a heavily textured TiO2 Table 3 Influence of Hough resolution on the accuracy of d-spacing determination from CHT analysis of a SADP of polycrystalline gold Resolution
Time (s)
Centre (x:y)
Max centre error (%)
Slopea
Rb
1 5 10 32 64
500 180 100 82 73
719.2:480.7 719.2:481.0 719.2:480.9 719.4:480.7 719.2:480.9
0.08 0.07 0.08 0.08 0.14
1.0079 1.0079 1.0084 1.0054 1.0062
0.99999 0.99999 0.99999 0.99999 0.99998
a
Slope of plot of CHT measured d-spacings versus actual d-spacings. Linear regression R-value of plot in footnote a.
b
Fig. 4. The CHT analysis superimposes a coloured mask on the SADP, showing the position and size of the rings: (a) polycrystalline gold: the SADP (left) and the mask (right) are shown separately for clarity; (b) textured TiO2 (anatase) showing that despite incomplete rings the CHT analysis has located the majority of Bragg peaks.
ARTICLE IN PRESS D.R.G. Mitchell / Ultramicroscopy 108 (2008) 367–374
(anatase) film [12], with two preferred orientational variants. The Bragg rings are incomplete and a beam stop outline runs through the pattern. Due to the incomplete diffraction rings it was necessary to use the rotational averaging facility during profile generation, to locate diffraction features which would not otherwise be intersected by a simple linear profile. The CHT analysis correctly located most of the principal Bragg peaks. The standard deviation on the average centre position (0.35%) was somewhat higher than that of the gold pattern under similar conditions (0.1%). The error in the d-spacing determined for the innermost diffraction ring (most sensitive to measurement error) was 1.6%, which compared with 0.7% for the innermost diffraction ring in the gold pattern. 6. Discussion Accurate measurement of electron diffraction patterns poses a number of problems. The uncertainty in the location of the pattern centre precludes accurate measurement of ring radii [9]. Centre of mass approaches to determining the pattern centre [3] are critically dependent upon a centro-symmetrical distribution of intensities in the patterns, and this symmetry may be lost due to the presence of a beam stop or an off-axis condition. Even with a known centre, other pattern effects, such as ring spottiness due to grain size/sampling effects and preferred orientation makes any radial intensity profile very sensitive to the positioning of the extraction line. Rotationally averaging a diffraction pattern about a centre can overcome these two effects [7]. However, this approach is critically dependent upon accurately knowing the centre, and errors in the pattern centre will lead to broadening of peaks in the intensity profile. Pattern centres can be crudely estimated via superposition of a circular feature onto a diffraction ring or by taking the average intercept of g vectors. However, the likely uncertainty of these approaches is probably of the order 75 pixels for a typical pattern captured at 1k 1k resolution. For a high index diffraction ring of radius E100 pixels, this leads to an error of 5%, in the d-spacing determination, which for many measurements is unsatisfactory. Distortion induced by the projector lens [1] can also compound these errors, due to the direction dependence of the measured d-spacing. Pattern measurement can also be time consuming, and so a readily implemented method of performing these measurements accurately and with limited user involvement is desirable. The approach of Hou [6] goes some way towards addressing this need. The approach has been to determine the pattern centre by iterating from a guessed starting position, to minimise the difference between the original image, and one rotated about the centre. This provides not only the pattern centre with high accuracy and precision, but also a measure of the pattern ellipsicity induced by the projector lens. Comparison of the centre determination of
373
the gold pattern (Fig. 4) using the CHT approach and that of Hou’s showed consistent sub-pixel agreement. The centre of mass approach was wildly inaccurate on this pattern. While the gold pattern was centro-symmetric, the weak ‘comet tail’ caused by the electronic shuttering of the CCD camera deflecting the intense transmitted beam off the CCD, was sufficient to shift the centre of mass away from the true pattern centre. As well as providing the pattern centre, the CHT automatically locates the diffraction rings and measures them with sub-pixel precision. Provided the pattern is calibrated beforehand, this yields a listing of d-spacings. Comparison of the measured d-spacings for gold with those expected showed excellent agreement (Table 3). Errors in pattern measurement are most acute for low index diffraction rings of small radius. In the case of the gold measurements in Table 3, the difference between the measured d-spacing and the actual d-spacing of the innermost diffraction ring was just 0.7%. Some of this error is due to the error in the camera length calibration. The radius of this ring was 106 pixels, and referring to the earlier discussion on the likely accuracy of crude estimates of the centre position, an uncertainty of 75 pixels in the centre position, corresponds to a potential error in the d-spacing of up to 5%. The results of the CHT are clearly quite superior. The CHT is also robust when dealing with patterns containing incomplete rings. Russ et al. [10] have shown that textured or single crystal patterns from asbestos polymorphs can be analysed using this method. The highly textured anatase pattern (Fig. 4b) was successfully analysed using the CHT method. The uncertainty in the location of the pattern centre was higher (0.35%) than that of the gold (0.1%), and the error in the d-spacing of the innermost diffraction ring (1.6%) was also higher than that of gold (0.7%). A perfectly circular feature in real space, will produce an intensity distribution in circular Hough space which is radially symmetrical, and contains a sharply defined maximum—the position of which defines the circular feature’s centre. Fig. 3b shows a symmetrical Hough space derived from such a circular feature. Circular Hough space is an accumulator space derived from overlapping circles, and the more pixels that are present on any given radius, the greater the peak value in the corresponding Hough space. If one considers two touching circles, the centres of which are at the 12 and 6 o’clock positions on a circular feature, the position of the overlap of these two circles (the clock face centre) is very sharply defined along the vector joining 12–6 o’clock. However, it is poorly defined at right angles to that vector (9–3 o’clock). If additional circles are added at the 3 and 9 o’clock positions, the centre becomes much more sharply defined, and improves with additional circles at other clock positions. The textured anatase pattern (Fig. 4b) consists of limited numbers of spot pairs which results in a lower accuracy of centre and radius determination compared with the polycrystalline gold ring pattern (Fig. 4b) for the
ARTICLE IN PRESS 374
D.R.G. Mitchell / Ultramicroscopy 108 (2008) 367–374
above reasons. Nevertheless, the performance of the CHT on the textured anatase pattern was satisfactory. The CHT is a computationally intensive method, and its implementation as a useable script for DM has necessitated a number of steps to ensure processing time is not excessive. The creation of a compiled plugin to handle the Hough processing resulted in a 50-fold increase in processing speed. The use of an intensity profile in which the user selects the radii of interest, rather than sampling the entire range of radii which may be present, results in a further 5–10-fold decrease in processing time. The inclusion of a resolution (‘Res’) parameter to allow sub-sampling of real space allows the user to tailor the processing time to keep analysis times reasonable. Significant sub-sampling can be used without accuracy adversely suffering (Tables 1 and 3). The default value (Res ¼ 5) results in a E2–3-fold reduction in processing time, compared with the maximum resolution. The result of these optimisations is that a typical 1k 1k diffraction pattern can be analysed with a processing time of just a few minutes, which certainly compares very favourably with fully manual approaches to measurement. The compiled plugin makes the CHT function available to the DM scripting language, and examples of its use are included in the instructions which can be downloaded with the package. This permits users to incorporate this functionality into their own scripts. Potential applications for this include automatic location/measurement of any circular features, such a convergent beam electron diffraction pattern discs, circular holes/pores/particles, etc. 7. Conclusions The circular Hough transform (CHT) has been successfully incorporated into a user-friendly software tool for automatic analysis of selected area electron patterns within Digital MicrographTM. This software tool can be freely downloaded from the Internet. The computationally intensive CHT method has necessitated the development of a compiled plugin to handle the key Hough transform processing, and a number of software optimisations to permit analysis of diffraction patterns in just a few minutes. The accuracy and robustness of the software has been tested with simulated diffraction patterns. The influence of projector lens distortion on the CHT analysis has been assessed, and the need to minimise any such distortion highlighted. Assessment of best case (polycrystalline gold) and worst case (textured anatase) experimental selected
area diffraction patterns, showed that the CHT software tool was robust and accurate, determining the pattern centres with a maximum uncertainty of 0.1% and 0.35%, respectively. Errors in the measurement of the d-spacing of the innermost diffraction ring (most sensitive to measurement errors), were 0.7% and 1.6%, respectively, although this error did include the error from the camera length calibration of the microscope. The incorporation a CHT function into the Digital MicrographTM scripting language will enable users to develop their own scripts for location/ measurement of circular features. Acknowledgements Dr. Bernhard Schaffer (Graz University of Technology, Austria) is thanked for provision of a C++ function to sort images into a list. This function was incorporated into the CircularHoughTransform.dll described here. Christoph Koch and Douglas Hague are thanked for contributing useful demonstration code of C++ and its compilation to the Digital Micrograph Script Database [5]. Dr. Bill Bertram (ANSTO) is thanked for mathematical assistance with algorithm development. Darren Attard (ANSTO) is thanked for testing this script. References [1] G.C. Capitani, P. Oleynikov, S. Hovmo¨ller, M. Mellini, Ultramicroscopy 106 (2006) 66. [2] Digital MicrographTM, /http://www.gatan.com/imaging/downloads. phpS, Gatan Inc., Pleasanton, CA, 2007. [3] DIFPACKTM, Diffraction analysis package for Digital MicrographTM, Gatan Inc., Pleasanton, CA, 2007. [4] D.R.G. Mitchell, B. Schaffer, Ultramicroscopy 103 (2005) 319. [5] Digital MicrographTM Script Database: /http://www.felmi-zfe. tugraz.at/dm_scripts/dmscript1.htmlS, 2007. [6] V. Hou, Diffraction Rings Distortion-Analysis, Digital Micrograph Script in Digital MicrographTM Script Database: /http://www.felmi-zfe. tugraz.at/dm_scripts/dmscript1.htmlS, 2006. [7] D.R.G. Mitchell, SADP Rotational Average, Digital Micrograph Script in Digital MicrographTM Script Database: /http://www.felmi-zfe. tugraz.at/dm_scripts/dmscript1.htmlS, 2004. [8] P.V.C. Hough, Method and means for recognizing complex patterns, US Patent, 3069654, 1962. [9] J.C. Russ, D.S. Bright, J.C. Russ, M. Hare, J. Comput. Assist. Microsc. 1 (1989) 3. [10] J.C. Russ, T. Taguchi, P.M. Peters, E. Chatfield, J.C. Russ, W.D. Stewart, Adv. X-Ray Anal. 32 (1989) 593. [11] J.C. Russ, The Image Processing Handbook, CRC Press (1996). [12] D.R.G. Mitchell, D.J. Attard, G. Triani, J. Cryst. Growth 285 (2005) 208.