CN103093503B - Based on the method for building up of the pulmonary parenchyma region surface model of CT image - Google Patents

Abstract

The present invention relates to a kind of method for building up of the pulmonary parenchyma surface model based on CT image, the method first utilizes often the open former figure of computing machine to CT lung scans image to carry out Gaussian smoothing operation, to eliminate image noise interference; Re-use three-dimensional Fast Marching algorithm and obtain pulmonary parenchyma region; Use Marching Cubes algorithm; The calculating of recycling surface curvature, determines the concrete form feature of pulmonary parenchyma, and removes the regions of non-interest determined, repairs, obtain the surface model of a smooth surface and complete lung's area-of-interest finally by marginalisation.The invention provides a kind of method of area of computer aided diagnostic imaging, the diagnosis difference that the difference due to subjective factors such as doctor's subjective experience, observation abilitys causes can be eliminated, and provide the higher reference identifying and diagnosing result of accuracy rate, thus diagnostic imaging is objectified more, improve efficiency and the accuracy of diagnosis.

Description

Based on the method for building up of the pulmonary parenchyma region surface model of CT image

Technical field

The present invention relates to a kind of disposal route of medical image, specifically a kind of method for building up of the pulmonary parenchyma region surface model based on CT image.

Background technology

Lung cancer is common primary lung tumors, and M & M rises year by year, becomes the malignant tumour that in global range, M & M is the highest at present.According to the investigation result display of the Ministry of Public Health in 2006, the lung cancer mortality of China increases obviously, has become the primary cancer cause of the death, has comparatively improved 75.77% the nineties, gets rid of the impact of age composition change, too increases 33.25%.The clinical manifestation of lung cancer is varied, and the early stage non-evident sympton of Most patients and sign, not easily find, during discovery, often arrive late period.If lung cancer in early days by Diagnosis and Treat, just correspondingly can improve cure rate and the survival rate of patient.

CT(Computer Tomography) computed tomography is technology the most frequently used in chest imaging, is widely used in chest pulmonary disease checkout and diagnosis.Lung cancer is show with the form of Lung neoplasm usually in CT image, and the various shapes of Lung neoplasm, not of uniform size, distributing position is not also fixed, and is easily connected with other tissue tights; Density and the pulmonary vascular of lung cancer are close, mostly show as the fine and close spot of circle or sub-circular in CT image, are only difficult to be distinguished with human eye.And in nearly hundred CT images that CT lung scans produces, tuberculous image only may have several.And the work analyzed all image datas is very uninteresting loaded down with trivial details, even if experienced doctor also can fail to pinpoint a disease in diagnosis or mistaken diagnosis unavoidably.Thus computer-aided diagnosis (Computer-Aided Diagnosis) has just become the important method addressed this problem.As a part for computer-aided diagnosis, to splitting fast and accurately of area-of-interest (Region ofInterest), be the primary link of work, be also an important foundation of carrying out next step diagnosis and detection classification simultaneously.

Doctor is when diagnosing by the disease of Thoracic CT scan image to patient, and usually interested is the region such as pulmonary parenchyma and lung inner disease foci, and namely this part region is called " area-of-interest " on CT scan image.Being partitioned into the area-of-interest in CT scan image exactly, to identification and the classification of lung tissue, and is all vital to the diagnosis etc. of lung cancer.But, also do not have a kind of computing machine that utilizes to carry out the method for building up of surface model to the lung's area-of-interest in CT scan image at present.

Summary of the invention

Object of the present invention is just to provide a kind of method for building up of the pulmonary parenchyma region surface model based on CT image, to solve the problem can not carrying out microcomputer modelling to the pulmonary parenchyma region in CT scan image, to meet the use needs of clinical research, Diagnosis and Treat.

The object of the present invention is achieved like this: a kind of method for building up of the pulmonary parenchyma region surface model based on CT image, comprises the following steps:

A, the smoothing operation of every pictures utilizing computing machine to sort in former figure to of CT lung scans group, to eliminate image noise interference;

B, acquisition pulmonary parenchyma region: on the former figure of this group after denoising, obtain initial pulmonary parenchyma region, recycle three-dimensional FastMarching algorithm, three-dimensional segmentation is carried out to the pulmonary parenchyma region obtained, and segmentation result is covered Hui Yuantu, obtain the original gray value in described pulmonary parenchyma region, and arrange by former figure sequence number;

C, surface modeling is carried out to pulmonary parenchyma region: utilize Marching Cubes algorithm, sequence picture b step handled well reads in calculator memory as 3 d-dem regular data field, to each cube (Cube) computing one by one in described data fields, draw equivalent point identical with described original gray value on each limit of each cube by linear interpolation, then simulate the contour surface in this cube with tri patch; All cubical contour surface in whole described data fields forms the initial surface model in described pulmonary parenchyma region;

D, determine the surface curvature of pulmonary parenchyma region initial surface model each point: the surface curvature calculating three vertex position places in the initial surface model of described pulmonary parenchyma region on each tri patch respectively, and the curvature value of described pulmonary parenchyma region initial surface model is greater than+1, be less than-1 and the part different colours of three curvature range between ± 1 indicate;

The determination of e, regions of non-interest and removing: statistical study is carried out to each point curvature value of the described pulmonary parenchyma region initial surface model calculated, determines the threshold value of a scope between 0.08-0.12 according to analysis result; Namely the part each point curvature value calculated in Step d being exceeded to this threshold value is defined as regions of non-interest; After the described regions of non-interest on the initial surface model of described pulmonary parenchyma region being removed, namely obtain the described pulmonary parenchyma region basic surface model of SOL;

F, after marginalisation repairing is carried out to described pulmonary parenchyma region basic surface model, the final surface model in a smooth surface and complete pulmonary parenchyma region can be obtained, be referred to as: " pulmonary parenchyma region surface model ".

The mode that marginalisation of the present invention is repaired uses the surface fracture region of radial basis function to described pulmonary parenchyma region basic surface model to repair.

The present invention is by drawing the 3-dimensional image in pulmonary parenchyma region, and by the removal of regions of non-interest and the repairing to damaged area, define the surface model of the complete lung's area-of-interest in surface, complete the information extracted lung's area-of-interest and contain from medical volume data, both the workload that doctor observes CT scan image had been decreased, again the image information of complexity and internal relation are shown in mode intuitively, thus assist physician carries out comprehensive and accurate analysis to pathology body and surrounding tissue, improve science and the accuracy of medical diagnosis.

The invention provides a kind of method of area of computer aided diagnostic imaging, the diagnosis difference that the difference due to subjective factors such as doctor's subjective experience, observation abilitys causes can be eliminated, and provide the higher reference identifying and diagnosing result of accuracy rate, thus diagnostic imaging is objectified more, improve efficiency and the accuracy of diagnosis.

Accompanying drawing explanation

Fig. 1 is the process flow diagram of modeling method of the present invention.

Fig. 2 is the triangular vector figure of gauging surface curvature.

Tu3Shi lung original CT scan picture.

Fig. 4 is the image in the pulmonary parenchyma region of extracting.

Fig. 5 is set up pulmonary parenchyma region initial surface illustraton of model.

Fig. 6 is the schematic diagram asking curvature on the initial surface model of pulmonary parenchyma region.

Fig. 7 is the pulmonary parenchyma region basic surface illustraton of model with local damage after removal regions of non-interest.

Fig. 8 is the final surface model figure in the smooth surface after repairing, complete pulmonary parenchyma region.

Embodiment

The software and hardware condition of embodiment of the present invention computing machine used is: dual-Core CPU E5800@3.20GHz, video card is NVIDIAGeForce GT430, internal memory 2.0GB, and operating system is WindowXP, and software-programming languages uses c++.

As shown in Figure 1, the implementing procedure of modeling method of the present invention is:

The first step, the smoothing operation of former figure to CT lung scans image: input the thoracic CT slices image that a group comprises the DICOM form of complete lung organ, one group of totally 564, former figure (quantity can have increase and decrease), size is 512 × 512 pixels (Fig. 3 is wherein one).Utilize Gaussian filter to the smoothing operation of the former figure of this group, to eliminate the image noise interference in each figure.Because medical image all contains a large amount of noises usually, impact can be brought on segmentation result, thus need to carry out denoising to image.

Second step, obtain pulmonary parenchyma region: on the former figure of this group after denoising, obtain initial pulmonary parenchyma region (Fig. 4), recycle three-dimensional FastMarching algorithm, carry out three-dimensional segmentation to the pulmonary parenchyma region obtained, concrete steps are: from the former figure of this group after denoising wherein one to have in the picture of complete pulmonary parenchyma selected Seed Points (X, a Y, Z), wherein X is the horizontal ordinate of Seed Points, and Y is the ordinate of Seed Points, and Z is selected picture is which is opened in this picture group sheet.Because the anatomical position of lung in thoracic cavity is relatively-stationary, so the position of pulmonary parenchyma is also relatively fixing in chest CT image, be positioned at the center of picture, computing machine can carry out cutting operation according to the position of Seed Points, and segmentation result saved as the picture of two-dimentional DICOM form and cover Hui Yuantu, obtain the original gray value (gray-scale value in the present embodiment is-700) in pulmonary parenchyma region.Because the order of section can have an impact to next step surface modeling, so will arrange by former figure sequence number the result be partitioned into.

3rd step, Marching Cubes algorithm is used to carry out surface modeling to pulmonary parenchyma region: the one group of sequence picture handled well by second step reads in calculator memory as 3 d-dem regular data field, set the gray-scale value of its contour surface as-700, computing machine can to each cube (Cube) computing one by one in described data fields, equivalent point identical with original gray value on each limit of each cube is drawn by linear interpolation, then the contour surface in this cube is simulated with a series of tri patch, expression is approached as in this cube one of contour surface.Cube computation on all former figure is also classified complete, namely builds the initial surface model (Fig. 5) in pulmonary parenchyma region.

4th step, determine the surface curvature of each point on the initial surface model of pulmonary parenchyma region, to determine the concrete form feature in pulmonary parenchyma region: the surface curvature calculating three vertex position places in described initial surface model on each tri patch respectively, and the curvature value of pulmonary parenchyma region initial surface model is greater than+1, be less than-1 and the part different colours of three curvature range between ± 1 indicate, to show in the picture better, and effects on surface curvature is facilitated to carry out statistical study.The present invention uses method of finite difference to ask maximum curvature to calculate, and result as shown in Figure 6.

Surface curvature on the initial surface model of pulmonary parenchyma region can symbolize the geometric attribute of body surface intuitively, than being easier to the concrete form feature determining object.Consider the style characteristic in pulmonary parenchyma region, therefore also can determine according to the size of surface curvature the regions of non-interest that needs to remove.In the present invention, the surface curvature of pulmonary parenchyma region initial surface model can calculate with reference to under type.

For three-dimensional body, surface normal curvature k in one direction _nit is the radius of a circle of cutting into slices closest to the normal in that direction.And the change of surface curvature is relevant with direction, but concerning a smooth surface, should meet:

$k_n=\left(\begin{array}{cc}s& t\end{array}\right)\left(\begin{array}{cc}e& f\\ f& g\end{array}\right)\left(\begin{array}{c}s\\ t\end{array}\right)=\left(\begin{array}{cc}s& t\end{array}\right)\mathrm{II}\left(\begin{array}{c}s\\ t\end{array}\right)$

Wherein, the unity-length vector that (s t) is local tangential plane, symmetric matrix II is Weingarten matrix or the second fundamental tensor, also can be obtained by the diagonalization of axis rotation method system:

$k_n=\left(\begin{array}{cc}{s}^{\′}& t^{\′}\end{array}\right)\left(\begin{array}{cc}{k}_1& 0\\ 0& k_2\end{array}\right)\left(\begin{array}{c}{s}^{\′}\\ t^{\′}\end{array}\right)=k_1{s}^{\′2}+k_2{t}^{\′2}$

Wherein, k ₁and k ₂for principal curvatures, s' and t' is principal direction.

The computing method of the surface curvature of surface model of the present invention are finite difference methods:

First, the surface curvature of each is calculated.The second fundamental tensor II is defined as follows:

$\mathrm{II}=\left(\begin{array}{cc}{D}_{n}n& D_{v}n\end{array}\right)=\left(\begin{array}{cc}\frac{{\∂}_n}{{\∂}_u}u& \frac{{\∂}_n}{{\∂}_v}u\\ \frac{{\∂}_n}{{\∂}_u}v& \frac{{\∂}_n}{{\∂}_v}v\end{array}\right)$

Wherein, (u is v) the vector of unit length of orthogonal coordinate system tangential direction, the vector in any direction in section is multiplied by (u is v) the inverse of the normal in this direction:

Ⅱs=D _sn

The vector of the inverse that it should be noted that normal inherently on section, it often has an ingredient on the s of direction, but in plane vertical with it, also have an ingredient (namely causing the surface of " torsion ").

Although this definition is only applicable to smooth surface, in the discrete case, we adopt finite difference.Such as, we define the difference of the normal on each summit on three direction vectors (edge) of a triangle (Fig. 2) and these directions, then have:

$\mathrm{II}\left(\begin{array}{c}{e}_{0}u\\ e_{0}v\end{array}\right)=\left(\begin{array}{c}(n_2-n_1)u\\ (n_2-n_1)v\end{array}\right)$

$\mathrm{II}\left(\begin{array}{c}{e}_{1}u\\ e_{1}v\end{array}\right)=\left(\begin{array}{c}(n_0-n_2)u\\ (n_0-n_2)v\end{array}\right)$

$\mathrm{II}\left(\begin{array}{c}{e}_{2}u\\ e_{2}v\end{array}\right)=\left(\begin{array}{c}(n_1-n_0)u\\ (n_1-n_0)v\end{array}\right)$

This provide the Linear Constraints of one group of the second fundamental tensor element, and this condition can use least square method to determine.It should be noted that this is estimated based on explication, unless this triangle self has the summit of three conllinear.

Secondly, be the conversion of coordinate system.Definition (u _f, v _f) be the curvature tensor in coordinate system face, suppose that each summit p has himself orthogonal coordinate system (u _p, v _p), this coordinate system is defined in the vertical plane of normal, obtains the formula of curvature tensor to apex coordinate thus.

Situation about will consider here is that the normal on plane and summit is equal, so (u _f, v _f) and (u _p, v _p) coplanar.

The Part I of the second fundamental tensor II can at (u _p, v _p) draw in coordinate system: p p

$e_p=\left(\begin{array}{cc}1& 0\end{array}\right)\left(\begin{array}{cc}{e}_p& f_p\\ f_p& g_p\end{array}\right)\left(\begin{array}{c}1\\ 0\end{array}\right)=u_p^T\mathrm{II}{u}_p$

E _pat (u _f, v _f) in can be expressed as:

$e_p=u_p^T\mathrm{II}{u}_p={\left(\begin{array}{cc}{u}_p& u_f\\ u_p& v_f\end{array}\right)}^T\mathrm{II}\left(\begin{array}{cc}{u}_p& u_f\\ u_p& v_f\end{array}\right)$

In like manner, draw:

$f_p=u_p^T\mathrm{II}{v}_p,g_p=v_p\mathrm{II}{v}_p$

When new and old coordinate not conllinear time, we can not simply by new u _pand v _pproject to old (u _f, v _f) in coordinate system.In general such projection is not unit length, this can cause the loss of the curvature when each coordinate transform (particularly mean curvature to be multiplied by press from both sides between normal cosine of an angle square).So it is coplanar by it and another coordinate system that first we rotate one of them coordinate system, doing so avoids cos ²the loss of θ curvature, adds the accuracy of valuation simultaneously.

5th step, the determination of regions of non-interest and removing: statistical study is carried out to each point curvature value of the above-mentioned pulmonary parenchyma region initial surface model calculated, the threshold value of a scope between 0.08-0.12 is determined according to analysis result, the threshold value value of the embodiment of the present invention is 0.12mm, and namely the region higher than this threshold value is defined as regions of non-interest.By analyzing the region that those may cause splitting error, as blood vessel or tracheae etc., because the radius of tracheae and blood vessel is relatively little, curvature is higher, so should remove these high curvature areas on the initial surface model of pulmonary parenchyma region.After removing these regions of non-interest, pulmonary parenchyma region basic surface model (Fig. 7) of SOL can be obtained.

6th step, marginalisation repairing is carried out to pulmonary parenchyma region basic surface model: above-mentioned pulmonary parenchyma region basic surface model is repaired, because after removing the regions of non-interest of higher curvature, formation pulmonary parenchyma region basic surface model having imperfect surface patches can be caused, so need to fill up the region lacked, the anatomical structure of surface model is made to seem smooth complete.Concrete mode uses the surface fracture region of radial basis function to pulmonary parenchyma region basic surface model to repair, and can obtain the final surface model in a smooth surface and complete pulmonary parenchyma region afterwards, namely pulmonary parenchyma surface model (Fig. 8).

Radial basis function is in fact look into/the process of extrapolation in one, constructs continuous print Implicitly function by a series of discrete sampling point heterogeneous, and only needing to carry out resampling to this function during reconstruction interpolation can go out damaged region.The surface patch that matching is left after removing operation is as implicit function G (x), and G (x) crosses point { x ₁, x ₂..., x _n.In order to obtain smooth matching, the energy of implicit function G (x) need be made minimum, can be measured by the quadratic sum of the second order local derviation to G (x), that is:

$E={\∫}_{\mathrm{\Ω}}(G_{\mathrm{xx}}^2+G_{\mathrm{yy}}^2+G_{\mathrm{zz}}^2+2G_{\mathrm{xy}}^2+2G_{\mathrm{xz}}^2+2G_{\mathrm{yz}}^2)\mathrm{d\Ω}$

Wherein, Ω presentation surface region, E is the degree of roughness measuring G (x).

By using radial basis function (radial basis functions) to represent curved surface, the minimization problem of energy can be solved:

$G\left(X\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_i\mathrm{\Φ}\left(\right||X-X_i|\left|\right)+P\left(X\right)$

Wherein, X=(x, y, z), Φ is radial basis function, λ _ireal coefficient, || x|| is Euclidean distance, { X _idiscrete point, P (X) is linear polynomial: P (X)=cX=c ₀+ c ₁x+c ₂y+c ₃z.

If G (X) is by all discrete point { X _i, then G (X _i)=0.In order to ensure G (X) non-zero, the outer surface pointed by the normal vector of discrete point some points can be increased.

Suppose always to have k outside surface sampled point, then formula G (X) can be written as AB=C, that is:

$\left[\begin{array}{cccccccc}{\mathrm{\Φ}}_{11}& {\mathrm{\Φ}}_{12}& ...& {\mathrm{\Φ}}_{1m}& 1& x_1& y_1& z_1\\ {\mathrm{\Φ}}_{21}& {\mathrm{\Φ}}_{22}& ...& {\mathrm{\Φ}}_{2m}& 1& x_2& y_2& z_2\\ .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .\\ {\mathrm{\Φ}}_{m1}& {\mathrm{\Φ}}_{m2}& ...& {\mathrm{\Φ}}_{\mathrm{mm}}& 1& x_m& y_m& z_m\\ 1& 1& ...& 1& 0& 0& 0& 0\\ x_1& x_2& ...& y_m& 0& 0& 0& 0\\ y_1& y_2& ...& y_m& 0& 0& 0& 0\\ z_1& z_2& ...& z_m& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\λ}}_1\\ {\mathrm{\λ}}_2\\ .\\ .\\ .\\ {\mathrm{\λ}}_m\\ c_0\\ c_1\\ c_2\\ c_3\end{array}\right]=\left[\begin{array}{c}{v}_1\\ v_2\\ .\\ .\\ .\\ v_m\\ 0\\ 0\\ 0\\ 0\end{array}\right]$

Wherein, Φ _ij=Φ (x _i-x _j), i, j ∈ [1, m], m=n+k.When sample is on the surface time, v _i=0; And work as v _iwhen=± 1, then represent this point respectively and be positioned at outside surface or surface.

When being mapped on medical image by implicit surfaces, we can check the implicit function value of each voxel, are worth for null representation border, negative value/and on the occasion of representing region, inside/outside portion respectively.

Utilize modeling method of the present invention, also can carry out the foundation of independent surface model to other area-of-interests such as Patients with Lung inner disease focis.

Claims (2)

1., based on a method for building up for the pulmonary parenchyma region surface model of CT image, it is characterized in that, comprise the following steps:

A, the one group of smoothing operation of every pictures of sorting in former figure utilizing computing machine to scan lung CT, to eliminate image noise interference;

B, acquisition pulmonary parenchyma region: on the former figure of this group after denoising, obtain initial pulmonary parenchyma region, recycle three-dimensional Fast Marching algorithm, three-dimensional segmentation is carried out to the pulmonary parenchyma region obtained, specifically: from the former figure of this group after denoising wherein one to have in the picture of complete pulmonary parenchyma a selected Seed Points (X, Y, Z), wherein X is the horizontal ordinate of Seed Points, Y is the ordinate of Seed Points, and Z is selected picture is which is opened in this picture group sheet; Computing machine carries out cutting operation according to the position of Seed Points, and segmentation result is covered Hui Yuantu, obtains the original gray value in described pulmonary parenchyma region, and arranges by former figure sequence number;

C, surface modeling is carried out to pulmonary parenchyma region: utilize Marching Cubes algorithm, sequence picture good for b step process is read in calculator memory as 3 d-dem regular data field, to each cube computing one by one in described data fields, draw equivalent point identical with described original gray value on each limit of each cube by linear interpolation, then simulate the contour surface in this cube with tri patch; All cubical contour surface in whole described data fields forms the initial surface model in described pulmonary parenchyma region;

D, determine the surface curvature of pulmonary parenchyma region initial surface model each point: the surface curvature calculating three vertex position places in the initial surface model of described pulmonary parenchyma region on each tri patch respectively, and the curvature value of described pulmonary parenchyma region initial surface model is greater than+1, be less than-1 and the part different colours of three curvature range between ± 1 indicate;

The determination of e, regions of non-interest and removing: statistical study is carried out to each point curvature value of the described pulmonary parenchyma region initial surface model calculated, determines the threshold value of a scope between 0.08-0.12 according to analysis result; Namely the part each point curvature value calculated in d step being exceeded to this threshold value is defined as regions of non-interest; After the described regions of non-interest on the initial surface model of described pulmonary parenchyma region being removed, namely obtain the described pulmonary parenchyma region basic surface model of SOL;

F, after marginalisation repairing is carried out to described pulmonary parenchyma region basic surface model, the final surface model in a smooth surface and complete pulmonary parenchyma region can be obtained.

2. according to the method for building up of the pulmonary parenchyma region surface model based on CT image described in claim 1, it is characterized in that, the mode that described marginalisation is repaired uses the surface fracture region of radial basis function to described pulmonary parenchyma region basic surface model to repair.