POLE FIGURE INVERSION USING FINITE ELEMENTS OVER RODRIGUES SPACE

Using finite elements over Rodrigues space, methods are developed for the formation and inversion of pole figures. The methods take advantage of the properties of Rodrigues space, particularly the fact that geodesics corresponding to pole figure projection paths are straight lines. Both discrete and continuous pole figure data may be inverted to obtain orientation distribution functions (ODFs) in Rodrigues space, and we include sample applications for both types of data.


INTRODUCTION
Since the early work of Roe (1965) and Bunge (1969), the analysis of texture, or preferred crystallographic orientation, has remained an active research area. The texture is quantified by an orientation distribution function (ODF) which often is regarded as the fundamental entity of quantitative texture analysis and whose determination is of prime importance. The measurement of pole density data using dif-*Corresponding author. Present Address: Department of M&AE, EBU2, UC San Diego, 9500 Gilman Dr., La Jolla, CA 92093-0411, USA, e-mail: nrbarton@ucsd.edu Software related to the methods presented in this paper will be made available at http:// anisotropy.mae.cornell.edu fraction methods provides a non-destructive route to determining the ODF of a material sample. While automated electron back-scatter diffraction (EBSD) methods now offer a viable alternative to determine an ODF (see Kunze et al., 1993;Weiland, 1994), the inversion of pole figure data from X-ray or neutron diffraction goniometers remains an important method of ODF determination. Several recent works (Bunge, 1982(Bunge, , 1985(Bunge, , 1999Wenk, 1985;Kocks et al., 1998) provide excellent summaries of available methodologies for this purpose. The basic task is the same for all approaches: evaluate a set of unknown coefficients that quantify an ODF such that the ODF is consistent with the pole figure data. For a specified sample direction, pole density contributions are obtained from all orientations that bring a given crystallographic direction or its negative into coincidence with the sample direction. These sets of contributing orientations are related by rotations about established axes and therefore, coincide with geodesics in the space of proper rotations. Thus the density of a pole figure for a particular sample direction is affected by the value of the ODF at all points along a pair of integration paths.
In all of the approaches, two characteristics of the representation must be specified: the parameterization of orientation space and the functions used to quantify the ODF over this parameterization. Even though the basic task of pole figure inversion is always the same, the specific methodologies to relate the coefficients of the ODF representation to the pole figure densities differ substantially. Important distinctions among the various approaches emerge in the context of dealing with the indeterminacies at the heart of this task and in relation to issues of sparse or contradictory pole figure data.
Indeterminacy in the evaluation of the ODF from pole figures arises from two sources. First, infinitely many pole figures are required to exactly resolve the ODF. However, the ODF can be evaluated approximately from limited pole figure data. Specifically, given that the ODF is represented with a function that is completely expressed using a finite number of parameters, such as a limited number of terms in a series expansion, those parameters can be determined from pole data. Ideally, the greater the amount of pole figure data, the better the approximation. For diffraction of polycrystalline materials, the number of useful poles is quite limited, so precision of the ODF in practice is restricted. Attempting to resolve the ODF to greater detail than the data allow introduces an ill-posedness, even when a finite number of parameters is employed.
A second source of indeterminacy enters in those cases where there exists antipodal symmetry in the pole figure data. This indeterminacy is not eliminated with increased numbers of pole figures as the experimental methods used to generate the pole figure data cannot distinguish between positive and negative along the same direction. As a result of this antipodal pole figure symmetry, the odd part of any functional description of the ODF does not influence the pole figure values. In the past, this motivated the complete exclusion of the odd terms in harmonic series expansions. Although this is now known to be too restrictive, the fact remains that the data allow only for determination of the even part of any function. Consequently, there remains an ambiguity in the odd part of the representation that, at least in principle, allows for multiple versions of the ODF that all render identical pole figures. Two factors intervene. First, additional constraints placed on the range of the ODF appear to mitigate this problem. In particular, acceptable ODFs have been obtained by forcing the ODF to be strictly non-negative rather than being comprised solely of even functions. Second, while the ODF itself is of general interest, the desired end usage is more often the ability to compute effective properties by taking orientational averages. For common crystal types, these averages are unaffected by the 'ambiguous' (or odd) part of the ODF. If the ODF reproduces the pole figure data correctly, then the orientational average taken with it correctly determines the average property. The method we present here does not attempt to resolve this issue of indeterminacy differently from other approaches previously reported (and summarized in the cited books). The representation of the ODF that we present has both odd and even parts in the functional form and may be constrained to be non-negative.
Other factors can further complicate the task of determining an ODF. Physical constraints in diffraction experiments can mean that data around the periphery of a pole figure are missing or are of poor quality. For some approaches, namely some implementations of the harmonic method, completion of the pole figure is required prior to performing the pole figure inversion to obtain the ODF. In the methodology presented here, pole figure completion is not required as the inversion may be performed using discrete values of the pole density together with an optimization method. However, it is possible to perform the inversion using pole figure data that have already been reduced to the form of smooth functions over the hemisphere. Examples of both methods are given.
Measurement errors may result in pole densities for different diffraction conditions (poles) that are not derivable from a common ODF and thus are inconsistent in this respect. The references cited above discuss these issues for a number of approaches, and we shall not reiterate that information here. However, in introducing a new methodology for determining the ODF from diffraction data, we recognize that we should address these issues at least conceptually to verify that the approach is capable of delivering an ODF in the face of inconsistent data. In rendering an ODF, the methods developed below arbitrate between conflicting data by minimizing a chosen objective function.
We begin by describing the new method in relation to existing ones. The most widely used procedures are broadly classified into two categories: series methods and discrete methods. Generally, with series methods the ODF is represented over an Euler angle domain with spherical harmonic functions; with discrete methods the ODF is represented using bins which partition an Euler angle domain. We discuss the new procedure from a perspective that is expanded somewhat from this classification using two defining characteristics for the representation of the ODF: the parameterization of orientation space and the interpolation of the ODF over the parameterization. There are two common choices for the orientation parameterization: Euler angles (of which there are several variants that get regular use) and angle-axis vectors (also having several variants, of which the Rodrigues parameters are the most common). Angle-axis representations may be considered a special case of quaternions (see Altmann, 1986), which recently have been discussed in the context of texture analysis by Schaeben et al. (2001b).
With regard to the representation of the ODF over a particular parameterization, we also consider two general categories: (1) bases comprised of a linear combination of functions each of which spans the entirety of orientation space (referred to here as global interpolants), and (2) bases comprised of a linear combination of functions that individually span only a subdomain of orientation space (referred to here as piecewise or local interpolants). In the latter case, the bases spanning the full orientation space are a patchwork of the functions over the subdomains (and are the traditional choice in finite element formulations). For perspective, we note that harmonic methods typically utilize global interpolants (harmonics) over an Euler angle parameterization of the orientation domain. The discrete methods employ piecewise interpolation (piecewise constant functions) typically over an Euler angle parameterization. Thus, from this perspective, both employ the same type of orientation parameterization (Euler angles), but differ in the choice of interpolation functions. The global interpolants in harmonic methods are continuous and have continuous derivatives. The piecewise interpolation in the discrete methods are constant within the bin and are discontinuous at the bin interfaces. Many other combinations of parameterization and basis functions are possible. We emphasize that the ODF can be adequately represented with a variety of combinations, and we are not claiming that any one of these is intrinsically superior. Our intent is to present a methodology for inversion of pole figure data using a different combination than has been reported previously and to demonstrate some of the advantages possessed by this combination.
The approach presented here is distinct from prior methodologies in several ways. It employs local interpolation functions over an angle-axis parameterization of orientation space. The interpolation functions are piecewise linear functions associated with tetrahedral finite elements. They display C continuity; the value of the ODF is continuous everywhere, but possesses discontinuous derivatives at the element interfaces. Within the elements, the interpolation has smooth derivatives, although with linear interpolation functions the second and higher derivatives are identically zero. The choice of orientation parameterization is based on the advantageous properties of the Rodrigues vectors. As discussed in the following sections, an extremely convenient property relates to the integration paths for pole projections: in Rodrigues space, these paths all fall along straight lines. There are further advantages, however, including: simple determination of the asymmetric domain (making for easy construction of the finite element mesh), well-behaved metric throughout the asymmetric domain, and straightforward enforcement of crystal symmetries. Due to these appealing properties, Rodrigues parameter space is finding more widespread use in the representation and analysis of preferred orientation (see, for example, Frank, 1988;Becker and Panchanadeeswaran, 1989;Randle, 1990;Heinz and Neumann, 1991;Neumann, 1991;Morawiec, 1995;Morawiec and Field, 1996;Morawiec, 1997;Adams and Olson, 1998). In fact, Morawiec and Field (1996) discuss in detail the geometric properties of Rodrigues space, noting the potential advantages of this space in the context of pole figure inversion. Using piecewise interpolants circumvents the difficult task of determining global interpolants in the Rodrigues parameter space that satisfy the symmetry constraints a priori.
Clearly, this representation of the ODF is based on a finite number of parameters (the number of degrees of freedom corresponds to the number of unconstrained nodal values in the mesh), so the first issue of indeterminacy is not a concern. A non-negativity constraint helps to deal with the second issue of indeterminacy -that associated with the inversion center and pole figure data having antipodal symmetry. The constraints are not always necessary and in some cases valid ODFs are obtained without their activation.
Before summarizing mathematical facts related to pole figures and orientations in Section 3, we briefly discuss in Section 2 the use of finite elements. Section 4 contains a description of a process by which the pole density may be determined given a finite element representation of the ODF in Rodrigues space. We then outline in Section 5 the formation of systems of equations for ODF determination from pole figure data. Example applications are also given in this section. Neutron diffraction methods were used to collect the data for the example applications, and we examine materials with both hexagonal and cubic crystal symmetry.

THE USE OF FINITE ELEMENTS
Volume, surface, and line finite elements are all used in the various aspects of the pole figure formation and inversion algorithms. Volume elements enter through the discretization of the Rodrigues parameterization of orientation space. Thus the basis functions for the ODF have local support, in contrast with the global support basis functions used in harmonic series expansion methods. As the regions of interest in Rodrigues space have planar boundaries but nontrivial geometric shapes, finite elements are well suited to the task of discretizing the space. The work of Dawson (1995, 2000) provides a thorough exposition of the use of finite elements over Rodrigues space, including a discussion of the symmetry conditions on the boundaries of the discretized region.
For the representation of pole figures, the unit hemisphere is discretized by surface finite elements. By discretizing in spherical coordinates, quadrilateral elements may be used to exactly map the hemisphere and satisfy antipodal symmetry conditions. Straight paths through Rodrigues space, corresponding to pole figure projections paths, are discretized using line finite elements with an individual line element spanning each volume element along the projection path. Integrals over the discretizations in various dimensions are all accomplished using Gaussian quadrature, taking into account the metrical properties of the spaces. As the results from pole figure path integrations are typically evaluated for a given finite element mesh and then stored, highly accurate quadrature may be used without concern for computation time. Whether the finite elements reside in a one, two, or three dimensional space, standard quadrature and mesh refinement methods may be applied.
With the ODF represented using finite elements over Rodrigues space, the pole figure inversion algorithm entails solving for nodal degrees of freedom. Thus, using the common terminology, the inversion method is ''direct'' in that the unknown quantities are values of the ODF at particular orientations. For another example of a direct pole figure inversion method, see the work of Schaeben and Siemes (1996), in which the pole figures and the ODF, represented in Euler angle space, are interpolated using piece-wise constant functions. The main contribution of the work presented here relates to the use of finite elements over Rodrigues space as a central component of the pole figure inversion method.

ORIENTATIONS AND POLE FIGURES
The orientation of a crystal lattice may be specified through the orientations of a set of crystallographic vectors, cˆi, that form a right-handed orthonormal basis. A proper rotation C then relates the lattice orientation to a reference orientation, with crystallographic vectors c i , through cˆj ¼ CÁc j or C ¼ cˆk c k . Parameterizing the operator C by its angle and axis of rotation, one obtains the following expression for C: C ¼ Rð, nÞ ¼ n n þ cosðÞðI À n nÞ þ sinðÞ skewðnÞ: ð1Þ The notation skew(n) represents the skew tensor defined so that skew(n) Á a ¼ n Â a for all a. 1 In the Rodrigues parameterization, the axis n is scaled by the angle of rotation according to resulting in the Rodrigues vector such that R R R r ðrÞ ¼ R(, n). Given two for the rotation C b followed by C a obeys the composition rule (Morawiec and Field, 1996) This relation is used in the development to follow and is intimately connected with geodesics being straight lines in Rodrigues space. The pole density function over sample directions y and diffraction plane normals h will be denoted P(h, y) with both h and y being unit vectors. For a particular h, the pole figure takes values at locations y on the unit sphere with antipodal symmetry. Pole figures may thus be represented on the upper hemisphere, S 2þ & < 3 , with locations on opposite sides of the equator taking the same value. Under ideal diffraction, as is assumed in classical texture analysis (Bunge, 1999), the pole figure at y receives contributions from the orientations for which the mapping C brings AE h or any symmetric equivalent into alignment with y. Given a crystal symmetry group H, symmetric equivalents of h are h (i) ¼ H (i) Á h with H ðiÞ 2 H. As H always contains the identity operator, h ¼ h ( j) for one of the j. For a symmetry group with n H operators, we have in which each of the P (h, y) represents a path integral. Denoting the image of h under the mapping C by hˆ¼ C Á h, in which A is the ODF, the probability density function for orientations defined on & SOð3Þ. Integration occurs over the set of orientations fC 2 jC Á h ¼ yg, denoted hˆk y, that bring h into coincidence with y. Note that Eq. (5) may be adjusted by a constant multiplying the integral in order to account for the chosen normalization methods for pole figures and ODFs. Given an ODF which satisfies the crystal symmetries, Eq. (4) collapses to two terms corresponding to AE h and no longer involves a summation over the symmetry operators. This will be achieved by defining A over a domain containing no symmetrically equivalent orientations on its interior, sometimes called an asymmetric domain. The asymmetric domain composed of orientations with smallest angle of rotation (closest to the origin for Rodrigues parameter space) is often called the fundamental region. Let fr denote the fundamental region and let @ fr , denote its boundary.
All orientations mapping h into y may be written as the operation C ¼ R R Rð , n Þ, which rotates h into y through the shortest possible angle, followed by a rotation C ¼ R R Rð , yÞ about y. To reduce the number of integrals in Eq. (4), we consider Eq. (5) in more detail. Let so that parameterizes a smooth curve in SO(3) for 0 2. We assume that A is piecewise continuous (integrable) and that A satisfies the symmetry relations of symmetry group H. The symmetry relations require A (C Á H (i) ) ¼ A(C) for all C 2 SOð3Þ and for all H ðiÞ 2 H. Given the above parameterization, Eq. (5) becomes with jdC/d j a constant in the geodesic parameterization. Now for Notice that ½Rð , n Þ Á H ðiÞ À1 Á h ðiÞ ¼ Rð , n Þ Á h ¼ y: As the members of H are isometries on SO(3), jdC ðiÞ ð Þ=dj ¼ jdCð Þ=d j and, using the fact that A satisfies the symmetry relations, Thus Eq. (4) reduces to given the appropriate conditions on A. Furthermore, if there exists an H ðiÞ 2 H such that h . This will be true if h is perpendicular to the axis of any binary symmetry. See Morawiec and Pospiech (1989) for a discussion of this and related properties of the pole density integration paths through orientation space. Despite the reduction in the number of terms in Eq. (4), crystal symmetries still play a crucial role in our algorithm for evaluating the path integrals through orientation space. The symmetry operators will be used to map the current integration segment into the fundamental region of Rodrigues space. This is equivalent to pre-applying a new symmetry operator every time the path crosses a boundary of the fundamental region (@ fr ). It is also equivalent to examining all of the paths in Eq. (4) and integrating over the portions falling inside fr .
FIGURE 1 Paths of the integrals P (h, y) and P (À h, y) through the cubic fundamental region of Rodrigues space corresponding to {123} diffraction planes for y parallel to (1, 0, 1). In the next section, we describe an algorithm for evaluating the pole density in terms of nodal degrees of freedom in the finite element discretization of the ODF. Figures 1 and 2 show illustrations of example paths through the fundamental region of Rodrigues space for materials with cubic crystal symmetry. In both figures, diffraction is from {123} planes. Figure 1 shows paths for a single y location and Fig. 2 shows paths for a collection of pole figure locations. The segments are straight and do not terminate inside the region. The paths are connected at the boundaries of fr by jumps according to the symmetry operators mapping opposing surfaces of the fundamental region into each other. These graphical representations of the paths help to illuminate the issues discussed in the following section.

PATH INTEGRALS IN RODRIGUES SPACE
The relationship between the ODF and a particular pole figure contribution, P (h, y), is central to both pole figure formation and pole figure inversion. This section develops that relationship for the case of the ODF represented using finite elements over Rodrigues space. For each P (h, y) of interest, we evaluate the path integral with the degrees of freedom of A as unknowns. The resulting expression for P (h, y) is a linear combination of the degrees of freedom of A, that is of the finite element nodal values. With the values P (h, y) expressed in terms of the degrees of freedom of A, we obtain a system of equations for pole figure values given any ODF represented using the chosen finite element basis functions. The information contained in this system of equations may then also be used for pole figure inversion.
The method employed here relies on parameterization of the path integrals through Rodrigues space, and the parameterization of the paths is discussed in Section 4.1. As an alternative, the integral in Eq. (5) may be evaluated without regard to the details of the chosen orientation space. Numerical quadrature could be used, with each quadrature point being mapped into the chosen orientation space and evaluated in terms of the coefficients in the representation of A. Such a method is conceptually simpler than explicitly parameterizing integration paths, and would work well for smooth basis functions. However, we choose to represent the ODF using finite elements in Rodrigues space, and the associated basis functions have derivatives which are discontinuous across element boundaries. To achieve greater accuracy in evaluating the integrals, we explicitly parameterize the paths of the integrals. Given the path parameterization, it is possible to match the end-points of the line finite elements along the integration path to the boundaries of the volume finite elements covering the fundamental region. As a result, the variation of A over each line finite element along the integration path is smooth and accuracy is improved.

Line Integral Parameterization in Rodrigues Space
Many of the results concerning the path integrals through the fundamental region of Rodrigues space are not particular to the finite element discretization, although they are used to advantage in the finite element algorithm. In this section, we review several key results.
To begin, let r denote the Rodrigues parameterization of C so that The parameter corresponds to Euclidean distance in Rodrigues space along the direction y. At any given location over the course of the path integration, one may use the crystal symmetries to map the path back into the fundamental region, fr , so that for i such that C (i) ( ) is in fr . If the Rodrigues parameterizations of the quantities in Eq. (12) exist individually, then where Note that one can form the Rodrigues parameterization of C (i) ( ) even if its parts do not individually have Rodrigues parameterizations due to the angles of the rotations involved.
Assuming for the purpose of illustration that r exists, we may use the composition rule for Rodrigues vectors (Eq. 3) and, from geometry, y Â (h Â y) ¼ h À (h Á y)y, to obtain (Morawiec and Pospiech, 1989;Kumar and Dawson, 2000) As neither of the vectorial portions of this expression depend on , this represents a straight line in Rodrigues space. It is in fact true that all geodesics are straight lines in Rodrigues space (Frank, 1988;Morawiec and Field, 1996). Lines through Rodrigues space parameterized by the form may always be returned to this form after composition with any constant Rodrigues vector, including the symmetry operators. In working with this form we will always have a 0 > 0 and will remain in a range such that the expression is bounded. After composition with a symmetry operator, both the reference location and the direction of the line generally change.
In our implementation, the integral is performed by tracing the path through the fundamental region, fr , and applying the appropriate symmetry operator when the path reaches @ fr in order to keep the parameterized path inside of fr . Alternatively, one could find all of the images of the integration path under the symmetry group and integrate over the segments that pass through the fundamental region. At the start of the integration, the reference location of a parameterized line is shifted into the fundamental region by selecting the H (i) that maps C into the fundamental region. Each time the path is mapped across @ fr , the parameterization in Eq. (16) is updated to reposition the reference location at the image point on @ fr . Throughout the reparameterizations, the argument of r( ) in Eq. (16) always describes a rotation of angle 2 arctan( ) about y. This fact allows one to relate lengths in the Rodrigues space to rotation angles about y, and thus to accumulate the contributions to the path integral.

Numerical Integration
Given the above results and a finite element discretization of fr in Rodrigues space (for a particular symmetry group), integration proceeds by tracking the integration path through each volume finite element in turn. Linear interpolation elements (4-noded tetrahedra) are used to discretize the orientation space. Given the linear interpolation, extreme values of A occur at nodes and constraining all nodes to have non-negative values of A enforces positivity over the whole domain. Linear elements also possess the advantage that inversion to map into the parent element space takes lines in the Rodrigues space into lines in the parent element space, simplifying the procedures for determining intersections of the integration paths with element boundaries. By considering the passage of the integration path through each finite element in turn, the integration procedure automatically adapts to the coarseness of the Rodrigues space finite element mesh.
Within a given finite element, A is interpolated using the element shape functions and the nodal point values associated with that element. The vector of all independent nodal values is represented by {A np }, and contains only non-redundant degrees of freedom (nodal points on the boundary are constrained to have the same value as nodes of symmetrically equivalent orientation on opposing boundaries). Given a pole figure location y, each pair of integrals corresponding to h and À h in Eq. (10) contributes to a system vector which maps the nodal point values of A into the value of the pole figure at y according to P A ðh, yÞ ¼ fmðh, yÞg T fA np g ð 17Þ with the notation on the right hand side indicating a vector dot product. To evaluate the contribution of one of the two terms in Eq. (10) to {m (h, y)}, integration over 2 follows the procedure: . Initialize the line parameterization so that it begins in fr . . Find the location in the finite element mesh.
. Move to the surface of the first element and start integration from that surface. . Integrate element by element to an accumulated rotation of 2, which brings the path back to the starting location (or to a symmetrically equivalent location if integration began on @ fr ). On crossing between volume elements: - Step across the interface between the volume elements into the new volume element. If this involves crossing @ fr , reparameterize so that the line is based at the boundary to be crossed and then pre-operate by the symmetry element that maps across the mesh to the opposite boundary. If a boundary of fr has not been crossed, check on the condition of the parameterization and reparameterize as necessary. -Find the angle of rotation about y needed to traverse the current element.
(a) Perform the (linear) isoparametric element inversion into the parent finite element space. (b) Determine the exit location from the element and the distance across the element along the integration path.
-Integrate in space using high order quadrature.
(a) At each quadrature point, use the element shape functions to find the contributions from the nodal values to A at the quadrature point. This involves mapping for the quadrature point to r using ¼ tan /2 and Eq. (16). (b) Assemble contributions to the system vector {m(h, y)}.
-Accumulate the contribution to the angle of rotation about y and check to see if the path integral is done. Figure 3 graphically depicts the passage of the path integral through a single finite element, and shows the quadrature point locations for an eight-point quadrature rule (Zienkiewicz, 1977).
Examples of the contributions from nodal values of A to P A (h, y) are shown in Figs. 4 and 5. The figures graphically depict the values in the m vectors that determine the given P A (h, y). The paths are shown as tubes through fr and the spheres drawn at the nodes are scaled according to the values of the contributions with the same scaling for both figures. The edges of finite elements through which the paths pass are also shown. For purposes of illustration, the values are those for the nonreduced nodal degrees of freedom. In forming the nodal contributions in the actual {m(h, y)}, values from symmetrically equivalent nodes are added. The integration path in Fig. 4 passes through the origin of the Rodrigues space and we see that the node at the origin thus has a relatively large value for its contribution to the pole figure location. As the h in Fig. 5 is aligned with a four-fold symmetry axis in the cubic sym-  metry group, the path overlaps itself four times in the fundamental region, reaching the starting point again after every /2 rotation. Some reduction in computational work would be achieved by taking advantage of such cases.

METHODS FOR POLE FIGURE FORMATION AND INVERSION
Equipped with the methods of the previous section, we are ready to form and invert pole figures. Eq. (17) may be used directly to obtain values at discrete locations on a pole figure from {A np }. For a single pole figure for which we desire values at a finite number of locations {y 1 ,. . . , y q }, the system matrix giving the values may be written where the A np j are the entries in {A np } for the given set of finite element basis functions. Each entry M ij in the system matrix corresponds to entry j out of {m (h, y i )}. The system matrix is generally not square. As an alternative, with pole figures discretized using finite elements, pole figure nodal values may be determined so that the pole figure best matches the values projected from the ODF over the whole of the surface of the hemisphere (S 2þ ). In this latter method, Eq. (17) is used at quadrature points in the mesh on S 2þ . Both approaches are used below for reconstructing pole figures following inversion of pole figure data to obtain the ODF.
We present two methods of pole figure inversion, one for discrete pole figure data and one for pole figure data defined over the whole of S 2þ . In the first method, we form the ODF by minimizing the difference between pole figure data at discrete locations and the values determined from the ODF using Eq. (17). With the nodal values of the ODF constrained to be non-negative and the objective function equal to the sum of squared differences between data and reconstructed pole figure values, a bound-constrained least squares problem is obtained. The second method also results in a bound-constrained least squares problem, but is based on finding the ODF that best matches pole figure functions defined over S 2þ . This case involves a surface integral over S 2þ , and, when evaluated numerically, results in the use of Eq. (17) at quadrature points in the mesh on S 2þ . Both of these two methods result in systems of equations that are linear in the unknown variables.
As input pole figure data must be scaled consistently, all input pole figure data are normalized so that integrating the data over S 2þ yields the surface area of S 2þ . Thus a uniform distribution on the hemisphere would have a value of unity and pole figure values are reported in multiples of a uniform distribution (MUD). Integration to scale discrete data is accomplished by first matching an interpolated pole figure function to the discrete data. Alternatively, one may scale discrete pole figure data using experimental measurements of the diffraction properties of the crystallographic planes.
Through Eq. (5), pole figures with uniform values of unity lead to an ODF also having a uniform value of unity. To be consistent with the pole figure normalization and Eq. (5), ODFs may be normalized so that R fr Ad ¼ R fr d (approximately 0.41 for materials obeying the cubic symmetry group). This integral condition may be added as an additional constraint in the least squares problem for determining the ODF, and this approach is used in one of the examples to follow. More commonly, A is normalized so that R fr Ad ¼ 1: Thus if A is uniform it has a value of A u ¼ 1=ð R fr dÞ. In the results presented here, all ODFs are plotted after being scaled so that R fr Ad ¼ 1: ODFs are normalized only after the reconstructed pole figures have been obtained.

Example Application -Pole Figure Inversion by Least Squares on Discrete Data
At a discrete set of locations y i and for families of diffraction planes with normals h () , let the pole figure data values be P d (h () , y i ). For simplicity of notation, assume that the y i are the same on all of the h () pole figures, as is the case in the example to follow. At the same pole figure locations, the pole density projected from the discretized ODF is Each row in the system of equations for the bound-constrained least squares problem thus consists of the vector {m(h () , y i )} and the number of unknowns equals the length of the {A np } vector.
For the problem to have a unique solution, there must be at least as many independent data measurements on the pole figures as there are independent degrees of freedom in the finite element discretization of A. The problem may be solved using standard methods for boundconstrained least squares problems. Given that the entries in the {m(h, y i )} are non-negative and the values in {A np } are also required to be non-negative, all of the P A (h, y i ) in Eq. (18) are necessarily non-negative. However, the least squares solution for {A np } from given pole figure data must be constrained in order to result in nonnegative values in {A np }.
As an example application of pole figure inversion using the method described above, we examine neutron diffraction data from the titanium alloy Ti-6A1-4V. The alloy is primarily at room temperature and we examine data only from the phase, which obeys the hexagonal symmetry group. Given the sharp transverse texture (c-axes aligned with the transverse direction), it is likely that the plate from which the material is taken was rolled near the transus temperature (Peters and Luetjering, 1980). The sharpness of the texture component partially motivates the use of data from this material, as many degrees of freedom are required in the vicinity of the texture component in order to capture the gradients in the ODF. Four pole figures were measured at the neutron diffraction facility in Chalk River, Canada. Data from these measurements are shown in Fig. 6. These and all other pole figures are displayed using the equal-area projection. Figure 7 depicts an ODF obtained by inversion of this pole figure data. The finite element mesh in this figure has 388 independent degrees of freedom. There are 4080 total measurements on the four pole figures, although not all of these are independent given the antipodal symmetry. Maximum and minimum singular values for the least squares system were approximately 4.0 and 0.08 respectively, so the system is well conditioned and there are no zero modes. That is, there are no modes in the shape of the ODF that do not directly influence the reconstructed values on the pole figures.
At the solution to the constrained problem, 278 of the constraints are active. Thus A has a value of zero at 278 independent nodes. The  Fig. 8, along with the differences between the reconstructed pole figures and the input data. For ease of comparison, scales are common between pole figures of the same type. In all of the results, including those in the following section, plots for different diffraction planes are not drawn on the same scale.
Pole figure inversion to obtain an ODF is possible with just one pole figure, and the results of this process using only the f2 1 1 1 12g pole figure data are shown in Fig. 9. The degree to which the ODF captures the f2 1 1 1 12g pole figure is slightly improved, but overall fidelity is, as one would expect, degraded. Figure 10 depicts an ODF obtained by pole figure inversion, using all four pole figures, onto a finer mesh. The system of equations, which now has 3080 independent degrees of freedom, remains well conditioned. At the solution, 2476 of the constraints are active. Figure 11 shows reconstructed and difference pole figures and we see FIGURE 8 Pole figures from the ODF in Fig. 7 and differences from input data. Pole figures are plotted on the same scale as in Figs. 6,9,and 11; differences are plotted on the same scale as in Figs. 9 and 11. that the ODF on the finer mesh more accurately captures the data. The residual for the new ODF is 0.61 times that of the ODF depicted in Fig.  7. Given the number of data measurements on each of the pole figures and the method used in this example, the more finely discretized ODF possesses too many degrees of freedom for its nodal values to be determined from fewer than four pole figures. As a check of these systems of equations and of the solution procedure, inversion of uniform pole figures results in ODFs that are uniform to ten significant figures. This method also has the property that accuracy is preserved after repeated projections and inversions. When the reconstructed pole figures are used as input data to the pole figure inversion algorithm, the new ODF is the same as that used to create the reconstructed pole figures. In a numerical verification of this fact, the differences in the nodal values of the ODFs are all less than 1.5 Â 10 À13 in magnitude. All of the systems considered here are well conditioned and in none of them does the ODF have modes which do not directly influence the values of the reconstructed pole figures.

Example Application -Pole Figure Inversion by Surface Integral Minimization
If the pole figure data is represented as a function over the surface S 2þ then the ODF may be determined by a surface integral minimization technique. For a set of pole figures with diffraction plane normals h () , the ODF minimizes  In all of the examples that follow, the P d h ðÞ are represented using an interpolation over the same finite element mesh of S 2þ as that which is used to numerically evaluate the surface integrals. These meshes are actually constructed in a spherical coordinate system, with the appropriate Jacobian for the mapping onto S 2þ included in the evaluation of the integrals. Nodal degrees of freedom on the surface mesh are constrained to obey the antipodal symmetry of pole figures. The unconstrained solution for {A np } solves where fm , q g fmðh ðÞ , y q Þg and the w q are weights at the quadrature points y q , including all surface element Jacobian contributions. Again, the bound-constrained problem may be solved using standard methods. We report the integrated relative error over each pole figure as and integrated error over each surface element according to for element S e . The examination of neutron diffraction data from an aluminum alloy serves as an example application of the surface integral minimization method. The material has a face-centered cubic crystal structure and thus obeys the cubic crystal symmetry group. Input data pole figures are shown in Fig. 12, with the data represented on a finite element mesh in which nodes are no more than 6 apart. Nodal values on the pole figures in Fig. 12 are obtained by finding the set that best matches the data obtained from the neutron diffraction measurements. The material exhibits a mixed rolling and cube texture, with the maximum ODF intensity for the cube component being greater than that for the rolling component.
FIGURE 13 ODF with 76 independent degrees of freedom, over the fundamental region for the cubic symmetry group. Values are in the range 0-9.5 and, as in Fig. 15, the level set is drawn at 2.55. The scale is the same as that in Fig. 15.   Fig. 13 has integrated relative error values (Eq. (21)) less than 0.103, 0.142 and 0.111, respectively. While the pole figure meshes each contain 526 independent degrees of freedom, the surface integrals for obtaining the ODF are evaluated with 2 Â 2 quadrature over each quadrilateral element, for a total of 2112 quadrature points over each pole figure mesh. Thus there are more than enough contributions to the system of 76 equations and the system is non-singular. The system has minimum and maximum eigenvalues of approximately 0.00143 and 0.281 so that, as in the examples in Section 5.1, the ODF does not possess modes that fail to contribute to the values of the pole figures.
Normalization of the pole figure data is quite accurate given that they are represented over the whole of S 2þ . It is, therefore, not surprising that the addition of the constraint R fr Ad ¼ R fr d to the least squares problem for ODF determination results in only minor changes in the resulting ODF. The norm of the difference between the integralconstrained ODF and the ODF shown in Fig. 13 is two orders of magnitude smaller than R fr d. To three significant figures, the integrated relative error values for the reconstructed pole figures are the same as those depicted in Fig. 14. Note that the ODF from the least square problem with the integral constraint is not the same as that obtained by normalizing the ODF obtained without the integral constraint. But for pole figures that are well normalized the two ODFs are not remarkably different. FIGURE 16 Pole figures from the ODF in Fig. 15 and elemental errors. Pole figures are plotting using the same scale as in Fig. 12 and element errors are plots using the same scale as in Fig. 14. FIGURE 15 ODF over the cubic fundamental region with 600 independent degrees of freedom. Values are in the range 0-12.5, though the scale is the same as that in Fig. 13. As in Fig. 13, the level set is drawn at 2.55.
An ODF computed on a mesh with 600 independent degrees of freedom yields integrated relative errors less than 0.0832, 0.0993, and 0.0907 for the {111}, {200}, and {220} pole figures, respectively. This more finely discretized ODF is depicted in Fig. 15, and Fig. 16 contains the reconstructed pole figures and the element errors. The numerical conditioning of the system of equations for this more finely discretized ODF is nearly as good as the conditioning of the system with 76 independent degrees of freedom.
As in the example in Section 5.1, the additional degrees of freedom in the ODF result in higher peak values and improved sharpness in the reconstructed pole figures. While ODF mesh refinement gives substantial improvement in the reconstructed pole figures for the titanium data (Figs. 8,11), the reconstructed pole figures in the current example do not approach the data pole figures as rapidly. This may be attributed to the noisiness of the data in Fig. 12. Neither smoothing nor sample symmetry operators have been applied to the data. The process of matching continuous pole figure functions to the neutron diffraction data may have exacerbated the fluctuations, but the fluctuations are apparent in the data obtained from the neutron diffraction facility.

CONCLUSIONS
The developments presented here allow one to obtain a finite element representation of the ODF over Rodrigues space directly from pole figure data. Input data may be represented over the surface of the pole figure or at discrete measurement locations. Using algorithms that exploit the properties of geodesics in Rodrigues space, namely that they are straight, path integrals for pole figure formation and inversion are evaluated accurately. By solving bound-constrained minimization problems, one obtains ODFs that are non-negative.
One can imagine several extensions of the current work, many of them drawing from methods and issues discussed in the literature. Adaptive mesh refinement in Rodrigues space might be used to improve the degree to which the ODF captures pole figure data. As in the wavelet representation work of Schaeben et al. (2001a), extra degrees of freedom might be added in order to capture sharp gradients in the density function. The potential utility of such localized mesh refinement is illustrated by the titanium alloy texture in Section 5.1. With the addition of degrees of freedom to the finite element representation, determinacy of the pole figure inversion system becomes an issue and various methods could be used to eliminate zero modes and to condition the problem. Potential heuristics might be similar to the entropy method employed in the work of Schaeben and Siemes (1996) or might involve penalization of gradients in the ODF through a diffusive contribution to the system. Methods that control gradients in the ODF would also be helpful in combating the effects of noisy input data such as that in Section 5.2. Other extensions could include consideration of data from incomplete pole figures or treatment of pole figures with overlapping contributions from multiphase materials. Diffraction experiments in which individual measurements are collected over sizable regions on the pole figure, such as those from time of flight neutron detectors, could be handled by associating data with surface integrals over appropriate regions on the pole figure.
Improvements in the efficiency of the implementation are also possible. For h that are aligned with one of the axes of the symmetry group, it would be possible to reduce the angular distance over which the path integral is evaluated. A more substantial time savings would be achieved by using sparse methods for the bound-constrained minimization problems. The direct methods currently in use work well for small problems, but storage and computational demands increase rapidly with system size. Standard packages for large sparse optimization with constraints, such as those discussed by Conn et al. (1992) or Murtagh and Saunders (1987), would considerably extend the size of problems solvable with reasonable effort on a given computational platform.