Applications of image processing in combustion research

Digital image processing has wide ranging applications in combustion research. The analysis of digital images is used in practically every scale of studying combustion phenomena from the scale of individual atoms to diagnosing and controlling large-scale combustors. Digital image processing is one of the fastest-growing scientific areas in the world today. From being able to reconstruct low-resolution grayscale images from transmitted signals, the capabilities have grown to enabling machines carrying out tasks that would normally require human vision, perception, and reasoning. Certain applications in combustion science benefit greatly from recent advances in image processing. Unfortunately, since the two fields - combustion and image processing research - stand relatively far from each other, the most recent results are often not known well enough in the areas where they may be applied with great benefits. This work aims to improve the accuracy and reliability of certain measurements in combustion science by selecting, adapting, and implementing the appropriate techniques originally developed in the image processing area. A number of specific applications were chosen that cover a wide range of physical scales of combustion phenomena, and specific image processing methodologies were proposed to improve or enable measurements in studying such phenomena. The selected applications include the description and quantification of combustion-derived carbon nanostructure, the three-dimensional optical diagnostics of combusting pulverized-coal particles and the optical flow velocimetry and quantitative radiation imaging of a pilot-scale oxy-coal flame. In the field of the structural analysis of soot, new structural parameters were derived and the extraction and fidelity of existing ones were improved. In the field of pulverized-coal combustion, the developed methodologies allow for studying the detailed mechanisms of particle combustion in three dimensions. At larger scales, the simultaneous measurement of flame velocity, spectral radiation, and pyrometric properties were realized.

APPLICATIONS OF IMAGE PROCESSING IN COMBUSTION RESEARCH by Pal Toth A dissertation submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Doctor of Philosophy Department of Chemical Engineering The University of Utah August 2014 Copyright © Pal Toth 2014 All Rights Reserved The Unive r si t y of Utah Gradua te School STATEMENT OF DISSERTATION APPROVAL The dissertation of ______________________Pal Toth has been approved by the following supervisory committee members: Eric G. Eddings Terry A. Ring JoAnn S. Lighty Arpad B. Palotas Ross T. Whitaker Chair 3/18/2014 Date Approved Member 3/18/2014 Date Approved Member 3/18/2014 Date Approved Member 3/18/2014 Date Approved Member 3/18/2014 Date Approved and by ___________________ Milind Deo___________________ , Chair/Dean of the Department/College/School o f ______________Chemical Engineering___________ and by David B. Kieda, Dean of The Graduate School. ABSTRACT Digital image processing has wide ranging applications in combustion research. The analysis of digital images is used in practically every scale of studying combustion phenom­ena from the scale of individual atoms to diagnosing and controlling large-scale combustors. Digital image processing is one of the fastest-growing scientific areas in the world today. From being able to reconstruct low-resolution grayscale images from transmitted signals, the capabilities have grown to enabling machines carrying out tasks that would normally require human vision, perception, and reasoning. Certain applications in combustion science benefit greatly from recent advances in image processing. Unfortunately, since the two fields - combustion and image processing research - stand relatively far from each other, the most recent results are often not known well enough in the areas where they may be applied with great benefits. This work aims to improve the accuracy and reliability of certain measurements in combustion science by selecting, adapting, and implementing the appropriate techniques originally developed in the image processing area. A num­ber of specific applications were chosen that cover a wide range of physical scales of combustion phenomena, and specific image processing methodologies were proposed to improve or enable measurements in studying such phenomena. The selected applications include the description and quantification of combustion-derived carbon nanostructure, the three-dimensional optical diagnostics of combusting pulverized-coal particles and the optical flow velocimetry and quantitative radiation imaging of a pilot-scale oxy-coal flame. In the field of the structural analysis of soot, new structural parameters were derived and the extraction and fidelity of existing ones were improved. In the field of pulverized-coal combustion, the developed methodologies allow for studying the detailed mechanisms of particle combustion in three dimensions. At larger scales, the simultaneous measurement of flame velocity, spectral radiation, and pyrometric properties were realized. CONTENTS ABSTRACT iii ACKNOWLEDGMENTS vii CHAPTERS 1. INTRODUCTION TO THE TECHNICAL CONTENT............................. 1 2. QUANTITATIVE DIFFERENTIATION OF POORLY ORDERED SOOT NANOSTRUCTURES: A SEMI-EMPIRICAL APPROACH. . . 4 2.1 Abstract ......................................................................................................................... 5 2.2 Introduction .................................................................................................................. 5 2.3 Materials and m e th o d s ............................................................................................... 6 2.4 Results and discussion.................................................................................................. 9 2.5 Conclusions.................................................................................................................... 12 2.6 Acknowlegdment........................................................................................................... 12 2.7 References...................................................................................................................... 12 3. A NOVEL FRAMEWORK FOR THE QUANTITATIVE ANALYSIS OF HIGH RESOLUTION TRANSMISSION ELECTRON MICROGRAPHS OF SOOT I. IMPROVED MEASUREMENT OF INTERLAYER SPACING ............................................................................. 13 3.1 Abstract ......................................................................................................................... 14 3.2 Introduction .................................................................................................................. 14 3.3 Materials and methods ................................................................................................ 15 3.4 Results and discussion.................................................................................................. 18 3.5 Conclusions.................................................................................................................... 24 3.6 Acknowlegdment ........................................................................................................... 24 3.7 References ....................................................................................................................... 24 4. A NOVEL FRAMEWORK FOR THE QUANTITATIVE ANALYSIS OF HIGH RESOLUTION TRANSMISSION ELECTRON MICROGRAPHS OF SOOT II. ROBUST MULTISCALE NANOSTRUCTURE QUANTIFICATION..................................................... 25 4.1 Abstract ......................................................................................................................... 26 4.2 Introduction .................................................................................................................. 26 4.3 Materials and methods ................................................................................................ 27 4.4 Applications - results and discussion........................................................................ 30 4.5 Conclusions.................................................................................................................... 38 4.6 Acknowlegdment........................................................................................................... 38 4.7 References...................................................................................................................... 38 5. AUTOMATED ANALYSIS OF HETEROGENEOUS CARBON NANOSTRUCTURES BY HIGH-RESOLUTION ELECTRON MICROSCOPY AND ON-LINE IMAGE PROCESSING........................ 39 5.1 Abstract ......................................................................................................................... 40 5.2 Introduction .................................................................................................................. 40 5.3 Materials and methods ................................................................................................ 41 5.4 Applications - results and discussion........................................................................ 43 5.5 Conclusion ....................................................................................................................... 48 5.6 Acknowlegdments......................................................................................................... 48 5.7 References ....................................................................................................................... 48 6. DETAILED INVESTIGATION OF SOOT NANOSTRUCTURE: EFFECT OF PR E S SU R E 50 6.1 Abstract ......................................................................................................................... 50 6.2 Introduction .................................................................................................................. 51 6.3 Materials and methods ................................................................................................ 52 6.4 Results and discussion .................................................................................................. 55 6.5 Conclusion ....................................................................................................................... 66 7. THREE-DIMENSIONAL STREAK IMAGING OF COMBUSTING COAL PARTICLES I. VELOCIMETRY 67 7.1 Abstract ......................................................................................................................... 67 7.2 Introduction .................................................................................................................. 67 7.3 Method description....................................................................................................... 69 7.4 Uncertainty in streak localization............................................................................. 75 7.5 Experimental.................................................................................................................. 84 7.6 Results of laboratory-scale experiments ................................................................. 87 7.7 Conclusion...................................................................................................................... 93 8. THREE-DIMENSIONAL STREAK IMAGING OF COMBUSTING COAL PARTICLES II. PYROMETRY 95 8.1 Abstract ......................................................................................................................... 95 8.2 Introduction .................................................................................................................. 95 8.3 Method description ....................................................................................................... 97 8.4 Uncertainty analysis .................................................................................................... 105 8.5 Experimental.................................................................................................................. 109 8.6 Experimental re su lts .................................................................................................... 109 8.7 Possible improvements ............................................................................................... 117 8.8 Conclusion......................................................................................................................117 9. THE POTENTIAL OF ON-LINE OPTICAL FLOW MEASUREMENT IN THE CONTROL AND MONITORING OF OXY-COAL FLAMES 120 v 9.1 Abstract ......................................................................................................................... 121 9.2 Introduction .................................................................................................................. 121 9.3 Materials and m e th o d s ............................................................................................... 122 9.4 Flame monitoring by O F V ........................................................................................ 124 9.5 Remarks ......................................................................................................................... 129 9.6 Conclusion......................................................................................................................132 9.7 Acknowlegdments.........................................................................................................132 9.8 References......................................................................................................................132 10. IMAGING OPTICAL DIAGNOSTICS OF 40 KW CO-AXIAL OXY-COAL FLAMES...............................................................................................134 10.1 Abstract ......................................................................................................................... 134 10.2 Introduction .................................................................................................................. 134 10.3 Materials and methods ................................................................................................ 135 10.4 Results and discussion .................................................................................................. 141 10.5 Conclusion ....................................................................................................................... 146 11. SUMMARY...................................................................................................................147 REFERENCES ...................................................................................................................150 vi ACKNOWLEDGMENTS This material is based upon work supported by the Department of Energy under Award Number DE-NT0005015. The views and opinions of the authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof. The author would like to thank Dana Overacker, David Wagner, Robert Cox, Tom Slowik, Prof. Philip Smith, and Mark Solum for the help they provided throughout this work. CHAPTER 1 INTRODUCTION TO THE TECHNICAL CONTENT This dissertation is a compilation of standalone articles. Each following chapter summa­rizes an approach to a distinct group of problems, along with experiments and experimental results, where applicable. Although it may seem that the problems that are discussed in Chapters 2-10 are at very different spatial scales, there is an underlying concept in the methodologies applied to solve these problems. The underlying concept in the methodologies is applying digital imaging, image process­ing and image analysis to study physical phenomena of a very different nature. As digital imaging is a well-defined concept, the methods applied are similar and are based on the same principles. This thesis therefore demonstrates how a mathematical concept can be applied in the study of diverse physical systems. Another, although more loose underlying similarity is found in the objects of studies themselves. Every chapter describes a specific problem that arose in the general scientific area of combustion research. Chapters 2-6 deal with soot formation in the combustion of liquid fuels (Diesel surrogates), while Chapters 7-10 investigate problems in pulverized-coal combustion. Chapter 2 is focused on the problem of quantifying soot nanostructure, as it is ob­served by high-resolution electron microscopy (HRTEM). HRTEM micrographs are two­dimensional intensity arrays (images), that convey information about the molecular struc­ture of soot particles. Interpreting this information in a quantitative fashion is a challenging task as the raw data must be heavily reduced in order to obtain meaningful scalar quantities that describe the three-dimensional structure. Chapter 2 presents an early attempt of this data reduction problem. As a result of this work, structural parameters (scalars) were developed that robustly quantify nanostructure. Chapters 3 and 4 present an advanced mathematical framework with which the struc­tural descriptors presented in Chapter 2 can be both generalized and their physical meanings 2 made more intuitive. The novelty of the approach presented in Chapters 3 and 4 lies in the application of filtering theory to the HRTEM micrographs. Image filtering is normally used to enhance, denoise, or reconstruct images or to achieve artistic effects. In these applications, image filtering is employed to extract physical information, like the spatial location and arrangement of molecules from the micrographs. The chapters up to Chapter 4 all deal with studying physical phenomena on the nanoscale, that is, on the scale of individual molecules. Chapters 5 and 6 build on Chapters 3 and 4, although in different ways. Chapter 5 presents software that applies the theory and implementation developed in Chapters 3 and 4 towards the integrated, automatic control of microscopes. The computational efficiency of the algorithms described in Chapter 4 allows for the utilization of the image filtering technique in real-time applications. In Chapter 6, the methodology presented in Chapters 3 and 4 is applied to extract physical information from "real" HRTEM micrographs and to study the effect of oxidation pressure on soot nanostructure. Chapters 5 and 6 can be viewed as practical applications of the methodology developed in the previous chapters. The methodology presented in Chapter 5 can extend the capability of the techniques described in Chapters 3 and 4 by being able to study physical structure on the microscale, that is, on the scale of soot aggregates. Chapters 7-10 deal with a different group of problems, namely the combustion of pulverized-coal. There are a number of physical and chemical aspects of coal combustion that have been widely studied in the literature and are important to know in order to be able to efficiently use coal as an energy resource. These aspects include the chemical kinetics of coal particle combustion, the fluid mechanics of pulverized-coal flames, and multiphase, multicomponent mass transfer in such systems. Chapters 7-10 focus on different parts of the problem, but the concept of increasing spatial size scales with each chapter is retained. Chapters 7 and 8 deal with the study of coal particle combustion, studied at the spatial scale of single particles. This is the scale of millimeters. The methodology described in these chapters allows for the nonintrusive optical measurement of the three-dimensional location, three-component velocity, temperature, and size of individual, combusting coal particles. The method is a combination of stereoscopic imaging, image processing, photometry, and pyrometry. The individual approaches developed in these different fields are combined in a novel way. Chapters 9 and 10 study pulverized-coal flames on the spatial scale of centimeters and meters. At this scale, it is not possible to observe individual coal particles combusting, but 3 their behavior can be deduced by observing the macroscopical characteristics of the flame. These characteristics are the radiative emission at different wavelengths and spatio-temporal properties like flame shape, flame length, and the apparent motion of flame structure. Chapter 9 describes a method with which a descriptor of flame velocity can be extracted from simple (although high-speed) photographs of flames. The method is tested on a 40 kW (nominally 120 kW) pilot-scale pulverized-coal flame. Chapter 10 utilizes an extension of the velocimetry method presented in Chapter 9, in which the capability of velocity measurement is extended by pyrometric capabilities that allow for the extraction of luminous fluxes in multiple wavebands, along with flame temperature and soot concentration. The combined capability of measuring these properties by using reasonably simple concepts of imaging, image processing, machine vision, and spectroscopy is novel. The method is applied to study the effect of oxygen injection and coal type on the already described 40 kW flame. A large part of the work (Chapters 2-5 and 9) has been submitted to and accepted by international scientific journals for publication. The remainder of the dissertation consists of manuscripts that are either under peer review or submission procedure at the time of the submission of this dissertation. For this reason, every chapter follows a journal article organization and the original text of the manuscripts is compiled as separate chapters of this dissertation. The following chapters have been published: • Chapter 2 - Quantitative differentiation of poorly ordered soot nanostructures: a semi-empirical approach. Published in Fuel, volume 99, pages 1-8 (2012), • Chapter 3 - A novel framework for the quantitative analysis of high resolution transmission electron micrographs of soot I. Improved measurement of interlayer spacing. Published in Combustion and Flame, volume 160, issue 5, pages 909-919 (2013), • Chapter 4 - A novel framework for the quantitative analysis of high resolution transmission electron micrographs of soot II. Robust multiscale nanostructure quan­tification. Published in Combustion and Flame, volume 160, issue 5, pages 920-939 (2013), • Chapter 5 - Automated analysis of heterogeneous carbon nanostructures by high-resolution electron microscopy and on-line image processing. Published in Ultrami­croscopy, volume 129, pages 53-62 (2013). • Chapter 9 - The potential of on-line optical flow measurement in the control and monitoring of oxy-coal flames CHAPTER 2 QUANTITATIVE DIFFERENTIATION OF POORLY ORDERED SOOT NANOSTRUCTURES: A SEMI-EMPIRICAL APPROACH Reprinted from Fuel, Vol. 99, P. Toth, A.B. Palotas, J. Lighty, C.A. Echavarria, Quan­titative differentiation of poorly ordered soot nanostructures: a semi-empirical approach, Pages 1-8, Copyright 2012, with permission from Elsevier. 5 Fuel 99 (2012) 1-8 ELSEVIER Contents lists available at SciVerse ScienceDirect Fuel jo u r n a l h om e p a g e : w w w .e ls e v ie r . c o m / lo c a te / fu e l Quantitative differentiation of poorly ordered soot nanostructures: A semi-empirical approach P. Totha'*, A.B. Palotasa, J. Lightyb, C.A. Echavarriab a Department o f Combustion Technology and Thermal Energy, University o f Miskolc, H3515 Miskolc-Egyetemvaros, Hungary b Department o f Chemical Engineering, University o f Utah, 50 S. Central Campus Drive, Salt Lake City, UT 84112-9203, United States A R T I C L E I NF O A B S T R A C T A novel, sensitive technique is presented yielding expressive semi-empirical order parameters (distance deviation parameter X and junction parameter t) describing the nanostructure of combustion derived carbonaceous materials. This method is purposively developed to enable the measurement of slight changes in structural properties of young soots originating from the combustion of different fuels - most of which are left undetected by conventional algorithms - by utilizing digital image processing of high resolution transmission electron micrographs (HRTEM images). The study has been motivated by a demand for better understanding the evolution of soot nanostructure. This paper presents the details of the proposed method, along with the comparison with the conventional methodology. Through sim­ulations, the correlation between previously defined order parameters and our new structural indices is demonstrated. As an example of practical utilization, the nanostructural evolution of different fuels with residence time and the extent of oxidation is shown by monitoring changes in X and t. The struc­tural development was analyzed using a specific sampling procedure: soots were collected at various heights above the burner surface. Published by Elsevier Ltd. Article history: Received 13 July 2010 Received in revised form 11 March 2012 Accepted 5 April 2012 Available online 21 April 2012 Keywords: Soot Nanostructure HRTEM Image processing 1. In tro d u ctio n Dark field lattice fringe imaging techniques utilizing high reso­lution transmission electron m icroscopes allow for the observation of soot structure a t the atomic levels. Several studies have been published discussing the properties of soot nanostructure based on qualitative [1-4] or quantitative [5-9] image analysis methods. The atomic-level investigations of combustion-derived carbon structures are generally initiated by either an a ttempt to determine the source of soot p o llu tio n [10,11] or, more universally, a m otive to gain a better insight on oxidation and graphite crystal formation. The latter incentive is important from further points of view, since the characteristics of soot nanostructure are believed to be corre­lated w ith macroscopic physical properties, e.g. oxidation reactivity and UV-VIS absorption [12] or radiant emissivity and extinction coefficients [13]. Others argue th a t the short or long range order of lattices has an effect on microstructural properties, porosity and surface area as well [14]. Much effort has been put in the attempts to representatively quantify structural information provided by HRTEM lattice images and several existing m ethodologies are known th a t may yield con­sistent data as a result. While the details of these algorithms are * Corresponding author. E-mail address: [e-mail] (P. Toth). 0016-2361/$ - see front m a tte r Published by Elsevier Ltd. http://dx.doi.org/10.1016/j.fuel.2012.04.013 quite different, the basic approaches and main steps of extracting physical information are generally the same. These methods usu­ally consist of the following: 1. the noise filtering of raw, digitized HRTEM images, 2. the classification of pixels as either objects of interest or background, 3. the measurement of some shape properties of these objects and the statistical evaluation of the data collected. Additionally, several studies utilize cluster or stack analysis in addition to these conventional steps in an attempt to quantify the order of objects and crystallinity of soot structure. HRTEM images are usually digitized as 8 bit grayscale bitmaps, where the value of each pixel - generally referred to as the gray­scale intensity value (/) - determines the corresponding pixel's probability of belonging to an object. Objects are believed to be graphene layers, also called fringes or lattices. Noise filtering basically means the removal of regional peaks in images caused by electrostatic effects, not carrying any structural meaning. Because of the characteristics of carbon soot nanostruc­tures - the definite range of fringe spacing and length has already been described by many authors who used different techniques [15-17] - the best results are usually obtained when using frequency domain filtering. Square-profiled bandpass frequency filters (referred to as ‘ideal' filters) were the most commonly used 6 P. Toth et al./Fuel 99 (2012) 1-8 filters in initial studies [5]. Sharma and co-workers [6] pointed out th a t important non-periodically occurring single layers are lost when using bandpass filtering with amorphous soot samples; their study proposed the use of low-pass ideal filters. For the detection process - the separation of objects and back­ground - basically all authors used binarization methods based on global thresholding. Choosing the appropriate method for the binarization step is crucial for obtaining consistent results, as the detected objects are the inputs of the next steps. Global threshold­ing means th a t a single 8 bit threshold value T is set before the detection procedure, and the values Bx,y o f the output binary image B are either 1 if Ix,y p T or 0 otherwise. The outputs of global thres­holding algorithms are therefore heavily determined by T, so it is not surprising th a t the m ethods used for assigning threshold values have been subjects of debate. Basically, there are authors w ho set T in order to obtain separated binary fringes in B [5] and authors w ho set T to get equal fractional coverage values for objects and non-ob­jects and used post-processing techniques to separate fringes after binarization [6,7]. These post-processing methods include the skel­etonization of fringes, detection of branchpoints, separation and reconnection algorithms based on complex criterion systems [6 -8 ]. From the structural data obtainable from HRTEM images, the most commonly measured parameters by recent studies are the interlayer distance, the fringe length and the fractional coverage of detected objects [5,6,8]. These are also the most easily verifiable by other techniques, such as X-ray diffractometry [18]. For the di­rect determination of the interlayer distance, to the best of our knowledge, all authors used the approximation proposed by Palotas and coworkers [5]. According to this study, the interlayer distance of two adjacent carbon fringes is the vector-vector distance of the two paralellized orientation vectors, where these vectors are obtained by averaging their orientations but keeping them at their original location, i.e. passing through the centroid of the fringe. Most recent publications place more emphasis on the classification o f ‘adjacent' and parallel groups of fringes, called ‘stacks'. Detailed algorithms for deciding an object's membership to stacks are pre­sented in [6,8]. If considering data describing the stacking oflayers, additional parameters arise, such as the diameters of stacks, the number of layers in stacks and other derived indices, e.g. the crys-tallinity index proposed byYang et al. [8]. Usingwhole distribution functions of the interlayer distance and the fringe orientation h can be expressive, because the amorphity of the texture - the deviation from perfect graphite structure - can be described by the statistical descriptors of these distributions. Shim and coworkers [7] intro­duced several interesting o rder parameters, such as the 2D nematic order parameter S2,N and the 2D polar o rder p arameter S2,P, defined as 2 S2A = 2cos2( h i ) - 1 (1) and S2,p = 1 - 2cos2(hi) (2) where h is a vector containing the angles between the orientation vectors of each fringe and a directional vector; the overline symbol means the arithmetic mean value of the elements of its argument. The directional vector is the reference vector for S2,N and the vector pointing from the concentric center to the center of mass of fringe i. For a perfectly ordered graphite structure, where all fringes are stacked together and their orientations are the same, the nematic order parameter is 1; for a perfectly concentric, onion-like carbon structure, the polar order parameter is 1. The aim of this research is to develop a method, with which even slight changes in amorphous soot nanostructures can be m onitored. Our previous experiences with conventional algorithms showed th a t the inherent uncertainty of these methods - mainly coming from the subjectivity of the b inarization process and the small sam­ple sizes due to the large amount of excluded data - makes the characterization of very similar structures difficult. 2. Materials and m e th o d s 2.1. Sample collection Soot samples were derived from the flames of ethylene and air, benzene and air and a surrogate fuel consisting of n-dodecane and m-xylene (9:1 V/V) and air. Samples were thermophoretically ta ­ken for HRTEM analysis. Ethylene and benzene were burned in air in a flat flame premixed burner under fuel rich conditions for studies of soot formation in aliphatic and aromatic flames [19]. This system consisted of a stainless steel chamber where fuel and air were properly mixed prior to entering the burner. The flame was stabilized over a tube bundle and shielded from atmo­spheric interference using a nitrogen shroud. A metallic mesh placed 3.5 cm above the burner surface distributed the flame uni­formly across the burner. The surrogate mixture and ethylene w ere also burned in air in a two stage premixed burner for studies of soot oxidation. In this system, soot was generated in either ethylene/air or surrogate/air fuel rich flames, which served as the first stage. The soot was then burned in a secondary, premixed burner under fuel lean condi­tions. Details of the premixing systems have been reported earlier [20]. Samples for HRTEM analysis were taken at different heights above the burner surface (HAB) using a thermophoretic probe commonly referred to as a frog tongue. A TEM grid holder w as a t­tached to a piston and compressed air at 60 psig was used to quickly insert the TEM grid into the flame. Multiple insertions were necessary to get a representative soot sample on the grid. The grid was oriented with the face parallel to the gas flow, so the distur­bance of the flame was minimal. Soot deposits on the grid because of the thermophoretic gradient between the cold grid and the hot flame, allowing freezing some heterogeneous reactions, avoiding changes in the soot morphology after the particles have been im­pacted upon the cold surface [21]. A summary of the collected samples, heights above burner and flame types is presented in Table 1. 2.2. Image acquisition HRTEM micrographs were produced using transmission elec­tron microscopes FEI, Models Tecnai F30 and F20 EFTEM operated under 200 keV accelerating potential. The images were digitized as 1024 x 1024 size images. The magnifications of these images w ere either 760,000x or 1,100,000x. 2.3. Image processing Our approach is fundamentally different from the previously described ones. Instead of trying to separate fringes based on hypothesized criteria and extract data from bitmaps, we tried to reduce the images to ‘networks' th a t can be topologically charac­terized. These topological parameters are correlated to the already defined order parameters later. For noise filtration, Gaussian bandpass filters have been used in the frequency domain. Gaussian filters are better than ideal filters if the harmonic ‘ringing' caused by the sharp changes in the transfer function may lead to false fringe detection, as in the case of high resolution soot TEM images. Despite the advantages of using a low-pass filter as proposed by Sharma and coworkers [6], bandpass filters w ere used, because of the additional intensity homogenizing 7 P. Toth et al./Fuel 99 (2 0 1 2 )1 -8 Table 1 Summary of collected samples. Sample name Fuel type Burner setup Heights above b u rn e r (mm) Benzene flame Benzene/air Premixed, fuel rich 5, 10, 15 Ethylene flame 1 Ethylene/air Two stage premixed, fuel rich/fuel lean 0.5, 2, 4 Ethylene flame 2 Ethylene/air Premixed, fuel rich 10, 15 Surrogate flame N-dodecane-m-xylene mixture/air Two stage premixed, fuel rich/fuel lean 1, 3, 5 properties of bandass Fourier filtering (inhomogeneous luminance, e.g. large patches of darker or brighter pixels are considered very low frequency noise). Frequency filters work in the frequency do­main of the image, obtained by the Fourier transform. This transfor­mation places the frequency data of the image from the spatial domain x, y into a new coordinate system called the frequency space u, v. This frequency space characterizes the different spatial frequencies of the image. In the following discussion, we use the convenient continuous function notation and we consider a cen­trally shifted frequency space, meaning th a t the lowest frequencies appear in the center of the frequency space and frequency increases as one m oves further away from this location, proportionally to the Euclidian distance from the central coordinate. In a coordinate sys­tem like this, ideal filters can be pictured as cylinders and Gaussian filters can be pictured as two-dimensional Gaussian surfaces, both centered around the origo. The following two equations describe ideal and Gaussian frequency filters, respectively. u2 + v2 6 rLP ! Oilp(u, v) = 1, else 0 OgMu, v) = exp 2 r 2„ (3) (4) Og,BP (u, v) = exp u2 + v2 2 r 2P 1 - exp u2 + v2 2 r HP (5) /■n/2 p m/2 -n/2 J-m/2 i-n/2 <■ m/2 exp u2 + v2 -n/2 -m/2 1 - exp 2 °IP u2 + v2 220r2 H P dudv = m n - nrH (6) (7) shifted coordinate system. Considering tha t most HRTEM images are digitized in m x m sizes, where m is a power of 2 to ease fre­quency filtering, therefore m = n , the integration of the left sides of (6) and (7) yield 2 p r 2Perf2 2 20L m2 - 2 n o 2HPerf2 (- PS- | = m2 - prH HP 2 20 H H (8) (9) It is apparent th a t the two equations show symmetry, therefore both 0 LP and 0HP can be determined by the implicit transcendent equation 0 rXP 2erf 0 xP (10) 2\/2oxp) The index LP refers to low-pass, since the filters allow low fre­quencies around u = 0, v =0. For both filters, there is a parameter controlling the amount of cut frequencies. The ideal filter can be pictured as a cylinder in the frequency space, centered around the u = 0, v = 0 coordinate, with a radius Vlp and a u nit height, there­fore the amount of frequency cut can be controlled by the radius Vlp. In the case of a two-dimensional Gaussian probability distribu­tion function, the parameter 0 LP would be the variance, ergo the ‘width' of the distribution. Thus it is easy to see by analogy, that olp will have the same rule for a Gaussian transfer function as the radius Vlp for an ideal filter. In the following we generalize the same idea for bandpass filters. A transfer function for a band­pass Gaussian filter can be generated by superimposing two Gauss­ian filters similar to the low-pass filter, with one - the high-pass component - being complemented. Thus, the complete function can be defined as follows: where the index HP refers to high-pass. The idea behind the deter­mination of o values is tha t the same extent of filtering can be achieved by equal ‘volume' ideal and Gaussian transfer functions. The ideal bandpass transfer function can be generally pictured as the difference of two cylinders; their radii Vlp and Vhp correspond to low- and high-pass cut-off frequencies, their heights are 1. Hence, to find the functions olp = frLp) and ohp = Kvhp) for an m x n size image, one can write The integration is carried out from -m /2 to m/2 (or n/2), be­cause the Fourier transform will generate a transformed image with the same dimensions as the original image, in a centrally where xp refers to either low-pass (Lp) or high-pass (Hp). In practi­cal situations, w hen Vxp is significantly lower than m, Eq. (10) can be simplified to Oxp rxP' (11) The relative error of this assumption is shown in Fig. 1. By choosing adequate values for Olp and Ohp it is possible to reduce the image only to the information having structural meaning, i.e. both high frequency and low frequency noise can be filtered. In our case, olp and ohp were determined a t cut-off frequencies so th a t the low-pass component filtered out p atterns repeating in less than 0.3 nm (the approximate shortest physically meaningful distance between carbon layers) and the high-pass component filtered out large patches with relative frequencies lower than 10. After frequency filtering, the inverse Fourier transformed image has been re-scaled and saturated to restore the original image for­mat. Since we used an adaptive thresholding technique, further intensity normalization of the restored image was not necessary. For the d etection process, w e used a local, adaptive thresholding technique. This method, along with the rest of the techniques are presented in matrix notation for convenience. Adaptive local threshold methods have been widely used in several image Fig. 1. The re la tive erro r of Eq. (11). 3 = n r n. 8 4 P. Toth et al./Fuel 99 (2 0 1 2 )1 -8 Fig. 2. The main step s o f th e n ew algorithm for th e d e te rm in a tio n o f X. Left: frequency filtered image. Center: binarized image (B) w ith th e thinned network indicated by lines (N). Right: th e grayscale image denote s th e distance transform of N, D and th e locations o f th e local maxima in D are indicated by lines. At th e se locations, th e values of D will be th e members o f th e s e t E, from w h ich X is computed. The padding in th e th ird image is necessary to eliminate errors caused by th e median filter a t image edges. Fig. 3. Examples o f artificially created binary fringe skeleton images (artificial versions of N). Top left: oriented, axially symme tric set. S2,n ~ 0.99. Top right: a shuffled version of th e oriented axial set. S2N « 0 .6 . Bottom left: oriented, concentrically symmetric set. S2,p « 0.99. Bottom right: a shuffled version of th e oriented polar set. S2,p « 0.6. processing areas, e.g. in document binarization and handwriting recognition [22,23]. ‘Local' m eans tha t the value o f Tis not constant for the entire image; instead, T is generalized as a matrix T with the same size as the image, thus a v ariation in the threshold values as a function of the location in the image is allowed. ‘Adaptive' means th a t T (or T) is not preliminarily set; its value is d ependent on o ther factors. For an input 8 bit grayscale image I, each entry TXy of the matrix T was calculated as the median value of the 1 x m sized neighborhood around the corresponding image pixel I^y, including Ix,y itself. As our images were digitized with an aspect ratio of 1, square-shaped median filters were used ( 1 = m). From our point of view, this binarization filter th a t produces the output binary image B, defined in the following way: Bxy = (12) has several advantages. First, median filtering has noise-removal ef­fects, making the edges of the detected objects ‘smoother', i.e. it Fig. 4. The re sults of th e stochastic simulations. Top: th e correlation o f X w ith S2,p. Bottom: th e correlation o f X w ith S2,p tends to remove spur pixels, which is beneficial in the next step. Second, this method considers regional changes in intensity rather than absolute values, thus less focused or less intense fringes can be detected as well. This feature explains the need for a frequency filter that does not produce harmonic ringing artifacts, thus the usage of Gaussian frequency filters. In the next step, the detected pattern is skeletonized using the algorithm presented Zhang and Suen [24]. This iterative algorithm removes edge pixels and leaves only the backbone of the structure, but does not tend to leave spurs and diagonal arms. The remaining 9 P. Toth et al./Fuel 99 (2 0 1 2 )1 -8 lines make a binary network of fringes N, which is the input for topological analysis. N only contains 1 and 0 values and 1 values form thin connected components. Two new structural parameters are proposed th a t can be mea­sured on these connected networks. The distance deviation p aram­e ter X (nm) is the standard deviation of the Euclidean distance set E, while the junction parameter t (1/nm) is the average number of branchpoints per network length unit. In the following, we define these parameters and discuss their physical significance. The calculation o f X can be achieved in a number of steps. First, the Euclidean distance transform o f the binary n etwork A = dist(N) is computed. The Euclidean distance transformed value of a point BXy in a binary set B is its Euclidean distance to the nearest point with the value :Bx,y, the negated value. Because N is a binary dis­crete network without blocks of adjacent 1 values, for coordinates (xo,yo), where N*0y0 = 1, DX0y0 is most likely a regional minimum of A. Note th a t these are the locations of the 1 valued pixels in N, thus it is not surprising th a t they belong to the regional minima of A - they denote the skeletons of the fringes in I. From the above mentioned, it is easy to see th a t points in the regional maxima in A correspond to half of the shortest distances between the centers of the nearest two objects; the connected lines formed by the adjacent points of the regional maxima of A are al­ways positioned halfways between objects in B. In the second step, values belonging to the set of regional maxima are extracted from A. We denote this set as E. Generally, there are two obvious ways of doing this: one could either use numerical gradient operators to find the locations of regional maxima or apply the same skeleton­ization routine to :N as one used to obtain N itself, so th a t the 1 valued pixels in the skeleton of :N will indicate the locations of the local maxima in A. Here the latter was used. Once the set of local maxima of A has been extracted by the above described method, the distance deviation parameter X can be readily computed by computing the standard deviation of the set of local maxima values. Therefore X expresses the amount of variation in the generalized interlayer distances. In other words, X quantifies how m uch the fringe distances deviate from the aver­age interlayer distance in an image, in nanometers. As can be seen, values in E are simply generalized versions of the already well-known interlayer spacing. The generalization originates in extending the applicability of the interlayer distance to fringes th a t are not parallel to each other and /or are curved. The advantages of a structural parameter like X are the signifi­cantly larger sample sizes and the convergent nature of the stan­dard deviation value. By using the original definition of interlayer distances, only one value is obtained from an adjacent pair of fringes - an analogous value in E can be extracted from every sin­gle local maximum in A. Ifthe images are sampled correctly, X will converge to a statistically representative value for the image or set of images coming from the same soot sample, thus the uncertainty in X will be negligible. The disadvantage is tha t although X is sim­ply the average deviation of the generalized interlayer distances, it is more difficult to grasp its direct physical meaning. X character­izes the texture in a complex way, as it is a parameter describing the uniformity of the lattice spacing distribution, the parallelity of the fringes and the order of their orientation as well.1 Its value is 0 for a perfectly graphitic texture containing only parallel, elongated fringes and increases with decreasing orderliness. A dem­onstration of the main steps of the algorithm for ease of understand­ing is presented in Fig. 2. The calculation of th e ju n c tio n parameter t is done by counting the branchpoints on the thinned networks and dividing this value Fig. 5. The verification of th e junction-finding algorithm. Left: a deta il of a HRTEM micrograph tak e n of an oriented Si single crystal. Right: The d ete cted network (lines) and branchpoints (dots). by the measured length o f the network. The detected branchpoints are filtered so th a t forks w ith an arm of length not exceeding a cer­tain threshold are not considered as branches. This w ay spur pixels and arms th a t are the results of the binarization of jagged contours are left neglected. The filtering algorithm works by checking the distance of each branchpoint from the nearest endpoint of an arm, and pixels connecting too closely spaced branch-end pairs are removed from the network. Branchpoints too close to each other are merged by binary dilations. Similar to distance deviation, t is an empirical parameter o f the orderliness ofcarbon texture as well: higher values indicate a m ore entangled structure, as the number of branchpoints increase with decreasing fringe lengths and overlapping layers. Its value is 0 for an ideal graphitic crystal with parallel fringes and no overlap­ping. A more entangled and amorphous structure will exhibit higher values for t. 3. Results and discussion To present the correlation of X with S2,p and S2,n, we have con­ducted several stochastic simulations by using artificially created fringe skeleton images. For the analysis of X = f(S2,jv), the images were constructed to contain short fringes with randomly posi­tioned centers of mass, but not allowing overlapping fringe skele­tons. Their orientations were also randomized, by using a random generator yielding uniformly distributed numbers in the range [-q,q]. For the simulation of concentrically symmetric structures, 1 Note th a t the same is true when using the conventional definition of interlayer spacing. Fig. 6. T h e v e rific a tio n o fth e ju n c tio n -fin d in g a lg o rithm w ith S i1 1 0 lattice images. The continuous line shows th e theore tic al value for Si 1 1 0 images. 5 10 6 P. Toth et al./Fuel 99 (2012) 1-8 a large number of fringe skeletons w ere rendered around a concen­tric center, equally spaced from each other and then perturbing their orientation angle by these random values. In this way, it is possible to observe how X changes with the order parameters S2,P and S2,N. Examples of these artificial images are presented in Fig. 3. By simplifying Eq. (2) by only considering concentric symmetry and if the number of fringes on each image is large, both S2,N and S2,P can be estimated by S2iX * E(2cos2(H ) - 1) (13) where H is a random variable of uniform distribution with values in the range [-q,q] and E() is the expected value of its argument. Solving (13) yields S „,sin(2q) Fig. 4 presents the results of the stochastic simulations. The tests were conducted with images representing gradually decreas­ing o rder parameters in each consecutive test. For each parameter, a total of 10 tests were evaluated, with images containing 400 fringes. It is apparent that the decrease in the order parameters resulted in an increase in X, but the correlation is not linear. The measured values of X increased monotonically in both cases, exhibiting steeply increasing periods at both ends of the curves. W ith our sets, changing the polar order parameter of the concentrically symmetric images induced less predictable changes in the values of the distance deviation parameter and the experimental data w as more scattered. The monotonity th a t can be observed in Fig. 4, along with the fact that X responds to many conditions - including these two order parameters - ensures that even slight changes in structure can be (14) representatively monitored. Another advantage is that X is sensitive Fig. 7. Some examples o f th e analyzed images. (A) Ethylene flame 1, HAB = 0.5 mm. (B) Ethylene flame 1, HAB = 4m m . (C) Surrogate flame, HAB = 1 mm. (D) Surrogate flame, HAB = 5 mm. (E) Benzene flame, HAB = 5 mm. (F) Benzene flame, HAB = 15 mm. 11 P. Toth et al./Fuel 99 (2012) 1-8 in the range of low 82? and S2N values. The downside is tha t it is impossible to gain information about specific structural parameters solely by observing changes in the value of X; it should rather be used as an intuitive empirical parameter of orderliness. For the verification of the junction finding algorithm, we used HRTEM images of Si lattices. At certain magnifications, Si single crystals look like perfect grids, where intense spots are believed to be Si atoms. For such arrangements and for infinitely large fields of view, the ideal value of t can be calculated as follows: u k!™ [l(k - 1)L + k(f- 1)Lj' (15) where k and l are the number of atoms th a t can be counted on the image horizontally and vertically and L « 0.2715 nm, half the lattice spacing of Si single crystals [25] (see Fig. 5 ). The horizontal axis illustrates the number of atoms (k x l) of the detail tha t has been cut from a large lattice image. Since the hori­zontal and vertical dimensions of these details were not the same, an average size of k* = -Jk7 \ (16) is shown. Fig. 6 shows the results of this verification. The continuous line represents the ideal, calculated value of t, while the circles mark the measured results. The errors of the algorithm are most com­monly caused by duplicate branchpoints, crystal imperfections, inaccurate network extraction and image noise, however the mag­nitude of errors is in an acceptable range. The above introduced parameters, X and t were calculated for samples originating from four fuels and with soot samples taken at different heights above burner (HAB). As HAB increases and the soot gets more mature, an increase in orderliness can be ex­pected, as the conventional analysis of [26] could demonstrate a slight shift in fringe lengths in the case of methane soots. An increasing orderliness would suggest a decrease in both X and t. A few examples of the analyzed images can be seen in Fig. 7 . From the real HRTEM images and from each sample, both X and t were computed. Because multiple images had been taken from the same soot sample, for the determination of X, images coming from the same sample were processed and sets of E (containing the generalized interlayer distances) were extracted. These sets were united and X was computed from these united sets, thus a single value for X was obtained for each soot sample, regardless of multiple images per sample. The number of distances extracted from one 1024 x 1024 image w as between 30,000 and 60,000, and several micrographs were analyzed for each sample. Separate val­ues of t were calculated for each micrograph, then average values were produced for different samples. The conventional analysis of these images was also conducted for the sake of comparison, using a methodology based on the tech­nique published by Palotas e t al. [5]. The measured parameters of the fringe images w ere the fringe lengths (nm) - measured as pro­posed by Vander Wal et al. [11], i.e. the length of the skeleton o f the fringe - and interlayer distances (nm), measured as proposed by Palotas e t al. [5]. The results of this statistical analysis are pre­sented in Fig. 8. As it is apparent from Fig. 8, conventional statistical analysis pro­duced results based on w hich the reliable differentiation of similar structures is difficult. The values of the interlayer distances show insignificant differences between samples and there is no apparent trend in the evolution of fringe lengths. Apart from this, the inher­ently high standard deviation of the results provided by conven­tional image analysis makes the differentiation impossible, based on mean values alone. Fig. 8. Results o f a conventional statistic al image analysis meth o d applied to our samples. Some datap o in ts w e re slightly displaced to improve visibility. For exact values regarding HAB see Table 1. The results for X and t obtained by the proposed method are presented in Fig. 9. The trends of the two parameters are quite obvious as both de­creased with increasing HABs for all four samples, meaning tha t an increasing magnitude of order could be detected in the soot structures. Interestingly enough, an even stronger correlation was found between X and t. Based on the observation of X, the highest degree of structural evolution happened in the soots from the ethylene flames, especially a t low HABs. The rather small change in X measured in surrogate soots can be explained by the fact th a t these soots were found to be fairly well ordered even at low HABs (while still being amorphous), however a huge drop in t indicates some extent of simplification in the texture (notice the highly developed concentric rings in Fig. 7D). Benzene flames showed a slight decrease in X and a significant drop in t as HAB increased, also suggesting a rather evenly spaced population of fringes at all HABs, with fringes getting more elongated. Note tha t the continuous lines in the plots are ju s t guides to the eye, as soot was only sampled at a few HABs per fuel. The uncertainty bounds on t were computed as the standard deviation of the separate val­ues for each image. Based on both structural parameters, soot from the surrogate fuel showed the highest structural order. This structural order w as still evolving between HABs 1 mm and 5 mm. Ethylene showed the lowest extent of orderliness, however be­tween HABs 0.5 mm and 4 mm its X values were quickly decreas­ing, indicating a very abruptly transforming soot structure. Based on X values, the orderliness of soot produced in the second ethyl­ene flame was less than th a t produced in the first ethylene flame. This can be explained by the premixing system - the first flame was produced by a two-stage premixer, thus a more perfect mix­ing of the fuel and preliminary oxidation happened, allowing quicker soot structure evolution. Despite the few sampling locations, it is easy to see trends in the evolution of the structure of soot. Different soots from different fuels start to age at different heights. The data suggests tha t there are zones of accelerating and decelerating evolution for each fuel and th a t the location of these zones depend on both the fuel and operating conditions. As an idea for future work we note th a t by using the above described sensitive structural parameters, it is pos­sible to map these regions in flames. 7 12 8 P. Toth et al./Fuel 0.12 o.5--------------*--------------J--------------1--------------■--------------*--------------*- 0 2.5 5 7.5 10 12.5 15 Height above burner, mm Fig. 9. The mea sured values o f X and t for th e fuel samples. Lines are guides to th e eye. 4. Conclusions A sensitive image analysis technique has b een developed for the representative quantification of soot n anostructure orderliness uti­lizing several image processing techniques. The technique proposed here is robust, since the evaluated number of elements in the mea­sured datasets is large, and produces consistently reproducible parameters with a resolution high enough to show differences in soot structure. The m ethod has been applied to soot samples of var­ious fuels, collected from different burner setups a t different heights above burner in order to monitor their n anostructural evolution. In all cases, decreasing tendencies were found in the values of the distance deviation and junction parameters as the heights above burner increased, indicating increasing structural orderliness in the soot. The resolute and accurate monitoring of the nanostructural evolution of sootis possible by the described m ethodology provided th a t HRTEM images were representative of the sample. Acknowledgment This work w as partially sponsored by the TAMOP-4.2.1.B-10/2/ KONV-2010-0001 project with support by the European Union, co­financed by the European Social Fund. 99 (2012) 1-8 References [1] Laurent P, Braekman-Danheux C, Rouzaud JN. Microtextural study of cokes from hydropyrolysis o f coals. Fuel 1995;74(2):201-7. [2] Auguie D, Oberlin M, Oberlin A, Hyvernat P. Microtexture of mesophase spheres as studied by high re solution conventional transmission electron microscopy (CTEM). Carbon 1 980;18:337-46. [3] Ishiguro T, Takatori Y, Akihama K. Microstructure o f diesel soot particles probed by ele ctron microscopy: first observation of in n e r core and o u te r shell. Combust Flame 1 997;108:231-4. [4] Vander Wal RL, Yezerets A, Currier NW, Kim DH, Wang CM. HRTEM study of diesel soot collected from diesel particula te filters. Carbon 2007;4 5 :7 0 -7 . [5] Palotas AB, Rainey LC, Feldermann CJ, Sarofim AF, Sande JBV. Soot m orphology: an application of image analysis in high-re solution transmission electron microscopy. Microsc Res Tech 1 9 96;33:266-78. [6] SharmaA, Kyotani T, Tomita A. A n ew quan titativ e approach for microstructural analysis of coal ch a r using HRTEM images. Fuel 1999;78:1203-12. [7] Shim HS, Hurta RH, Yang NY. A methodology for analysis of 002 lattice fringe images and its application to combustion-derived carbons. Carbon 2 000;38:29-45. [8] He YangJ, Hui Cheng S, Wang X, Zhang Z, Rong Liu X, Hua Tang G. Quantitative analysis of micro stru c tu re of carbon materials by HRTEM. Trans Nonferrous Met Soc China 2006;1 6 :7 9 6 -8 0 3 . [9] Rouzaud JN, Clinard C. Quantitative high-re solution transmission electron microscopy: a promising tool for carbon materia ls characterization. Fuel Process Technol 2 0 0 6 ;7 7 -7 8 :2 2 9 -3 5 . [10] Palotas AB, Rainey LC, SarofimAF, SandeJBV. Whe re did th a t soot come from? Chemtech 1 998;28:24-30. [11] Vander Wal RL, Bryg VM, Hays MD. Fingerprinting soot (towa rds source identification): physical s tru c tu re and chemical composition. Aerosol Sci 2 010;41:108-17. [12] Alfe M, Apicella B, Barbella R, Rouzaud JN, Tregrossi A, Ciajolo A. S tru ctu re - p rope rty relationship in n an o stru c tu res of young and m a tu re soot in premixed flames. Proc Combust Inst 2 0 09;32:697-704. [13] Shaddix CR, Palotas AB, Megaridis CM, Choi MY, Yang NY. Soot graphitic order in laminar diffusion flames and a large-scale JP-8 pool fire. In t J Heat Mass Transfer 2 0 05;48:3604-14. [14] Zhu W, Miser DE, Chan WG, Hajaligol MR. HRTEM investigation o f some commercially available furnace ca rbon blacks. Carbon 2 004;42:1841-5. [15] Henson R, Reynolds W. Lattice p a ram e te r changes in irradiated graphite. Carbon 1 9 65;3:277-87. [16] Tidjani M, Lachter J, Kabre T, Bragg R. Structural d isorder induced in graphite by grinding. Carbon 1 9 86;24:447-9. [17] Aladekomo J, Bragg R. S tructural transformations induced in g raphite by grinding: analysis of 0 0 2 X-ray diffraction line profiles. Carbon 1990;28:897-906. [18] Aso H, Matsuoka K, Sharma A, Tomita A. Structural analysis o f PVC and PFA carbons p repared a t 5 0 0 -1 0 0 0 0C based on elemental composition, XRD, and HRTEM. Carbon 2 0 04;42:2963-73. [19] Echavarria CA, Sarofim AF, LightyJS, D'Anna A. Modeling and m ea su rem en ts of size distributions in premixed eth y len e and b enz ene flames. Proc Combust Inst 2 009;32:705-11. [20] Echavarria CA. Evolution of soot size d istribution during soot formation and soot ox id atio n -fragm en tatio n in premixed flames: experimental and modeling study. PhD Thesis. University o f Utah; 2010. [21] Dobbins RA, Megaridis CM. Morphology of flame-generated soot as d ete rmined by th e rmophoretic sampling. Langmuir 1 987;3:254-9. [22] Sauvola J, Pietikainen M. Adaptive do cum en t image binarization. Pattern Recogn 2 0 00;33:225-36. [23] Gatos B, Pratikakis I, P erantonis S. Adaptive degraded docum en t image binarization. Pattern Recogn 2006;3 9 :3 1 7 -2 7 . [24] Zhang TY, Suen CY. A fast parallel algorithm for thin n in g digital patterns. Image Process Comput Vision 1 984;27:236-9. [25] O'mara WC. Handbook of semiconductor silicon technology. William Andrew Publishing; 1990. [26] Alfe M, Apicella B, RouzaudJN, Tregrossi A, Ciajolo A. The effect of temperature on soot properties in premixed methane flames. Combust Flame 2010;157:1959-65. CHAPTER 3 A NOVEL FRAMEWORK FOR THE QUANTITATIVE ANALYSIS OF HIGH RESOLUTION TRANSMISSION ELECTRON MICROGRAPHS OF SOOT I. IMPROVED MEASUREMENT OF INTERLAYER SPACING Reprinted from Combustion and Flame, Vol. 160, Pal Toth, Arpad B. Palotas, Eric G. Eddings, Ross T. Whitaker, JoAnn S. Lighty, A novel framework for the quantitative anal- ysis of high resolution transmission electron micrographs of soot I. Improved measurement of interlayer spacing, Pages 909-919, Copyright 2013, with permission from Elsevier. 14 Combustion and Flame 160 (2013) 909-919 Combustion and Flame A novel framework for the quantitative analysis of high resolution transmission electron micrographs of soot I. Improved measurement of interlayer spacing Pal Totha,*\ Arpad B. Palotasb, Eric G. Eddingsa, Ross T. Whitakerc, JoAnn S. Lightya a Department of Chemical Engineering, University of Utah, 50 S. Central Campus Drive, Salt Lake City, UT 84112-9203, United States b Department of Combustion Technology and Thermal Energy, University of Miskolc, H3515 Miskolc-Egyetemvaros, Hungary c School of Computing, University of Utah, 3893 Warnock Engineering Building, Salt Lake City, UT 841112-9205, United States A R T I C L E I NF O A B S T R A C T The reliable and reproducible quantitative image analysis of digital micrographs from high resolution transmission electron microscopy (HRTEM) of soot has been an area of interest since the early nineties. Since the resolution of HRTEM images is usually sufficient to carry out structural measurements at the atomic level, the information obtained from these images is very valuable as it potentially yields insight into very specific soot oxidation processes; however, extracting physically meaningful, reliable, accurate and statistically robust data from HRTEM images is not an easy process. Data extraction is hindered by the presence of overlapping structures, varying focus, contrast, illumination levels and noise in the images. In this paper a novel image analysis framework is presented to address these issues and explore the possibility of the extraction of high-fidelity structural data from HRTEM soot images. Emphasis is on the analysis of images of mostly amorphous, poorly ordered soot structures, as these are the most diffi­cult to analyze. Published by Elsevier Inc. on behalf of The Combustion Institute. Article history: Received 27 October 2012 Received in revised form 8 January 2013 Accepted 8 January 2013 Available online 6 February 2013 Keywords: Soot Nanostructure HRTEM 1. Introduction The digital image analysis of high resolution transmission elec­tron microscopy (HRTEM) images of soot allows for the character­ization of the carbon structure at the atomic level. The nanoscale observation of soot structure is motivated by either an incentive to determine the source of soot pollution [1 ,2 ] or to better under­stand combustion processes [3 ]. A typical soot HRTEM micrograph shows an overlapping pattern of periodically occurring dark and bright lines of varying orientation and contrast. These lines are also called fringes and they are understood to be projections of graph­ene layers formed by phase contrast imaging principles. The mean­ing of dark and bright fringes are not obvious, as depending on the imaging conditions, they can either indicate carbon atoms or the spaces between them. Although not typical, the meaning of bright and dark fringes can be also dynamically reversed depending on imaging conditions and sample thickness. This phenomenon will be referred to as phase inversion in this paper. Also, it is practically important to satisfy some general conditions for the reliable mea­surement of geometric properties of fringe images: first, fringe contrast should be maximized and second, image regions where * Corresponding author. E-mail address: [e-mail] (P. Toth). fringe contrast changes abruptly should not be used for geometric analysis [4]. Image processing and analysis methods aiming to ex­tract structural information from HRTEM images of soot have been developed and applied since the mid-nineties [5 -1 0 ]. These meth­ods differ in details, however the basic procedure consists of the following steps: 1. The pre-filtering of the digital micrograph. This step usually consists of frequency filtering; i.e., the removal of unwanted frequencies in the image. Since carbon layers can only appear in a well-known frequency band determined by physically meaningful values of their interlayer spacing, bandpass fre­quency filtering is ideal for noise reduction in HRTEM soot images. 2. The detection of separate fringes. This step is usually an image binarization procedure, meaning that the initially gray scale image is transformed to a binary image in which fringes are indicated as 1 values and background is indicated by 0 values. Until lately, in most cases this transformation has been a global, non-adaptive binarization process; i.e., a single pixel intensity threshold has been set to determine whether a particular pixel belongs to a fringe or background. Authors have started to report results obtained by methods utilizing adaptive binariza-tion [10]. The outcome of this step is a binary image in which individual objects (fringes) can be detected and labeled. 0010-2180/$ - see front matter Published by Elsevier Inc. on behalf of The Combustion Institute. http://dx.doi.org/10.1016/j.combustflame.2013.01.002 15 3. The postprocessing of binary objects. In this step, the labeled binary fringes are processed further. The postprocessing tech­niques vary from author to author. Some use geometric criteria for selecting valid fringe candidates [5], some use skeletoniza­tion algorithms to reduce fringes to curves or line segments [6 -1 0 ] and some implement fringe separation/reconnection logic [7]. 4. The data extraction step. In this step geometric information is extracted from the postprocessed binary fringe objects. Geo­metric data can include fringe lengths [6 -1 0 ], fringe tortuosity [9,10], fringe separation (interlayer distance) [5 -1 0 ] and fringe orientation [5,8], among others. The obtained geometrical data is in the form ofvectors or sets of the mentioned properties. Each value in each set corresponds to a particular detected fringe in the analyzed micrograph. After extrac­tion, these sets can be statistically described by constructing prob­ability distribution functions (PDFs), specifying mean and standard deviation values and so forth. Structural order can be characterized by quantifying the symmetries and deviations in fringe orientation or fringe spacing values [6,8]. The common drawbacks of method­ologies that are based on the steps described above are the usage of subjectively set image processing parameters, results that are sen­sitive to these parameters and the insufficient amount of structural data extracted due to the oversimplification of images (the detec­tion and artificial separation of individual fringes). Because of all these reasons, it is understood that the quantitative characteriza­tion and differentiation of real soot HRTEM images can be a diffi­cult problem, especially when trying to differentiate poorly ordered, highly amorphous and/or only slightly different samples. In fact, the quantitative image analysis of amorphous soot samples has been considered an unsolved problem since the publication of the first papers in this area, especially when on aims to measure interlayer spacing statistics in soot structures (e.g., see insufficient number of data points in [11] or insignificant differences in statis­tics in [1 2 ]). Recently, we have developed an image analysis proce­dure with which distances can be measured between curved graphene layers as well, thus increasing the fidelity of the obtained distributions; however, this method does not avoid the subjectivity of binarization techniques and can only be used to obtain semi-quantitative descriptors [13]. 2. Materials and methods In this section a detailed description of the proposed image pro­cessing and analysis methodology is given along with the details and origins of the image sets used for the validation and verifica­tion of the proposed algorithm. 2.1. Image analysis procedure - overview The image processing methodology proposed in this work is completely different from the methods discussed in Section 1. The proposed method has been developed to address and over­come all the specific issues of the standard methods and it utilizes recent advances and ideas from the image analysis and signal pro­cessing literature. Similar algorithms have been applied in HRTEM crystallography[14,15]. While these published algorithms are sim­ilar, they are not applicable to images of more complex nanostruc­tures. The basic objectives of the proposed methodology are the following: 1. To extract as much structural data as possible from a single micrograph. A higher number of data points means more accu­rate and reliable statistical evaluation of the structure. Being 910 P. Toth et a l f Combustion i able to extract the highest amount of structural information from a single image has significant importance in cases where the availability of samples or micrographs is limited. Since typ­ical procedures for soot sampling deposit abundant quantities of soot particles on TEM grids, the scarcity of samples is usually not a problem when analyzing real soot. However, depending on the heterogeneity of the soot particles, acquiring and analyz­ing a number of micrographs that is sufficient for obtaining robust statistics can be a time-consuming process. In these cases, being able to extract as much information as possible from each micrograph not only increases the fidelity of the results, but has practical and economic importance as well. 2. To extract information only from reliable image regions in order to minimize measurement uncertainty and errors. 3. To measure structural parameters as accurately as possible. 4. To minimize the number and effect of subjectively chosen and set image processing parameters that are purely technical and do not hold any real physical meaning (e.g., threshold values, filter kernel sizes and parameters for fringe detection logic - see [5] or [11] for typically arising parameters). 5. To be able to handle typical aberrations and artifacts present in HRTEM micrographs; e.g., noise, inhomogeneous illumination and phase inversion. To achieve these objectives, one must consider a different mod­el for soot HRTEM images than the standard model. By image mod­el we refer to a mindset that is used to interpret these images. The standard model assumes that the micrographs contain identifiable or detectable objects; - i.e., imaged physical bodies (atoms or lay­ers of atoms) with well defined boundaries and the standard meth­ods aim to locate these bodies and their boundaries in the images. Such a model inherently leads to large amounts of eliminated data, as it is only interested in the detected bodies themselves (which are well defined and quantized) and their geometric properties. In­stead, our approach interprets the micrographs as continuous (at least at the image level) projections of patterns in electromagnetic fields1 and tries to analyze the patterns evolved in the projections. This image model is also more realistic from the quantum-physical point of view. Mathematically, images are interpreted as superposi­tions of two-dimensional patches of sinusoidal patterns with varying phase, frequency, amplitude and orientation corrupted by noise and contrast inhomogeneity. The proposed parameters to use for struc­tural characterization are exactly these (and the maps of these) sig­nal properties. It is easy to see that, following the proposed image model, there is no need for a detection step; - i.e., the labeling of well-defined binary fringes. Information can be extracted at the na­tive resolution of the image, meaning that every single pixel in the image will yield a set of structural parameters. The approach results in an abundant flow of information, which contributes to the statis­tical robustness of the extracted data. It is also possible to evaluate certain image regions based on their quality (signal strength) and only extract structural information from reliable regions. Since the approach is a spectral technique and is not based on pixel-level manipulations, sub-pixel accuracy can be achieved in the measure­ment of structural parameters. The upper limit of pointwise accuracy of these measurement is only set by the Nyquist sampling criterion. In the following, a detailed description of the methodology and ap­proach is given. 2.2. Image analysis procedure - dev elopment Three structural parameters are proposed in this study. These are the local orientation h, local wavelength k and local modulation Flame 160(2013) 909-919 1 In other words, signals. 16 P. Toth e t a!./Combustion and Flame 160 (2013) 909-919 strength i of the signal. These describe the sinusoidal pattern in a specific location in the image, but since the sinusoidal patterns are in fact the interpretations of fringes, the proposed structural prop­erties are simply the continuous generalized versions of the al­ready used and published fringe structural parameters. More precisely, local orientation is a generalized version of fringe orien­tation and local wavelength is the generalized version of the inter­layer spacing. The local modulation strength is a new parameter that describes and quantifies the reliability and local anisotropy of the neighborhood around a specific pixel. To extract these gen­eralized structural parameters, a localized frequency filtering ap­proach is applied based on filtering the images with a set of Gabor filters, also called a Gabor filter bank. Several similar algo­rithms have been proposed in the field of medical image process­ing and computer vision for applications like image registration [16,17]. The general approach is a hybrid spatial/frequency domain filtering method and it combines the advantages (sub-pixel resolu­tion, spectral representation and localization) of both [18]. Gabor filters are quadrature filters that can be fine-tuned to specific wavelengths, scales and orientations. Filtering (mathematically realized by two-dimensional convolution) the image with a Gabor filter yields a filtered image, also called a response. In this re­sponse, locations where the wavelength and orientation of the sinusoidal pattern in the original image were the closest to the wavelength and orientation of the Gabor filter yield the highest re­sponse value. Thus, by repeatedly filtering the image with a set of differently tuned Gabor filters and recording the responses, it is possible to extract local wavelength and orientation values at each image pixel. The term quadrature refers to the phase-independent response magnitudes of Gabor filtering. Since the phase of a sinu­soidal pattern is directly connected to its local intensity, it is easy to see that Gabor filtering analysis is immune to the phase inver­sion phenomenon mentioned in Section 1. In other words, it is not important whether dark or bright pixels indicate atoms or spaces between atoms - the responses will be the same in both dark and bright regions. Gabor filters also achieve the theoretical lowest localization uncertainty, thus they are ideal candidates for applications where one aims to minimize measurement uncertain­ties [19,20]. The analytical form of a two-dimensional Gabor filter in the spatial domain is the following: g(x, y ) = 1 exp 2 r 2 2 r2 x02 y 'M /2ftix'\ exp j = W (x ,y )S (x ,y ) (1) By definition, the two-dimensional (x, y) spatial Gabor filter is a superposition of a Gaussian window W(x, y) and a complex sinu­soidal S(x , y ). r x is the scale parameter of the Gaussian window in the x direction and r y is the scale parameter of the Gaussian window in the y direction. i is the complex unit vector and g x is the wavelength of the complex sinusoidal (in other words, the wavelength to which the filter is tuned to). The rotated coordinates x0 and y 0 can be expressed as x' = x cos (g „ )+ y sin(g„) y = y cos(gh) - x sin(g„) (2) (3) where g e is the orientation angle to which the filter is tuned to. gk and ge are not the same as the local wavelength k and orientation h of the image at a particular pixel, thus the different notation. The frequency domain form of the Gabor filter is a Gaussian surface. It can be written as the Fourier transform of g(x, y): G(u, v) = F [g (x ,y)| = exp r^ (2p + g ,u ')2 A 0 .5 (v < r y f ) (4) where u and v are the horizontal and vertical frequency coordinates, respectively. u 0 and v0 are rotated frequency coordinates and are ob­tained similarly to x0 and y 0. F denotes Fourier transformation with respect to a centrally shifted frequency domain (u, v). A quadrature Gabor filter pair along with its frequency domain representation is shown in Fig. 1 . The filter demonstrated in Fig. 1 is only a single member of the filter bank that is used to extract wavelengths and orientations from the image. The design of the whole filterbank is the most prac­tical to carry out in the frequency domain. The general idea behind the filter bank design is to cover and sample the relevant frequency bands of the image with Gabor filters. The filter bank is therefore a collection of sets of filter parameters ge, gk, r x, r y and by the design of the filter bank, one refers to the determination of these parame­ters. The determination of ge and gk is an easier process, since typ­ical HRTEM soot images contain orientations in the full range [ -p /2 , p/2[2 and wavelengths corresponding to a well defined fre­quency band derived from valid interlayer spacing values. One thing to note is that while the orientation scale is linear, the frequency (and therefore wavelength) scale is not, thus one must sample wave­lengths exponentially (lower frequencies/higher wavelengths must be sampled more densely). Also note that due to the direction ambi­guity of orientation angles, e values only cover a half circle (a range of p). In other words, an orientation angle value of e1 means exactly the same orientation as 01 ± lp, where l is an integer. One needs to specify a number of filters ne to cover the range of e and another number n to cover the range of k. Once these numbers are specified, the orienta­tion and wavelength values to sample can be calculated as p p g e; = - 2 + j nTe g k.k = kminq q = (5) (6) (7) where j and k are indices of positive integers between [1, nh] and [1, nk], respectively; 9j is the jth orientation, kk is the kth wavelength, kmin is the lowest wavelength (highest frequency) of interest and kmax is the highest wavelength (lowest frequency) of interest. The permutation of the sets of gh,j and gk,k obtained in this way maps out the filter center frequencies in the frequency domain. After these center frequencies are found, one still needs to specify the scale parameters r x and r y of each filter. The idea behind the spec­ification of the scale parameters is that the value of adjacent filters can be specified at a point equally far away from their center fre­quencies. Since there are n x nk possible permutations of sets of gej and gk,k, ne x values of r XtJ-,k and r y jk must be specified. By fix­ing the filter values qh and pk at intermittent points between their central frequencies corresponding to the overlapping values in the h (angular) and k (radial) directions in the frequency domain, it is apparent that the scale parameter values can be determined by the following relationships [2 1 ]: r x.i.k - I 2p (-1-------2_\ r y.j.k = (1 + q)(ln (pk))1/2 0.5 2p i ta^2% u -; 2(ln(pe)) (8) (9) Note that since the frequency responses of the Gabor filters are normalized, both ph and pk must be less than 1. If this condition is satisfied, the above equations yield real values for scale parameters. Here the symbols represent an interval that is open from the left and closed from the right. 1 2g k2 17 P. Toth et al/Combustion and Flame 160 (2013) 909-919 Fig. 1. From left to right: the real component of a Gabor filter in the spatial domain, the imaginary component (quadrature pair) of a Gabor filter in the spatial domain, the Fourier transformed Gabor filter in the frequency domain. The parameters used are: go = p/3, gk = 0.35 nm, ox = 0.25, Gy = 0.5. This filter would give the highest response at locations where fringes with a spacing of 0.35 nm oriented at p/3 radians are present in an image with a size of 200 x 200 pixels and a nm per pixel ratio of 5/150. The frequency domain filter is a Gaussian surface centered around frequencies (u0, v>). The Euclidean distance of this central point from the origo is exactly 1/g> Note that the filter shown here is enlarged for visualization purposes. Typical filters used for data extraction fit inside a 64 x 64 kernel window. It is also worth noting that the scale parameters are only functions of gk,k and not g ej, which is understandable since the orientation range is sampled linearly. Figure 2 demonstrates the design and construction of a typical Gabor filter bank for soot HRTEM image analysis. Note that while the filter bank essentially contains param­eter sets of go, gk, rx, Gy, to calculate these sets, the set nA, n e, p x, po must be determined. Therefore, for the balance of this paper the set nk, no, pk, po will be referred to as the filter bank parameters. To extract orientations and wavelengths, the HRTEM image is filtered with each filter in the constructed filter bank. The response values for each image pixel are recorded. Intuitively, one can deter­mine the orientation and wavelength at a specific pixel by identi­fying the filter which gave the highest response at that pixel. The orientation and wavelength of the particular pixel can be esti­mated as the orientation and wavelength to which the identified best filter is tuned to. While this method is intuitive and graphical, the resulting orientations and wavelengths will be heavily quan­tized depending on the sampling density of the filter bank. In order to obtain continuous and smooth orientation and wavelength maps, an interpolation method is used in this study. The basic idea of this technique is that the intermediate responses of the filter bank can be approximated by using the individual frequency re­sponses of the filters. The sum of individual filter responses, weighted by the discrete response values obtained by filtering, very closely approximates the continuous response surface one would obtain by filtering with a very densely spaced filter bank. In other words, the ideally high-resolution response surface can be approximated by a finite number of filter responses. The errors introduced by this approximation have been quantified by Perona in 1991 [22]. It has been shown that the error of the approximation quickly approaches zero as the number of filters increases. For our application, the error is practically negligible (approximately 1%) in o and k if one uses 15 or more filters per orientation and wave­length range. Mathematically, to find o and k at each pixel, one needs to find the global maximum of the interpolated response surface obtained by the following: Ro _ (R e (I0 go) + Im(10 go) \1/2| R(u, v ) = ^ 2 Go(u, v)Ro o=1 (umax, Vmax) _ argmax[R(u, v)| p 1 0 _ 2 + tan 1 (Vmax/umax) k _ K a x + vmax) ^ 1 R (umax; vmax) Rideal(umax; v max) (10) (11) (12) (13) (14) (15) Fig. 2. A demonstration of the design and construction of a Gabor filter bank. Left: HRTEM image of a real soot sample. Middle: the power spectrum (magnitude of Fourier transformation, plotted on a logarithmic scale) of the soot image. Circles denote the selected frequency band for further analysis. This frequency band corresponds to wavelengths between 0.3 and 0.6 nm. These values are chosen based on a priori knowledge of the physically meaningful interlayer spacing values for soot or can be selected automatically by locating the frequency band containing the highest spectral energy (as it is demonstrated by the bright ring around the origo). Right: gray scale values show the superimposed Gabor filters in the frequency domain. Crosses denote the center frequencies of the Gabor filters. The parameters of the filter bank are: no =10, nx = 10, po = 0.05, pk= 0.05. Note that the filter bank shown in this figure was not the filter bank that was used to extract the data that will be presented in Section 3.3. 912 18 where R(u, v) is the interpolated (continuous) response surface at a particular pixel (X0, y 0), o is the linear index of filters in the filter bank, Go is the filter frequency response of the oth filter as it is given by Eq. (4), Ro is the discrete response magnitude at a particular pixel (X0, y0) of the oth filter and (umax, vmax) is the frequency pair which maximizes R(u, v). I refers to the gray scale image itself, 0 is the symbol of convolution, Re() means the real part of its arguments and Im() means the imaginary part of its argument. The ideal max­imum filter response Rideai(u, v) refers to the theoretical maximum response at a location in the continuous interpolated response sur­face. This value is calculated by convolving an artificial image tem­plate containing a sinusoidal pattern oriented and spaced appropriately to a specific filter with maximum contrast and inten­sity (thus effectively a two-dimensional square-wave signal). This normalization scales the modulation strength between 0 and 1 and therefore makes it a good descriptor for comparing different images taken under similar imaging conditions. Note that it is pos­sible to scale the modulation strength by the maximum modulation strength found in the particular image, making it an intensity-invariant structural parameter and allowing for the comparison of images taken under different imaging conditions. To find the global maximum, a local direct search algorithm has been applied [23]. Global optimization algorithms are usually computationally inten­sive and difficult to implement, but the global maximum can be reliably found by local maximum search algorithms if the initial guess for the maximum is very close to the global maximum. Assuming a smooth and unimodal response surface, the center fre­quency of the filter that yielded the highest discrete response is a very good initial guess in this problem. The direct search algorithm has been implemented by using Matlab's fminsearch function. Fig­ure 3 demonstrates the maximum search step in the continuous re­sponse surface. Thus the complete algorithm can be summarized as follows: P. Toth et al./Combustion i 2.3. Soot sampling Graphite materials (Union Carbide) from three different fuels (anthracene, bifluorenyl and p-terphenyl) were previously synthe­sized by high temperature (3 0 0 0 °C) pyrolysis. The sample prepa­ration procedure has been reported earlier [1]. In brief, the samples were ground in an agate mortar and pestle and then ultra-sonically dispersed in ethanol. The suspension was deposited drop-wise onto a copper TEM grid coated with a lacey carbon film. HRTEM micrographs of the carbon particles were obtained using a transmission electron microscope operated at 2 0 0 keV. Amor­phous soot samples were taken using thermophoretic sampling methods previously described in [3 ]. In brief, benzene was burned in air under fuel rich conditions in a premixed flat flame burner that was initially developed for studying aliphatic and aromatic flames. This system consisted of a stainless steel chamber where fuel and air were properly mixed prior to entering the burner. The flame was stabilized over a tube bundle and was shielded from atmospheric interference using a nitrogen shroud. A metallic mesh placed 3.5 cm above the burner surface stabilized and distributed the flame uniformly across the burner. Samples for HRTEM analysis were taken at different heights above the burner surface (HAB) using a thermophoretic probe commonly referred to as a frog ton­gue. A TEM grid holder was attached to a piston and compressed air at 60 psig was used to quickly insert and remove the TEM grid from the flame. Multiple insertions were necessary to get a repre­sentative soot sample on the grid. The grid was oriented with the face parallel to the gas flow, so that the disturbance of the flame was minimal. Soot deposits on the grid due to the thermophoretic gradient between the cold grid and the hot flame, which also al­lows freezing of some heterogeneous reactions, as well as avoiding changes in the soot morphology after the particles have impacted upon the cold surface. Flame 160 (2 0 1 3 )9 0 9 -9 1 9 913 1. Load the image. 2. Specify the wavelength band of interest either manually or automatically based on the power spectrum. 3. Calculate filter parameters by using Eqs. (5 )-(9 ) and con­struct the filter bank. 4. Filter the image with each filter in the filter bank and record the discrete responses at each pixel. 5. At each pixel, calculate the frequencies that maximize the continuous interpolated response surface using optimiza­tion and store these frequencies using Eqs. (1 0 )-(1 2 ). 6. Calculate structural parameters h, k and 1 at each pixel using Eqs. (1 3 )-(1 5 ) and store. While filtering can be realized by numerical convolution using computationally inexpensive Fast Fourier Transforms (FFTs), the whole algorithm represents a huge computational load due to steps 5 and 6. To overcome this problem and keep the computa­tion times under reasonable limits, steps 5 and 6 have been implemented as parallel algorithms. Under a Matlab environment using a modern quad-core computer computations take 10­60 min, depending on the image and selected filter bank. Process­ing the typical HRTEM image shown in Fig. 2 took 17 min. If in a specific application computation time is an issue, the algorithm can be further optimized by only using discrete filter responses in the close neighborhood of the filter with the maximum dis­crete response to compute R(u , v) in the frequency domain, man­ually specifying regions of interests (ROIs) or outsourcing parts of the computation to graphical processing units (GPUs) that are commonly found in modern personal computers. Figure 4 demon­strates typical extracted fields of the three derived structural parameters. 3. Results and discussion To demonstrate the capabilities ofthe proposed algorithm, tests with both artificial and real images have been carried out. Testing with artificial images is important, since it can achieve what test­ing with real data cannot: verification of the method by compari­son with accurately known ground truth. To demonstrate that the proposed method is capable of reproducing already published results - with higher resolution and accuracy - results from a pa­per presenting standard structural analysis of different graphites [1 ] have been compared to results obtained by the method pre­sented here. In addition to demonstrate the methods capability of analyzing and differentiating poorly ordered soot structures, re­sults for images of immature soot obtained by standard image analysis methods and the method proposed here are compared [5 ]. 3.1. Artificial images Since the method proposed here is practically a signal analysis routine with a number of features specifically designed for soot HRTEM image analysis, artificial images can be constructed to val­idate the algorithm and quantify its accuracy. Images showing a sinusoidal pattern with varying orientation and wavelength have been generated by taking the sine of an underlying polynomial phase surface. The theoretical wavelength can be calculated as the inverse of the phase gradient magnitude, while the theoretical orientation can be calculated as the orientation angle of the gradi­ent. After generating the images, they have been artificially cor­rupted by zero-mean Gaussian noise. The standard deviation of the noise has been increased in steps in order to evaluate the algo- 19 914 P. Toth et al./Combustion and Flame 160 (2013) 909-919 280 300 320 340 -2 0 2 2 2.5 3 x U u Fig. 3. Two different locations in the real HRTEM soot image shown in full size in Fig. 2 are shown in the first column. WhiteX symbols indicate the analyzed pixel. The second column shows the reconstructed continuous response surfaces for both pixels. Gray dots indicate the center frequencies of the filters. The grayscale values show the values of the response surfaces. The third column shows the results of the maximum search algorithm. These plots show the vicinities of the maximum to improve visibility. The initial guesses for the maximum locations are indicated by the / symbol. The refined maximum found is shown by the ★ symbol. Notice that the angular coordinate of the maximum corresponds to a direction normal to the local orientation (see Eq. (13)) and the radial coordinate corresponds to the reciprocal local wavelength (see Eq. (14)). 250 300 350 400 450 250 300 350 400 450 250 300 350 400 450 X X X Fig. 4 . Extracted modulation strength (i), orientation (h) and wavelength (k) maps from a detail of the real soot HRTEM image shown in Fig. 2, from left to right, respectively. In the first image, hues of orange indicate the value of the modulation strength parameter - the more orange a pixel is, the highest its i parameter. Notice how the map of i highlights the best defined and most anisotropic regions. These regions correspond to well-imaged crystalline regions developing short-range order. In the second image, the orientation angles of each pixel are plotted as blue line segments. The orientation of these lines indicate h. Note that h is only accurate and meaningful in unimodally oriented areas - these are the areas with high i values. Other areas return a h value corresponding to the spectrally strongest orientation angle. The third image shows a color coded plot of k. The values indicate the distances of fringes at each pixel in the direction that gave the highest response value. If there is only one well defined direction, the value of k is the generalized interlayer spacing in nanometers at the particular pixel. The visual interpretation of the map of k of a real soot image is not easy. A more easily understandable representation will be presented in Fig. 5, Section 3.1. 20 P. Toth e t al./Combustion and Flame 160 (2013) 909-919 Fig. 5. Top left: an artificial image constructed based on an underlying polynomial phase surface. The image has been corrupted by zero-mean Gaussian noise. The noise level of this particular image is -2.5 decibels. Top middle: ground truth local wavelength map (k) of the artificial image computed as the gradient magnitude of the noise-free underlying phase. Top right: ground truth local orientation angle map (h) of the artificial image. Bottom left: extracted modulation strength map (i). Higher values mean more reliable locations in the image in terms of wavelength and orientation measurement. Bottom middle: measured local wavelength map. Bottom right: measured local orientation angle map. Note that the upper right corner is inaccurate in both the measured k and h maps, but the values of i corresponding to this region are also low, indicating less reliable estimation. 915 px and p0 Noise level, dB Fig. 6. Results of error quantification tests. Top left: the relative mean error in k as a function of n and n> Denser sampling in the frequency domain yields lower errors, although it is worth noting that the accuracy of estimations of k has been found insensitive to nh in the tested parameter range. Top right: the relative mean error in h as a function of nh and nk. Unlike k, estimations of h have been found to be sensitive to both nh and nk in the tested parameter range. Bottom left: the effect of ph and pk on the relative mean error in h and k. It has been found that there are ideal values for these parameters with which the error is minimal. Overall, the effect of the two p parameters has not effected error values significantly. Bottom right: the effect of noise levels on the accuracy of estimations. The mean error in k stayed under 1%, while the mean error in h stayed under 2% for typical noise levels. In the high noise region, noise effected errors the most among the evaluated conditions. Note that in the range of typical noise levels the algorithm is practically immune to noise, as the errors introduced by noise are negligible. Excessive noise levels resulted in inaccurate estimations. It is not advised to use the proposed algorithm for images with noise levels above 5 decibels. rithm's capability to handle image noise. Noise levels are given in where LN is the noise level in decibels, r N is the standard deviation decibels calculated by the following relationship: of the Gaussian noise and max(/) is the maximum signal strength r r (the highest pixel value in the image). A logarithmic noise range be- Ln = 10log ----- (16) tween - 2 0 and +10 decibels has been evaluated. Soot HRTEM N [max(/)J ' ' 21 916 P. Toth et al./Combustion and Flame 160 (2013) 909-919 Bifluorenyl Anthracene p-terphenyl 0.32 0.34 0.36 0.38 0.32 0.34 0.36 0.38 0.32 0.34 0.36 0.38 interlayer spacing, nm interlayer spacing, nm interlayer spacing, nm Fig. 7. Results of graphite structure analysis. Top left: HRTEM image of graphite originating from burning bifluorenyl. Top middle: HRTEM image of graphite originating from burning anthracene. Top right: HRTEM image of graphite originating from burning p-terphenyl. Bottom left: interlayer spacing PDF of bifluorenyl graphite. Bottom middle: interlayer spacing PDF of anthracene graphite. Bottom right: interlayer spacing PDF of p-terphenyl graphite. In the PDF graphs the results from Ref. [1]are compared with the results obtained by our proposed algorithm. images obtained by modern microscopes have noise levels typically between - 1 5 and - 5 decibels. Figure 5 shows a generated artificial image along with its wavelength and orientation fields and the re­sults of analysis carried out by the algorithm proposed here. A series of tests have been conducted to map out the errors and sensitivities of the method. Since the wavelength ranges are approximated based on a priori physical knowledge and the filter scales are dependent on filter spacing, the proposed algorithm has only four independent parameters: nh, nk, ph and pk (see Sec­tion 2 .1 ). It will be shown that within a reasonably wide range, the results are accurate and insensitive to variations in these parameters. The most important parameters are the number of fil­ters in both the k and h direction, as these determine the accuracy of the estimation of the continuous response surface by discrete fil­ter responses (see Section 2 .1 ) at a particular pixel. A set of tests have been carried out by keeping ph and pk constant (0.9) and try­ing all combinations of nh and nk values between 10 and 30. The ef­fect of changing p h and pk has also been evaluated by keeping nh and nk constant (both 20) and varying both p h and p k between 0.65 and 0.97 in steps of 0.02. Since the two p values effect the shape (aspect ratio) of each filter simultaneously, there is no real need to evaluate their different combinations. A noise level of - 2 .5 decibel has been used to test the algorithms sensitivity to in­put parameters. To evaluate the algorithms sensitivity to noise, a series of images have been generated by using the same underlying pattern but with increasing noise levels. These images have been processed and analyzed by the following filter bank parameters: n e = 20, nk = 20, p e = 0.9, pk = 0.9. The algorithm has been evaluated based on error norms. Each resulting h and k map has been com­pared to the theoretical maps and their relative errors have been computed. A single scalar error measure has been computed for each test for both h and k by taking the mean of the error matrices. Figure 6 shows the results of these tests. To summarize the results ofvalidation tests it can be stated that under typical image conditions, the mean error in both h and k is under 1%. If the frequency domain is sufficiently densely sampled by filters the results are insensitive to the particular choices in nh and nk. Our results agree well with the reportings of Perona [2 2 ]; namely, it is advised to use at least 15 filters for both h and k. It is also worth mentioning that for most soot images, p h and p k should be set so that the filter bank parameters result in circular filters; i.e., the aspect ratio r x/ r y should be approximately 1 for each filter, as this will result in the lowest error. 3.2. Analysis of graphite The results presented in this section are reproductions ofthe re­sults of Palotas et al. published in 1998 [1]. The original study fo­cused on interlayer measurements in HRTEM images of graphites from different sources. The results of interlayer spacing measure­ments for graphites from bifluorenyl, anthracene and p-terphenyl (among others) were presented and the methodology for generat­ing these graphites can be found in Section 2.3 or in more detail in the original publication. The original image processing algorithm used in [1] is a standard routine based on frequency filtering, glo­bal binarization and fringe detection. The interlayer distances were measured as the distances between parallel fringe pairs, with ori­entations being made equal to the mean of the two orientation val­ues of the fringe pair. Details of this method can be found in [5]. The interlayer spacing values in three graphite HRTEM images have been measured by the algorithm proposed in this paper. Since graphitic structures contain long, parallel carbon layers, measuring 22 P. Toth e t al./Combustion and Flame 160 (2013) 909-919 Statistics of the measurements of graphite interlayer spacing compared to the statistics of Palotas et al. [1]. Notice the significantly richer datasets obtained by the proposed method. Sample Mean (nm) Stddev (nm) Mode (nm) Comp. time Results from Ref. [ 1 ] Bifluorenyl Anthracene p-Terphenyl Proposed algorithm Bifluorenyl Anthracene p-Terphenyl 557 451 468 170,274 74,978 15,320 0.3355 0.3404 0.3481 0.3357 0.3401 0.3482 0.0049 0.0079 0.0093 0.0063 0.0066 0.0101 0.3322 0.3316 0.3479 0.3316 0.334 0.3443 3s 3s 3s 5 min 5 min 5 min 917 Table 1 # Of data Fig. 8. Top row: HRTEM images of benzene soot sampled at different HABs (5 mm, 10 mm and 15 mm from left to right, respectively). Bottom left: interlayer spacing PDFs obtained by the proposed method and a standard method [5]. Bottom middle: orientation angle PDFs obtained by the proposed method and a standard method [5]. Bottom right: modulation strength PDFs obtained by the proposed method. Notice the high-fidelity, high-resolution interlayer PDFs in the case of the proposed method. The rich datasets allowed for setting the bins by steps of 0.0035 nm in the case of interlayer distances and pi/100 radians in the case of orientation angles. For the standard algorithm, ten evenly spaced bins have been defined in the shown ranges. the orientation maps was not of practical importance. This section aims to demonstrate that the generalized interlayer spacing values (k) obtained by the method presented in this paper yield practi­cally the same distribution as the interlayer spacing values ob­tained by the standard basic structural unit analysis algorithm. The following filter bank parameters have been used to obtain the data: n k = 1 5 , ne = 30, pk = 0.99, pe = 0.9. The values of the ob­tained maps of i spanned a range of 0 -0 .2 . Pixels with i values above 0.06 have been used as reliable pixels for the measurement of k. Figure 7 shows the analyzed graphite images along with the measured distributions of k compared to the interlayer distribu­tions first presented in [1]. The information obtained by using the method proposed here implies almost exactly the same physical structures as the infor­mation obtained by previous authors. The only difference is in the accuracy and robustness of estimation. Table 1 shows the ex­tracted statistics of the measurements. Our method produced three orders of magnitude richer datasets on average compared to the re­sults of Palotas et al. [1]. Notice that the number of reliable data points was lower for p-terphenyl graphite than for the other two. It is easy to understand Table 2 The statistical evaluation of interlayer distance data obtained from the benzene images. This table compares the results of the proposed method with a standard algorithm described in [5]. Notice the lower standard deviations obtained by our method and the trend in the modes of the interlayer distances. There is no apparent trend in the results obtained by the standard method. On average, the datasets obtained by the proposed methods contained three orders of magnitude more information at a cost in computation time, which was approximately a 100 times longer for the proposed algorithm. HAB # Of data Mean (nm) Stddev (nm) Mode (nm) Comp. time Algorithm from Ref. [5] 5 mm 31 0.5273 0.0531 0.5222 12 s 10mm 56 0.4427 0.0598 0.4056 13 s 15 mm 104 0.4377 0.0751 0.4056 13 s Proposed algorithm 5 mm 78,259 0.4071 0.0184 0.402 14 min 10 mm 66,739 0.4001 0.0209 0.3949 15 min 15mm 168,504 0.4103 0.0243 0.3843 17 min why by looking at the image of p-terphenyl graphite - there were less non-overlapping regions with good contrast. Table 1 demon­strates very good agreement between the extracted mean and 23 918 P. Toth et al./Combustion and Flame 160 (2013) 909-919 450 500 550 600 0.33 0.34 0.35 0.36 interlayer spacing, nm 350 400 450 0.3 0.4 0.5 interlayer spacing, nm Fig. 9. The effect of different methods for incorporating i in the obtained interlayer spacing distributions. In the left column, micrographs of the bifluorenyl graphite (top) and benzene soot (bottom) are shown with isocontours of i overlain. The three shades of gray of the contours correspond to the thresholds shown in the legends of the right column. In the right column, resulting interlayer spacing distributions are shown for the bifluorenyl graphite sample (top) and benzene soot (bottom). Three distributions are shown obtained by thresholding i by using the threshold values shown in the legends and a fourth one (denoted by red crosses) obtained by weighing k with i. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) mode values of the interlayer spacing distributions. The slight dif­ferences in the obtained standard deviation values can be ex­plained by the inherent volatile nature of variance statistics; in other words, the number of data points extracted by the conven­tional technique was insufficient for the accurate estimation of the standard deviation of the distributions. Supposedly, introduc­ing additional samples to the conventional measurement would converge the standard deviation estimates toward those obtained by the proposed method. 3.3. Analysis o f amorphous soot To demonstrate the proposed methods applicability to the anal­ysis of highly amorphous soot samples, HRTEM images of benzene soot sampled at different HABs have been processed. Typically, the compaction of the layered structure is expected during the matu­ration or oxidation of black carbon or soot particles [12,24], result­ing in slight shifts in the interlayer spacing distributions towards lower values. Benzene soot images have been processed with the following filter bank parameters: nx = 1 5 , no = 30, p x = 0.97, po = 0.97. The modulation strength values ranged from 0 to 0.6 - a reliability threshold of 0.15 has been selected. Figure 8 shows the results of the measurements. It is clear from Fig. 8 that the proposed method was able to identify a slight shift in interlayer distances towards shorter dis­tances. This shift corresponds to the maturation and compaction of soot nanostructure. Interestingly, the standard deviation of the interlayer spacing data also increased as HAB/residence time in­creased. The standard method failed to identify any trends. It is worth noting that in the case of completely amorphous structure, the standard method produced significantly overestimated inter­layer distances, which was caused by the lack of truly parallel fringes. As the number of data points extracted by the standard method increased, its results started to get closer to the results ob­tained by the new method. Orientation angle PDFs obtained by the standard method showed rough consistency with the results of the proposed method. While PDFs of k and o can provide physical information on the structure, the PDFs of i can be used as finger­printing tools. The PDFs of i are indicative of the overall crystalline order in the structure - for similar images, the more they are skewed towards higher values, the higher the structural order is. PDFs of i show a trend of increasing orderliness in benzene soots as HAB increases, as would be expected (notice the increasing frac­tion of pixels above the threshold i value - the locations of the peaks are irrelevant, since low i values represent image back­ground and unreliable image regions; their corresponding pixels are excluded from further calculations). The statistical summary of the results is shown in Table 2 . 3.4. Choosing the threshold o f i Strictly speaking, the modulation strength parameter i is a descriptor without exact physical meaning. i is basically a normal­ized convolution product of an image detail representing a carbon sub-structure with a best-fit filter kernel modeling a structural primitive. In other words, i is a scalar that represents the similar­ity between the local structure and the image ofan idealized build­ing block of carbon layers. i therefore carries no absolute quantitative information,3 but is instead used to classify image re­gions based on their overall quality. It is easy to understand that ac­tual values of i depend not only on the local contrast, but on the local morphology of the carbon layers as well. In this section a num­ber of simple strategies are proposed on how to use i in order to only include high-fidelity information in the obtained interlayer spacing and orientation distributions. 3 i can be arbitrarily scaled - only its distribution is important. 24 P. Toth et al./Combustion and Flame 160 (2013) 909 -9 1 9 919 The two simplest ways to incorporate i a s a reliability measure are thresholding (used in Sections 3 .2 and 3 .3 ) and weighing. As discussed above, thresholding is a procedure that defines a critical value of i , under which pixels are considered as unreliable and are excluded from further calculations. The threshold values presented in Sections 3 .2 and 3 .3 were set manually, by overlaying the scalar fields of i on the micrographs and choosing thresholds that de­fined and enclosed visually appealing regions. Interestingly, the manually chosen thresholds roughly corresponded to 30% of the highest modulation strength values observed in the analyzed micrographs, although it is expected that this rule of thumb does not generalize very well. Another method of determining a feasible threshold value is to use a priori information regarding the interlayer statistics and opti­mize the threshold of i , such that the obtained interlayer spacing distributions best approximate the expected outcome. This route is only recommended when the a priori information can be regarded as highly accurate; e.g., in the case of graphite. Obviously, a single value should be chosen for a complete set of images, so that the re­sults are comparable. Instead of setting strict rules on the threshold of i , the modula­tion strength can be used as a soft measure of reliability as well. In a Bayesian sense, each interlayer spacing data point (observation) can be given a reliability factor or weight - the value of i at the same location. Instead of filtering out a number of less reliable pix­els and building a distribution of the remaining, one can build weighted distributions, to which less reliable pixels contribute less than more reliable ones. Figure 9 illustrates the strategies dis­cussed in this section. Since graphite generally has elongated and straight carbon lay­ers, it is understandable that distributions of i extracted from micrographs of graphite are mostly bimodal - one peak represents the background, where i is practically zero and the other the struc­ture, where i is closely maximum. Therefore, in the case of graph­ite samples it is not surprising that choosing any value of i as a threshold that is above the value representing the background suf­fices. In the case of benzene soot, changing the threshold has a greater effect (although not significantly, provided that at least the background is excluded by thresholding), however there are no exact methods of determining the threshold value. Since many factors have an effect on the actual values of i , using the visual (manual) approach discussed above is recommended. Naturally, for similar structures, the same threshold should be used to facili­tate semi-quantitative comparison. Notably, in most practical situ­ations, the modes and means of the interlayer spacing distributions are insensitive to the threshold of i , provided that the threshold is within a ‘reasonable' range - the only effected parameters are usu­ally the variances. If one aims to avoid subjective thresholding completely, using weighted distributions is recommended. In all the cases studied in this paper, the weighted distributions were very close to the ones utilizing manually set thresholds. 4. Conclusion A novel image processing framework for the analysis of amor­phous soot HRTEM images has been designed, developed, tested and evaluated in this paper. The proposed method is completely different from all previously published approaches and is capable of extracting approximately 1000 times more structural informa­tion from the same amount of micrographs than standard meth­ods. The method has been purposely developed to be able to extract the most accurate and reliable structural information from soot HRTEM images. To achieve this objective, the method has been designed to be practically immune to image noise, phase inversion phenomena and to yield the lowest localization uncer­tainty that is theoretically possible. Unlike standard methods, the proposed algorithm provides structural information at native im­age resolution; i.e., every image pixel yields a set of structural parameters. The method has been tested on artificially created images, real micrographs of graphite and real micrographs of amorphous soot. It has been found that the method provides infor­mation consistent with results obtainable by standard algorithms for graphitic structures; however, it is still able to provide high-fidelity data in cases where most standard techniques fail - specif­ically in the case of amorphous soot samples. As a demonstration of the technique's capabilities, we have applied the methodology to laboratory collected benzene soot. Our analysis revealed an in­crease in the degree of order (compactness) in soot structure with oxidation (maturation). The finding is consistent with literature data, however the supporting data is orders of magnitude more robust. Acknowledgments This work was partially sponsored by the TAM0P-4.2.1.B-10/2/ K0NV-2010-0001 Project with support by the European Union, co­financed by the European Social Fund. The authors would like to thank Carlos Andres Echavarria at the University of Utah for pro­viding some HRTEM images. References [1] A.B. Palotas, L.C. Rainey, A.F. Sarofim, J.B.V. Sande, R.C. Flagan, Chemtech 28 (1998) 24-30. [2] L.C. Rainey, A.B. Palotas, A.F. Sarofim, J.B.V. Sande, Applied Occupational and Environmental Hygiene 11 (1996) 777-781. [3] C.A. Echavarria, Evolution of Soot Size Distribution During Soot Formation and Soot Oxidation-Fragmentation in Premixed Flames: Experimental and Modeling Study, Ph.D. Thesis, University of Utah, 2010. [4] M.J. Hytch, T. Plamann, Ultramicroscopy 87 (2001) 199-212. [5] A.B. Palotas, L.C. Rainey, C.J. Feldermann, A.F. Sarofim, J.B.V. Sande, Microscopy Research and Technique 33 (1996) 266-278. [6] J. Yang, S. Cheng, X. Wang, Z. Zhang, X. Liu, G. Tang, Transaction of Nonferrous Metals Society of China 16 (2006) 796-803. [7] A. Sharma, T. Kyotani, A. Tomita, Fuel 78 (1999) 1203-1212. [8] H.S. Shim, R.H. Hurt, N.Y.C. Yang, Carbon 38 (2000) 29-45. [9] J.N. Rouzaud, C. Clinard, Fuel Processing Technology 77-78 (2002) 229-235. [10] K. Yehliu, R.L.V. der Wal, A.L. Boehman, Combustion and Flame 158 (2011) 1837-1851. [11] K. Yehliu, R.L.V. der Wal, A.L. Boehman, Carbon 49 (2011) 4256-4268. [12] C.R. Shaddix, A.B. Palotas, C.M. Megaridis, M.Y. Choi, N.Y.C. Yang, International Journal of Heat and Mass Transfer 48 (2005) 3604-3614. [13] P. Toth, A.B. Palotas, J. Lighty, C.A. Echavarria, Fuel 99 (2012) 1-8. [14] M.J. Hytch, Scanning Microscopy 11 (1997) 53-66. [15] M.J. Hytch, E. Snoeck, R. Kilaas, Ultramicroscopy 74 (1998) 131-146. [16] M.I. Elbakary, M.K. Sundareshan, Pattern Recognition Letters 26 (2005) 2154­2173. [17] M.I. Elbakary, M.K. Sundareshan, Image and Vision Computing 25 (2007) 663­670. [18] N. Bonnet, Micron 35 (2004) 635-653. [19] D. Gabor, Journal of the Institution of Electrical Engineers III: Radio and Communication Engineering 93 (1946) 429-441. [20] G.H. Granlund, H. Knutsson, Signal Processing for Computer Vision, first ed., Kluwer Academic Publishers, Dordrecht, Netherlands, 1994. [21] J. Ilonen, J.K. Kamarainen, H. Kalviainen, Efficient Computation of Gabor Features, Technical Report, Lappeenranta University of Technology, Finland, 2005. [22] P. Perona, in: IEEE Computer Society Conference on Computer Vision and Pattern Recognition, IEEE Computer Society 1991, pp. 222-227. [23] J.C. Lagarias,J.A. Reeds, M.H. Wright, P.E. Wright, SIAMJournal of Optimization 9 (1998) 112-147. [24] X. Zhang, A. Dukhan, I. Kantorovich, E. Bar-Ziv, A.F. Sarofim, in: Twenty-sixth Symposium (International) on Combustion, The Combustion Institute, 1996 pp. 3111-3118. CHAPTER 4 A NOVEL FRAMEWORK FOR THE QUANTITATIVE ANALYSIS OF HIGH RESOLUTION TRANSMISSION ELECTRON MICROGRAPHS OF SOOT II. ROBUST MULTISCALE NANOSTRUCTURE QUANTIFICATION Reprinted from Combustion and Flame, Vol. 160, Pal Toth, Arpad B. Palotas, Eric G. Eddings, Ross T. Whitaker, JoAnn S. Lighty, A novel framework for the quantitative analysis of high resolution transmission electron micrographs of soot II. Robust multiscale nanostructure quantification, Pages 920-932, Copyright 2013, with permission from Elsevier. 26 Combustion and Flame 160 (2013) 920-932 Combustion and Flame A novel framework for the quantitative analysis of high resolution transmission electron micrographs of soot II. Robust multiscale nanostructure quantification Pal Totha,*\ Arpad B. Palotasb, Eric G. Eddingsa, Ross T. Whitakerc, Joann S. Lightya a Department o f Chemical Engineering, University o f Utah, 50 S. Central Campus Drive, Salt Lake City, UT 84112-9203, United States b Department o f Combustion Technology and Thermal Energy, University o f Miskolc, H3515 Miskolc-Egyetemvaros, Hungary c School o f Computing, University o f Utah, 3893 Warnock Engineering Building, Salt Lake City, UT 841112-9205, United States A R T I C L E I NFO A B S T R A C T The quantitative characterization of mostly amorphous soot structures is a difficult problem. High reso­lution electron microscopy is a tool capable of providing structural information related to the crystalline order in soot; however, well-defined and exhaustive structural parameters are needed for quantification. The typical observable field of view and the insufficient amount of structural information extractable from a single electron micrograph pose another problem in obtaining reliable statistical description. This paper has two objectives: first, to show that the already developed and published structural descriptors can be united by introducing a general model for the characterization of molecular order and second, to extend the general filtering approach presented in Part I of this study to allow for the efficient extraction of such general parameters. The computational background is described with automatic, real-time future applications in mind. Published by Elsevier Inc. on behalf of The Combustion Institute. Article history: Available online 12 February 2013 Keywords: Soot Nanostructure HRTEM 1. Introduction Graphite, an allotrope of carbon has a hexagonal lattice struc­ture with an atom spacing of 0.142 nm and interlayer spacing of 0.335 nm. The layers o f the hexagonal carbon structure are referred to as graphene layers. These graphene layers can be imaged and quantified by using high resolution transmission electron micros­copy (HRTEM). In a HRTEM micrograph, graphene layers appear as dark or bright linear patterns, also called fringes. Soot is a product o f the pyrolysis of carbonaceous materials and is generally considered as amorphous carbon [1]. Despite being amorphous, in most cases soot nanostructures show some degree of crystalline order typically in the form of graphite microcrystals (mesophasic crystalline units or clusters exhibiting short-range or­der in the form of parallel graphene layers, also called stacks), par­tial fullerenic (graphene layers in a concentrically symmetric, onion-like structure) or partial graphitic (longer range parallelity of layers) order [2-4]. Soot nanostructure is effected by combus­tion conditions, the thermal environment [3] and also by the com­busted fuel type [5-7]. Further, in this paper we use the term ‘nanostructure' to refer to the degree of crystalline order in soot. Soot n anostructure quantification has been an area of interest since * Corresponding author. E-mail address: [e-mail] (P. Toth). the early nineties [8]. Being able to quantitatively describe soot nanostructure is important from two perspectives. First, it has been hypothesized th a t soot nanostructure has an effect on soot reactivity. Correlations are available showing the interdependence of oxidation kinetics and nanostructure [2,3,6­11]. It is understood th a t the reactivity of edge site carbon atoms (at the edges of graphene layers) is higher than th a t of carbon atoms in the basal plane (surrounded by other carbon atoms in the graphene layer) [6,12]. Since the transformation of combus­tion- derived carbon has b een found to ultimately lead towards gra­phitic structure as maturation proceeds (e.g., in the case of biodiesel soot [13] or in large-scale pool fires [14]), the appearance and growth of graphitic clusters can eventually cause passivation. The ordering of the initially amorphous structure also results in the compaction of the carbon structure [15] and decreasing macro-scopical surface area for surface reactions. These processes may have an effect of the reactivity of soot formed in high-temperature flames as well. Similar studies are available on coal and coal char nanostructure [16,17] and carbon blacks [18]. Second, since soot nanostructure is affected by the fuel source [5-7], quantitative descriptors may be used as forensic tools [19­21]. Nanostructural parameters from HRTEM micrographs have also been combined w ith chemical composition information to ob­tain a dataset more suitable for source identification [21]. Along with X-ray diffractometry (XRD), digital processing and analysis of digitized HRTEM micrographs is a technique capable 0010-2180/$ - see front m a tte r Published by Elsevier Inc. on b eha lf o f The Combustion Institute. http://dx.doi.org/10.1016/j.combustflame.2013.01.003 27 P. Toth et al./Combustion and Flame 160 (2013) 920-932 921 of providing quantitative soot nanostructure information. Since XRD structural parameters are derived based on assumed crystal models and are averaging in nature (references to published XRD analysis techniques are given in [22]), and due to the fact that HRTEM allows for the direct obser