CN101071506A - System and method for local pulmonary structure classification for computer-aided nodule detection - Google Patents

Abstract

A method for classifying pulmonary structures in digitized images includes providing approximate target structure locations of one or more target structures in a digitized 3-dimensional (3D) image, fitting an anisotropic Gaussian model about said approximate target locations to generate more precise 3D target models and center locations of said one or more target structures, warping each said 3D target model into a 3D sphere, constructing a bounding manifold about each said warped 3D sphere, and identifying clusters on said bounding manifold wherein said one or more target structures are classified.

Description

Be used for the local pulmonary textural classification to carry out the system and method for area of computer aided nodule detection

The cross reference of relevant U. S. application

The U.S. Provisional Application No.60/761 that the application requires people such as Bahlmann to submit on January 25th, 2006, the right of priority of 927 " Local Pulmonary Structure Classification for Computer-Aided NoduleDetection (the local pulmonary textural classification is in order to the area of computer aided nodule detection) ", its content is introduced in here as a reference.

Technical field

The present invention is directed to the classification of partial structurtes type in digitized medical images.

Background technology

Only just had 160,000 people of surpassing to die from lung cancer in the U.S. last year.Though prevention lung cancer the best way is a non-smoking, early detection is to improve the key of patient's prognosis.When detecting cancer in early days and undergoing surgery, non-small cell lung cancer I phase patient's five-year survival rate is 60%-80%.Yet the five-year survival rate that the patient of motionless operation faces only is 10%.1.

Imaging technique (for example computer tomography (CT) scanning) provides non-invasive, sensitive early detection method.In Thoracic CT scan, the computer aided detection of lung tubercle and diagnosis (CAD) have been reduced the possibility of personal error for more effective and standardization diagnostic procedure.In CT scan, the lung tubercle presents the compact block of different shape and size.They may for example organize blood vessel with other or pleura separates or adhesion.

Recently, many technology have been proposed in thin layer (thin-slice) CT, tubercle detected automatically and to classify, comprise that region growing and automatic threshold are definite, template matches with Gauss (Gaussian) lung nodule model, use the 3D tubercle to select and the squelch filtrator, tubercle coupling, and variable-geometry and intensity template.Yet the shortcoming of these present situations of CAD system is to distinguish tubercle and other compact texture (for example blood vessel).Because employed circular hypothesis in most systems, crooked vascular and point of crossing thereof often are mistakenly detected as tubercle, thereby cause (FP) case of false positive (falsepositive).

In order to reduce the quantity of this class FP, proposed two kinds of solutions in the past: the filtrator based on correlativity is followed the tracks of by the vascular axis that analysis was provided based on Hai Sen (Hessian) to carry out vascular tree reconstruction and utilization in order to the area-of-interest that strengthens fuzzy shape analysis.The shortcoming of last method comprises its ineffective activity.Employed under study for action simple structure template can not be handled the vascular shape and the topology of many complexity.On the other hand, very big although back one method can be handled more irregular structure calculated amount.

Summary of the invention

Exemplary embodiments of the present invention as described herein generally comprises and is used for the method and system of local structure type (for example tubercle, vascular and point of crossing) being classified in Thoracic CT scan.This classification is useful for the computer aided detection (CAD) of lung tubercle, and can be used as the aftertreatment ingredient of any lung CAD system so that reduce the false positive (FP) that is produced by vascular and point of crossing.This classification therefore supposition provides positive candidate by this CAD system or from radiologist's report, thereby concentrates on the problem that reduces FP.In this scheme, thereby classification results provides a kind of eliminating to improve overall performance by the false-positive effective means that vascular and point of crossing cause.

Convert the classification of various 3D topological structures to simple more 2D data clusters classification according to the method for the embodiment of the invention, to its available in the literature more general and solution flexibly, and it more is applicable to visual.Except the benefit on calculating, this method has following advantage, promptly more general and analytical technology inventory (inventory) and more intuitive visual possibility flexibly, this with radiologist's possible reciprocation environment in be useful.

Given tubercle candidate makes anisotropic Gaussian be fit to data at first steadily.Resulting Gauss center and spreading parameter are used to affine standardized data territory, so that the anisotropy ellipsoid of match is bent to the isotropy spheroid of fixed size.A kind of automated process can extract the 3D spherical flow shape of the suitable border surface that contains target structure.Carrying out yardstick by data-driven entropy minimization method selects.Adopt EM cluster, the directional statistics that for example utilizes automatic modulus to estimate and utilize the technology such as hierarchical clustering of the Bhattacharyya distance of revising to come analysis stream shape at high strength group (cluster) corresponding to outstanding structure.The type of lung structure is determined on high strength group's estimate amount dominance ground: tubercle (0), adhesion tubercle (1), vascular (2), point of crossing (＞3).Utilize directional statistics, especially multivariate to hold Gauss's modeling (multivariate wrappedGaussian modeling), method has been expanded the Gauss curve fitting method according to an embodiment of the invention, comprises automatic modulus selection.

Outside the scope of the CAD of lung, sorting technique can be used to provide the significant information of blood vessel structure in various fields (for example angiography) according to an embodiment of the invention.These local programs are more flexible and effective than prior art, and will help to improve the accuracy of common lung's CAD system.In addition, selected volume of interest (VOI) expression has useful visual ability in the modeling part, for example helps the 2D border manifold (boundingmanifold) of the interactive expansion of user (radiologist).

In this field, must the issue a certificate qualitative examination to selected chest CT image examples of favourable classification.Algorithm can be classified to the example of tubercle, adhesion tubercle, vascular and vascular point of crossing steadily according to an embodiment of the invention.

According to an aspect of the present invention, a kind of method of lung's structure being classified at digitized image of being used for is provided, be included in the approximate target structure position that one or more target structures are provided in digitizing three-dimensional (3D) image, the center that the match anisotropic Gaussian model produces more accurate 3D target model and described one or more target structures around each described approximate target position, make each described 3D target model bend to the 3D spheroid, around each described bending 3D spheroid, make up border manifold, and the identification group wherein classifies to described one or more target structures on described border manifold.

According to a further aspect in the invention, digitized image comprise with 3 d grid on the corresponding a plurality of intensity in some territory.

According to a further aspect in the invention, the match anisotropic Gaussian model comprises employing Gauss yardstick-spatial mean value drift analysis and selects based on Jansen-Shannon (Jensen-Shannon) auto-bandwidth of divergence around approximate target position, produce the 3D ellipsoid model of described target structure, the ellipsoidal center of wherein said 3D is corresponding with the center and the covariance of described Gauss model with size.

According to a further aspect in the invention, crooked described 3D target model comprises the described 3D ellipsoid of affine standardization, and wherein zoom direction and coefficient obtain from the structure covariance of described anisotropic Gaussian model.

According to a further aspect in the invention, make up border manifold and comprise that further the 3D unfolded surface with bending spheroid becomes 2D to represent, and determine the radius of suitable border manifold.

According to a further aspect in the invention, with the 3D unfolded surface of bending spheroid become 2D represent to comprise surface transformation with described bending spheroid become spherical coordinates (θ, ), wherein  ∈ [π, π] and θ ∈ [π, π].

According to a further aspect in the invention, the radius of determining suitable border manifold comprises: the spherical flow shape that makes up a plurality of different radiis around described bending spheroid; Each spherical flow shape is launched into 2D to be represented; Will be in the intensity level distribution standardization on the spherical flow shape of each described expansion; Be the spherical flow shape calculating strength entropy of each described expansion, wherein intensity level is counted as probable value, and wherein entropy distributes and is defined; And find and make the described entropy minimized radius that distributes.

According to a further aspect in the invention, the identification group comprises that the use expectation maximization makes vector variable Θ=(θ ₁..., θ _F) ^TMultivariate hold Gaussian distribution N _W ^PMixing (Θ) $N_w\left(\mathrm{\&Theta;}\right)=\underset{p=1}{\overset{p}{\mathrm{\&Sigma;}}}{c}_p{N}_w^p\left(\mathrm{\&Theta;}\right)$ Be fit to by the outstanding object of described border manifold, described border manifold is followed the minimum description length criterion; Mixed components probability c wherein _pIn expectation maximization, estimated, wherein each the dimension in, θ _iSatisfy θ=χ _w=χ mod 2 π ∈ (π, π], N _W ^P(Θ) satisfy $N_w^p\left(\mathrm{\&Theta;}\right)=\underset{{\mathrm{<figure-callout\; id="1"\; label="k"\; filenames="A20071008604100231.png"\; state="\{\{state\}\}">k</figure-callout>}}_1=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\&CenterDot;\&CenterDot;\&CenterDot;\underset{k_F=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{N}^p(\mathrm{\&Theta;}+2\mathrm{\&pi;}{k}_1{\mathrm{<figure-callout\; id="1"\; label="e"\; filenames="A20071008604100231.png"\; state="\{\{state\}\}">e</figure-callout>}}_1+\&CenterDot;\&CenterDot;\&CenterDot;+{2\mathrm{\&pi;k}}_F{e}_F)$ , e wherein _k=(0 ..., 0,1,0 ..., 0) ^TBe to be k Euclid's basis vector, the k item be designated as 1 and other place be designated as 0, the wherein estimation of mixed components p

With

Be from sample set X={ θ in expectation maximization ⁽¹⁾..., 9 ^(M)In based on direction mean value ${\left({\widehat{\mathrm{\&mu;}}}_{\mathrm{\&theta;}}\right)}_f=\arg \left(\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\exp \left(i{\mathrm{\&theta;}}_f^{\left(m\right)}\right)\right)$ And covariance ${\widehat{\mathrm{\&Sigma;}}}_{\mathrm{\&theta;}}=\frac{1}{M-1}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&Theta;}}^{{\left(m\right)}^{\&prime;}}{\mathrm{\&Theta;}}^{{\left(m\right)}^{\&prime;T}}$ Obtain, wherein ${\mathrm{\&theta;}}^{{\left(m\right)}^{\&prime;}}=({\mathrm{\&theta;}}_f^{\left(m\right)}-{\left({\widehat{\mathrm{\&mu;}}}_{\mathrm{\&theta;}}\right)}_f)\mathrm{<figure-callout\; id="2"\; label="mod"\; filenames="A20071008604100221.png,A20071008604100231.png"\; state="\{\{state\}\}">mod</figure-callout>}2\mathrm{\&pi;}$ , and wherein observe (observation) X directly from 2D unfolded image I (θ, ) the middle extraction, the wherein (θ of each sampling _m,  _m) ∈ (π, π] * (π, π] occurrence number and corresponding image array value I (θ _m,  _m) be set pro rata.

According to a further aspect in the invention, described method comprises uses the cohesion hierarchical clustering that the group in the mutual preset distance is merged, and is used for multivariate and holds right the equaling of Gaussian distribution $\frac{1}{8}(({\mathrm{\&mu;}}_2-{\mathrm{\&mu;}}_1)\mathrm{<figure-callout\; id="2"\; label="mos"\; filenames="A20071008604100221.png,A20071008604100231.png"\; state="\{\{state\}\}">mos</figure-callout>}2\mathrm{\&pi;}{)}^T{\left(\frac{{\mathrm{\&Sigma;}}_1+{\mathrm{\&Sigma;}}_2}{2}\right)}^{-1}$ $\left(({\mathrm{\&mu;}}_2-{\mathrm{\&mu;}}_1)\mathrm{<figure-callout\; id="2"\; label="mod"\; filenames="A20071008604100221.png,A20071008604100231.png"\; state="\{\{state\}\}">mod</figure-callout>}2\mathrm{\&pi;}\right)+\frac{1}{2}\ln \frac{|{\mathrm{\&Sigma;}}_1+{\mathrm{\&Sigma;}}_2|}{\sqrt{\left|{\mathrm{\&Sigma;}}_1\right|\left|{\mathrm{\&Sigma;}}_2\right|}}$ Distance metric, μ wherein ₁And μ ₂Be the right mean value of Gaussian distribution, and ∑ ₁And ∑ ₂It is its variance separately.

According to a further aspect in the invention, the class of lung structure is determined that by the quantity that holds the gaussian component group relevant with target structure wherein Gu Li tubercle has 0 group, and the adhesion tubercle has 2 groups, and vascular has 4 groups, and the intersecting blood vessels point has 6 or multigroup more.

According to a further aspect in the invention, provide a kind of computer-readable program storage device, visibly embodied by the executable instruction repertorie of computing machine to be implemented in the method step of in the digitized image lung structure being classified.

Description of drawings

Fig. 1 (a)-(c) illustrates the method for lung structure classification according to an embodiment of the invention.

Fig. 2 (a)-(g) illustrates the effect of expansion ellipsoid and the corresponding image intensity entropy of histogram of different radii r according to an embodiment of the invention.

Fig. 3 illustrates the cluster (clustering) that has directional data according to an embodiment of the invention.

Fig. 4 (a)-(d) and Fig. 5 (a)-(d) illustrate the illustrative example of the lung structure sorting technique of the one embodiment of the invention that is used for Thoracic CT scan.

Fig. 6 (a)-(b) illustrates the example of directional data according to an embodiment of the invention.

Fig. 7 is the process flow diagram of sorting technique according to an embodiment of the invention.

Fig. 8 is the block diagram that is used to carry out the typical computer of sorting technique according to one embodiment of the invention.

Embodiment

The typical embodiment of the present invention as described herein generally comprises and is used for the system and method for local structure type being classified in chest scan.Therefore, although the present invention admits of various modifications and replaceable form, its specific embodiments is represented by giving an example in the accompanying drawings and here will be described in detail.Yet, should be understood that this is not intended that the invention be limited to particular forms disclosed, and opposite, the present invention has covered all modifications, equivalence and the replacement scheme that drops in the spirit and scope of the invention.

As used herein, term " image " refers to the multidimensional data of being made up of discrete pictorial element (for example, the voxel of the pixel of 2-D image and 3-D image).For example, image can be by computer tomography, magnetic resonance imaging, ultrasonic or the medical image to main body (subject) that any other medical image system known in those skilled in the art is gathered.Image also can for example be provided by the non-medical environment as these remote sensing systems, electronics microscopy etc.Although image can be regarded as from R ³To the function of R, yet method of the present invention is not restricted to these images, and can be applied to the image of any dimension, for example 2-D picture and 3-D volume.For 2 or 3 d image for, the territory of image typically 2 or the 3 d image rectangular array, wherein each pixel or voxel can be addressed with reference to one group 2 or 3 mutually orthogonal axles.Suitably, term used herein " numeral " or " digitized " refer at the numeral of obtaining system by numeral or obtaining from the conversion of analog image or the image or the volume of digitized format.

Categorizing system comprises according to an embodiment of the invention: (1) is used for the module of anisotropic Gaussian match, and (2) are at the structure of the 2D of lung structure boundary stream shape, and (3) are to the sane cluster analysis (cluster analysis) of described stream shape.(2) part is used based on the data-driven yardstick of entropy minimization and is selected.(3) part is used statistical analysis technique, for example with the cluster based on expectation maximization (EM) of automatic modulus selection (automatic mode numberselection), directional data modeling and based on the hierarchical clustering of Bhattacharyya apart from variable.High strength group's quantity will directly determine the lung structure class in this analysis.Different with other spherical method of rebuilding as vascular tree and so on, this method allows the lung structure inspection that localizes flexibly.

In the nodule detection application background, incorrect detection and vascular that cut apart and vessel branchings structure are represented false positive (FP) case.Sorting technique is refused all so non-nodular structures according to an embodiment of the invention, and as accessory substance, infers the kind of the lung structure type in the research, i.e. tubercle, adhesion tubercle, vascular or vascular point of crossing.And the solution of classifying according to an embodiment of the invention can play posttreatment filter in lung's CAD system, to reduce vascular and the caused FP in point of crossing.

The flow chart of lung's sorting technique is shown among Fig. 7 according to an embodiment of the invention.The apparent position of method supposition lung structure for example is provided from CAD system, radiologist's handbook reading matter or report etc. according to an embodiment of the invention.One presses (one-click) nodule segmentation algorithm can be used to the location and cut apart target structure, and described target structure comprises tubercle, adhesion tubercle, vascular or vascular point of crossing.Referring now to this figure,, both having fixed on step 71 provides the tubercle position candidate as this semi-automatic initialization of cutting apart solution.In step 72, make anisotropic Gaussian model be fit to target structure intensity, form the more accurate 3D ellipsoid model of target.In step 73, these ellipsoids are bent to the 3D spheroid.In step 74, make up border manifold from bending 3D spheroid.According to one embodiment of the invention, this structure comprises that the 3D unfolded surface with spheroid becomes the 2D spherical coordinates to represent, then determines the radius of suitable border manifold.In step 75, carry out the cluster analysis of border manifold, then carry out aftertreatment in step 76.The details of these steps will be described below.

Referring to step 72, according to an embodiment of the invention algorithm based on: use Gauss's yardstick-spatial mean value drift to analyze and select to make steadily anisotropic intensity model to be fit to data based on Gauss based on the Jansen-Shannon auto-bandwidth of divergence.These technology from coarse CAD or manually initialization accurate estimation to pinwheel is provided.From the target structure border, extract 3D ellipsoid stream shape.Yet the ellipsoid match is not unessential usually, alleviates this task by the selection to local segmentation of structures, and described partial structurtes are cut apart according to Gauss's mean parameter and covariance and provided the center of tubercle and the accurate estimation of ellipsoidal shape.The robustness of this estimation also allows non-tubercle zone (for example interested vascular and vascular point of crossing/branch) cut apart.

In order to simplify mathematical notation, in step 73, (VOI) carries out affine standardization to original volume of interest.This comprises makes the VOI bending convert the isotropy spheroid of fixed size to the anisotropy ellipsoid that will cut apart, and it is set at the center of VOI.Affine standardized parameter, that is to say that zoom direction and coefficient can be from by directly obtaining the Eigenvalue Analysis of cutting apart the estimated structure covariance of module.

Fig. 1 (a)-(c) illustrates exemplary according to an embodiment of the invention lung structure classification.Fig. 1 (a) represents original volume of interest (VOI) and the tubercle candidate of cutting apart, uses the nodular structure of ellipsoid match, here is vascular.The ellipsoid match obtains from the anisotropic Gaussian fitting module.The affine standardization of the VOI that Fig. 1 (b) expression is original, wherein ellipsoid is bent to the isotropy spheroid.Fig. 1 (c) is illustrated in apart from r _BoundThe border manifold of structure is cut in punishment, and it is launched into the 2D image and by spherical polar coordinates θ and  parametrization.Obtain the gradation of image value by the trilinearly interpolated method.

Refer again to Fig. 7, in step 74, determine the kind of lung structure type in the research by suitable stream shape being carried out cluster analysis, described suitable stream shape is calculated from the borderline region of target structure.Spherical flow shape makes up from affine standardized 3D rendering according to an embodiment of the invention.Geometrically, purpose is to represent to exceed a little the spherical layer of target structure border surface, makes to contain information relevant for the outstanding object (object) by the surface.The supposition in initial VOI of its shape is an ellipsoid, and is especially proportional with the ellipsoid that is obtained from cut apart based on Gauss's anisotropy.Therefore, in affine standardized expression, also with by central point (a _Bound, b _Bound) and radius r _BoundThe isotropy ball shape that limits is corresponding.Although central point is identical with the ellipsoidal central point of cutting apart, radius of a ball r _BoundWill determine with the data-driven version of being explained hereinafter.

Suppose a fixing r _Bound, border manifold represent can from Descartes (x, y, z) be converted into spherical coordinates (θ, ).Here, θ refers to the position angle, and  refers to polar angle.The result be affine standardized ellipsoid launch expression be 2D image array I (θ, ).Figure (1c) illustrates the result of lung structure example.It should be noted that opposite with general custom, polar angle comprises Interval at interval _=2 π (replacing π), promptly  ∈ [π, π] causes the double appearance of Descartes's voxel.The reason of introducing this redundancy be the following cluster that will introduce aspect two parameters at its corresponding Interval at interval _θAnd Interval _On require I (θ, cyclic behaviour ), just, I (θ+Interval _θ, +Interval _)=I (θ, ).Under the spherical coordinates situation, if Interval _=π, then this obviously can not realize.

According to one embodiment of present invention, on the basis of the entropy of intensity distributions, use data-driven method to determine suitable radius r _BoundIn order to promote this method, consider Fig. 2 (a)-(f), every width of cloth figure all illustrate have different radii r (θ, the ellipsoid that ) launches in the territory represent, as shown in FIG..Fig. 2 (g) expression for radius r ∈ 1 ..., 32} is about corresponding image histogram entropy E that image intensity calculated _rTypical image entropy can basis

(θ ) calculates from image intensity I.After suitably standardization was carried out in distribution to image intensity value, the stream shape image of expansion was counted as the 2D likelihood function.Then, directly come the calculating strength entropy with the standardization intensity level that is interpreted as probable value.The selection of radius comprises automatically selectes radius, and it is the most special to make that the high strength group shows in corresponding stream shape owing to outstanding structure.This less have lower entropy with comparing than the image of long radius by several stream shape images of being formed of group as shown in Fig. 2 (d) with having, and its reason is following argument intuitively.Shown in Fig. 2 (a)-(b), small radii makes corresponding border ellipsoid pass target structure inside, causes having the high entropy of more flat likelihood.On the other hand, shown in Fig. 2 (e)-(f),, also cause high entropy than long radius owing to nigh other " non-target " structure occurs.Therefore, suitable radius r _BoundForm entropy distribution E _rThe part minimize.At this on the one hand, selected r _BoundBe positioned at positive difference coefficient

Appearance place for the first time, just,

$r_{\mathrm{bound}}=\underset{r}{\min }\left\{r\right|E_{r+1}>E_r\}.$

After converting the part of 3D lung structure to the 2D image, can use through abundant research, effectively and be easy to visual 2D image analysis technology.As can being seen from Fig. 1 (c), border manifold contains the lung structure valuable information of classifying.In fact, high strength group's quantity discloses the type of lung structure, and it is equal to the quantity of the outstanding object that passes through the border that limited.The classification of one embodiment of the invention is based upon on the basis of this observation, wants to carry out following territory hypothesis:

0 proximity structure that group's indication is not connected in the border manifold, therefore, the structure of cutting apart is corresponding with isolated tubercle;

1 group in the border manifold indicates the single connection to adhesion structure, and this stems under many circumstances and is attached to the more tubercle of macrostructure (for example lung wall etc.);

2 groups indicate two connections, and this is the most common for blood vessel; And

＞3 groups indicate the vascular point of crossing.

According to one embodiment of present invention, high strength group's quantity is discerned by cluster analysis, carries out in the step 75 of Fig. 7.The basis of the cluster strategy of one embodiment of the present of invention is based on the Gauss curve fitting of expectation maximization (EM).Except standard EM Gauss cluster characteristic, the clustering algorithm of one embodiment of the present of invention need be reflected in (θ, ) continuity that occurs at 2D border manifold edge of image place in the territory.Particularly, (θ ) is represented corresponding to directional data by parameterized border manifold by spherical angle variable.For the direction as shown data, referring to the simplicity of illustration of Fig. 3.(θ, ) suitable clustering algorithm should regain single group in the territory in direction.Yet, replacing the direction modeling with linear, each in three kinds of observable structures all will form independently group.And the clustering algorithm of one embodiment of the invention should be able to automatically be determined modulus.

Directional data can be visualized as hyperspherical lip-deep point with bidimensional on the circumference of circle.Fig. 6 (a)-(b) illustrates the example of directional data.The situation that produces directional data is as mentioned below.For circle variable θ, addition " a+b " becomes " (a+b) mod2 π ", here angle with at interval (π, π] expression.It should be noted, under the situation of this supposition, mod (asking mould) operator also be mapped to (π, π].Suppose variable μ _θ ^cAnd V _θ ^cThe round counterpart of expression mean value and variance.Using $\widetilde{\mathrm{\&theta;}}=(\mathrm{\&theta;}-v)\mathrm{<figure-callout\; id="2"\; label="mod"\; filenames="A20071008604100221.png,A20071008604100231.png"\; state="\{\{state\}\}">mod</figure-callout>}2\mathrm{\&pi;}$ Under the situation of the zero direction drift of expression, to μ _θ ^cAnd V _θ ^cReasonable definition should remain unchanged.The circle variable unchangeability should for

${\mathrm{\&mu;}}_{\widetilde{\mathrm{\&theta;}}}^c=({\mathrm{\&mu;}}_{\mathrm{\&theta;}}^c-v)\mathrm{<figure-callout\; id="2"\; label="mod"\; filenames="A20071008604100221.png,A20071008604100231.png"\; state="\{\{state\}\}">mod</figure-callout>}2\mathrm{\&pi;},$

$V_{\widetilde{\mathrm{\&theta;}}}^c=V_{\mathrm{\&theta;}}^c.$

Yet, can prove easily that from the example of Fig. 6 desired unchangeability is violated.Fig. 6 represents the simple collection that circle is observed, and it is corresponding to ν=-0.5 π, Θ in Fig. 6 (a)=0.1 π, and 0.25,0.6s}, and among Fig. 6 (b) $\widetilde{\mathrm{\&Theta;}}=\{0.6\mathrm{\&pi;},0.7\mathrm{\&pi;},-0.9\mathrm{\&pi;}\}$ Observe for these, the unbiasedness maximum likelihood (ML) of mean value and variance is estimated and can be calculated as ${\widehat{\mathrm{\&mu;}}}_{\mathrm{\&theta;}}=0.3\mathrm{\&pi;},$ ${\widehat{\mathrm{\&sigma;}}}_{\mathrm{\&theta;}}\&ap;0.26\mathrm{\&pi;},$ ${\widehat{\mathrm{\&mu;}}}_{\widetilde{\mathrm{\&theta;}}}\&ap;0.13\mathrm{\&pi;}$ With ${\widehat{\mathrm{\&sigma;}}}_{\widetilde{\mathrm{\&theta;}}}\&ap;0.90\mathrm{\&pi;},$ This has obviously violated the drift unchangeability.In the drawings, the value of mean value and variance uses the position of bullet and the length of subsidiary arc (accompanying arc) to illustrate respectively.Therefore, for the circle data, the linear definition of mean value and variance depends on zero direction greatly, and this is disproportionate behavior and requirement for suitable processing.

In order to handle this situation, suppose that round stochastic variable θ has PDF _p(θ).Under the situation of conformance with standard statistical property, PDF should satisfy p (θ) 〉=0 He $\underset{-\mathrm{\&pi;}}{\overset{\mathrm{\&pi;}}{\&Integral;}}p\left(\mathrm{\&theta;}\right)\mathrm{d\&theta;}=1$ Variable θ represents plural e ^{I θ}, and adopt and pass through ${\mathrm{\&rho;}}_{\mathrm{\&theta;}}\exp \left(i{\mathrm{\&mu;}}_{\mathrm{\&theta;}}^c\right)=E\left[\exp \left(\mathrm{i\&theta;}\right)\right]$ The round mean direction μ that limits _θ ^cWith circle variance V _θ ^cSymbol, wherein $V_{\mathrm{\&theta;}}^c=1-{\mathrm{\&rho;}}_{\mathrm{\&theta;}}.$ Parameter p _θBe called composition length.Figuratively speaking, μ _θ ^cBe e ^{I θ}Expect phase place and ρ _θBe e ^{I θ}Desired length. $V_{\mathrm{\&theta;}}^c\&Element;\left[\mathrm{0,1}\right]$ Measure the amount of deviation (dispersion).Can show that desired drift unchangeability is satisfied in these definition of mean value and variance, and can be used as the suitable counterpart of linear average and variance.

The statistical modeling that many models are used for directional data has been proposed.According to one embodiment of present invention, use multivariate to hold Gaussian distribution, it is the expansion that holds Gaussian distribution.The circumference of the circle of the Gaussian distribution N (x) of the variable x on the line can " hold (the wrapped) " unit's of wearing radius.Just, hold variable θ=χ _w=mod2 π ∈ (π, π] hold Gaussian distribution N _w(θ) be

$N_w\left(\mathrm{\&theta;}\right)=\underset{k=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}N(\mathrm{\&theta;}+2\mathrm{k\&pi;}).$

Vector variable Θ=(θ ₁..., θ _F) ^TMultivariate hold Gaussian distribution and can be defined as similarly:

$\left(1\right),N_w\left(\mathrm{\&Theta;}\right)=\underset{{\mathrm{<figure-callout\; id="1"\; label="k"\; filenames="A20071008604100231.png"\; state="\{\{state\}\}">k</figure-callout>}}_1=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\&CenterDot;\&CenterDot;\&CenterDot;\underset{k_F=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}N(\mathrm{\&Theta;}+2\mathrm{\&pi;}{k}_1{\mathrm{<figure-callout\; id="1"\; label="e"\; filenames="A20071008604100231.png"\; state="\{\{state\}\}">e</figure-callout>}}_1+\&CenterDot;\&CenterDot;\&CenterDot;+2\mathrm{\&pi;}{k}_F{e}_F),$

E wherein _k=(0 ..., 0,1,0 ..., 0) ^TBe k Euclid's basis vector, k item element be designated as 1 and other place be designated as 0.Fig. 3 illustrates the example that the bidimensional multivariate holds Gauss.

Show that the suitable little variance in the assigned direction variable can be utilized

$\left(2\right),{\left({\widehat{\mathrm{\&mu;}}}_{\mathrm{\&theta;}}\right)}_f=\arg \left(\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\exp \left(i{\mathrm{\&theta;}}_f^{\left(m\right)}\right)\right)$

With

$\left(3\right),{\widehat{\mathrm{\&Sigma;}}}_{\mathrm{\&theta;}}=\frac{1}{M-1}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&Theta;}}^{{\left(m\right)}^{\&prime;}}{\mathrm{\&Theta;}}^{{\left(m\right)}^{\&prime;T}}$

From sample set X={ θ ⁽¹⁾..., θ ^(M)Middle accurate mean value and the covariance estimation that obtains equation (1)

With

, wherein ${\mathrm{\&theta;}}^{{\left(m\right)}^{\&prime;}}=({\mathrm{\&theta;}}_f^{\left(m\right)}-{\left({\widehat{\mathrm{\&mu;}}}_{\mathrm{\&theta;}}\right)}_f)\mathrm{<figure-callout\; id="2"\; label="mod"\; filenames="A20071008604100221.png,A20071008604100231.png"\; state="\{\{state\}\}">mod</figure-callout>}2\mathrm{\&pi;}$ , i ²=-1, and " arg " is the phase place of plural number.In order to simplify, suppose cycle and the θ of 2 π for all the dimension f among the Θ all recessively _f∈ (π, π] scope.

Expectation maximization (EM) algorithm is a class statistical method that is used for finding at probability statistics model the maximal possibility estimation of parameter, and wherein model depends on unobservable potential variable.EM replaces between carry out desired (E) step and maximization (M) step, the wherein expectation of desired step by comprising that potential variable calculates likelihood as their observed arrivals, maximization steps will be by maximizing the maximal possibility estimation of calculating parameter in the expectation likelihood that the E step is found.The parameter that finds in step M then is used to begin another E step, and repeats said process.The EM algorithm will improve initial estimation θ in the mode of iteration ₀The estimation θ new with structure ₁..., θ _n

If y represents the fragmentary data be made up of the value of observable variable, and x represents the data that lack, and x and y lump together and just form a complete data collection so.Suppose p be have by vectorial θ the joint probability distribution function of complete data of given parameter: p (y, x| θ).This function provides the complete data likelihood.Then, utilize the Bayes rule, what the condition of given unobservable variable distributed is contemplated to be

$p\left(x\right|y,\mathrm{\&theta;})=\frac{p(y,x|\mathrm{\&theta;})}{p\left(y\right|\mathrm{\&theta;})}=\frac{p\left(y\right|x,\mathrm{\&theta;})p\left(x\right|\mathrm{\&theta;})}{\&Integral;p\left(y\right|\hat{x},\mathrm{\&theta;})p\left(\hat{x}\right|\mathrm{\&theta;})d\hat{x}}.$

This formulate only require about given unobservable data observation likelihood p (y|x, θ) and the knowledge of the Probability p of unobservable data (x| θ).From θ _nObtain θ _N+1Independent maximization steps be:

${\mathrm{\&theta;}}_{n=1}=\arg \underset{\mathrm{\&theta;}}{\max }{E}_x\left[\log p(y,x|\mathrm{\&theta;})\right|y],$

E wherein _xθ was fixed on θ during the condition that [] is illustrated in x distributed _nThe time log p (y, x| θ) conditional expectation.Log-likelihood log p (y, x| θ) often is used to substitute real likelihood p (y, x| θ), and its reason is that it has produced the formula that is more prone to, but still obtains its maximal value at the some place identical with likelihood.In other words, given viewed variable under the parameter value situation formerly, θ _N+1It is the value of conditional expectation (E) maximization (M) that makes the complete data log-likelihood.Usually, by forming Lagrange (Lagrangian) function of log-likelihood, obtain maximal value with respect to mean value and covariance estimation derivative then.

In the EM clustering algorithm, rule, linear Gauss model can hold Gauss model with above-mentioned multivariate and replace.Especially, EQ. (1) and EQ. (2) and (3) alternative initial linear equivalent in E and M step respectively on the other hand on the one hand.

According to one embodiment of the invention, can obtain the model that modulus is selected from the limited mixing of adopting Gauss's estimation and the correction EM clustering algorithm of following the Model Selection of minimum description length (MDL) criterion, use so that the quantity of the component in mixing minimizes.Usually, the input of EM clustering algorithm is the sample set X={ (θ that observes ₁,  ₁) ..., (θ _M,  _M), and current data be 2D (image) matrix I (θ, ).In order to overcome this incompatibility, directly (θ is extracted in ), wherein (the θ of each sampling from I to observe X _m,  _m) ∈ (π, π] * (π, π] occurrence number and corresponding image array value I (θ _m,  _m) be set pro rata.

One of shape of the outstanding structure of reality in border manifold and oval gaussian shape be at once, takes place with an aspect of Gauss EM cluster correlation.In this case, the match of expectation EM algorithm has this structure of one group of gaussian component.This effect will obviously have a negative impact to classification, and wherein the quantity of component has played integral action.

Refer again to Fig. 7, according to one embodiment of present invention,, use aftertreatment and merge suitable component in step 76.Especially, this aftertreatment of one embodiment of the present of invention can be counted as second cluster analysis, and described second cluster analysis is merged into single group to the gaussian component collective row analysis of all EM matches and with subclass, up to a certain yardstick.The well known in the art a kind of technology that is used for this situation is a cohesion hierarchical clustering (agglomerative hierarchical clustering).In hierarchical clustering, the cluster space is represented with the distance of its element.In this case, element is that multivariate holds Gaussian function, and the descriptive statistics symbol is used to geometric configuration.Suitable (with computable analysis) of Gaussian distribution, statistical distance was measured is the Bhattacharyya distance

$D_{\mathrm{Bhatt}}({\mathrm{\&mu;}}_1,{\mathrm{\&Sigma;}}_1,{\mathrm{\&mu;}}_2,{\mathrm{\&Sigma;}}_2)=\frac{1}{8}{({\mathrm{\&mu;}}_2-{\mathrm{\&mu;}}_1)}^T{\left(\frac{{\mathrm{\&Sigma;}}_1+{\mathrm{\&Sigma;}}_2}{2}\right)}^{-1}({\mathrm{\&mu;}}_2-{\mathrm{\&mu;}}_1)+\frac{1}{2}1n\frac{|{\mathrm{\&Sigma;}}_1+{\mathrm{\&Sigma;}}_2|}{\sqrt{\left|{\mathrm{\&Sigma;}}_1\right|\left|{\mathrm{\&Sigma;}}_2\right|}}.$

Yet, D _BhattThere is not to consider to hold Gauss's direction character.Thereby, D has been proposed according to one embodiment of present invention _BhattThe correction modification, " hold Bhattacharyya distance ":

$D_{\mathrm{Bhatt}}^w({\mathrm{\&mu;}}_1,{\mathrm{\&Sigma;}}_1,{\mathrm{\&mu;}}_2,{\mathrm{\&Sigma;}}_2)=\frac{1}{8}{\left(({\mathrm{\&mu;}}_2-{\mathrm{\&mu;}}_1)\mathrm{<figure-callout\; id="2"\; label="mod"\; filenames="A20071008604100221.png,A20071008604100231.png"\; state="\{\{state\}\}">mod</figure-callout>}2\mathrm{\&pi;}\right)}^T{\left(\frac{{\mathrm{\&Sigma;}}_1+{\mathrm{\&Sigma;}}_2}{2}\right)}^{-1}\left(({\mathrm{\&mu;}}_2-{\mathrm{\&mu;}}_1)\mathrm{<figure-callout\; id="2"\; label="mod"\; filenames="A20071008604100221.png,A20071008604100231.png"\; state="\{\{state\}\}">mod</figure-callout>}2\mathrm{\&pi;}\right)+\frac{1}{2}1n\frac{|{\mathrm{\&Sigma;}}_1+{\mathrm{\&Sigma;}}_2|}{\sqrt{\left|{\mathrm{\&Sigma;}}_1\right|\left|{\mathrm{\&Sigma;}}_2\right|}}.$

At last, the quantity that holds the gaussian component group has determined the class of lung structure: the isolated tubercle of 0 representative, and 2 * 1=2 represents the adhesion tubercle, and 2 * 2=4 represents vascular, and＞2 * 3=6 represents fat pipe point of crossing.It is owing to the interval of the twice in polar coordinates , as discussed above that coefficient gets 2.

The limitation of this method of the embodiment of the invention is the following fact, promptly in that (θ, ) yardstick depends on the position in the territory.A replacement scheme according to an embodiment of the invention will be: distribute with von Mises-Fisher and can avoid this problem to the direction data modeling.Form is taked in vonMises-Fisher distribution for d dimension unit random vector:

f(x|μ，k)＝c _d(k)exp(kμ ^Tx)，

‖ μ ‖=1 wherein, k 〉=0, and d 〉=2.Generalized constant c _d(k) by

$c_d\left(k\right)=\frac{k^{d/2-1}}{{\left(2\mathrm{\&pi;}\right)}^{d/2}{I}_{d/2-1}{\left(k\right)}^{\&prime;}}$

Provide, wherein I _r() expression r rank first kind correction Bessel (Bessel) function.Distribute and come parametrization with mean direction (mean direction) μ and lumped parameter (concentration parameter) k, it characterizes (x| μ, k) vector of unit length that is extracted is how to concentrate on forcefully around the mean direction μ according to f.Yet the parameter estimation that von Mises-Fisher distributes generally includes to be found the solution the implicit equation of Bessel (Bessel) function ratio, and there are not analytic solution in it.

According to one embodiment of present invention, carry out qualitative test for the lung structure classification that is proposed.Figure 4 and 5 have been showed the classification diagram of chest CT image, all have two examples to represent to each class " tubercle ", " adhesion tubercle ", " vascular ", " vascular point of crossing ".

Fig. 4 (a)-(d) and 5 (a)-(d) have described the illustrative example of the lung structure sorting technique of the one embodiment of the invention that is used for Thoracic CT scan.Each row is cut apart and is confirmed corresponding with an example.About being adhered to the tubercle of lung wall, and the row 1 and 2 of Fig. 5 has shown blood vessel structure to preceding two row of Fig. 4 about tubercle object, back two row, and row 3 and 4 is vascular point of crossing.Row (a) illustrate the CT VOI in three quadrature xsects.Segmentation result illustrates with ellipsoid.Row (b) have been represented the affine standardization of initial VOI, make the 3D ellipsoid be bent to spheroid.Row (c) have shown in (θ, the constructed border manifold of ) launching in the territory.Yet be noted that and introduced the added strength threshold value.Adopt this step as quick, the simple and easy means of eliminating the low-intensity structure, described low-intensity structure may be disturbed Gauss EM cluster.The graphic presentation of row in (d) by result based on the algorithm fitted Gaussian mixture model of EM.Dotted ellipse shape is corresponding with the cluster gaussian component based on EM, and the solid line ellipse has been described the group through aftertreatment.

As shown in row (a), 3D is cut apart the tubercle of the isolated and adhesion of divided ownership as shown in Figure 4, also can cut apart false positive blood vessel and vascular point of crossing as shown in Figure 5.In row (d), the border manifold image is converted into sampled data collection X, as described in the 2.2.1 part.In addition, row (d) shown based on EM hold Gauss's clustering result, just, the mean value of k component and covariance with dashed lines are oval to be represented.Be noted that (θ, ) continuity of the boundary in territory in the row 3,4 of Fig. 4 and row 3,4 especially at Fig. 5.For visual purpose, comprised the diagram of hierarchical clustering aftertreatment.The group who obtains from this aftertreatment uses k ₂The oval expression of solid line, its central point and expansion are corresponding to holding mean value and covariance that Gauss's mean value calculation draw from all at one in the group of aftertreatment.What note is that if group's radix is lower, then this diagram can form the ellipse of degeneration, for example in the row 2 of Fig. 5.From number of components k ₂Infer the structure class, can verify that current classification has provided correct option to all 8 examples.Other situation has obtained similar results.

It is worthy of note the limitation of this classification, this can cause misclassification in some cases.(in the θ-Wei of 2D image, become big unworthily after launching corresponding to the structure at the limit place of θ=0 and =π) at stream shape 3D spheroid.This situation can be compared with the phenomenon from cartography, and wherein the north and south poles zone occupies bigger zone at ratio on 2D Mercator (Mercator) the projection world map on the spherical world globe of 3D.In the illustrated embodiment, this behavior can be observed in the row 4 of Fig. 5 in the above, wherein the high-strength structure of  ≈ π spread all over θ gamut (π, π].Thereby suggestion is careful when especially being reached a conclusion in these limit zones from the scaling relation of the stream shape of launching.In fact this is the shortcoming that holds Gauss's modeling, especially launches.In this, it should be noted that this phenomenon has been avoided in von Mises-Fisher modeling, because hypothesis is not launched.

It should be understood that the present invention can realize with the process of hardware, software, firmware, specific purposes and the various forms of combination thereof.In one embodiment, the present invention can be implemented as the tangible application program that is embodied on the computer-readable program memory device with form of software.This application program can be uploaded to and be carried out by the machine that comprises any suitable structure.

Fig. 8 is the block diagram that is used to carry out the typical computer system of sorting technique according to an embodiment of the invention.Referring now to Fig. 8,, is used to realize that computer system 81 of the present invention especially can comprise CPU (central processing unit) (CPU) 82, storer 83 and I/O (I/O) interface 84.Computer system 81 is general by I/O interface 84 and display 85 and various input equipment 86 (for example mouse and keyboard) coupling.Support circuit can comprise circuit as cache memory, power supply, clock circuit and communication bus and so on.Storer 83 can comprise random access memory (RAM), ROM (read-only memory) (ROM), disc driver, tape drive etc. and combination thereof.The present invention can be used as the routine 87 that is stored in the storer 83 and realizes, and carries out to handle the signal from signal source 88 by CPU82.Similarly, computer system 81 is general-purpose computing systems, becomes dedicated computer system when carrying out routine 87 of the present invention.

Computer system 81 also comprises operating system and micro-instruction code.Various process as described herein and function can be by the part of the micro-instruction code of operating system execution or the part of application program (or and combination).In addition, other various peripherals can be connected to computer platform, for example auxiliary data storage device and printing device.

Will be further understood that so the mode that is programmed according to the present invention, the actual connection between the system unit (or process steps) can be different because more described in the accompanying drawings composition system units and method step can be realized with form of software.Given instruction of the present invention provided here, those of ordinary skill in the related art can imagine these and similarly embodiment or configuration of the present invention.

Although the present invention is described in detail with reference to preferred embodiment, those skilled in the art recognize because of this, in not departing from the spirit and scope of the present invention of setting forth as appending claims, can carry out various modifications and replacement to it.

Claims (24)

1, a kind of method of lung structure being classified at digitized image of being used for may further comprise the steps:

The approximate target structure position of one or more target structures is provided in digitizing three-dimensional (3D) image;

Match anisotropic Gaussian model around described approximate target position is to produce the center of more accurate 3D target model and described one or more target structures;

Make each described 3D target model bend to the 3D spheroid;

Around each described bending 3D spheroid, make up border manifold; And

The identification group wherein classifies to described one or more target structures on described border manifold.

2, the method for claim 1, wherein said digitized image comprise with 3 d grid on the corresponding a plurality of intensity in some territory.

3, the method for claim 1, wherein the match anisotropic Gaussian model comprises around approximate target position: use Gauss's metric space average drifting analysis and select to produce the 3D ellipsoid model of described target structure based on the Jansen-Shannon auto-bandwidth of divergence, the ellipsoidal center of wherein said 3D is corresponding with the center and the covariance of described Gauss model with size.

4, method as claimed in claim 3, wherein crooked described 3D target model comprises the described 3D ellipsoid of affine standardization, wherein zoom direction and coefficient are to obtain from the structure covariance of described anisotropic Gaussian model.

5, the method for claim 1 wherein makes up border manifold and comprises that further the 3D unfolded surface with bending spheroid becomes 2D to represent, and determines the radius of suitable border manifold.

6, method as claimed in claim 5, wherein the 3D unfolded surface of bending spheroid is become 2D represent to comprise surface transformation with described bending spheroid become spherical coordinates (θ, ), wherein  ∈ [π, π] and θ ∈ [π, π].

7, method as claimed in claim 5 determines that wherein the radius of suitable border manifold comprises: the spherical flow shape that makes up a plurality of different radiis around described bending spheroid; Each spherical flow shape is launched into 2D to be represented; To the distribution carrying out of the intensity level on the spherical flow shape of each described expansion standardization; To the spherical flow shape calculating strength entropy of each described expansion, wherein intensity level is counted as probable value, and wherein entropy distributes and is defined; And find and make the described entropy minimized radius that distributes.

8, the method for claim 1, wherein the identification group comprises that the use expectation maximization makes vector variable Θ=(θ ₁..., θ _F) ^TMultivariate hold Gaussian distribution N _w ^pMixing (Θ) $N_w\left(\mathrm{\&Theta;}\right)=\underset{p=1}{\overset{p}{\mathrm{\&Sigma;}}}{c}_p{N}_w^p\left(\mathrm{\&Theta;}\right)$ Be fit to by the outstanding object of described border manifold, described border manifold is followed the minimum description length criterion; Mixed components probability c wherein _pIn expectation maximization, estimated, wherein each the dimension in, θ _iSatisfy θ=x _w=xmod 2 π ∈ (π, π], N _w ^p(Θ) satisfy $N_w^p\left(\mathrm{\&Theta;}\right)=\underset{k_1=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\&CenterDot;\&CenterDot;\&CenterDot;\underset{k_F=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{N}^p(\mathrm{\&Theta;}+2\mathrm{\&pi;}{k}_1{e}_1+...+2\mathrm{\&pi;}{k}_F{e}_F)$ , e wherein _k=(0 ..., 0,1,0 ..., 0) ^TBe k Euclid's basis vector, the k item be designated as 1 and other place be designated as 0, the wherein estimation of mixed components p

With

Be from sample set X={ θ in expectation maximization ⁽¹⁾..., θ ^(M)In based on direction mean value ${\left({\widehat{\mathrm{\&mu;}}}_{\mathrm{\&theta;}}\right)}_f=\arg \left(\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\exp \left(i{\mathrm{\&theta;}}_f^{\left(m\right)}\right)\right)$ And covariance ${\widehat{\mathrm{\&Sigma;}}}_{\mathrm{\&theta;}}=\frac{1}{M-1}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&Theta;}}^{{\left(m\right)}^{\&prime;}}{\mathrm{\&Theta;}}^{{\left(m\right)}^{\&prime;T}}$ Obtain, wherein ${\mathrm{\&theta;}}^{{\left(m\right)}^{\&prime;}}=({\mathrm{\&theta;}}_f^{\left(m\right)}-{\left({\widehat{\mathrm{\&mu;}}}_{\mathrm{\&theta;}}\right)}_f)\mathrm{mod}2\mathrm{\&pi;},$ And directly (θ ) extracts, wherein (the θ of each sampling from 2D unfolded image I wherein to observe X _m,  _m) ∈ (π, π] * (π, π] occurrence number and corresponding image array value I (θ _m,  _m) be set pro rata.

9, the method for claim 1 wherein further comprises and uses the cohesion hierarchical clustering that the group in the mutual preset distance is merged, and is used for multivariate and holds right the equaling of Gaussian distribution $\frac{1}{8}{\left(({\mathrm{\&mu;}}_2-{\mathrm{\&mu;}}_1)\mathrm{mod}2\mathrm{\&pi;}\right)}^T{\left(\frac{{\mathrm{\&Sigma;}}_1+{\mathrm{\&Sigma;}}_2}{2}\right)}^{-1}$ $\left(({\mathrm{\&mu;}}_2-{\mathrm{\&mu;}}_1)\mathrm{mod}2\mathrm{\&pi;}\right)+\frac{1}{2}\ln \frac{|{\mathrm{\&Sigma;}}_1+{\mathrm{\&Sigma;}}_2|}{\sqrt{\left|{\mathrm{\&Sigma;}}_1\right|\left|{\mathrm{\&Sigma;}}_2\right|}}$ Distance metric, μ wherein ₁And μ ₂Be the right mean value of Gaussian distribution, and ∑ ₁And ∑ ₂It is its variance separately.

10, the method for claim 1, wherein the class of lung structure is determined by the quantity that holds the gaussian component group relevant with target structure, and wherein Gu Li tubercle has 0 group, and the adhesion tubercle has 2 groups, and vascular has 4 groups, and there are 6 or multigroup more in the vascular point of crossing.

11, a kind of method of lung structure being classified at digitized image of being used for may further comprise the steps:

The target position of one or more 3D spheroids is provided in digitizing three-dimensional (3D) image, described image comprise with 3 d grid on the corresponding a plurality of intensity in some territory, each 3D spheroid is all represented target structure in described image;

Around described 3D spheroid, make up the spherical flow shape of a plurality of different radiis;

To each described spherical flow shape calculating strength entropy, wherein intensity level is counted as probable value, and wherein entropy distributes and is defined;

Find to make the described entropy minimized radius that distributes, wherein saidly minimize radius limited boundary stream shape;

The unfolded surface of border manifold is become 2D spherical coordinates (θ, ) expression, wherein  ∈ [π, π] and θ ∈ [π, π];

Use expectation maximization to come match vector variable Θ=(θ ₁..., θ _F) ^TMultivariate hold Gaussian distribution N _wMixing (Θ) $N_w\left(\mathrm{\&Theta;}\right)=\underset{p=1}{\overset{p}{\mathrm{\&Sigma;}}}{c}_p{N}_w^p\left(\mathrm{\&Theta;}\right)$ , mixed components probability c wherein _pIn expectation maximization, estimated, wherein θ _i=(θ _i,  _i) be group by the outstanding target structure of described border manifold, and wherein lung structure is classified by a plurality of outstanding groups.

12, method as claimed in claim 11 further comprises the intensity distributions standardization that makes on each of described a plurality of spherical flow shapes.

13, method as claimed in claim 11 wherein provides the target position of one or more 3D spheroids to comprise: the approximate target structure position that described one or more target structures are provided in described digitizing three-dimensional (3D) image; The center that the match anisotropic Gaussian model produces more accurate 3D ellipsoid target model and described ellipsoid model around described approximate target position; And with the affine 3D spheroid that is standardized into of described ellipsoid model.

14, method as claimed in claim 11, wherein each the dimension in, θ _iSatisfy θ=x _w=x mod 2 π ∈ (π, π], N _w(Θ) satisfy $N_w\left(\mathrm{\&Theta;}\right)=\underset{k_1=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\&CenterDot;\&CenterDot;\&CenterDot;\underset{k_F=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}N(\mathrm{\&Theta;}+2\mathrm{\&pi;}{k}_1{e}_1+\&CenterDot;\&CenterDot;\&CenterDot;+2\mathrm{\&pi;}{k}_F{e}_F)$ , e wherein _k=(0 ..., 0,1,0 ..., 0) ^TBe k Euclid's basis vector, the k item be designated as 1 and other place be designated as 0, the wherein estimation of mixed components p

With

Be from sample set X={ θ in expectation maximization ⁽¹⁾..., θ ^(M)In based on direction mean value ${\left({\widehat{\mathrm{\&mu;}}}_{\mathrm{\&theta;}}\right)}_f=\arg \left(\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\exp \left(i{\mathrm{\&theta;}}_f^{\left(m\right)}\right)\right)$ And covariance ${\widehat{\mathrm{\&Sigma;}}}_{\mathrm{\&theta;}}=\frac{1}{M-1}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&Theta;}}^{{\left(m\right)}^{\&prime;}}{\mathrm{\&Theta;}}^{{\left(m\right)}^{\&prime;T}}$ Obtain, wherein ${\mathrm{\&theta;}}^{{\left(m\right)}^{\&prime;}}=({\mathrm{\&theta;}}_f^{\left(m\right)}-{\left({\widehat{\mathrm{\&mu;}}}_{\mathrm{\&theta;}}\right)}_f)\mathrm{mod}2\mathrm{\&pi;}$ , and directly (θ ) extracts, wherein (the θ of each sampling from 2D unfolded image I wherein to observe X _m,  _m) ∈ (π, π] * (π, π] occurrence number and corresponding image array value I (θ _m,  _m) be set pro rata.

15, computer-readable program storage device visibly embodies by the executable instruction repertorie of computing machine, implements to be used for the method step of lung structure being classified at digitized image, said method comprising the steps of:

The approximate target structure position of one or more target structures is provided in digitizing three-dimensional (3D) image;

The center that the match anisotropic Gaussian model produces more accurate 3D target model and described one or more target structures around described approximate target position;

Make each described 3D target model bend to the 3D spheroid;

Around each described bending 3D spheroid, make up border manifold; And

The identification group wherein classifies to described one or more target structures on described border manifold.

16, computer-readable program storage device as claimed in claim 15, wherein said digitized image comprise with 3 d grid on the corresponding a plurality of intensity in some territory.

17, computer-readable program storage device as claimed in claim 15, wherein the match anisotropic Gaussian model comprises around approximate target position: use Gauss's metric space average drifting analysis and select to produce the 3D ellipsoid model of described target structure based on the Jansen-Shannon auto-bandwidth of divergence, the ellipsoidal center of wherein said 3D is corresponding with the center and the covariance of described Gauss model with size.

18, computer-readable program storage device as claimed in claim 17, wherein crooked described 3D target model comprises the described 3D ellipsoid of affine standardization, wherein zoom direction and coefficient are to obtain from the structure covariance of described anisotropic Gaussian model.

19, computer-readable program storage device as claimed in claim 15 wherein makes up border manifold and comprises that further the 3D unfolded surface with bending spheroid becomes 2D to represent, and determines the radius of suitable border manifold.

20, computer-readable program storage device as claimed in claim 19, wherein the 3D unfolded surface of bending spheroid is become 2D represent to comprise surface transformation with described bending spheroid become spherical coordinates (θ, ),  ∈ [π wherein, π] and θ ∈ [π, π].

21, computer-readable program storage device as claimed in claim 19 determines that wherein the radius of suitable border manifold comprises: the spherical flow shape that makes up a plurality of different radiis around described bending spheroid; Each spherical flow shape is launched into 2D to be represented; Will be in the intensity level distribution standardization on each described expansion spherical flow shape; To the spherical flow shape calculating strength entropy of each described expansion, wherein intensity level is counted as probable value, and wherein entropy distributes and is defined; And find and make the described entropy minimized radius that distributes.

22, computer-readable program storage device as claimed in claim 15, wherein

The identification group comprises makes vector variable Θ=(θ ₁..., θ _F) ^TMultivariate hold Gaussian distribution N _w ^pMixing (Θ) $N_w\left(\mathrm{\&Theta;}\right)=\underset{p=1}{\overset{p}{\mathrm{\&Sigma;}}}{c}_p{N}_w^p\left(\mathrm{\&Theta;}\right)$ Be fit to by the outstanding object of described border manifold, described border manifold is followed the minimum description length criterion; Mixed components probability c wherein _pIn expectation maximization, estimated, wherein each the dimension in, θ _iSatisfy θ=x _w=x mod 2 π ∈ (π, π], N _w ^p(Θ) satisfy $N_w^p\left(\mathrm{\&Theta;}\right)=\underset{k_1=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\&CenterDot;\&CenterDot;\&CenterDot;\underset{k_F=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{N}^p(\mathrm{\&Theta;}+2\mathrm{\&pi;}{k}_1{e}_1+...+2\mathrm{\&pi;}{k}_F{e}_F)$ , e wherein _k=(0 ..., 0,1,0 ..., 0) ^TBe k Euclid's basis vector, the k item be designated as 1 and other place be designated as 0, the wherein estimation of mixed components p

With

Be from sample set X={ θ in expectation maximization ⁽¹⁾..., θ ^(M)In based on direction mean value ${\left({\widehat{\mathrm{\&mu;}}}_{\mathrm{\&theta;}}\right)}_f=\arg \left(\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\exp \left(i{\mathrm{\&theta;}}_f^{\left(m\right)}\right)\right)$ And covariance ${\widehat{\mathrm{\&Sigma;}}}_{\mathrm{\&theta;}}=\frac{1}{M-1}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&Theta;}}^{{\left(m\right)}^{\&prime;}}{\mathrm{\&Theta;}}^{{\left(m\right)}^{\&prime;T}}$ Obtain, wherein ${\mathrm{\&theta;}}^{{\left(m\right)}^{\&prime;}}=({\mathrm{\&theta;}}_f^{\left(m\right)}-{\left({\widehat{\mathrm{\&mu;}}}_{\mathrm{\&theta;}}\right)}_f)\mathrm{mod}2\mathrm{\&pi;},$ And directly (θ ) extracts, wherein (the θ of each sampling from 2D unfolded image I wherein to observe X _m,  _m) ∈ (π, π] * (π, π] occurrence number and corresponding image array value I (θ _m,  _m) be set pro rata.

23, computer-readable program storage device as claimed in claim 15, described method comprise that further use cohesion hierarchical clustering merges the group in the mutual preset distance, is used for multivariate and holds right the equaling of Gaussian distribution $\frac{1}{8}{\left(({\mathrm{\&mu;}}_2-{\mathrm{\&mu;}}_1)\mathrm{mod}2\mathrm{\&pi;}\right)}^T{\left(\frac{{\mathrm{\&Sigma;}}_1+{\mathrm{\&Sigma;}}_2}{2}\right)}^{-1}$ $\left(({\mathrm{\&mu;}}_2-{\mathrm{\&mu;}}_1)\mathrm{mod}2\mathrm{\&pi;}\right)+\frac{1}{2}\ln \frac{|{\mathrm{\&Sigma;}}_1+{\mathrm{\&Sigma;}}_2|}{\sqrt{\left|{\mathrm{\&Sigma;}}_1\right|\left|{\mathrm{\&Sigma;}}_2\right|}}$ Distance metric, μ wherein ₁And μ ₂Be the right mean value of Gaussian distribution, and ∑ ₁And ∑ ₂It is its variance separately.

24, computer-readable program storage device as claimed in claim 15, wherein the class of lung structure determines that by the quantity that holds the gaussian component group relevant with target structure wherein Gu Li tubercle has 0 group, and the adhesion tubercle has 2 groups, vascular has 4 groups, and there are 6 or multigroup more in the vascular point of crossing.

Images