U.S. patent application number 13/228517 was filed with the patent office on 2012-03-15 for methods and systems for producing an implant.
This patent application is currently assigned to CASE WESTERN RESERVE UNIVERSITY. Invention is credited to H. DAVID DEAN, ALEXANDROS T. MOULLAS, ROBERT A. RATCHESON, KRISHNAMOORTHY SUBRAMANYAN.
Application Number | 20120063655 13/228517 |
Document ID | / |
Family ID | 27386671 |
Filed Date | 2012-03-15 |
United States Patent
Application |
20120063655 |
Kind Code |
A1 |
DEAN; H. DAVID ; et
al. |
March 15, 2012 |
METHODS AND SYSTEMS FOR PRODUCING AN IMPLANT
Abstract
A computer implemented method for determining the 3-dimensional
shape of an implant to be implanted into a subject includes
obtaining a computer readable image including a defective portion
and a non-defective portion of tissue in the subject, superimposing
on the image a shape to span the defective portion, and determining
the 3-dimensional shape of the implant based on the shape that
spans the defective portion.
Inventors: |
DEAN; H. DAVID; (SHAKER
HEIGHTS, OH) ; SUBRAMANYAN; KRISHNAMOORTHY;
(PALATINE, IL) ; MOULLAS; ALEXANDROS T.;
(THESSALONIKI, GR) ; RATCHESON; ROBERT A.; (SHAKER
HEIGHTS, OH) |
Assignee: |
CASE WESTERN RESERVE
UNIVERSITY
CLEVELAND
OH
|
Family ID: |
27386671 |
Appl. No.: |
13/228517 |
Filed: |
September 9, 2011 |
Related U.S. Patent Documents
|
|
|
|
|
|
Application
Number |
Filing Date |
Patent Number |
|
|
10089467 |
Mar 27, 2002 |
|
|
|
PCT/US2000/022053 |
Aug 11, 2000 |
|
|
|
13228517 |
|
|
|
|
60148393 |
Aug 11, 1999 |
|
|
|
60148275 |
Aug 11, 1999 |
|
|
|
60148277 |
Aug 11, 1999 |
|
|
|
Current U.S.
Class: |
382/128 |
Current CPC
Class: |
B33Y 50/00 20141201;
G06T 17/00 20130101; G06K 9/34 20130101; G16H 50/50 20180101; G06T
7/0014 20130101; G06K 2209/05 20130101; G16H 30/40 20180101; B33Y
80/00 20141201; G16H 20/40 20180101; A61F 2002/2835 20130101; G06T
2207/30008 20130101; A61F 2/28 20130101; A61F 2/30942 20130101;
A61B 6/505 20130101; G06K 9/6251 20130101 |
Class at
Publication: |
382/128 |
International
Class: |
G06K 9/00 20060101
G06K009/00 |
Claims
1. A computer implemented method for determining the 3-dimensional
shape of an implant to be implanted into a subject, the method
comprising: obtaining a computer readable image including a
defective portion and a non-defective portion of tissue in the
subject; superimposing on the image a shape to span the defective
portion; and determining the 3-dimensional shape of the implant
based on the shape that spans the defective portion.
2. The method as set forth in claim 1, wherein the shape that spans
the defective portion is drawn on the image as the image is
displayed by a computer.
3. The method as set forth in claim 1, wherein the shape that spans
the defective portion is a right-left mirror image of a portion of
the subject's anatomy that includes the defective portion.
4. The method as set forth in claim 3, wherein the determining the
3-dimensional shape of the implant is determined as a function of
respective shapes of the defective portion and the right-left
mirror image.
5. The method as set forth in claim 1, wherein the shape that spans
the defective portion in the subject is a template representing one
of corresponding normal tissue in another subject and an average
shape of corresponding normal tissue.
6. The method as set forth in claim 5, wherein the determining the
3-dimensional shape of the implant is determined as a function of
respective shapes of the defective portion and the template.
7. The method as set forth in claim 1, wherein the superimposing on
the image the shape to span the defective portion includes: mapping
highly curved portions of surfaces of at least one of: the subject
and another individual, the subject and a first average, and two or
more individuals other than the subject to produce a second
average; and establishing correspondence between the mapped highly
curved portions.
8. The method as set forth in claim 1, wherein the computer
readable image consists of slices of the defective portion and the
non-defective portion.
9. The method as set forth in claim 1, wherein the computer
readable image consists of scan lines of the defective portion and
the non-defective portion.
10. The method as set forth in claim 1, wherein the computer
readable image consists of voxels of the defective portion and the
non-defective portion.
11. The method as set forth in claim 1, further including:
determining a position for optimal adjacency between the implant
and the non-defective portion of the subject.
12. A system for determining a 3-dimensional shape of an implant to
be implanted into a subject, the system comprising: a
computer-readable medium having data representing an image
including a defective portion and a non-defective portion of tissue
in the subject; a graphical user interface operable to superimpose
on the image a shape to span the defective portion; and a logic
configured to determine the 3-dimensional shape of the implant
based on the shape that spans the defective portion.
13. The system as set forth in claim 12, wherein the graphical user
interface is further operable for a user to draw the shape that
spans the defective portion onto the image.
14. The system as set forth in claim 12, wherein the shape that
spans the defective portion is a right-left mirror image of a
portion of the subject's anatomy that includes the defective
portion.
15. The system as set forth in claim 14, wherein the determining
the 3-dimensional shape of the implant is determined as a function
of respective shapes of the defective portion and the right-left
mirror image.
16. The system as set forth in claim 12, wherein the shape that
spans the defective portion in the subject is a template
representing one of corresponding normal tissue in another subject
and an average shape of corresponding normal tissue.
17. The system as set forth in claim 16, wherein the determining
the 3-dimensional shape of the implant is determined as a function
of respective shapes of the defective portion and the template.
18. The system as set forth in claim 12, further comprising: a
second logic configured to map highly curved portions of surfaces
of at least one of: the subject and another individual, the subject
and a first average, and two or more individuals other than the
subject to produce a second average; and establish correspondence
between the mapped highly curved portions.
19. The system as set forth in claim 12, wherein the image consists
of slices of the defective portion and the non-defective
portion.
20. The system as set forth in claim 12, wherein the image consists
of scan lines of the defective portion and the non-defective
portion.
21. The system as set forth in claim 12, wherein the image consists
of voxels of the defective portion and the non-defective
portion.
22. The system as set forth in claim 12, wherein the logic is
further configured to determine optimal adjacency between the
implant and the non-defective portion of the subject.
Description
[0001] This application is a continuation of U.S. application Ser.
No. 10/089,467, filed Mar. 27, 2002, which is a U.S. National Stage
Entry of International Application PCT/US2000/22053, filed Aug. 11,
2000, which claims the benefit of U.S. Provisional Application Nos.
60/148,275, 60/148,277, and 60/148,393, which were all filed Aug.
11, 1999. All prior applications are incorporated by reference.
BACKGROUND OF THE INVENTION
[0002] The present invention relates to fabricating implants to
replace bony structures. More specifically, the invention relates
to a system and methodology for fabricating a "drop in" replacement
for a particular segment of missing bony structure, in which the
implant fits precisely within the contours of the missing segment
and thus minimizes complications during the surgical procedure to
install the implant, and subsequently during the healing process.
It will be appreciated, however, that the invention is also
amenable to other like applications.
[0003] Various systems and methods of fabricating prosthetic
implants are known in the prior art. Examples of such prior systems
and methods include U.S. Pat. Nos. 4,436,684; 5,274,565; 5,357,429;
5,554,190; 5,741,215; and 5,768,134. Each of these patents,
however, suffer from many disadvantages that have collectively
limited the usefulness of their methods and implants to the
relevant field.
Neural Network Segmentation
[0004] Although clinical visualization of three dimensional organ
surfaces embedded in an imaged volume is now common, segmentation
techniques are often ineffective and time consuming Current three
dimensional (3D) surface segmentation is often dependent on the
contrast and resolution of two dimensional (2D) anatomical
boundaries identified by an operator proceeding slice by slice
through the image volume. Existing automated 3D surface and volume
segmentation algorithms, such as region growing, are most commonly
based on voxel intensity, nearest neighbor similarity, and edge
strength. An active contours energy minimization approach attempts
to automate 2D outline detection. An extension to 3D suggests
regions to be segmented by internally enveloping them with three
dimensional bubbles. Another method based explicitly on high level
volumetric shape features utilizes prior information obtained from
previously imaged anatomical objects. This approach is based on
prior morphological knowledge and a highly integrated set of model
based grouping operations. The use of prior information allows the
last method to be more adaptive to an expected surface when it
encounters noise due to artifact or non-descript structures such as
in thin areas. However, detection of anatomical surfaces as complex
as the skull or soft tissue face with prior information alone is
computationally expensive. It is likely manual work will be
necessary where the signal is too low, especially in regions with
high artifact noise or where artifact-interrupted surface
morphology may obfuscate a match with database information. The
relationship of artifacts to surface morphology is dependent on
both surface proximity to the artifact-generating structure and
patient pose. To overcome these limitations, a third approach,
combines semi-automatic image segmentation procedures with
intelligent scissors. A user steers the image segmentation via live
wire and live lane tools through the image volume.
[0005] Current slice-based segmentation methods tend to be time
consuming. Imprecision correlates with the numbers of decisions the
operator is forced to make. These are concentrated in areas of poor
contrast, artifact, or insufficient geometry for the user to
recognize local aspects of the patient's anatomy (i.e., the
sufficiency of anatomical cues are heavily dependent on structural
size and degree of curvature as well as how patient pose affects
the representation in tomographic images). Difficulty in making
segmentation decisions based on these three (3) types of
information may be compounded by a user environment where an image
volume is presented as a series of tomographic or projection images
that lose resolution when reformatted. In all cases these two
dimensional views require anticipation of how the segmented
outlines will combine to form the 3D surface that the user wishes
to render.
[0006] Deformable curve fitting for edge detection seems less
adaptable than manual segmentation for determining the location of
anatomical edges. Curve fitting runs into trouble when edges vary
in signal strength, requiring multiple initialization and guidance
steps to obtain best-fit. An operator familiar with the salient
anatomy quickly adapts to provide continuity between areas of
material inhomogeneities (i.e., varying pixel intensity) that may
signal a break to an edge detection algorithm. It is possible that
these edges can be spanned by the use of geometrically ideal shapes
or those based on prior information. However, these simple models
are less adaptable, or less decisive, than manual segmentation.
[0007] Neural networks, applied to volume image segmentation, are
capable of modeling structural components in the image through the
processing of the gray scale and spatial information recorded in
volume images. They take advantage of the multidimensional nature
of medical images. Neural network classification can be either
supervised or unsupervised. To date, supervised neural networks
have been preferred, despite their dependence on user interaction
to train or correct their results. Unsupervised neural networks do
not rely on pre-labeled samples. Rather, they require careful
selection of input events that, when clustered, are useful to the
user.
Fitting a Deformable Template to a Segmented Surface
[0008] There is debate as to whether comparison of three
dimensional (3D) anatomical shapes should be done on the basis of
whole objects or elementally. One conventional method for 3D
anatomical surface parsing begins with identification of a
specialized class of space curves, which are referred to as "ridge
curves."
[0009] A consensus ridge curve definition would cite a maximum
first principal curvature along sharp edges as allowing their
representation as a space curve. The curvature found along ridge
curves is regular enough that a class of curvature extremal
landmarks can be readily detected. These are referred to in the
literature as Type II landmarks, which may be extracted manually or
automatically. These landmarks may form the initial seed points on
which to register a wire frame (i.e., ridge curve and "geodesic"
lines) plate. "Geodesic" space curves, a shortest path traversed on
a graphical manifold, are useful as a means to link Type II
landmarks. The resulting fitted ridge curve and "geodesic"
wireframe are used as a means to tile out and parameterize 3D
surface images. This toolkit allows parameterized (labeled) surface
extractions, multi-surface registration, and averaging for
cephalometric comparisons of the skull, soft tissue face, and
cerebral ventricles.
[0010] In regards to parameterized 3D surface extraction, the
toolkit facilitates an initially manual superimposition of a ridge
curve-based deformable template to the perspectively rendered
surface of interest at the Type II landmarks. However, there is
only an approximate tie between the surface tile points bounded by
ridge and "geodesic" space curves, and the surface voxels
originally identified (segmented) in the raw CT or MR scan volume
image data. The "geodesic" criteria are less appropriate and these
lines have been renamed "tiling curves." Additionally, the entire
surface is encoded as set of B-spline space curves, in order to
facilitate whole surface averaging algorithms, (i.e., to average
multiple surfaces space curve by space curve). The toolkit links
biological labeling (homology mapping) and geometric surface
subdivision, resulting in a homology mapped (tile name)
parameterized (u,v space of tile) surface.
[0011] The primary basis for positioning multiple biological
surfaces for conventional averaging is a specialized class of space
curves referred to as "ridge curves". One definition of ridge
curves describes them as biologically significant space curves
found along sharp edges, corresponding to first principal curvature
maxima. These sharp edges are the most bent part of the surface.
Type II landmarks, identified as recurrent curvature maxima
(bends), are found along ridge curve space curves. They span these
landmarks with what are called "geodesics" in order to create three
or four sided tiles ("tiling curves").
[0012] It has been proposed to add tile parameters (i.e., the u,v
coordinates and tile name) to the xyz point coordinates found in 3D
skull surface images to produce a feature-preserving average or
morphometric comparison of homologous landmark. A geometric
subdivision (ridge curve and geodesic tile boundaries) of the
surface links biological labels (homology mapping) to tile grid (u,
v) parameterization of the extracted points (x, y, z). The toolkit
begins superimposition of the ridge curve and geodesic (tiling
curve) deformable wireframe template with the identification of
Type II landmarks. Using a thin plate spline, the ridge curve and
geodesic template wireframe is warped to the surface's Type II
landmarks. Superimposition of the interceding ridge curve and
geodesic wire segments to best position proceeds automatically.
Averaging of Multiple Tiled Out Specimens
[0013] When too much of the patient is missing from the image, it
is useful to match the average to the patient to obtain a
reasonable shape to fill the defect. Unlike images of individuals,
averages are usually relatively smooth and symmetric (right to
left). However, estimating the shape of a defect of large size with
an average of too few subjects is also not helpful. In this case,
it may be better to use the image of a single relative.
[0014] Homology-based 3D surface averaging requires appropriate
superimposition of sample member surface points within a shared
reference frame. The toolkit obtains this position in a series of
global warping and local unwarping maneuvers. First an average of
all sample specimen's Type II landmarks is obtained. Each specimen
is then globally warped to these landmarks and then locally
unwarped. A rough average of the ridge and geodesic curve sample is
produced. Next, because the space curve sample members are fit to
the resulting rough average space curve, perpendicular bisectors
may be positioned equidistantly along the rough average. An average
of the intercepted points appropriately samples the shape of each
space curve on the entire surface. The u,v parameters determine the
number of points along either side of the four (4) sided tiles and,
therefore, contain u x v internal space curves, which are averaged
in the same manner as the boundary curves. This method insures that
these sampling planes (perpendicular bisectors) are never parallel
to the sampled data.
[0015] The present invention provides a new and improved apparatus
and method which overcomes the above-referenced problems and
others.
SUMMARY OF THE INVENTION
[0016] A method for determining a shape of a medical device to be
implanted into a subject produces an image including a defective
portion and a non-defective portion of a surface of a tissue of
interest included in the subject. The tissue of interest is
segmented within the image. A template, representing a normative
shape of an external anatomical surface of the tissue of interest,
is superimposed to span the defective portion. An external shape of
an implant, is determined as a function of respective shapes of the
defective portion as seen in the template, for repairing the
defective portion.
[0017] In accordance with one aspect of the invention, a volumetric
image of the tissue of interest is produced. In accordance with a
more limited aspect of the invention, one of a CT and MR image is
produced.
[0018] In accordance with another aspect of the invention, a
position for seating the implant into the defective portion is
determined. The implant is tapered as a function of a seating
strategy and a space available for the implant.
[0019] In accordance with another aspect of the invention, the
template is superimposed to the surface of the tissue of interest
via warping.
[0020] In accordance with another aspect of the invention, the
template is warped to an external surface of the non-defective
portion of the tissue of interest.
[0021] In accordance with another aspect of the invention, a ridge
and tiling curve homology map is produced of another anatomical
surface. The average normative shape is produced from control
surface image representations of the anatomical surface of
interest.
[0022] In accordance with a more limited aspect of the invention,
the normative shape is determined from a mirror image of the tissue
of interest.
[0023] In accordance with another aspect of the invention, the
normative shape is modified as a function of a shape of the
non-defective portion of the anatomical surface of interest.
[0024] A system determines a shape of a medical device to be
implanted into a subject. An imaging device produces an image
including a defective portion and a non-defective portion of a
tissue of interest included in the subject. A segmenting means
segments the tissue of interest. A template spans the defective
portion. The template represents a normative shape of an external
anatomical surface of the tissue of interest. A determining means
determines an external shape of an implant, as a function of
respective shapes of the defective portion and the template, for
repairing the defective portion.
[0025] In accordance with one aspect of the invention, the image is
a volumetric image of the anatomical surface of interest.
[0026] In accordance with a more limited aspect of the invention,
the volumetric image is one of a CT and MR image.
[0027] In accordance with another aspect of the invention, the
template specifies a seated edge, which is fixed to the subject,
and a surface shape, which is not affixed to the subject.
[0028] In accordance with another aspect of the invention, the
template is warped to a deformed bone.
[0029] In accordance with another aspect of the invention, the
template is warped to an external surface of the non-defective
portion of the tissue of interest.
[0030] In accordance with another aspect of the invention, the
normative shape of the template is determined from an additional
tissue representative of the hard tissue of interest.
[0031] In accordance with another aspect of the invention, the
normative shape of the template is determined as a function of an
average shape of the tissue of interest.
[0032] A method for repairing a defect in a tissue of interest
included in a subject produces a volumetric image showing a
defective portion, which includes the defect, and a non-defective
portion of the hard tissue of interest within the subject. The
tissue of interest is segmented from the image. A template, having
an average shape of the tissue of interest, is warped over the
defective and non-defective portions. A shape of the implant is
determined, as function of respective shapes of the defective
portion and the template. The implant is inserted into the
defective portion for repairing the defect.
[0033] In accordance with one aspect of the invention, the implant
is cast.
[0034] In accordance with another aspect of the invention, a mirror
of another tissue of interest, representing a bilateral mirror
image of the tissue of interest, is homology mapped. The normative
shape is determined from the mirrored other section of tissue.
[0035] One advantage of the present invention is that it presents
an unsupervised neural network algorithm (i.e., self organizing
feature map, or SOFM) for automated and accurate 3D surface
segmentation, the latter being a necessity for computer generation
of cranial implants, diagnostic cephalometrics and image guided
surgery.
[0036] Another advantage of the present invention is that user
interaction in the SOFM algorithm is limited to determining which
segmented surfaces represent the anatomy desired by the user.
[0037] Another advantage of the present invention is that it
rapidly assembles automatically segmented, and anatomically valid
(i.e., no correction required), surfaces into a surface of interest
for rendering of 3D surfaces found in volume images.
[0038] Another advantage of the present invention is that it saves
time over manually identifying each pixel or currently available
semi-automatic methods.
[0039] Another advantage of the present invention is that it
locates continuous surfaces, even in areas of thin structure or
imaging artifact (e.g., metallic artifact in x-ray or CT).
[0040] Another advantage of the present invention is that it
combines spatial and intensity information to track 3D surfaces for
segmentation.
[0041] Another advantage of the present invention is that it
provides a means for anatomically valid segmentation (because 3D
surfaces are tracked oblique to the scan plane, irrespective of
patient pose (orientation of scan planes are determined relative
the patient anatomy) in the scanner).
[0042] Another advantage of the present invention is that it
provides a warping method that gives precise, accurate and
efficient ridge curve-based template registration for extraction of
a parameterized surface.
[0043] Another advantage of the present invention is that a fitted
template allows matching of mirrored (filling defect with right or
left opposing side data) or average (normative) surfaces, and
production of averages. The template establishes the homology for
these three types of surface matching (i.e., mirroring, one-to-one
fit, averaging).
[0044] Another advantage of the present invention is that it
provides a reduced time for the warping method that is largely due
to an automated (simulated annealing) search for ridge and tiling
curves in the volume image.
[0045] Another advantage of the present invention is that it
provides faster and more reproducible and, consequently, more
reliable automated ridge and tiling curve fitting, especially for
internal tile points which were previously generated from an error
map between an ideal surface and the true surface.
[0046] Another advantage of the present invention is that it
provides a new basis for determining tiling curve and internal tile
point location by establishing more reliable homology mapping of
these features.
[0047] Another advantage of the present invention is that it
provides a template fitting and surface averaging method that will
work well on the entire external surface of the body, as well as
many internal organs.
[0048] Another advantage of the present invention is that it
encodes the entire surface as series of B-spline space curves,
including the internal tile points.
[0049] Another advantage of the present invention is that all
surface operations (e.g., registration or averaging) may be
accomplished with series of homologous space curves, rather than
whole surfaces, thereby allowing for computationally less expensive
serial or parallel processing of space curve sets.
[0050] Another advantage of the present invention is that it
improves the homology assignment of both tiling curves and internal
tile points through new surface extraction techniques.
[0051] Another advantage of the present invention is that it uses
average 3D surface images to model surfaces in patient images for
rapid prototyping of prosthetic implants.
[0052] Still further advantages of the present invention will
become apparent to those of ordinary skill in the art upon reading
and understanding the following detailed description of the
preferred embodiments.
BRIEF DESCRIPTION OF THE DRAWINGS
[0053] The invention may take form in various components and
arrangements of components, and in various steps and arrangements
of steps. The drawings are only for purposes of illustrating a
preferred embodiment and are not to be construed as limiting the
invention.
[0054] FIG. 1A illustrates a bone of interest according to the
present invention;
[0055] FIG. 2 illustrates a flowchart according to the present
invention;
[0056] FIG. 3 illustrates a local negative feedback which excites
neuron drop off allowing coordinated input of simulated neuronal
field to a neural network;
[0057] FIG. 4 illustrates event inputs to a Network;
[0058] FIG. 5 illustrates a current implementation of a SOFM neural
network configuration shown as flow of event space image element
input to processor layer of neural network;
[0059] FIG. 6 illustrates Gaussian and Zero crossing curves;
[0060] FIG. 7 illustrates processor network configurations
according to the present invention;
[0061] FIGS. 8 and 9 illustrate flowcharts of the segmentation
process;
[0062] FIG. 10 illustrates multiple cluster back-projection from
map nodes;
[0063] FIG. 11 illustrates a SOFM-based image segmentation and
surface generation interface;
[0064] FIGS. 12A-12E illustrate a feature cluster map, back
projected pixels and triangulated surface rendering; FIGS. 12A and
12B display the resulting triangulated surface rendered soft tissue
face and skull. FIGS. 12C and 12D present the voxel back projection
from the selected feature cluster nodes; FIG. 12E presents the
feature cluster map;
[0065] FIGS. 13A and 13B illustrates a self organizing feature map
formation process in which an initial random processor location
evolves into the intended grid-like structure as ordering is
obtained by the neural network;
[0066] FIG. 14 illustrates both sessions of manual and SOFM
segmented surfaces for skull surfaces;
[0067] FIG. 15 illustrates both sessions of manual and SOFM
segmented surfaces for soft tissue face surfaces;
[0068] FIGS. 16A and 16B illustrate plots of accumulated
differences across all five data sets per slice the soft tissue
face surface and the skull surface, respectively;
[0069] FIG. 17 illustrates orthogonal projections of 2D images;
[0070] FIG. 18 illustrates a new ridge curve-based deformable
template wireframes superimposed on a soft tissue face graphical
manifold, with Type II landmarks, ridge curve arcs, and tiling
curves indicated;
[0071] FIGS. 19A and 19B illustrate Type II landmarks of the soft
tissue face and boney skull templates;
[0072] FIGS. 20A and 20B illustrate a CT slice highlighting
segmented surface voxels of interest which are indexed to vertices
(points) in an adjacent rendered manifold surface;
[0073] FIGS. 21A and 21B illustrate a candidate wireframe prior to
the warp at lower left and the result of the thin plate spline warp
of the wireframe onto the Type II landmarks, respectively;
[0074] FIG. 22A illustrates a space curve C from a sweeping plane
positioned at each control point, with tangent vector {right arrow
over (t)}, along a normal section of segmented voxel surface S;
[0075] FIG. 22B illustrates a surface normal curvature k.sub.n in
direction {right arrow over (t)}, where k.sub.n attains principal
curvature values k.sub.1 and k.sub.2;
[0076] FIGS. 23A and 23B illustrate a plane constructed from the
two endpoint landmarks and an object's center of mass;
[0077] FIG. 24 illustrates a trace along the candidate template's
control point to the image surface;
[0078] FIG. 25A illustrates a cross section on the top of the skull
and tiling curve control point search planes intersecting that
image surface;
[0079] FIG. 25B illustrates a series of search plane
iterations;
[0080] FIG. 26A illustrates rotating search planes;
[0081] FIG. 26B illustrates search planes for one ridge curve
arrayed along a warped ridge curve-based deformable template space
curve;
[0082] FIG. 26C illustrates a final ridge curve registration found
by the simulated annealing algorithm;
[0083] FIG. 26D illustrates a candidate tiling curve and its
associated search planes;
[0084] FIG. 26E illustrates the fitted control points;
[0085] FIG. 27 illustrates an entire final ridge and tiling curve
wireframe superimposed on the image surface;
[0086] FIG. 28 illustrates a flowchart for minimizing the cost
function;
[0087] FIG. 29 illustrates a temperature cooling schedule from an
initial temperature 100 to a final temperature, -10;
[0088] FIGS. 30A, 30B, and 30C illustrates an overview of a ridge
curve-based deformable template registration process in SASE;
[0089] FIG. 31 illustrates three orthogonal projection views
showing direct differences between two surface extractions;
[0090] FIGS. 32A and 32B illustrate the soft tissue face surfaces
and boney skulls, respectively, obtained via different methods
versus the toolkit on the right side;
[0091] FIGS. 33A and 33B illustrate color coded difference images
of the skull and soft tissue face surface extraction comparisons
with the segmented voxel surface and color coded difference images
of SASE skull and soft tissue face surface extractions,
respectively;
[0092] FIGS. 34A and 34B illustrate a tile point density;
[0093] FIG. 35 illustrates the toolkit and SASE surface
extractions;
[0094] FIG. 36A illustrates an alignment of u and v internal tile
curves in a 4 sided tile;
[0095] FIG. 36B illustrates an alignment of u and v internal tile
curves in a 3 sided tile;
[0096] FIG. 37A illustrates a soft tissue face ridge curve-based
deformable template definition;
[0097] FIG. 37B illustrates a boney skull deformable template
definition;
[0098] FIG. 38A illustrates a translation of two shapes into a
shared two dimensional common coordinate system;
[0099] FIG. 38B illustrates scaling adjustments between the two
shapes;
[0100] FIG. 38C illustrates the rotation, in order to obtain the
Procrustes distance (square root of sum of squared) between
homologous point landmarks pairs: A . . . D, a . . . d;
[0101] FIG. 39 illustrates a space curve averaging procedure;
[0102] FIG. 40A illustrates a soft tissue face average surface
generation;
[0103] FIG. 40B illustrates a skull average surface generation;
[0104] FIG. 41A illustrates a soft tissue face surface average
comparison of various methods;
[0105] FIG. 41B illustrates a boney skull surface average
comparison of various methods;
[0106] FIG. 42A illustrates a surface average tile-by-tile color
coded error maps for skull and soft tissue face surface
extractions;
[0107] FIG. 42B illustrates an SSA surface average tile-by-tile
color coded error maps for skull and soft tissue face surface
extraction;
[0108] FIG. 43A illustrates a tile point density for soft tissue
face average;
[0109] FIG. 43B illustrates an SSA method tile point density;
[0110] FIG. 43C illustrates a boney skull surface average;
[0111] FIG. 43D illustrates an SSA method skull tile point density;
and
[0112] FIG. 44 illustrates an inter-method average surface
continuity.
DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS
[0113] With reference to FIGS. 1, 1A, and 2, a subject 10, or
patient, is imaged using an imaging device 12 in a step A.
Preferably, the image is a volumetric image obtained using a
computerized tomography ("CT") scanner according to a preferred
protocol. However, other images, obtained using other imaging
techniques (e.g., magnetic resonance) are also contemplated. Once
the data is captured from the CT scanner 12, it is stored in a
processor 14 (preferably, a transportable storage medium) in a step
B (although it could be transferred over a network). Then, the
processor 14 segments the data in a step C for extracting a region
of the image including a target tissue (e.g., liver, lung, or hard
tissue such as bone) 16 of interest within the subject; the
extracted region 16 includes two (2) portions (i.e., one portion
having a defect 18 and one portion 20 without the defect).
Preferably, the non-defective portion 20 substantially surrounds
the defective portion. An external surface of the extracted region
16 of the image is mapped in a step D. A decision is made in a step
E whether to register the mapped data of the non-defective portion
of the extracted region to a memory device in the processor 14,
which includes normative (average) data for various bones as a
function of various factors (e.g., race, sex, age, etc) of the
subject. If the data is transmitted to the memory device, the data
representing points on the mapped surface of the non-defective
portion of the extracted region are incorporated into averaged data
for the respective bone based on the various factors of the subject
in a step F.
[0114] After the image is segmented and mapped, a template 24 is
superimposed in a step G, to an external surface of a normative
shape of the target tissue 16 (i.e., without the defect 18).
Preferably the template is superimposed via warping. However, a
best fit process is also contemplated. The template specifies a
"seated" edge, which is fixed to the subject (patient), and the
remaining surface shape, which is not affixed to the subject.
[0115] It is to be understood that the normative shape represents a
desired or average shape of the target tissue. Preferably, the
normative shape is determined in a bilateral structure from a
substantially mirror image of the target tissue (e.g., a defect is
present in a left side of a skull and complete data is available
for a right (control) side of a skull). However, if too much of the
target tissue is missing on the affected side, or the defect
crosses the midline in a way that prevents establishing the
midline, appropriate normative average surfaces (which, as
described below, are already homology mapped) are used. In this
case, the normative average data is brought in and registered to
the homology mapped normative surface surrounding the defect. Once
the template is warped, points on the template are mapped, in a
step H, to points on the external surface of the normative shape of
the bone of interest 16.
[0116] A shape of an implant 26 to be "dropped in" the defect 18 is
determined, in a step I, as a function of a difference between the
mapped points on the external surface of the target tissue and the
external surface of the template. The implant 26 is seated (fit) in
the defect 18 of the target tissue 16 in a step J. A fit of the
implant (seating and shape relative to surrounding anatomy) is
analyzed in order to produce an image description of the best shape
for the implant on computer, not manually. That image description
is saved in industry standard format, such as stereolithography
(STL), and printed on a 3D rendering device (e.g.,
stereolithographic printer). The fit is then accomplished using a
tapering scheme based on the seating strategy and space available
for the implant.
[0117] If the implant is not produced in biocompatible material, it
is cast in an implantable material. It is to be understood that the
present methods work with inert casting materials (e.g.,
bioplastics such as Polymethymethacrylate) or future tissue
engineered materials.
[0118] A stereolithography manufacturing method is contemplated for
forming the either the implant master shape or the actual
implant.
I. Segmentation
[0119] As described below, the step C of segmenting the image is
preferably executed using self organizing feature maps ("SOFM"),
which are used for computationally encoding neural maps.
[0120] The SOFM incorporates a modified Hebbian rule. The resulting
mapping of data ordering resembles observed neurological patterns.
Since then, the SOFM has found widespread use in many practical
engineering applications. The SOFM derives primarily from
observations of two biologically equivalent computational
hypotheses defined as:
[0121] (1) The Hebbian Rule: "A synapse is strengthened or
stabilized when pre- and post synaptic activations are correlated."
Change in synaptic strength is defined as:
s.sub.t+1=s.sub.t+.alpha..eta. (1)
where s, is a connection strength at time t, .eta. is an output (IP
in FIG. 4) equal to s.sub.t.times..xi., .xi. is an input pattern,
and .alpha. is a factor, a user-selected value, that is made to
decrease for each iteration, to ensure it reaches zero. This, in
turn, ensures appropriate synaptic strength formation and system
stabilization.
[0122] (2) Local Feedback: Each neuron laterally influences the
neighboring neurons to either an excitatory or inhibitory state.
This results in a competition which enables the formation of a map.
FIG. 3 is an example of an excitation function. The SOFM based 3D
surface segmentation presented here extends the concepts of ordered
feature maps. In the neural network, each normalized input feature
vector is scaled or weighted by multiplication with a user-assigned
constant vector. This process begins with an event generation step,
which results in a pattern .xi.={.xi..sub.1, .xi..sub.2, . . . ,
.xi..sub.n} that is sent to a "resonator" from which potential
connection are made. The weights determine feature connection
strengths. The scaling results in a vector xi=c*F, with c the
feature connection strengths and F the input vector are shown in
FIG. 4. Processor weights adapt to input stimulation and in so
doing the neural network organizes an ordered map.
[0123] The unit i has connection weights, .mu..sub.i1, .mu..sub.i2,
. . . , .mu..sub.in, which are expressed as a column vector
mi=[.mu..sub.i1, .mu..sub.i2, . . . , .mu..sub.in].sup.T. The
discriminant function of a unit is defined as:
.eta. i = j = 1 n .mu. ij .xi. j = m i T .xi. ( 2 )
##EQU00001##
[0124] For unit k after selecting a winner from each iteration, the
weight vector is updated as:
m i ( t + 1 ) = m i ( t ) + .alpha. ( .xi. - m i ( t ) ) m i ( t )
+ .alpha. ( .xi. - m i ( t ) ) N ( 3 ) ##EQU00002##
where t is the discrete point at which events are sequentially sent
to the network. .parallel.o.parallel..sub.N is the normalization
process, ensuring the weight vectors are bounded. The network's
selectivity improves as more input vectors are presented. Network
formation requires no initial correlation or order between any of
the voxel input data elements. Hence, no initialization or operator
assistance is required during feature cluster formation. Iterative
steps lead array units to represent the topological order of
different combinations of input image features.
[0125] 1. Voxel Input Data Elements (VIDE)
[0126] One dimensional feature maps with voxel intensity input data
have been extended to include six (6) voxel input data elements
(i.e., patterns) in a two dimensional processor and feature map for
interactive surface segmentation. The image feature vector is
generated for each input voxel, its elements are: (1) A distance
from volume centroid to voxel location in three dimensional space;
(2) The histogram equalized or normalized voxel intensity; (3) The
Gaussian smoothed voxel intensity; 4) The difference of Gaussians;
5) The Laplacian zero crossings; and, 6) Invariant image
moments.
[0127] The calculation of VIDE parameters is explained below. In
general, intensity information is incorporated to provide voxel
differentiation. Spatial distance parameters provide positional
information, such as where a surface is located. Gaussian
parameters used in self organizing map-based segmentation provide
additional adaptive qualities to the 3D surface detection operation
by smoothing artifact noise and compensating for scan protocol
sequellae (resolution, pose, and beam strength). The "difference of
Gaussian" parameter extends smoothing across a detected surface by
introducing effects from local and distant voxels (i.e., the
maximum distance between voxels is that of the surface detected but
not elsewhere in the image volume). The "zero crossing" parameter
provides transitional information important for some kinds of edge
detection. The "moment invariant" parameters provide a positional
and scale invariant signature of a surface in each region of the
image volume; a computationally intensive parameter, moment
invariants are initially set to 10 weight (Table I), and raised
here only in small, user-defined, regions (e.g., the margins of
thin structures) for fine tuning.
[0128] This set of voxel input data element vectors have normalized
values of: spatial location via distance transform, voxel density
with respect to nearest neighbors via Gaussian and difference of
Gaussians, and boundary changes captured via Laplacian zero
crossings. The two Gaussian input elements prevent feature
inclusion of the spurious background noise commonly associated with
3D images. These parameters are computed during a preprocessing and
feature-extraction procedure (see FIG. 5). All VIDEs are normalized
to a unit range of 0 to 1.0. This bounds the neural network output
(i.e., second layer produced from VIDE input).
[0129] The voxel input data element vectors are defined as row
vectors F.sub.n=[E.sub.dn,V.sub.n,g,g.sub.d,Z.sub.c,m] with six (6)
elements, shown in Table I, where n varies from 1 to the total
number of input image voxels. The weighting factors are defined for
each element in order to prefer one element over the others. The
resultant VIDE vector .xi. is defined as:
.xi.=.kappa.F.sub.n=[.kappa..sub.1E.sub.dn,.kappa..sub.2V.sub.v,.kappa..-
sub.3g,.kappa..sub.4g.sub.d,.kappa..sub.5Z.sub.c,.kappa..sub.6m]
(4)
where .kappa..sub.i defines the percentage of preference assigned
for each VIDE. The differential preference of input VIDE results in
associations between neural network processors and the voxel
clusters. Table 1 presents the default preferences for each element
during experiments.
[0130] A. Input element 1: Normalized Euclidean Distance
(E.sub.dn)
[0131] The radial distance from the image volume centroid is
referred to as the moment of inertia. The image volume centroid is
(C.sub.x,C.sub.y,C.sub.z) as determined from:
C X = ( z = 0 d y = 0 h x = 0 w xV ( x , y , z ) ) / ( z = 0 d y =
0 h x = 0 w x ) , C y = ( z = 0 d y = 0 h x = 0 w yV ( x , y , z )
) / ( z = 0 d y = 0 h x = 0 w y ) , C z = ( z = 0 d y = 0 h x = 0 w
zV ( x , y , z ) ) / ( z = 0 d y = 0 h x = 0 w z ) ( 5 )
##EQU00003##
where V(x,y,z) is voxel intensity, d,h,w are number of slices,
slice height, and slice width. Euclidean distances are computed
from the centroid to each input voxel as E.sub.d= {square root over
((x-C.sub.x).sup.2+(y-C.sub.y).sup.2+(z-C.sub.z).sup.2))}{square
root over
((x-C.sub.x).sup.2+(y-C.sub.y).sup.2+(z-C.sub.z).sup.2))}{square
root over
((x-C.sub.x).sup.2+(y-C.sub.y).sup.2+(z-C.sub.z).sup.2))}. This
distance is normalized as
E.sub.dn=E.sub.d/.parallel.E.sub.d.parallel..sub.max.
[0132] B. Input element 2: Normalized Voxel Density (V.sub.n)
[0133] This value depends on a prior choice of gray scale histogram
normalization or equalization via a lookup table procedure. The
value V.sub.n=f(V(x,y,z))/.parallel.f(V).parallel..sub.max is
calculated from histogram mapping and normalized with the maximum
value obtained from the entire volume. The best results are
obtained via gray scale equalization.
[0134] C. Input Element 3: Gaussian (g)
[0135] A 3D Gaussian is defined as:
G ( x , y , z ) = 1 2 .pi. .sigma. 2 - ( x 2 + y 2 + z 2 ) / ( 2
.sigma. 2 ) ( 6 ) ##EQU00004##
where .sigma. defines the spread of the Gaussian curve as shown in
FIG. 6. Use of two Gaussian features allows the neural network to
produce feature clusters in noisy regions. Pre-computed numbers are
used to normalize the 3D Gaussian. A kernel is computed from this
Gaussian function with pre-specified kernel size ranging from 2 to
10. The smoothed voxel intensity results from the convolution of
the kernel over the present voxel and its immediate three (3)
dimensional neighbors.
[0136] D. Input Element 4: The Difference of Gaussians
(g.sub.d)
[0137] This is an extension of the above Gaussian input element
with dual .sigma..sub.1,.sigma..sub.2 values, shown separately in
FIG. 6. This is defined as:
G.sub.x(x,y,z)=.parallel.G.sub.1(x,y,z)-G(x,y,z).parallel. (7)
where
G 1 ( x , y , z ) = 1 2 .pi. .sigma. 1 2 - ( x 2 + y 2 + z 2 ) / (
2 .sigma. 1 2 ) ##EQU00005##
and G(x,y,z) is defined in equation 6. As in Equation 6,
pre-computed numbers are used to normalize this Gaussian function,
A kernel is computed from G.sub.d(x,y,z) and g.sub.d is (2x)
obtained by convolving this kernel with the present voxel and its
immediate neighbors (i.e., the 6 directly adjacent voxels).
[0138] E. Input Element 5: Laplacian Zero Crossing (Z.sub.c)
[0139] Laplacian Zero Crossings are presented by as:
Z.sub.c(x,y,z)= {square root over
((Z.sub.c(x)+Z.sub.c(y)+Z.sub.c(z)))}{square root over
((Z.sub.c(x)+Z.sub.c(y)+Z.sub.c(z)))}{square root over
((Z.sub.c(x)+Z.sub.c(y)+Z.sub.c(z)))} (8)
this value is obtained from the Laplacian computed in each axis in
the image volume, where:
Z c ( x ) = ( ( m = d - d m = d - d V ( x , y + n , z + m ) sgn ( n
) ) / L ) 2 , Z c ( y ) = ( ( m = d - d m = d - d V ( x + p , y , z
+ m ) sgn ( n ) ) / L ) 2 , Z c ( z ) = ( ( m = d - d m = d - d V (
x + p , y + n , z ) sgn ( n ) ) / L ) 2 ##EQU00006##
and, where sgn(n) is the sign of the scalar x. L is equal to
(2d+1)d where d is a predefined length; and V is voxel
intensity.
[0140] F. Input Element 6: Moment Invariant (m.sub.i)
[0141] The moment invariant is an image statistic that is invariant
to rotation, translation, and scale in a region. Moment invariants
are derived from the definitions of moments, centralized moments,
and normalized central moments. The image moments are defined
as:
m pqr = x y z x p y q z r f ( x , y , z ) ( 9 ) ##EQU00007##
where p,q,r.epsilon.{0, 1, . . . , l}. If f(x,y,z) is a piece-wise
continuous function with bounded support in a 3D image, then the
bounding is set by an m.times.n.times.o voxel array. It is
positioned at the selected voxel and bounded by length l. The
centralized moment .mu..sub.pqr is defined as:
.mu. pqr = x y z ( x - x _ ) p ( y - y _ ) q ( z - z _ ) r f ( x ,
y , z ) ( 10 ) ##EQU00008##
where x=m.sub.100/m.sub.000, y=m.sub.010/m.sub.000,
z=m.sub.001/m.sub.000 (i.e., centroid coordinates) and the point
(x,y,z) is called the image region centroid. The normalized central
moment is defined as: .eta..sub.pqr=.mu..sub.pqr/.mu.'.sub.000
where .gamma.=(p+q+r)/2+1. From the normalized central moments,
three (3) orthogonal invariant moments are defined as:
.phi..sub.1=.eta..sub.200.eta..sub.020.eta..sub.002,
.phi..sub.2=2.eta..sub.110.eta..sub.011.eta..sub.101,
.phi..sub.3=.eta..sub.020.eta..sup.2.sub.101+.eta..sup.2.sub.011.eta..su-
b.200+.eta..sup.2.sub.110+.eta..sub.002 (11)
[0142] These three orthogonal moment invariants were selected
because they are irreducible and distinct. When approaching tissue
boundaries the moment invariants may provide information that
complements the zero crossing data. The strength of the zero
crossing 3D surface signal is expected to be more dependent on the
object's pose (i.e., relation of surface orientation to scan
orientation) than the moment invariants.
[0143] 2. Back projection of Feature Clusters
[0144] In the preferred embodiment, the SOFM neural network
utilizes 4.times.4 to 16.times.16 processor network configurations
in either a rectangular or a hexagonal grid formation (see FIG. 7).
In a rectangular grid, each processor has four connections to its
neighbors. In a hexagonal grid, each processor unit has six
connections to its neighbors. The input layer, translates each
voxel into a vector containing the VIDEs defined in Equation 4. The
input layer has direct connection to all processors in the
formative network configuration and is associated with a weight as
shown in FIG. 5. A flowchart of the complete surface segmentation
process is shown in FIG. 8.
[0145] A. Preprocessing to Remove Background Noise and Head
Holder
[0146] The head holder and background noise are separated from the
head object in a two step procedure. First, a grayscale threshold
is manually set; next, a volumetric region growing method is
applied. This removes voxels associated with the head bolder,
background noise, and most metallic artifact outside the head; it
also limits the voxels sent to the neural network to those relevant
to the 3D surfaces internal to the head.
[0147] B. Randomization of Input Voxels
[0148] The randomization of the VIDE vectors in the input event
space is necessary in order to obtain an ordered map. Each voxel
intensity in the image volume is stored along with its position
(i.e., x, y, z coordinates) in a randomly selected entry in a hash
table (i.e., a register queue length N). The random position
R.sub.p in the hash table is computed as
R.sub.P=hashTable[length.times.randNumber], where the hash table is
initialized sequentially from 0 to the total number of voxels N.
The random number is generated in a range between 0.0 and 1.0. The
length of the hash table starts with value N. Subsequently, N
decrements for each position computed and is swapped with entries
for that position's current length in the hash table. This
procedure ensures randomization of feature vectors in the input
element event space and, subsequently, the map formation by
adapting weights in the self organization process. The position of
the intensity in the hash table entry is replaced during map
formation with the processor identity; processors represent feature
clusters of potential anatomical surfaces.
[0149] C. Self-Organizing Feature Map Algorithm
[0150] With reference to FIGS. 8 and 9, weights are randomly
initialized around 0.5.+-.x,x.epsilon.{0.0, . . . , 0.1}, in a step
K, with default initial learning rate and neighborhood excitation
parameters, set by the user.
[0151] A pattern .xi., a VIDE vector, is presented, and the network
output is evaluated in a step L.
[0152] The unit {i,j} with the minimum output is selected and
designate as a winner:
|P.sub.i,j|.sub.winner=.parallel..xi.-W.sub.i,j.parallel..sub.m- in
in a step M.
[0153] Weights are updated, in a step N, using a learning rule
where neighbor excitement levels are a function of self and
neighbor excitation rates
W i , j ( t + 1 ) winner = W i , j ( t ) + .alpha. i , j ( t ) [
.xi. - W i , j ( t ) ] , W m , n ( t + 1 ) = { W m , n ( t ) +
.beta. i , j ( t ) [ .xi. - W i , j ] , , if ( ( m , n ) .di-elect
cons. N i , j ( t ) W m , n ( t ) , ##EQU00009##
where N.sub.i,j(t) is a neighborhood of {i,j} at time t;
.alpha..sub.i,j(t), the learning rate for the winner processor;
.beta..sub.i,j(t), is the rate at which its excitation influences
its neighbors.
[0154] The value of .alpha..sub.i,j(t) and .beta..sub.i,j(t) is
decreased, and the neighborhood, N.sub.i,j(t) is shrunk in a step
O.
[0155] Steps K-O are repeated for all the VIDE vectors:
W.sub.i,j=[w.sub.1,w.sub.2,w.sub.3, . . . ,w.sub.6]i=1,2, . . .
,Ij=1,2, . . . ,II.times.II.epsilon.{4 . . .
16}Q.sub.i,j=.parallel..kappa.F.sub.n-W.sub.i,j.parallel.
.xi.=.kappa.F.sub.n'
n.epsilon.{1 . . .
N}NQ.sub.i,j.alpha..sub.i,j(t).beta..sub.i,j(t).beta..sub.i,j(t).gamma..s-
ub.i,jid.epsilon.{0 . . . 255}W.sub.i,jN.
[0156] 3. Interactive Image Segmentation and Rendering
[0157] The SOFM 3D surface segmentation program integrates the SOFM
algorithm, volumetric image reading, slice image display and
editing, and graphical surface rendering of output feature cluster
voxel data (i.e., obtained 3D surfaces).
[0158] The 2D SOFM processor map displayed in the window (FIGS. 10,
11, and 12) allows the user to select voxel clusters associated
with each processor. Clicking on a displayed node activates back
projection of the associated voxels to a 3D graphics window.
Successive feature clusters associated with desired surfaces may be
accumulated, saved as a group, and back-projected to the graphics
display (FIG. 10). At this stage (i.e., selection of map nodes) no
surface rendering is provided.
[0159] Our tests were with 3D head CT data. The voxel display of
these data allows the user to quickly identify the processor nodes
associated with skull, soft tissue face, oral, or cervical
vertebral surfaces. After having identified the group of processors
associated with the desired 3D surface, the voxels are triangulated
for 3D surface rendering.
[0160] A. Self Organization (Ordering) Process
[0161] During the automated image segmentation process an ordering,
from an initially random configuration of processors to a state of
equilibrium of the self organizing feature map must occur. The
weight vector of a processor is given by W.sub.i,j=[w.sub.1,
w.sub.2, w.sub.3 . . . w.sub.6], where i=1, 2, . . . , I, j=1, 2, .
. . , I, and I.times.I, is the number of processors assigned to the
network configuration, where I.epsilon.{4 . . . 16}. A maximum of
256 processors were chosen, a condition which insures that the
identity of a processor is always found in the range of 0 to 255.
This allows each processor to be stored within one byte, thereby
saving memory for further computation.
[0162] A VIDE is generated from each voxel in the input volume
event space hash table and is sent to the processor layer. Once a
processor captures a voxel, the cluster identity is placed back in
the hash table register replacing the input voxel intensity cell
datum. Later it will be displayed as part of a surface with all
voxels belonging to this cluster. Voxel intensity information is
obtained from the original input image volume via corresponding
hash table registers. A processor's computed output is:
Q.sub.i,j=.parallel..kappa.F.sub.n'-W.sub.i,j.parallel.(12)
where .xi.=.kappa.F.sub.n' is a VIDE vector, n.epsilon.{1 . . . N}
and N is the total number of voxels. Among the processors a winner
is identified. The distance measure Q.sub.i,j determines if it is
lower than the initially chosen value. If so, the processor
identity changes from its initial state; the new processor will
continue to claim the voxel as long as no other processor claims it
is the "winner." The winner successively strengthens its weight
vector by a fraction of the input element values sampled by the
.alpha..sub.i,j(t), a self excitation parameter. The value of
.alpha..sub.i,j(t) is set and then decreased successively,
eventually resulting in no change to the weight vector. This
ensures that the network will reach an equilibrium. Another
parameter, .beta..sub.i,j(t), influences neighborhood processor
excitation and determines how much the winner can excite or inhibit
its neighbor's weights. .beta..sub.i,j(t) is assigned a fixed value
that decreases over time to avoid wide fluctuations in the
organization process. Using neighborhood Gaussian function curves
(See FIG. 3), defined by the parameter .gamma..sub.i,j, the number
of neighboring processors is successively decreased. This number is
initialized with a default value at the beginning. The "winner"
processor takes the current input voxel to its cluster along with
its spatial coordinates. This update replaces the processor's
identity, id.epsilon.{0 . . . 255}, with the voxel density
currently in the register.
[0163] Prior to the mapping process the weight vectors, W.sub.i,j
of each processor, are initialized with random values. Therefore,
the initial topology of the axis spanned by connection strength
between any two of the VIDE elements is unorganized. This is shown
in FIG. 13A. The processor's physical connections remain constant
throughout the process, (i.e., processor 1's connections to
processor 2 remain the same). However, the mapping topology with
respect to each processor's connection strength changes throughout
this process. The random starting topology, shown in FIG. 13A, is
drawn by plotting each processor's weight vector, chosen two weight
elements at a time, onto two dimensional axes. The initial pattern
has no meaning, however the final pattern obtained presents
processors (feature clusters) containing voxels representing 3D
surfaces and similar 3D surfaces will be shown in the same SOFM map
region. We have chosen voxel intensity and Euclidean distance and
their corresponding members in processor weight vectors to plot the
map topology in FIGS. 13, 10, and 11.
[0164] After N input VIDEs are fed through the input layer,
producing each processor's final weight vector, the network adopts,
approximately, the original "flattened" processor grid topology.
The axis of the grid is spanned by processor connection strength
values. This grid is a 2D projection image of the feature cluster
VIDE hyperspace. An example is seen in FIG. 13B.
[0165] SOFM feature vector adaptation results from two basic
mechanisms: (1) Identifying a "winner" and allowing it to modify
it's connection strengths; (2) allowing the winner to selectively
excite/inhibit its neighbors. The first mechanism ensures
connection strengths that reflect input event values. The second
mechanism ensures that related processors, especially those
representing anatomical surfaces, are mapped adjacently. The
observed surfaces, displayed as point clusters from a set of
neighboring processor nodes, illustrate the second mechanism in
FIG. 10.
[0166] The voxels associated with the winners are clustered and
classified as separate groups. These clusters are identified
through the processor's identity, a unique value which is stored in
the hash table during the self organization process. The resulting
clusters capture continuous anatomical surfaces [FIG. 10] that
exist in the 3D image volume as an intrinsic model of the VIDE
event space. This approach differs from physically based models
(prior information) or user-defined (algorithm assisted or manual
selection) segmentation methods that rely on information extrinsic
to the image volume. Our method adaptively locates surfaces
captured in the image volume using voxel data intrinsic to the
image volume. We involve the user to verify and determine the
surface identity. This introduces supervised learning to the
unsupervised neural network algorithm.
[0167] B. Surface Rendering
[0168] Two different types of surface rendering methods (i.e.,
triangulation of voxels associated with selected feature clusters)
have been implemented. The first approach is primarily a contouring
of the outermost voxels. Assuming the volume contains a thin
boundary surface, all cells on the surface are classified as useful
vertices. Each slice is divided into four quadrants, and each
quadrant is swept by 256 lines from it's outer most volume boundary
to the slice image centroid. In each ray, the first encounter with
the boundary surface is defined as a boundary vertex. The remaining
sweeps obtain subsequent vertices in each slice, resulting in a
closed curve that traverses all of the detected slice vertices.
Since the fine sweeps are produced at half degree angular rotation
intervals, the points generated are aligned in the correct
neighborhood. Next, four sided polygons are generated from all
quadruples of neighboring vertices in a pair of consecutive closed
curves. These polygons are degenerated to two triangular faces.
This results in a complete triangulated surface of the segmented
image object. This procedure is extremely fast and well suited to
display simple surfaces such as the soft tissue face (see FIG.
12A). However, views of objects with complex topology, such as the
skull, provide unsatisfactory results.
[0169] A second approach is taken to rendering complex surfaces
such as the skull. We used a known method of decomposing volumetric
data into tetrahedra. This approach takes two adjacent slices and
traverses neighboring voxels eight at a time, treating them as a
tetrahedron. The decomposition of tetrahedra produces the oriented
cycles necessary to represent triangular faces of an isosurface.
This approach yields, for example, a skull surface (see FIG. 12B).
We also compared the Doi and Koide tetrahedral decomposition to the
marching cubes algorithm. Both algorithms use bi-linear
interpolation to compute partial voxel values. Computation of a
complete skull surface rendering took approximately 3 to 4 minutes
as opposed to 2 to 3 minutes, on an SGI Octane workstation.
[0170] C. Combining Multiple Surfaces and Filling Holes
[0171] The user can combine the feature clusters associated with
multiple SOFM processors. If unexpected holes are seen in the
surface, perhaps due to partial voluming "drop out" of thin
structures, the volume image region immediately around the missing
area may be selected as a new region of interest (ROI). A new high
resolution processor map of voxel clusters is generated for this
zoomed-in ROI, allowing careful surface assembly in regions of thin
bone or skin.
[0172] D. Partial Volume Error
[0173] This error results in surface "drop out" and is a well known
cause of surface inaccuracy in 3D CT visualization of thin bone
segments. The problem becomes more difficult where bone voxel
intensities are reduced to tissue levels. Region growing and manual
intensity thresholding techniques for isosurface generation cannot
overcome these errors. We find that selecting a region presenting
thin bone segments and assigning higher preferences to VIDEs 5 and
6 (i.e., Laplacian zero crossing and Invariant image moments)
allows a local SOFM map to capture more of the desired voxels
representing surfaces formed by thin bone. For example, the
computationally expensive moment invariant was set at 10 weighting
on this initial pass, then elevated to 40-60% in user-defined
regions around thin bone structures such as the infra-orbital
foramina and orbital walls (FIG. 14).
[0174] 4. SOFM Precision Tests and Results
[0175] The SOFM program was tested with five 3D CT data sets
(B2537, B2621, B3037, B3095, B3195) provided by a recall of
subjects. All subjects are above 60 years in age, female, have
dental fillings, and often have large oral prostheses (dentures).
One operator (AM) segmented all the images twice by two methods:
(1) slice-based methods (i.e., intensity thresholding, manual pixel
selection, and 2D and 3D region growing) in the AUTO program, and
(2) Self Organizing Feature Map. Both the slice-based and SOFM
tools provide orthogonally reformatted views of the original image
and segmented pixels. More than a day passed between each
segmentation session. Intra-operator variance (precision), using
either segmentation method, was measured between the two
segmentation attempts. In manual segmentation, skull and soft
tissue face segmentations were separate activities. In SOFM
experiments, the same ordered map was used to produce the skull and
soft tissue surfaces in a single attempt.
[0176] The original, unsegmented CT data sets are
512.times.512.times.115 in dimension with 1.904 mm slice thickness
and 0.5 mm pixel resolution. This scan protocol was chosen to
provide accurate source data for boney skull and soft tissue face
isosurface images. This data was reduced to 256.times.256.times.115
in dimension with 1.904 mm slice thickness and 1 mm pixel
resolution. This resolution is sufficient in terms of the resulting
surface segmentation and saves memory resources on the
computer.
[0177] Five measures were used to determine both intra-operator
within-method variance and between-method variance in the resulting
segmented images. The former measures how reproducible the
segmentation is; low variance implies high precision. The latter
measures how different the results of the two methods are. We
expect the differences to be local to regions of thin structures or
metallic artifact where qualitative inspection suggests the
greatest benefit from SOFM segmentation. These five measures are:
(1) Visual comparison of the skull and soft tissue face surfaces
obtained; (2) Difference of the overall volume of the segmented
image; (3) Two dimensional difference images computed slice by
slice for each slice image within the volume; (4) Three dimensional
difference images viewed from three orthogonal eye points; and (5)
Procrustes superimposition of manually located, highly reliable,
single point anatomical landmarks. The last measure is similar to a
qualitative visual determination of the completeness of anatomical
features seen on the segmented surface. Clearer features will
result in more easily, and more repeatably, detected anatomical
landmark coordinates by trained workers.
[0178] A. Test 1: Visual Comparisons: The left hand side of FIGS.
14 and 15 show both manual slice-based segmentation attempts. The
right hand side shows both SOFM segmentation results. SOFM
segmentation of soft tissue face recovers thin tissue structures
around the nostrils, lips, and ears, producing an improved image.
There is also dramatic improvement in reducing the effect of the
metallic (dental filling) artifact around the lips. Note the
variable posterior sagging of inframandibular soft tissue due to
the affects of gravity on the supine individual within the head
holder. SOFM skull segmentation (FIG. 14) improves areas of thin
bone around the infraorbital and mental foramina and the thin
posterior and medial orbital walls. These visual results indicate
that the slice-based segmentation does not allow as effective a
search for these surfaces.
[0179] B. Test 2: Volume Differences: Table II presents volumetric
differences between the soft tissue face obtained between an
average of multiple attempts and each attempt in both manual and
SOFM methods and between the two method averages. Table III
presents the same data for the segmented skull image volumes. The
difference image volumes were computed by multiplying the total
number of voxels in each segmented image by voxel volume (1.904
mm3) SOFM segmentation results show less difference between
sessions. Given the additional detail it is not surprising that in
comparison with manual slice-based segmentation, SOFM segmentation
results in an average of 3.9% more volume for skull segmentation
and 0.5% more for soft tissue face segmentation. The higher
difference for skull segmentation likely indicates greater
reduction of partial volume error.
[0180] C. Test 3: 2D Slice Image Difference Data: This test
presents 2D slice by slice difference profiles for each image
volume using either segmentation method. These difference profiles
result from calculating the absolute difference between
corresponding pixels in the homologous segmented slices. The
accumulated difference across all five data sets are plotted per
slice in FIG. 16A for the soft tissue face surface and FIG. 16B for
the skull surface. Regions where the slices contain thin bone and
dental filling artifact require manual pixel selection during
manual slice-based segmentation in AUTO. These are areas where the
user must make many decisions that are not well replicated on
subsequent segmentation attempts. Manual slice-based segmentation
results in higher intra-operator variability. Regions with thin
soft tissue or bone, or dental artifact, contributed most to the
difference plots (FIG. 16) between the slice-based manual and SOFM
segmentation results, this suggests that the latter was more the
reliable method.
[0181] D. Test 4: 3D Difference Images: This test consists of 3D
difference images computed between the volume image resulting from
the two segmentation methods. The 3D difference images, seen in
three orthogonal projections of 2D images in FIG. 17, graphically
illustrate the affect of highly variable human decision making in
the regions of metallic artifact and thin structures during manual
slice-based segmentation. Note that thin areas of the skull
surround the pneumatic sinuses (i.e., frontal, maxillary, ethmoid,
sphenoid, and mastoid). Areas of thin soft tissue face structures
(especially, the nostrils, eyelids, ears, and lips) are more fully
extracted by SOFM than with slice-based manual methods.
[0182] E. Test 5: Anatomic Landmark Localization: The user was
asked to pick reliable anatomic landmarks twice on the soft tissue
face and skull surface segmented either manually or by SOFM. These
landmarks include 21 on the soft tissue face surface (Table IV) and
26 on the skull surface (Table V). The error found between landmark
locations in these two sessions was determined via Procrustes fit.
The error reported here is the average value across all 5
superimpositions. We contend that the improved quality of the
imaged surface, especially in areas of artifact or thin structures,
reduces intra-operator variance (i.e., raises precision) between
landmarks picked on two attempts of either the soft tissue face or
skull surfaces segmented via the SOFM method; the less ambiguous
SOFM surface renderings allow our skilled operator to place these
landmarks more precisely.
[0183] F. Time of Segmentation: Skull and soft tissue face
segmentation sessions (AM) were reduced from approximately 3 to
approximately 0.5 hours. In either method, the majority of time is
spent correcting dental artifact or partial volume error. We are
aware that the skull is most often segmented at the CT console by
intensity thresholding alone. This method may save time but results
in the poor segmentation of thin structures or those near metallic
artifact.
II. Warping
[0184] As described below, the step G of warping the template to an
external surface of a normative shape of the bone of interest is
preferably executed using a Simulated Annealing-based Surface
Extraction (SASE) process.
[0185] In order to save time, improve precision (repeatability),
and increase accuracy of the extracted anatomical surface, a
ridge-curve based deformable template superimposition strategy
combines simulated annealing of a ridge and tiling curve wireframe
and subsequent surface normal based identification of surface tile
points. The annealing algorithm searches the volume image and
displays the result to the operator on a graphical manifold
rendering of the surface of interest. Unlike energy minimization
algorithms, which require a good initial parameter estimation to
avoid converging to a misleading local minimum, annealing methods
exploit stochastic search methods to locate a global minimum.
Several stochastic methods form a basis for simulated annealing.
Space curve fitting of 3D surfaces is a goal not uncommon in the
field of computer vision. It is merely object recognition followed
by shape registration.
[0186] B-spline space curve encoding of labeled surfaces has
application beyond surface averaging and craniometrics or image
co-registration. We seek highly accurate surface extractions to
model cranial prosthetic implants for the skull. Homology mapped,
B-spline space curve encoded surfaces prove useful in morphometric,
biomechanical, and deformable model assessment of many organ
systems. Additionally, parameterized surface extractions,
especially of the external body surfaces, have application in
animation.
[0187] 1. Ridge Curve-Based Deformable Template Creation
[0188] We chose to use B-spline surface encoding because it
provides a smooth surface approximation, reduced data set, and
easier assignment of labels to various components (see Appendix 2
for B-spline curve and surface construction). We have
re-implemented skull and soft tissue face ridge curve-based
deformable templates in this fashion. FIG. 18 shows one of these
new ridge curve-based deformable template wireframes superimposed
on a soft tissue face graphical manifold, with Type II landmarks,
ridge curve arcs, and tiling curves indicated.
[0189] Tiling Curves: Since ridge curves rarely intersect,
additional space curves are necessary to tile the surface with a
single enclosing wireframe. In order to identify a unique
demarcating line within a search space we define tiling curves as
space curves that: (a) traverse a low curvature area (i.e., flat
surface); and (b) are weighted to pass close to a plane passing
through the surface between the tiling curve end points and the
center of mass of the entire image surface (See Sections II.2.C and
II.2.G).
[0190] Initially, we approximate tiling curves as straight lines
between endpoint landmarks. In FIG. 18, B-spline ridge and tiling
curves are indicated by dotted lines while their control points are
shown as small black cross-hatches. The Type II landmarks are shown
as large black cross-hatches.
[0191] As one draws either a ridge or tiling curve, a corresponding
B-spline curve is fit to the surface via global curve interpolation
to a set of control points sampled from the user drawn space curve.
The B-spline interpolated 3D points collected are re-sampled and
indexed to vertices on the graphical manifold surface and assigned
a control point vector, P.sub.1. In the case of ridge curves, our
resampling maintains the curvature information along ridge curves
when determining how to space the control points. On the other
hand, tiling curve control points are equi-distantly spaced. The
curve smoothness, set by the operator, further determines
redistribution of control points, if necessary. Knot-vectors and
knotspacing require a B-spline interpolation function assignment
for each space curve, by the operator during the template creation
process.
[0192] FIGS. 19A and 19B display the Type II landmarks and Appendix
1A and 1B list the non-landmark components of the soft tissue face
and boney skull templates, respectively. The soft tissue face
template consists of 56 Type II landmarks, 103 space curves (54
ridge curves and 49 tiling curves), and 44 surface tiles [FIG. 19A,
Table VI]. The skull template consists of 54 Type II landmarks, 104
space curves (40 ridge curves and 64 tiling curves), and 47 surface
tiles [FIG. 19B, Table VII].
[0193] 2. Ridge Curve-Based Deformable Template Superimposition
[0194] A. Graphical Surface and Volume Image Input: The 3D surface
input we use results from segmented volume images. The segmented
voxels are rendered as a graphical manifold. Segmentation occurs
via an unsupervised neural network program. The graphical manifold
is built with sub-sampled voxels to provide an interactive context
to the operator.
[0195] The simulated annealing ridge curve-based deformable
template superimposition algorithm evaluates features found in the
volume image, not on the graphical manifold surface. The
unsubsampled volume image provides additional cues for best
locating ridge or tiling curves. The result of the search is
displayed on the graphical manifold surface where the operator can
make manual corrections to algorithmic insufficiencies. The CT
slice shown in FIG. 20A highlights the segmented surface voxels of
interest which are indexed to the vertices (points) in the adjacent
rendered manifold surface (FIG. 20B).
[0196] During loading of the graphical manifold surface data, a 3D
lookup table (i.e., geometric hash-table) of spatial x,y,z and
curvature metrics (Listed in Table 8) corresponding to each vertex
is generated. This provides the simulated annealing based surface
extraction complete indexing between the rendered graphical
manifold and the source volume image.
[0197] B. Candidate Wireframe: The first step in the Simulated
Annealing-based Surface Extraction (SASE) process is the operator's
manual location of the Type II landmarks on the graphical manifold
surface. These landmarks attach the ridge curve-based deformable
template to the graphical manifold surface via a thin plate spline
warp. FIG. 21A shows a candidate wireframe prior to the warp at
lower left. FIG. 21B shows the result of the thin plate spline warp
of the wireframe onto the Type II landmarks. Some of the landmarks
(large cross-hatches) are visible. Note how the candidate wireframe
undulates in and out of the surface.
[0198] C. Volume Image Features: Simultaneous to loading of the
graphical manifold data and lookup table generation, a set of
surface geometry features are generated for use in the simulated
annealing procedures. A method is known for computing principal
directions and principal curvature values (FIG. 22) for any point
on the segmented image surface directly from the input volume image
data. In addition, the principal normal and tangent vectors are
derived at the same time. These latter parameters, in addition to
providing the basis for the simulated annealing process cost
functions, are also used to bound the search space in the volume
image linked to the graphical manifold surface displayed for the
operator.
[0199] The search for volume image principal curvature normal
vectors and principal curvature tangent vectors occurs initially at
a point, I, a control point found on the B-spline space curve,
defined at location {x.sub.0,y.sub.0,z.sub.0} to have normal vector
{right arrow over (N)} on the volume image surface S defined as
I(x,y,z)=I(x.sub.0,y.sub.0,z.sub.0). I lies in the same direction
as the curvature gradient (i.e., changing surface normals) of the
surface voxel data:
.gradient.I=(I.sub.x,I.sub.y,I.sub.z)|(
N=(.gradient.I)/.parallel..gradient.I.parallel.) (EQ 13)
[0200] In order to evaluate the local surface curvature we consider
a space curve C from a sweeping plane positioned at each control
point, with tangent vector {right arrow over (t)}, along a normal
section of segmented voxel surface S (FIG. 22A). The normal vector
at any point along C is N=(n.sub.x,n.sub.y,n.sub.z).
Differentiation of .gradient.I.t=0 yields surface normal curvature
k.sub.n in direction {right arrow over (t)}, where k.sub.n attains
principal curvature values k.sub.1 and k.sub.2. Principal
directions {right arrow over (e)}.sub.1 and {right arrow over
(e)}.sub.2 show the path of the principal curvatures (FIG.
22B).
[0201] On any surface, the local curvatures vary between the
principal curvatures k.sub.1 and k.sub.2. Therefore, the curvature
parameters we use combine these measures:
K=k.sub.1k.sub.2 (EQ 14)
M = 1 2 { k 1 + k 2 ) ( EQ 15 ) ##EQU00010##
where K is Gaussian curvature, independent of principal directions
and M is mean curvature changes sign with opposing principal
directions.
[0202] Twelve feature dimensions are used for simulated annealing
registration of ridge or tiling curves (Table VIII). The 3rd order
image derivatives are obtained from image derivatives:
H = [ I xx I xy I xz I yz I yy I yz I zx I zy I zz ] ( EQ 16 )
##EQU00011##
[0203] The Partial derivatives of H along each axis are
H.sub.x,H.sub.y,H.sub.z:
H x = [ I xxx I xyz I xzx I yxx I yyx I yzx I zxx I zyx I zzx ] , H
y = [ I xxy I xyy I xzy I yxy I yyy I yzy I zxy I zyy I zzy ] , H z
= [ I xxz I xyz I xzz I yxz I yyz I yzz I zxz I zyz I zzz ] ( EQ 17
) ##EQU00012##
[0204] Among these 27 third order partial derivatives, ten
derivatives are chosen, I.sub.xxx, I.sub.xxy, I.sub.xxz, I.sub.xyy,
I.sub.xyz, I.sub.xzz, I.sub.yyy, I.sub.yyz, I.sub.yzz, I.sub.zzz,
preferring derivatives along the three axes. Essentially, the
maxima of these ten third order partial derivatives provide an
indication of sharp geometric transition. A pre-determined
inter-knot spacing restricts any overlap of search space between
any two control points. The Type II landmarks are used for ridge
curve arc-separation and to anchor tiling curves.
[0205] Features for Ridge Curve Arcs: We use Gaussian and mean
curvatures, and ten third order partial image derivatives to obtain
maximum curvature prescribing control points along the surface
between ridge curve arc endpoints (i.e., Type II landmarks).
[0206] Features for Tiling Curves: We use the same set of features
to position tiling curve control points along lines of minimum
curvature. In addition, we use a thirteenth feature, a deviation
measure, d.sub.m. It is found between the candidate control point
normal and the nearest ray within a plane constructed from the two
endpoint landmarks and the object's center of mass (FIG. 23). The
deviation measure helps limit the tiling curve search space so that
no boundaries overlap.
[0207] D. Initial Control Point Location: When the template is
warped to operator located landmarks on a patient image the
candidate space curves connecting the Type II landmarks and their
control points may not lie on the image surface. The true surface
is located using a trace along the candidate control point normal.
Thus, the initial candidate control point is assured to be on the
image surface. FIG. 24 illustrates a trace along the candidate
template's control point to the image surface. The annealing
process starts with the first point found on the surface as an
initial candidate control point.
[0208] E. Search Plane: A search plane is anchored to the candidate
wireframe (i.e., ridge curve-based deformable template) and rotates
about the axis of that point's surface normal in the volume image.
Note that when ridge or tiling curves sit on a two sided tile
boundary, these cut planes rotate a full 360 degrees (FIG. 25B).
The sweep rotation about each control point normal is set to
increment a default 45 degrees (operator chooses between 10 and 90
degrees). If search in one direction fails, an opposite direction
search is attempted. This is important where a ridge arc is at the
surface margin.
[0209] F. Random Candidate Control Point Generation: The principal
directions constrain the search space for each control point along
the B-spline space curves. The maximum principal curvature bounds
the space in which it is profitable to search for optimal B-spline
space curve control points. The search space for each B-spline
control point during ridge curve or tiling curve superimposition is
defined from the point normals and tangents bounded by the
principal directions. FIG. 25A shows a cross section on the top of
the skull and tiling curve control point search planes intersecting
that image surface. FIG. 25B shows a series of search plane
iterations.
[0210] Candidate control points are evaluated during each search
plane iteration. Search plane rotations select a randomly located
set of M (chosen by operator, default is 2) candidate control
points for k search iterations. Random candidate positions are
selected from the space curves intersected by the search plane
rotations through the image surface. Their features (Table VIII)
are computed and evaluated by the cost function in comparison to
the current candidate control point position, .zeta.. The actual
candidate locations on the image surface chosen for evaluation via
variables r and s using reflection coefficient parameter
.alpha.:
r = ( .alpha. N 2 ) ( ( N + 1 ) + N - 1.0 ( EQ 18 ) s = ( .alpha. N
2 ) ( ( N + 1 ) - 1.0 ( EQ 19 ) ##EQU00013##
where N is the number of dimensions (i.e., number of image features
in the case presented). Randomly distributed position vectors
(I.sub.p) are computed by using r and s as:
Ip=S(hrand(r),wrand(s)) (EQ 20)
where S is the search space spanned by width, w, and height, h,
computed from the normal and tangent at the template B-spline space
curve control point. The rand(r),rand(s) are the coordinates with
the search plane. Each position within the search space plane
carries a single vector for the three coordinate axes, which is
indexed to a segmented surface voxel in the volume image space.
[0211] The probable distribution of positions are subsequently used
in finding the best ridge or tiling curve control point candidate.
The values for r and s are decreased after each iteration:
r=r-(r/k),s=s-(s/k) (EQ 21)
and successively decreased, as k, the current iteration, approaches
the maximum number of iterations, k.sub.max.
[0212] G. Simulated Annealing Cost Functions: The cost function
evaluates the candidates for B-spline control points intersected by
the search plane via simulated annealing. The procedure used here
is based on the Metropolis algorithm.
[0213] Ridge Curve Cost Function: The value of the ridge curve
search space cost function is computed from the twelve features
(Table VIII):
E ( .zeta. ) = features = 0 12 F .zeta. i - 2 ( EQ 22 )
##EQU00014##
where .zeta. is a candidate space curve control point. Since the
ridge curve requires maximum curvature values, the inverse of these
squared distance feature is minimized during simulated
annealing.
[0214] Tiling Curve Cost Function: For tiling curves, the cost
function is defined as:
E ( .zeta. ) = features = 0 13 F .zeta. i - 2 ( EQ 23 )
##EQU00015##
The curvature-based features, as well as the deviation measure
(d.sub.m) are minimized directly as squared distance features for
tiling curve control point identification (Table VIII).
[0215] H. Illustration of Control Point Location: The rotating
search planes are intended to identify optimal control point
locations to assign to the ridge or tiling curve (FIG. 26A). For
example, search planes for one ridge curve (i.e.,
menton_to.sub.--1_gonion) are arrayed along the warped ridge
curve-based deformable template space curve in FIG. 26B. Here these
search planes pass along the candidate ridge curve control point
normal superiorly to find intersection with the voxels representing
the base of the mandible. FIG. 26C, shows the final ridge curve
registration found by the simulated annealing algorithm. Similarly,
FIG. 26D shows a candidate tiling curve (i.e., r_ASSZ_to_r_PBP) and
its associated search planes. FIG. 26E shows the fitted control
points. FIG. 27 shows the entire final ridge and tiling curve
wireframe superimposed on the image surface.
[0216] I. Simulated Annealing Algorithm Description: Four key
elements are considered in our simulated annealing algorithm
implementation: 1) The configuration definition; 2) Random
candidate generation mechanism (i.e., definition of neighborhoods
on the configuration space); 3) A cost function; and 4) A cooling
schedule. In our implementation, we have first defined the
configuration space with respect to the curvature features
necessary to locate ridge or tiling curves on an image surface
(Section II.2.C). Second, candidate control points are randomly
generated from the search planes (Section II.2.F). Third, standard
cost functions for location of ridge and tiling curves have been
defined (Section II.2.G). Finally, the annealing procedure is
described here.
[0217] The annealing algorithm (FIG. 28) attempts to minimize the
cost function E(.zeta.) where .zeta. is a continuously modified
candidate B-spline space curve control point, which lies in a
search plane intersecting the surface. A set of cost functions are
constructed from the sequence of template control point
deformations .zeta..sup.1, .zeta..sup.2, to .zeta..sup.k starting
from .zeta..sup.0 at the initial ridge curve-based deformable
template control point estimate. For each deformation, we have a
set of P candidates .zeta..sub.0.sup.k, .zeta..sub.1.sup.k, . . . ,
.zeta..sub.p.sup.k within multiple search plane in current
iteration, k. From the initial candidates, the best candidate
{circumflex over (.zeta.)}.sup.k is selected.
{circumflex over (.zeta.)}.sup.k=.zeta. of
.parallel.E.sub.j(.zeta..sup.k).parallel..sub.min (EQ 24)
as k approaches the maximum number of iterations the cost function
E(.zeta.) reaches a stable control point .zeta..sup.s.
Simultaneously, the temperature, T.sub.k, a component of the
annealing procedure, is set to decrease monotonically. This
parameter limits the search area and determines the probability of
determining the best candidate between iterations. As the
temperature lowers, the search area and number of probable
candidates, narrows. At each iteration, k, the algorithm determines
a new value of the B-spline control point deformation parameters,
based on the previous value (i.e, .zeta..sup.k-1) and compares it
with the present value.
[0218] .zeta..sup.k is determined with probability p.sub.k-1, which
are defined as:
.zeta. k = { .zeta. ^ .di-elect cons. p k - 1 .zeta. k - 1 p k - 1
( EQ 25 ) pk - 1 = min ( exp { - [ E ( .zeta. ^ ) - E ( .zeta. k -
1 ) ] T k - 1 } , 1 ) ( EQ 26 ) ##EQU00016##
[0219] T.sub.k is a monotonically decreasing temperature sequence,
E({circumflex over (.zeta.)}) is the cost function computed at the
current position I, and E(.zeta..sup.k-1) is the cost function
computed in the preceding iteration.
[0220] The temperature schedule T.sub.k [34]:
T k = T 0 ( T f T 0 ) k k max ( EQ 27 ) ##EQU00017##
where T.sub.0 is initial temperature, T.sub.f is final temperature,
and k.sub.max is the number of maximum iterations allowed. FIG. 29
presents the temperature cooling schedule from the initial
temperature 100 to the final temperature, -10. We set the number of
iterations to 120 for our initial testing (see Section 4).
[0221] J. Final B-spline Ridge and Tiling Curve Registration: FIG.
30 presents an overview of the ridge curve-based deformable
template registration process in SASE. Manual editing of control
point location is provided in the event that the simulated
annealing algorithm fails to register a ridge curve arc to the
operator's satisfaction. The operator may also locally vary the
search parameters (i.e., sweep rotation degree, search plane height
and width, or annealing temperature).
[0222] 3. Surface Extraction
[0223] The registered B-spline control points are interpolated to
produce a set of tile boundary points (i.e., smooth ridge or tiling
curves). We set the default number of tile boundary points for each
space curve to 50, assuming this would be adequate in all cases.
Each tile, therefore, consists of 2500 parametric points. Our soft
tissue face ridge curve-based deformable template consists of 44
surface tiles [Tables VI and VII], 106376 tile points. Our skull
ridge curve-based deformable template consists of 47 surface tiles
[Tables VI and VII], 113334 tile points. Note that some of the
tiles are constructed with three bounding curves and the fourth
side is a degenerate endpoint of the first and third sides. Since
they are redundant, we can reduce the number of tile points encoded
somewhat from the fully quadrilateral tile case of 110000 (soft
tissue face) and 117,500 (skull).
[0224] A. Space Curve Registration To Volume Image Surface: The
trace of I.sub.n={I.sub.nx,I.sub.ny,I.sub.nz} provides the path
along each tile boundary point normal to the voxel representing the
segmented surface in the volume image. The position of the detected
volume image surface point along the directed normal replaces the
computed space curve point. FIG. 30A shows tiling space curve
points registered onto the volume image surface via reference to
the linked graphical manifold image. The surface normal at each
point is seen as a single perpendicular line.
[0225] For each computed tile boundary point in a ridge curve or
geodesic B-spline, a nearest point is located on the segmented
surface voxels. The surface normal vectors, at current curve point
I, are taken as the average of the local normal and interpolated
end point normals:
I.sub.nx=(((I.sub.x-I.sub.cx)/xsize)+(e1.sub.nx(1-p)+e2.sub.nxp))/2
(EQ 28)
I.sub.ny=(((I.sub.y-I.sub.cy)/ysize)+(e1.sub.ny(1-p)+e2.sub.nyp))/2
(EQ 29)
I.sub.nz=(((I.sub.z-I.sub.cz)/zsize)+(e1.sub.nz(1-p)+e2.sub.nzp))/2
(EQ 30)
where xsize, ysize and zsize are the width, height and slices in
the volume image. I.sub.cx, I.sub.cy, I.sub.cz is the midpoint of
each coordinate, computed during the loading of a graphical
manifold produced from the external tile boundary points. The
I.sub.cx, I.sub.cy, and I.sub.cz are:
I cx = 1 T i = 1 T V ix , I cy = 1 T i = 1 T V iy , I cz = 1 T i =
1 T V iz ( EQ 31 ) ##EQU00018##
where e1.sub.n and e2.sub.n are the tiling space curve end point
(i.e., Type II landmarks) normals, and p is the distance ratio of
the current position I to the nearest space curve end point.
[0226] B. Interior Tile B-spline Space Curves Extraction: The
internal surface tile point vertices are determined via bi-linear
Coons surface interpolation (equation 16). The Coons surface is
converted into set of interior tile B-spline space curves along u
and v parametric directions. The normals of each surface tile point
are computed (equations 16, 28, 29, and 30). These normals are used
to index the tile points to the volume image surface voxels (FIG.
30B). We compute the normal at each parametric surface tile point
by directing a trace from the u, v location calculated from the
Coon's surface to the voxels representing the segmented surface in
the volume image. The detected volume image surface point replaces
the Coons surface tile point. FIG. 30C shows the extracted surface
as it is displayed on the graphical manifold for the operator's
approval.
[0227] 4. Precision and Accuracy Tests and Results
[0228] The SASE 3D surface extraction program was written in C in
the Silicon Graphics Unix (IRIX) environment and tested on a SGI
Octane workstation (R10000 CPU, 256 Mbyte RAM). User interaction
occurs via X11/Motif protocols.
[0229] We tested the SASE program with five 3D CT data sets (B2537,
B2558, B2621, B3095, B3195) provided by a recall of subjects to the
Bolton-Brush Growth Study Center at Case Western Reserve
University. All subjects were females above 60 years in age. The
boney skull and soft tissue face surfaces were segmented in a
separate program. We measured the precision of our ridge
curve-based deformable template registration and surface extraction
by repetition. We compared both our surface extraction and that
obtained by the toolkit to the segmented surface voxels in the
volume image to test inter-method accuracy.
[0230] A. Precision Tests: Intra-operator precision was tested by
repeating the SASE template registrations and surface extractions.
Direct differences between the two surface extractions are
projected into three orthogonal views (FIG. 31). Although,
differences were minimal, most were localized to the oral and
orbital regions. The oral region imprecision may be due to dental
artifact. The orbital region imprecision may be due to surface
ambiguity in thin bone structures. Operator manual effort in
correcting ridge curve extraction varied requiring between 10%
(soft tissue face) and 30% (skull) of the total surface extraction
time in the SASE method. No tiling curve correction was necessary.
We also computed the square root mean of sum of squared distance
between all tile points (Table IX).
[0231] B. Inter-Method Accuracy Tests: 1. Visual Comparison: The
left side of FIG. 32A shows the soft tissue face surfaces obtained
via the SASE methods versus the toolkit on the right side. FIG. 32B
presents the same comparison for the boney skull.
[0232] 2. Computed Differences: The same segmented surface voxel
volume image data were used for template registration and surface
extraction by both the toolkit and the SASE program. This provides
a baseline to compare how accurately the extracted surfaces match
the volume image surface.
[0233] The toolkit generates an initial graphical manifold via the
Wrapper and then converts it to Alligator/Winged Edge format. The
Toolkit's Alligator/Winged Edge graphical manifold is subsampled,
first, during surface rendering from the segmented voxel step,
introducing an error on the order of .+-.1 to .+-.3 mm for a
constant sub-sampling rate of 2.times.2.times.2. Second, a
wireframe is nudged into best position on this graphical manifold
by interpolation. This step increases the error at each point by
.+-.1 to .+-.3 mm depending on the variable level of parametric
smoothing that the program automatically applies. The SASE surface
extraction program indexes (links) the vertices of a less
sub-sampled graphical manifold, used for operator interaction and
approval, to the original segmented volume image voxels, where
wireframe superimposition operations begin.
[0234] In our second inter-method test, the tile points in surfaces
extracted by both methods were mapped to the segmented voxels in
their shared source volume image space. The distances between the
extracted surfaces and the nearest surface voxel in the volume
image along the surface tile point normals were computed. Table X
lists the square root of sum of squared distances for each method
for all five CT data sets and both surface extraction methods. The
surfaces, deviated from the original image voxels by an average of
6 mm for both the skull and soft tissue face surfaces. The SASE
method average error was 0.47 mm for the soft tissue face and 0.46
mm for the boney skull (i.e., pixel resolution). FIG. 33A shows
color coded difference images of the skull and soft tissue face
surface extraction comparisons with the segmented voxel surface.
Similarly, FIG. 33B shows the same color coded difference images of
SASE skull and soft tissue face surface extractions. The error is
scaled from lower difference in blue to larger differences in red.
Note that despite higher sampling density the internal tile points
are as true as the tile borders in the SASE extractions. Table XI
presents the same results tile by tile. Note that the larger tiles
have the greatest error in the toolkit extractions, especially on
the calvarium.
[0235] C. Other Inter-method Differences: 1. Tile Point Density:
FIG. 34 shows the tile point density. The toolkit extracts fewer
points than the SASE program.
[0236] 2. Surface Patch Continuity: Continuity of surface
curvatures across surface tile boundaries is a critical rendering
issue for parameterized ridge curve-based deformable template
surface extraction. The toolkit interpolates surface normals
providing smooth boundary transitions (i.e., no seam). FIG. 35
displays the toolkit and SASE surface extractions. The SASE program
avoids interpolation, using the unaltered, detected segmented voxel
surface normals.
[0237] 3. Tendency of Toolkit to Enlarge Ridge Features: As seen in
FIG. 32, the toolkit tends to produce a bowing (see chin cleft,
reduced nasal aperture, and reduced orbital aperture) of some ridge
curves. This is primarily due to the interpolation process for
surface tile boundary points beside Type II landmarks.
[0238] 4. Time Taken: Skull surface extraction with the toolkit for
an experienced operator requires 2-3 hours. In the SASE interface,
the operator may extract a skull surface in less than one hour,
including the manual effort in correcting 10 to 30% of the
automated ridge curve extractions.
III. Averaging
[0239] 1. Parametric Surfaces and Surface Extractions
[0240] A. Surface Parameterization
[0241] We use a semi-automatic surface parameterization program,
Simulated Annealing Surface Extraction (SASE), to homology map the
segmented volume image surface. The SASE surface parameterization
occurs via superimposition of a ridge curve-based deformable
template to each sample member's surface of interest. The first
step in the SASE method is manual location of Type II landmarks on
a graphical manifold representation of the segmented surface
voxels. Second, the candidate ridge curve-based deformable template
is warped to these landmarks as in the toolkit. Finally, ridge and
tiling curve superimposition occurs via a simulated annealing
algorithm which primarily evaluates surface curvature information
found in the unsubsampled volume image data. The superimposed ridge
curve-based deformable template defines a homology mapped
parametric surface that is encoded as a series of B-spline space
curves. The toolkit does not use volume image data to superimpose
ridge curve-based deformable templates, relying instead on the
operator to manually nudge the ridge curve deformable template onto
a graphical manifold.
[0242] B. B-Spline Ridge and Tiling Curve Wireframe
[0243] A B-spline space curve may be recorded as a set of 3D
{Q.sub.k},(k=1, . . . , n) points incrementally stacked along a
sweep in 3D space represented as a pth degree non-rational
B-spline. By assigning a parameter value, u.sub.k, to each
{Q.sub.k}, with an appropriate knot vector U={a, . . . a,
u.sub.p+1, . . . ,u.sub.m-p-1, b, . . . , b}, where a and b are
initial and final conditions (i.e., 0 . . . 1), and m is an
ascending function, a B-spline space curve is:
Q k = C ( u k ) = i = 0 n N i , p ( u k ) P i ( EQ 32 )
##EQU00019##
where N.sub.i,p(u.sub.k) is the pth degree B-spline basis function.
These basis functions define the non-periodic and non-uniform knot
vector, U, and control points, P.sub.1. We set the degree of
interpolation components as 3, and order 4 (degree+1) for all space
curves. B-spline basis functions of the 3rd degree are defined
as:
N i , 0 ( u ) = { 1 , ... , u , < u < u i + 1 0 , ... , ( u
.ltoreq. u , ) or ( u .gtoreq. u i + 1 ) ( EQ 33 ) N i , 3 ( u ) =
u - u i u i + 3 - u i N i , 2 ( u ) + u i + 4 - u u i + 4 - u i + 1
N i + 1 , 2 ( u ) ( EQ 34 ) ##EQU00020##
from the initial condition, these interpolant variables are
computed with the knot vectors.
[0244] We use a chord length method to select each knot vector and
a linear averaging method to determine knot spacing along the ridge
curve and tiling curve B-spline space curves. Knot density is
pre-determined by the operator. Since each axis is interpolated
separately, as in the 2 dimensional case, the x, y, and z axes can
be treated separately to encode each space curve.
[0245] C. B-Spline Surface Tile Grids
[0246] Based on ridge curve and geodesic tile (tiling curve)
boundaries alone, a surface tile is defined as a bi-linearly
blended Coons surface. We extend this definition. We take the
internal tile point grid and tie it to the surface voxels segmented
from the original volume. Next, we encode these points as row and
column B-spline curves in u, v space (FIG. 36A). This improved
homology mapping prevents overlapping of internal tile points
during averaging. Our surface tiles are either three or four sided.
In the case of a three sided tile, the fourth side is generated by
placing the same number of second side points, with each point
filled by a common degenerate value from the intersection of the
first and third boundary space curves (FIG. 36B).
[0247] Our template uses a bi-linear blending function to define
location for interior tile points:
Q ( u , w ) = [ 1 - u u 1 ] [ - P ( 0 , 0 ) - P ( 0 , 1 ) P ( 0 , w
) - P ( 1 , 0 ) - P ( 1 , 1 ) P ( 1 , w ) P ( u , 0 ) P ( u , 1 ) 0
] [ 1 - w w 1 ] ( EQ 35 ) ##EQU00021##
where the functions ((1-u), u, (1-w), w) (i.e., blending functions)
blend the boundary space curves to produce the preliminary internal
shape of the surface as a cosine function; P is the parametric
function representing the interior tile points with respect to
parametric values u and w. The surface points generated are fitted
to the original image surface via a surface normal trace to the
surface voxels segmented from the volume image. The trace modifies
the shape of the tile by registering each final surface tile point
to the original volume image surface (i.e., segmented voxels). The
operator sees the result projected onto a graphical manifold
surface. The surface tile points are recorded as rows of B-spline
space curves spanning the two pairs of tile bounding space curves
(FIG. 36).
[0248] The caption to FIG. 37 lists the Type II landmarks found in
both soft tissue face and boney skull surface templates. Our
current soft tissue face template consists of 56 landmarks, 103
space curves (64 ridge curves and 49 geodesics), and 48 surface
tiles (FIG. 37A). Our current skull template consists of 54
landmarks, 104 space curves (40 ridge curves and 64 geodesics), and
47 surface tiles.
[0249] 2. Computational Tools for Average Generation
[0250] As with the toolkit, the first step in our surface averaging
algorithm is the generation of an average Type II landmark
configuration. We use the Procrustes superimposition method to
produce an average Type II landmark configuration. The ridge and
tiling space curve B-splines making up each sample member are
globally warped to this average landmark shape and locally
unwarped.
[0251] A. Procrustes Superimposition for Production of Average Type
II Landmark Configuration
[0252] The Procrustes superimposition method is applied to two
landmark coordinate configurations at a time. The two homologous
shapes are referred to as landmark configurations X.sub.1 and
X.sub.2, are pxk matrices with k-dimensional coordinates of the p
vertices of Type II landmarks. Since the geometric operations are
linear, this could be extended to include any number of
configurations, from X.sub.2 to X.sub.N. The order in which
Procrustes Fitting occurs has no effect on the shape of the final
configuration. A Procrustes fit attempts to minimize the sum of the
squared distances between corresponding vertices of both
configurations.
[0253] Translation: The coordinates of X.sub.1 and X.sub.2 are
translated to a shared origin, thereby k dimensions are reduced
from the overall pk dimensions. This is achieved by subtracting the
respective means of each axis from the X.sub.1 and X.sub.2
matrices. X.sub.1 and X.sub.2 are translated to the origin by
pre-multiplying X by (I-P), where P is a p.times.p identity matrix,
and is a matrix with every element equal to 1/p. Therefore I-P
is:
[ 1 - 1 p - 1 p - 1 p - 1 p - 1 p 1 - 1 p - 1 p - 1 p - 1 p - 1 p 1
- 1 p - 1 p - 1 p - 1 p - 1 p 1 - 1 p ] ( EQ 36 ) ##EQU00022##
and the translation of X.sub.1 and X.sub.2 to a shared origin is
computed as X.sub.1'=(1-P)X.sub.1 and X.sub.2'=(I-P)X.sub.2. This
translation is shown in FIG. 38A.
[0254] Scaling: A scalar value, called "Centroid Size," S.sub.1 of
X.sub.1, is calculated by taking the square root of the sum of
squared distances of each landmark to the centroid. Centroid size
is computed for each axis separately and divided for each control
point P:
S.sub.1= {square root over
(trace((I-P)X.sub.1X.sub.1'(I-P)))}{square root over
(trace((I-P)X.sub.1X.sub.1'(I-P)))}(EQ 37)
[0255] Similarly, S.sub.2 of X.sub.2 is calculated. The translated
and scaled versions
X 1 '' = 1 S 1 X 1 ' and X 2 '' = 1 S 2 X 2 ' ##EQU00023##
are seen in FIG. 38B.
[0256] Rotation: The translated and scaled configuration is next
rotated to obtain a minimum root sum of square distance between
corresponding points. The first step is to rotate all
configurations about the first configuration. Next the "rough"
average (FIG. 38) position of all the rotated configurations is
determined. Finally, all configurations are rotated to the rough
average position to obtain the final average configuration. The
configuration X.sub.1'' is rotated to overlay on configuration
X.sub.2'' so that the root sum of squared distances between the
corresponding landmarks is minimized. First the X.sub.2''
configuration is rotated by transposing it and pre-multiplying it
by a H, a k.times.k matrix. Therefore, the X.sub.2'' configuration
to rotated X.sub.1'' is
{circumflex over (X)}.sub.2'=X.sub.2''H (EQ 38)
Where {circumflex over (X)}'.sub.2 is rotated X.sub.2''
configuration, t is the transpose operator and H is obtained by
H=VSU' where V and U are obtained from the singular value
decomposition of the product of transpose of X.sub.1'' and
X.sub.2'' as:
X.sub.1''X.sub.2''=U.SIGMA.V' (EQ 39)
S is the matrix of diagonal elements of .+-.1 corresponding to the
signs of the .SIGMA. matrix in equation 39. The average
configuration X.sub.A is the average of X.sub.1'',{circumflex over
(X)}.sub.2. The rotation is shown in FIG. 38C. For the multiple
configuration average, we first compute the rough average X'.sub.A
by rotating the configurations to the first configuration, one at a
time:
X.sub.2'={X.sub.1'',X.sub.2''},X.sub.3'={X.sub.1'',X.sub.3''} . . .
X'.sub.N={X.sub.1'',X.sub.N''} (EQ 40)
where N is the number of configurations to be superimposed.
X'.sub.A is computed as:
X A = 1 N { X 1 '' + X 2 ' ^ + ... + X N ' ^ } ( EQ 41 )
##EQU00024##
[0257] To produce the final Type II landmark and tile boundary
B-spline control point average, we rotate each configuration to the
rough average configuration, X.sub.A, one at a time:
{circumflex over (X)}.sub.1={X.sub.A',X.sub.1''},{circumflex over
(X)}.sub.2={X.sub.A',X.sub.2''} . . . {circumflex over
(X)}.sub.N={X.sub.A',X.sub.N''} (EQ 42)
The consensus average landmark configuration X.sub.A is computed
from {circumflex over (X)}.sub.1, {circumflex over (X)}.sub.2, . .
. , {circumflex over (X)}.sub.N as:
X A = 1 N { X ^ 1 + X ^ 2 + ... + X ^ N } ( EQ 43 )
##EQU00025##
[0258] B. Thin Plate Splines
[0259] We use a thin plate spline to initially position all sample
member surfaces at the average Type II landmark configurations. The
first landmark transformations matched homologous regions of two
and three dimensional wireframe drawings of organisms. A spatial
interpolation routine, widely used in geology and metallurgy, the
thin plate spline, was suitably modified and applied to this
problem. A thin plate spline warp results in a bending energy
matrix describing the shape deformation of the anchor points (i.e.,
control points and Type II landmarks). We use the thin plate spline
to provide a global warp of sample member's ridge, tiling curve,
and internal tile 13-spline points to the rough average Type II
landmark configuration.
[0260] Although the warp is based only on the Type II landmarks
shared by all sample members, the appropriate position of all
surface tile points is achieved via the mapping function
f(R.sup.2-R.sup.2) such that f(S.sub.1)=S.sub.A, based on the Type
II landmarks. This function, a bi-harmonic equation U(r)=r.sup.2
log r.sup.2, is a generalized solution. The basis function, r=
{square root over (x.sup.2+y.sup.2)} is limited to the two
dimensional case. We extend it to three dimensional surfaces as
U(r)=r, where r= {square root over (x.sup.2+y.sup.2+z.sup.2)}. A 3D
thin plate spline function can be generalized and interpreted
visually as
{(r,s,i,.delta.(r,s,t)).epsilon.R.sup.4.parallel.r,s,t.epsilon.R}.
If S.sub.A is our consensus average Type II landmark reference
configuration, with k=3 dimensions and p vertices, (i.e., a
k.times.p matrix), the matrix S.sub.A and function matrix K are
defined as:
S A i = [ 1 1 1 ... 1 x 1 x 2 x 3 ... x p y 1 y 2 y 3 ... y p z 1 z
2 z 3 ... z p ] , 4 .times. p ( EQ 44 ) K = [ 0 U ( r 12 ) ... u (
r 1 p ) U ( r 21 ) 0 ... U ( r 2 p ) ... ... ... ... U ( r p 1 ) U
( r p 2 ) ... 0 ] , p .times. p ( EQ 45 ) ##EQU00026##
[0261] U(r.sub.i,j) is defined in 3D as the distance from the
i.sup.th to j.sup.th the point in the reference configuration
as:
U(r.sub.i,j)+ {square root over
((x.sub.i-x.sub.j)+(y.sub.i-y.sub.j)+(z.sub.i-z.sub.j))}{square
root over
((x.sub.i-x.sub.j)+(y.sub.i-y.sub.j)+(z.sub.i-z.sub.j))}{square
root over ((x.sub.i-x.sub.j)+(y.sub.i-y.sub.j)+(z.sub.i-z.sub.j))}
(EQ 46)
[0262] Next, a projection operator L is defined as:
L = [ K S A S A i O ] , ( p + 4 ) .times. ( p + 4 ) ( EQ 47 )
##EQU00027##
O is a 4.times.4 null matrix. We use these mapping functions to
globally warp a series of ridge curve based deformable template
(i.e., ridge, tiling, and internal tile B-spline space curves)
wireframes to the rough average Type II landmark configuration via
a global warp. Once in position, each space curve is locally
unwarped to begin the averaging process.
[0263] In general, thin plate spline warping is defined with a set
of homologous points Y, as L.sup.-1Y, and the transformation is
defined:
S ( X ni ' ) = A [ 1 x y z ] ' ni + j P n j U ( r j ) , n = 1 ... N
, i = 1 ... I ( EQ 48 ) ##EQU00028##
[0264] 3. Average Surface Generation
[0265] Our space curve averaging is a two step process. The first
step begins following the local unwarping of space curves that were
globally thin plate spline warped to an average landmark
configuration. Then all sample members are locally unwarped and
averaged, B-spline sample by B-spline sample (i.e., each set of
homologous ridge, tiling, or internal tile curves).
[0266] A. Average Type II Landmark Configuration
[0267] We use the Procrustes Superimposition method to align the
sample member landmark configurations and compute an average
configuration from them. Procrustes fitting removes variation due
to position, orientation, and scale prior to our calculating a
simple average. In equation 45 we see that the computed average
configuration is set to X.sub.A.
[0268] B. Average Space Curve Generation
[0269] This is accomplished algorithmically as follows: Let the
configurations S.sub.n where n varies from 1 to N wireframe space
curves (i.e., The wireframe consists of the tile boundary B-spline
space curves). X.sub.A is an average landmark configuration
(Equation 45). Each sample configuration S.sub.n is warped to
X.sub.A landmarks. We refer to this sample of globally warped
configurations as S.sub.n. Next, each member of S.sub.n is locally
unwarped to the original configuration S.sub.n for each space
curve. The locally unwarped and assembled space curves are shown in
FIG. 39 drawn between hypothetical Type II landmarks A and B. A
"rough" average curve is produced as a direct average of the
control point coordinates (note: this is possible for B-splines but
not for Bezier splines where the control points may not lie on the
encoded space curve). The rough average curve is then re-splined.
The re-splined rough average curve is displayed as a dashed line in
FIG. 39. The second step is the generation of final average space
curves. Sampling planes (perpendicular bi-sectors) are positioned
at regularly spaced intervals along the rough average curve. A
control point for each sampling plane is produced from the average
of intercepted points. Finally, a B-spline space curve is generated
from the average intercept points (FIG. 39, thick dark line).
[0270] The space curve averaging is done separately for each type
of space curve (i.e., ridge, tiling, and internal tile curves).
Ridge curves are further divided into arc-segments between Type II
landmarks for averaging. The complete set of average space curves
represents the average surface.
[0271] C. Internal Tile Curve Averaging
[0272] Internal tile curves have either u or v orientation (see
FIG. 36). Therefore each internal tile point belongs to two
B-splines. We compute the average u and v internal tile B-spline
space curves separately. The average of these two points is then
substituted back to internal tile vertices, in order to prepare a
triangulated (graphical manifold) surface. Algorithmically,
B-spline internal tile curve segments in the u and v direction may
be represented as matrices:
S i ( C u ) = [ C 1 , 1 C 1 , 2 C 1 , v C 2 , 1 C 2 , 2 C 2 , v C u
, 1 C u , 2 C u , v ] , S i ( C v ) = [ C 1 , 1 C 1 , 2 C 1 , u C 2
, 1 C 2 , 2 C 2 , u C v , 1 C v , 2 C v , u ] ( EQ 49 )
##EQU00029##
The computed interior surface tile control points thus are the
average of the above two matrices, S.sub.t(C.sub.u) and
S.sub.t(C.sub.v):
S i ( C ) = S i ( C u ) + S i ( C v ) i = [ C 1 , 1 C 1 , 2 C 1 , v
C 2 , 1 C 2 , 2 C 2 , v C v , 1 C v , 2 C u , v ] = [ C 1 , 1 C 1 ,
2 C 1 , u C 2 , 1 C 2 , 2 C 2 , u C v , 1 C v , 2 C v , u ] ( EQ 50
) ##EQU00030##
where t varies from 1 to the number of tiles in the whole sample
surface.
[0273] 4. Precision and Inter-Method Tests and Results
[0274] The SSA program was written in C in the Silicon Graphics
Unix (IRIX) environment. It was tested on a SGI Octane workstation
(R10000 CPU, 256 M byte RAM). The program interface integrates
Procrustes superimposition, thin plate spline warp, and surface
averaging methods with graphical manifold surface rendering
operator feed back. User interaction occurs via X11/Motif
protocols.
[0275] To test the SSA program, we averaged soft tissue face
surfaces (FIG. 40A) and skull surfaces (FIG. 40B) from segmented
from 3D CT volume images. These five 3D CT data sets (B2537, B2621,
B3037, B3095, B3195) were provided by a recall of subjects. All
subjects are above 60 years in age, female, have dental fillings,
and often have large oral prostheses (dentures). The soft tissue
face and skull surface segmentation of the volume images was
performed in a separate program, SOFM. The ridge curve-based
deformable template superimposition and surface parameterization
was done in an another program, SASE.
[0276] A. Intra-Operator Precision Test
[0277] We conducted a precision test of the SSA methods between two
sessions by one operator. The operator introduces variation through
initial seeding and subsequent fine tuning of the ridge and tiling
curve template superimposition. The overall root mean square error
between the surface tile points of two averaging sessions, of the
boney skull surface average was 0.0 mm and for the soft tissue face
surface average was 0.0 mm. This is complete agreement between
sessions. This intra-operator reproducibility result (i.e.,
precision) provides evidence that the SSA method is reliable.
[0278] B. Inter-Method Tests
[0279] We next compared the surface sampling used to produce
average surfaces by the SSA program versus that of the toolkit
using the same segmented voxel data obtained from the SOFM program.
We hypothesize that improved homology of tiling and internal tile
curves in the SSA method will reduce the variance of the sample
members about the average.
[0280] i. Visual Comparison
[0281] We parameterized all 10 surfaces with both the SASE and the
toolkit for the soft tissue face (FIG. 41A) and the boney skull
(FIG. 41B). Note that the toolkit sample and average surface images
appear to bow along every ridge curve edge (e.g., orbits and nasal
regions). Under the chin, we note that this problem has exaggerated
a notch not seen in the SSA average. This bulging also results in
smaller orbital and nasal apertures relative to the shape seen in
the original segmented surface voxels. This and the smoother
surface may be an artifact of splining a lesser number of points
than in the SSA method.
[0282] ii. Inter-Method Comparison of Average Surface Warped to
Members' Original Segmented Voxel Surface
[0283] We thin plate spline warped the average surface to the Type
II landmark configuration of it's member surfaces. Then, we traced
each fitted average surface tile point to the segmented surface
image. Both averaging methods (i.e., toolkit and SSA) used surfaces
extracted from the same segmented volume image, providing a common
basis for comparison between the two methods. The computed square
root of sum of squared differences are listed in Table XII. Note,
the number of tile points in the toolkit surface extractions is
well below that used in the SSA method. Given this difference, we
observed an average of 6 mm distance for the toolkit averaging
method for all points in the soft tissue face and skull surfaces,
whereas the average SSA distance were 0.448 mm in the case of the
soft tissue average and 0.419 mm in the case of the skull average.
Table XIII presents these results tile by tile. FIG. 42A displays
color coded difference images of toolkit skull and soft tissue face
surface averages with their source segmented voxel surfaces.
Similarly FIG. 42B shows the same color coded difference images of
the SSA skull and soft tissue face surface averages. The same range
of difference is scaled from lowest difference in blue to largest
difference in red in both images. Note that despite higher sampling
density the internal tile points are as true as the tile borders in
the SSA extractions.
[0284] iii. Results
[0285] In FIG. 43A, we see the toolkit average surface prior to the
final step of normal interpolation and averaging. FIG. 43B presents
the same surfaces afterward. FIG. 43C presents the final SSA
average surface; the SSA algorithm does not average surface
normals. Note that seams are present but they less disturb the
continuity of surface curvature between the tile surfaces. In both
cases these are graphical rendering artifacts, not an indication of
inaccuracy in surface parameterization or a non-shape preserving
average.
[0286] 5. Conclusion
[0287] One reason that the SSA average surface results better
preserve surface curvature continuity, and reduce variance about
the average, is the higher sampling density (FIG. 44). The
toolkit's graphical manifold is subsampled at four stages. Surface
rendering of the segmented voxels introduces an error of .+-.1 to
.+-.3 mm assuming a sub-sampling rate of 2.times.2.times.2.
[0288] Second, a wireframe is manually nudged into best position on
this graphical manifold and then interpolated. This step increases
the error at each point by .+-.1 to .+-.3 mm depending on the
variable level of parametric smoothing that the program
automatically applies. Third, the surface is sparsely encoded. This
may eliminate some otherwise useful information from the patient
surface image.
[0289] Finally when generating the average, only the surface tile
boundaries (i.e., ridge curve and geodesic lines) are truly
averaged in the toolkit. The average interior tile points are
assigned, not found. Overall, the errors accumulated at each
surface tile point vary from 3 to 6 mm.
[0290] It appears to us that the SSA method improves tiling and
internal tile curve homology assignment. The SSA method better
preserves overall shape of 3D surfaces during averaging because of
(1) the improved homology assignment of internal tile points, and
(2) the extension of the original space curve averaging method to
the entire surface.
[0291] The utility of average surface images for clinical diagnosis
needs validation. Their use for boney prosthetic design is
apparent. These methods were originally developed to quantify
diagnostic morphometric differences of clinical populations,
however we expect these data have applications at every stage of
patient care, as well as other deformable model applications
including animation.
[0292] The invention has been described with reference to the
preferred embodiment. Obviously, modifications and alterations will
occur to others upon reading and understanding the preceding
detailed description. It is intended that the invention be
construed as including all such modifications and alterations
insofar as they come within the scope of the appended claims or the
equivalents thereof.
* * * * *