Topologyfaithful nonparametric estimation and tracking of bulk interface networks
Abstract
The main focus of this paper is a nonparametric filtering technique for the estimation of interface geometry in bulk materials obtainable from modern imaging measurements. The filtering methodology relies on an assumed hierarchy of topological features present in a typical interface network, such as foam interfaces and grain boundary networks in polycrystalline materials. Each type of topological feature is treated in order of rank in the hierarchy, with the lowerlevel feature being filtered subject to the positional constraints imposed by the higherlevel features. Such a scheme is an alternative to existing surface smoothing/estimation techniques in microstructural materials science, in which the explicit treatment of different elements of the network topology is absent, or at best arbitrarily parameterized. We describe the ramifications of this technique in the usual microstructural applications in which the computation of important physical quantities is predicated on the precise estimation of the interface features. As an additional application, we describe a novel fronttracking algorithm for quantifying the transport of such interfaces.
1 Introduction
The morphology of surfaces and interfaces has garnered great interest in many fields of scientific and engineering research. Such studies have implications in applied physics, materials science, biology, pharmacology, chemical engineering and computer vision [1, 2, 3]. A vast part of this research is predicated on the proper imaging of interfaces in the medium of interest. Interfacial networks are composed of twodimensional, possibly curved interfaces that separate two distinct regions of homogeneous matter, such as gas in bubble foams, or phases or crystalline orientations in solids. We use in this paper language relevant to interfaces in polycrystalline materials but the methods described are equally applicable in other fields by straightforward adaptation of the terminology. The three dimensional entities with more or less uniform crystalline characteristics henceforth will be referred to as ‘grains’.
A particular type of grain boundary can be specified by five parameters on the mesoscale where ‘meso’ refers to a length scale that is large compared to interatomic distances but small compared to a typical grain size. Among the several possible parameterizations; we choose the set of three specifying the relative crystal orientations of the grains, and two specifying the local normal direction relative to the crystal axes in one of the grains. The normal direction in the other crystal frame can be computed from these five parameters. This parameterization ignores a microscopic relative translation on the atomic scale and thereby atomiclevel faceting of the interface, a feature addressed explicitly in molecular statics and dynamics simulations. The set of these five parameters is said to specify the grain boundary character[4, 5, 6, 7]. Note that the character between two grains can vary over the twodimensional boundary between them because, while the misorientation is fixed, the local normal typically varies significantly over a curved grain surface. Similar characterizations can be made for triple lines (two misorientations and a tangent line) and quad points (three misorientations). Finally, we note that crystal symmetry is typically exploited to reduce these specifications to unique ‘fundamental zones’ that span physically distinct ranges of orientations or misorientations.
Whether from a basic or applied science viewpoint, the importance of characterizing grain boundaries in this manner cannot be overstated. In polycrystalline materials, the local interfacial energy density and mobility are known to be sensitive to the five grain boundary parameters at each location [8, 9, 7, 10]. It also informs applications like grain boundary engineering whose eventual goal is to precisely manipulate bulk material properties through the tuning of the grain boundary character distributions [11, 12, 4, 13]. Further, it is wellknown that the topological elements of a grain boundary network like triple junctions and quad points are hotbeds of activity with respect to precipitate diffusion [14, 15, 16] and strain accumulation [17, 18]. Real grain boundary networks are usually the starting point for atomistic and continuum simulations of microstructure evolution, the physics of which is most difficult to model at triple lines and quad points.
All these applications are predicated upon measurements of the various topological features of a grain boundary network, which are inevitably subject to noise, whether through experimental resolution or image gridding. Given the generally accepted assumption of mesoscopically smooth interfaces, this necessitates the use of a smoothing estimator prior to any further analysis. Owing to the diverse roles of topological elements such as triple lines and quad points in microstructure phenomenology, an important motivation for this novel filtering technique and other recent ones [19] is to give them their due importance through explicit treatment.
Other factors motivating this work are:

Unlike voxelized images of most everyday objects, there exists no general intuition for the shape of a grain in a sample, and therefore a grain boundary. In the former case, iterative smoothing algorithms such as Laplace and Taubin smoothing [20] yield an acceptable result that is partially helped along by the user’s advance knowledge of the object in question. However these methods can suffer from under or oversmoothing if the number of iterations or step size are not chosen properly.

Explicit modeling techniques [21] more often than not belie the sheer variety in the observed structure of grain boundaries and network topologies.

Existing nonparametric techniques [22] require the use of a smoothing window of a userdefined size.
The methodology described here internally optimizes a compromise between fidelity to the input data points and a constrained Laplacian smoothing. An objective function is minimized with respect to this compromise. The algorithm requires no user input in terms of filter parameters, only that the connectivity of the nodes be specified in advance, in the form of a graph. We distinguish the type of kernel resulting from graphconnectedness to a given node from a fixedsize window centered on that node since the former, which we rely upon, does not take into consideration the physical distance between neighboring nodes, and only keeps track of the connectivity.
While the grains in polycrystals can take on essentially arbitrary shapes, the topological features typically encountered, and which are explicitly dealt with here, can be demonstrated with children’s (or adult’s) building blocks. Place two, say, cubic blocks (grains) together on a surface with edges aligned (blocks 1 and 2). The two blocks meet at a two dimensional interface (grain boundary 12). Now place a third block on the same surface so that it forms boundaries with both blocks 1 and 2 (boundaries 13 and 23). These boundaries meet at triple line 123. Now, place block 4 on top of these three so as to form boundaries with all three below (boundaries 14, 24, and 34). One now has new triple lines 124, 134, and 234. Furthermore, triple line 123 now terminates at a quadruple point, 1234, where all four grains meet. Unless one makes special alignments to again align edges, these are the topological features that will characterize an extended group of similarly stacked blocks.
In this paper, all line junctions of interfaces in a network are referred to by the generic term ‘triple line’ in allusion to the fact that energetically stable junctions in 3D are shared between exactly three interfaces. The incidental existence of a ‘tuple line’ in a polycrystalline material where interfaces intersect is known to be energetically unstable, forcing the interface topology to deform to a lower energy configuration [23, 24]. Likewise, a node of termination of such triple lines is referred to as a ‘quad point’ irrespective of the value of , alluding to the fact that is the physically stable configuration in bulk materials. The mathematical machinery developed in this paper is as appropriate for contrived interface networks that deviate from this this topological rule as for digital images of real bulk microstructure, in which these deviations are almost never observed.
We first describe the topological hierarchy in general terms and then address the interface estimation procedure, which is a modification of Laplacian smoothing of a set of meshed surface points. This is followed by the application of the estimation algorithm to pixelated versions of easily parameterized geometric primitives, in particular circles, spheres and cylinders. Postsmoothing errors are quantified in terms of estimated sizes of these primitives as well as estimated normals for specific geometries. We then address specific cases of interest in mesoscale materials science: two and threedimensional grain boundary networks. The former finds relevance in the study of thin films and the latter in that of bulk material behavior (most prominently in the computation of grain boundary character distribution plots, a common characterization of materials microstructure). We demonstrate how the user is freed from the largely intuitive choices of smoothing parameters that is characteristic of iterative or windowed techniques. We then describe in some detail the applicability of this surface estimation algorithm to finite element methods in materials science as well as interface velocity estimation, which is a new capability made possible with data obtained from modern nondestructive imaging techniques. A new nonparameteric algorithm to achieve the latter is described.
2 General formalism
Consider a set of noisy sample points in dimensional Cartesian coordinates that denote an imaged grain boundary. A subset of these points is tagged as a ‘perimeter’ that samples the edges of the grain boundary feature, with the same grid resolution as the interior. For example, in three dimensions, these points could represent a twodimensional boundary including all edge points, or a onedimensional triple line with its terminating quad points. We also specify a connectivity for every point in , described by a graph Laplacian matrix :
(1) 
where is the number of points connected to and is an indicator function that is if point is connected to point and otherwise. We require that all remain constrained to their initial positions while the are smoothed, all the while adhering to the same node connectivity. We denote this smoothing operation notionally by .
As a general rule, we enter points into our hierarchy such that all are at one level above all . Notationally the hierarchy level or ‘rank’ is denoted by a function such that ; the sole purpose of being to distinguish points of different ranks and the actual returned value being a matter of choice.
Figure 1 visualizes two common systems with different hierarchy sizes. Keeping in mind that in an interface network in dimensional space there exist in general objects of dimensionality , we define the rank function , where corresponds to the lowestdimensional object in the network to which belongs. For example, a triple point in a 2dimensional image is assigned a rank of because it is a zerodimensional object, while a grain boundary interior point has a rank of .
Type of  

Triple point  0  2 
boundary interior  1  1 
Type of  

Quad point  0  3 
Triple line  1  2 
Boundary interior  2  1 
Tables 1 and 2 show the general rule that for a given network, . We note that the feature of a topological element that decides its rank is its dimensionality rather than its name. For example, if a quad line existed in a network for which (i.e. intersection of four grain boundary surfaces) its rank would be . Based upon these definitions, the smoothing algorithm for a set of interface points in a dimensional network is as follows:
In summary, points of rank are smoothed while holding in place all previously smoothed connected points of rank , with highestrank points essentially undergoing unconstrained smoothing (since is initially an empty set). If the highest rank elements have as do quad points when or triple points when , then one can skip ‘smoothing’ them altogether. This scheme gives the aforementioned topological features their due importance relative to one another. The prerequisite of having sample points labeled according to topological feature is readily achievable by nearest neighborbased clustering algorithms [25], which is in fact the very information represented in the NodeType dataset in a typical DREAM.3D microstructure file.
3 Constrained smoothing
is based on a nonparametric regression that involves penalizing, in Cartesian componentbycomponent fashion, the displacement between each estimated smoothed point and its unsmoothed neighbors. If of initial points are mobile (), a measure of the nearest neighbor fluctuations of each Cartesian component of can be estimated with , where represents a vector of only the that require smoothing, is a modified graph Laplacian operator expressing the connectivity of the mobile nodes and denotes constants that are determined from the remaining . in fact specifies the Dirichlet boundary conditions to Laplace’s equation. Specific examples of and are described presently. In the case of no constraints, , is an empty set and is the full graph Laplacian introduced in Equation (1). performs simultaneous filtering of each component by negotiating a tradeoff between fidelity to the raw data and minimization of fluctuations between smoothed neighbors through a scalar control parameter . A control function is defined to this end:
(2)  
Here represents the array corresponding to that is further along in the smoothing process. At the extreme values of and , the minimizer of respectively favors complete data fidelity () and complete Laplacesmoothing (). We further define an objective function that penalizes fluctuations between each smoothed point and its nearest unsmoothed neighbors based on the connectivity specified in the full graph Laplacian :
(3) 
Crucially, we require that the minimizer of be reached by always satisfying the optimality condition of the control function in (2) with respect to and therefore indirectly through variation of the parameter alone. This makes the smoothing operation on the a onedimensional minimization in that can easily be achieved by a binary search in the interval . Briefly, the objective function is the actual quantity being minimized in the regression, but the path taken in the objective function landscape is decided by the control function .
We define the matrix of unsmoothed starting points as having rows and columns where is the dimensionality of the points. Similarly we define the identicallysized matrix as the solution resulting from applying to . We define and as the diagonal and adjacency matrices respectively of . We rely on the following intermediate definitions to obtain the reduced Laplacian and constant matrices:

Let the integer set denote the indices of the points that remain fixed (i.e. or equivalently )

If for an integer , is an integer set such that , then let , i.e. the complement of with respect to .

Let the submatrix of a matrix formed by:

the rows whose indices are in be denoted by .

the rows and columns whose indices are in be denoted by .

then the reduced Laplacian and constant matrix are defined:
(4)  
(5)  
For example, if points in dimensions are to be smoothed in which the are connected sequentially with and to be fixed, then:
(6) 
Equation (6) denotes the full smoothed solution including the constrained points, with as defined earlier. If and respectively satisfy the freeboundary and constrained Laplace equations, then is it clear that and . The smoothing problem is stated more compactly as the following optimization problem:
(7)  
(8)  
We note from Equation (8), which is derived from the minimizer of the control function in Equation (2), that is an matrix and that in Equation (7) the argument of the trace operator is a symmetric matrix with nonnegative eigenvalues (the case of zero eigenvalues implies that the sample points are flattened in at least one dimension). Significantly, the trace of this matrix and therefore the objective function itself represents the aggregate squared Euclidean distance of each node from its unsmoothed neighbors. The equivalence of Equations (3) and (7) is seen in the simple onedimensional smoothing example (): . Each column in in the matrix corresponds to one such Cartesian component in the sample frame of reference and each element is of the form . This is precisely the argument of the operation in (3). The objective function in (7) represents the operation in (3) being performed simultaneously on all Cartesian components. Minimizing them simultaneously is completely equivalent to minimizing the (reference frameinvariant) trace of the matrix .
The algorithm is finally given by:
The resulting surface consisting of the smoothed points with the preserved original connectivity is nonparametric. Qualitatively, the algorithm attempts to determine the least jagged surface passing in between the sample points, thus maintaining data fidelity. This precludes a major problem in applying iterative Laplacelike techniques, that of over or undersmoothing. For the applications of in the remainder of this text, the threshold value of was taken to be .
We point out that the smoothing scheme outlined in this section allows users to define the components of the algorithm for specific requirements. For instance, users of mesh smoothing algorithms like finite element method (FEM) might want to explicitly incorporate mesh quality metrics into the objective function as an alternative to remeshing. This in turn may well decide the manner of stepping in the objective function landscape and therefore shape the control function. Our smoothing paradigm permits the flexibility of userdefined objective and control functions, all the while heeding the hierarchy between the components of the grain boundary network. We address the applicability of smoothed meshes obtained from the objective function in Equation (7) to finite element applications in Section 8.
4 ing a digitized curve in a plane
We describe as a first demonstration of the problem of filtering a set of pixelated points representative of a curve in a plane. We list the following general properties of such a set of points:

Every point has at least one nearest neighbor in at least one direction on this integer grid. This is characteristic of discretized sampling of continuous curves and surfaces in general.
Such points can be obtained from pixelated images in standard formats by first generating a phase field (for instance a field of unique integers assigned to each grain), taking the magnitude of the gradient of this field and binarizing it. A morphological ‘skeletonizing’ operation can then be applied to this binarized field [29, 30, 31]. This same technology is used in the field of biometrics, for example, to thin down highresolution images of fingerprints to features of singlepixel thickness for further analysis. For nfHEDM images, one can collect directly the voxel (volume pixel) edges that border two different grains, as decided by some segmentation criterion. In our example the coordinates of the sample points are integers on a square grid. The image of the curve has been skeletonized to ensure that each sample point has no more than two of the immediate squaregrid neighbors belonging to the pixelated line (interior points have two neighbors and the terminal points have one). The results of constrained smoothing on such a pixelated curve is shown in Figure 2 with the perimeter points fixed at the unsmoothed grid point locations.
In Figure 2(a) the unfiltered sample points are taken to be the pixel centers. That the final estimated curve passes almost through most of the pixel corners is a consequence of the algorithm attempting to choose a smooth path between the unfiltered sample points, much like an oldfashioned wooden spline negotiating a smooth path between its fixed control points. Given a fixed boundary condition, the skeletonization of the original pixelated curve ensures that the smooth estimate stays within a pixel width of the original smooth curve. This is also seen in the smooth estimates of the twodimensional crosssections of grains in Figure 6(d).
5 ing known shapes
In this section the hierarchical algorithm is applied to open surfaces that are geometry primitives in two and three dimensions. Simple parameterizations for these primitives provide a means of comparison with a smoothed solution on a pointtopoint basis. We describe trends in the errors for different primitives as a function of voxel density. Focusing attention on open surfaces allows us to simulate smoothing in the presence of topological features characteristic of a grain boundary network. The three primitives chosen are a 2D circle, a 3D sphere, and a 3D cylinder. Surface slices from these primitives were characterized by the following:

The gridding resolution was chosen in terms of the number of voxels per unit length, .

Circle: A semicircular arc of unit radius, whose endpoints were reset to unit radius after discretization. These endpoints were held constrained during smoothing.

Sphere: A square patch spanning in two mutually perpendicular directions cut out from the surface of a sphere of radius units, with the edges of the square treated as triple lines and the vertices as quad points. The quad points alone were constrained to lie on the sphere, while the others were subject to discretization on a cubic lattice.

Cylinder: A rectangular patch cut from the surface of a cylinder of radius units, parallel to its axis and spanning along the azimuth. The edges of the rectangle were treated as triple lines and the vertices as quad points. The quad points were constrained to remain on the surface of the cylinder.
The choice of the sphere and cylinder radii are indicative of the typical size of a grain from earlier nfHEDM measurements [32]. The quality of smoothing was expressed as the error in the estimated radius of the primitive in question. Specific examples smoothing on these primitives are shown in Figure 3.
The fidelity of the final smoothed result to the original primitive was quantified in terms of the pointtopoint difference in radii of the original and smoothed surfaces. In the case of the sphere and the cylinder, the first spherical and cylindrical polar coordinates are respectively used (i.e and ). Shown in Figure 4 are trends in the estimated error and its standard deviation , taken over the smoothed mesh nodes, for all three primitives. If the unit of length is taken to be a millimeter, the relative error in the region between and is particularly relevant for techniques like nfHEDM since they correspond to a pixel size range of to , which brackets the known experimental resolution [33]. For comparison, the radii of the spherical and cylindrical patches were chosen to be .
Primitive  

Circle  
Sphere  
Cylinder 
The power law behavior in the error trends was tested for each primitive as is seen in the straight line fits in Figure 4. The coefficients and of the estimated power law are listed in Table 3. Of particular interest is which is seen to lie close to for all three primitives. This is simply explained by the fact that the length error is , the size of one voxel. The lowered fit quality for the sphere and cylinder is attributed to the difficulty in obtaining a perfect stairstepped mesh for these primitives. We further note that the relative error for each primitive is around a fraction of a percent () at the spatial resolution of nfHEDM ( or pixels per millimeter).
Another smoothing quality metric that is easily calculated for two dimensions is the error in the local normal of a curve. This is particularly relevant to surface imaging applications. We estimate the local deviation about the known normal of the sections of a pixelated straight line that has been ed with its endpoints held at their true positions. This deviation is determined for different inclinations of the original line to the Cartesian grid. A schematic and results for the inclination range of to is shown in Figure 5. It is seen that while the mean deviation varies with the inclination of the original line, it falls within a few degrees of the actual normal which justifies the use of in calculations of grain boundary character distribution [34, 5, 6].
6 Results: twodimensional microstructure smoothing
The constrained smoothing is next demonstrated on a real twodimensional microstructure imaged with nfHEDM. The sample points of each grain boundary are expressed in integer units of suitable inplane step size (in this case, ) and are classified as belonging to a grain interior, grain boundary or triple point depending on the number of unique grains represented in the neighborhood . The optimization is performed simultaneously over and and therefore is an matrix. The results on a section of microstructure are shown in Figure 6.
7 Results: threedimensional microstructure smoothing
Hierarchical smoothing is demonstrated on select grains of a wellordered threedimensional microstructure measured with nfHEDM and is compared to the results of Laplace smoothing. The prerequisite node connectivity on the grain surfaces was obtained by first segmenting the microstructure into its constituent grains and then triangulating the faces of the cubic voxels along grain boundaries [25]. The result of this operation is a ‘quickanddirty’ Delaunay mesh on stepped grain surfaces characteristic of discrete sampling. A few important points about this bookkeeping process that inform the subsequent smoothing are:

The grain surface nodes are unambiguously classified into their topological types i.e. boundary interiors, triple lines or quad points. It is worth mentioning that the preprocessing described above is external to the hierarchical smoothing algorithm itself and as such is not the focus of this work.

The meshing and bookkeeping is done in such a manner as to ensure that nodes along triple lines and at quad points are shared between the neighboring topological features.
Figures 7, 8 and 9 show comparisons of hierarchical smoothing with Laplace smoothing in which the ease of movement of the triple points is enhanced by assigning them a greater Laplace smoothing parameter . The values of for the parent volumes for each of these grains were determined by trial and error (as a user would have to do) by visually minimizing the distortion from the original squaregridded grain. This is a highly inefficient process that is not required in hierarchical smoothing.
8 Mesh quality
We address the suitability of output for finite element applications, which are predicated on the availability of surface meshes with reasonably isotropic mesh elements (equilateral, in the case of triangular) in order to avoid errors from piecewise linear interpolation. While there exist other sophisticated methods of quantifying the quality of a mesh element [35], we implement a simple metric for triangular elements that tests their closeness to an equilateral triangle [36]. The quality of a triangle of area and side lengths , and is defined to be:
(9) 
which gives for an equilateral triangle. Shown in Figure 10 are elementwise quality plots of select grains in a 3dimensional volume. Figure 11 histograms the mesh quality over all surface elements in the entire 424 grain volume. The flattening of some of the mesh elements at the triple junctions of the 3dimensional grains in Figures 7, 8 and 9 is evidenced by the slight peak in the distribution in Figure 11 at low qualities ().
9 Front tracking
We describe as a final application of surface ing a methodology to infer directly the velocity field on a sampled interface on a nodetonode basis given the two snapshots of the interface before and after its differential migration, what we call the fronttracking problem. We present this novel computation as one riding the current wave of synchrotronbased nondestructive measurement techniques such as HEDM and DCT, and whose eventual applicability is predicated upon the estimation of the geometry of mesoscopically smooth grain boundaries. It is generalized insofar as to account for changes in the interface geometry, such as size (or equivalently number of sample points), local curvature, surface reorientation and torsion. An example of of a pysical measurement that would serve as input to this algorithm would be a grain coarsening experiment by ex situ annealing [32] in which bulk grain boundary structure is measured before and after a shortterm annealing regimen. A satisfactory quantification of the intermediate stages of boundary migration in this manner opens up the possibility of studying the influence of local curvaure and its rate of change on coarsening. The emphasis on differential migration is significant because it permits us to make a qualitative assumption about the interface transport: that the sample nodes of the interface take the straightest and simplest path in the intervening space to reach the configuration described by the target set of points. A schematic of such transport is shown in Figure 12 for a variety of morphological distortions.
We state our problem in terms of two Cartesian point sets and in 3D space ( and in number respectively). These two sets represent an interface before and after differential migration. Note that may not be equal to . In practice this could be because the size of the grain boundary changed during transport or the sampling of the interfaces was done at different resolutions, or a combination of both.

We seek a set of displacement vectors , one for each such that the set mimics the configuration of the target surface as closely as possible, with as uniform a sampling as posible.

Under these conditions, it is assumed that in the time of transport (usually obtained from experiment) of the intermediate positions of the boundary nodes are:
(10) 
The velocity field over the nodes is then simply given by:
We seek to compute the node displacements since the velocities are estimated to first order with knowledge of the time of transport .
We demonstrate how this problem is solvable using linear optimization for which there exist a variety of numerical software packages. We calculate each displacement as a weighted sum of all possible displacements from to each :
(11) 
The unknown weights are determined by the following linear optimization problem in the unknowns :
(12) 
subject to the constraints:
(13)  
(14)  
(15) 
The objective function in Equation (12) attempts to find the smallest possible aggregate weighted squared displacement between the two point sets. The positivity constsaint (13) ensures directionality of the displacements towards the target surface. Constraint (14) ensures that the actually terminate on the target surface and not before or after, while the equiweight constraint (15) adjusts the displacements so that the target surface is as uniformly sampled by the point set as can be managed. The inferred displacements of the sample nodes of a meshed surface in 3D is shown in Figure 13. Additional examples can be seen in the supplementary material.
It is wellknown, however, that the true quantity of interest is the local normal velocity of the interface. This is easily estimated for each surface mesh element from the inferred velocities of its constituent nodes. For example, if the nodes of a triangular surface mesh element of inclination have inferred displacements then the normal displacement of the element is estimated as the projection of the displacement of the element centroid along the direction of inclination:
(16) 
In closing, we note the following points:

The linear optimization problem always has a unique solution and therefore the front tracking algorithm never fails. It falls on the user to decide whether the interface data satisfy the requirement of differential migration. For interface migration over extended periods of time in which the final interface is dramatically different from the initial interface in geometry or topology, it may very well be that this configuration was not achieved through migration along straight lines in actual physical samples.

It is observed in Figure 13 that the edge nodes of the initial interface are not mapped exactly to the edges of the final interface; in fact they are mapped to a location decidedly in the interior of the target interface. This is characteristic of node transport between point sets of different sizes: . In order to retain topological consistency one may adopt a convention of hierarchical transport similar to the smooth surface estimator described earlier, in which the linear optimization algorithm is applied separately to subsets of the point sets. For example, quad points in the initial surface are explicitly transported to quad points in the target surface, triple line nodes to triple line nodes and grain boundary interior nodes to grain boundary interior nodes. The local normal displacement can be estimated as before.
10 Summary and discussion
A new smoothing estimator for interface networks was demonstrated and compared to the performance of an established but generic smoothing algorithm in current use. The requisite inputs to the algorithm are:

Discretely sampled polycrystal interface data, particularly measured by modern highresolution synchrotronbased experiments such as nfHEDM.

Topological characterization of the sample points as belonging to grain boundary interiors, triple lines or quad junctions. This is easily achieved by existing software packages like DREAM.3D [25].

A mesh of the discretized sample points, also achieved quite simply by the QuickMesh feature of DREAM.3D.
The new algorithm organizes the topological elements of the network into an hierarchy depending on which elements physically border other elements and smooths each element set constraining its bordering elements to their fixed positions. This treatment gives the physically relevant higherorder topolological elements like triple lines and quad points their due consideration and retains the geometric discontinuities along a grain surface resulting from their existence. By its very design the smooth estimator returns a curve that passes in between the original noisy sample points, which is highly significant for pixelated images obtained from digital measurement techniques. Reorienting the sampling grid still ensures that the estimated surface lies within a pixel width of the true interface. Further, all elements belonging to a particular hierarchy rank in the entire volume are smoothed simultaneously so that they are ready to be used as Dirichlet boundary conditions for elements of lower rank that connect to them. The method is completely nonparametric, permitting the automated smoothing of imaged bulk structures and does not suffer from userrelated effects like over, undersmoothing or fixedsize window artifacts. Repeated applications of the smoothing on the same point set results in better smooth approximations with decreasing the extent of waviness along the smoothed surface.
The technique is predicated on the nearest neighbor connectivity of the surface nodes being known in advance, which is easily achievable by existing algorithms and is already implemented in opensource microstructure software packages [25]. The additional requirement of labeling surface nodes as belonging to either grain boundary interiors, triple lines or quad points is also achieved by these software packages. The case of faceted grain boundary interiors is likewise handled by appropriate labeling of this kind, for purely bookkeeping purposes. This is a task external to the smoothing itself and as such is outside the scope of this work.
The relative errors in the smoothing of known shapes are demonstrated to be a fraction of a percent for typical resolutions of microstructure imaging techniques. The estimated normals in the case of two dimensions were found to be within thresholds characteristic of bin sizes used in plots of grain boundary character distribution (GBCD). The ability to handle data of different dimensions in a generalized manner allows this technique to be used for surface experiments such as optical metallography and EBSD, as well as 3D bulk techniques like nfHEDM. Computation of the mesh quality on a 3dimensional grain boundary network consisting of 424 grains revealed that the overwhelming majority of mesh elements have a quality above a comfortable as required by simple applications. If need be one may remesh the output of the smoothing routine subject to junction constraints in order to obtain a more uniform meshing of the grain boundary interiors.
We further detailed a new computational scheme to quantify the local displacement of a grain boundary using mathematical optimization techniques. The estimated transport of interfaces, coupled with the nonparamatric estimation of interface geometry using , potentially sets the stage for direct study of the influence of local curvature and its rate of change on interface migration in real materials, something which has not yet been achieved in three dimensions for real polycrystals.
In summary, hierarchical smoothing as a surface estimation technique, apart from yielding aesthetically pleasing results, finds applications in the computation of grain boundary character distribution (GBCD), finiteelement modeling of microstructure and the study of influence of curvature on interface migration in polycrystals, the latter uniquely in conjunction with modern synchrotronbased nondestructive measurement techniques. As a further implementationbased consideration, the separate treatment of the different topological features suggests parallelizeability in order to reduce computation time.
11 Acknowledgements
The authors would like to thank David Menasche and Anthony Rollett at Carnegie Mellon University for their valuable input. This research was supported by NSF Award DMR1105173.
12 References
References
 [1] MJ McCready, Eleni Vassiliadou, and TJ Hanratty. Computer simulation of turbulent mass transfer at a mobile interface. AIChE Journal, 32(7):1108–1115, 1986.
 [2] Tim McInerney and Demetri Terzopoulos. A dynamic finite element surface model for segmentation and tracking in multidimensional medical images with application to cardiac 4d image analysis. Computerized Medical Imaging and Graphics, 19(1):69–83, 1995.
 [3] H. et al. Lyklema. Fundamentals of Interface and Colloid Science. Academic Press, 2005.
 [4] David M Saylor, Bassem S El Dasher, Anthony D Rollett, and Gregory S Rohrer. Distribution of grain boundaries in aluminum as a function of five macroscopic parameters. Acta Materialia, 52(12):3649 – 3655, 2004.
 [5] A. Khorashadizadeh, D. Raabe, S. Zaefferer, G. S. Rohrer, A. D. Rollett, and M. Winning. Fiveparameter grain boundary analysis by 3d ebsd of an ultra fine grained cuzr alloy processed by equal channel angular pressing. Advanced Engineering Materials, 13(4):237–244, 2011.
 [6] Sutatch Ratanaphan, Yohan Yoon, and GregoryS. Rohrer. The five parameter grain boundary character distribution of polycrystalline silicon. Journal of Materials Science, 49(14):4938–4945, 2014.
 [7] Sutatch Ratanaphan, David L. Olmsted, Vasily V. Bulatov, Elizabeth A. Holm, Anthony D. Rollett, and Gregory S. Rohrer. Grain boundary energies in bodycentered cubic metals. Acta Materialia, 88:346 – 354, 2015.
 [8] David L. Olmsted, Stephen M. Foiles, and Elizabeth A. Holm. Survey of computed grain boundary properties in facecentered cubic metals: I. grain boundary energy. Acta Materialia, 57(13):3694 – 3703, 2009.
 [9] David L. Olmsted, Elizabeth A. Holm, and Stephen M. Foiles. Survey of computed grain boundary properties in facecentered cubic metalsâii: Grain boundary mobility. Acta Materialia, 57(13):3704 – 3713, 2009.
 [10] Vasily V. Bulatov, Bryan W. Reed, and Mukul Kumar. Grain boundary energy function for fcc metals. Acta Materialia, 65:161 – 175, 2014.
 [11] David M. Saylor, Adam Morawiec, and Gregory S. Rohrer. The relative free energies of grain boundaries in magnesia as a function of five macroscopic parameters. Acta Materialia, 51(13):3675 – 3686, 2003.
 [12] David M. Saylor, Adam Morawiec, and Gregory S. Rohrer. Distribution of grain boundaries in magnesia as a function of five macroscopic parameters. Acta Materialia, 51(13):3663 – 3674, 2003.
 [13] Jia Li, Shen J. Dillon, and Gregory S. Rohrer. Relative grain boundary area and energy distributions in nickel. Acta Materialia, 57(14):4304 – 4311, 2009.
 [14] Alexander H. King. Triple lines in materials science and engineering. Scripta Materialia, 62(12):889 – 893, 2010. Viewpoint set no. 46 âTriple Linesâ.
 [15] B. Bokstein, V. Ivanov, O. Oreshina, A. Peteline, and S. Peteline. Direct experimental observation of accelerated zn diffusion along triple junctions in al. Materials Science and Engineering: A, 302(1):151 – 153, 2001.
 [16] T.P. Swiler, V. Tikare, and E.A. Holm. Heterogeneous diffusion effects in polycrystalline microstructures. Materials Science and Engineering: A, 238(1):85 – 93, 1997. Microstructure Evolution in Bulk Phases F.
 [17] A D Rollett, R A Lebensohn, M Groeber, Y Choi, J Li, and G S Rohrer. Stress hot spots in viscoplastic deformation of polycrystals. Modelling and Simulation in Materials Science and Engineering, 18(7):074005, 2010.
 [18] J.L.W. Carter, N. Zhou, J.M. Sosa, P.A. Shade, A.L. Pilchak, M.W. Kuper, Y. Wang, H.L. Fraser, M.D. Uchic, and M.J. Mills. Characterization of Strain Accumulation at Grain Boundaries of NickelBased Superalloys, pages 43–52. John Wiley & Sons, Inc., 2012.
 [19] SB Lee, G S Rohrer, and A D Rollett. Threedimensional digital approximations of grain boundary networks in polycrystals. Modelling and Simulation in Materials Science and Engineering, 22(2):025017, 2014.
 [20] G. Taubin. Curve and surface smoothing without shrinkage. In Proceedings of the Fifth International Conference on Computer Vision, ICCV ’95, pages 852–, Washington, DC, USA, 1995. IEEE Computer Society.
 [21] Jun Wang and Zeyun Yu. A novel method for surface mesh smoothing: Applications in biomedical modeling. In BrettW. Clark, editor, Proceedings of the 18th International Meshing Roundtable, pages 195–210. Springer Berlin Heidelberg, 2009.
 [22] E J Lieberman, A D Rollett, R A Lebensohn, and E M Kober. Calculation of grain boundary normals directly from 3d microstructure images. Modelling and Simulation in Materials Science and Engineering, 23(3):035005, 2015.
 [23] L.A. BarralesMora. 2D and 3D Grain Growth Modeling and Simulation. Cuvillier, 2008.
 [24] G. Gottstein and L.S. Shvindlerman. Grain Boundary Migration in Metals: Thermodynamics, Kinetics, Applications, Second Edition. Materials Science & Technology. CRC Press, 2011.
 [25] Michael Groeber and Michael Jackson. Dream.3d: A digital representation environment for the analysis of microstructure in 3d. Integrating Materials and Manufacturing Innovation, 3(1):5, 2014.
 [26] Adam J. Schwartz, Mukul Kumar, Brent L. Adams, and David P. Field, editors. Electron Backscatter Diffraction in Materials Science. Springer Science + Business Media, 2009.
 [27] U. Lienert, S.F. Li, C.M. Hefferan, J. Lind, R.M. Suter, J.V. Bernier, N.R. Barton, M.C. Brandes, M.J. Mills, M.P. Miller, B. Jakobsen, and W. Pantleon. Highenergy diffraction microscopy at the advanced photon source. JOM, 63(7):70–77, 2011.
 [28] S. F. Li and R. M. Suter. Adaptive reconstruction method for threedimensional orientation imaging. Journal of Applied Crystallography, 46(2):512–524, Apr 2013.
 [29] Robert M. Haralick and Linda G. Shapiro. Computer and Robot Vision. AddisonWesley Longman Publishing Co., Inc., Boston, MA, USA, 1st edition, 1992.
 [30] Louisa Lam, SeongWhan Lee, and Ching Y. Suen. Thinning methodologiesa comprehensive survey. IEEE Trans. Pattern Anal. Mach. Intell., 14(9):869–885, September 1992.
 [31] T. Yung Kong and Azriel Rosenfeld, editors. Topological Algorithms for Digital Image Processing. Elsevier Science Inc., New York, NY, USA, 1996.
 [32] CM Hefferan, SF Li, J Lind, U Lienert, AD Rollett, P Wynblatt, and RM Suter. Statistics of high purity nickel microstructure from high energy xray diffraction microscopy. Computers, Materials, & Continua, 14(3):209–220, 2010.
 [33] C. M. Hefferan, S. F. Li, J. Lind, and R. M. Suter. Tests of microstructure reconstruction by forward modeling of high energy xray diffraction microscopy data. Powder Diffraction, 25:132–137, 6 2010.
 [34] S D Sintay and A D Rollett. Testing the accuracy of microstructure reconstruction in three dimensions using phantoms. Modelling and Simulation in Materials Science and Engineering, 20(7):075005, 2012.
 [35] Pascal J. Frey and Houman Borouchaki. Surface mesh quality evaluation. International Journal for Numerical Methods in Engineering, 45(1):101–118, 1999.
 [36] Randolph E. Bank. PLTMG: a software package for solving elliptic partial differential equations, volume 7 of Frontiers in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1990. Users’ guide 6.0.
 [37] MATLAB. version 8.5 (R2015a). The MathWorks Inc., Natick, Massachusetts, 2015.