CN101278316A - System and method for automatic segmentation of vessels in breast MR sequences - Google Patents

Abstract

A method for segmenting digitized images includes providing (61) a digitized image, selecting (62) a point with a median enhancement greater than a predefined threshold, wherein a contrast enhancing agent was applied to the subject matter of said digitized image prior to acquisition of said image, defining (63) a shape matrix for the selected point in said image from moments of the intensities in a window of points about said selected point, calculating (64) eigenvalues of said shape matrix, determining (65) an eccentricity of a structure underlying said point from said eigenvalues, and segmenting (67) said image based on said eccentricity values, wherein the steps of defining a shape matrix, calculating eigenvalues of said shape matrix, and determining the eccentricity of the underlying structure are repeated (66) for all points in said image.

Description

Automatically the system and method for cutting apart the vascular in the breast MR sequences

The cross reference of relevant U. S. application

The application requires people's such as Hermosillo, on February 1st, 2006 submitted to, sequence number is 60/764, people's such as 122 U.S. Provisional Application " Automatic segmentation of vesselsin breast MR sequences as a false positive elimination techniquefor automatic lesion detection and segmentation using the shapetensor " and Hermosillo, on August 2nd, 2005 submitted to, sequence number is 60/704, the right of priority of 930 U.S. Provisional Application " Method for automatic extraction ofimage structure based on the second order geometric moment ", the content of these two U.S. Provisional Applications is incorporated herein by reference.

Technical field

The present invention be directed to cutting apart of digitized medical image.

Background technology

The MR sequence that contrast strengthens is the strong diagnostic tool that is used for detecting the damage of breast.Usually, strengthen the suspicious enhancing zone that obtains in the thing with respect to (pre-contrast) before strengthening in back (post contrast) the acquisition thing by identification and begin diagnosis.Therefore, making this process automation is the action that computer-aided detection system need be carried out.The difficulty of this system is the following fact, and promptly except damage, many non-suspect structure also strengthen in the image after enhancing.Major part in these structures is a vascular.The main type of the false positive structure that vascular occurs when being be enhanced regional after automatically damage check being injection of contrast medium.

The shadow that dynamically subtracts that strengthens the back t1 weighted image is performed as the part of agreement routinely and utilizes Magnetic resonance imaging (MRI) assessment injury of breast.Because damage comprises high vascular distribution usually, so the perfusion of contrast preparation makes damage seem brighter than background, and therefore this mode is quite responsive.Automatically cutting apart damage can provide accurate measurement result automatically for the radiologist, and makes these measurement results become more consistent between the reader.If not cutting apart the fact that penetrates vascular because the vascular that is attached to those damages causes, then the region growing partitioning algorithm or even simple thresholding (thresholding) can be used to cut apart those damages.Therefore, remove vascular and can make that to cut apart task easy.On the other hand, the automatic detection requirement of damage will damage the ability that makes a distinction with various types of normal configuration of also utilizing contrast preparation to strengthen.These normal configuration comprise parenchymal tissue, vascular, inframammillary zone and the heart zone on every side of breast.It is interesting the automated process etc. of blood vessel structure to be cut apart in research and development in the mode of similar CT or MR angiography etc.Document about this theme is very abundant, these documents have not only been described automated process but also described semi-automatic method, these methods covered very wide scope model, suppose and technology.In the clinical workflow environment, the extraction blood vessel structure should be fully automatically and require the only computing time in several seconds.A kind of carry out good, can easily utilize clinical data to verify and the technology implemented easily comprises the use square, in the research document, this is seldom reported.The previous method based on square comprises the vascular in the infrared image of the skin that uses invariant moments (moment invariant) to extract and characterize LASER HEATING, use geometric moment to come to extract blood vessel structure, and characterize this vascular and calculate multiresolution (multi-resolution) the square wave filter that is used for extracting linear structure from the big 2D image of noise from big CT data centralization.

It is different among the method that document proposes to use geometric moment to extract picture structure.Usually, the binary image that obtains behind thresholding (binarized image) is gone up and is calculated moment of inertia.The problem of doing like this is that threshold value is difficult to choose usually, and may not allow to detect [, because low threshold value will cause less vascular (tending to have than low-intensity) and adjacent structure to merge.Another problem of thresholding is that structure becomes " pixelation (pixelized) ", promptly manifests sharp edges, and these sharp edges make the true shape out of true of the calculating of its shape with respect to understructure.

The alternative of thresholding is to use the image intensity function f to calculate square as density function.But, noise (SN) than low zone in, be difficult to set up the structure that detects elongation about the threshold value of the eccentricity of fitted ellipse.For example, Fig. 1 (a) has described from the MIP of the sub-block (sub-volume) that real image extracted that centers on the vascular connecting portion.Top line shows the initial voxel value of using the arest neighbors interpolation.Middle row shows the bianry image that obtains behind manual thresholding.Threshold value is conditioned catches two vasculars, and this is to be difficult to the task of realization automatically.The pixelation effect of thresholding clearly, this has influenced the accuracy of shape description symbols.The third line shows the same sub-block body of the more complicated interpolation scheme of use.

Summary of the invention

One exemplary embodiment of the present invention described herein generally includes the system and method that is used for detecting automatically bright tubular structure, and based on the geometric moment that is used for from image extraction tubular structure the vascular of breast MR sequences is cut apart automatically.According to the method for the embodiment of the invention eigenwert, and be in harmonious proportion and carry out thresholding to image, wherein under the situation of low-down noise (SN) ratio, can restore structure reliably based on the shape tensor.As method based on the eigenwert of average structure tensor, according to the method for the embodiment of the invention and do not rely on single order image derivative, perhaps as method based on Hai Sai (Hessian) eigenwert, method according to the embodiment of the invention does not rely on second order image derivative yet, and avoided based on the method for Hai Sai or structure tensor intrinsic level and smooth to output.Method according to the embodiment of the invention can be moved fast, and each sequence only needs several seconds.The breast MR sequences of crossing based on the motion correction of test result shows, method according to the embodiment of the invention is cut apart vascular reliably, keep damage intact simultaneously, and all be better than differential technique aspect susceptibility and the bearing accuracy, and select parameter more insensitive yardstick.

According to aspect of the present invention, a kind of method that is used to cut apart digitized image is provided, this method comprises: digitized image is provided, and this digitized image comprises a plurality of intensity corresponding to the some territory on the 3D grid; Define the form matrix of described institute reconnaissance according to the square of the intensity in the some window of the institute's reconnaissance in described image; Calculate the eigenwert of described form matrix; Determine the eccentricity of the structure under described point according to described eigenwert; With cut apart described image based on described eccentricity value, wherein in the described image repeat the step of the eccentricity of the step of described definition form matrix, the step of eigenwert of calculating described form matrix and definite understructure a little.

According to another aspect of the present invention, the intermediate value that institute's reconnaissance has greater than predetermined threshold strengthens, and wherein the main substance (subjectmatter) to described digitized image applies contrast-enhancing agent before gathering described image.

According to another aspect of the present invention, be normalized in the predetermined scope with the difference of the intermediate value that strengthens preceding enhancing image and with described difference by the intermediate value of getting described contrast image, calculate intermediate value and strengthen.

According to another aspect of the present invention, form matrix S _αBe defined as

$S_{\mathrm{\α}}=\left(\begin{array}{ccc}{\mathrm{\μ}}_{\mathrm{xx},\mathrm{\α}}& {\mathrm{\μ}}_{\mathrm{xy},\mathrm{\α}}& {\mathrm{\μ}}_{\mathrm{xz},\mathrm{\α}}\\ {\mathrm{\μ}}_{\mathrm{xy},\mathrm{\α}}& {\mathrm{\μ}}_{\mathrm{yy},\mathrm{\α}}& {\mathrm{\μ}}_{\mathrm{yz},\mathrm{\α}}\\ {\mathrm{\μ}}_{\mathrm{xz},\mathrm{\α}}& {\mathrm{\μ}}_{\mathrm{yz},\mathrm{\α}}& {\mathrm{\μ}}_{\mathrm{zz},\mathrm{\α}}\end{array}\right),$

Wherein

${\mathrm{\&mu;}}_{\mathrm{xx},\mathrm{\&alpha;}}=\frac{{\mathrm{<figure-callout\; id="2"\; label="m"\; filenames="A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{\mathrm{2,0,0},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{\mathrm{1,0,0},\mathrm{\&alpha;}}^2}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$

${\mathrm{\&mu;}}_{\mathrm{xy},\mathrm{\&alpha;}}=\frac{{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">\; <figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>\; </figure-callout>}}_{\mathrm{1,1,0},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{\mathrm{1,0,0},\mathrm{\&alpha;}}{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{\mathrm{0,1,0},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$ ${\mathrm{\&mu;}}_{\mathrm{yy},\mathrm{\&alpha;}}=\frac{{\mathrm{<figure-callout\; id="2"\; label="m"\; filenames="A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{0,\mathrm{2,0},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{\mathrm{0,1,0},\mathrm{\&alpha;}}^2}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$

${\mathrm{\&mu;}}_{\mathrm{xz},\mathrm{\&alpha;}}=\frac{{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">\; <figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>\; </figure-callout>}}_{\mathrm{1,0},1,\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{\mathrm{1,0,0},\mathrm{\&alpha;}}{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{\mathrm{0,0,1},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$ ${\mathrm{\&mu;}}_{\mathrm{yz},\mathrm{\&alpha;}}=\frac{{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">\; <figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>\; </figure-callout>}}_{\mathrm{0,1,1},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{\mathrm{0,1,0},\mathrm{\&alpha;}}{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{\mathrm{0,0,1},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$ ${\mathrm{\&mu;}}_{\mathrm{zz},\mathrm{\&alpha;}}=\frac{{\mathrm{<figure-callout\; id="2"\; label="m"\; filenames="A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{0,0,2,\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{\mathrm{0,0,1},\mathrm{\&alpha;}}^2}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$

Square m wherein _{P, q, r, α}Be defined as

$m_{p,q,r,\mathrm{\&alpha;}}(x_0,y_0,z_0)={\&Integral;}_{R^3}{(x-x_0)}^p{(y-y_0)}^q{(z-z_0)}^{r}f{(x,y,z)}^{\mathrm{\&alpha;}}w(x-x_0,y-y_0,z-z_0)\mathrm{dxdydz},$

Wherein w is the window function with tight support (contact support) p, q, r μ 0 and α μ 1.

According to another aspect of the present invention, by on the limited neighborhood of each point and calculate integration.

According to another aspect of the present invention, window function is defined by following formula:

V wherein _x, v _y, v _zBe the picture point spacing, N _x, N _y, N _zBe defined nonnegative integer, wherein window size comprises maximum diameter interested.

According to another aspect of the present invention, this method comprises uses the described square of arest neighbors interpolation calculation, and proofreaies and correct described form matrix according to following formula

${\hat{S}}_{\mathrm{\&alpha;}}+\frac{1}{12}\left(\begin{array}{ccc}{v}_x^2& 0& 0\\ 0& v_{\mathrm{<figure-callout\; id="2"\; label="y"\; filenames="A20068003676700201.png"\; state="\{\{state\}\}">y</figure-callout>}}^2& 0\\ 0& 0& v_{\mathrm{<figure-callout\; id="2"\; label="z"\; filenames="A20068003676700201.png"\; state="\{\{state\}\}">z</figure-callout>}}^2\end{array}\right),$

V wherein _x, v _y, v _zIt is the picture point spacing.

According to another aspect of the present invention, this method comprises uses Tri linear interpolation to calculate described square.

According to another aspect of the present invention, α=1, and proofread and correct described form matrix according to following formula

${\hat{S}}_{\mathrm{\&alpha;}}+\frac{1}{6}\left(\begin{array}{ccc}{v}_x^2& 0& 0\\ 0& v_{\mathrm{<figure-callout\; id="2"\; label="y"\; filenames="A20068003676700201.png"\; state="\{\{state\}\}">y</figure-callout>}}^2& 0\\ 0& 0& v_{\mathrm{<figure-callout\; id="2"\; label="z"\; filenames="A20068003676700201.png"\; state="\{\{state\}\}">z</figure-callout>}}^2\end{array}\right),$

V wherein _x, v _y, v _zIt is the picture point spacing.

According to another aspect of the present invention, computer-readable program storage device is provided, this program storage device has visibly comprised the programmed instruction that can be carried out the method step that is used to cut apart digitized image by computer run.

Description of drawings

Fig. 1 (a) show according to the embodiment of the invention from MIP around the sub-block that true picture extracted of vascular connecting portion.

Fig. 1 (b) shows according to the simulation vascular of the embodiment of the invention and in the detection that does not have to use under the situation of thresholding moment of inertia.

Fig. 2 shows the basis function that is used for the 1D linear interpolation according to the embodiment of the invention.

Fig. 3 (a)-(c) shows cutting apart according to the macrolesion of the embodiment of the invention.

Fig. 4 (a)-(c) shows cutting apart according to a plurality of little damages of the embodiment of the invention.

Fig. 5 shows use shape tensor the cutting apart the blood vessel structure among the breast MRI according to the embodiment of the invention.

Fig. 6 shows the process flow diagram based on the dividing method of square according to the embodiment of the invention.

Fig. 7 is the block diagram based on the illustrative computer system of the dividing method of square that is used to implement according to the embodiment of the invention.

Embodiment

One exemplary embodiment of the present invention described here generally includes the system and method that is used for detecting automatically bright tubular structure, and the vascular in the breast MR sequences is cut apart automatically.According to the method for the embodiment of the invention eigenwert based on the shape tensor.This method can with compare based on the method for the eigenwert of average Hai Sai with based on the method for the eigenwert of average structure tensor.Match the structured descriptor that can be regarded as second order according to the sea that second derivative defines.Similarly, described structure tensor is the structured descriptor of single order.Described shape tensor can be regarded as the structured descriptor of zeroth order.

As used herein such, term " image " refers to the multidimensional data of being made up of discrete pictorial element (for example, the voxel of the pixel of 2D image and 3D rendering).Image can be for example by computed tomography, Magnetic resonance imaging, ultrasonic imaging or well known to a person skilled in the art the medical image of the main body that any other medical image system is collected.This image also can be provided by non-medical environment, such as being provided by remote sense system, electron microscopy etc.Although can think that image is from R ³To the function of R, but method of the present invention is not limited to such image, and can be applied to the image of any dimension, for example 2D picture or 3D volume.For 2 dimension or 3 d images, the territory of image is 2 dimensions or 3 rectangular arrays of tieing up normally, and wherein each pixel or voxel can be with reference to one group 2 or 3 mutually orthogonal incompatible addressing of axle.Term used herein " numeral " or " digitized " will be represented by the digital collection system or by the numeral obtained from the conversion of analog image or the image or the volume (as suitable) of digitized format.

Method according to the embodiment of the invention is exerted one's influence to image intensity by the second order geometric moment of (bright) structure of calculating lower floor.Method can be used to still also can not have application process under the situation of this threshold value by to strengthening the binary image that the image applications threshold value is obtained after the initial enhancing.The eigenwert of second order geometric moment is the classical tool of shape characterization in the object identification.But these eigenwerts never are used as the wave filter that is used to extract picture structure.Given bianry image consider the little sub-block (its size is relevant with structures of interest) around each pixel, and the shape tensor is defined as the second moment of the position of bright voxel with respect to the center of this sub-block in that place.For center pixel bright and with the enough approaching voxel in center of the shape of lower floor, calculate the eigenwert of shape tensor, and will be worth λ ₁-λ ₂/ (λ ₁+ λ ₂) compose and give filter response, wherein λ ₂＞λ ₁It is eigenvalue of maximum.

According to the embodiment of the invention, how much 3D squares can be defined as:

$m_{p,q,r,\mathrm{\&alpha;}}(x_0,y_0,z_0)={\&Integral;}_{R^3}{(x-x_0)}^p{(y-y_0)}^q{(z-z_0)}^{r}f{(x,y,z)}^{\mathrm{\&alpha;}}w(x-x_0,y-y_0,z-z_0)\mathrm{dxdydz},$

Wherein w is positive and symmetrical window function, and this window function has tight support p, q, r μ 0 and the α μ 1 that the location is provided.The shape tensor on α rank is defined as according to these squares

$S_{\mathrm{\&alpha;}}=\left(\begin{array}{ccc}{\mathrm{\&mu;}}_{\mathrm{xx},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{xy},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{xz},\mathrm{\&alpha;}}\\ {\mathrm{\&mu;}}_{\mathrm{xy},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{yy},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{yz},\mathrm{\&alpha;}}\\ {\mathrm{\&mu;}}_{\mathrm{xz},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{yz},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{zz},\mathrm{\&alpha;}}\end{array}\right),$

Wherein

${\mathrm{\&mu;}}_{\mathrm{xx},\mathrm{\&alpha;}}=\frac{{\mathrm{<figure-callout\; id="2"\; label="m"\; filenames="A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{\mathrm{2,0,0},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{\mathrm{1,0,0},\mathrm{\&alpha;}}^2}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2}$

${\mathrm{\&mu;}}_{\mathrm{xy},\mathrm{\&alpha;}}=\frac{{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">\; <figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>\; </figure-callout>}}_{\mathrm{1,1,0},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{\mathrm{1,0,0},\mathrm{\&alpha;}}{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{\mathrm{0,1,0},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2}$ ${\mathrm{\&mu;}}_{\mathrm{yy},\mathrm{\&alpha;}}=\frac{{\mathrm{<figure-callout\; id="2"\; label="m"\; filenames="A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{0,\mathrm{2,0},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{\mathrm{0,1,0},\mathrm{\&alpha;}}^2}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2}$

${\mathrm{\&mu;}}_{\mathrm{xz},\mathrm{\&alpha;}}=\frac{{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">\; <figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>\; </figure-callout>}}_{\mathrm{1,0},1,\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{\mathrm{1,0,0},\mathrm{\&alpha;}}{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{\mathrm{0,0,1},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2}$ ${\mathrm{\&mu;}}_{\mathrm{yz},\mathrm{\&alpha;}}=\frac{{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">\; <figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>\; </figure-callout>}}_{\mathrm{0,1,1},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{\mathrm{0,1,0},\mathrm{\&alpha;}}{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{\mathrm{0,0,1},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2}$ ${\mathrm{\&mu;}}_{\mathrm{zz},\mathrm{\&alpha;}}=\frac{{\mathrm{<figure-callout\; id="2"\; label="m"\; filenames="A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{0,0,2,\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{{\mathrm{<figure-callout\; id="1"\; label="m"\; filenames="A20068003676700191.png,A20068003676700201.png"\; state="\{\{state\}\}">m</figure-callout>}}_{\mathrm{0,0,1},\mathrm{\&alpha;}}^2}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2}.$

This matrix is symmetrical, so its all eigenwerts are real numbers.Making three eigenwerts is λ ₃＞λ ₂＞λ ₁μ 0, and then filter response can be defined as

For wire or the cylindrical-shaped structure such as vascular, C _Lineλ 1.

According to embodiments of the invention, based on S _αThe eigenwert 0[λ of (α＞＞1) ₁[λ ₂[λ ₃Calculate the eccentricity of the shape of lower floor.Along with α becomes big, higher intensity level is given more importance, acts on like that thereby almost be similar to thresholding.As shown in the simulated experiment of Fig. 1 (b), the high value of α can be dealt with low-down SN ratio, and the synthetic tubular structure that wherein has the even noise of increase utilizes the shape tensor of classical (classic) inertial matrix and α=15 detected.

Fig. 1 (b) shows simulation vascular and detection thereof, this detections utilization do not have thresholding standard moment of inertia and utilize the shape tensor of α=15.Row from left to right illustrate: (1) is the center section of synthetic volume initially, (2) its maximum intensity projection (MIP), (3) removed the MIP of the volume of vascular by standard square method, (4) MIP by the detected vascular of square method, (5) use the shape tensor of α=15 to remove the MIP of volume of vascular and the MIP that (6) use the detected vascular of shape tensor of α=15.The increase level of the even noise of six line display additivitys, thus following SN given ratio from top to bottom respectively: (1) 56.3, (2) 36.7, (3) 20.4, (4) 11.6, (5) 5.5 and (6) 0.8dB.For every kind of algorithm, be identical between the threshold value on the excentricity of shape is expert at.In all cases, for S ₁₅, examination criteria is $\frac{{\mathrm{\&lambda;}}_3}{{\mathrm{\&lambda;}}_2}>15,$ And for corresponding to S ₁Inertial matrix, examination criteria is $\frac{{\mathrm{\&lambda;}}_3}{{\mathrm{\&lambda;}}_2}>2.$ Under actual conditions, have been noted that this improved detection performance.

In fact, above-mentioned integration usually by on the limited neighborhood around each voxel and replace because f is only known at the voxel location place.Can be at all experiment supposition, mapping function is provided by following formula:

V wherein _x, v _y, v _zBe the image voxel spacing, and N _x, N _y, N _zBe defined nonnegative integer, make window size comprise maximum diameter interested.Then, given image is considered little sub-block and definition around each pixel

${\hat{m}}_{p,q,r,\mathrm{\&alpha;}}=\underset{i=1}{\overset{{2N}_x}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{{2N}_y}{\mathrm{\&Sigma;}}}\underset{k=0}{\overset{{2N}_z}{\mathrm{\&Sigma;}}}{\left({\mathrm{iv}}_x\right)}^p{\left({\mathrm{jv}}_y\right)}^q{\left({\mathrm{kv}}_z\right)}^r{\mathrm{\&rho;}}_{\mathrm{ijk}}^{\mathrm{\&alpha;}},$

ρ wherein _IjkIt is image value corresponding to the voxel place of index i, j, k.The eigenwert 0[λ of compute matrix ₁[λ ₂[λ ₃

${\hat{S}}_{\mathrm{\&alpha;}}=\left(\begin{array}{ccc}{\widehat{\mathrm{\&mu;}}}_{\mathrm{xx},\mathrm{\&alpha;}}& {\widehat{\mathrm{\&mu;}}}_{\mathrm{xy},\mathrm{\&alpha;}}& {\widehat{\mathrm{\&mu;}}}_{\mathrm{xz},\mathrm{\&alpha;}}\\ {\widehat{\mathrm{\&mu;}}}_{\mathrm{xy},\mathrm{\&alpha;}}& {\widehat{\mathrm{\&mu;}}}_{\mathrm{yy},\mathrm{\&alpha;}}& {\widehat{\mathrm{\&mu;}}}_{\mathrm{yz},\mathrm{\&alpha;}}\\ {\widehat{\mathrm{\&mu;}}}_{\mathrm{xz},\mathrm{\&alpha;}}& {\widehat{\mathrm{\&mu;}}}_{\mathrm{yz},\mathrm{\&alpha;}}& {\widehat{\mathrm{\&mu;}}}_{\mathrm{zz},\mathrm{\&alpha;}}\end{array}\right),$

As above calculated value wherein

But be to use the summation square to calculate.The excentricity of understructure or length growth rate can be measured ε=(λ by the excentricity of classics ₃-λ ₂)/(λ ₃+ λ ₂) measure, this measures the value of getting between 0 to 1, perhaps gets ratio λ simply ₃/ λ ₂(supposition λ ₂＞0)

Owing to do not suppose the differentiability of image intensity function f based on the method for square, thus can use simple interpolation scheme, such as arest neighbors interpolation or Tri linear interpolation, with the integration of the function that calculates interpolation, rather than on the voxel value and.Can expect, particularly under the situation of Tri linear interpolation, use these principal value of integrals to have accuracy preferably.Use equation

${\&Integral;}_{(i-1/2)v_x}^{(i+1/2)v_x}\mathrm{dx}=v_x,$

${\&Integral;}_{(i-1/2)v_x}^{(i+1/2)v_x}\mathrm{xdx}=v_x^{2}i,$

${\&Integral;}_{(i-1/2)v_x}^{(i+1/2)v_x}{x}^2\mathrm{dx}=v_x^3({\mathrm{<figure-callout\; id="2"\; label="i"\; filenames="A20068003676700201.png"\; state="\{\{state\}\}">i</figure-callout>}}^2+\frac{1}{12}),$

Can see, for arest neighbors interpolation integration, top matrix

Should replace by following formula

${\hat{S}}_{\mathrm{\&alpha;}}+\frac{1}{12}\left(\begin{array}{ccc}{\mathrm{<figure-callout\; id="2"\; label="v\; x"\; filenames="A20068003676700201.png"\; state="\{\{state\}\}">v</figure-callout>}}_x^{}& 0& 0\\ 0& {\mathrm{<figure-callout\; id="2"\; label="v\; y"\; filenames="A20068003676700201.png"\; state="\{\{state\}\}">v</figure-callout>}}_y^{}& 0\\ 0& 0& {\mathrm{<figure-callout\; id="2"\; label="v\; z"\; filenames="A20068003676700201.png"\; state="\{\{state\}\}">v</figure-callout>}}_z^{}\end{array}\right).$

Under the situation of Tri linear interpolation, function f is by ∑ _{I, j, k}ρ _Ijkg _IjkProvide, wherein i, j, k are the index of image voxel, ρ _IjkBe the image value of voxel, and

Then, write out $\underset{\mathrm{xyz}}{\&Integral;}(\&CenterDot;)\&equiv;{\&Integral;}_{(k-1)v_z}^{(k+1)v_z}{\&Integral;}_{(j-1)v_y}^{(j+1)v_y}{\&Integral;}_{(i-1)v_x}^{(i+1)v_x}(\&CenterDot;)\mathrm{dxdydz}:$

$\underset{\mathrm{xyz}}{\&Integral;}{g}_{\mathrm{ijk}}=v_x{v}_y{v}_z,$

${\&Integral;}_{\mathrm{xyz}}{\mathrm{xg}}_{\mathrm{ijk}}={\mathrm{iv}}_x^2{v}_y{v}_z,$

${\&Integral;}_{\mathrm{xyz}}{\mathrm{xyg}}_{\mathrm{ijk}}={\mathrm{ijv}}_x^2{v}_y^2{v}_z,$

${\&Integral;}_{\mathrm{xyz}}{x}^2{g}_{\mathrm{ijk}}=({\mathrm{<figure-callout\; id="2"\; label="i"\; filenames="A20068003676700201.png"\; state="\{\{state\}\}">i</figure-callout>}}^2+\frac{1}{6})v_x^3{v}_y{v}_z,$

Make, for the Tri linear interpolation under the situation of α=1, matrix

Should replace by following formula

${\hat{S}}_{\mathrm{\&alpha;}}+\frac{1}{6}\left(\begin{array}{ccc}{v}_x^2& 0& 0\\ 0& {\mathrm{<figure-callout\; id="2"\; label="v\; y"\; filenames="A20068003676700201.png"\; state="\{\{state\}\}">v</figure-callout>}}_y^{}& 0\\ 0& 0& {\mathrm{<figure-callout\; id="2"\; label="v\; z"\; filenames="A20068003676700201.png"\; state="\{\{state\}\}">v</figure-callout>}}_z^{}\end{array}\right).$

Under the situation of the general shape tensor (α＞1) that uses Tri linear interpolation, this situation becomes more complicated, wherein f ^αBy (∑ _{I, j, k}ρ _Ijkg _Ijk) ^αProvide.Although corresponding integration still can enclosed calculate, complexity significantly increases.According to embodiments of the invention,, should note as under the situation of above-mentioned α=1, calculating g in order to calculate corresponding shape tensor _Ijk ^αSquare die on.Continue, above-mentioned square can use not too directly but the method that can be generalized to α＞1 obtains.Under the 1D situation, can be thus completed, can directly be generalized to the situation of 2D and 3D.Suppose for k＜i-2 or k＞i+2, ρ _k=0, then obtain

${\&Integral;}_{(i-2)v_x}^{(i+2)v_x}f\left(x\right)\mathrm{dx}=\underset{k=i-1}{\overset{i+1}{\mathrm{\&Sigma;}}}{\&Integral;}_{(k-1)v_x}^{(k+1)v_x}{\mathrm{\&rho;}}_k{g}_k\left(x\right)\mathrm{dx}$

$={\&Integral;}_{(i-2)v_z}^{(i-1)v_z}{\mathrm{\&rho;}}_{i-1}{g}_{i-1}\left(x\right)\mathrm{dx}+{\&Integral;}_{(i-1)v_z}^{i{v}_x}({\mathrm{\&rho;}}_{i-1}{g}_i\left(x\right)+{\mathrm{\&rho;}}_i{g}_i\left(x\right))\mathrm{dx}$

$+{\&Integral;}_{i{v}_z}^{(i+1)v_z}({\mathrm{\&rho;}}_i{g}_i\left(x\right)+{\mathrm{\&rho;}}_{i+1}{g}_{i+1}\left(x\right))\mathrm{dx}+{\&Integral;}_{(i+1)v_x}^{(i+2)v_x}{\mathrm{\&rho;}}_{i+1}{g}_{i+1}\left(x\right)\mathrm{dx}$

$=v_x({\mathrm{\&rho;}}_{i-1}+{\mathrm{\&rho;}}_i+{\mathrm{\&rho;}}_{i+1}).$

Four top integrations can obtain according to the piecewise linear basis function of 3 shown in Fig. 2.With reference to this figure, the first basis function g _I-1Be defined within territory (i-2) v _xTo iv _xOn, the second basis function g _iBe defined within territory (i-1) v _xTo (i+1) v _xOn, and the 3rd function g _I+1Be defined within territory iv _xTo (i+2) v _xOn.

The method of calculating integration can be generalized to α＞1.For example, can calculate

${\&Integral;}_{(i-2)v_x}^{(i+2)v_x}f{\left(x\right)}^2\mathrm{dx}={\&Integral;}_{(i-2)v_x}^{(i-1)v_x}{\left({\mathrm{\&rho;}}_{i-1}{g}_{i-1}\left(x\right)\right)}^2\mathrm{dx}+{\&Integral;}_{(i-1)v_x}^{{\mathrm{iv}}_x}{({\mathrm{\&rho;}}_{i-1}{g}_i\left(x\right)+{\mathrm{\&rho;}}_i{g}_i\left(x\right))}^2\mathrm{dx}$

$+{\&Integral;}_{i{v}_x}^{(i+1)v_x}{({\mathrm{\&rho;}}_i{g}_i\left(x\right)+{\mathrm{\&rho;}}_{i+1}{g}_{i+1}\left(x\right))}^2\mathrm{dx}+{\&Integral;}_{(i+1)v_z}^{(i+2)v_z}{\left({\mathrm{\&rho;}}_{i+1}{g}_{i+1}\left(x\right)\right)}^2\mathrm{dx}$

$=(\frac{2}{3}{\mathrm{\&rho;}}_{i-1}^2+\frac{1}{3}{\mathrm{\&rho;}}_{i-1}{\mathrm{\&rho;}}_i+\frac{2}{3}{\mathrm{\&rho;}}_i^2+\frac{1}{3}{\mathrm{\&rho;}}_i{\mathrm{\&rho;}}_{i+1}+\frac{2}{3}{\mathrm{\&rho;}}_{i+1}^2)v_x$

Similarly,

${\&Integral;}_{(i-2)v_x}^{(i+2)v_x}\mathrm{xf}{\left(x\right)}^2\mathrm{dx}=({\mathrm{i\&rho;}}_i^2+(i-1){\mathrm{\&rho;}}_{i-1}^2+(\frac{1}{2}i-\frac{1}{4}){\mathrm{\&rho;}}_{i-1}{\mathrm{\&rho;}}_i+(\frac{1}{2}i+\frac{1}{4}){\mathrm{\&rho;}}_i{\mathrm{\&rho;}}_{i+1}+(i+1){\mathrm{\&rho;}}_{i+1}^2){\mathrm{<figure-callout\; id="2"\; label="v\; x"\; filenames="A20068003676700201.png"\; state="\{\{state\}\}">v</figure-callout>}}_x^{},$

And

${\&Integral;}_{(i-2)v_x}^{(i+2)v_x}{x}^{2}f{\left(x\right)}^2\mathrm{dx}=\left(\begin{array}{c}(\frac{11}{15}+\frac{4}{3}i+\frac{2}{3}{i}^2){\mathrm{\&rho;}}_{i+1}^2+(\frac{1}{10}+\frac{1}{3}i+\frac{1}{3}{i}^2){\mathrm{\&rho;}}_i{\mathrm{\&rho;}}_{i+1}\\ +(\frac{1}{15}+\frac{2}{3}{i}^2){\mathrm{\&rho;}}_i^2+(\frac{1}{10}-\frac{1}{3}i+\frac{1}{3}{i}^2){\mathrm{\&rho;}}_{i-1}{\mathrm{\&rho;}}_i+(\frac{11}{15}-\frac{4}{3}i+\frac{2}{3}{i}^2){\mathrm{\&rho;}}_{i-1}^2\end{array}\right){\mathrm{<figure-callout\; id="3"\; label="v\; x"\; filenames="A20068003676700211.png"\; state="\{\{state\}\}">v</figure-callout>}}_x^{}$

Although may find general formula at 3D situation and given α＞1, improve for potential degree of accuracy, final polynomial complexity is quite high.Under the 2D situation, four top integrations become 16 integrations, and become 64 integrations in 3D.

Method according to the embodiment of the invention is tested on the breast MR dynamic sequence of crossing above 100 motion corrections.The result who is obtained shows, can cut apart vascular reliably, keeps damage intact simultaneously.According to the embodiment of the invention, on the moving window of fixed measure, calculate square, but only consider that intermediate value strengthens the point that is higher than given threshold value.It is enough low that this threshold value can be chosen, so that detect even little vascular.This is not difficult to be provided with because and do not rely on calculating, and just quicken whole process by handling less voxel.Deduct the value that obtains thing before strengthening and calculate the intermediate value enhancing by will strengthen intermediate value that the back obtains thing in each image voxel.By using affine function, make final enhancing in scope [0,200] then with this difference standardization.Fig. 3 (a)-(c), Fig. 4 (a)-(c) and Fig. 5 show described result's several typical example.

Fig. 3 (a)-(c) shows cutting apart of macrolesion, and Fig. 4 (a)-(c) shows cutting apart of a plurality of little damages.For these figure, picture (panel) (a) shows after the initial enhancing of thresholding and strengthens image, and picture (b) shows detected vascular, and picture (c) shows the damage that removes vascular.

Fig. 5 shows shape tensor the cutting apart the blood vessel structure among the breast MRI of using α=6.These three row show the orthogonal view of same patient.First row illustrates the initial MIP that intermediate value strengthens.Second row illustrates the equal volume that removes vascular automatically.The third line illustrates the MIP of the vascular that removes separately.Notice that the diameter vascular very different with enhanced level correctly cut apart.Make λ by the eigenwert of getting the shape tensor ₃/ λ ₂Detection is carried out in＞3 place.

In each figure of these figure, notice that how little vascular is all correctly cut apart, and how little sphere structure all remains intact.As further checking, under 40 kinds of situations, extract blood vessel structure, 75 places damage altogether that three radiologist have observed these structures and mark according to method of the present invention.Vascular is correctly cut apart in all cases, and all damages that are labeled all remain intact.

Figure 6 illustrates process flow diagram based on the dividing method of square according to the embodiment of the invention.With reference now to this figure,,, provides and want divided image in step 61.At calculating the shape tensor in the image as at the determined voxel of step 62, the intermediate value contrast of this voxel strengthens above predetermined threshold value.In step 63, in the window of the fixed measure that centers on selected voxel, calculate square according to its definition shape tensor.In step 64, calculate the eigenwert of this shape tensor, and, determine the eccentricity of understructure in step 65.In step 66, this process circulation, processed up to each voxel.In step 67, based on the eccentricity split image of deriving by the shape tensor.

Extract local shape information based on the method for square can with compare based on the method for higher-order image derivative.For example, gradient square tensor (GST, Gradient SquareTensor) (or structure tensor) has been proposed as the robust method of estimating the partial structurtes dimension.Therefore this method is based on first order derivative and can be called the structured descriptor of single order.The eigenwert of Hai Sai also provides local image structure information, and set point etc. the principal curvatures of level (isolevel).Hai Sai and principal curvatures define according to second derivative, and thereby can be called the structured descriptor of second order.The shape tensor can be counted as the structured descriptor of zeroth order.The shape tensor is based on integration, and thereby with compare based on the method for single order or second derivative, the shape tensor has very sane attribute for noise.In addition, do not need to suppose any differentiability on the image function, this has simplified modeling.Problem based on the method for shape tensor is that connecting portion is not detected.In addition, need to understand preferably and determine whether and to come the computational geometry shape attribute according to the eigenwert of the shape tensor of α＞1.

Should be understood that the present invention can make up with various forms of hardware, software, firmware, dedicated process or its realizes.In one embodiment, the present invention can be implemented as the application program that visibly is comprised on the computer-readable program memory device with software.This application program can be uploaded in the machine that comprises any appropriate system framework and by this machine run.

Fig. 7 is the block diagram based on the illustrative computer system of the dividing method of square that is used to implement according to the embodiment of the invention.With reference now to Fig. 7,, is used to implement computer system 71 of the present invention and wherein can comprises CPU (central processing unit) (CPU) 72, storer 73 and I/O (I/O) interface 74.Computer system 71 is coupled to display 75 and various input equipment 76 by I/O interface 74 usually, all mouses in this way of input equipment and keyboard.Support circuit can comprise circuit and communication bus such as high-speed cache, power supply, clock circuit.Storer 73 can comprise random-access memory (ram), ROM (read-only memory) (ROM), hard disk drive, tape drive etc. or its combination.The present invention may be implemented as the routine 77 that is stored in the storer 73 and handles signal from signal source 78 by CPU 72 operations.Similarly, computer system 71 is general-purpose computing systems, and when operation routine 77 of the present invention, this general-purpose computing system becomes dedicated computer system.

Computer system 71 also comprises operating system and micro-instruction code.Various process described herein and function can be the parts of micro-instruction code, or via the part (or its combination) of the application program of operating system.In addition, various other peripherals can be connected to computer platform, such as being connected to additional data storage device and printing device.

Should be further understood that, because some shown in the accompanying drawing are formed system units and method step can be realized with software, so the actual connection between these system units (or process steps) can depend on the mode of the present invention's programming and difference.The given instruction of the present invention that provides here, those skilled in the art can imagine these and similarly embodiment or configuration of the present invention.

Although describe the present invention in detail with reference to preferred embodiment, it will be understood by those skilled in the art that under the situation of the spirit and scope of the present invention of in not departing from the claims of enclosing, being set forth and to carry out various modifications and replacement to it.

Claims (20)

1. method of cutting apart digitized image, this method may further comprise the steps:

Digitized image is provided, and this digitized image comprises a plurality of intensity corresponding to the some territory on the 3D grid;

Define the form matrix of described institute reconnaissance according to the square of the intensity in the some window of the institute's reconnaissance in described image;

Calculate the eigenwert of described form matrix;

According to described eigenwert determine the structure under described point eccentricity and

Cut apart described image based on described eccentricity value, wherein in the described image re-define the step of the eccentricity of the step of form matrix, the step of eigenwert of calculating described form matrix and definite understructure a little.

2. method according to claim 1, wherein, the intermediate value that the reconnaissance of described institute has greater than predetermined threshold strengthens, and wherein the main substance to described digitized image applies contrast-enhancing agent before gathering described image.

3. method according to claim 2 wherein, is normalized in the preset range with the difference of the intermediate value that strengthens preceding enhancing image and with described difference by the intermediate value of getting described contrast image, calculates described intermediate value and strengthens.

4. method according to claim 1, wherein, described form matrix S _αBe defined as

$S_{\mathrm{\&alpha;}}=\left(\begin{array}{ccc}{\mathrm{\&mu;}}_{\mathrm{xx},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{xy},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{xz},\mathrm{\&alpha;}}\\ {\mathrm{\&mu;}}_{\mathrm{xy},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{yy},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{yz},\mathrm{\&alpha;}}\\ {\mathrm{\&mu;}}_{\mathrm{xz},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{yz},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{zz},\mathrm{\&alpha;}}\end{array}\right),$

Wherein,

${\mathrm{\&mu;}}_{\mathrm{xx},\mathrm{\&alpha;}}=\frac{m_{\mathrm{2,0,0},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{m_{\mathrm{1,0,0},\mathrm{\&alpha;}}^2}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$

${\mathrm{\&mu;}}_{\mathrm{xy},\mathrm{\&alpha;}}=\frac{m_{\mathrm{1,1,0},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{m_{\mathrm{1,0,0},\mathrm{\&alpha;}}{m}_{\mathrm{0,1,0},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$ ${\mathrm{\&mu;}}_{\mathrm{yy},\mathrm{\&alpha;}}=\frac{m_{0,\mathrm{2,0},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{m_{\mathrm{0,1,0},\mathrm{\&alpha;}}^2}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$

${\mathrm{\&mu;}}_{\mathrm{xz},\mathrm{\&alpha;}}=\frac{m_{\mathrm{1,0},1,\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{m_{\mathrm{1,0,0},\mathrm{\&alpha;}}{m}_{\mathrm{0,0,1},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$ ${\mathrm{\&mu;}}_{\mathrm{yz},\mathrm{\&alpha;}}=\frac{m_{\mathrm{0,1,1},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{m_{\mathrm{0,1,0},\mathrm{\&alpha;}}{m}_{\mathrm{0,0,1},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$ ${\mathrm{\&mu;}}_{\mathrm{zz},\mathrm{\&alpha;}}=\frac{m_{0,0,2,\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{m_{\mathrm{0,0,1},\mathrm{\&alpha;}}^2}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$

Square m wherein _{P, q, r, α}Be defined as

$m_{p,q,r,\mathrm{\&alpha;}}(x_0,y_0,z_0)={\&Integral;}_{R^3}{(x-x_0)}^p{(y-y_0)}^q{(z-z_0)}^{r}f{(x,y,z)}^{\mathrm{\&alpha;}}w(x-x_0,y-y_0,z-z_0)\mathrm{dxdydz},$

Wherein w is the window function with tight support p, q, r μ 0 and α μ 1.

5. method according to claim 4, wherein, by on the limited neighborhood of each point and calculate integration.

6. method according to claim 4, wherein, described window function is defined by following formula

Wherein, v _x, v _y, v _zBe the picture point spacing, N _x, N _y, N _zBe defined nonnegative integer, wherein window size comprises maximum diameter interested.

7. method according to claim 4 further comprises and uses the described square of arest neighbors interpolation calculation, and proofreaies and correct described form matrix according to following formula

${\hat{S}}_{\mathrm{\&alpha;}}+\frac{1}{12}\left(\begin{array}{ccc}{v}_x^2& 0& 0\\ 0& v_y^2& 0\\ 0& 0& v_z^2\end{array}\right),$

V wherein _x, v _y, v _zIt is the picture point spacing.

8. method according to claim 4 further comprises and uses Tri linear interpolation to calculate described square.

9. method according to claim 8, wherein, α=1, and proofread and correct described form matrix according to following formula

${\hat{S}}_{\mathrm{\&alpha;}}+\frac{1}{6}\left(\begin{array}{ccc}{v}_x^2& 0& 0\\ 0& v_y^2& 0\\ 0& 0& v_z^2\end{array}\right),$

V wherein _x, v _y, v _zIt is the picture point spacing.

10. method of cutting apart digitized image, this method may further comprise the steps:

Digitized image is provided, and this digitized image comprises a plurality of intensity corresponding to the some territory on the 3D grid;

Define the form matrix of described institute reconnaissance, wherein said form matrix S according to the square of the intensity in the some window of the institute's reconnaissance in the described image _αBe defined as

$S_{\mathrm{\&alpha;}}=\left(\begin{array}{ccc}{\mathrm{\&mu;}}_{\mathrm{xx},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{xy},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{xz},\mathrm{\&alpha;}}\\ {\mathrm{\&mu;}}_{\mathrm{xy},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{yy},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{yz},\mathrm{\&alpha;}}\\ {\mathrm{\&mu;}}_{\mathrm{xz},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{yz},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{zz},\mathrm{\&alpha;}}\end{array}\right),$

Wherein,

${\mathrm{\&mu;}}_{\mathrm{xx},\mathrm{\&alpha;}}=\frac{m_{\mathrm{2,0,0},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{m_{\mathrm{1,0,0},\mathrm{\&alpha;}}^2}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$

${\mathrm{\&mu;}}_{\mathrm{xy},\mathrm{\&alpha;}}=\frac{m_{\mathrm{1,1,0},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{m_{\mathrm{1,0,0},\mathrm{\&alpha;}}{m}_{\mathrm{0,1,0},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$ ${\mathrm{\&mu;}}_{\mathrm{yy},\mathrm{\&alpha;}}=\frac{m_{0,\mathrm{2,0},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{m_{\mathrm{0,1,0},\mathrm{\&alpha;}}^2}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$

${\mathrm{\&mu;}}_{\mathrm{xz},\mathrm{\&alpha;}}=\frac{m_{\mathrm{1,0},1,\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{m_{\mathrm{1,0,0},\mathrm{\&alpha;}}{m}_{\mathrm{0,0,1},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$ ${\mathrm{\&mu;}}_{\mathrm{yz},\mathrm{\&alpha;}}=\frac{m_{\mathrm{0,1,1},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{m_{\mathrm{0,1,0},\mathrm{\&alpha;}}{m}_{\mathrm{0,0,1},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$ ${\mathrm{\&mu;}}_{\mathrm{zz},\mathrm{\&alpha;}}=\frac{m_{0,0,2,\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{m_{\mathrm{0,0,1},\mathrm{\&alpha;}}^2}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$

Square m wherein _{P, q, r, α}Be defined as

$m_{p,q,r,\mathrm{\&alpha;}}(x_0,y_0,z_0)={\&Integral;}_{R^3}{(x-x_0)}^p{(y-y_0)}^q{(z-z_0)}^{r}f{(x,y,z)}^{\mathrm{\&alpha;}}w(x-x_0,y-y_0,z-z_0)\mathrm{dxdydz},$

Wherein w is the window function with tight support p, q, r μ 0 and α μ 1;

Calculate the eigenwert of described form matrix; With

Determine the eccentricity of the structure under described point according to described eigenwert.

11. method according to claim 10, further comprise in the described image re-define the step of the eccentricity of the step of form matrix, the step of eigenwert of calculating described form matrix and definite understructure a little, and cut apart described image based on described eccentricity value.

12. a computer-readable program storage device, it visibly comprises the programmed instruction of being carried out the method step that is used to cut apart digitized image by computer run, and the method comprising the steps of:

Digitized image is provided, and this digitized image comprises a plurality of intensity corresponding to the some territory on the 3D grid;

Define the form matrix of described institute reconnaissance according to the square of the intensity in the some window of the institute's reconnaissance in described image;

Calculate the eigenwert of described form matrix;

According to described eigenwert determine the structure under described point eccentricity and

Cut apart described image based on described eccentricity value, wherein in the described image re-define the step of the eccentricity of the step of form matrix, the step of eigenwert of calculating described form matrix and definite understructure a little.

13. computer-readable program storage device according to claim 12, wherein, the intermediate value that the reconnaissance of described institute has greater than predetermined threshold strengthens, and wherein the main substance to described digitized image applies contrast-enhancing agent before gathering described image.

14. computer-readable program storage device according to claim 13 wherein, is normalized in the preset range with the difference of the intermediate value that strengthens preceding enhancing image and with described difference by the intermediate value of getting described contrast image, calculates described intermediate value and strengthens.

15. computer-readable program storage device according to claim 12, wherein, described form matrix S _αBe defined as

$S_{\mathrm{\&alpha;}}=\left(\begin{array}{ccc}{\mathrm{\&mu;}}_{\mathrm{xx},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{xy},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{xz},\mathrm{\&alpha;}}\\ {\mathrm{\&mu;}}_{\mathrm{xy},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{yy},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{yz},\mathrm{\&alpha;}}\\ {\mathrm{\&mu;}}_{\mathrm{xz},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{yz},\mathrm{\&alpha;}}& {\mathrm{\&mu;}}_{\mathrm{zz},\mathrm{\&alpha;}}\end{array}\right),$

Wherein,

${\mathrm{\&mu;}}_{\mathrm{xx},\mathrm{\&alpha;}}=\frac{m_{\mathrm{2,0,0},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{m_{\mathrm{1,0,0},\mathrm{\&alpha;}}^2}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$

${\mathrm{\&mu;}}_{\mathrm{xy},\mathrm{\&alpha;}}=\frac{m_{\mathrm{1,1,0},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{m_{\mathrm{1,0,0},\mathrm{\&alpha;}}{m}_{\mathrm{0,1,0},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$ ${\mathrm{\&mu;}}_{\mathrm{yy},\mathrm{\&alpha;}}=\frac{m_{0,\mathrm{2,0},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{m_{\mathrm{0,1,0},\mathrm{\&alpha;}}^2}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$

${\mathrm{\&mu;}}_{\mathrm{xz},\mathrm{\&alpha;}}=\frac{m_{\mathrm{1,0},1,\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{m_{\mathrm{1,0,0},\mathrm{\&alpha;}}{m}_{\mathrm{0,0,1},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$ ${\mathrm{\&mu;}}_{\mathrm{yz},\mathrm{\&alpha;}}=\frac{m_{\mathrm{0,1,1},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{m_{\mathrm{0,1,0},\mathrm{\&alpha;}}{m}_{\mathrm{0,0,1},\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$ ${\mathrm{\&mu;}}_{\mathrm{zz},\mathrm{\&alpha;}}=\frac{m_{0,0,2,\mathrm{\&alpha;}}}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}}-\frac{m_{\mathrm{0,0,1},\mathrm{\&alpha;}}^2}{m_{\mathrm{0,0,0},\mathrm{\&alpha;}}^2},$

Square m wherein _{P, q, r, α}Be defined as

$m_{p,q,r,\mathrm{\&alpha;}}(x_0,y_0,z_0)={\&Integral;}_{R^3}{(x-x_0)}^p{(y-y_0)}^q{(z-z_0)}^{r}f{(x,y,z)}^{\mathrm{\&alpha;}}w(x-x_0,y-y_0,z-z_0)\mathrm{dxdydz},$

Wherein w is the window function with tight support p, q, r μ 0 and α μ 1.

16. computer-readable program storage device according to claim 15, wherein, by on the limited neighborhood of each point and calculate described integration.

17. computer-readable program storage device according to claim 15, wherein, described window function is defined by following formula

V wherein _x, v _y, v _zBe the picture point spacing, N _x, N _y, N _zBe defined nonnegative integer, wherein window size comprises maximum diameter interested.

Use the described square of arest neighbors interpolation calculation and proofread and correct described form matrix 18. computer-readable program storage device according to claim 15, described method further comprise according to following formula

${\hat{S}}_{\mathrm{\&alpha;}}+\frac{1}{12}\left(\begin{array}{ccc}{v}_x^2& 0& 0\\ 0& v_y^2& 0\\ 0& 0& v_z^2\end{array}\right),$

V wherein _x, v _y, v _zIt is the picture point spacing.

19. further comprising, computer-readable program storage device according to claim 15, described method use Tri linear interpolation to calculate described square.

20. computer-readable program storage device according to claim 19, wherein, α=1, and proofread and correct described form matrix according to following formula

${\hat{S}}_{\mathrm{\&alpha;}}+\frac{1}{6}\left(\begin{array}{ccc}{v}_x^2& 0& 0\\ 0& v_y^2& 0\\ 0& 0& v_z^2\end{array}\right),$

V wherein _x, v _y, v _zIt is the picture point spacing.

Images