JP5132559B2 - Digital image segmentation method and computer-readable program storage device - Google Patents

Description

US Patent Application No. 60 / 764,122, "Automatic segmentation of vessels in breast MR sequences as a false positive elimination technique for automatic lesion" by Hermosillo et al., Filed Feb. 1, 2006. detection and segmentation using the shape tensor "and US Provisional Application No. 60 / 704,930," Method for automatic extraction of image structure based on the second order geometric moment "by Hermosillo et al. Both of which are hereby incorporated by reference.

The present invention is directed to segmentation of digital medical images. In particular, the present invention relates to a digital image segmentation method and a computer readable program storage device for performing the method.

  A contrast agent applied MR sequence is a high performance diagnostic tool for detecting lesions in the breast. In general, this diagnosis first identifies suspicious enhancement regions in the acquired image after application of the contrast agent (post-contrast) as compared to the acquired image before application of the contrast agent (pre-contrast). Therefore, to automate this process requires the execution of a computer-aided detection system. The problem with this system is that in the image after contrast agent application, some unsuspecting structure shadows are created in addition to the lesion. Most of these structures are blood vessels. A blood vessel is the main type of false positive structure that appears when a lesion is automatically detected as an area to be emphasized after contrast agent injection.

  Dynamic subtraction of TI-weighted images after contrast agent application is usually performed as part of a protocol for examining breast lesions by nuclear magnetic resonance imaging (MRI). Because the lesion usually contains a high vascular distribution, the perfusion of the contrast agent makes the lesion appear brighter than the background, so this physiotherapy is extremely sensitive. By automatically segmenting the lesion, the radiologist can make accurate automatic measurements, and such measurements are consistent for all viewers. If the vessel leading to the lesion does not even segment to the vessel, a region-enhanced segmentation algorithm or simpler thresholding can be used to segment such a lesion. Thus, segmentation can be facilitated by removing blood vessels. On the other hand, in order to automatically detect a lesion, it is necessary to have an ability to distinguish between a lesion and various types of normal structures that are also emphasized by a contrast medium. Such structures include breast parenchyma, blood vessels, areas under the nipple, and areas surrounding the heart. The development of automated methods for segmenting vascular structures in physical therapy such as CT and MR angiography has attracted attention. There is a great deal of literature on this subject describing both automatic and semi-automatic methods, including a wide variety of models, assumptions and techniques. In the context of a clinical workflow, the extraction of vascular structures must be fully automatic and require only a little computation time. One technique that works well, can be easily verified with clinical data, and is easy to implement requires the use of moments that are rarely reported in the research literature. Previous methods based on moments use moment invariant to extract and characterize blood vessels in infrared images of laser-heated skin, extract vessel structures from large CT datasets and characterize blood vessels There are multi-resolution computational moment filters for using geometric moments and for extracting linear structures from extremely noisy two-dimensional images.

  The use of geometric moments to extract the image structure depends on the methods proposed in the literature. In many cases, the moment of inertia is calculated based on a binary image obtained after thresholding. A problem associated with this is that the threshold is generally difficult to select, and if the threshold is low, a thin blood vessel that tends to have low brightness is assimilated with a nearby structure, so that a thin blood vessel cannot be detected. Another problem associated with thresholding is that the structure is “pixelated”, ie, a sharp edge is generated, which causes the calculation of the shape of the structure relative to the true shape of the structure under the edge. Is inaccurate.

  An alternative to thresholding is to calculate the moment using the image brightness function f as the density function. However, in regions where the signal-to-noise (SN) ratio is low, it becomes difficult to set a threshold for the eccentricity of the ellipse that is applied to detect elongated structures. For example, FIG. 1A shows a MIP of a subvolume extracted from a real image around a blood vessel branch. The top row shows the original voxel values using nearest-neighbor interpolation. The middle row shows the binary image obtained after manual thresholding. The threshold is adjusted to acquire both vessels, which is a fairly difficult task to do automatically. The impact of thresholding pixelation is obvious, which affects the accuracy of shape descriptors. The third row shows the same subvolume using a more advanced interpolation scheme.

  The exemplary embodiments of the invention described herein generally provide a high intensity tubular structure for automatic segmentation of blood vessels in a chest MR sequence based on geometric moments for extracting the tubular structure from an image. And a method and system for automatically detecting its function. The method according to an embodiment of the present invention eliminates thresholding of the image by reliably restoring the structure with a very low signal-to-noise (SN) ratio based on the eigenvalue of the shape tensor. The method according to an embodiment of the present invention does not depend on the first order image derivative, such as a method based on the eigenvalue of the average structure tensor, or the second order image derivative, such as a method based on the Hessian eigenvalue, and the Hessian or An output smoothing circuit inherent in the structure tensor based approach is avoided. The method according to an embodiment of the invention can be performed quickly and only takes a few seconds per sequence. Examining the results based on the motion compensated chest MR sequence, the method according to one embodiment of the present invention performs a method that reliably segments the blood vessel while leaving the lesion intact, and that differs in both sensitivity and localization accuracy, It can be seen that there is not much influence from the scale selection parameter.

According to the present invention, providing a digital image including a plurality of intensities corresponding to a plurality of point domains on a three-dimensional grid; and a shape matrix of one selected point in the image. a step of defining the luminance moment in the window of a plurality of points around the selected points, and calculating the eigenvalues of the shape matrix, the eccentricity of the structure underlying the points from the eigenvalues Determining the rate; segmenting the image based on the eccentricity value; defining the shape matrix; calculating eigenvalues of the shape matrix; and eccentricity of the substructure determining the rate may be repeated for all points in the image, wherein the selected point, defined in advance An intermediate enhancement greater than the value, the contrast agent, segmentation method for a digital image that is applied to the contents of said digital image before the acquisition of the image is provided.

  According to yet another aspect of the present invention, the selected point has a median enhancement that is greater than a predetermined threshold, wherein the contrast agent is the digital image prior to acquisition of the image. Applied to the content.

  According to still another aspect of the present invention, the intermediate enhancement determines a difference between an intermediate value of the image to which the contrast agent is applied and an intermediate value of the image before the contrast agent is applied, and determines the difference in advance. It is calculated by normalizing to be within the specified range.

  According to yet another aspect of the present invention, the shape matrix Sα is defined as follows.

here,

The moments m _{p, q, r,} α are defined as follows:

w is a window function with compact support p, q, rμ0 and αμ1.

  According to yet another aspect of the invention, the integral is calculated by the sum of the entire finite neighborhood around each point.

  According to yet another aspect of the invention, the window function is defined as follows.

Here, v _x , v _y , and v _z are image point intervals, N _x , N _y , and N _z are non-negative defined integers, and the window size is configured by the maximum diameter of interest.

  According to yet another aspect of the invention, the method includes calculating the moment using nearest neighbor interpolation and correcting the shape matrix according to the following equation:

v _x , v _y , v _z are image point intervals.

  According to yet another aspect of the invention, the method includes calculating the moment using trilinear interpolation.

  According to yet another aspect of the invention, α = 1 and the shape matrix is corrected according to the following equation:

v _x , v _y , v _z are image point intervals.

In accordance with the present invention, there is provided a computer readable program storage device tangibly recorded with a program of instructions readable by a computer and executable by a computer to perform the method steps for segmenting a digital image. .
That is, according to the present invention, there is provided a computer readable program storage device recording a program of instructions readable by a computer and executable by a computer to perform the method steps for segmenting a digital image. Providing a digital image including a plurality of intensities corresponding to a plurality of point domains on a three-dimensional grid; and a shape matrix of one selected point in the image around the selected point Defining from moments of brightness within a window of a plurality of points; calculating an eigenvalue of the shape matrix; determining an eccentricity of a structure underlying the point from the eigenvalue; and the eccentricity Segmenting the image based on values; Including the step of defining a shape matrix, calculating eigenvalues of the shape matrix, and determining an eccentricity of a substructure are repeated for all points in the image and the selected point Is proposed a computer readable program storage device having an intermediate enhancement greater than a predetermined threshold and in which a contrast agent is applied to the contents of the digital image prior to acquisition of the image .

  Exemplary embodiments of the invention described herein generally include systems and methods for automatically detecting high brightness tubular structures and automatically segmenting blood vessels in a chest MR sequence. The method according to an embodiment of the invention is based on the eigenvalues of a shape tensor. This method can be compared to a method based on the mean Hessian eigenvalue and a method based on the mean structure tensor eigenvalue. A Hessian is defined from a second derivative and can be regarded as a structure descriptor of order 2. Similarly, a structure tensor is a degree 1 structure descriptor. The shape tensor can be regarded as a structure descriptor of degree 0.

As used herein, the term “image” refers to multidimensional data composed of discrete image elements (eg, pixels of a two-dimensional image and voxels of a three-dimensional image). Image is, for example, computer tomography, nuclear magnetic resonance imaging, an ultrasonic or subject of the medical image acquired by another image diagnostic systems known to those skilled in the art. Images may also be provided from non-medical situations such as remote sensing systems and electron microscopes. Images can be regarded from R ³ and functions to ^{R (function from R 3 to R} ) , the methods of the present invention is not limited to such images, any dimension image (e.g., two-dimensional image Or three-dimensional volume). For a two-dimensional or three-dimensional image, the domain of the image is typically a two-dimensional or three-dimensional rectangular array, and each pixel or voxel is addressed relative to a set of two or three mutually orthogonal axes. Can be specified. The terms “digital” and “digitizing”, as used herein, refer to an image or volume in a digital or digitized format obtained by conversion from a digital image acquisition system or analog image, as appropriate.

A method according to an embodiment of the present invention processes image luminance by calculating a second-order geometric moment of a base (high luminance) structure. Although the method can be applied to a binary image obtained by applying a threshold value to an image after the initial application of a contrast agent, the method can be applied without this threshold value. The eigenvalue of the secondary geometric moment is a classic tool for characterizing shapes in object recognition. However, such an eigenvalue has not been used as a filter for extracting an image structure until now. Given a binary image, a small subvolume around each pixel (whose size is associated with a given structure) is considered, and the shape tensor is in that location relative to the center of the subvolume. It is defined as the second moment of the position of the high brightness voxel. For voxels whose center pixel is bright and close enough to the center of the base shape, the eigenvalues of the shape tensor are calculated and the filter response is assigned the value λ ₁ -λ ₂ / (λ ₁ + λ ₂ ). Here, λ ₂ > λ ₁ is the maximum eigenvalue.

  According to one embodiment of the present invention, a geometric three-dimensional moment can be defined as follows:

  Where w is a positive and symmetric window function with a compact support that provides localization p, q, rμ0 and αμ1. The shape tensor of order α is defined as follows from the viewpoint of these moments.

here,

This matrix is symmetric, so all its eigenvalues are real. Assuming that the _three eigenvalues are λ ₃ > λ ₂ > λ ₁ μ0, the filter response can be defined as follows.

In the case of a linear or cylindrical structure such as a blood vessel, C _line λ1.

According to an embodiment of the invention, the eccentricity of the basic shape is calculated based on the eigenvalue of Sα 0 [λ ₁ [λ ₂ [λ _3] . Here, α >> 1. The greater α is, the greater importance is given to the higher luminance value, which acts almost like threshold processing. When the value of α is high, the S / N ratio becomes extremely low as shown in the simulated experiment of FIG. 1 (b), and the synthetic tubular structure to which uniform noise is added becomes a conventional matrix of inertia (matrix of inertia). ) And α = 15 shape tensor.

FIG. 1B shows a blood vessel simulated with a standard moment of inertia without threshold processing and a shape tensor of α = 15, and its detection. The columns are, from left to right, (1) the central slice of the original composite volume, (2) its maximum intensity projection image (MIP), (3) the blood vessel MIP with the blood vessel removed by the standard moment method, (4 ) Shows the MIP of the blood vessel detected by the moment method, (5) the MIP of the volume from which the blood vessel has been removed using the shape tensor of α = 15, and (6) the MIP of the blood vessel detected using the shape tensor of α = 15. . The six rows represent progressively higher levels of additive uniform noise, from top to bottom (1) 56.3, (2) 36.7, (3) 20.4, (4) 11.6, ( 5) The SN ratio of 5.5 and (6) 0.8 dB is shown. The shape eccentricity threshold is the same for the entire row for each algorithm. In all cases, the detection criterion, if the inertia matrix corresponds to S ₁₅ is _{(λ 3 / λ 2)>} 15, when corresponding to S ₁ was _{(λ 3 / λ 2)>} 2. This improvement in detection performance was seen in actual cases.

  In practice, the above integral is generally replaced by a sum over a finite neighborhood around each voxel, since only f is known at the voxel location. For all experiments, it can be assumed that the localization function is given by:

Where v _x , v _y , v _z are voxel intervals of the image, and N _x , N _y , N _z are defined integers that are not negative such that the window size is composed of the maximum diameter of interest. is there. Next, for a given image, consider a small subvolume around each pixel and define:

Here, ρ _ijk is the value of the image at the voxel corresponding to the index i, j, k.

The eigenvalues 0 [λ ₁ [λ ₂ [λ ₃ ] of the above matrix are calculated, where the values ^ μ... Are the same as described above but using the total moment. The eccentricity or elongation of the substructure is determined according to a typical eccentricity criterion ε = (λ ₃ −λ ₂ ) / (λ ₃ + λ ₂ ) taking a value of 0 to 1 or when λ ₂ > 0. It can be measured simply by the ratio λ ₃ / λ ₂ .

  Since the moment-based method does not assume the differentiability of the image brightness function f, a simple interpolation method such as nearest neighbor or tri-linear interpolation is used instead of the sum of voxel values. The integral of the interpolation function can be calculated. In particular, in the case of trilinear interpolation, higher accuracy can be expected by using these integration values.

  In the case of nearest neighbor interpolation, it can be seen that the following equation must be used to replace the matrix ^ Sα with the following equation:

  In the case of trilinear interpolation, the function f is then given by:

Here, i, j, and k are indices of the image voxel, ρ _ijk is an image value in the voxel, and g _ijk is as follows.

  Next, it can be described as follows.

  Therefore, in the trilinear interpolation when α = 1, the matrix Sα should be replaced by the following equation.

The situation becomes more complicated for general shape tensors using trilinear interpolation (α> 1), where fα is given by (Σ _{i, j, k} ρ _ijk g _ijk ) α.

The corresponding integral can be calculated in closed form, but is considerably more complicated. Note that according to one embodiment of the present invention, the calculation of the moment of gα _ijk is no longer valid for calculating the corresponding shape tensor, as in the case of α = 1 described above. Initially, the aforementioned moment can be determined using a method that is not very straightforward but can be generalized to α> 1. This can be done in a one-dimensional case, and the two-dimensional and three-dimensional cases are simply generalizations of it. Assuming that ρ _k = 0 for k <i−2 or k> i + 2, the following equation is obtained.

The above four integrals can be obtained from the three piecewise-linear basis functions shown in FIG. Referring to the figure, the first basis function g _i-1 is defined from domain (i-2) v _x to iv _x , and the second basis function g _i is defined from domain (i-1) v _x ( i + 1) v _x and the third function g _{i + 1} is defined from domain iv _x to (i + 2) v _x .

  This integral calculation method can be generalized to α> 1. For example, it can be calculated as follows.

Similarly,

Also

  In the three-dimensional case, it may be possible to find a general expression when α> 1, but the resulting polynomial is quite complex for accuracy improvement. In the two-dimensional case, the above four integrals are 16 integrals, and in three dimensions, there are 64 integrals.

  The method according to one embodiment of the present invention was tested with over 100 motion corrected breast MR dynamic sequences. From the obtained results, it can be seen that the blood vessels can be reliably segmented while leaving the lesion intact. According to one embodiment of the invention, moments are calculated over a fixed size sliding window, but only points where the intermediate enhancement is above a predetermined threshold are considered. This threshold can be selected low enough to detect even small blood vessels. The calculation does not depend on the threshold value, so it is not difficult to set. The intermediate enhancement is calculated by subtracting the value of the acquired image before applying the contrast agent from the intermediate value of the acquired image after applying the contrast agent in each image voxel. This difference is then normalized by applying an affine function so that the resulting enhancement is in the range [0,200]. 3 (a)-(c), 4 (a)-(c) and FIG. 5 show some representative examples of the results.

  3 (a)-(c) show segmentation of large lesions, and FIGS. 4 (a)-(c) show segmentation of multiple small lesions. In both of these figures, panel (a) shows the image after thresholding the first contrast agent application, panel (b) shows the detected blood vessel and panel (c) shows the blood vessel removed. Showing lesions.

FIG. 5 shows segmentation of the vasculature in chest MRI using a shape tensor with α = 6. Three columns show orthogonal views of the same patient. The first line shows the original MIP with intermediate enhancement. The second row shows the same volume with the blood vessels automatically removed. The third line shows the MIP of only the removed blood vessel. Note that completely different diameters and enhancement levels of blood vessels are accurately segmented. Detection was performed by taking a place where the eigenvalue of the shape tensor was (λ ₃ / λ ₂ )> 3.

  Note that in each of these figures, the thin blood vessels are also accurately segmented and the small spherical structures are left intact. As further verification, the method according to one embodiment of the present invention extracted vascular structures in 40 cases investigated by three radiologists who marked a total of 75 lesions. In all cases, the blood vessels were correctly segmented and all marked lesions remained intact.

  FIG. 6 shows a flowchart of a method for moment-based segmentation according to an embodiment of the present invention. Referring to the figure, in step 61, an image to be segmented is provided. According to the determination of step 62, the shape tensor of the voxel in the image where the intermediate contrast enhancement exceeds a predefined threshold is calculated. At step 63, the moments used to define the shape tensor are calculated over a fixed size window around the selected voxel. In step 64, the eigenvalue of the shape tensor is calculated, and in step 65, the eccentricity of the substructure is determined. The process loops at step 66 until all voxels have been processed. At step 67, the image is segmented based on the eccentricity value obtained from the shape tensor.

  The moment-based method of extracting local shape information can be compared with a method based on higher order image derivatives. For example, Gradient Square Tensor (GST), or structural tensor, has been proposed as a robust method for evaluating local structural dimensions. This is based on the first derivative and is therefore sometimes referred to as a degree 1 structure descriptor. The Hessian eigenvalue also provides local image structure information and the principal curvature of the isolevel at a given point. The Hessian and principal curvature are defined from the second derivative and are therefore sometimes referred to as degree 2 structure descriptors. The shape tensor can be regarded as a structure descriptor of degree 0. It is based on integration and thus has the property that it is extremely resistant to noise compared to methods based on either the first or second derivative. Furthermore, the image function need not assume differentiability that simplifies modeling. The problem with the shape tensor based method is that the joint is not detected. A better understanding is also needed to determine if the geometric shape characteristics can be calculated from the eigenvalues of the shape tensor with α> 1.

  It should be understood that the present invention can be implemented in various forms of hardware, software, firmware, dedicated processes, or combinations thereof. In one embodiment, the present invention can be implemented in software as a tangible application program embodied on a computer readable program storage device. The application program may be uploaded and executed on a device that includes any suitable architecture.

  FIG. 7 is a block diagram of an exemplary computer system for implementing a moment-based segmentation method according to one embodiment of the invention. Referring to FIG. 7, a computer system 71 for implementing the present invention may include a central processing unit (CPU) 72, a memory 73, and an input / output (I / O) interface 74, among others. The computer system 71 is generally connected to a display device 75 and various input devices 76 such as a mouse and a keyboard via an input / output interface 74. The support circuit can include circuits such as a cache, a power supply, a clock circuit, and a communication bus. Memory 73 may include random access memory (RAM), read only memory (ROM), disk drive, tape drive, etc., or a combination thereof. The present invention can be implemented as a routine 77 that is stored in the memory 73 and executed by the CPU 72 to process signals from the signal source 78. Accordingly, the computer system 71 is a general-purpose computer system that becomes a dedicated computer system when executing the routine 77 of the present invention.

  The computer system 71 also includes an operating system and microinstruction code. The various processes and functions described herein may be part of microinstruction code or part of an application program (or a combination thereof) that is executed via an operating system. In addition, various other peripheral devices may be connected to the computer platform such as an additional data storage device and a printing device.

  Further, since some of the system components and method steps shown in the accompanying drawings can be implemented in software, the actual connections between the system components (or process steps) can depend on the manner in which the invention is programmed. Please understand that it may be different. From the teachings of the invention provided herein, one of ordinary skill in the related art will be able to review these and similar embodiments or configurations of the invention.

  Although the invention has been described in detail in the preferred embodiments, those skilled in the art will recognize that various modifications and changes may be made to these embodiments without departing from the spirit and scope of the invention as set forth in the appended claims. It should be understood that alternatives can be made.

(A) is a diagram showing the MIP of a subvolume extracted from a real image around a blood vessel branch according to an embodiment of the present invention, and (b) is a moment of inertia without threshold processing according to an embodiment of the present invention. Of a simulated blood vessel and its detection FIG. 6 shows basis functions used for 1D linear interpolation according to an embodiment of the present invention. Diagram showing segmentation of large lesions according to one embodiment of the invention FIG. 4 illustrates segmentation of multiple small lesions according to one embodiment of the present invention. FIG. 4 illustrates segmentation of vascular structures in chest MRI using shape tensors according to one embodiment of the present invention. Flowchart of a method for moment-based segmentation according to an embodiment of the invention 1 is a block diagram of an exemplary computer system that implements a method of moment-based segmentation according to one embodiment of the invention.

Explanation of symbols

75 Display device 76 Input device 77 Memory 78 Signal source

Claims (18)

Providing a digital image including a plurality of intensities corresponding to a plurality of point domains on a three-dimensional grid;
Defining a shape matrix of one selected point in the image from a moment of intensity in a window of points around the selected point;
Calculating eigenvalues of the shape matrix;
Determining the eccentricity of the structure underlying the point from the eigenvalue;
Segmenting the image based on the eccentricity value;
A step of defining the shape matrix, calculating the eigenvalues of the shape matrix, determining the eccentricity of the substructure is, is repeated for all points in the image,
The selected point has an intermediate enhancement greater than a predetermined threshold and a contrast agent is applied to the contents of the digital image prior to acquisition of the image Method. The intermediate enhancement obtains a difference between an intermediate value of an image to which the contrast agent is applied and an intermediate value of an image before the application of the contrast agent, and normalizes the difference to be within a predetermined range. the method of claim 1 which is calculated by. The shape matrix Sα is defined as follows:

here,

And moments m _{p, q, r,} α are defined as follows:

2. A method according to claim 1, wherein w is a window function with compact stands p, q, r [mu] 0 and [alpha] [mu] 1. 4. The method of claim 3 , wherein the integral is calculated by a sum over a finite neighborhood around each point. The window function is defined by

v _x, v _y, v _z is the image point _{spacing, N x, N y, N} z is an integer defined nonnegative The method of claim 3 window size including the maximum diameter interest. Calculating the moment using nearest neighbor interpolation and correcting the shape matrix according to the following equation:

v _x, v _y, v _z The method of claim 3 is the image point spacing. The method of claim 3 , further comprising calculating the moment using trilinear interpolation. α = 1, and correcting the shape matrix according to the following equation:

v _x, v _y, v _z A method according to claim 7 is an image point spacing. Providing a digital image including a plurality of intensities corresponding to a plurality of point domains on a three-dimensional grid;
Defining a shape matrix of one selected point in the image from a moment of intensity in a window of points around the selected point;
Calculating eigenvalues of the shape matrix;
Determining the eccentricity of the structure underlying the point from the eigenvalue,
The shape matrix Sα is defined as follows:

here,

The moment m _{p, q, r,} α is defined as:

w is a compact table p, q, Ri window function der with rμ0 and αμ1,
The selected point has an intermediate enhancement greater than a predetermined threshold and a contrast agent is applied to the contents of the digital image prior to acquisition of the image Method. Defining a shape matrix; calculating eigenvalues of the shape matrix; determining an eccentricity of a substructure of all points in the image; and segmenting the image based on the eccentricity value The method of claim 9, further comprising repeating the step of converting. To perform the method steps for segmenting a digital image, a computer readable program storage device which records a program of instructions executable by read-possible and and computer by a computer,
Providing a digital image including a plurality of intensities corresponding to a plurality of point domains on a three-dimensional grid;
Defining a shape matrix of one selected point in the image from a moment of intensity in a window of points around the selected point;
Calculating eigenvalues of the shape matrix;
Determining the eccentricity of the structure underlying the point from the eigenvalue;
Segmenting the image based on the eccentricity value;
A step of defining the shape matrix, calculating the eigenvalues of the shape matrix, the step of determining the eccentricity of the substructure, is repeated for all points in the image,
The selected point has an intermediate enhancement greater than a predetermined threshold and a contrast agent is applied to the content of the digital image prior to acquisition of the image
Computer readable program storage device. The intermediate enhancement is calculated by obtaining a difference between an intermediate value of the image to which the contrast agent is applied and an intermediate value of the image before the contrast agent is applied, and normalizing the difference within a predetermined range. A computer readable program storage device according to claim 11 . The shape matrix Sα is defined as follows:

here

The moment m _{p, q, r,} α is defined as:

12. A computer readable program storage device according to claim 11 , wherein w is a window function having compact stands p, q, r [mu] 0 and [alpha] [mu] 1. 14. The computer readable program storage device of claim 13 , wherein the integral is calculated by a sum over a finite neighborhood around each point. The window function is defined by the following equation:

v _x, v _y, v _z is the image point _{spacing, N x, N y, N} z is an integer defined nonnegative, window size according to claim 13, which is constituted by the maximum diameter interest Computer readable program storage device. The method further includes calculating the moment using nearest neighbor interpolation and correcting the shape matrix according to the following equation:

v _x, v _y, v _z is a computer readable program storage device of claim 13 is an image point spacing. The computer readable program storage device of claim 13 , further comprising calculating the moment using trilinear interpolation. α = 1, and the shape matrix is corrected according to the following equation:

v _x, v _y, v _z is a computer readable program storage device of claim 17 is an image point spacing.

Images