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

Abstract

The invention relates to a method for establishing a lung parenchyma surface model based on CT images. The method is to use a computer to perform Gaussian smoothing operations on each original image of a lung CT scan image to eliminate image noise interference; and then use three-dimensional Fast Marching The algorithm obtains the lung parenchyma area; uses the Marching Cubes algorithm; then uses the calculation of the surface curvature to determine the specific morphological characteristics of the lung parenchyma, and removes the determined non-interesting area, and finally obtains a smooth and complete surface through marginalization repair A surface model of the lung region of interest. The invention provides a method for computer-aided image diagnosis, which can eliminate the difference in diagnosis caused by subjective factors such as doctor's subjective experience and observation ability, and provide a reference identification and diagnosis result with high accuracy, so that image diagnosis It is more objective and improves the efficiency and accuracy of diagnosis.

Description

Method of establishing surface model of lung parenchyma based on CT images

technical field

The invention relates to a medical image processing method, in particular to a method for establishing a surface model of a lung parenchyma region based on a CT image.

Background technique

Lung cancer is a common primary lung tumor, and its morbidity and mortality are increasing year by year. It has become the malignant tumor with the highest morbidity and mortality worldwide. According to the survey results of the Ministry of Health in 2006, the mortality rate of lung cancer in my country has increased significantly, and it has become the number one cause of cancer death, which has increased by 75.77% compared with the 1990s. Excluding the influence of age structure changes, it has also increased by 33.25%. The clinical manifestations of lung cancer are various, and most patients have no obvious symptoms and signs in the early stage, so it is difficult to find out, and they often have reached the late stage when they are discovered. If lung cancer can be diagnosed and treated at an early stage, the cure rate and survival rate of patients can be improved accordingly.

CT (Computer Tomography) is the most commonly used technique in chest imaging and is widely used in the detection and diagnosis of chest and lung diseases. Lung cancer is usually manifested in the form of pulmonary nodules in CT images, and pulmonary nodules have various shapes, sizes, and distribution positions are not fixed, and are easily connected with other tissues; the density of lung cancer is similar to that of pulmonary blood vessels. Most of them appear as circular or nearly circular dense spots in CT images, which are difficult to distinguish by human eyes alone. Moreover, among nearly a hundred CT images generated by a lung CT scan, there may be only a few images with nodules. The work of analyzing all image data is very boring and tedious, and even experienced doctors will inevitably miss or misdiagnose. Therefore, computer-aided diagnosis (Computer-Aided Diagnosis) has become an important method to solve this problem. As a part of computer-aided diagnosis, the rapid and accurate segmentation of the Region of Interest is the first step in the work, and it is also an important basis for the next step of detection and diagnostic classification.

When doctors diagnose patients' diseases through chest CT scan images, they are usually interested in areas such as lung parenchyma and intrapulmonary lesions, which are called "regions of interest" on the CT scan images. Accurately segmenting regions of interest in CT scan images, identifying and classifying lung tissue, and diagnosing lung cancer are all crucial. However, there is currently no method for establishing a surface model of a lung region of interest in a CT scan image using a computer.

Contents of the invention

The object of the present invention is to provide a method for establishing a surface model of the lung parenchyma region based on CT images to solve the problem that computer modeling cannot be performed on the lung parenchyma region in the CT scan image, so as to meet the needs of clinical research, diagnosis and treatment. need.

The object of the present invention is achieved in that a method for establishing a surface model of a lung parenchyma region based on a CT image comprises the following steps:

a. Use a computer to perform a smoothing operation on each picture in a group of sorted original pictures of lung CT scans to eliminate image noise interference;

b. Obtain the lung parenchyma area: obtain the initial lung parenchyma area on the original image after denoising processing, and then use the 3D FastMarching algorithm to perform 3D segmentation on the obtained lung parenchyma area, and overlay the segmentation result back to the original image , to obtain the original gray value of the lung parenchyma region, and arrange them according to the serial number of the original image;

c. Perform surface modeling on the lung parenchyma area: use the Marching Cubes algorithm to read the sorted images processed in step b into the computer memory as a three-dimensional discrete regular data field, and calculate each cube (Cube) in the data field one by one processing, by linear interpolation to obtain the isovalue points on each side of each cube that are the same as the original gray value, and then use the triangular patch to fit the isosurface in the cube; all the data fields in the entire the isosurfaces of the cube constitute the initial surface model of the lung parenchymal region;

d. Determine the surface curvature of each point of the initial surface model of the lung parenchyma region: respectively calculate the surface curvatures at the three vertices on each triangle patch in the initial surface model of the lung parenchyma region, and convert the initial surface of the lung parenchyma region to The parts of the model whose curvature values are greater than +1, less than -1 and within ±1 are marked with different colors;

e. Determination and removal of non-interest regions: perform statistical analysis on the calculated curvature values of each point of the initial surface model of the lung parenchyma region, and determine a threshold value in the range of 0.08-0.12 according to the analysis results; for d The portion where the curvature value of each point calculated in the step is higher than the threshold is determined as a non-interest area; after removing the non-interest area on the initial surface model of the lung parenchyma region, the area with local defects is obtained. a basic surface model of the lung parenchyma region;

f. After marginalizing and repairing the basic surface model of the lung parenchyma region, a smooth and complete final surface model of the lung parenchyma region can be obtained, which is called "the surface model of the lung parenchyma region".

The method of repairing the marginalization in the present invention is to use radial basis function to repair the surface damage area of the basic surface model of the lung parenchyma region.

The present invention forms a surface model of the lung region of interest with a complete surface through the three-dimensional image rendering of the lung parenchyma region, and removes the non-interest region and repairs the damaged region, and completes the extraction from the medical volume data. The information contained in the area of interest in the lungs not only reduces the workload of doctors observing CT scan images, but also displays complex image information and internal relationships in an intuitive way, thereby assisting doctors to comprehensively and accurately detect diseased objects and surrounding tissues. The analysis can improve the scientificity and accuracy of medical diagnosis.

The invention provides a method for computer-aided image diagnosis, which can eliminate the difference in diagnosis caused by subjective factors such as doctors' subjective experience and observation ability, and provide a reference identification and diagnosis result with high accuracy, so that image diagnosis It is more objective and improves the efficiency and accuracy of diagnosis.

Description of drawings

Fig. 1 is a flowchart of the modeling method of the present invention.

Figure 2 is a triangular vector diagram for calculating surface curvature.

Figure 3 is the original CT scan image of the lungs.

Figure 4 is an image of the extracted lung parenchyma region.

Fig. 5 is a diagram of the established initial surface model of the lung parenchyma region.

Fig. 6 is a schematic diagram of calculating the curvature on the initial surface model of the lung parenchyma region.

Fig. 7 is a diagram of the basic surface model of the lung parenchyma region with local damage after removal of the non-interest region.

Fig. 8 is a final surface model diagram of the repaired smooth surface and complete lung parenchyma region.

Detailed ways

The software and hardware condition of computer used in the embodiment of the present invention is: Dual-Core CPU E5800 3.20GHz, graphics card is NVIDIAGeForce GT430, memory 2.0GB, operating system is WindowXP, software programming language is c++.

As shown in Figure 1, the implementation process of the modeling method of the present invention is:

The first step is to perform a smoothing operation on the original image of the lung CT scan image: input a set of chest CT slice images in DICOM format containing complete lung organs, a set of original images with a total of 564 pieces (the number can be increased or decreased), the size It is 512×512 pixels (Figure 3 is one of them). The Gaussian filter is used to smooth the group of original images to eliminate the image noise interference in each image. Because medical images usually contain a lot of noise, which will affect the segmentation results, it is necessary to denoise the image.

The second step is to obtain the lung parenchyma area: obtain the initial lung parenchyma area on the original image after denoising processing (Figure 4), and then use the 3D FastMarching algorithm to perform 3D segmentation on the obtained lung parenchyma area. The specific steps is: select a seed point (X, Y, Z) from one of the original images with complete lung parenchyma after denoising, where X is the abscissa of the seed point, and Y is the The ordinate, Z is the number of the selected picture in this group of pictures. Since the anatomical position of the lung in the thoracic cavity is relatively fixed, the position of the lung parenchyma in the chest CT image is also relatively fixed. It is located in the center of the picture. The computer will perform the segmentation operation according to the position of the seed point, and save the segmentation result as The image in the two-dimensional DICOM format was overlaid back to the original image to obtain the original gray value of the lung parenchyma region (the gray value in this embodiment is -700). Since the order of slices will affect the next step of surface modeling, the segmented results should be arranged according to the serial number of the original image.

The third step is to use the Marching Cubes algorithm to model the surface of the lung parenchyma: read a group of sorted pictures processed in the second step into the computer memory as a three-dimensional discrete regular data field, and set the gray value of the isosurface to be -700, the computer will calculate and process each cube (Cube) in the data field one by one, obtain the equivalent point on each side of each cube with the same value as the original gray value through linear interpolation, and then use a series of triangles The isosurface in the cube is fitted by the patch as an approximate representation of the isosurface in the cube. All the cubes on the original image were calculated and classified, and the initial surface model of the lung parenchyma region was constructed (Fig. 5).

The fourth step is to determine the surface curvature of each point on the initial surface model of the lung parenchyma region to determine the specific morphological characteristics of the lung parenchyma region: respectively calculate the surface curvature at the three vertices on each triangular patch in the initial surface model , and the parts of the three curvature ranges whose curvature values are greater than +1, less than -1 and between ±1 of the initial surface model of the lung parenchyma region are marked with different colors to better display in the image and facilitate the identification of Statistical analysis of surface curvature. The present invention uses the finite difference method to calculate the maximum curvature, and the result is shown in FIG. 6 .

The surface curvature on the initial surface model of the lung parenchyma region can intuitively represent the geometric properties of the object surface, and it is relatively easy to determine the specific morphological characteristics of the object. Considering the shape characteristics of the lung parenchyma region, the non-interest region to be removed can also be determined according to the magnitude of the surface curvature. In the present invention, the surface curvature of the initial surface model of the lung parenchyma region can be calculated with reference to the following manner.

For a three-dimensional object, the surface normal curvature kn in a certain direction _is the radius of the circle closest to the normal slice in that direction. The change of surface curvature is related to the direction, but for a smooth surface, it should satisfy:

${\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; =</span>\left(\begin{array}{cc}\mathrm{<span\; class="notranslate">\; the\; s</span>}& \mathrm{<span\; class="notranslate">\; t</span>}\end{array}\right)\left(\begin{array}{cc}\mathrm{<span\; class="notranslate">\; e</span>}& \mathrm{<span\; class="notranslate">\; f</span>}\\ \mathrm{<span\; class="notranslate">\; f</span>}& \mathrm{<span\; class="notranslate">\; g</span>}\end{array}\right)\left(\begin{array}{c}\mathrm{<span\; class="notranslate">\; the\; s</span>}\\ \mathrm{<span\; class="notranslate">\; t</span>}\end{array}\right)<span\; class="notranslate">\; =</span>\left(\begin{array}{cc}\mathrm{<span\; class="notranslate">\; the\; s</span>}& \mathrm{<span\; class="notranslate">\; t</span>}\end{array}\right)\mathrm{<span\; class="notranslate">\; II</span>}\left(\begin{array}{c}\mathrm{<span\; class="notranslate">\; the\; s</span>}\\ \mathrm{<span\; class="notranslate">\; t</span>}\end{array}\right)$

Among them, (s t) is the unit length vector of the local tangent plane, and the symmetric matrix II is the Weingarten matrix or the second basic tensor, which can also be obtained by diagonalizing the rotating local coordinate system:

${\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; =</span>\left(\begin{array}{cc}{\mathrm{<span\; class="notranslate">\; the\; s</span>}}^{<span\; class="notranslate">\; \&prime;</span>}& {\mathrm{<span\; class="notranslate">\; t</span>}}^{<span\; class="notranslate">\; \&prime;</span>}\end{array}\right)\left(\begin{array}{cc}{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}& \mathrm{<span\; class="notranslate">\; 0</span>}\\ \mathrm{<span\; class="notranslate">\; 0</span>}& {\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}\end{array}\right)\left(\begin{array}{c}{\mathrm{<span\; class="notranslate">\; the\; s</span>}}^{<span\; class="notranslate">\; \&prime;</span>}\\ {\mathrm{<span\; class="notranslate">\; t</span>}}^{<span\; class="notranslate">\; \&prime;</span>}\end{array}\right)<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; the\; s</span>}}^{<span\; class="notranslate">\; \&prime;</span>\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; t</span>}}^{<span\; class="notranslate">\; \&prime;</span>\mathrm{<span\; class="notranslate">\; 2</span>}}$

Among them, k ₁ and k ₂ are the main curvatures, and s' and t' are the main directions.

The calculation method of the surface curvature of the surface model of the present invention is a finite difference method:

First, calculate the surface curvature for each face. The second basic tensor II is defined as follows:

$\mathrm{<span\; class="notranslate">\; II</span>}<span\; class="notranslate">\; =</span>\left(\begin{array}{cc}{\mathrm{<span\; class="notranslate">\; D.</span>}}_{\mathrm{<span\; class="notranslate">\; no</span>}}\mathrm{<span\; class="notranslate">\; no</span>}& {\mathrm{<span\; class="notranslate">\; D.</span>}}_{\mathrm{<span\; class="notranslate">\; v</span>}}\mathrm{<span\; class="notranslate">\; no</span>}\end{array}\right)<span\; class="notranslate">\; =</span>\left(\begin{array}{cc}\frac{{<span\; class="notranslate">\; \&PartialD;</span>}_{\mathrm{<span\; class="notranslate">\; no</span>}}}{{<span\; class="notranslate">\; \&PartialD;</span>}_{\mathrm{<span\; class="notranslate">\; u</span>}}}\mathrm{<span\; class="notranslate">\; u</span>}& \frac{{<span\; class="notranslate">\; \&PartialD;</span>}_{\mathrm{<span\; class="notranslate">\; no</span>}}}{{<span\; class="notranslate">\; \&PartialD;</span>}_{\mathrm{<span\; class="notranslate">\; v</span>}}}\mathrm{<span\; class="notranslate">\; u</span>}\\ \frac{{<span\; class="notranslate">\; \&PartialD;</span>}_{\mathrm{<span\; class="notranslate">\; no</span>}}}{{<span\; class="notranslate">\; \&PartialD;</span>}_{\mathrm{<span\; class="notranslate">\; u</span>}}}\mathrm{<span\; class="notranslate">\; v</span>}& \frac{{<span\; class="notranslate">\; \&PartialD;</span>}_{\mathrm{<span\; class="notranslate">\; no</span>}}}{{<span\; class="notranslate">\; \&PartialD;</span>}_{\mathrm{<span\; class="notranslate">\; v</span>}}}\mathrm{<span\; class="notranslate">\; v</span>}\end{array}\right)$

Among them, (u, v) is the unit vector in the tangent direction of the orthogonal coordinate system, and the vector in any direction of the tangent plane multiplied by (u, v) is the reciprocal of the normal in this direction:

Ⅱs=D _s n

Note that the inverse of the normal is itself a vector on the tangent plane, which tends to have a component in the direction s, but also a component in the plane perpendicular to it (i.e. the surface causing the "twist") .

Although this definition only applies to smooth surfaces, in the discrete case we employ finite differences. For example, we define the three vector directions (edges) of a triangle (Fig. 2) and the difference of the normals of each vertex in these directions, then:

$\mathrm{<span\; class="notranslate">\; II</span>}\left(\begin{array}{c}{\mathrm{<span\; class="notranslate">\; e</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}\mathrm{<span\; class="notranslate">\; u</span>}\\ {\mathrm{<span\; class="notranslate">\; e</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}\mathrm{<span\; class="notranslate">\; v</span>}\end{array}\right)<span\; class="notranslate">\; =</span>\left(\begin{array}{c}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; no</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; no</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span>\mathrm{<span\; class="notranslate">\; u</span>}\\ <span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; no</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; no</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span>\mathrm{<span\; class="notranslate">\; v</span>}\end{array}\right)$

$\mathrm{<span\; class="notranslate">\; II</span>}\left(\begin{array}{c}{\mathrm{<span\; class="notranslate">\; e</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}\mathrm{<span\; class="notranslate">\; u</span>}\\ {\mathrm{<span\; class="notranslate">\; e</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}\mathrm{<span\; class="notranslate">\; v</span>}\end{array}\right)<span\; class="notranslate">\; =</span>\left(\begin{array}{c}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; no</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; no</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span>\mathrm{<span\; class="notranslate">\; u</span>}\\ <span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; no</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; no</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span>\mathrm{<span\; class="notranslate">\; v</span>}\end{array}\right)$

$\mathrm{<span\; class="notranslate">\; II</span>}\left(\begin{array}{c}{\mathrm{<span\; class="notranslate">\; e</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; u</span>}\\ {\mathrm{<span\; class="notranslate">\; e</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; v</span>}\end{array}\right)<span\; class="notranslate">\; =</span>\left(\begin{array}{c}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; no</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; no</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span>\mathrm{<span\; class="notranslate">\; u</span>}\\ <span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; no</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; no</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span>\mathrm{<span\; class="notranslate">\; v</span>}\end{array}\right)$

This provides a set of linear constraints on the elements of the second base tensor, which can be determined using least squares. Note that this estimate is based on exact definitions, unless the triangle itself has three collinear vertices.

Secondly, it is the conversion of the coordinate system. Define (u _f , v _f ) as the curvature tensor of the coordinate system surface, assuming that each vertex p has its own orthogonal coordinate system (u _p , v _p ), which is defined on the vertical plane of the normal, From this one obtains the formula for curvature tensor to vertex coordinates.

The situation to be considered here is that the normals of the plane and the vertex are equal, so (u _f ,v _f ) and (u _p ,v _p ) are coplanar.

The first part of the second elementary tensor II can be drawn in the (u _p , v _p ) coordinate system: p p

${\mathrm{<span\; class="notranslate">\; e</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}<span\; class="notranslate">\; =</span>\left(\begin{array}{cc}\mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}\end{array}\right)\left(\begin{array}{cc}{\mathrm{<span\; class="notranslate">\; e</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}& {\mathrm{<span\; class="notranslate">\; f</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}\\ {\mathrm{<span\; class="notranslate">\; f</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}& {\mathrm{<span\; class="notranslate">\; g</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}\end{array}\right)\left(\begin{array}{c}\mathrm{<span\; class="notranslate">\; 1</span>}\\ \mathrm{<span\; class="notranslate">\; 0</span>}\end{array}\right)<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}\mathrm{<span\; class="notranslate">\; II</span>}{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}$

e _p in (u _f , v _f ) can be expressed as:

${\mathrm{<span\; class="notranslate">\; e</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}\mathrm{<span\; class="notranslate">\; II</span>}{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}<span\; class="notranslate">\; =</span>{\left(\begin{array}{cc}{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}& {\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; f</span>}}\\ {\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}& {\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; f</span>}}\end{array}\right)}^{\mathrm{<span\; class="notranslate">\; T</span>}}\mathrm{<span\; class="notranslate">\; II</span>}\left(\begin{array}{cc}{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}& {\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; f</span>}}\\ {\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}& {\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; f</span>}}\end{array}\right)$

Similarly, it is concluded that:

${\mathrm{<span\; class="notranslate">\; f</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}\mathrm{<span\; class="notranslate">\; II</span>}{\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; g</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}\mathrm{<span\; class="notranslate">\; II</span>}{\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}}$

We cannot simply project the new u _p and v _p into the old (u _f ,v _f ) coordinate system when the old and new coordinates are not collinear. In general such projections are not unit length, which results in a loss of curvature for each coordinate transformation (in particular the mean curvature is multiplied by the square of the cosine of the angle between the normals). Therefore, we first rotate one of the coordinate systems so that it is coplanar with the other coordinate system, which avoids the loss of cos ² θ curvature and increases the accuracy of the estimation.

The fifth step is to determine and remove the non-interest area: perform statistical analysis on the curvature values of each point of the initial surface model of the lung parenchyma area calculated above, and determine a threshold value ranging from 0.08 to 0.12 according to the analysis results. In the embodiment of the invention, the threshold value is 0.12mm, and the region higher than this threshold value is determined as the non-interest region. By analyzing those areas that may cause segmentation errors, such as blood vessels or trachea, since the trachea and blood vessels have relatively small radii and high curvature, these high-curvature areas should be removed from the initial surface model of the lung parenchyma. After removing these non-interest regions, the basic surface model of the lung parenchyma with local defects can be obtained (Fig. 7).

The sixth step is to repair the edge of the basic surface model of the lung parenchyma region: the above basic surface model of the lung parenchyma region is repaired because the removal of the non-interest region with high curvature will cause the basic surface model of the lung parenchyma region to There is formation of incomplete surface fragments, so missing areas need to be filled to make the anatomy of the surface model appear smooth and complete. The specific method is to use the radial basis function to repair the surface damaged area of the basic surface model of the lung parenchyma region, and then obtain a smooth and complete final surface model of the lung parenchyma region, that is, the lung parenchyma surface model (Fig. 8) .

The radial basis function is essentially a process of interpolation/extrapolation. A continuous implicit function is constructed through a series of non-uniform discrete sampling points. During reconstruction, only the function needs to be re-sampled to interpolate the damaged area. . Fit the remaining surface patch after the removal operation as an implicit function G(x), and G(x) passes through the points {x ₁ ,x ₂ ,…,x _n }. In order to obtain a smooth fit, the energy of the implicit function G(x) needs to be minimized, which can be measured by the sum of the squares of the second-order partial derivatives of G(x), namely:

$\mathrm{<span\; class="notranslate">\; E.</span>}<span\; class="notranslate">\; =</span>{<span\; class="notranslate">\; \&Integral;</span>}_{\mathrm{<span\; class="notranslate">\; \&Omega;</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; G</span>}}_{\mathrm{<span\; class="notranslate">\; xx</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; G</span>}}_{\mathrm{<span\; class="notranslate">\; yy</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; G</span>}}_{\mathrm{<span\; class="notranslate">\; zz</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 2</span>}{\mathrm{<span\; class="notranslate">\; G</span>}}_{\mathrm{<span\; class="notranslate">\; xy</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 2</span>}{\mathrm{<span\; class="notranslate">\; G</span>}}_{\mathrm{<span\; class="notranslate">\; xz</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 2</span>}{\mathrm{<span\; class="notranslate">\; G</span>}}_{\mathrm{<span\; class="notranslate">\; yz</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span>\mathrm{<span\; class="notranslate">\; d\&Omega;</span>}$

where Ω represents the surface area and E is the roughness of the measurement G(x).

The energy minimization problem can be solved by representing surfaces using radial basis functions:

$\mathrm{<span\; class="notranslate">\; G</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; no</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}\mathrm{<span\; class="notranslate">\; \&Phi;</span>}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; P</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span>$

Among them, X=(x, y, z), Φ is the radial basis function, λ _i is the real coefficient, ||x|| is the Euclidean distance, {X _i } is a discrete point, P(X) is a linear polynomial : P(X)=cX=c ₀ +c ₁ x+c ₂ y+c ₃ z.

If G(X) passes through all discrete points {X _i }, then G(X _i )=0. To ensure that G(X) is non-zero, points can be added at the outer surface to which the normal vectors of discrete points point.

Assuming that there are a total of k sampling points on the outer surface, the formula G(X) can be written as A·B=C, namely:

$\left[\begin{array}{cccccccc}{\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; 11</span>}}& {\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; 12</span>}}& <span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span>& {\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}\mathrm{<span\; class="notranslate">\; m</span>}}& \mathrm{<span\; class="notranslate">\; 1</span>}& {\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}& {\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}& {\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}\\ {\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; twenty\; one</span>}}& {\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; twenty\; two</span>}}& <span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span>& {\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}\mathrm{<span\; class="notranslate">\; m</span>}}& \mathrm{<span\; class="notranslate">\; 1</span>}& {\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}& {\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}& {\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}\\ <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>\\ <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>\\ <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>& <span\; class="notranslate">\; .</span>\\ {\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}& {\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}\mathrm{<span\; class="notranslate">\; 2</span>}}& <span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span>& {\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; mm</span>}}& \mathrm{<span\; class="notranslate">\; 1</span>}& {\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}& {\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}& {\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}\\ \mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; 1</span>}& <span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span>& \mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}\\ {\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}& {\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}& <span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span>& {\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}\\ {\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}& {\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}& <span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span>& {\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}\\ {\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}& {\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}& <span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span>& {\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}\end{array}\right]\left[\begin{array}{c}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}\\ {\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}\\ <span\; class="notranslate">\; .</span>\\ <span\; class="notranslate">\; .</span>\\ <span\; class="notranslate">\; .</span>\\ {\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}\\ {\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}\\ {\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}\\ {\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}\\ {\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 3</span>}}\end{array}\right]<span\; class="notranslate">\; =</span>\left[\begin{array}{c}{\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}\\ {\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}\\ <span\; class="notranslate">\; .</span>\\ <span\; class="notranslate">\; .</span>\\ <span\; class="notranslate">\; .</span>\\ {\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}\\ \mathrm{<span\; class="notranslate">\; 0</span>}\\ \mathrm{<span\; class="notranslate">\; 0</span>}\\ \mathrm{<span\; class="notranslate">\; 0</span>}\\ \mathrm{<span\; class="notranslate">\; 0</span>}\end{array}\right]$

Among them, Φ _ij =Φ( _xi -x _j ), i,j∈[1,m], m=n+k. When the sampling point is on the surface, v _i =0; and when v _i =±1, it means that the point is located outside or inside the surface, respectively.

When mapping an implicit surface onto a medical image, we can check the value of the implicit function at each voxel, where a value of zero indicates the boundary, and a negative/positive value indicates the inner/outer region, respectively.

Using the modeling method of the present invention, a separate surface model can also be established for other regions of interest such as lesions in the patient's lung.

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.