|
Medical Patent Abstract
A method for determining the location, shape and orientation of
a tumor in a medical image includes finding a plurality of spatial
extrema .mu. of a D-dimensional spatial signal f for a set of bandwidths
H by performing mean shift-based gradient-ascent iterations for
a set of bandwidths H and then determining a D-dimensional spread
and orientation of the signal about each extrema .mu. by estimating
a covariance .SIGMA. of the signal f for each extrema .mu.. The
optimal estimate of .mu. and .SIGMA. is determined by performing
a Jensen-Shannon divergence on the full set of .mu. and .SIGMA..
Medical Patent Claims
What is claimed is:
1. A method for analyzing three-dimensional structures in medical
images, in order to determine the location, shape and orientation
of a tumor, said method comprising the steps of: finding at least
one spatial extrema .mu. of a D-dimensional spatial signal f for
a set of bandwidths H; estimating a D-dimensional spread and orientation
of the signal about each extrema .mu. by estimating a covariance
.SIGMA. of the signal f for each extrema .mu.; and finding a most
stable estimate of .mu. and .SIGMA. from the at least one extrema
.mu. and covariances .SIGMA. to find from the signal f an optimal
estimate of a target structure wherein the spatial extrema .mu.
of signal f are found by solving an equation defined by .function..ident..intg..mu..times..times..PHI..function..mu..times..fu-
nction..mu..times.d.mu..intg..PHI..function..mu..times..function..mu..time-
s.d.mu. ##EQU00006## where m(x; H) is an extended fixed-bandwidth
mean shift vector, x is a spatial location corresponding to a signal
measurement, and .PHI.(x-u;H) is a Gaussian kernel for a bandwidth
H.
2. The method of claim 1, wherein .PHI.(x-.mu.; H) can be defined
as .function..times..function..mu. ##EQU00007## wherein D.sup.2(x,.mu.;
H)=(x-.mu.).sup.T H.sup.-1(x-.mu.) and H.sup.-1 is a weighted harmonic
mean of the bandwidth matrices, .function..times..function..times.
##EQU00008## and wherein the weights can be defined as w.sub.i .function..times..function..times..function..mu..times..times..function..-
times..function..mu. ##EQU00009## and can be normalized to unity.
3. The method of claim 1, wherein finding the spatial extrema .mu.
of signal f comprises the steps of: making an estimate of an extrema
and evaluating the extended fixed-bandwidth mean shift vector for
this extrema, wherein y.sub.1 is used to denote the first term of
m(x; H) for the initial estimate of .mu..sub.1; evaluating the extended
fixed-bandwidth mean shift vector m(x; H), replacing the second
term with y.sub.1 and replacing the initial estimate of .mu..sub.1
with the previous evaluation of m(x; H); repeating said evaluations
of the extended fixed-bandwidth mean shift vector m(x; H), each
time each time replacing the second term with the first term from
the previous iteration, and evaluating the first term on the previous
evaluation of m (x; H), until a value of .mu..sub.k is found for
which the extended mean shift vector m (x; H) is sufficiently close
to zero; and partitioning data points of the signal by grouping
said data points that converge into the same extrema.
4. The method of claim 1, wherein the step of estimating the covariance
comprises the steps of defining the covariance by A.SIGMA.=B where
A=(m.sub.1; . . . ; m.sub.t.sub.u).sup.T H.sup.-T is a D.times.D
dimensional positive definite matrix and B=(b.sub.i; . . . ; b.sub.t.sub.u).sup.T,
with .SIGMA.H.sup.-1m.sub.i.apprxeq.b.sub.i.ident..mu.-x.sub.i-m.sub.i;
and evaluating the equation .SIGMA.*=U.sub.p.SIGMA..sub.P.sup.-1U.sub.{tilde
over (Q)}.SIGMA..sub.{tilde over (Q)}U.sub.{tilde over (Q)}.sup.T.SIGMA..sub.P.sup.-1U.sub.P
.sup.T; where A.sup.TA=U.sub.P.SIGMA..sub.P.sup.2U.sub.P.sup.T and
B.sup.TB.ident.Q with {tilde over (Q)}=.SIGMA..sub.PU.sub.P.sup.TQU.sub.P.SIGMA..sub.P=U.sub.{tilde
over (Q)}.SIGMA..sub.{tilde over (Q)}.sup.2U.sub.{tilde over (Q)}.sup.T
are symmetric Schur decompositions.
5. The method of claim 1, wherein the most stable estimate of .mu.
and .SIGMA. is found by calculating, for a neighborhood a about
each bandwidth value h, a Jensen-Shannon divergence defined by .function..times..times..times..times..times..times..SIGMA..times..times.-
.SIGMA..times..times..times..mu..mu..times..times..SIGMA..times..mu..mu.
##EQU00010## .times..times..mu..times..times..times..mu..times.
##EQU00010.2##
6. A method for analyzing three-dimensional structures in medical
images, in order to determine the shape and orientation of a tumor
whose location has been marked, said method comprising the steps
of: sampling a set of starting points from a neighborhood of a marker;
estimating an extrema .mu. by performing mean shift-based gradient-ascent
iterations for a set of bandwidths H; performing a regular sampling
of points in a neighborhood of the extrema estimate .mu.; estimating
a spread and orientation about the extrema .mu. by estimating a
covariance .SIGMA. of the sampling of points of the extrema .mu.;
and estimating the stability of .mu. and .SIGMA. by calculating
a Jensen-Shannon divergence.
7. A method for determining the location, shape and orientation
of a tumor in a medical image, said method comprising the steps
of: defining an extended fixed-bandwidth mean shift vector m(x;
H) by an equation defined by .function..ident..intg..times..times..PHI..function..times..function..tim-
es.d.intg..PHI..function..times..function..times.d ##EQU00011##
wherein x is a spatial location corresponding to a signal measurement,
and wherein .PHI.(x-.mu.; H) can be defined as .function..times..function..mu.
##EQU00012## where D.sup.2(x,.mu.; H)=(x-.mu.).sup.T H.sup.-1(x-.mu.)
and H.sup.-1 is a weighted harmonic mean of the bandwidth matrices,
.function..times..function..times. ##EQU00013## and where the weights
can be defined as .function..times..function..times..function..mu..times..times..function..-
times..function..mu. ##EQU00014## and can be normalized to unity
.PHI.(x-.mu.; H); making an estimate of an extrema and evaluating
the extended fixed-bandwidth mean shift vector for this extrema,
wherein y.sub.1 is used to denote the first term of m(x; H) for
the initial estimate of .mu..sub.1; evaluating the extended fixed-bandwidth
mean shift vector m(x; H), replacing the second term with y.sub.1
and replacing the initial estimate of .mu..sub.1 with the previous
evaluation of m(x; H); repeating said evaluations of the extended
fixed-bandwidth mean shift vector m(x; H), each time each time replacing
the second term with the first term from the previous iteration,
and evaluating the first term on the previous evaluation of m (x;
H), until a value of .mu..sub.k is found for which the extended
mean shift vector m (x; H) is sufficiently close to zero; partitioning
data points of the signal by grouping said data points that converge
into the same extrema; estimating the D-dimensional spread and orientation
of the signal by defining the covariance as A.SIGMA.=B where A=(m.sub.1;
. . . ; m.sub.t.sub.u).sup.T H.sup.-T is a D.times.D dimensional
positive definite matrix and B=(b.sub.i; . . . ; b.sub.t.sub.u).sup.T,
with .SIGMA.H.sup.-1m.sub.i.apprxeq.b.sub.i.ident..mu.-x.sub.i-m.sub.i
and evaluating the equation .SIGMA..times..SIGMA..times..times..SIGMA..times..times..SIGMA..times.
##EQU00015## .times..times..SIGMA..times..times..times..times..times..times..ident..ti-
mes..times. ##EQU00015.2## .SIGMA..times..times..times..SIGMA..times..SIGMA..times..times..times.
##EQU00015.3## symmetric Schur decompositions; and finding a most
stable estimate of .mu. and .SIGMA. from the at least one extrema
.mu. and covariances .SIGMA. by calculating, for a neighborhood
a about each bandwidth value h, a Jensen-Shannon divergence defined
by .function..times..times..times..times..times..times..SIGMA..times..times.-
.SIGMA..times..times..times..mu..mu..times..times..SIGMA..times..mu..mu.
##EQU00016## where .mu..times..times..times..mu. ##EQU00017## and
finding an extrema of said divergences, to find from the signal
f an optimal estimate of a target structure.
8. A program storage device readable by a computer, tangibly embodying
a program of instructions executable by the computer to perform
the method steps for analyzing three-dimensional structures in medical
images, in order to determine the location, shape and orientation
of a tumor, said method steps comprising: finding at least one spatial
extrema .mu. of a D-dimensional spatial signal f for a set of bandwidths
H; estimating a D-dimensional spread and orientation of the signal
about each extrema .mu. by estimating a covariance .SIGMA. of the
signal f for each extrema .mu.; and finding a most stable estimate
of .mu. and .SIGMA. from the at least one extrema .mu. and covariances
.SIGMA. to find from the signal f an optimal estimate of a target
structure wherein the spatial extrema .mu. of signal f are found
by solving an equation defined by .function..ident..intg..times..times..PHI..function..times..function..tim-
es.d.intg..PHI..function..times..function..times.d ##EQU00018##
where m(x; H) is an extended fixed-bandwidth mean shift vector,
x is a spatial location corresponding to a signal measurement, and
.PHI.(x-u; H) is a Gaussian kernel for a bandwidth H.
9. The computer readable program storage device of claim 8, wherein
.PHI.(x-.mu.; H) can be defined as .function..times..function..mu.
##EQU00019## wherein D.sup.2(x,.mu.; H)=(x-.mu.).sup.T H.sup.-1
is a weighted harmonic mean of the bandwidth matrices, .function..times..function..times..function..mu..times..times..function..-
times..function..mu. ##EQU00020## and wherein the weights can be
defined as .function..times..times..function..times..function..mu..times..times..-
times..function..times..function..mu. ##EQU00021## and can be normalized
to unity.
10. The computer readable program storage device of claim 8, the
method steps for finding the spatial extrema .mu. of signal f comprises:
making an estimate of an extrema and evaluating the extended fixed-bandwidth
mean shift vector for this extrema, wherein y.sub.1 is used to denote
the first term of m(x; H) for the initial estimate of .mu..sub.1;
evaluating the extended fixed-bandwidth mean shift vector m(x; H),
replacing the second term with y.sub.1 and replacing the initial
estimate of .mu..sub.1 with the previous evaluation of m(x; H);
repeating said evaluations of the extended fixed-bandwidth mean
shift vector m(x; H), each time each time replacing the second term
with the first term from the previous iteration, and evaluating
the first term on the previous evaluation of m (x; H), until a value
of .mu..sub.k is found for which the extended mean shift vector
m (x; H) is sufficiently close to zero; and partitioning data points
of the signal by grouping said data points that converge into the
same extrema.
11. The computer readable program storage device of claim 8, wherein
the method steps for estimating the covariance comprises defining
the covariance by A.SIGMA.=B where A=(m.sub.1; . . . ; m.sub.t.sub.u).sup.T
H.sup.-T is a D.times.D dimensional positive definite matrix and
B=(b.sub.i; . . . ; b.sub.t.sub.u).sup.T, with .SIGMA.H.sup.-1m.sub.i.apprxeq.b.sub.i.ident..mu.-x.sub.i-m.sub.i;
and evaluating the equation .SIGMA..times..SIGMA..times..times..SIGMA..times..times..SIGMA..times..ti-
mes..times..SIGMA..times..times..times..times..times..times..ident..times.-
.times..times..times..SIGMA..times..times..times..times..times..SIGMA..tim-
es..SIGMA..times. ##EQU00022## are symmetric Schur decompositions.
12. The computer readable program storage device of claim 8, wherein
the method steps for finding the most stable estimate of .mu. and
.SIGMA. comprise calculating, for a neighborhood a about each bandwidth
value h, a Jensen-Shannon divergence defined by .function..times..times..times..times..times..times..SIGMA..times..times.-
.times..SIGMA..times..times..times..mu..mu..times..times..times..SIGMA..ti-
mes..mu..mu. ##EQU00023## where .mu..times..times..times..mu. ##EQU00024##
and finding an extrema of said divergences.
Medical Patent Description
BACKGROUND OF THE INVENTION
A problem in the volumetric medical image analysis is to characterize
the 3D local structure of tumors across various scales, because
the size and shape of tumors varies largely in practice. Such underlying
scales of tumors also provide useful clinical information, correlating
highly with probability of malignancy. There are a number of previously
proposed approaches addressing this problem. However, these prior
art approaches are prone to be sensitive to signal noises and their
accuracy degrades when the target shapes differ largely from an
isolated Gaussian. In the medical domain, these constraints are
too strong since many tumors appear as irregular shapes within noisy
background signals.
SUMMARY OF THE INVENTION
This invention is directed to methods for the robust estimation
of the covariance matrix that describes the spread and 3-dimensional
(3D) orientation of the structure of interest. Given an input signal,
a mean shift-based gradient ascent is performed for an extended
mean-shift vector using all the available data points for each analysis
bandwidth. The data points that converge to the same point are grouped
together, forming a set of local structure candidates. These convergence
candidates can be interpreted as spatial extrema of the signal.
Then, for each candidate, the underlying scale is determined by
estimating a covariance matrix by a constrained least-squares method.
Finally, for each candidate, a stability test is performed across
the analysis scales, resulting in an optimal scale estimate for
each local target. As a result, one can find a signal's local scales
that can vary spatially.
These methods utilize an algorithm, referred to herein as a mean
shift-based bandwidth selection algorithm, for analyzing general
discretized signals and estimating fully parameterized covariance
matrices. With the methods of the invention, it becomes possible
to address the problem of representing local structures in images.
The robust mode and scale selection of the mean shift-based bandwidth
selection algorithm together with the consideration of the fully
parameterized covariance matrix also enables one to estimate a tumor's
scale in more flexible and robust manner, mitigating the aforementioned
shortcomings of the previous methods. Since many target objects
in the medical domain possess complex 3D structures, the methods
of the invention can be deployed for a number of different application
scenarios.
BRIEF DESCRIPTION OF THE DRAWINGS
FIG. 1 is a flow chart depicting a method of a preferred embodiment
of the invention.
FIG. 2 is a flow chart depicting a method of another preferred
embodiment of the invention.
FIGS. 3-4 are images depicting results of analyses performed using
a preferred method of the invention.
FIG. 5 depicts a preferred embodiment of a computer system for
implementing a preferred method of the invention.
DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT
1. Fixed-bandwidth Mean Shift Vector for Continuous Signals
A medical image can be represented by a D-dimensional continuous
signal f:R.sup.D.fwdarw.R evaluated at n d-dimensional points x.sub.i,
and the uncertainty associated with each point x.sub.i can be represented
by a D.times.D matrix H.sub.i, for an i.di-elect cons.1, . . . n.
The matrices H.sub.i are referred to as bandwidth matrices. The
signal can have one or more extrema. An extrema of the signal can
be associated with a location of a tumor or other target object.
Referring now to FIG. 1, a first step in the analysis of the medical
image represented by f is to determine, at step 11, the spatial
extrema .mu. of the signal. To find the extrema, one can first define
a function, m (x; H), where x is a spatial location corresponding
to a signal measurement and H is the corresponding bandwidth, referred
to herein as an extended mean shift vector, by
.function..ident..intg..mu..times..times..PHI..function..mu..times..functi-
on..mu..times.d.mu..intg..PHI..function..mu..times..function..mu..times.d.-
mu. ##EQU00001## where a Gaussian kernel .PHI.(x-.mu.; H) can be
defined as
.function..times..function..mu. ##EQU00002## with D.sup.2(x,.mu.;
H)=(x-.mu.).sup.T H.sup.-1(x-.mu.) and H.sup.-1 is a weighted harmonic
mean of the bandwidth matrices,
.function..times..function..times. ##EQU00003## The weights can
be defined as
.function..times..function..times..function..mu..times..times..function..t-
imes..function..mu. ##EQU00004## and can be normalized to unity.
Eq. (3) can be used to locate spatial extrema .mu. of f given a
fixed analysis bandwidth H as follows. First, make an estimate of
an extrema, .mu..sub.1, and then evaluate m.sub.1 (x; H) for this
extrema from Eq. (3). If y.sub.1 is used to denote the result of
the first term of Eq. (3) for the initial estimate of .mu..sub.1,
then, for the next iteration of Eq. (3), x is replaced with y.sub.1
and .mu..sub.1 is replaced with m.sub.1 (x; H), denoted as .mu..sub.2.
This process can be repeated, each time replacing the second term
of Eq. (3) with the result of the first term from the previous iteration,
and evaluating the first term on the previous evaluation of m (x;
H). For each iteration k of Eq. (3), the resulting difference will
converge to zero. The value of .mu..sub.k for which the extended
mean shift vector m (x; H) is sufficiently close to zero can be
taken as an extrema of the signal f. The data space of the signal
can be partitioned by grouping data points that converge into the
same extrema.
2. Constrained Least-squares Solution of Covariance Matrices
The next step, step 12 of FIG. 1, is to estimate the D-dimensional
spread and orientation of the tumors whose center .mu. as a spatial
extremum was found in step 11. The geometrical information of a
D-dimensional local surface can be characterized by a covariance
matrix .SIGMA. estimated at the extrema.
The covariance .SIGMA. can be defined by the equation m(x).apprxeq.-H(.SIGMA.+H).sup.-1(x-.mu.)
when f can be approximated by a Gaussian. This can be rewritten
in the following simple form, .SIGMA.H.sup.-1m.sub.i.apprxeq.b.sub.i.ident..mu.-x.sub.i-m.sub.i
(4)
Considering all the trajectory points {x.sub.i: i=1, . . . , t.sub.u}
that converge to an extremum .mu., an over-complete set of linear
equations can be contructed: A.SIGMA..apprxeq.B, (5) where A=(m.sub.1;
. . . ; m.sub.t.sub.u).sup.TH.sup.-T, B=(b.sub.i; . . . ; b.sub.t.sub.u).sup.T,
and where .SIGMA. is a symmetric, positive definite matrix in R.sup.D.times.D.
The covariance can be estimated by a constrained least-squares solution
of Eq. (5). This solution yields the following closed form, .SIGMA.*=U.sub.p.SIGMA..sub.P.sup.-1U.sub.{tilde
over (Q)}.SIGMA..sub.{tilde over (Q)}U.sub.{tilde over (Q)}.sup.T.SIGMA..sub.P.sup.-1U.sub.P
.sup.T; (6) where the solution involves the following symmetric
Schur decompositions: A.sup.TA=U.sub.P.SIGMA..sub.P.sup.2U.sub.P.sup.T
and B.sup.TB.ident.Q with {tilde over (Q)}=.SIGMA..sub.PU.sub.P.sup.TQU.sub.P.SIGMA..sub.P=U.sub.{tilde
over (Q)}.SIGMA..sub.{tilde over (Q)}.sup.2U.sub.{tilde over (Q)}.sup.T.
This closed form can be found by determining the unique minimizer
for an area, g(Y).ident..parallel.AY-BY.sup.-T.parallel..sub.F.sup.2,
where .SIGMA.=YY.sup.T. 3. Scale Selection
The above two steps can result in pairs of center location and
covariance estimates {.mu..sub.h; .SIGMA..sub.h} for each analysis
bandwidth H. The next step, step 13, concerns finding the optimal
estimate of the target structures analyzed across a range of bandwidths.
This optimal estimate can be found by using a form of the Jensen-Shannon
divergence,
.function..times..times..times..times..times..times..SIGMA..times..times..-
SIGMA..times..times..times..mu..mu..times..times..SIGMA..times..mu..mu..ti-
mes..times..times..times..mu..times..times..times..mu. ##EQU00005##
Given a neighborhood parameter a, the divergence can be computed
for each analysis bandwidth h. The extremum of the divergences JS(h)
across the bandwidths can provide a final scale estimate that is
most stable over a range of scales.
The stability test described requires the set of analysis bandwidths
a priori. In one embodiment of the invention, H=hI and h is varied
with a constant step. In order to achieve higher performance for
the scale selection, it is preferred to have more densely distributed
analysis bandwidths. However, such dense sampling can prohibitively
enlarge the search space, especially when a fully parameterized
H is considered.
4. Algorithm for Local Multi-scale Analysis
In some application scenarios, the task to be solved is to represent
the scale of local structure whose rough location is provided by
another means. An example is the structural analysis of tumors whose
locations in a volumetric image are provided manually by radiologist.
The simplest strategy in such a case is to perform the mean shift
iteration only from the given marker point. The convergence point
serves as the tumor center estimate and all the trajectory points
from the marker are used to estimate the scale. This naive strategy
can fail when the provided locations are contaminated by uncertainties
and when the iteration converges too soon, forcing the Eq.(4) to
be under-complete and rank-deficient. These issues can be addressed
by the following steps, depicted in FIG. 2. First, at step 21, consider
a set of starting points sampled from the neighborhood of the marker.
At step 12, after performing mean shift iterations, the point to
which most starting points converged serves as a location estimate
.mu.. Next, at step 13, a regular sampling around the estimate .mu.
can be performed. The scale estimate .SIGMA. can be given by solving,
at step 14, Eq. (4) using all the trajectory points that converged
to .mu.. The same stability test of Eq. (7) can be used for the
final estimate at step 15.
6. EXAMPLES
A 3D domain implementation of the local multi-scale analysis algorithm
described in Section 5 is evaluated with high-resolution computerized
tomography (HRCT) images of 14 patients displaying pulmonary tumors.
A total of 44 analysis scales with 0.25 interval h=(0.25.sup.2;
. . . ; 11.sup.2) and a=1 were used. The rough location of the tumors
were provided. As a pre-process, volumes of interest of size 32.times.32.times.32
are extracted using the markers. FIGS. 3 and 4 show examples of
the resulting center and part-solid nodules whose geometrical shapes
are more deviated from the simple Gaussian structure. The correct
estimation of the tumor locations, spreads, and 3D orientations
for these difficult cases demonstrates the effectiveness of the
methods of the invention.
It is to be understood that the present invention can be implemented
in various forms of hardware, software, firmware, special purpose
processes, or a combination thereof. In one embodiment, the present
invention can be implemented in software as an application program
tangible embodied on a computer readable program storage device.
The application program can be uploaded to, and executed by, a machine
comprising any suitable architecture.
Referring now to FIG. 5, according to an embodiment of the present
invention, a computer system 101 for implementing the present invention
can comprise, inter alia, a central processing unit (CPU) 102, a
memory 103 and an input/output (I/O) interface 104. The computer
system 101 is generally coupled through the I/O interface 104 to
a display 105 and various input devices 106 such as a mouse and
a keyboard. The support circuits can include circuits such as cache,
power supplies, clock circuits, and a communication bus. The memory
103 can include random access memory (RAM), read only memory (ROM),
disk drive, tape drive, etc., or a combinations thereof. The present
invention can be implemented as a routine 107 that is stored in
memory 103 and executed by the CPU 102 to process the signal from
the signal source 108. As such, the computer system 101 is a general
purpose computer system that becomes a specific purpose computer
system when executing the routine 107 of the present invention.
The computer system 101 also includes an operating system and micro
instruction code. The various processes and functions described
herein can either be part of the micro instruction code or part
of the application program (or combination thereof) which is executed
via the operating system. In addition, various other peripheral
devices can be connected to the computer platform such as an additional
data storage device and a printing device.
It is to be further understood that, because some of the constituent
system components and method steps depicted in the accompanying
figures can be implemented in software, the actual connections between
the systems components (or the process steps) may differ depending
upon the manner in which the present invention is programmed. Given
the teachings of the present invention provided herein, one of ordinary
skill in the related art will be able to contemplate these and similar
implementations or configurations of the present invention.
While the present invention has been described in detail with reference
to a preferred embodiment, those skilled in the art will appreciate
that various modifications and substitutions can be made thereto
without departing from the spirit and scope of the invention as
set forth in the appended claims. |