Abstract

Existing physical artifacts including sculpture, mechanical parts, and anatomical structures are commonly acquired by modern surface and volumetric scanning technologies for archival, visualization, and diagnostic purposes. While the native representations for such data are largely sufficient for visualization purposes, more advanced field simulation currently requires extensive manual conversions into simplified surface and volume meshes compatible with the traditional finite element analysis pipeline. These conversions are tedious, error-prone, and require expertise in the mesh construction process. We propose to automate field simulation on acquired artifacts by bypassing the difficult geometric and topological meshing problems using a meshfree paradigm based on approximate distance fields computed from the native acquired data. To achieve complete automation, significant advances are required in the efficient computation of approximate distance from acquired data, surface and volume integration over the geometric domain, and visualization of the computational results on the acquired geometry. We are exploring this analysis process and investigating, developing, and implementing the necessary data structures and algorithms. The successful outcome of the proposed research will be an automated process for performing rapid field analysis on physical artifacts in situ.

An example of the Scan and Solve process in 2D is shown above. The process starts with a scan of the specimen. The resulting image is then processed to produce a binary image. The Euclidean distance transform is then applied to compute a distance map containing the distance to the boundary at every pixel. Samples of distance are taken from this map and a grid of bi-cubic B-splines is fit to these samples. A normalized trimming operation is performed to produce the distance field (shown partially) for the fixed boundary condition. Integration over the interior and boundary of the domain is performed to produce a linear system. The linear system is solved and the principal stress s1 is derived from the computed displacement field.

Distance-Based Meshfree Analysis

The proposed approach for automation uses meshfree technology described in detail elsewhere, but for our purposes it is sufficient to have a spatial function, ω, that is zero with non-vanishing gradient where the boundary conditions are specified. Inside the problem domain, ω must also be smooth, and non-zero. Central to the meshfree analysis process is the concept of a solution structure; a function that combines the specified boundary conditions with sufficient degrees of freedom to approximate the solution to the governing differential equation. For example, a problem with convective boundary conditions uses a solution structure that includes basis functions, ω, and derivatives of the basis functions in the direction of the gradient of ω. Incorporating a Ritz formulation, a linear system is produced with the known matrix and vector components assigned using volume and surface integrals computed over the interior and boundary of the domain. The solution of the system yields the vector of coefficients for the basis functions. These coefficients are then substituted back into the basis functions in the solution structure to yield the spatial distribution of the field.

ω plays a key role in the meshfree analysis process; appearing in the solution structure, and acting as a point membership classification function to facilitate computation of the necessary integrals. To serve these purposes, ω must take on a value of zero with non-vanishing gradient at the boundary while being non-zero everywhere within the domain. Further, because the solution structure includes products of ω with the basis functions, as well as derivatives of basis functions in the direction of the gradient of ω, it must be sufficiently differentiable that the basis functions are able to accurately model the distribution of the physical field without artifacts due to non-differentiability in ω. In addition, the solution structure and the integrals requires extensive evaluation of ω and its derivatives. Our representation for ω must therefore maximize efficiency. Euclidean distance to the boundary (hereafter, Euclidean distance) satisfies these requirements everywhere, except at the medial surfaces where it is not differentiable.

 

a. The zero set for the exact distance field computed using the Euclidean Distance Transform.	b. The zero set for an approximation represented on a 50x50x200 grid of tri-linear B-splines.	c. The zero set for an approximation represented on a 30x30x150 grid of tri-cubic B-splines.	d. The zero set for an approximation represented on a 10x10x50 grid of tri-cubic B-splines.
Figures show approximations with varying accuracies achieved using different grids of B-splines to approximate the Euclidean distance field. Each figure shows the zero set of the approximation extracted using a polygonization algorithm. The source CT data is the left femur from the Visible Human Male.

Computational Support

Distance Computation For voxel input, Euclidean distance can be rapidly computed by applying a Euclidean distance transform (EDT) to binary voxel data from segmenting the greyscale input. The EDT algorithm propagates distance values from boundary voxels to all other voxels in the domain. The resulting distance map holds at each voxel, the distance to the nearest boundary voxel. Further, distance values at voxels outside the domain are assigned a negative value to ensure the basis functions fit to the samples will pass through zero at the boundary and to encourage approximation of a gradient with magnitude 1. For other geometric input, spatial organizational schemes, or even the graphics processor, might be employed to minimize the computational overhead.

Numerical Integration The formation of the linear system requires numerical integration over the interior and boundary of the domain. Computation of both types of integrals can be abstracted to the parameterization of space to allocate quadrature points and weights over unparameterized geometry while guaranteeing integral results with a desired degree of accuracy. Further, these integrals must be computed over the individual supports of the basis functions being used to represent the solution field. Development of the volumetric integration algorithms is ongoing, but testing has shown that w may be used as a point membership classification function to support the allocation of quadrature points. Tensor-product Gaussian quadrature rules can be applied for supports fully interior to the domain, and subdivision and appropriate change of coordinate systems can be used for for supports intersecting the boundary. The computation of surface integrals is straightforward, given a representation of the geometric boundary in the form of a mesh or parametric surface. However, with geometric input in the form of segmented CT image slices or as ω, no parameterized representation of the boundary is available. Instead, methods must be developed to extract the zero set of ω by reconstruction, or to map points and weights from the points near the boundary (e.g. centroids of boundary voxels) to the zero set of ω.

Visualization Following solution of the linear system, and substitution of the coefficients into the solution structure, the results must be presented to the user for display. Like integration, the display process requires the parameterization of the geometric boundary to generate points where the solution field can be sampled, mapped to color, and projected for display. The display algorithms must be sufficiently fast to enable real-time interaction with the user.

[1] M. Freytag, V. Shapiro, and I. Tsukanov, Field Modeling with Sampled Distances, Computer Aided Design, Volume 38, Issue 2, February 2006, Pages 87-100.