## Mathematical Formulation

As an introduction, we begin with the derivation of surface estimation problem in two dimensions. The goal is to simultaneously estimate the interface between two materials and their densities, ¡d and ¡\. Thus we have a background with density p0 and collection of solid objects with density p1. We denote the (open) set of points in those regions as the closure of that set, the surface, as S.

The projection of a 2D signal f (x, y) produces a sinogram given by the radon transform as

-to J—to whereR0x = xcos(0) + ysin(0)isarotationandprojectionofapointx = (x, y) onto the imaging plane associated with 0. The 3D formulation is the same, except that the signal f (x, y, z) produces a collection of images. We denote the projection of the model, which includes estimates of the objects and the background, as p(s, 0). For this work we denote the angles associated with a discrete set of projections as 01,...,0N and denote the domain of each projection as S = s1; ...sM. Our strategy is to find p0, and p1 by maximizing the likelihood.

If we assume the projection measurements are corrupted by independent noise, the log likelihood of a collection of measurements for a specific shape and density estimate is the probability of those measurements conditional on the model, ln P (p(s1, 01), p(s2, 01),p(sm , 0N )|S, Po, P1)

We call the negative log likelihood the error and denote it ^data. Normally, the probability density of a measurement is parameterized by the ideal value, which gives

Edata =

where E(p^j, p, j) = — ln P(p,j, p^j) is the error associated with a particular point in the radon space, and pit j = p(sj, 0,).In the case of independent Gaussian noise, E is a quadratic, and the log likelihood results in a weighted least-squares in the radon space. For all of our results, we use a Gaussian noise model. Next we apply the object model, shown in Fig. 8.17, to the reconstruction of f. If we let g(x, y) be a binary inside-outside function on then we have the following approximation to f (x, y):

Density=Pj Density=p(

Density=Pj Density=p( Figure 8.17: The model is the interface between two densities, which are projected onto the imaging plane to create p(s, 0j).

Applying the radon transform to the model and substituting for p gives

Edata = J2 J2 E ( S, + ¡1 - ¡>] x - Sj Pij) , (8.37)

where K(Sj, 0i) is the projection of the background—it depends on the geometry of the region over which the data is taken and is independent of the surface estimate. For some applications we know that ¡}0 = 0, and the term ¡}0K is zero. The integral over O results from integrating g over the entire domain.

The proposed strategy is to alternately (i.e. separately) update the shape of the surface model and the density parameters. For the surface shape, a gradient descent minimization approach describes the deformation of the surface, with respect to an evolution parameter t, as it progressively improves its fit to the Figure 8.18: The reconstruction strategy starts with an initial surface estimate and iteratively modifies its shape and the associated density parameters to achieve a good fit to the input data.

sinogram data. The incremental change in the likelihood is dEdr = fs £ £ drtEp ^dx=I £ £E' ppj d-drdx'

where E' = dE/d p, which, for Gaussian noise, is simply the difference between p and p. Next we must formulate dp/dt, which, by the transport equation, is d% = [ft - fto] d f x - sj)dx

where n is an outward pointing surface normal and v(x) is the velocity of the surface at the point x. The derivative of Edata with respect to surface motion is therefore dE f N M

dEdata dt p N M

= [ßi - ßo] / > > E'p j, pj 8R x - sj )n(x) ■ v(x) dx. (8.40)

Note that the integral over dx and the 8 functional serve merely to associate sj in the ith scan with the appropriate x point. If the samples in each projection are sufficiently dense, we can approximate the sum over j as an integral over the image domain, and thus for every x on the surface there is a mapping back into the ith projection. We denote this point si(x). This gives a closed-form expression for the derivative of the derivative of Edata in terms of the surface velocity, dEdata = [ft - ft0] f N ei(x)n(x) • v(x)dx, (8.41)

where ei(x) = E'(1p(si(x), 0i), p(sj,(x), 0i)) is the derivative of the error associated with the point si(x) in the ith projection. The result shown in Eq. (8.41) does not make any specific assumptions about the surface shape or its representation. Thus, this equation could be mapped onto any set of shape parameters by inserting the derivative of a surface point with respect to those parameters. Of course one would have to compute the surface integral, and methods for solving such equations on parametric models (in the context of range data) are described in .

For this work we are interested in free-form deformations, where each point on the surface can move independently from the rest. If we let xt represent the velocity of a point on the surface, the gradient descent surface free-form surface Figure 8.19: The model expands or contracts based on the difference in the sinograms between the projected model and the measured data.

motion is xt = - ^Ex1 = (¡0 - ¡¡1) E ei(x)n(x). (8.42)

Thus, at a point x e S, the ith projection has the effect of causing the surface to expand or contract according to the difference between the projected model values and the measured data at the point si(x), the projection of x (Fig. 8.19). The surface motion prescribed by a collection of projections is the sum of motions from the individual projections. In the case of continuous set of angles, the surface motion at a point is proportional to the sinusoidal line integral on the error sinogram, which is e(s, 0).

### 8.6.2.1 Density Parameter Estimation

The density parameters also affect the error term in Eq. (8.37). We propose to update the estimate of the surface model iteratively, and at each iteration we re-estimate the quantities ¡0 and ¡1 in such a way that the energy Edata is minimized. Treating ^ as fixed, Eq. (8.37) has two unknowns, ¡0 and ¡1, which are computed from the following system:

In the case of a Gaussian noise model, (8.43) is a linear system. Because of variations in instrumentation, the contrast levels of images taken at different angles can vary. In such cases we estimate sets of such parameters, i.e., ¡io(0i) and p1(0i) for i = 1,..., N.

To extend the domain to higher dimensions, we have x e R™, and S c 1R™-1 and the mapping si : IR™ ^ S models the projective geometry of the imaging system (e.g. orthographic, cone beam, or fan beam). Otherwise, the formulation is the same as in 2D.

One important consideration is to model more complex models of density. If po and p1 are smooth, scalar functions defined over the space in which the surface model deforms and g is a binary function, the density model is f (x) = Ao(x) + (Mx) - Ao(x)) g(x, y). (8.44)

The first variation of the boundary is simply dX N

Note that this formulation is different from that of Yu et al. , who address the problem of reconstruction from noisy tomographic data using a single density function f with a smoothing term that interacts with a set of deformable edge models T. The edges models are surfaces, represented using level sets. In that case the variational framework for deforming T requires differentiation of f across the edge, precisely where the proposed model exhibits (intentionally) a discontinuity.