An inverse problem for Voronoi diagrams: A simplified model of non‐destructive testing with ultrasonic arrays

In this paper, we study the inverse problem of recovering the spatially varying material properties of a solid polycrystalline object from ultrasonic travel time measurements taken between pairs of points lying on the domain boundary. We consider a medium of constant density in which the orientation of the material's lattice structure varies in a piecewise constant manner, generating locally anisotropic regions in which the wave speed varies according to the incident wave direction and the material's known slowness curve. This particular problem is inspired by current challenges faced by the ultrasonic non‐destructive testing of polycrystalline solids. We model the geometry of the material using Voronoi tessellations and study two simplified inverse problems where we ignore wave refraction. In the first problem, the Voronoi geometry itself and the orientations associated to each region are unknowns. We solve this nonsmooth, nonconvex optimisation problem using a multistart non‐linear least squares method. Good reconstructions are achieved, but the method is shown to be sensitive to the addition of noise. The second problem considers the reconstruction of the orientations on a fixed square mesh. This is a smooth optimisation problem but with a much larger number of degrees of freedom. We prove that the orientations can be determined uniquely given enough boundary measurements and provide a numerical method that is more stable with respect to the addition of noise.

carried out using a single transducer (as opposed to the arrays of sensors available in seismic and medical imaging), and so until recently, development of tomographic capabilities has been somewhat neglected. However, with the rapid uptake of ultrasonic phased arrays (which are capable of simultaneously transmitting and receiving ultrasound across multiple elements 5 ), tomographic inversion is now becoming feasible [6][7][8][9] and presents a potential solution to the challenge of imaging defects embedded in complex media. There are many cases in industry where, due to manufacturing conditions or design, heterogeneous and locally anisotropic microstructures can arise [10][11][12] and defect detection can prove difficult. One particularly challenging example of this occurs in austenitic steel welds where, due to the application of extreme thermal gradients during the manufacturing process, the primary axes of the microstructural crystals align with the direction of heat flow. 13,14 This results in a highly scattering, refractive and locally anisotropic medium where contributions to the scattered field by small defects are often obscured. However, it has been shown that a priori knowledge of such microstructures can be used to compensate for these effects and better focus the defect scattering energy to produce more reliable flaw reconstructions. [7][8][9]15 Inspired by this particular problem in the ultrasonic NDE community, we propose a new approach to extracting information on the spatially varying material properties of this class of locally anisotropic media from travel time information collected on the boundary of the object. This specific problem has previously been tackled using Markov chain Monte Carlo (MCMC) methods, [7][8][9] where an ensemble of samples estimating the posterior probability density function is produced and its associated moments are used to describe the probability distribution of the model given the observed data. These MCMC methods are generally applicable to most high-dimensional inverse problems but are usually associated with prohibitive computational costs from the iterative runs of forward models over large and complex parameter spaces, which are in some cases inefficiently sampled. 16 In this paper, we instead cast the problem as a deterministic optimisation problem, to encourage more efficient sampling of the parameter space and allow better scaling of higher dimensional problems in future works. A global, multistart, non-linear least squares method is implemented to deal with the nonconvex objective function, and the sensitivity of the method to additive noise is examined. Moreover, we introduce a new geometric inverse problem, (IP1), which we believe is a challenging, interesting mathematical problem in its own right due to its nonsmooth and nonconvex nature.

Outline of the paper
In Section 2, we state two inverse problems, (IP1) and a simplified version (IP2). These can be viewed as reduced models for the non-destructive testing of metals using ultrasound. We show how they are derived in Section 3. We study the first inverse problem (IP1) theoretically in Section 4.1 (properties of the objective function, non-uniqueness of solutions) and numerically in Section 4.3. In Sections 5.1 and 5.2, we develop the theory and a numerical method for the simplified inverse problem (IP2).

STATEMENT OF THE INVERSE PROBLEMS
In this section, we give a precise statement of the geometric inverse problems that we study in this paper. These problems are inspired by the NDE of metals using ultrasonic arrays. In Section 3, we discuss the relation between these mathematical problems and the physical NDE problem.

The Voronoi inverse problem
Let Ω ⊂ R 2 be open, bounded and convex. Given X = (x 1 , … , x n ) ∈ Ω n , the Voronoi diagram generated by X 17,18 is the collection of sets {V k (X)} n k=1 defined by Each Voronoi cell V k (X) is the intersection of a finite number of half-planes with Ω and is therefore the intersection of a convex polygon with Ω. Observe that V ∩ V k = ∅ if j ≠ k and that the collection of Voronoi cells {V k (X)} n k=1 forms a polygonal partition of Ω (up to a set of measure zero). Sometimes, we simply write V k instead of V k (X) if it is not important to emphasise the generator set X. An edge of the Voronoi diagram is a line segment of the form V ∩ V k such that length (V ∩ V k ) > 0.
We introduce the notion of an oriented Voronoi diagram, which is a Voronoi diagram {V k (X)} n k=1 with an orientation k ∈ [0, 2 ] assigned to each Voronoi cell. Let Θ = ( 1 , … , n ) be the collection of all orientations in a diagram. We say that (X, Θ) generates the oriented Voronoi diagram {(V k (X), k )} n k=1 .
Given two points a, b ∈ Ω, let a,b ∶ [0, 1] → Ω, a,b (s) = (1 − s)a + sb be the line joining a to b. Let a, b ∈ [0, 2 ) be the direction of the line a, b , defined by Below, we will restrict our attention to lines a, b that do not contain any edges of the Voronoi diagram, which means that Note that this condition does not preclude a, b from intersecting edges, only from containing edges, and it is a generic condition satisfied by 'almost every' line a, b . Let S ∶ R → (0, ∞) be continuously differentiable and /2-periodic. We define the travel time T(a, b; X, Θ) of any line a, b satisfying Equation (2) by Note that the kth term in the sum is non-zero if and only if the line a, b intersects the Voronoi cell V k , and 'most' terms in the sum are zero. We refer to S as the slowness function since, if we think of T(a, b) as a time, then S has units of inverse velocity. We can also interpret T(a, b) as the 'length' of the chord a, b of the Voronoi diagram with respect to the weight function S. Indeed, if S takes the value 1 at every point, then T(a, b) = length( a,b ). Observe that the weight S( a, b − k ) depends on the relative angle between the line a, b and the orientation k of Voronoi cell k. The physically realistic assumption that S is ( /2)-periodic, and hence -periodic, ensures that T(a, b) = T(b, a).
This equivalent definition of T(a, b) also requires a, b to satisfy Equation (2) since is not defined on the edges of the Voronoi diagram. The problem is that, while the piecewise constant function is defined almost everywhere (with respect to the two-dimensional Lebesgue measure), in Equation (4), we are integrating along a one-dimensional set.
We extend the definition of T(a, b) to those lines a, b that violate Equation (2) by defining This choice ensures that the travel time of each line a, b is as short as possible and so seems natural from a modelling point of view. Other choices are of course possible. In Section 4, we will study the following inverse problem: fix n, M ∈ N. Given pairs of boundary points In other words, we want to determine (or fit) an oriented Voronoi diagram from the travel times between boundary points, or from the 'length' of chords of the Voronoi diagram. We refer to Equation (3) as the forward problem and (IP1) as the inverse problem. A necessary (but not sufficient) condition for the given travel times {T i } M i=1 to be compatible with an oriented Voronoi diagram (i.e., a necessary condition for the infimum in [IP1] to be zero) is that We also assume that S is a nonconstant function; otherwise, T(a i , b i ; X, Θ) is independent of (X, Θ), and every oriented Voronoi diagram satisfies (IP1).

The fixed-grid inverse problem
In this section, we state a simplified version of (IP1). Let Ω = (0, L 1 ) × (0, L 2 ) be a rectangular domain. We simplify the inverse problem (IP1) by taking the Voronoi diagram {V k (X)} n k=1 to be a given square grid and by solving (IP1) just for Θ. To be precise, let h > 0, N = ⌊L ∕h⌋, j ∈ {1, 2}. For l j ∈ {1, 2, … , N j }, define x l 1 l 2 ∈ Ω by In Section 5, we will study the following inverse problem: fix h, M > 0. Given pairs of boundary points In other words, we want to determine the orientations of the grid cells from the travel times between boundary points.

DERIVATION OF THE INVERSE PROBLEMS
In this section, we derive the inverse problems (IP1), (IP2) from the more general setting of travel time tomography in anisotropic media. Our physical motivation is the determination of metal microstructure and defects from the travel times of ultrasonic waves. 5,8,9,19 We follow the notation of Benmansour et al. 20 Let Ω ⊂ R 3 be an open, bounded, convex set representing a metal sample, a ∈ Ω be the location of an ultrasound transmitter, and b ∈ Ω be the location of an ultrasound receiver. The time it takes for an ultrasound wave to travel from a to b is where the minimisation is over all curves ∶ [0, 1] → Ω (ultrasound rays) joining the transmitter a to the receiver b, and is an anisotropic metric. Here, ∶ Ω × S 2 → (0, ∞) has units of inverse velocity, and (x, ) is the slowness of a ray passing through point x in direction . The travel time tomography problem is to approximate an unknown metric from a collection of travel time measurements {t (a i , b i )}, a i , b i ∈ Ω. As observed, for example, in Benmansour et al, 20 the map  → t (a, b) is the pointwise infimum of a family of linear functions of and hence is concave. Consequently, the travel time t depends in a smooth (Lipschitz continuous) way on the metric . Unfortunately, the inverse problem (IP1) does not enjoy the same smoothness. A further challenge is that the metric is anisotropic. If we wish to solve the inverse problem numerically using a gradient method, then we need to compute an element of the superdifferential of the concave function  → t (a, b). This can be done for the isotropic case (where (x, ) is independent of ) using for example the supergradient marching algorithm of Benmansour et al, 20 which applies automatic differentiation to a fast-marching algorithm. 21 One approach to the anisotropic case could be to generalise Benmansour et al 20 to anisotropic metrics, using an anisotropic fast-marching algorithm such as those presented in other studies, 8,22,23 but we do not pursue this here.
For simplicity, we will restrict our attention to two dimensions, Ω ⊂ R 2 and ∶ Ω × S 1 → (0, ∞). We write elements of S 1 as = (cos , sin ), where ∈ [0, 2 ). Since our motivation is the non-destructive testing of metals, we assume that has the form for some orientation function : Ω → [0, 2 ] and /2-periodic slowness function S ∶ R → (0, ∞). Here, (x) is the orientation of the metal lattice at point x (we assume that the metal is single-phase, which means that the lattice type is the same at every point, and polycrystalline, which means that the orientation of the lattice can vary from point to point). Waves in metals move at different speeds in different directions, due to the preferred directions of the lattice, which is why the slowness (x, ) depends on the relative angle between the lattice orientation and the ray direction .
In this paper, we work with two different models for . In Section 4, we assume that is a piecewise constant function with respect to a Voronoi diagram {V k } n k=1 : where V k is the characteristic function of the Voronoi cell V k and k ∈ [0, 2 ]. This choice of has a strong physical motivation. The microstructure of polycrystalline metals consists of grains, which are regions of (almost) constant lattice orientation, and grains often resemble Voronoi cells. This resemblance is not coincidental; grains would form a perfect Voronoi diagram under the highly idealised cooling conditions where a molten metal starts solidifying from several points (the Voronoi generators) at exactly the same time and solidifies at exactly the same rate in every direction. Consequently, it is common to model grains as Voronoi cells. [24][25][26][27] The typical grain size in steels is smaller than the wavelength of ultrasound, and so in this paper, we use the term grain somewhat imprecisely to mean an averaged grain cluster of a size that can be resolved by ultrasound waves. For the choice of given in Equation (9), the inverse problem is to determine the unknown oriented Voronoi diagram {(V k , k )} n k=1 from a collection of known travel times {t (a i , b i )}, a i , b i ∈ Ω. In Section 5, we assume that Ω = (0, L 1 ) × (0, L 2 ) is rectangular and that is a piecewise constant function with respect to a square grid: where h > 0 is the grid size, l j ∈ {1, … , L j /h}, l 1 l 2 ∈ [0, 2 ]. This choice of is motivated by electron backscatter diffraction (EBSD) measurements, where the lattice orientation is measured at each point in a grid on a polished metal surface. Grouping together, grid cells with the same orientation result in the grains referred above. For this choice of , the inverse problem is to determine the unknown orientations { l 1 l 2 } from a collection of known travel times {t ( Note that we take the grid to be known since each grid cell corresponds to an EBSD pixel. To arrive at the inverse problems (IP1), (IP2), we make one more simplifying assumption. The travel time (or forward problem) t (a, b) defined in Equation (7) is difficult to evaluate since it involves solving a minimisation problem for the optimal ray . Instead, we simplify the forward problem by replacing the optimal ray with the straight ray a,b (s) = (1 − s)a + sb joining the transmitter a to the receiver b. Therefore, we replace the true travel time t (a, b) with the approximate travel time This approximation is equivalent to assuming that the inhomogeneity and anisotropy of the material is small. While this is not always the case for materials of interest in ultrasonic tomography problems, it leads to a challenging, new mathematical problem that is interesting in its own right. By substituting the expressions for and from Equations (8) and (9) into Equation (11), we finally arrive at the inverse problem (IP1), as we now demonstrate. Given a, b ∈ Ω, let a, b be the direction of the line joining a to b, that is, This is exactly the travel time function (3) for the inverse problem (IP1). The advantage of the straight-ray approximation is that the forward problem (12) is very easy to evaluate. The disadvantage, however, as we shall see in the following section, is that the travel time now depends in a nonsmooth way on the metric (on the generators of the Voronoi diagram). A similar calculation to Equation (12) shows that substituting the expressions for and from Equations (8) and (10) into Equation (11) gives the inverse problem (IP2). This concludes the derivation of the inverse problems.
Finally, a quick remark about the isotropic case, where (x, ) = (x) is independent of . In this case, Equation (11) is simply the Radon transform of along the line a, b . It follows that the infinite family of travel times {T (a, b) : a ∈ Ω, b ∈ Ω} uniquely determines . This is because is determined by its Fourier transform̂, which is determined by its Radon transform: recall that if L( , t) = {t + s ⟂ ∶ s ∈ R} is the line in R 2 with unit normal ∈ S 1 and signed distance t ∈ R from the origin, then (where we have extended to R 2 by zero outside Ω). Therefore, the inverse problem has a unique solution if we know the travel time between every pair of boundary points (admittedly an unrealistic requirement in practice). Computation (13) fails for anisotropic metrics, however, and so we cannot use the Radon transform approach in this paper.

SOLUTION OF THE VORONOI INVERSE PROBLEM
In this section, we study the inverse problem (IP1).

Theory
First, we introduce some notation. For n ∈ N, we let  n (Ω) ⊂ Ω n denote the set of admissible generators of a Voronoi diagram of n points: where T was defined in Section 2.1. The goal of inverse problem (IP1) is to minimise F, that is, to find the oriented Voronoi diagram (X, Θ) giving the best least squares fit of . The challenge is that F is nonconvex and nonsmooth (discontinuous) and the graph of F has flat regions where F is constant, as we illustrate in the following examples. In each example, we use the same slowness function, where c 1 > c 2 > 0 are constants. For appropriate choices of c 1 and c 2 (c 1 = 5621 ms −1 , c 2 = 540 ms −1 ) this is a good fit for the slowness function of austenitic steel.
. We show that F is not continuous at the point ) . 1∕2)).
Therefore, F is not continuous at (X 0 , Θ 0 ). Other points of discontinuity include for all s ∈ (0, 2), y ∈ (0, 1). In general, F is continuous everywhere except at nongeneric points (X, Θ) such that the Voronoi diagram generated by X has an edge contained in the ray Example 4.2 (F is not lower semicontinuous). From the previous example, we can also see that F is not lower semicontinuous. By Equation (5), Therefore, F is not lower semicontinuous at (X 0 , Θ 0 ). (T on the other hand is lower semicontinuous, thanks to the choice (5).) (i) Number of constraints M is less than the number of degrees of freedom 3n. This is the case in Example 4.1, where we have M = 1 rays and 3n = 6 degrees of freedom. For all t ∈ (0, 2), ∈ [0, 2 ], F(X t , (3 ∕2, )) = 0 (since the ray a 1 ,b 1 only intersects the cell V 1 (X t )) and therefore (X t , (3 /2, )) are all minimisers (and there are many others). Note also that F(X t , Θ 0 ) =constant >0 for all t ∈ (− 1, 0) and so the energy landscape can also be flat in regions where it is not minimal. (ii) Lack of injectivity of S. For example, in Example 4.1, where S is the /2-periodic function defined in Equation (15), other minimisers include F(X t , (0, )), F(X t , ( /2, )) and F(X t , ( , )) for all t ∈ (0, 2), ∈ [0, 2 ]. (iii) Small grains. 'Small' Voronoi cells that are not intersected by any of the rays (a i , b i ) also give rise to a flat energy landscape and non-uniqueness of minimisers. For example, let Ω = (0, 1) × (0, 1), n = 2, M = 9 (so that the number of constraints is greater than the number of degrees of freedom), and . .
Then, F(X, (0, 2 )) = 0 for all X sufficiently close to X 0 and for all 2 ∈ [0, 2 ] since none of the rays a i ,b i intersect V 2 . Therefore, F has infinitely many minimisers.
Remark 4.4 (Relaxation and regularisation). In general, establishing the existence of global minimisers of F is challenging since F is not lower semicontinuous and its domain is not compact, and so we cannot apply the Weierstrass Extreme Value Theorem. 28, p. 3 The set  n (Ω) is not compact since you can find a sequence (X m ) m∈N ⊂  n (Ω) converging to X ∈ Ω n but with x i = x j for some i, j, that is, the Voronoi generators x m i , x m collide in the limit. The lack of compactness of  n (Ω) could be circumvented by replacing it by the compact set where the Voronoi generators are not allowed to be closer than distance > 0. The lack of lower semicontinuity of F could be circumvented by replacing it with its relaxation This is lower semicontinuous by definition and hence has a minimiser in the compact set  n (Ω) × [0, 2 ] n . It agrees with F at all generic points (X, Θ) such that the Voronoi diagram generated by X does not have an edge contained in a ray a i ,b i . Alternatively, we could regularise the objective function F. We do not take either of these approaches however in this paper. In practice, since 'most' points in the domain of F are generic, then a gradient method appears to perform fairly well at numerically minimising F, see Section 4.3.

The gradient of F
We demonstrate how to compute the gradient of F at generic points where it is differentiable. It is straightforward to compute F/ Θ. The challenge is to compute F/ X since this depends on how the lengths of the intersections a i ,b i ∩V k (X) vary as X varies. Fix k ∈ {1, … , n} and assume that a, b ∈ Ω and X ∈  n (Ω) are such that a, b intersects V k (X) in exactly two points and that neither of these points is a vertex of V k (X). Define p 1 to be the first point where a, b intersects V k (X) and p 2 to be the second point, that is, (p 1 − a) · (b − a) < (p 2 − a) · (b − a). For j ∈ {1, 2}, let p j belong to the edge E km ∶= V k ∩ V m for some m j ∈ {1, … , n} ∖ {k}, m 1 ≠ m 2 , and let = x m − x k , which is normal to edge E km . Let a, b be normal to a, b . We will show how to compute the derivative with respect to X of It is easy to see that since p 1 ∕ x l = 0 unless l ∈ {k, m 1 } and p 2 ∕ x l = 0 unless l ∈ {k, m 2 }. Therefore, we just need to compute the four derivatives We will do this for the case = 1; the case = 2 is analogous. Let q = (x k + x m 1 )∕2, which lies on the line containing the edge E km 1 . Since p 1 lies on the edge E km 1 and on the ray a, b , we have where we have taken 1 , a, b to be column vectors. By differentiating this linear system for p 1 , we see that in order to compute p 1 / x k and p 1 ∕ x m 1 , it is sufficient to compute These can be computed easily from the identities It is now straightforward to combine these expression to find an (admittedly rather lengthy) expression for F/ X, which we omit due to space restrictions. In the following section, we will use the gradient of F for numerical inversion.

Numerical inversion
In this section, we solve the inverse problem (IP1) numerically. In all of our numerical experiments we consider a square metal sample with side length L = 0.06, that is, Ω = (0, 0.06) × (0, 0.06). We assume that a series of m ultrasound transducer elements are placed on each side of the square. These elements are indexed from 1 to 4 m starting at the bottom left corner and continuing clockwise. Each element acts as both a transmitter and receiver and is paired with every other element excluding those on the same edge of the square (since we do not want ultrasound rays a, b to lie along edges of the square). This gives rise to a set of M = 6 m 2 transmitter-receiver pairs In this section, we consider four distinct numerical experiments, which were chosen to illustrate different aspects of the optimisation problem. In Examples 4.5 and 4.6, m = 2 transducer elements are placed on each edge of the sample, giving rise to M = 24 transmitter-receiver pairs (a i , b i ) with coordinates where k 1 = 0.015, k 2 = 0.045. The number of cells n to be recovered and the number of initial guesses used in the multistart optimisation framework are then varied. In Example 4.7, a more complex geometry with n = 10 cells is studied, and the number of elements m placed on each side of the boundary is necessarily increased to constrain the problem. These three examples demonstrate the method's ability to recover the Voronoi diagram in the absence of system noise. The effect of noise is then examined in Example 4.8, and the subsequent findings motivate the need to consider the alternative inverse problem (IP2).

Example 4.5.
In this example, the objective is to recover the oriented Voronoi diagram with n = 3 generators shown in Figure 1A Figure 1B), where the slowness function is as defined in Equation (15). In this example and the following two examples, we use noise-free synthetic data for the travel times, that is, we take T i to be the travel times (3) for the true oriented Voronoi diagram. To solve the optimisation problem (IP1), MATLAB's lsqnonlin (non-linear least squares method) is used and a multistart framework is implemented to cope with the nonconvex objective function. To begin, 20 arbitrary initial guesses are selected, and we pass lsqnonlin the true gradient of F (computed in the previous section). A plot of the estimated oriented Voronoi diagram is shown in Figure 1B   Example 4.6. This example examines the effect of the number of initial guesses used in the multistart optimisation framework. This time we take n = 5 Voronoi cells, and the number of transmitter-receiver pairs is the same as in the previous example. Figure 2A shows the true geometry, and Figure 2B-F shows the Voronoi diagrams recovered using a multistart adaptation of MATLAB's lsqnonlin function with {1, 25, 50, 100, 150} random initial guesses, respectively (only the best result is shown in each case). Visually, Figure 2D-F appears to reconstruct the true Voronoi diagram very well. Table 1 records the number of initial guesses and the associated L 2 -errors || true − computed || L 2 (Ω) . It can be observed that as the number of initial guesses is increased, the L 2 -errors decrease.
Example 4.7. In Examples 4.5 and 4.6, it was shown that (X, ) can be accurately recovered for n = 3, 5 with just m = 2 elements placed on each side of the sample. Here, we consider a more complex geometry where the number of cells is increased to n = 10. Since this example has 3n = 30 degrees of freedom, more transmitter and receiver pairs are required to recover (X, Θ) uniquely. We study the cases where m = 2, 3, 4, 5 elements are placed on each side of the sample (i.e., M = 24, 54, 96, 150 transmitter-receiver pairs, respectively). The optimisation problem is solved using MATLAB's lsqnonlin function within a multistart framework, implemented with 50 random initial guesses. Figure 3A depicts the true oriented Voronoi diagram, and Figure 3B-E depicts the computed oriented Voronoi diagrams for the cases M = 24, 54, 96, 150, respectively. The transducer elements are depicted in Figure 3B-E as green stars. It can be seen in Figure 3B that for m = 2 (where the number of unknowns 3n is greater than the number of constraints M = 6m 2 ), the algorithm fails to recover the true Voronoi diagram. Figure 3C-D shows better recovery but still fail to correctly capture some aspects of the geometry. In this case, we require m = 5 (M = 150) to successfully recover a Voronoi diagram with a satisfactory L 2 -error (see Table 2).   Figures 1B, 2F and 3E), but we now study the cases where noise is added to the travel time data in the following sense. We define noisy travel times T(a, b) = T(a, b)(1 + ), where ∼ N(0, 1) is a normally distributed random variable with mean 0 and variance 1, and dictates the level of noise added (the cases where = 0, 0.01, 0.02, 0.05 are studied here). In Figure 4, the first row shows the true oriented Voronoi diagrams; the second shows the computed Voronoi diagrams as recovered in Examples 4.5-4.7 and the third, fourth and fifth rows show the computed Voronoi diagrams in the case where 1%, 2% and 5% noise has been added to the true travel times, respectively. In Table 3, the respective L 2 -errors are recorded, and it can be observed that the problem is very sensitive to the addition of noise. For example, in the case where n = 3, the L 2 -error with 5% noise is more than double the L 2 -error observed when there is no noise. To summarise, we have shown in this section that Voronoi diagrams with n = 3, 5, 10, generators can be recovered very well provided there are enough transmitter-receiver pairs and enough initial guesses for the multistart optimisation method. However, in the final example, it was shown that the inverse problem (IP1) is very sensitive to the addition of noise. This motivates the examination of the simplified inverse problem (IP2) in the next section.

SOLUTION OF THE FIXED-GRID INVERSE PROBLEM
In this section, we study the inverse problem (IP2).   guesses used in the optimisation algorithm, the number of transducer elements on each side of square and the L 2 -error with and without noise added to the true travel times

Theory
Let Ω = (0, L 1 ) × (0, L 2 ), h > 0. Define X h as in Section 2.2 to be the generators of a grid Voronoi diagram with grid size h. Let N h = ⌊L 1 ∕h⌋ · ⌊L 2 ∕h⌋ be the number of grid cells. Given M pairs of boundary points where T was defined in Section 2.1. The goal of inverse problem (IP2) is to minimise F h , that is, to find the grid orientations Θ giving the best least squares fit of . Unlike in Section 4, the objective function here is continuously differentiable (since the Voronoi diagram is fixed and since S is continuously differentiable), and its domain is compact. Therefore, the existence of a minimiser of F h is trivial. Fixing the Voronoi diagram to be a grid has a disadvantage however from the modelling point of view; in order to represent well the (approximately polygonal) grain microstructure of a metal (see Section 3), we require many more degrees of freedom when we use a grid than if we use a Voronoi diagram. Therefore, the price to pay for a smooth objective function is a larger number of degrees of freedom.
The following theorem shows that the travel times uniquely determine the grid orientations if we have a sufficient number of transmitters and receivers. Recall that each Voronoi cell in the square grid is denoted by V l 1 l 2 with the bottom left cell being V 11 , the cell above that being V 12 and with this first column of cells continuing up to V 1N 2 . Similarly, to the right of cell V 11 is cell V 21 and at the end of this first row is cell V N 1 1 . We will need to further develop our earlier notation for the following theorem. Denote the four edges of each Voronoi cell by V (1)  (0, h∕2) 11 so that each of the rays will have a distinct direction a 1 ,b 1 and a 1 ,b 2 . From Equation (12), then, Suppose that there are two orientations 11 and̄1 1 that give rise to the same travel times so that T(a 1 , b 1 ; 11 ) = T(a 1 , b 1 ;̄1 1 ) and T(a 1 , b 2 ; 11 ) = T(a 1 , b 2 ;̄1 1 ), then From the symmetry of S the first of these equations implies that a 1 ,b 1 − 11 = a 1 ,b 1 −̄1 1 or a 1 ,b 1 − 11 = ∕2 − ( a 1 ,b 1 −̄1 1 ). That is 11 =̄1 1 or 11 = 2 a 1 ,b 1 − ∕2 −̄1 1 . Similarly, the second equation implies 11 =̄1 1 or Since a 1 ,b 1 ≠ a 1 ,b 2 then we must have 11 =̄1 1 . Hence, the orientation is uniquely determined in cell V 11 . Let us turn our attention now to cell V 12 lying immediately above cell V 11 . Consider now the two straight rays going from a 1 = (h, 0) ∈ V (2) 12 , which will therefore have distinct directions a 1 ,b 1 and a 1 ,b 2 . From Equation (12), then, Suppose that there are two orientations 12 and̄1 2 that give rise to the same travel times so that T(a 1 , b 1 ; 12 ) = T(a 1 , b 1 ;̄1 2 ) and T(a 1 , b 2 ; 12 ) = T(a 1 , b 2 ;̄1 2 ). Since the lengths and 11 are known, then we have As before, we can use the symmetries of S to show that 12 =̄1 2 and hence the orientation is uniquely determined in cell V 12 . It is clear that one can continue in this fashion up the first column of cells and uniquely determine each of their orientations. We turn our attention then to the determination of the orientations in the second column of cells. Consider the two straight rays which start from a 1 = (2h, 0) ∈ V (2) 21 , intersect with the pointsb 1 21 , before arriving at b 1 , b 2 ∈ Ω (3) . As above, it is clear that these two rays will have distinct directions a 1 ,b 1 and a 1 ,b 2 . From Equation (12), then, Note that given the sequential nature in which the orientations are being determined, all parameters in the latter summation in both of these equations are known. Suppose that there are two orientations 21 and̄2 1 that give rise to the same travel times so that T(a 1 , b 1 ; 21 ) = T(a 1 , b 1 ;̄2 1 ) and T(a 1 , b 2 ; 21 ) = T(a 1 , b 2 ;̄2 1 ). This implies S( a 1 ,b 1 − 21 ) = S( a 1 ,b 1 −̄2 1 ) and S( a 1 ,b 2 − 21 ) = S( a 1 ,b 2 −̄2 1 ).
As before, we can use the symmetries of S to show that 21 =̄2 1 and hence the orientation is uniquely determined in cell V 21 . We can continue in this manner up the second column; however, when we get to cell V 2N 2 , the endpoints of the rays now lie on the upper edge of cell V 1N 2 and so b 1 , b 2 ∈ V (1)  The methodology in the proof of this theorem can form the basis of a computer algorithm that, while demonstrating the correctness of the proof, shows that it is not a practical algorithm for recovering the cell orientations; a very small amount of noise in the travel time data highlights its instability. Let us consider a domain with L 1 = L 2 = 0.1, N 1 = N 2 = 10 and the orientations in each cell randomly assigned from a uniform distribution with range (0, /2). An example Voronoi diagram with the colours representing the orientations is shown in Figure 5A. The travel times T(a, b) for the rays are calculated, and then, noise is added according to T(a, b) = T(a, b)(1 + ) where ∼ N(0, 1) is a normally distributed random variable with mean 0 and variance 1, and dictates the level of noise. The noisy travel time dataT(a, b) is then used in the above algorithm to recover the orientations in each cell. In Figure 5B, the recovered orientations are shown when there is no noise added ( = 0). Comparing these two plots visually, it can be seen that there is a perfect recovery of the cell orientations. Figure 5C then shows the recovery when 0.001% noise is added ( = 0.00001) and a visual comparison shows that around 75% of the orientations are recovered. As more noise is added, then Figure 5D with 0.01% noise shows a recovery of around 50% of the orientations, Figure 5E with 0.1% noise shows a recovery of around 40%, and finally, Figure 5F with 1% noise the recovered orientation map is very poor indeed.
So it is very clear that while this approach does show that with a sufficient number of travel times one can uniquely determine the orientations of this grid Voronoi diagram, the algorithm is inherently unstable and does not present a practical solution to this inverse problem. This instability arises for a number of reasons not least of which is that it is a recursive algorithm and errors in the orientation in one cell are compounded by affecting the calculations of the cell orientations that follow. In addition, as the algorithm uses the inverse of the slowness curve S −1 ( ), then if is close to one of the turning points of S( ), then any small rounding errors are magnified. Also, when recovering the orientation of a cell that is far from the starting point of the ray, the two rays that traverse this cell are almost parallel and very close    together. They therefore sample the slowness curve S( ) in close proximity and this can again produce significant errors in the determination of the cell orientation.

Numerical inversion
In this section, we present a more practical numerical method for solving the inverse problem (IP2) that is less sensitive to noise.

Example 5.2.
We consider the same square domain Ω and slowness function S as in Example 4.5. We place m = 10 transducer elements on each side of the square (shown as green stars in Figure 6B-O) to give M = 6m 2 = 600 transmitter-receiver pairs (see Section 4.3). We define (synthetic, noise-free) travel times {T i } 600 i=1 to be the true travel times given by Equation (3) for the oriented Voronoi diagram with n = 10 cells shown in Figure 6A. Given a grid size h > 0, our goal is to minimise F h (defined in Equation 16) to recover the N h grid orientations l 1 l 2 . In this example, will illustrate the effect of the grid size. To minimise F h , we use MATLAB's lsqnonlin (non-linear least squares) function within a multistart optimisation framework, with 100 random initial guesses. Figure 6B-O shows the reconstructed grid orientations with a varying number of grid cells N h ∈ {3 2 , … , 16 2 }. As the number of grid cells N h increases, the recovered grid Voronoi diagrams more closely resemble the true Voronoi diagram shown in Figure 6A. The L 2 -errors || true − grid || L 2 (Ω) are given in Table 4, and although this shows a generally decreasing trend as the number of grid cells increases, the improvement is small. This can be attributed to the fact that as the number of grid cells increases, the number of degrees of freedom N h of the problem also increases, and it becomes more difficult to find the global minimum of (IP2). Example 5.3 (Effect of noise). In this example, we study the effect of additive noise in the travel time data. We consider noise of the form T(a, b) = T(a, b)(1+ ) (given in Section 5.1) in the context of Example 5.2 with N h = 100. Figure 7A shows the true Voronoi diagram, and Figure 7B-E shows the computed Voronoi diagrams with 0%, 1%, 2% and 5% noise ( = 0, 0.01, 0.02, 0.05). The L 2 -error for 0% noise is 1.44 × 10 −4 , and the L 2 -error for 1%, 2% and 5% noises are 1.51 × 10 −4 , 1.49 × 10 −4 , 1.92 × 10 −4 , respectively. We see that this numerical method is much more robust to noise that the exact method given in Section 5.1.

CONCLUSIONS
This paper presents two industrially relevant inverse problems from a deterministic perspective. The first, (IP1), considers the inversion of an oriented Voronoi diagram from ultrasonic travel time measurements made on its boundary. A non-trivial calculation of the gradient of the objective function with respect to the generators of the Voronoi diagram is given. By combining this with a multistart non-linear least squares optimisation method, it is shown that, given enough transmitter-receiver pairs on the boundary of the domain, we can faithfully reconstruct the locally anisotropic metric in the absence of noise. The second inverse problem, (IP2), studies reconstruction of the orientations only, over a fixed regular mesh. In this case, it is shown that given enough data measurements on the boundary, the orientations can be reconstructed exactly and uniquely. However, due to steep gradients present in the slowness curve and cumulative error propagation, the method is inherently unstable. To remedy this, we presented a numerical method that is inexact but less sensitive to noise.
Note in both cases, to avoid running a nested optimisation in order to find the shortest ray path between two boundary points, the simplifying assumption that the rays travel in straight lines has been made. Of course, refraction is important to model wave front propagation accurately in highly anisotropic media, and it would be interesting to extend our methods to include ray-bending. Another area for further investigation is the role of regularisation in solving the inverse problem; due to the nonsmooth underlying metric we wish to reconstruct, it is not immediately clear what type of regularisation should be applied.