Mutual information based registration of 3D-image volumes on GPU using novel accelerated methods of histogram computation

Description

This application claims the benefit of filing of U.S. Provisional Patent Application Ser. No. 60/639,446, filed Dec. 23, 2004, entitled “Mutual Information Based Registration of 3d-image Volumes on GPU Using Novel Accelerated Methods of Histogram Computation,” and U.S. Provisional Patent Application Ser. No. 60/641,686, filed Jan. 6, 2005, and also entitled “Mutual Information Based Registration of 3d-image Volumes on GPU Using Novel Accelerated Methods of Histogram Computation,” the teachings of both of which are incorporated herein by reference.

BACKGROUND OF THE INVENTION

The invention relates to digital data processing and, particularly, to registering 3D image volumes. The invention has application, by way of non-limiting example, in medical imaging, e.g., in aligning diagnostic images acquired by computed tomography (CT) and/or magnetic resonance imaging (MRI). Other applications include microscopy and geosciences, to name but a few fields.

A common problem in three-dimensional (3D) imaging is, given two different 3D image volumes of the same object or of two different but corresponding objects, to compute a 3D transformation (e.g., rotation, translation, scaling) that brings the two images into alignment, so that after transformation, corresponding points in the two images are coincident or at least close to each other. This problem arises, for example, in aligning multiple 3D images (e.g., of a patient's head) generated by Computed Tomography (CT) and/or Magnetic Resonance Imaging (MRI).

A widely accepted standard procedure for solving this problem is to search in a space of transformations for the transformation that maximizes the “mutual information” of the two images. (See, by way of non-limiting example, Viola, “Alignment by Maximization of Mutual Information,” PhD thesis, Massachusetts Institute of Technology (MIT), Cambridge, Mass., March 1995, also referred to as A.I. Technical Report No. 1548, dated June 1995, of the MIT Artificial Intelligence Laboratory, the teachings of which are incorporated herein by reference).

Thus, given two digital images A and B with range [0 . . . N−1], i.e.

• • A: ³→[0 . . . N−1], B: ³→[0 . . . N−1]

then for a given transformation T: ³→³, the Mutual Information is defined by

H ⁡ ( A , B ∘ T ) = - ∑ a = 0 N ⁢ ⁢ ∑ b = 0 N ⁢ ⁢ p A , B ∘ T ⁡ ( a , b ) ⁢ log ⁡ ( p A , B ∘ T ⁡ ( a , b ) )

Here P_A,B∘T(a,b) denotes the probability that a pixel value a in image A coincides with a pixel value b in the transformed image B∘T. In order to compute this probability, essentially a 2D histogram has to be computed, namely:

p A , B ∘ T ⁡ ( a , b ) = card ⁡ ( { x ∈ G ❘ A ⁡ ( x ) = a ⋀ B ⁡ ( T ⁡ ( x ) ) = b } ) ⁢ 1 card ⁢ ⁢ ⁢ G

Here card denotes the number of elements in a set. The set G of locations x over which the histogram is computed may be chosen to consist of all pixels centres of the reference image A, or of any other possibly smaller set of points.

With this background, an object of this invention is to provide improved methods and apparatus for digital data processing.

A more particular object, by way of example, is to provide improved methods and apparatus for registering 3D image volumes.

Yet another more particular object, again, by way of example, is to provide improved such methods and apparatus as can be implemented at lower cost. And, a related aspect of the invention is to provide such methods and apparatus as can be implemented using standard, off-the-shelf components.

SUMMARY OF THE INVENTION

The foregoing are among the objects attained by the invention which provides, in some aspects, methods and apparatus for registering 3D image volumes by utilizing a graphics processing unit (GPU), to reduce computation time, e.g., in comparison to conventional implementations on a central processing unit (CPU).

Related aspects of the invention provide such methods and apparatus as utilize a GPU to facilitate maximizing the “mutual information” of the two images volumes being registered.

Yet still further related aspects of the invention provide such methods and apparatus as utilize the GPU to compute a histogram that facilitates determining a probability that a pixel value in one of the images being registered coincides with a pixel value in a transform of the other image.

Further aspects of the invention provide such methods and apparatus as utilize a two-pass approach to perform data scattering that supports histogram computation on the GPU. That approach includes the steps of, according to some aspects of the invention, (i) “drawing” a polygon covering a two-component render target having as many pixels as there are samples, (ii) interpreting that render target as a buffer of vertex positions and, for each sample, rendering into a new render target, a single vertex using coordinates from that buffer.

Related aspects of the invention provide such methods and apparatus in which the aforementioned two-step approach includes setting up the two-component render target with as many pixels as there are samples. The locations x of the samples can be stored in a texture residing on the GPU. This texture can have the same resolution as the render target. In addition, the images being aligned can be stored in (3D) textures also residing on the GPU.

A further related aspect of the invention provides such methods and apparatus in which the “drawing” step includes drawing a simple polygon covering the whole render target. For each pixel the location x of corresponding sample point is looked up. The intensities A(x) and B(T(x)) are then looked up and written into the render target.

A still further related aspect of the invention provides such methods and apparatus in which the interpreting step comprises (i) interpreting the two-component floating point render target as a buffer containing vertex positions, (ii) setting up the new render target to be of size N×N and initializing it to zeros, and (iii) for each sample, rendering a single vertex using the coordinates from the vertex buffer. Additive blending can be enabled and for each sample the corresponding output pixel is incremented by one

According to related aspects of the invention, after all samples have been rendered, the N×N render target represents the aforementioned histogram. This can be read back, for example, into main memory, and mutual information computed (e.g., on the CPU) as described above.

Still further aspects of the invention provide methods and apparatus as described above in which data sets that are too large to fit entirely into texture memory are broken down into smaller “bricks,” which are processed independently. According to related aspects of the invention, the resulting histograms are combined in a further step. According to still further related aspects of the invention, if more samples are taken than fit into a single texture, multiple passes are rendered, accumulating into the same histogram render target.

According to still other aspects of the invention, as an alternative to storing the sample locations in a texture, they are computed on-the-fly in the GPU's vertex or pixel shader.

Methods and apparatus according to the foregoing aspects of the invention can be used in connection with rigid or affine transformations, as well as for any other transformation(s) that can be evaluated on a GPU.

Methods and apparatus according to these aspects of the invention can parallel those described above except that, for example, in the drawing step, instead of rendering the intensities A(x) and B(T(x)) into a two-channel render target, the value A(x)*N+B(T(x)) is rendered into a single channel render target.

Methods and apparatus according to these aspects of the invention can further parallel those described above except that, by way of further example, in the interpreting step, a ID histogram is computed over the render target.

Methods and apparatus according to these aspects of the invention can still further parallel those described above except that, by way of still further example, the resultant ID histogram is interpreted as a 2D histogram using the following definition:

p A , B ∘ T ⁡ ( a , b ) = p A , B ∘ T 1 ⁢ D ⁡ ( a * N + b )

These and other aspects of the invention are evident in the text that follows.

BRIEF DESCRIPTION OF THE ILLUSTRATED EMBODIMENT

A more complete understanding of the invention may be attained by reference to the drawings, in which:

FIG. 1 depicts an environment in which methods according to the invention can be practiced, as well as apparatus and systems according thereto;

FIG. 2 is a flow chart depicting a method according to the invention for computing histograms used in determining the mutual information in two images.

DETAILED DESCRIPTION OF THE ILLUSTRATED EMBODIMENT

Described below are methods, apparatus and systems according to the invention that, inter alia, register three-dimensional (3D) image volumes by utilizing a graphics processing unit (GPU, also called “graphics card”) to facilitate maximizing the “mutual information” of the two images volumes and, thereby, reducing computation time necessary for such registration.

Apparatus and System Architecture

Such methods can be implemented in an imaging and digital processing environment of the type disclosed in commonly assigned U.S. Pat. No. 7,120,283, entitled Improved Methods and Apparatus for Back-Projection and Forward-Projection, the teachings of which are incorporated herein by reference. Referring to FIG. 1, that environment comprises a computer aided tomography system 10—although methods according to the invention can be implemented in other environments, as well.

The system 10 includes a conventional image acquisition apparatus 12 that generates projection images 14 of a volume 18 of an object 16. In the illustrated environment, this is accomplished in the conventional manner, e.g., by illuminating the object 16 with radiation from a source 22 and detecting with detector 20 (such as a charged-coupled device or other sensor array), radiation not absorbed by the object 16. Generally, multiple projection images are obtained, e.g., at different respective angles, in order that a three-dimensional representation of the object can be reconstructed, e.g., using the principles of computed tomography, computed tomosynthesis, and so forth.

Digital data processor 26 can be configured to perform that reconstruction (though such functionality is not critical to practice of the invention). In the illustrated embodiment, the digital data processor 26 forms part of a workstation that, additionally, registers pairs of such image volumes produced by such reconstruction. Alternatively, or in addition, digital data processor 26 can register, for example, an image volume produced by such reconstruction with an image volume generated by other medical imaging apparatus, e.g., other CT systems (not shown), magnetic resonance imaging systems (not shown), positron emission systems (not shown), to name a few. Still further, by way of example, digital data processor 26 can register pairs of image volumes generated by such other medical imaging apparatus.

More generally, although FIG. 1 shows a computer aided tomography system, it will be appreciated that digital data processor 26 and the techniques discussed herein are more broadly applicable to registering image volumes for a full range of medical and non-medical imaging applications and techniques, such as, by way of non-limiting example, geo-sciences, or microscopy, computer tomographs, magnetic resonance imaging devices, or confocal microscopes, to name just a few. Moreover, although digital data processor 26 is shown here has a stand-alone device, in other embodiments, it may be configured differently, e.g., as an embedded processor, or otherwise.

Regardless, it will be appreciated that methods as described herein, when implemented on such a digital data processor, form apparatus according to the invention. Further, when implemented on systems that include a digital data processor and other apparatus, e.g., image acquisition apparatus, such methods form systems according to the invention.

As shown in the drawing, the digital data processor 26 includes a central processing unit (CPU) 30, dynamic memory (RAM) 32, and I/O section 34, all of the type conventionally known the art. The digital data processor 26 may be coupled, via I/O section 34, with acquisition device 12, as well as with a monitor or other presentation device 28 on which the registered (and unregistered) images may be displayed. I/O section 34 can also provide coupling to other apparatus, e.g., other medical imaging systems, that generate image volumes requiring registration.

Illustrated digital data processor 26 also includes a graphical processing unit (GPU) 36 that is coupled to the CPU 30, through which it can access the other elements of the digital data processor 26, as shown. The GPU 36 includes pixel shaders 36a, vertex shaders 36b and textures 36c, all of the conventional type known in the art. The GPU serves, in the illustrated embodiment, as a coprocessor, operating with and under the control of the CPU 30 to facilitate registering three-dimensional (3D) image volumes in accord with the teachings hereof. Other embodiments of the invention employ multiple GPUs for this purpose, each responsible for a respective portion of processing described below. Still other embodiments use no GPU at all, relying on other types of processing elements (such as array processors, game engine processors, and so forth) to provide or supplement such processing, all in accord with the teachings hereof.

In embodiments which use GPUs, e.g., as described above, preferred such devices are of the variety having programmable vertex shaders and programmable pixel shaders and are commercially available in the marketplace from ATI Research (for example, the Radeon™ 9700 processor), Nvidia (for example, the GeForce™ FX and Quadro® processors) and/or otherwise compatible with vertex shader version 3.0, ARB_vertex_shader extension (for Approach 1 described below) and/or the ARB_IMAGING_EXTENSION (for Approach 2 described below).

Methods

Given two different 3D image volumes of the same object or of two different but corresponding objects, digital data processor 26 computes a 3D transformation (e.g. rotation, translation, scaling) that brings the two images into alignment, so that after transformation, corresponding points in the two images are coincident or at least close to each other. As noted above, those images may be from volume reconstructions performed by processor 26 itself and/or by other medical imaging apparatus. To this end, CPU 30 executes software that provides for such image volume alignment by searching a space of transformations for the transformation that maximizes the “mutual information” of the two images, for example, in accord with the general techniques disclosed by Viola, “Alignment by Maximization of Mutual Information,” supra (the teachings of which are incorporated herein by reference), as modified in accord with the teachings hereof. Thus, that software executes in accord with a general methodology that may be described as follows.

Given two digital images A and B with range [0 . . . N−1], i.e.,

• • A: ³→[0 . . . N−1], B: ³→[0 . . . N−1]

then for a given transformation T: ³→³, the Mutual Information is defined by

H ⁡ ( A , B ∘ T ) = - ∑ a = 0 N ⁢ ⁢ ∑ b = 0 N ⁢ ⁢ p A , B ∘ T ⁡ ( a , b ) ⁢ log ⁡ ( p A , B ∘ T ⁡ ( a , b ) )

Here p_A,B∘T(a,b) denotes the probability that a pixel value a in image A coincides with a pixel value b in the transformed image B∘T.

In order to compute that probability, digital data processor 26 computes, with the aid of GPU 36, a 2D histogram that can be expressed as follows:

p A , B ∘ T ⁡ ( a , b ) = card ⁡ ( { x ∈ G ❘ A ⁡ ( x ) = a ⋀ B ⁡ ( T ⁡ ( x ) ) = b } ) ⁢ 1 card ⁢ ⁢ ⁢ G

In the foregoing expressions,

• • G denotes a set of locations within the images over which the probability and the histogram are to be computed. The set G may be chosen to consist of all pixels (and specifically, for example, all pixel centres) of the reference image A, or of any other possibly smaller set of points.
  • card denotes the number of elements in the set,
  • N is the maximum pixel intensity value in the images A, B being aligned (and, more precisely, among the locations within set G). More generally, N can represent any other maximum “binning” value for histogram generation. Thus, in the illustrated embodiment, N can be 256 where A and B are full-range gray scale images. However, to minimize processing resources and/or time, lower values of N can be used as well, and
  • x denotes locations within the set.

In the illustrated embodiment, the aforementioned histogram is computed in accord with any of the approaches described below. These may be better understood by reference to the flowcharts of FIG. 2. In the discussion that follows, it will be appreciated that the trans-formation T (shown in the drawings as B∘T and in the text below as B(T(x)) is not limited to rigid or affine transformations, but could be any transformation that can be evaluated on a GPU. Moreover, it will be appreciated that, although the methods described below are applied to computation of a histogram for two images, one of which is transformed, those methods are equally applicable to computation of histograms for non-transformed images.

Approach 1

In order to efficiently compute the histogram, digital data processor 26 makes use of the capabilities of GPU 36 (which can be of the types indicated above or other modern GPU adaptable for use in accord with the teachings hereof). Using the highly parallelized pixel processing units native to the GPU, data necessary for histogram computation can be efficiently gathered via texture lookups. However, computing the histograms requires data scattering—i.e., for each pair of input values, the respective bin count of the histogram must be incremented. Since data scattering is not supported directly by typical GPU's, the two-pass approach shown in FIG. 2 and described below is utilized:

1. A two-component (or two-channel) floating point render target is set up with as many pixels as there are samples; see, step 40. The locations x of the samples are stored in a texture residing on the GPU; see, step 42. This texture has the same resolution as the render target. In addition, the images A and B to be aligned are stored in (3D) textures also residing on the GPU; see, step 44.

2. In a first pass, a simple polygon covering the whole render target is drawn; see, step 46. For each pixel in the render target, the location x of corresponding sample point is looked up; see, step 46a. Then the intensities A(x) and B(T(x)) are looked up and written into that pixel of the render target; see, step 46b.

3. In a second pass, the histogram is rendered; see, step 48. To do so, the two-component floating point render target is interpreted as a buffer containing vertex positions. In step 48a, a new render target of size N×N is setup and initialized with 0. For each sample a single vertex is rendered using the coordinates from the vertex buffer. Additive blending is enabled and for each sample the corresponding output pixel is incremented by one. See step 48b.

4. After all samples have been rendered, the N×N render target represents the histogram that was to be computed. The histogram is read back into main memory and mutual information is computed on the CPU according to the equations above. See step 50.

In computing the histogram by the foregoing approach, if the data sets are too large to fit entirely into texture memory, then they are broken down into smaller “bricks”, which are processed independently. Then the resulting histograms are combined in a further step. Also if there are more samples than would fit into a single texture, e.g., in steps 42 and 44, then multiple passes are rendered, accumulating into the same histogram render target.

As an alternative to storing the sample locations in a texture in step 42, they can be computed on-the-fly in the vertex or pixel shader, e.g., in step 46a. In case of a regular sampling, this can be more efficient.

In order to convert a render target into vertex data, in step 48, the illustrated embodiment utilizes functionality exposed, for example, in DirectX 9.0 graphics© 2002 Microsoft Corp with the vertex shader version 3.0 and in OpenGL2003 with the ARB_vertex_shader extension, or their equivalents.

Approach 2

This alternative approach does not require the vertex data functionality but instead a ID histogram functionality, as defined in the ARB_IMAGING_EXTENSION, which is part of OpenGL 1.2. However, since not all GPUs necessarily implement this extension in hardware, this approach is not necessarily efficient on all GPUs that support the extension.

Approach 2 can be implemented identically to Approach 1 (described above), with the following differences:

In step 40, a render target of only one channel is set up.

In step 46b, instead of rendering the intensities A(x) and B(T(x)) into a two-channel output, the value A(x)*N+B(T(x)) is rendered into a single channel render target.

In steps 48a-48b, a 1D histogram is computed over the render target. The size of the histogram must be N*N.

In step 50, the resulting render target—a 1D histogram—is the interpreted as a 2D histogram using the following definition:

p A , B ∘ T ⁡ ( a , b ) = p A , B ∘ T 1 ⁢ D ⁡ ( a * N + b )

A further understanding of the foregoing may be attained by reference to the following, the teachings of which are incorporated herein by reference, and the citation of which is not intended as, nor shall it be deemed, an admission that any or all of them is/are prior art: P. A. Viola. Alignment by Maximization of Mutual Information. PhD thesis, Massachusetts Institute of Technology (MIT), Cambridge, Mass., March 1995. http://citeseer.ist.psu.edu/viola95alignment.html; MICROSOFT CORP. 2002. DirectX 9.0 graphics, Dec. Available at http://msdn.microsoft.com/directx; OpenGL2003: OpenGL ARB; OpenGL Extension ARB_vertex_shader, 2003, http://oss.sgi.com/projects/ogl-sample/registry/ARB/vertex_shader.txt.

Described above are methods and apparatus meeting the objects set forth herein. Those skilled in the art will appreciate that the embodiments described above are merely examples of the invention and that other embodiments, incorporating changes thereto, fall within the scope of the invention. Thus, by way of non-limiting example, it will be appreciated that the methods described herein can be used to calculate histograms not only for use in determining mutual information, but also for determining entropy and related expressions. In view of the foregoing, what we claim is: