## Implementation

The level set equations are solved by finite differences on a discrete grid, i.e. a volume. This raises several important issues in the implementation. These issues are the choice of numerical approximations to the PDE, efficient and accurate schemes for representing the volume, and mechanisms for computing the sinogram-based deformation in Eq. (8.47).

where D = is the deviation from flatness [99].

8.6.4.1 Numerical Schemes

Osher et al. [30] have proposed an up-wind method for solving equations of the form = ■ v, of which = |V\$| ^iei(x), from Eq. (8.47), is an example. The up-wind scheme utilizes one-sided derivatives in the computation of |V where the direction of the derivative depends, point-by-point, on the sign of the speed term i ei(x). With strictly regulated time steps, this scheme avoids overshooting (ringing) and instability.

Under normal circumstances, the curvature term, which is a directional diffusion, does not suffer from overshooting; it can be computed directly from first-and second-order derivatives of \$ using central difference schemes. However, we have found that central differences do introduce instabilities when computing flows that rely on quantities other than the mean curvature. Therefore, we use the method of differences of normals [101,102] in lieu of central differences. The strategy is to compute normalized gradients at staggered grid points and take the difference of these staggered normals to get centrally located approximations to N (as in Fig. 8.20). The normal projection operator n <g> n is computed with gradient estimates from central differences. The resulting curvatures are

Figure 8.20: The shape matrix B is computed by using the differences of staggered normals.

treated as speed terms (motion in the normal direction), and the associated gradient magnitude is computed using the up-wind scheme.

### 8.6.4.2 Sparse-Field Method

The computational burden associated with solving the 3D, second-order, nonlinear level set PDE is significant. For this reason several papers [34,35] have proposed narrow-band methods, which compute solutions only for a relatively small set of pixels in the vicinity of k level set. The authors [36] have proposed a sparse-field algorithm, which uses an approximation to the distance transform and makes it feasible to recompute the neighborhood of the level set model at each time step. It computes updates on a band of grid points, called the active set, that is one point wide. Several layers around this active set are updated in such a way as to maintain a neighborhood in order to calculate derivatives. The position of the surface model is determined by the set of active points and their values.

The tomographic surface reconstruction problem entails an additional computational burden, because the measured data must be compared to the projected model at each iteration. Specifically, computing pj can be a major bottleneck. Computing this term requires recomputing the sinogram of the surface/object model as it moves. In the worst case, we would reproject the entire model every iteration.

To address this computational concern, we have developed the method of incremental projection updates (IPU). Rather than fully recompute p at every iteration, we maintain a current running version of p and update it to reflect the changes in the model as it deforms. Changes in the model are computed only on a small set of grid points in the volume, and therefore the update time is proportional to the area of the surface, rather than the size of the volume it encloses.

The IPU strategy works with the the sparse-field algorithm as follows. At each iteration, the sparse-field algorithm updates only the active layer (one voxel wide) and modifies the set of active grid points as the surface moves. The incremental projection update strategy takes advantage of this to selectively update

□ Area contributing to s

■ Area contributing to s

■ Area contributing to s

□ Area contributing to s

■ Area contributing to s

■ Area contributing to s

Figure 8.21: A weighting coefficient for each voxel determines the portions of the discrete sinogram influenced by incremental changes to a grid point.

the model projection to reflect those changes. At each iteration, the amount of change in an active point's value determines the motion of that particular surface point as well as the percentage of the surrounding voxel that is either inside or outside of the surface. By the linearity of projection, we can map these changes in the object shape, computed at grid points along the surface boundary, back into the sinogram space and thereby incrementally update the sinogram. Note that each 3D grid point has a weighting coefficient (these are precomputed and fixed), which is determined by its geometric mapping of the surrounding voxel back into the sinogram, as in Fig. 8.21. In this way the IPU method maintains subvoxel accuracy at a relatively low computational cost.

### 8.6.4.4 Initialization

The deformable model fitting approach requires an initial model, i.e. \$ (x, t = 0). This initial model should be obtained using the "best" information available prior to the surface fitting. In some cases this will mean thresholding a grayscale reconstruction, such as FBP, knowing that it has artifacts. In practice the initial surface estimate is impacted by the reconstruction method and the choice of threshold, and because we perform a local minimization, these choices can affect the final result. Fortunately, the proposed formulation is moderately robust with respect to the initial model, and our results show that the method works well under a range of reasonable initialization strategies.

Specimen

Contrast i-

Specimen

Contrast i-

Figure 8.22: (a) Transmission electron microscopy is used to image very small specimens that have been set apart from the substrate by a contrast agent. (b) TEM imaging technology provides projections over a limited set of angles.

Figure 8.22: (a) Transmission electron microscopy is used to image very small specimens that have been set apart from the substrate by a contrast agent. (b) TEM imaging technology provides projections over a limited set of angles.