MULTISCALE MORPHOLOGICAL TOPOLOGY CORRECTION OF

CORTICAL SURFACES

FIELD OF THE INVENTION

The present invention relates generally to the field, of medical imaging and associated data processing. More specifically, the invention relates to improved techniques for multiscale volume-based retrospective topology correction in cortical surface reconstruction.

BACKGROUND OF THE INVENTION

Medical imaging techniques, and particularly magnetic resonance imaging (MRI), is a valuable tool in studying both the structure and the function of human cortex. It is also valuable to reconstruct a representation of the cortical surface from an acquired image. The reconstructed surface is the boundary between foreground voxels in the image, which correspond to bone and tissue, and background voxels in the image. Although a normal human cortex is highly convoluted, the cortical surface (with the brain stem artificially closed) is topologically equivalent to a sphere. However, due to various artifacts in the acquired images, the reconstructed cortical surface may contain various topological defects such as handles and cavities. Consequently, researchers have developed various techniques for reconstructing cortical surfaces that are topologically correct.

Techniques for generating topologically correct cortical surfaces fall into two broad classes: deformable model methods and retrospective topology correction methods. Deformable model methods begin with a topologically correct surface model. The surface geometry is then deformed without altering the topology. This approach is performed as an integral part of the image segmentation process that distinguishes foreground from background voxels. Retrospective topology correction methods take as given a segmented image with foreground and background voxels identified. In a post-processing step after image segmentation, the topology is corrected using a topology correction technique that eliminates all topological defects until the topology is equivalent to that of a sphere.

Retrospective topology correction techniques may be divided into two types: surface-based approaches and volume-based approaches. In a surface-based approach (see Fig. 5), the cortical surface model 502 is first extracted from the segmented image volume 500 and then the topology of the surface is corrected 504 until it is homotopic to a sphere, e.g., Guskov[l], Wood[10], and Jaume[4,5]. In a volume-based approach (see Fig. 6), the topology of the segmented image volume 600 is first corrected until it is homotopic to a ball 602 and then the surface of the topologically correct volume is extracted 604, e.g., Shattuck[9], Han[2,3], Segonne[7,8], and Kriegeskorte[6]. The techniques of Shattuck and Han are examples of multiscale volume-based approaches that gradually correct topology errors at increasing length scales. At each length scale, topological defects (e.g., handles) are eliminated using a cost function that depends on the scale. In all existing techniques, topological defects are eliminated either from the foreground voxels, the background voxels, or both. In Shattuck, handle elimination cost is measured in terms of the number of removed voxels. In Han, the cost of handle elimination is mainly measured by the distance-to-surface metric. A graph-based method is then used to determine the final handle cuts attempting to minimize the number of removal voxels at a specific scale.

[1] Guskov, I. and Wood, Z. Topological noise removal. In Graphics Interface 2001.

[2] Han, X., Xu ₅ C, Braga-Neto, U. and Prince, J. L. Topology Correction in Brain Cortex Segmentation Using a Multiscale, Graph-Based Algorithm. IEEE Trans. Medical Imaging, 21 (2002) 109-121.

[3] Han, X., Xu, C, Braga-Neto, U. and Prince, J. L. Graph-based topology correction for brain cortex segmentation, in Proc. XVIIth Int. Conf. Information Processing in Medical Imaging, June 2001.

[4] Jaume, S., Rondao, P., Macq, B. Open Topology: A Toolkit for Brain Isosurface Correction, MICCAI Open Source Software Workshop, Oct 2005.

[5] Jaume, S. Topology Simplification Algorithm for the Segmentation of Medical Scans, Ph.D. Dissertation, University of Louvain (Belgium), Feb 2004.

[6] Kriegeskorte, N. and Goebel R. An Efficient Algorithm for Topologically Correct Segmentation of the Cortical Sheet in Anatomical MR Volumes. Neurolmage, 14 (2001) 329- 346.

[7] Segonne, F., Grimson, E., and Fischl, B. Topological correction of subcortical segmentation. In Proc. MICCAF03, Montreal, November 2003, 695-702.

[8] Segonne, F. Segmentation of Medical Images under Topological Constraints, PhD dissertation, Dec. 2005.

[9] Shattuck, D. W. and Leahy, R. M. Automated Graph-based Analysis and Correction of Cortical Volume Topology. IEEE Trans. Medical Imaging, 11 (2001) 1167-1177.

[10] Wood, Z., Hoppe, H., Desbrun, M. and Schroder, P. Removing Excess Topology. ACM Trans. Graphics, 23 (2004) 190-208.

SUMMARY OF THE INVENTION

In one aspect, the present invention provides a computer-implemented method for improved multiscale volume-based retrospective topology correction. In contrast with prior methods which correct topological defects on at most two objects, namely, the foreground and/or background voxels, the present method corrects topological defects on three objects, namely, white matter, grey matter, and background voxels. The method includes obtaining a segmented volumetric image comprising white matter (W), grey matter (G) and background (B) voxels. The image may be preprocessed if necessary to ensure that W consists of one connected component and contains no cavities. At each of a sequence of increasing scales, a topology correction filter is successively applied to fill tunnels and cut handles of W. When W is nomotopic to a ball, the topology correction is completed and a cortical surface model is extracted from the topology corrected volumetric image. The topology correction technique preferably comprises a morphological topology correction (MTC) filter that associates a larger correction cost to topology corrections at locations having a larger surface-like property (e.g., "flat" and "wide" regions).

BRIEFDESCRIPTION OF THE DRAWINGS

Figs. IA-B show 3D images illustrating topological corrections and surface-likeness properties of objects according to the teachings of the invention.

Figs. 2A-B show 3D images illustrating various point types according to the teachings of the invention.

Figs. 3 A-C illustrate a dilation procedure according to the teachings of the invention.

Figs. 4A-C show images illustrating regions of a cortex before and after topology correction according to the teachings of the invention.

Fig. 5 is a schematic diagram illustrating a surface-based approach to retrospective topology correction.

Fig. 6 is a schematic diagram illustrating a volume-based approach to retrospective topology correction.

Fig. 7 is a flow chart illustrating steps of a multiscale volume-based approach to retrospective topology correction according to the teachings of the invention.

Fig. 8 is a schematic diagram of an object illustrating the skeleton of a fill according to the teachings of the invention.

Fig. 9 is a schematic diagram of an object illustrating the thickness of a skeleton according to the teachings of the invention.

Fig. 10 is a schematic diagram of an object illustrating how simultaneously removing all residual points may create new handles in some rare situations according to the teachings of the invention.

DETAILED DESCRIPTION

Preferred embodiments of the present invention will now be described with reference to various drawing figures. These specific embodiments contain various details for the purposes of illustration only. Those skilled in the art will appreciate that the general principles of the present invention are not limited by these particulars, and many variations and alterations are possible and evident from these examples and associated teachings.

Notation and Definitions

As is customary in the art, a 3D binary image is defined by quadruple (Z ³ , n, n ^c , F), where Z ³ is the 3D cubic lattice representing all voxels (also called points) in the image, F C Z ³ represents the set of foreground voxels, n is an integer representing the adjacency in F ₅ and n ^c is an integer

representing the adjacency in F ^c = Z ³ - F, i.e., in the complement of F. Three types of adjacency commonly used are 6- adjacency, 18- adjacency, and 26-adjacency. Two voxels are 6-adjacent if they share a face, they are 18-adjacent if they share a face or an edge, and 26-adjacent if they share a face, an edge, or a corner. We use X generically to denote either F or F ^c . X is normally in m-adjacency, where m is 26 or 6. An n-neighbor of a point p is a point that is n-adjacent to p. The set of n-neighbors of a point p is denoted as N _n (p). Let N(p) denote N ₂₆ (P) U {p}. A set of points S C Z ³ is n-connected if S cannot be partitioned into two subsets that are not n-adjacent to each other. Topologically compatible adjacencies/connectivities of F and F ^c are (6, 26), (6, 18), (18, 6) and (26, 6).

The central concept in digital topology is the definition of simple point. A point in a binary image (Z ³ , n, n ^c , F) is simple if it can be added to or removed from F without changing the topology of both F and F ^c , i.e. without changing the number of connected components, cavities and handles of both F and F ^c . A simple point p can be characterized by two topological numbers T _n (p, F) and T _π o(p, F ^c ). In the case where m=26, the topological number T _m (x, X) may be defined to be the number of connected components in N ₂₆ (X) H X ; when m=6, we may define T _m (x, X) to be the number of connected components UiN _1S (X) H X that are adjacent to x. Intuitively, T _m (x, X) can be seen as the number of "ports" at which the point x is connected to the set of points X in m-adjacency. T _n c(p, F ^c ) — 1 can also be seen as the number of tunnels in the image ^λ N(p) = (N(p), n, n ^c , (N(p) O F) — {p}). A point p is called a simple point if and only if T _n (ρ, F) = 1 and T _n c(p, F ^c ) = l.

Another important concept is the multisimple point. A point p is multisimple relative to F (or F ^c ) if and only if it can be added to or removed from F (or F ^c ) without changing the number of handles and cavities of F (or F ^c ). Splitting and merging connected components in F (or F ^c ) are allowed. Let T ⁺ _n (p, F) and T ⁺ _n c(x, F ^c ), respectively, denote the number of foreground and background components in the whole image excluding p that are adjacent to p. Then p is multisimple relative to F if and only if T _n c(p, F ^c ) = 1 and T ⁺ _n (p, F) = T _n (p, F); p is multisimple relative to F ^c if and only if T _n (p, F) = 1 and T ⁺ _n c(p, F ^c ) = T _n (p, F ^c ).

Overview

A general outline of a method for reconstructing a topologically correct cortical surface model is shown in Fig. 7. In an initial step 700, a segmented volumetric image is obtained, e.g., from an

MRI device. The image is assumed to be segmented into white matter (W), grey matter (G) and background (B) voxels. Using 6-adjacency for G and B and 26-adjacency for W, the segmented image is then preprocessed if necessary in step 702 to ensure that W consists of a single connected components and contain no cavities.

Traditionally, volume-based methods of topology correction of the white matter are twofold in that there are two basic types of tunnel filling: filling the tunnels of the white matter or filling the tunnels of the complement of the white matter. Note that the second is. equivalent to cutting the handles of the white matter. Whenever a tunnel of an object (the white matter or its complement) is filled, the points used to fill the tunnel are always from the complement of the object.

Brain MRI segmentation, however, is usually able to separate the grey matter from the rest of the complement of the white matter and there is certain information provided by the prior segmentation that is not used by the traditional twofold methods. When a tunnel of the white matter is filled, there are generally three possibilities on the composition of the points used to fill the tunnel: the points are only from G; the points are only from B; the points are from both G and B. Considering the radiological property of Tl-weighted brain MRI (the average gray level of air, cerebrospinal fluid, grey matter, and white matter are in the ascending order) and the layered organization of W, G, and B regions, it is reasonable to assume that the points from B have less possibility of actually belonging to white matter than points from G. In other words, it is reasonable to prefer to use exclusively the points from G to fill the tunnels in W. Points in B are used to fill a tunnel in W only when the tunnel is passed through by one or more handles of B. Based on this rationale, a threefold topology correction method has been invented that involves three types of tunnel filling: filling the tunnels of the union of B and G using points from W; filling tunnels of W using points from G; filling tunnels of the union of G and W using points from B.

At each of a sequence of increasing scales (e.g., scales s=0, I ₅ 2, etc.), a topology correction filter is successively applied to cut handles and fill tunnels in W, as shown in step 704. The topology correction filter preferably comprises a morphological topology correction (MTC) filter that associates a larger correction cost to topology corrections (e.g., "cuts") at locations having a larger surface-like property (e.g., "flat" and "wide" regions). In other words, the technique

exploits the shape characteristic of surface-likeness of white matter as well as the layered organization of the W, G, and B regions.

In one embodiment, the following loop corrects the topology of W in a multiscale manner starting from scales of 0. At the end, W is homotopic to a ball. An object is homotopic to a ball if it can be transformed into a point through repeated removal of simple points.

1. Initialize scale s to s=0.

2. At each scale s: a) Apply MTC filter at scale s to fill tunnels of the union of G and B (i.e., the complement of W) using points from W. The points used to fill the tunnels are moved to G. Put it in other way, this step cuts the handles of W at scales s using points from G. The loop terminates if now W is homotopic to a ball. b) Apply MTC filter at scale s to fill the tunnels in the union of G and W (i.e. the complement of B) using points from B. The points used to fill the tunnels are moved to G. In other words, this step cuts the handles of B at scale s using points from G. c) Apply MTC filter at scale s to fill the tunnels of W using points from G. The points used to fill the tunnels are moved to W. d) Fill cavities (if any) in W. In extremely rare situations, step c) may create cavities in W. The loop terminates if now W is homotopic to a ball.

3. If W is not homotopic to a ball, increment s by 1 and go to step 2.

Whenever W is homotopic to a ball, the topology correction is completed, and in step 706 a cortical surface model is extracted from the topology corrected volumetric image, e.g., using a well-known topologically consistent marching cubes isosurface algorithm to generate the triangulated surface representation of the cortical surfaces. When 1 and 0 are used to respectively represent the foreground (in 26-adjacency) and the background (in 6-adjacency) object, an isovalue less than 0.25 should be used to avoid topological paradoxes.

Morphological topology correction (MTC) filter

In preferred embodiments, the topology correction filter seeks to preserve the surface-likeness of white matter. To preferentially preserve the shape of surface-like objects, the correction cost of handle cut B in Fig. IA is larger than that of cut A because the object is "wider" and more like a

surface at B than at A. Similarly, in the right object in Fig. IB the tunnel 100 is preferentially filled (i.e., the associate background handle cut) instead of cutting the foreground handle.

A typical application of the MTC filter is to fill the tunnels of object P using points from a set Q at a specific scale s. Here Q is a subset of the complement of P. Put it in another way, the MRC filter cuts the handles of the complement of P by moving a subset of points from Q to P. The adjacency of P and its complement should be compatible. Instances of P and Q and their adjacencies are:

1) P is the union of B and G and Q is W. P is in 6 adjacency. P ^c and Q are in 26 adjacency.

2) P is the union of G and W and Q is B. P is in 26 adjacency. P ^c and Q are in 6 adjacency.

3) P is the W and Q is G. P is in 26 adjacency and P ^c and Q are in 6 adjacency.

A minimal set of connected points f C Q used to fill a tunnel (or multiple tunnels simultaneously) is referred to as a fill. It can be also referred to as a cut in the handle-cutting perspective. In other words, a fill for the object P is a cut for the complement of P. A characteristic of a fill f is that it contains no multisimple points relative to P ^c — f (where P ^c represents the complement of P). A fill must be minimal in that removal of any proper subset of the fill will lead to filling of fewer or no tunnels.

The cost of filling tunnels with a fill f is measured as follows. First, let k be the skeleton off. A fill f is a 3D surface- like object and its skeleton is the loci of the centers of all maximal discs in f. For example, the skeleton of a fill is delineated with thick black lines in Fig. 8. In many other topology correction methods, the cost of a fill f corresponds to the radius of the largest maximal disc in f. A preferred embodiment of this invention uses a metric called diameter of a fill and takes it as the cost of the fill. A diameter of a fill is defined as the maximum sum of the distance between two skeleton ends and the radii of the discs centered at the two ends. For example, the diameter of the fill in Fig. 8 is the distance of the path between C and D passing through the skeleton between A and B. The middle point E in the path is referred to as the center of the fill. Since f is a subset of P° and removal of f cuts the correspondent handle in P ^c , the diameter of the fill f can be seen as a metric that measures the wideness of the P ^c where f locates.

To fill a tunnel, there may be many alternative fills. A preferred embodiment of the invention uses those with smaller diameters (i.e., less costs) and follows a multiscale approach (i.e., it first

uses fills with the smallest diameter to eliminate tunnels at the lowest cost and gradually eliminates tunnels with greater costs using fills with increasingly larger diameters).

To fill tunnels at scale s of P using points from Q ₅ MTC first performs iterative shape and topology preserving geodesic erosion on P ^c with respect to P ^c — Q to find centers of fills with diameters being 2s or 2s+l. After the centers of the fills at scale s (if any) are found, the fills are extracted from Q using three stages of geodesic dilations with topology control.

Iterative shape and topology preserving geodesic erosion (ISTPGE)

In a scenario of applying MTC filter to fill tunnels of P using points from Q, an operation of iterative shape and topology preserving geodesic erosion (ISTPGE) is first used to find the centers of the fills. Under certain constraints, ISTPGE iteratively erodes P ^c and the result curve skeleton points are centers of the candidate fills.

To clearly describe ISTPGE, we first define various point types (see Figs. 2A and 2B). First, we extend the definition of simple point to thicksimple point. A point p G X is a thicksimple point if p is a simple point and removal of p from X does not create new non-simple points in its 26- neighborhood N ₂ O(Jp)- For a digital object, the thickness of its skeleton can be 1 unit or 2 units. For example, the skeleton 902 of the object 900 in Fig. 9 is in 2 units thick. In this case, the points in the skeleton are still simple and the definition of thicksimple points are used to identify those points in the skeleton of 2 units thick (thick skeleton). If a point is not thicksimple, then it should be in the thick skeleton of the object.

We define- a border point p £ X as a point 6-adjacent to X ^c . Any points in X that are not border points are referred to as interior points of X . If a border point p G X is not a simple point, it should be in a topology-preserving surface skeleton of X . If T _ra c(p, X ^c ) is greater than 1, then p is referred to as a surface point; otherwise, p is referred to as a curve point. Note that in this context p is not a simple point, which means that either T _m (p, X)≠l or/and T _m c(p, X ² ) ≠l. Now that T _m c(p, X ⁰ ) ≤l and p is not a simple point, T _m (p, X)≠l must be true. So there are two cases: T _m (p, X)=O or T _1n (P, X)>1. The first case means that p is an isolated point, which can be seen as a degenerated curve segment. The second case means that there are more than one connected components in p's neighborhood. There are three subcases in the second case: p is actually a curve junction point when T _m (p, X)>2; p is a point between a curve and a surface edge; p is

really a curve point. We can refer to p as curve point without distinction for all these cases. Figs. 2A and 2B illustrate various point types in 26-adjacency (left object) and 6-adjacency (right object). The point types illustrated are curve points 1, thick curve points 2, surface points 3, thick surface points 4, thicksimple points 5, and curve end points 6.

If a border point p is a simple point but not a thicksimple point, then it should be in the object's topology-preserving thick surface skeleton. For example, all points in the skeleton in Fig. 9 are simple but not thicksimple. It is critical to distinguish thicksimple points from non-thicksimple points for identification of thick surface skeletons of X. If removal of a simple point p from X creates a new surface point in N ₂₆ (P) ₅ then we refer to p as a thick surface point; otherwise we refer to the non-thicksimple point p as a thick curve point. A thicksimple point p may be in one of the following situations:

1. At the boundary of a volume portion of X;

2. At the edge of a thin or thick surface portion of X 's surface skeleton;

3. In a thick curve portion of X's surface skeleton (only when m = 26);

4. At the end of a curve or a thick curve portion of X. In 26-adjacency, a point p is a curve end if it is thicksimple and all points in N ₂ 6(p) are 26-adjacent to each other. In 6-adjacency, a point p is a curve end if it is thicksimple and N ₆ (p) either contains exclusively curve or thick curve points, or contains one thicksimple point as well as curve or thick curve points.

When applied on a binary object X, a conventional morphological erosion technique removes all border points without preserving topology and shape of X. If we sequentially remove border points that are simple at the ^' moment of removal, the topology of X is preserved but the shape of X may still be significantly changed. For example, the object in Fig. 9 can be eroded to a single point by removing only simple points. In addition to checking whether a border point is simple to impose topology preservation, ISTPGE preserves object shape by thinning thick curve portions and eroding edge points of surface portions and boundary points of volume portions. At a specific scale s, ISTPGE is applied s iterations on X with respect to a mask set M, each which consists of following steps.

1. Identify border points and classify them into thicksimple points, thick surface points, thick curve points, surface points, curve points, and curve ends.

2. Remove thick curve points and thicksimple points that are not curve ends if they are simple at the moment of removal until there exists no such point. An additional condition to remove a point is that the point cannot belong to the mask set M.

Geodesic dilation with topology control

After ISTPGE is applied on P ^c with respect to P ^c -Q at a specific scale, the set of all curve points that are in Q, denoted R, is identified in the erosion result P ^Cf . (Note: P ^c ' is the result of applying ISTPGE on P ^c with respect to P ^c -Q) Let B ₀ = P ^C/ - R. We refer to R and B ₀ as residual set and body set respectively. Each residual point is the center of a candidate cut for cutting handles of P ^c (i.e., a candidate fill for filling tunnels of P). Removal of R from P ^c ' may decrease the number of handles in P ⁰ ', which is exactly desirable. However, it may also disconnect P ^c ', and R may be not minimal with respect to cutting handles in P ^c '. By several stages of geodesic dilation with topology control, we recover a subset of residual points and find the final cuts for cutting handles of P ^c (i.e., filling tunnels of P).

If removal of a residual point p from P ⁰ ' decreases the number of handles in P ^c \ p is referred to as a handle residual point; otherwise, p is referred to as a finger residual point. Simultaneously removing all residual points may actually create new handles in some rare situations. For example, in the object (in 6-adjacency) in Fig. 10, points 7 to point 30 are all identified as residual point at first. But removal of the residual point 7 and 30 simultaneously creates a new tunnel. For any residual point p, if T _m (p, P°-R)>1 (m is the adjacency of P ^c -R), then p is moved out of R into Bo. In this way, it is guaranteed that removal of all residual points won't create new handles. According to the mechanism of ISTPGE, the object P ^c is in the same scale of wideness at each location of handle residual points. Every handle residual point is associated with the same cost of handle cutting and the dilation procedure will prefer those in the middle of the connected components of handle residual points.

The following three dilation stages iterate in the same manner. The dilation in every stage involves a seed set S and a condition set C . In each iteration of every stage, any points in C that are adjacent to S are marked first and then are recovered if they satisfy some additional conditions. Each stage of dilation terminates when no more points can be recovered. Note that whenever a point is recovered, B ₀ and S are augmented with the recovered point. C is also

updated with the removal of the point. Also that after a point is removed from R, it is not referred to as a residual point any more.

In stage I ₅ S = B ₀ , C = R, and finger residual points are identified and recovered. In each iteration, a marked point p is identified as a finger residual point and recovered if T _m (p, R ⁺ ) ≤ 1, and p is adjacent to a body component Bi with d(B{) = 1 and is multisimple relative to Bo at the moment of recovery. Here a body component is a connected component in the body set B ₀ and d(Bi) denotes the degree of a body component Bi and is defined as d(Bi) = ^j =i N T _m (rj, Bi), where η represents a residual point and N represents the total number of residual points. d(Bi) can be seen as the number of ports at which Bj is connected to residual points. R ⁺ refers to the union of R and all body components with degree greater than 1.

In stage 2 (see Figs. 3 A-C), a minimal set of handle residual points R ^* C R is identified and any points in R - R ^* are recovered. The condition set C = R and the seed set S is initialized as the set of points in B ₀ that are not finger residual points identified in stage 1. In each iteration of stage 2, a marked point is recovered if it is multisimple relative to Bo at the moment of recovery. Figs. 3A-C give an illustration of stage 1 and stage 2. In stage 3, S = Bo, C = P ^c — Bo, and a marked point is recovered if it is multisimple relative to'Bo at the moment of recovery. Eventually, the original object P ^c is recovered except for the points in the cuts that cut the handles of P ^c (i.e., the fills that fill the tunnels of P).

Since a recovered point has to be multisimple, the number of handles in P ^c ' is not changed. Merging of body components may happen and is desirable. In rare situations, however, there might be more than one connected component in the result body set and we only keep the largest component in the case of using points from W to fill tunnels of the union of G and B

Figs. 4A-C demonstrate the behavior of the method on eliminating a handle in the white matter and a resulting cortical surface after topology correction. The handle in the white matter is removed by filling the associated tunnel 400. The method can detect that the white matter handle is wider than the associated tunnel (i.e., the gray matter handle), despite the fact that the former is thinner than the latter. Therefore, the algorithm chose to fill the tunnel instead of cutting the white matter so that the shape of the white matter is better maintained.