WO2000045326A1  Method and apparatus for processing images with curves  Google Patents
Method and apparatus for processing images with curves Download PDFInfo
 Publication number
 WO2000045326A1 WO2000045326A1 PCT/US2000/001631 US0001631W WO0045326A1 WO 2000045326 A1 WO2000045326 A1 WO 2000045326A1 US 0001631 W US0001631 W US 0001631W WO 0045326 A1 WO0045326 A1 WO 0045326A1
 Authority
 WO
 WIPO (PCT)
 Prior art keywords
 curve
 image data
 computer
 cause
 module configured
 Prior art date
Links
 238000004422 calculation algorithm Methods 0 abstract claims description 32
 230000000875 corresponding Effects 0 abstract claims description 24
 210000003484 anatomy Anatomy 0 claims description 15
 238000004519 manufacturing process Methods 0 claims 8
 238000000034 methods Methods 0 description 8
 210000004556 Brain Anatomy 0 description 6
 230000001054 cortical Effects 0 description 4
 241000282414 Homo sapiens Species 0 description 3
 238000002583 angiography Methods 0 description 3
 230000037361 pathway Effects 0 description 3
 210000001772 Blood Platelets Anatomy 0 description 2
 230000000996 additive Effects 0 description 2
 239000000654 additives Substances 0 description 2
 238000002591 computed tomography Methods 0 description 2
 210000004884 grey matter Anatomy 0 description 2
 238000006011 modification Methods 0 description 2
 230000004048 modification Effects 0 description 2
 230000001423 neocortical Effects 0 description 2
 238000005457 optimization Methods 0 description 2
 230000035945 sensitivity Effects 0 description 2
 210000004885 white matter Anatomy 0 description 2
 210000004720 Cerebrum Anatomy 0 description 1
 208000004981 Coronary Disease Diseases 0 description 1
 240000008570 Digitaria exilis Species 0 description 1
 206010061216 Infarction Diseases 0 description 1
 241000124008 Mammalia Species 0 description 1
 210000000478 Neocortex Anatomy 0 description 1
 230000001174 ascending Effects 0 description 1
 230000003935 attention Effects 0 description 1
 230000017531 blood circulation Effects 0 description 1
 238000004590 computer program Methods 0 description 1
 201000008739 coronary artery diseases Diseases 0 description 1
 230000002596 correlated Effects 0 description 1
 238000000354 decomposition Methods 0 description 1
 238000002059 diagnostic imaging Methods 0 description 1
 238000009826 distribution Methods 0 description 1
 238000003708 edge detection Methods 0 description 1
 230000000694 effects Effects 0 description 1
 230000001747 exhibited Effects 0 description 1
 239000000284 extracts Substances 0 description 1
 238000003384 imaging method Methods 0 description 1
 230000036039 immunity Effects 0 description 1
 238000005304 joining Methods 0 description 1
 238000002595 magnetic resonance imaging Methods 0 description 1
 239000000203 mixtures Substances 0 description 1
 230000000877 morphologic Effects 0 description 1
 210000000653 nervous system Anatomy 0 description 1
 238000005381 potential energy Methods 0 description 1
 230000001603 reducing Effects 0 description 1
 230000002104 routine Effects 0 description 1
 239000007787 solids Substances 0 description 1
 230000001629 suppression Effects 0 description 1
 238000007514 turning Methods 0 description 1
 230000002792 vascular Effects 0 description 1
Classifications

 G—PHYSICS
 G06—COMPUTING; CALCULATING; COUNTING
 G06K—RECOGNITION OF DATA; PRESENTATION OF DATA; RECORD CARRIERS; HANDLING RECORD CARRIERS
 G06K9/00—Methods or arrangements for reading or recognising printed or written characters or for recognising patterns, e.g. fingerprints
 G06K9/62—Methods or arrangements for recognition using electronic means
 G06K9/6201—Matching; Proximity measures
 G06K9/6202—Comparing pixel values or logical combinations thereof, or feature values having positional relevance, e.g. template matching
 G06K9/6203—Shifting or otherwise transforming the patterns to accommodate for positional errors
 G06K9/6204—Matching of contours
Abstract
Description
METHOD AND APPARATUS FOR PROCESSING IMAGES WITH CURVES
RELATED APPLICATION
This patent application claims priority to U.S. provisional patent application 60/117,591 filed January 27, 1999, which is herein incorporated by reference in its entirety.
GOVERNMENT RIGHTS
This work was supported in part by the following U.S. Government grants: NIH grant 50567; NSF grant BIR 9424264; and NSF grant R01 NS3536802. The U.S. Government may have certain rights in the invention.
BACKGROUND OF THE INVENTION
The present invention relates to image processing systems and methods, and more particularly to image processing systems that process images with curves. Curves in an image provide information regarding underlying structures in the image. Curves that are apparent in brain images, for example, provide physicians with useful landmarks that can be used for treating patients. Pronounced gross morphological features of the cerebral hemisphere in mammals are known to manifest a diverse and complex arrangement of the sulcal fissures visible throughout the cortical surface of a mammalian brain. Some atlases of human anatomy catalog curved features comprising major sulci and gyri in images.
Computational metrics defined by cortical geometry such as geodesic length has attracted the attention of the neuroscience community. These special types of curved structures in the anatomy provide information for studying, among other things, the role of wiring length in the general layout of the nervous system. Despite their anatomic and functional significance curves in anatomical images such as the gyri, sulci, and many stable cortical pathways consistently appearing in the images often exhibit pronounced variability in size and configuration.
Conventional approaches for locating and characterizing curves in images are lacking because of limitations such as computational inefficiency, noise sensitivity, and an inability to accurately model curves of interest in images. There is, therefore, a need for an image processing system that overcomes the limitations of the conventional techniques.
SUMMARY OF THE INVENTION
The present invention provides a methodology for processing images with curves. Additional features and advantages of the invention will be set forth in the description which follows, and in part, will be apparent from the description, or may be learned by practicing the invention. The objectives and other advantages of the invention will be realized and obtained by the method and apparatus particularly pointed out in the written description and the claims hereof as well as in the appended drawings.
To achieve these and other advantages and in accordance with the purpose of the invention, as embodied and broadly described, a system according to the invention identifies image data points defining a curve. The method comprises the steps of determining a start point and an end point for the curve, establishing a search space that includes at least the start point, the end point, and other image data elements comprising the curve, and searching the search space using a dynamic programming algorithm to locate image data elements corresponding to the curve.
Another embodiment consistent with the present invention identifies image data points defining a curve. The method comprises the steps of determining a start point and an end point for the curve, generating a model of the curve, establishing a search space that includes at least the start point, the end point, and other image data elements comprising the curve, and searching the search space using a dynamic programming algorithm and the model for the curve to locate image data elements corresponding to the curve.
Yet another embodiment consistent with the present invention matches a first curve to a second curve. The method comprises the steps of identifying a first curve, identifying a second curve, generating a higher order distance measure for comparing the first curve and the second curve, and matching the first curve to the second curve using the higher order distance measure. Both the foregoing general description and the following detailed description are exemplary and explanatory and are intended to provide further explanation of the invention as claimed.
DESCRIPTION OF THE FIGURES
The accompanying drawings, which are incorporated in and constitute part of the specification, illustrate a presently preferred embodiment of the invention and together with the general description given above and detailed description of the preferred embodiment given below, serve to explain the principles of the invention.
Fig. 1 is an embodiment of an apparatus for processing curves in images consistent with the present invention;
Fig. 2 is a flow diagram for identifying a target curve in image data consistent with the present invention;
Fig. 3 is a flow diagram for identifying a target curve in image data using a curve model consistent with the present invention;
Fig. 4 is a flow diagram of a method for curve matching consistent with the present invention;
Fig. 5 is a schematic of a method for curve matching;
Fig. 6 is a schematic of a method for curve matching consistent with the present invention; and
Fig. 7 is an embodiment of an apparatus for processing curves in medical images consistent with the present invention.
DETAILED DESCRIPTION OF THE INVENTION
A system is disclosed which locates and matches curves in an image. Fig. 1 is an embodiment of an apparatus for processing curves in images consistent with the present invention. Curve processor 100 includes curve processor software 102 and curve processor hardware 104. Curve processor 100 receives image data 106 for processing. Image data 106 is, for example, any image that contains image data elements, at least some of which representing curved structural features. This would include medical images (e.g., computed tomography and magnetic resonance images) of the brain and other anatomical structures. Curve processor hardware 104 includes at least one hardware processing element suitable for processing image data. Curve processor software 102 is computer program code for processing image data 106 using curve processor hardware 104. The image processing operations described in greater detail below are performed by curve processor software 102 and/or curve processor hardware 104.
One skilled in the art will recognize that a parallel computer platform having multiple processors is also a suitable hardware platform for the present invention, including, but not limited to, parallel machines and workstations with multiple processors. Curve processor 100 can be a single computer, or several computers can be connected through a communications network to create a logical curve processor with the operations of curve processor software 102 and curve processor hardware 104 distributed across the computers. The functions performed by curve processor 100 can also be executed by hardware and/or software associated with medical devices such as, for example, a surgical navigation device. Moreover, the computational operations described herein can be performed by a general purpose computer, or other computational devices in a standalone configuration or with the operations distributed among several devices in a distributed network, including, but not limited to, networks linked by the Internet. In an embodiment consistent with the present invention, curve processor 100 includes a user interface for operator input (not shown).
In an embodiment consistent with the present invention, curve processor 100 identifies curves in image data 106. Geodesies are one example of the type of curves that can appear in image data 106. Two important neocortical curves are length minimizing geodesic curves and the curves determined by the sulcal fissures. For identifying these types of curves, dynamic programming is adapted for optimization on triangulated surfaces. Fig. 2 is a flow diagram for identifying a target curve in image data 106. To identify the points in an image comprising a geodesic, two image data elements are identified corresponding to the start point, s, and the end point, t, of the curve (step 200). Points s and t can be identified manually by an operator using a user interface to curve processor 100. Alternatively, points s and t can be identified using image processing algorithms, such as filters and edge detection routines known in the art. Once s and t are defined, there are many possible paths through the image data to connect s and t, therefore, initially there are many possible paths in the image data that are candidates for the curve being identified by curve processor 100. These many possible paths comprise a search space (step 202). An embodiment consistent with the present invention uses dynamic programming to efficiently determine which of the many possible paths through the image data corresponds to a geodesic (step 204).
In an embodiment of the dynamic programming approach, curve processor 200 generates a twicedifferentiable representation of the neocortical surface at the interface of gray and white matter supporting tangent, curvature, and torsion calculations. A cortical surface Mis constructed from a triangulated graph placed at the interface of the white and gray matter neocortex. For each point ? e M, the surface is preferably expressed as the graph of the function z(x, y) =βx, y), such that it is locally quadratically approximated by flp) =f_{x}(p) =f_{y}(p) ^{~} 0, with
f ( ^{v} x, ' y " ) ' = — ( V f xx x ^{2} + 2  xy xy ^{1} + f yy y^{1 2}) /,'
with the 2 x 2 shape operator defined via the curvatures as
( f xx f xy
\ f xy f yy,
The maximum and minimum eigenvalues κ_{t}, κ_{2} are the principal curvatures of Mat/?, with the unit vector directions t, and t_{2} in which these extreme values occur called the principal directions. For example, the surface created consists of 15,000 triangles, with each neighborhood of size 812 neighbors.
One set of curved anatomical features suitable for identification by curve processor 100 are the deepest paths in the valleys of the sulci are called the fundus beds. Fundus beds resemble crest lines corresponding to points where the maximal absolute principal curvature has a local maximum. Crest lines and curves are the loci of points x e M where the maximal absolute principal curvature κ_{m∞}(x) has a local maximum, expressed mathematically as V κ_{mαI}(x) • t_{m}(x) = 0 for t_{m} the principal unit direction corresponding to the maximal principal curvature. To find image data elements corresponding to curves in a way that reduces sensitivity to noise in an image, the problem of tracking such curves is posed as a control/optimization problem. Following the dynamic programming approach, curve processor 100 searches for curves that pass through regions of curvature joining the start and end points (s, t) in the surface. For noise immunity, curve processor 100 preferably uses a sequentially additive energy associated with candidate curves and uses dynamic programming for its minimization. To apply dynamic programming, a cost function such as /α_{(}. _{r)}(κ_{mαr}(x)  Kfdx, with ST preferably assigned the largest maximal curvature on the surface, and minimize over all such paths on the triangulated graph representation of the surface. In images of the human anatomy using this cost function in conjunction with dynamic programming generates curves that agree with those generated by an anatomist attempting to outline similar anatomical structures.
Using dynamic programming, curve processor 100 transforms image data 106 into a search space, a finite state space S of size S = N and a graph representing the possible transitions between states. Curve processor 100 preferably computes the shortest paths between the specified initial states s and the final state t. For example, when the optimal path has no more than K nodes, the total number of paths of length K between points s and t are of the order N*. If the cost is additive over the length of the path, dynamic programming reduces the complexity of the search algorithm to order of KN ^{2}. Using c^{k}(x_{k}, x_{k+]}) to denote the cost incurred for the transition from state x_{k} e S to x_{k+]} e S at each time k, curve processor 100 selects the path π which minimizes the expected cost
JXx_{0}) = ∑ c ^^ ,
with the path π identified with the sequence of arcs (s =j],j_{2}), (JiJi), —, (J_{Λ}._{/}. , ^{=} *)■> from the initial node s of the path to the terminal node t. In an embodiment consistent with the present invention, suppression of k dependence in c(i, j) means the cost is independent of time. The graph is constructed to that it preferably has no negative arc lengths, so that c(i, j) > 0, and arcs of infinite cost c(i, j) = ∞ signifies that there is no arc from node i to nodej. An optimal path need not have more than N arcs (number of nodes in the graph) and hence takes no more than N moves. Curve processor 100 seeks to find the optimal path through the search space in N moves. Degenerate moves from a node i to itself with cost c(i, i) = 0 are also permitted. The degenerate moves signify that the length of the path may be less than N. Thus, denoting the optimal cost of getting from node i to node t in (Ν  i) moves as J_{k}(i), i e S, k = 0, 1 , , Ν  1, then the optimal Νlength path J_{0}(i) from i to t is given by the final step of the following algorithm, with J_{Ν}_, (i) = c^{N"}'(i, t), and J (t) = min, = 1 ... N {c^{k} (i, j) + J_{k+1} (j)}, k = 0,l, ..., N  2, i e S.
Therefore, following a dynamic programming approach consistent with the present invention, curve processor 100 generates a triangulated graph from image data 106 corresponding to a brain surface. The curve start and end points, s and t, are nodes in the triangulated graph. The triangulated graph is used to reduce the search space operated on by the dynamic programming algorithm to operate on the image data elements that represent the most likely components of the geodesic. This approach allows for more efficient curve identification.
The steps of the dynamic programming curve identification method for curve identification consistent with the present invention are now set forth in greater detail when applied to locating image curves corresponding to geodesies. Geodesies on the continuum surface correspond to length minimizing curves restricted to the surface. The triangulated graph M discussed above is adapted to have indices /, j, the index nodes in the graph, and x_{y} x_{r} e l^{3} in their positions in M^{3}. Accordingly, given a two dimensional triangulation {x. e M}, a platelet P, of point x, is the set of triangles (with index triples (j , , j_{2} , j_{3}) specifying their vertices) sharing x, as a common vertex
^{p}xΛό_{1}' J_{2}'
A path α with length L(α) on the triangulated surface M is sequence of arcs (j,, j_{2})> 0z» Js). •••» Oki k •■•> GNI N) such thaty, eP_{lk}_„ V&, with k=Nl {a) k Σ = l J_{k}'J J + 1 '
where d(i, j) = _{\}j ( x  ] ) ^{2} + ( ^{2}  x^{2} ^{2} )^^{2} + (x^x,^{3})^{2}
For s and t, nodes on the triangulated surface M, consider the collection of all paths α(s, t) e P_{st}(M) connecting (s, t), and define the discrete geodesic
ά(s, t) =^αeP_{Sft}(Λf) L(α) = rain L(β) \. βeP_{s}__{t}(W)
Given a curve passing through some point on the graph, the curve passes through one of its neighbors (analogous to being in the tangent space for the continuum representation). This structure of the triangulated graph allows for the pruning of the state space at each discrete time thereby reducing the complexity of the algorithm. Instead of searching over the whole state space S at time k, curve processor 100 examines the restricted state space S_{k} c S of points in the platelet which terminates at point t in at most N  k moves. Accordingly, the state spaces and costs functions are preferably represented as follows:
{i\ieP_{t}} , S_{k}={i\ieP., je S_{ktl}} , c^{k} (i, j) =d(i, j) , jeP. with c^{k} ( i , j ) = ∞ for j <£ P_{i} , with
J (i) = c (i,t), ieS ,
J_{k}(i) = (i, j) + J_{k+1}(j) . ieS_{k},
k = 0,l,...,N2.
Then the following Algorithm 1 with J_{0}(s) defined in the algorithm produces a shortest path (but not necessarily a unique path) between nodes s and t with length given by J_{0}(s): ά = argmin L {cx) α ( s, t) eP_{s t} (M)
ALGORITHM 1. Initialize: J_{k}(i) < ∞ i ≠ t, for all k, S_{N} < t , J_{k}(t) < 0: for k <— N  1 down to 0 do
S_{k}  { I j e P_{p}j e S_{t+1}}, set c^{k}(i, j),j e S_{M}, i <= S_{k}.
J_{k} ( i) = min { c^{k} ( i, j ) + J_{k +l} ( j ) } , ±e S_{k} .
The above Algorithm 1 produces the shortest curve between points on the triangulated graph M. Since the curve is constrained to the lattice points of the graph, the computed shortest path may be different from the true geodesic. For dealing with this issue, the search for the shortest path is extended to the edges of the triangles by dividing the edges of the triangle into equal parts. To each original site x„ assign as its neighbors all the sites which sit on the edges opposite to x, in triangles which have x, as a vertex. To the new set of sites which belong to the triangle edges, assign as neighbors, all sites lying on the four opposite edges to that edge thus creating new neighborhoods for each point. As refinement increases, the discrete geodesic α. approximates the true geodesic arbitrarily closely.
Another method consistent with the present invention identifies and tracks curves in image data 106 using dynamic programming and a curve model (fig. 3). As described above, a start point and an end point for the target curve is determined (step 300). Curve processor 100 generates a model of the geometry of the curve being tracked (step 302). Curve processor 100 establishes a search space that includes the start point, the end point, and other image data elements comprising the target curve (step 304). Curve processor 100 uses the curve model to focus the search performed by the dynamic programming algorithm in the search space (306).
One model consistent with the present invention for curve tracking is a Frenet representation model. A Frenet curve representation models a curve's speed, torsion, and curvature. For each curve x(s) parametized by arc length s, curve processor 100 associates the orthonormal tangent field T, the normal field N, and the binormal field R. Using the Frenet representation, if x() is a unit speed curve with curvature and torsion fields K, τthen T'(s) = κN(s), N'(s) =  κT(s) + τB(s), B'(s) =  τN(s), with curvature
and torsion
For the prior distribution, curve processor 100 associates the potential energy with the mean lengths, curvatures, and torsions 1 , K , τ with the Nlength paths assuming the simplest model of constant inverse variabilities η, β, γ across the entire fundus:
Ν
[7( 1 , K , τ ) = η ∑ T . 2
£ = 1 ^{' β} A∑ C K , K, Y . = l J = l ^{(}χ  χ
Parameters for the curve model can be estimated using example images by fitting polynomic splines to the surface nodes from which speeds, curvatures, and torsions are calculated. Empirical estimates of the mean lengths, curvature and torsion functions of the fundus beds can be computed from images according to the formulas:
1 M i . m  τ .
T& ^{1} k i K
M Σ — 1 ^ k ' m=l ^{k} M [
1 ^{M} M ^{N}N W _{m=}ι _{k = 1}
To synthesize the fundus curves, the random process ( , κ_{κ}, τ_{κ}, k = 1 , ..., N) is preferably generated to be Gaussian, with the discretized Frenet equations solved 0/45326 PCT/USM/Ol^{β}l
11 sequentially. Suitable translation and rotation is chosen to match the start points of the curve removing ambiguity of the special Euclidean group. One skilled in the art will recognize that there are several other models for representing curves within the scope of the present invention.
After modeling the curve to be tracked, the curve model is incorporated into the dynamic programming algorithm. The dynamic programming algorithm and curve model are integrated preferably according to the following approach. The cost of a candidate curve (s, t) e P(s, t) is defined as J (s,t) (κ_{max}(x)  κ)^{2}dx, with K assigned the largest maximal curvature of the surface. For finding gyral crowns, use the extrema of negative curvatures. For a surface symmetrical about a crest line, i.e., one where t_{m} is perpendicular to the crest line, minimizing the above function gives (κ_{max} (x)  K) V κ_{max} ^{•} t_{m} = 0 implying V κ_{max} ^{•} t_{m} = 0, which is precisely the equation for the crest line. For regions of the sulcus where the basin is flat, κ_{max} is constant, and the minimizer of the above functional produces shortest paths through these regions.
Curve processor 100 incorporates the prior shape of fundi into the dynamic programming algorithm by associating with a candidate path a(s, t) e the set of node positions with principal maximal curvatures of the surface through which the path traverses on the surface as (x^ K_{max}(xJ), k = 1 , ..., N, and the sequence of lengths, curvatures, and torsions (l_{h} κ_{h} τ_{k}), k = 1, ..., N of the curve. Using Bayes law, the posterior potential of a discrete path (s, t) through the triangulated graph connecting (s, t) becomes
K ( X k '  K "/ I + ( \\K" mmaaxx ( ^{vv} X " JJ. ++ 11 J ' '  K) /^{2} ^{ 2} k Σ
. ( \K max { X k, )  K max ( * X k. + .1 ) ' ) '^{'} 00/45326
^{■}12 where Δ  + l k k + 1 and σ^{2} is a weight determining the relative contribution of each term, and curve processor 100 preferably uses the corrected trapezoid rule for the integration. Dynamic programming preferably maximizes over the space of curves (s, t) e R „_{j}(M). Therefore, given s and t the start and end point of the curves α(s, t) e P _{s t}(M), then the Algorithm 1 , with the cost for transition c^{k}(i, j) = ∞, j <f .P_{i} and forj e 3^{>} _{i},
x ^{)}  )^{2 +} (^{K}maxk _{+} lM^{2}
X , J ■ r
+ \ max k ^{■} K max ( ^{%} X k, .X Δ
produces the maximum a posteriori estimator of H( ),
ά ( s, t) , t ( M) H (α) =
Another operation performed by curve processor 100 is curve matching. Curves in image data provide useful landmark for image processing applications such as image registration. In image registration, corresponding features in two or more images are correlated to provide a map that relates image data elements among the registered images. Curve matching provides a tool for generating a coarse image data correlation in the registration process by correlating corresponding curves appearing in the images to be registered.
Fig. 4 is a flow diagram of a method for curve matching consistent with the present invention. Curve processor 100 identifies a first curve in a first image (step 400). Curve processor 100 identifies a second curve in a second image (step 402). Note also that the first and second curves can also be in the same image. Curve processor 100 then generates a higher order distance measure for comparing the two curves, such a distance measure using a Frenet representation of the curves (step 404). Curve processor 100 matches the first cirve to the second curve using the distance measure (step 406). The foregoing steps are described in greater detail below.
Given a curve α(t), t e [0,1], the Frenet equations provide the differential geometric characterization of α in three dimensions (3D) based on its speed, curvature, torsion, and the orthogonal frame (T, N, B):
Here T is the unit tangent vector field, N is the unit normal vector field, B is the unit binormal field on α, and the speed, curvature, and torsion parameters are given by: ^{:} α', K = α' x c"/α'^{3}, τ = (α' x α") • α' " /α' x α". F(t) describes the flow of the orthogonal frame through its tangent space.
Given two curves α(t) and β(t) and their orthogonal frames F_{α}(t), F_{p}(f) parameterized on the unit interval [0, 1], define the set of diffeomorphisms the index set [0, 1] to itself as
Φ = [φ: [0, l] ^{«} [0, l]]
A distance p(α, β; φ) is defined using the Frobenius norm between any two 3^{χ}3 matrices A and B by: trace (A  B) (A  B)^{T}.
The distance p between the curves and β is: p(α, β; Φ) = P trace (F_{α}(t)  F_{p}(t))(F_{B}(t)  F_{p}(t))^{T}dt
= 2 J ^(v.CtKCt)  v_{p}(φ(t))κ_{p}(φ(t)))^{2}dt
+ 2 J ' (v„(t)T_{e}(t)  v_{p}(φ(t))T_{p}(φ(t)))^{2}dt
Due to the properties of the Frenet representation, the distance p is inherently invariant to spatial position and orientation. The problem of matching curves across different brains is now defined as one of finding the particular diffeomorphism that minimizes the above distance. WOOΘ/45326 PCT/USOO/01631
14
Given two curves α(t), β(t), t e [0, 1] and the set of diffeomorphism Φ = lφ:φ(α(t)) = β(t)J, curve processor 100 matches the two curves by finding the diffeomorphism that minimizes the above distance with a penalty on the velocity field:
Φ .α, β ) argmin^_{φ} v„ v ( t ) ^{l}dt + p (α, β ; φ )
Introducing the simple normsquare function   v  v_{a} \  ^{2} in the matching provides control over local scale. For generating the lowest cost diffeomorphism, the cost function p^{Δ} is the discrete approximation to the distance. Curve processor 100 preferably uses a bipartite matching algorithm to reduce the complexity of the problem [see, e.g., Sedgewick, "Algorithms," (1983)]. Representing the curve as a linear array of straight line segments and assuming that curvature and torsion are piecewise constant, the distance assigned to the correspondence becomes: p^{Δ}( , β; φ) = a ∑ [ v j ) κ_{a} ( j )

where a, b, and c are coefficients picked by an anatomist or determined by curve processor 100 to adjust the weight of the matching based on the speed, curvature, or torsion terms. Choosing to match based on speed (a ≠ 0, b = c = 0) emphasizes uniform stretching of the curves. Matching based on curvature (b ≠ 0, a = c = 0) emphasizes the turning of the curves in the plane. Matching based on torsion (c ≠ 0, a = b = 0) emphasizes twisting of the curves out of the plane. Matching based on a weighted combination of these criteria is also possible. In the bipartite matching process implemented by curve processor 100, the target curve α is preferably sampled with n equally spaced points and the template curve β is sampled with m points (m is preferably equal to or approximately equal to N(n  1), where N is the spacing between the samples in the target). The neighborhood of each point i in the target curve is defined to be al points j e β such that [i  j] < N. We now have a bipartite weighted graph in which there are two distinct set of nodes (samples), and all edges in the graph connect two samples i e α, j e β where j is defined to be a neighbor of i. The weights on the edges are defined by: (i, j) = a[v_{α}(i)  v_{p}(i)]^{2} + b[v_{«}(i)κ_{β}(i) vp(j)κ_{p}G)]^{2} + c[v_{β}(i)τ„(i)  v_{p}α)τ_{p}0)]^{2}.
Let i represent one of the n samples in the target. For each such sample, there is an associated cost array, where cost[i, j] = w(i, j) if the edge exists (j is in the neighborhood of i). cost[i, j] = ∞ if the edge does not exist. There is also a prefer[i, k] array for each sample i in the target, which contains the indices of the template samples j sorted by their costs in ascending order. For example, if min.cost[i,j] = j' then prefer [i, 1] = j', i.e., if j' is the template point (among all the template points j) for which the cost function for the target point i is minimized, then j' becomes the first point on point i's preference list. Curve processor 100 also keeps track of how far down each point in the target has progressed in its preference list. This is handled by the index[i] array, initialized to 1. The current match in the template assigned to point i is stored in the array match[i]. The algorithm proceeds as follows:
1. Initialize: index[i] = 1 Vi e 1, ... n, match [1] = prefer [1, 1]
2. for i = 2 to n do repeat isLegal = FALSE match[i] = preferfi, index[i]] if matchfi] > match[i  1] isLegal = TRUE else find all p and q such that
*p ≥ index [i  1] and q ≥ index [i] where (p, q) ≠(match [i  1], matchfi]) *match[i] < prefer [i  1], p], prefer [i, q] < matchfi  1 ] and δ,l(p) = cost[i 1, match[i  1]]  cost[i  1, prefer [i  1, p] δ.(q) = cost[i, match [i]]  cost[i, prefer[i, q]]) Δ(p, q) = δ,._{1}(p) + δ,(q) p* = argmiiyHp, q) q* = argmin_{p}Δ(p, q) index[i  1] = p* index [i] = q* if i ≠ 1 i = i  1 WαOO/45326 PCT/USOO/01631
16 until isLegal = TRUE.
The matching algorithm assigns each sample in the target the match with the lowest possible cost. If the match does not violate the diffeomorphism, the inner loop terminates. If it does, the matching algorithm finds the next match that does not violate the diffeomorphism for the last two points and maintains the lowest cost constraint. It then steps back to check if the new match violates the diffeomorphism for the previous points.
Fig. 5 is a schematic of a method for curve matching using only a single parameter, arc length. Fig. 5 includes two gaussian curves to be matched, target curve 502 and template curve 504. The curves have been generated with different means and standard deviations. The curves are representative of curves that appear in brain images. Template curve 504 was evaluated at points more than two standard deviations from the mean to create a relatively flat portion of the curve to the right of the peak. Target curve 502 has thirteen equally spaced landmark points 506am.
After executing a matching algorithm comparing the arc length between points in target curve 502 and template curve 504, points corresponding to target image 502 points 506am are shown on template curve 504 as points 508am. The limitations of using arc length alone for curve matching are apparent from the correlation of target points 506am to template points 508am. For example, some points to the left of the peak in target curve 502 are matched with points to the right of the peak in template curve 504. Also note that some points along the descending slope of target image 502 are matched to points along the flat region of template curve 504, even though target curve 502 does not have a flat region.
The limitations in the curve matching in Fig. 5 can be addressed by including higher order information in the matching algorithm. As described in greater detail above, this higher order information includes, for example, using a distance measure that includes information about a curve's curvature and torsion. Fig. 6 shows the result of applying a method for curve matching consistent with the present invention using higher order information (e.g., curvature and torsion) in the matching process. 17
Target curve 602 has landmark points 606am. Target curve 602 and template curve 604 are gaussian curves as described for Fig. 5. After applying the matching algorithm using a distance measure that includes curvature and torsion, points 608am in template curve 604 were matched to points 606am in the target curve, respectively. In contrast to the matching shown in Fig. 5, Fig. 6 shows that points matched in the target and template curves correspond to the same side of the curve peak. Also, no points in the target curve are mapped to the flat region of the template. Thus, using this higher order in the matching distance measure can improve matching.
Locating and matching curves in image data is useful in many image processing applications including medical image processing. An embodiment consistent with the present invention generates vascular connections by generating curves in imagery of any dimension from intensity imagery such as CTangiography or MRangiography. Since it is based on dynamic programming it can generate a globally optimal decomposition of the solid volume into arterial pathways. Fig. 7 is a block diagram of an embodiment of a curve processor for medical images consistent with the present invention. Image data 706 is medical image data containing curve features such as CTangiography. Since curve processor 100 can process three dimensional imagery, by using the three dimensional geometry of the images, curve processor 100 can process data independently from the viewing orientation and is therefore not limited by the viewer. Curve processor 100 extracts curves from image data 706 and generates data 704 containing curves representing arterial pathways 702.
Another embodiment calculates a structural representation of the blood flow through a threedimensional curve in angiographic images by calculating the local average intensity in the arterial enclosed region providing an automated diagnostic tool for measuring flow versus arclength. Another embodiment of curve processor 100 provides a threedimensional, real time viewing instrument of arbitrary geometry vasculature. Given the structural description of the orientation of vessels and their extent, curve processor 100 generates view of the structure by removing the obscuring effects of twodimensional projections. Using dynamic programming curve processor 100 generates multiple arterial projections and handles bifurcations. Moreover, curve processor 100 can produce arterial flow diagrams illustrating blockage as a function of WO00/45326 PCT/USOO 
18 position along a structural representation. Such a description provides an analyzer for arterial blockage and ballooning such as is associated with coronary disease and arterial infarctions.
While the disclosed system and method is useful for medical imaging systems used for noninvasive exploration of human anatomy, for example, computed tomography and magnetic resonance imaging, this invention can also be used on images acquired from other imaging modalities. Furthermore, application of the present invention is not limited to anatomical images. This invention also applies to nonanatomical images, including, but not limited to, satellite imagery, photographs, radar images, and images acquired from multiple sources.
It will be apparent to those skilled in the art that various modifications and variations can be made to the disclosed embodiments of the present invention without departing from the spirit or scope of the invention. Thus it is intended that the present invention cover the modifications and variations of this invention provided they come within the scope of the appended claims and their equivalents.
Claims
Priority Applications (2)
Application Number  Priority Date  Filing Date  Title 

US11759199P true  19990127  19990127  
US60/117,591  19990127 
Applications Claiming Priority (1)
Application Number  Priority Date  Filing Date  Title 

AU26247/00A AU2624700A (en)  19990127  20000127  Method and apparatus for processing images with curves 
Publications (1)
Publication Number  Publication Date 

WO2000045326A1 true WO2000045326A1 (en)  20000803 
Family
ID=22373754
Family Applications (1)
Application Number  Title  Priority Date  Filing Date 

PCT/US2000/001631 WO2000045326A1 (en)  19990127  20000127  Method and apparatus for processing images with curves 
Country Status (2)
Country  Link 

AU (1)  AU2624700A (en) 
WO (1)  WO2000045326A1 (en) 
Cited By (11)
Publication number  Priority date  Publication date  Assignee  Title 

DE10319546A1 (en) *  20030430  20041125  Siemens Ag  Automatic vein structure anomaly detection procedure recognises vein structure in tomography data, prepares skeletal path and classifies by comparison with characteristic values 
EP1280096A3 (en) *  20010723  20050810  Siemens Aktiengesellschaft  Method for matching two digital images consisting of polygonal stroke sequences 
WO2008141046A1 (en) *  20070509  20081120  Edsa Micro Corporation  Systems and methods for creation of a schematic user interface for monitoring and predicting the real time health, reliability and performace of an electrical power system 
US8131401B2 (en)  20060719  20120306  Power Analytics Corporation  Realtime stability indexing for intelligent energy monitoring and management of electrical power network system 
US8155943B2 (en)  20071012  20120410  Power Analytics Corporation  Systems and methods for automatically converting CAD drawing files into intelligent objects with database connectivity for the design, analysis, and simulation of electrical power systems 
US8170856B2 (en)  20060412  20120501  Power Analytics Corporation  Systems and methods for realtime advanced visualization for predicting the health, reliability and performance of an electrical power system 
US8229722B2 (en)  20070516  20120724  Power Analytics Corporation  Electrical power system modeling, design, analysis, and reporting via a clientserver application framework 
US8253723B2 (en)  20050622  20120828  Koninklijke Philips Electronics N.V.  Method to visualize cutplanes for curved elongated structures 
US8775934B2 (en)  20060719  20140708  Power Analytics Corporation  Systems and methods for creation of a schematic user interface for monitoring and predicting the realtime health, reliability and performance of an electrical power system 
US9092593B2 (en)  20070925  20150728  Power Analytics Corporation  Systems and methods for intuitive modeling of complex networks in a digital environment 
US9557723B2 (en)  20060719  20170131  Power Analytics Corporation  Realtime predictive systems for intelligent energy monitoring and management of electrical power networks 
Citations (4)
Publication number  Priority date  Publication date  Assignee  Title 

US5734739A (en) *  19940531  19980331  University Of Washington  Method for determining the contour of an in vivo organ using multiple image frames of the organ 
US5768413A (en) *  19951004  19980616  Arch Development Corp.  Method and apparatus for segmenting images using stochastically deformable contours 
US5841958A (en) *  19950119  19981124  Nec Research Institute, Inc.  Bipartite matching 
US5898797A (en) *  19940415  19990427  Base Ten Systems, Inc.  Image processing system and method 

2000
 20000127 AU AU26247/00A patent/AU2624700A/en not_active Abandoned
 20000127 WO PCT/US2000/001631 patent/WO2000045326A1/en active Application Filing
Patent Citations (4)
Publication number  Priority date  Publication date  Assignee  Title 

US5898797A (en) *  19940415  19990427  Base Ten Systems, Inc.  Image processing system and method 
US5734739A (en) *  19940531  19980331  University Of Washington  Method for determining the contour of an in vivo organ using multiple image frames of the organ 
US5841958A (en) *  19950119  19981124  Nec Research Institute, Inc.  Bipartite matching 
US5768413A (en) *  19951004  19980616  Arch Development Corp.  Method and apparatus for segmenting images using stochastically deformable contours 
Cited By (13)
Publication number  Priority date  Publication date  Assignee  Title 

EP1280096A3 (en) *  20010723  20050810  Siemens Aktiengesellschaft  Method for matching two digital images consisting of polygonal stroke sequences 
DE10319546A1 (en) *  20030430  20041125  Siemens Ag  Automatic vein structure anomaly detection procedure recognises vein structure in tomography data, prepares skeletal path and classifies by comparison with characteristic values 
US7546154B2 (en) *  20030430  20090609  Siemens Aktiengesellschaft  Method and apparatus for automatic detection of anomalies in vessel structures 
US8253723B2 (en)  20050622  20120828  Koninklijke Philips Electronics N.V.  Method to visualize cutplanes for curved elongated structures 
US8170856B2 (en)  20060412  20120501  Power Analytics Corporation  Systems and methods for realtime advanced visualization for predicting the health, reliability and performance of an electrical power system 
US8131401B2 (en)  20060719  20120306  Power Analytics Corporation  Realtime stability indexing for intelligent energy monitoring and management of electrical power network system 
US8775934B2 (en)  20060719  20140708  Power Analytics Corporation  Systems and methods for creation of a schematic user interface for monitoring and predicting the realtime health, reliability and performance of an electrical power system 
US9557723B2 (en)  20060719  20170131  Power Analytics Corporation  Realtime predictive systems for intelligent energy monitoring and management of electrical power networks 
WO2008141046A1 (en) *  20070509  20081120  Edsa Micro Corporation  Systems and methods for creation of a schematic user interface for monitoring and predicting the real time health, reliability and performace of an electrical power system 
AU2008251610B2 (en) *  20070509  20130523  Edsa Micro Corporation  Systems and methods for creation of a schematic user interface for monitoring and predicting the real time health, reliability and performace of an electrical power system 
US8229722B2 (en)  20070516  20120724  Power Analytics Corporation  Electrical power system modeling, design, analysis, and reporting via a clientserver application framework 
US9092593B2 (en)  20070925  20150728  Power Analytics Corporation  Systems and methods for intuitive modeling of complex networks in a digital environment 
US8155943B2 (en)  20071012  20120410  Power Analytics Corporation  Systems and methods for automatically converting CAD drawing files into intelligent objects with database connectivity for the design, analysis, and simulation of electrical power systems 
Also Published As
Publication number  Publication date 

AU2624700A (en)  20000818 
Similar Documents
Publication  Publication Date  Title 

Van den Elsen et al.  Automatic registration of CT and MR brain images using correlation of geometrical features  
Ng et al.  Medical image segmentation using kmeans clustering and improved watershed algorithm  
Shen et al.  MRI fuzzy segmentation of brain tissue using neighborhood attraction with neuralnetwork optimization  
Maintz et al.  Evaluation of ridge seeking operators for multimodality medical image matching  
Hartley  Projective reconstruction and invariants from multiple images  
Bankman  Handbook of medical image processing and analysis  
Stefanescu et al.  Grid powered nonlinear image registration with locally adaptive regularization  
Pizer et al.  Segmentation, registration, and measurement of shape variation via image object shape  
Chakraborty et al.  Deformable boundary finding in medical images by integrating gradient and region information  
Koss et al.  Abdominal organ segmentation using texture transforms and a hopfield neural network  
Kirbas et al.  Vessel extraction techniques and algorithms: a survey  
US6754374B1 (en)  Method and apparatus for processing images with regions representing target objects  
Schmitt et al.  New methods for the computerassisted 3D reconstruction of neurons from confocal image stacks  
US6785409B1 (en)  Segmentation method and apparatus for medical images using diffusion propagation, pixel classification, and mathematical morphology  
Li et al.  Fusing images with different focuses using support vector machines  
Paragios et al.  Nonrigid registration using distance functions  
Davatzikos et al.  Hierarchical active shape models, using the wavelet transform  
Van Ginneken et al.  Segmentation of anatomical structures in chest radiographs using supervised methods: a comparative study on a public database  
Lorenzen et al.  Multimodal image set registration and atlas formation  
US20020141626A1 (en)  Automated registration of 3D medical scans of similar anatomical structures  
Wyawahare et al.  Image registration techniques: an overview  
Pitiot et al.  Expert knowledgeguided segmentation system for brain MRI  
Zoroofi et al.  Automated segmentation of acetabulum and femoral head from 3D CT images  
EP2284795A2 (en)  Quantitative analysis, visualization and motion correction in dynamic processes  
Kim et al.  A fully automatic vertebra segmentation method using 3D deformable fences 
Legal Events
Date  Code  Title  Description 

AK  Designated states 
Kind code of ref document: A1 Designated state(s): AE AL AM AT AU AZ BA BB BG BR BY CA CH CN CR CU CZ DE DK DM EE ES FI GB GD GE GH GM HR HU ID IL IN IS JP KE KG KP KR KZ LC LK LR LS LT LU LV MA MD MG MK MN MW MX NO NZ PL PT RO RU SD SE SG SI SK SL TJ TM TR TT TZ UA UG UZ VN YU ZA ZW 

AL  Designated countries for regional patents 
Kind code of ref document: A1 Designated state(s): GH GM KE LS MW SD SL SZ TZ UG ZW AM AZ BY KG KZ MD RU TJ TM AT BE CH CY DE DK ES FI FR GB GR IE IT LU MC NL PT SE BF BJ CF CG CI CM GA GN GW ML MR NE SN TD TG 

121  Ep: the epo has been informed by wipo that ep was designated in this application  
DFPE  Request for preliminary examination filed prior to expiration of 19th month from priority date (pct application filed before 20040101)  
REG  Reference to national code 
Ref country code: DE Ref legal event code: 8642 

122  Ep: pct application nonentry in european phase 