CN106373168A - Medical image based segmentation and 3D reconstruction method and 3D printing system - Google Patents

Abstract

The invention discloses a medical image based segmentation and 3D reconstruction method and a 3D printing system. The method comprises segmenting a medical image and carrying out 3D reconstruction. The 3D printing system carries out 3D printing on a reconstruction model by utilizing a fused decomposition modeling (FDM) method. In the image segmentation part, a unique segmentation algorithm based on medical image (DICOM) preprocessing is provided, and a segmentation task is completed effectively. In the reconstruction part, an MC (Marching Cubes) algorithm is optimized, and the model in the invention is improved greatly. In an interactive interface, a visualization tool enables interaction between a user and the model to certain extent. Finally, a 3D printer prints the obtained model and obtains a simulated entity of bones, and doctors can touch and make interaction with the simulated entity.

Description

A medical image-based segmentation and three-dimensional reconstruction method, 3D printing system

technical field

The invention relates to the field of medical image processing and its three-dimensional reconstruction technology and 3D printing, in particular to a medical image-based segmentation and three-dimensional reconstruction method and a 3D printing system.

Background technique

At present, research institutions in many countries in the world are devoting themselves to research and development in this field, and have developed some clinical-oriented, simple-function 3D visualization systems for medical images. Its main features are as follows: First, it mainly runs on mainframes and workstations (such as SGI, SUN and other graphics workstations), or as supporting software for some special equipment (such as CT, MRI, etc.), which is expensive and not suitable for my country's national conditions , and most of them only generate two-dimensional images, and rarely use three-dimensional software; instead, they are basically foreign commercial software, and the three-dimensional volume rendering technology for lesions and human tissues has not been fully mastered in China. It should also be noted that although these systems provide only limited clinical functionality, they have already played an important role in clinical research. It is precisely because medical image processing and its three-dimensional reconstruction technology play a huge role in promoting clinical medicine, so it has been valued by many countries in the world, making medical image processing and its visualization technology one of the hot spots in biomedical engineering research, also known as Visualization in scientific computing is one of the most interesting and fastest growing fields. Although researchers at home and abroad have achieved fruitful results in this field, there are still many places worthy of further exploration and research because of its complex and mysterious research objects and the broad fields and subjects involved in the research. In addition, the traditional 3D model can only be projected on the computer screen, which limits the intuitive experience of the human body's 3D vision and causes some information loss. A new technology that has emerged recently: 3D printing, can overcome these limitations and produce a touchable three-dimensional model, thereby more accurately simulating the patient's physical condition and lesion degree. Not only is it beneficial for doctors to diagnose more accurately, but it can also perform simulation operations before clinical operations, so as to reduce as many operational errors as possible in real operations. In the future, 3D printing based on image processing and 3D reconstruction will become a new and important hot direction in the medical field.

Image segmentation is an important part of the field of computer vision (Computer Vision), and its goal is to separate the objects in the image from the background, so as to better lay the foundation for subsequent detection and recognition. In the past few decades, scholars in this field have proposed many segmentation algorithms: among them, the snake algorithm proposed by Kass continuously approaches the segmentation target through dynamic contour lines; the level set algorithm proposed by Sethian and Osher evolves the low-dimensional closed curve The problem is transformed into the implicit equation of the overall evolution of the level set function in the high-dimensional space to solve. Many people integrate other information in the basic level set method. For example, Mr. Wang from the Institute of Automation, Chinese Academy of The segmentation achieved better results; Graph Cuts proposed by Boykov in MICCAI and ICCV in 2000 and 2001 respectively is also a very classic segmentation algorithm, which regards the pixels in the image as the nodes of the graph, and the pixels The relationship between represents the weight of the edge in the graph, and the boundary of the object is obtained by the minimum cut solution. After that, Boykov and his students further studied the energy optimization algorithm, color image segmentation algorithm (Grab Cuts) and so on on this basis.

And 3D reconstruction is an important research direction of Computer Graphics. Its method is often used in 3D design, CAD modeling, and is usually closely related to some content in computer vision. Reconstruction algorithms are mainly divided into surface reconstruction and volume reconstruction. Since this article only focuses on surface reconstruction, the volume reconstruction algorithm will not be introduced in detail. Surface reconstruction is to obtain the surface structure of the three-dimensional model through algorithms, the most commonly used of which is the MC (Marching Cube) algorithm. This is an efficient surface reconstruction algorithm. The triangular patches are obtained through the classification of volume data, and finally the adjacent triangular patches are connected to obtain the reconstruction result. On top of the classic MC algorithm, Ju Tao proposed an interactive topology modification to correct some errors in the model; Jiantao Pu further matched the results obtained by MC through the puzzle game, and successfully obtained the reconstructed structure of the lung trachea. In addition, there is also the Poisson reconstruction proposed by Kazhdan, which can reconstruct the directed point set into a three-dimensional model of watertight structure through Sutokes formula.

Contents of the invention

The technical problem to be solved by the present invention is how to perform optimized image segmentation and three-dimensional reconstruction of medical images, and at the same time, how to 3D print the model to improve the degree of wide application in the medical field.

To solve the above technical problems, the present invention provides a medical image-based segmentation and three-dimensional reconstruction method, comprising the following steps: processing redundant information in the original medical image, removing information within the upper and lower threshold ranges of pixel values; Increase the contrast of the image, extract the boundary information between the bones and the background in the above image, and obtain a gradient map; use the canny algorithm to detect the edge of the image after Gaussian smoothing operation on the gradient map, determine the contour and add contour restrictions; segment the above image The obtained images are reconstructed in 3D by MC algorithm.

Further, the method also includes the following operation steps of visualization and graphical interface:

VTK is used to realize the visualization of the 3D model after the above 3D reconstruction, and Qt is used as GUI;

And, displaying a two-dimensional image and processing interaction with a user through the first graphical interface;

The generated three-dimensional model is displayed and the interaction with the user is processed through the second graphical interface.

Furthermore, when performing three-dimensional reconstruction by the MC algorithm, the following optimization steps are included:

Use triangular strips to connect scattered triangular patches to reduce memory usage;

According to the coordinates of three vertices that need to be stored in a single triangular patch, N triangles need to store 3N points,

If these triangular faces are connected into a triangular strip, N+2 points are needed.

Furthermore, when performing three-dimensional reconstruction by the MC algorithm, the following optimization steps are included:

Reduce the time and space complexity by reducing the number of triangle faces;

First, determine which type each triangle belongs to, and then decide which triangle to delete according to the attributes of each type of triangle;

Second, relocate the resulting holes and fill them in until the user-set reduction ratio is reached.

Furthermore, when performing three-dimensional reconstruction by the MC algorithm, the following optimization steps are included:

Use connectivity detection to filter discrete noise regions and remove impurities;

First, arrange all the triangular faces in sequence, and then use the first triangular face as the root of the tree;

Second, gradually find all connected triangles as a whole;

Finally, the entire model is divided into several disjoint regions, and the number of triangles in each region is calculated, and the two largest regions are selected.

Furthermore, when performing three-dimensional reconstruction by the MC algorithm, the following optimization steps are included:

Detect holes and refine the model by filling holes;

First, traverse all the edges. If an edge is not used by two or more triangular patches at the same time, after storing the edges, get the set of such single edges, and find out from the set that can be connected end to end The subset of is the hollow contour;

Then, triangulate the area to form a new triangular patch and fill the holes;

Finally, by calculating the angles between all triangular patches and adjacent patches, finding the abnormal patches among them, forming a set of abnormal patches, and then searching for the closed areas, these depressed areas can be automatically located, and then through The method of adding triangular faces in this area finally completes the repair of the model.

Based on the above method, the present invention also provides a system for segmentation and three-dimensional reconstruction based on medical images, including:

Image segmentation module and 3D reconstruction module,

The image segmentation module is used to process redundant information in the original medical image, and remove information within the upper and lower threshold ranges of pixel values; by increasing the contrast of the image, extract the boundary information between the bone and the background in the above image , to obtain a gradient map; after the Gaussian smoothing operation of the gradient map, the canny algorithm is used to detect the edge of the image, and then the contour is determined and the contour limit is added;

The 3D reconstruction module is used to receive the processing result of the above-mentioned image segmentation module, display the generated 3D model through the second graphical interface and process the interaction with the user.

Furthermore, the method further includes a 3D printing module, configured to receive the processing result of the 3D reconstruction module, and construct the reconstruction model in the 3D reconstruction module by layer-by-layer printing.

Based on the above, the present invention provides a visualization and a system with a graphical interface, including the segmentation and three-dimensional reconstruction system, and also includes a display,

The display is configured as:

VTK is used to realize the visualization of the 3D model after the above 3D reconstruction, and Qt is used as GUI;

And, displaying a two-dimensional image and processing interaction with a user through the first graphical interface;

The generated three-dimensional model is displayed and the interaction with the user is processed through the second graphical interface.

Based on the above, the present invention also provides a 3D printing system, which adopts the segmentation and three-dimensional reconstruction method to obtain the reconstruction model, which is characterized in that the reconstruction model is 3D printed by FDM fused deposition manufacturing method, wherein the Used in the FDM method: one of ABS acrylonitrile-butadiene-styrene copolymer and/or PLA biodegradable plastic polylactic acid

Beneficial effects of the present invention:

In the image segmentation part, the present invention proposes its own unique segmentation algorithm based on medical file (DICOM) preprocessing, and effectively completes the segmentation task. For the reconstruction part, the moving cube (Marching Cubes, MC) algorithm is optimized, so that the model in the present invention is finally greatly improved. For interactive interfaces, visualization tools enable a certain level of interaction between the user and the model. Finally, the obtained model is printed by a 3D printer, and finally a simulated body of the skeleton is obtained, so that the doctor can perform touch interaction.

Description of drawings

Fig. 1 is a schematic flow chart of the medical image-based segmentation and three-dimensional reconstruction method of the present invention;

Figure 2(a)-Figure 2(b) is a comparison of the effects of filtering discrete noise areas;

Figure 3(a)-Figure 3(b) is a comparison of the effect after filling holes;

Figure 4(a)-Figure 4(d) is a comparison diagram of the effect after hole filling; Figure 4(c)-Figure 4(d) is a comparison diagram of the effect of different segmentation algorithms;

Figure 5(a)-Figure 5(c) is a schematic diagram of the effect of the Graph Cuts improved algorithm;

Figure 6(a)-Figure 6(b) are the segmentation results using auto-graph cuts, and a schematic diagram of the segmentation effect on adjacent objects;

Figure 7(a)-Figure 7(b) are the comparison between the original image and the effect image after removing redundant information;

Fig. 8(a)-Fig. 8(b) are respectively the display effect after obtaining the gradient and the effect of stretching the gradient range for hours;

Figure 9 is a schematic diagram of the result after the canny operation;

Figure 10(a)-Figure 10(b) are schematic diagrams of the outermost contour and the segmentation effect after projecting the result back to the original image, respectively;

Figure 11 Comparison of the effect of the segmentation method;

Fig. 12 is a schematic structural diagram of the medical image-based segmentation and three-dimensional reconstruction system of the present invention.

detailed description

The principles of the disclosure will now be described with reference to some example embodiments. It can be understood that these embodiments are described only for the purpose of illustrating and helping those skilled in the art to understand and implement the present disclosure, rather than suggesting any limitation to the scope of the present disclosure. The disclosure described herein may be implemented in various ways other than those described below.

As used herein, the term "comprising" and its variations may be understood as open-ended terms meaning "including but not limited to". The term "based on" may be understood as "based at least in part on". The term "one embodiment" can be read as "at least one embodiment". The term "another embodiment" may be understood as "at least one other embodiment".

CT (computed tomography) refers to computerized tomography, which emits X-rays through computer-controlled instruments, and then uses highly sensitive instruments to measure the human body according to the absorption and transmission rates of X-rays by different tissues of the human body. Input the measurement data into the electronic computer for processing, and then the cross-sectional image of the inspected part of the human body can be displayed, so that the doctor can know the symptoms and degree of the patient's disease without cutting the wound.

The absorption and passage of X-rays will be reflected on the image by the computer for doctors to diagnose and treat. The most commonly used one is DICOM (Digital Imaging and Communication in Medicine). DICOM files are international standards for medical digital imaging and communication, and are currently the only It is a specification strictly followed by Everbright medical imaging equipment manufacturers. It mainly consists of DICOM file header and pixel data. The pixel data part stores 12-bit grayscale image information, namely -2048~2048. Since ordinary computers can only display 28, that is, 255 pixel depths, certain conversions are required to display on our screens. In this paper, we use a new algorithm to process CT images and display them on the screen.

The current main segmentation methods include threshold segmentation, snake contour model, Level Set model, Chan-Vese segmentation, Graph Cut and other methods. Since the threshold segmentation is the simplest, the author mainly introduces the following segmentation methods.

2.2.1 Snake model

The snake model must first give an initial curve near the area of interest, then define the energy function, and then continue to converge, so that the curve deforms in the image and continuously approaches the target contour. The Snake model proposed by Kass et al. consists of a set of control points:

v(s) = [x(s), y(s)], s ∈ [0, 1]

Then, an energy function (reflecting the relationship between energy and contour) is defined on these control points:

${\mathrm{<span\; class="notranslate">\; E.</span>}}_{\mathrm{<span\; class="notranslate">\; t</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; t</span>}\mathrm{<span\; class="notranslate">\; a</span>}\mathrm{<span\; class="notranslate">\; l</span>}}<span\; class="notranslate">\; =</span>{<span\; class="notranslate">\; \&Integral;</span>}_{\mathrm{<span\; class="notranslate">\; the\; s</span>}}<span\; class="notranslate">\; \{</span>\mathrm{<span\; class="notranslate">\; \&alpha;</span>}<span\; class="notranslate">\; |</span>\frac{<span\; class="notranslate">\; \&part;</span>}{<span\; class="notranslate">\; \&part;</span>\mathrm{<span\; class="notranslate">\; the\; s</span>}}\mathrm{<span\; class="notranslate">\; v</span>}{<span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; \&beta;</span>}<span\; class="notranslate">\; |</span>\frac{{<span\; class="notranslate">\; \&part;</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{<span\; class="notranslate">\; \&part;</span>{\mathrm{<span\; class="notranslate">\; the\; s</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}\mathrm{<span\; class="notranslate">\; v</span>}{<span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; E.</span>}}_{\mathrm{<span\; class="notranslate">\; e</span>}\mathrm{<span\; class="notranslate">\; x</span>}\mathrm{<span\; class="notranslate">\; t</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; \}</span>\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; the\; s</span>}$

The first term is called elastic energy, the second term is bending energy, and the third term is external energy. As for the local image information around the control point, gradient features are often used, also known as image force. Finally, the final segmentation contour is obtained by solving the minimum value of the energy function Etotal(v).

${\mathrm{<span\; class="notranslate">\; E.</span>}}_{\mathrm{<span\; class="notranslate">\; e</span>}\mathrm{<span\; class="notranslate">\; x</span>}\mathrm{<span\; class="notranslate">\; t</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; s</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">\; v</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; \&dtri;</span>\mathrm{<span\; class="notranslate">\; I</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; )</span>{<span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}$

2.2.2 Level Set Function algorithm

The Level Set method raises some low-dimensional calculations to a higher dimension, and regards the N-dimensional description as a level of N+1 dimensions. For example, a circle on a two-dimensional plane, such as x2+y2=1, can be regarded as the level 1 of the binary function f(x, y)=x2+y2. Therefore, when calculating the change of this circle, you can first find The change of f(x,y), and then seek its 1 level set.

The level set method implicitly expresses the planar closed curve as an equivalent curve of the continuous function surface φ(x, y, t) with the same function value. Usually, the target curve is implicitly expressed in the zero level set function {φ(x,y,t)}, that is, the zero level set corresponding to time t:

$\left\{\begin{array}{c}\mathrm{<span\; class="notranslate">\; C</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; p</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; )</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">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; \&phi;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; \}</span>\\ \mathrm{<span\; class="notranslate">\; C</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; p</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</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">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; \&phi;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; \}</span>\end{array}\right.<span\; class="notranslate">\; ,</span>$

Suppose the plane closed curve used for evolution is C(p,t)=(x(p,t),y(p,t)), p is any parameterized variable, and t is time. Assuming that the inward normal vector of the curve is N and the curvature is k, the evolution of the curve along its normal vector direction can be expressed by the following partial differential equation:

$\frac{<span\; class="notranslate">\; \&part;</span>\mathrm{<span\; class="notranslate">\; C</span>}}{<span\; class="notranslate">\; \&part;</span>\mathrm{<span\; class="notranslate">\; t</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; V</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; )</span>\mathrm{<span\; class="notranslate">\; N</span>}$

Depend on Totally differentiating with respect to t, we get while the inward normal vector Organized to get:

This is the equation for curve evolution with level sets. Using the level set method to realize the active contour model has the following advantages: the evolution curve can naturally change the topological structure with the evolution of φ, and can split, merge, form sharp corners, etc.; since φ always remains a complete function during the evolution process, Therefore, it is easy to realize approximate numerical calculation; the level set method can be extended to the evolution of high-dimensional surfaces, which simplifies the complexity of three-dimensional segmentation theory and application. However, there are also many problems in the level set function, such as the choice of initialization. Since the energy function defined by the level set is not necessarily a convex function, the result is often a local minimum rather than a global minimum. Unless the initial selection position is very good, you can make the function converge to the global extremum to get better results. At the same time, the convergence process of LSF is often accompanied by a large number of iterations. This process is very time-consuming. It takes about one minute for a slightly larger image (512X 512) to iterate 1000 times. This speed is necessary for thousands of medical Graphics are not satisfactory.

2.2.3 Chan–Vese contour model

This method uses the geometric characteristics of the contour curve to establish the energy function of the contour curve movement (deformation). By minimizing this energy function, the contour curve gradually approaches the target boundary in the image, and uses the level set function to transform the contour curve motion equation into a solution Numerical partial differential equation problems. Among these methods, the Chan-Vese model proposed by Chan and Vese based on the Mumford-Shah segmentation model is a research hotspot. This model treats the original image as a simple form consisting of discontinuous sets and piecewise constant images, and the stopping function no longer depends on the local gradient of the image, but the global information of homogeneous regions.

The mathematical expression of the Mumford-Shah model is as follows: Let u(x,y) be the image function defined on the region Ω, C is the boundary of the image, divide the image into several homogeneous regions, and obtain the segmented image μ0(x,y ). The M-S model is to find the correct smooth image boundary C0, so that the error between the obtained segmented image μ0MS(x, y) and the original image u(x, y) is minimized, that is, the following energy function is minimized:,

$<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; C</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}^{\mathrm{<span\; class="notranslate">\; m</span>}\mathrm{<span\; class="notranslate">\; S</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; minF</span>}}^{\mathrm{<span\; class="notranslate">\; m</span>}\mathrm{<span\; class="notranslate">\; S</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; C</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; \&mu;</span>}\mathrm{<span\; class="notranslate">\; L</span>}\mathrm{<span\; class="notranslate">\; e</span>}\mathrm{<span\; class="notranslate">\; no</span>}\mathrm{<span\; class="notranslate">\; g</span>}\mathrm{<span\; class="notranslate">\; t</span>}\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; C</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; \&lambda;</span>}\underset{\mathrm{<span\; class="notranslate">\; \&Omega;</span>}}{<span\; class="notranslate">\; \&Integral;</span><span\; class="notranslate">\; \&Integral;</span>}<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; u</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">\; 2</span>}}\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; x</span>}\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; +</span>\underset{\mathrm{<span\; class="notranslate">\; \&Omega;</span>}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; C</span>}}{<span\; class="notranslate">\; \&Integral;</span><span\; class="notranslate">\; \&Integral;</span>}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; \&dtri;</span>{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}{<span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; x</span>}\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; the\; y</span>}$

Among them, u is the known noise-containing image, μ ₀ is the segmented image, μ, λ, ν are the weight coefficients, Length(C) represents the one-dimensional Hausdorff measure of the boundary curve C, and the image μ _{0MS is} this energy function the minimum solution of .

Based on MS theory, Chan and Vese proposed the Chan-Vese model in 2001. The mathematical description of the model is: Let the active contour line C divide the image I defined on the area Ω into two parts, which are respectively recorded as inside (C) and outside (C), and c1 and c2 are the curves inside and outside respectively. The average value of image grayscale, let the energy functional function of the image be:

$\begin{array}{c}\mathrm{<span\; class="notranslate">\; E.</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; \&phi;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\\ \mathrm{<span\; class="notranslate">\; \&mu;</span>}{<span\; class="notranslate">\; \&Integral;</span>}_{\mathrm{<span\; class="notranslate">\; \&Omega;</span>}}\mathrm{<span\; class="notranslate">\; \&delta;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&phi;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; \&dtri;</span>\mathrm{<span\; class="notranslate">\; \&phi;</span>}<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; x</span>}\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; +</span>{<span\; class="notranslate">\; \&Integral;</span>}_{\mathrm{<span\; class="notranslate">\; \&Omega;</span>}}\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&phi;</span>}<span\; class="notranslate">\; )</span>\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; x</span>}\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; \&Integral;</span><span\; class="notranslate">\; \&Integral;</span><span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; I</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{<span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&phi;</span>}<span\; class="notranslate">\; )</span>\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; x</span>}\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; \&Integral;</span><span\; class="notranslate">\; \&Integral;</span><span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; I</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>\\ {\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{<span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&phi;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span>\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; x</span>}\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; the\; y</span>}\end{array}$

length(C) is the length of the boundary curve C, Aera(C) is the area of the inner area of the curve C, μ, v≥0, λ1, λ2>0 are the weight coefficients, the first two items are called "smooth items", the control curve Maintain a certain degree of smoothness during the evolution process; the latter two items are called "fitting items", which mainly move the segmentation curve to the edge of the image to minimize the fitting error. The position of the final segmentation contour line C and the unknown constants c1, c2 are obtained by minimizing the above energy functional. In order to obtain the expression of the energy functional under the level set, the level set function φ can be defined as follows:

$\begin{array}{c}\mathrm{<span\; class="notranslate">\; E.</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; C</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; \&mu;</span>}\mathrm{<span\; class="notranslate">\; L</span>}\mathrm{<span\; class="notranslate">\; e</span>}\mathrm{<span\; class="notranslate">\; no</span>}\mathrm{<span\; class="notranslate">\; g</span>}\mathrm{<span\; class="notranslate">\; t</span>}\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; C</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; .</span>\mathrm{<span\; class="notranslate">\; A</span>}\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; e</span>}\mathrm{<span\; class="notranslate">\; a</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; no</span>}\mathrm{<span\; class="notranslate">\; the\; s</span>}\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; e</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; C</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; \&Integral;</span>{<span\; class="notranslate">\; \&Integral;</span>}_{\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; no</span>}\mathrm{<span\; class="notranslate">\; the\; s</span>}\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; e</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; C</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; I</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}{<span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; x</span>}\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; the\; y</span>}\\ <span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; \&Integral;</span>{<span\; class="notranslate">\; \&Integral;</span>}_{\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; u</span>}\mathrm{<span\; class="notranslate">\; t</span>}\mathrm{<span\; class="notranslate">\; the\; s</span>}\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; e</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; C</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; I</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; c</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{<span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; x</span>}\mathrm{<span\; class="notranslate">\; d</span>}\mathrm{<span\; class="notranslate">\; the\; y</span>}\end{array}$

According to the level set solution method described above, the numerical solution of the equation is finally obtained.

2.2.4 Graph Cuts Algorithm

Different from several previous algorithms, Graph Cuts proposed by Boykov adopts a completely different theory, that is, the method of maximum flow and minimum cut in graph theory. Graph theory is an important branch of mathematics, which first appeared in Euler's 1736 classic Königsberg bridge problem. Since its appearance, it has been deeply studied by scientists and applied to various fields on a large scale. Since this article focuses on the segmentation problem in vision, it mainly focuses on the application of graph theory in computer vision, especially the use of graph cuts in image segmentation. In some literatures, images are optimally partitioned into K parts such that the maximum cut among segmented objects is minimized. In this formulation, however, the segmentation is more biased towards small objects. Other literature attempts to solve the problem of segmentation by normalizing the cost function of cut. After optimization of the results, a non-deterministic polynomial hard problem (NP-hard) is obtained. Graph cuts have been used in other literature to minimize energy functions, which have been applied to image inpainting, 3D object reconstruction, stereo reconstruction, and many other problems in computer vision.

Let me introduce the basic principle of Graph Cuts: an undirected graph G=<V, E> is defined as a collection of nodes (vertex V) and undirected edges (E) connecting these vertices. Each pixel of an image is regarded as a node, and the connection between pixels is regarded as a set of undirected edges. In the graph, each edge e∈E is assigned a non-negative weight we. There are two special nodes called endpoints. Every cut is an edge A subset of , so that the endpoints can be separated. Generally speaking, the cost of the cut set is defined through joint optimization, that is, the cost of the edges it divides:

$<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; C</span>}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; e</span>}<span\; class="notranslate">\; \&Element;</span>\mathrm{<span\; class="notranslate">\; C</span>}}{<span\; class="notranslate">\; \&Sigma;</span>}{\mathrm{<span\; class="notranslate">\; w</span>}}_{\mathrm{<span\; class="notranslate">\; e</span>}}$

The cost function used by Boykov is very similar to the conclusions drawn in MAP-MRF. Any set of data elements is defined as P, and some adjacent areas are named point sets N by some unordered pairs {p,q}. For example the set P contains pixels for a 2D grid or volume data for a 3D grid, while, in a standard 8-bit, 26-bit contiguous system, N could contain unordered pairs of adjacent pixels. A=(A1,...,Ap,...,A|p|) is a two-dimensional vector, and its internal member Ap is assigned to each pixel p of the point set P. Each Ap can be both an "object" and a "background". Vector A defines a segmentation structure. Soft constraints based on boundary and region properties are described by a cost function. The function E(A) is as follows:

E(A)＝λ·R(A)+B(A)

Where, R(A)=∑ _p R _P (A _P )/B(A)=Σ _{{p, q}} B _{{p, q}} σ(A _P , A _q ),

And: σ(A _P , A _q )=1 (if Ap, Aq are not equal), otherwise 0.

Among them, the coefficient λ>=0 defines the importance of the region item R(A) relative to the boundary item B(A). The region term R(A) represents the penalty term for assigning pixels to "object" and "background", corresponding to RP("object"), RP("background"). For example, RP(·) reflects the degree to which the intensity of the pixel P fits the object or background intensity model.

B(A) contains the boundary terms that split A. The coefficient B{p,q}>=0 can be understood as a penalty item for disconnection of p and q. Normally, when the pixels p and q are similar, B{p,q} is very large; when the intensity difference is large, B{p,q} is close to 0. The penalty term B{p,q} can be p, A decreasing function of the q-distance. The cost B{p,q} can be based on intensity gradient, Laplacian zero connection, gradient direction and some other conditions. As for the weights of each side, we use Boykov's method in this paper, namely:

$\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\underset{\mathrm{<span\; class="notranslate">\; p</span>}<span\; class="notranslate">\; \&Element;</span>\mathrm{<span\; class="notranslate">\; P</span>}}{\mathrm{<span\; class="notranslate">\; max</span>}}\underset{\mathrm{<span\; class="notranslate">\; q</span>}<span\; class="notranslate">\; :</span><span\; class="notranslate">\; \{</span>\mathrm{<span\; class="notranslate">\; p</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; q</span>}<span\; class="notranslate">\; \}</span><span\; class="notranslate">\; \&Element;</span>\mathrm{<span\; class="notranslate">\; N</span>}}{<span\; class="notranslate">\; \&Sigma;</span>}{\mathrm{<span\; class="notranslate">\; B</span>}}_{<span\; class="notranslate">\; \{</span>\mathrm{<span\; class="notranslate">\; p</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; q</span>}<span\; class="notranslate">\; \}</span>}$

In this way, the present application constructs a complete graph and obtains its energy function. The next job is to quickly calculate the minimum cut (maximum flow) in the graph to obtain the segmentation result. Another feature of Graph Cuts is that when users find that the segmentation results are not accurate enough, they can manually modify the seed points to improve the accuracy. Moreover, the algorithm does not need to recalculate the maximum flow and minimum cut of the entire graph from scratch, but adds certain constants to the region items in the original calculation process, which greatly improves the calculation efficiency and achieves almost real-time interactive processing.

2.3 Classic algorithm improvement and new algorithm

2.3.1 Improvements to Graph Cuts

Graph Cuts is an interactive segmentation method as described in 2.2.4. For multi-target segmentation, we need to manually specify the desired segmentation target, which is understandable, and sacrifice a certain amount of time efficiency in exchange for higher segmentation accuracy. The situation is also necessary. However, there are often thousands of medical images, especially CT images. It is unimaginable to manually calibrate such a large number of images one by one. Therefore, we can no longer stay in the simple interactive semi-automatic segmentation method, but need to explore a new fully automatic Graph Cuts algorithm.

Here we first observe medical images. Since the density of most bones is much greater than other parts of the human body, it can be guaranteed that most of the pixels above a certain value will be the target we want to segment---bones, so that we can Find the seed point of the segmentation target through the corresponding threshold, as the hard limit of Graph Cuts. Similarly, after a lot of experimental observations, we found that all pixels smaller than a certain value are basically not the areas we want to segment. Therefore, we can set another threshold, and all pixels smaller than this value are regarded as the seed points of the background, which is also the hard limit part of Graph Cut. Then the final segmentation result is obtained through the traditional Graph Cut algorithm. Figure 5(a) is the original picture. In Figure 5(b), the bright spots are the seed points of the object area, and the gray points are the seed points of the background area. It can be seen that all the seed points are basically marked correctly. Figure 5(c) is the result obtained after segmentation, and the segmentation effect is not bad.

2.3.2 A new method based on image preprocessing and Canny operator

Although the Graph Cuts method designed in this application can segment bones well, its effect is often unsatisfactory when some edges are blurred, there are too many objects to be detected, and the distance between them is too close. It is very likely that the complete edge cannot be detected or the adjacent objects are merged into one, the effect is shown in Figure 6(b). Therefore, the present application proposes a new method to segment bones more accurately.

Since the original medical image is a 12-bit grayscale image, it contains a lot of useless redundant information, and it cannot be displayed on our normal display screen. So here we take a series of operations to remove redundant information, highlight the contrast between the object and the background, and finally get a better segmentation result.

1. Removal of redundant information

First, bones are much denser than other tissue parts of the body. Through a large number of experiments and careful observations, we found that most of the weak boundaries between bones and other tissues are distributed within a certain range of pixels, that is, pixels that are too large or too small can be identified as "bones" and "background". Therefore, we use a piecewise function to retain weak boundary information and remove some information that is no longer needed:

$\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\left\{\begin{array}{cc}\mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; \&le;</span>\mathrm{<span\; class="notranslate">\; 50</span>}\\ \mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 50</span>}& \mathrm{<span\; class="notranslate">\; 50</span>}<span\; class="notranslate">\; <</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; \&le;</span>\mathrm{<span\; class="notranslate">\; 305</span>}\\ \mathrm{<span\; class="notranslate">\; 255</span>}& \mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ></span>\mathrm{<span\; class="notranslate">\; 305</span>}\end{array}\right.$

Parts smaller than 50 are automatically zeroed, that is, information with small pixel values is ignored, because these pixels do not belong to bones and are redundant information; and weak edges are mostly distributed between 50 and 305 pixels, so we keep This part of the information; pixels larger than 305 can basically be identified as bones, so set all of them to the maximum value. The processing effect of the piecewise function and the overall fitting effect are shown in Figure 7(b), which shows that it increases the contrast between the object and the background.

2. Find the gradient

The image just now can highlight the contrast between the object and the background, but what we need is the edge information of the bones, so here we can further process the image, and only extract the boundary information between the bones and the background. Therefore, we use the method of obtaining the image gradient to obtain the edge information. However, the size of the gradient only occupies a small part of the 8-bit image, so here we stretch the value range of the image to further increase the contrast of the image. The effect is shown in Figure 8(b).

3. Gaussian smoothing and canny edge detection

On the basis of the gradient map obtained above, we first perform Gaussian smoothing. Because in the end we hope to get the complete object outline, that is, the closed point set. We need to smooth the image, because smoother smoothing produces rounder edges. Figure 9 The result after canny operation.

Then use the canny algorithm to detect the edge of the image. The main process of Canny is to first use the sobel difference operator to obtain the x and y direction derivatives of the grayscale image, and to obtain the gradient size and direction of each point of the image. Then set the high and low thresholds, the gradient greater than the high threshold is a strong edge pixel, and the gradient greater than the low threshold is a potentially weaker edge point. Non-maximum suppression is performed along the gradient direction among the strong edge points left after the first screening, and the adjacent weak edge points are found along the strong edge points after the second screening to obtain the final edge.

4. Find contours and add constraints

On the basis of the above canny edge, we need to find out the closed contour. If there are concentric contours, the outermost contour should be chosen as the splitting edge instead of the inner contour. Finally, we need to add some restrictions to the contour, such as the area of the contour should be greater than a certain fixed value, so as to filter some noise contours outside the main contour, so that the obtained results are more accurate. Figure 10(b) is the effect after the edge is displayed on the original image.

2.4 Experimental comparison

Among the methods described above, we choose the Chan-Vese algorithm, the improved Graph Cuts algorithm (auto-Graph Cuts), and the algorithm (IP+Canny) that combines image preprocessing and Canny operator proposed by ourselves for comparison. test. Here select different parts of human bones (toes, soles, leg bones and pelvis). One thing to note is that in auto-Graph Cuts, we use the method of obtaining the maximum flow that was later improved by Yuri et al., which effectively improves the time and space efficiency of the program. The specific comparison results are shown in Figure 11.

The following conclusions can be drawn from the above:

●The Chan-Vese algorithm is not suitable for obvious boundary segmentation, and it tends to converge on weak boundaries (the junction of flesh and air in medical images), rather than the result we want. Therefore, the segmentation results of the first three pictures are very poor, and the last one has a weak boundary between bones and muscles, but the segmentation effect is at the forefront of all algorithms. However, its algorithm is the most time-consuming because it requires multiple iterations;

● Auto-Graph Cut performs well in most cases, and its time complexity is the lowest. Of course, this does not calculate the time to provide the template, and the optimized maximum flow calculation method is used here, and this algorithm has been applied to actual products. However, for weak boundaries, its resolution ability is still weak, and for objects that are closer together, it tends to classify these objects into one category, thus automatically ignoring the weak boundary in the middle;

●IP+Canny algorithm, that is, a new algorithm proposed by ourselves. Its segmentation effect is stronger than auto-Graph Cuts when segmenting weak boundaries, and it is much better than Chan-Vese when segmenting strong boundaries. And the running time is also less, which meets our requirements for fast segmentation. And it is more suitable for the segmentation of the whole body skeleton.

3. About 3D model reconstruction

The 2D image sequence generated by medical imaging equipment already contains 3D information of human tissues and organs. The analysis and processing of 2D tomographic images can realize the segmentation and extraction of human organs and lesions, and reconstruct the 3D models of tissues, organs and lesions. , realizing the sectioning, window opening and three-dimensional display of the model will provide doctors with an intuitive technical means to help them accurately determine the location, shape and spatial relationship between the lesion and the surrounding tissues and organs, so as to design a precise Improve the effectiveness and accuracy of diagnosis and treatment.

3.1 Visualization of 3D spatial data field

1. Data type

Algorithms for 3D spatial visualization are closely related to data. The data type includes two levels of meaning: one is the type of the data itself. In the visualization of scientific computing, there are three types of data: scalar, vector and tensor; it is the type of data distribution and connection relationship. The objects of three-dimensional spatial data visualization include both the results of computer science calculations and the measurement data of measuring instruments. In scientific computing, the characteristics of the research object are often described by a set of equations, and generally only the numerical solutions of these equations can be obtained. For this reason, it is necessary to discretize the space defined by them, discretize into solid units, surface units, line segments or grid points, and then use the numerical solution method to obtain the function values at these discrete units. Therefore, the resulting data of scientific computing tends to be discrete rather than continuous. As for spatial measurement data, such as meteorological monitoring data, CT or MRI scan data of the human body, etc., they are usually also discrete. The reason is that it is difficult for people to obtain continuous measurement data in space. Therefore, the visualized objects are generally spatially discrete three-dimensional data.

The connection relationship between discrete data in three-dimensional space can be divided into three types: structured data, unstructured data and unstructured mixed data. Structured data refers to discrete spatial data that is logically organized into a three-dimensional array, and each element has its own layer number, row number, and column number. Structured data can be divided into regular grid structured data and irregular grid structured data. Regular grid structured data is distributed on a three-dimensional grid of points consisting of cubes or cuboids. The regularized volume data we refer to generally refers to structured regular volume data. Regular volume data is scalar data or vector data defined on a three-dimensional grid.

2. Definition and characteristics of medical volume data

Medical volume data generally refers to scalar or vector structured data defined on three-dimensional spatial grids. These grids are usually orthogonal grids. Data are generally defined on grid nodes, and eight adjacent grid nodes form a cube. Each small grid is the basic unit of volume data. It is assumed that the medical image sequence obtained by CT is put into a three-dimensional coordinate system after image registration and image interpolation, and each image is placed on a plane parallel to the XY plane, and the coordinates in the Z direction of the plane are set by CT when the image is generated. Relevant parameters determine that when the distance between two images is greater than the resolution of a single image, interpolation is performed between the two images, thus forming the above-mentioned volume data set, in which the volume data value is the original data value of the image.

Medical volume data is discrete volume data with optimal structure and regularity. It has the following characteristics: the data field corresponding to the image generated by CT is uniform and regular, and the value of the original data represents the nature of human tissue at a certain point (such as attenuation coefficient, proton density, etc.). The vertical resolution of the data field is generally lower than the horizontal resolution. For existing equipment, a two-dimensional resolution of 512X512 can generally be obtained, but the number of tomographic images per unit length is limited. Volume data is huge. If a group of 1000 tomographic image sequences of 512X512 are stored in 12 gray levels, it will take about 500MB. The image contains noise. The performance of imaging equipment and the user's operation level will affect the image quality, which will inevitably bring noise. In the visualization process, these specific characteristics of medical images must be considered, and corresponding measures should be taken in the processing process and specific algorithms to obtain the best results. good effect.

The reconstruction algorithms of medical images are mainly divided into two categories: volume reconstruction and surface reconstruction.

Volume reconstruction is a method of directly projecting voxels to a display plane, which is called a volume rendering method based on volume data, also known as a direct rendering method. The central idea of volume rendering technology is to specify an opacity (Opacity) for each voxel, and consider the transmission, emission and reflection of each voxel to light. The transmission of light depends on the opacity of the voxel; the emission of light depends on the objectness of the voxel (Objeemess), the greater the objectness, the stronger the emitted light; The angle relationship of the incident light. The most common method is the ray-casting method. Its main principle is to start from each pixel on the screen and emit a ray along the set viewpoint direction, and this ray passes through the three-dimensional data field. Select several equidistant sampling points along this ray, and perform trilinear interpolation from the color and opacity values of the eight voxels nearest to a certain sampling point to obtain the opacity and color values of the sampling point. Calculate the color value and opacity value of all sampling points on the ray, so as to calculate the color value at the pixel point on the screen. Ray casting experiments are also carried out in our experiments.

It can be seen that it can not only see the human bone structure, but also the muscle tissue wrapped around the bone. Although volume reconstruction has the ability to reproduce the internal structure of objects, it is not suitable for our goal. First of all, we only need to reconstruct human bones to prepare for 3D printing, so we don’t need to see the specific structure of other tissues and organs in the human body; secondly, body reconstruction requires extremely high hardware, which takes up a lot of memory and slow processing speed. Computers are not up to its operational requirements and are therefore unsuitable for real-time interaction with users. So in the end we adopted the method of surface reconstruction.

Surface reconstruction, compared with volume reconstruction, occupies less memory and has high rendering efficiency, and can basically meet our basic requirements for reconstructing human bones. Among the several methods belonging to surface reconstruction, we choose the most classic MC (MarchingCubes) algorithm as the basic algorithm, and then make certain improvements on the basis of it, so as to achieve better results.

Among the surface reconstruction algorithms, the isosurface reconstruction algorithm is the most classic. The purpose of our surface reconstruction is to construct a three-dimensional geometric model of the corresponding tissue or organ using the regions extracted by segmentation. The construction of isosurface is one of the commonly used methods to restore the 3D geometric model of objects from volume data. If we regard the volume data as a sampling set of a certain physical attribute in a certain spatial region, and the values at non-sampling points are estimated by interpolation of adjacent sampling points, then the values of all points with the same value in the spatial region The set will define one or more surfaces, called isosurfaces. Because different substances have different physical properties, appropriate values can be chosen to define isosurfaces that represent the interface of different substances. That is, an isosurface defined with appropriate values can represent the surface of a substance. An isosurface is a collection of all points with the same value in space, which can be expressed as {(x, y, z), f(x, y, z)=c} where C is a constant. However, not every voxel has an isosurface. When the corners of a voxel are greater than C or less than C, there is no isosurface inside. Value surface, we call such voxels as boundary voxels. The part of the isosurface in a boundary voxel is called the isosurface patch of the voxel. The isosurface is a cubic surface, and the intersection line between it and the boundary voxel surface is a hyperbola, and this hyperbola is only composed of the Determined by the four corner points on the surface. These isosurfaces have isotopological consistency, that is, they can form a continuous surface without holes or suspended surfaces (except at the boundary of volume data). Because for any two voxels whose boundaries are coplanar, if there is an intersection line between the isosurface and their common surface, the intersection line is the intersection line between the isovalue patch and the common surface in the two boundary voxels, that is to say, The two isosurfaces match completely, so it can be considered that the isosurface is a continuous surface composed of many isosurfaces.

The basic assumption of the MC algorithm is that the data field along the side of the cube changes in a continuous linear shape, that is to say, if the two vertices of a side are greater than or less than the value of the isosurface, only one point on the side is The intersection of this edge with the isosurface. Determining the distribution of cubic voxel isosurfaces is the basis of the algorithm.

First of all, we can regard the processed image slice data as composed of some grid points, which represent density values. Read out two slices at a time to form a layer (Layer). The eight points corresponding to the upper and lower sides of the two slices form a Cube, also called Cell, Voxel, etc. The 8 vertices of the cube are composed of 4 pixels on adjacent layers, and these 8 pixels form a cube. We call this cube a voxel. In order to determine the subdivision method of the isosurface in the voxel, a threshold value is required for the isosurface, and then the eight vertices of the voxel are classified to determine whether the vertex is located in the isosurface or in the isovalue surface; and then determine the subdivision mode of the isosurface according to the vertex classification results. Vertex classification rules are:

If the data value of the vertex is greater than the value of the isosurface, it is defined that the vertex is inside the isosurface, which is recorded as "1"; if the vertex density value < the threshold value, it is set as Outside (1).

If the data value of the vertex is less than the value of the iso-surface, it is defined that the vertex is outside the iso-surface, recorded as "0". Vertex density value ≥ threshold value, Inside(0).

First determine which voxels the isosurface passes through, and then determine how the isosurface intersects the voxels. When the values of some pixels in a voxel are greater than the threshold, while other pixels are less than the threshold, then the isosurface must pass through this voxel. If the values of 8 pixels in a voxel are all less than the threshold or all are greater than the threshold, then the The voxel does not intersect the isosurface, and the isosurface does not pass through the voxel. When a voxel intersects the isosurface, there must be some pixel values greater than the threshold and some less than the threshold. Each pixel has two states, either greater than the threshold or less than the threshold to determine the voxel containing the isosurface. For a voxel whose 8 corner points are all 1 or 0, it belongs to the "0" structure and no isosurface passes through the voxel. When there is a corner point marked as 1, it is the No. 1 structure. We use a triangular patch to represent the isosurface, which divides the corner point and other seven corner points into two parts. For the remaining configurations, multiple triangles will be generated. flag (i, j, k) = 0 (1) 256 cases. Therefore, there are 256 combined states in total. Each combination corresponds to a case where the isosurface intersects the voxel. Because 8 points have a symmetrical relationship, 256 combinations can be simplified to 15 situations. Each relationship corresponds to how the isosurface intersects the voxel. After knowing how the isosurface intersects the voxel, the intersection point between the isosurface and the cube edge can be obtained. The patch formed by these intersection points is a part of the isosurface . When all the voxels intersecting the isosurface are found and the corresponding intersecting surfaces are obtained, the isosurface is also obtained. Finally, connecting these isosurfaces constitutes the reconstructed object.

A brief review of its algorithm process, the main steps are as follows:

(1) Construct an index table of 256 intersecting relationships according to the symmetrical relationship. The table indicates which edge of the voxel the isosurface intersects with;

(2) Extract adjacent 8 of the adjacent two-layer pictures to form a voxel and number these 8 pixels;

(3) Determine whether the pixel is 1 or 0 according to the comparison between each pixel and the threshold;

(4) The 01 string formed by these 8 pixels is formed into an 8-bit index value;

(5) use the index value to look up the corresponding relationship with the pixels in the index table on the top, and find the point with each side of the cube;

(6) Use the intersection points to form a triangular facet or a polygonal facet;

Traverse all voxels of the 3D image, and execute (2) to (6) repeatedly.

Improvement of MC Algorithm

After MC reconstruction, we can already get a fairly good reconstruction effect. However, there are still many flaws in the reconstructed 3D model, so we need to make a series of optimizations to get better results.

1. Use triangular strips to connect scattered triangular faces

The 3D model is composed of many triangular faces, but too many triangular faces will occupy a large amount of memory, which is a huge challenge for the PCs we usually use. Therefore, here we use the method of connecting triangle patches into triangle strips to reduce memory usage. As shown in the figure above, a single triangle patch needs to store the coordinates of three vertices, then N triangles need to store 3N points. And if these triangle faces are connected into a triangle belt, only N+2 points are needed. Save 2/3 of the memory. This is a huge improvement for our skeletal models that often have hundreds of thousands or even millions of triangular faces.

2. Reduce the number of triangles

Although we generate triangle strips by connecting triangle patches, the memory usage is greatly reduced. However, it is worth noting that we hope that the final 3D model can be successfully 3D printed into a solid body, but the accuracy of 3D printers is limited (generally 0.1mm), which means that although we have a very fine 3D model, but 3D printers cannot recognize such fine structures. Therefore, adhering to the idea of "not seeking the best, only seeking the most suitable", we greatly reduce the triangular faces here. The benefits of doing so also include reducing time and space complexity, which is beneficial to subsequent operations.

Here we use the method of William (Schroeder W J, Zarge J A, Lorensen W E. Decimation of triangle meshes[C]//ACM Siggraph Computer Graphics. ACM, 1992,26(2):65-70.) to first determine each triangle Which type of patch belongs to, and then according to the attributes of each type of triangular patch, such as the distance between the simple plane and the ground, the distance between the boundary type and the edge, etc., it is determined which triangle to delete, and then the generated hole is positioned and filled. These holes are trimmed until the user-set reduction ratio is reached. It is the effect diagram achieved by different reduction ratios. However, when the setting value of this method is too high, it is easy to change the topology of the model, and even produce holes. The solution is to try to avoid setting too high a reduction ratio, while filling holes in the next section.

3. Filter discrete noise regions

We can find that in addition to the bone part we want, there are some other impurities in the picture, such as: CT detection bed, some noise left over from the segmentation stage. These parts are all parts we don't need, so we need to filter them out through the algorithm. Here we mainly use connectivity detection, first arrange all the triangles into a sequence, then take the first triangle as the root, and then gradually find all the connected triangles as a whole. Eventually the whole model is divided into several disjoint regions. Then calculate the number of triangles in each area, and select the two largest areas to get the bone model we want, as shown in Figure 2(a)-Figure 2(b), the connection parts in the model, Fig. 2(a) is to mark different parts, and Fig. 2(b) is the filtered result.

4. Fill holes

After we get the regions we want, we start refining our model and reducing the holes in it. We first need to detect holes, here we use this algorithm: first traverse all edges, if an edge is not used by two or more triangles at the same time, then store this edge. Afterwards, we obtained this kind of unilateral set, and found the subset that can be connected end to end, which is the empty outline. This area is then triangulated to form a new triangular patch to fill the hole. Specifically, as shown in Figure 3(a)-Figure 3(b), it is a comparison picture before and after hole filling. However, we found that there are still some suspected holes in the model, but they have not been detected or filled. By meshing the triangular facets, we found that these areas are not holes, but some concave structures, which may be related to the insufficient fineness of the segmentation results when we segmented. After a lot of experiments, we found that the angles between the triangular patches at the edges of these concave regions are very different. Therefore, by calculating the angles between all triangular patches and adjacent patches, finding out the abnormal patches, forming a set of abnormal patches, and then searching for the closed areas, we can automatically locate these depressed areas. Afterwards, by adding triangular patches in this area, the model is finally repaired. The details are shown in Figure 4(a)-Figure 4(b). Figure 4(a) is the detected concave region, and Figure 4(b) distinguishes different concave regions.

Effect of Different Segmentation Algorithms on Reconstruction Effect

After completing the proposal and improvement of the reconstruction algorithm, we began to combine the previous segmentation results for the final 3D reconstruction. We further process the previous image segmentation, set all the segmented target area pixels to "1", and all the background area pixel values to "0", and extract the target area to obtain point cloud data.

Since Chan-Vese is not effective for bone data segmentation, we mainly compare auto-GraphCuts and our own IP+Canny segmentation results for reconstruction. After reconstruction, the results can be seen in Figure 4(c)-Figure 4(d) As shown, Figure 4(c) is the reconstruction result of auto-Graph Cuts, and Figure 4(d) is the reconstruction result after IP+Canny segmentation. We can find that the reconstruction results of auto-Graph Cuts are very bad, and the middle fault is serious, while the overall effect of the reconstruction results of IP+Canny is good. We zoomed in on some areas of the two to observe their internal structures in detail.

Visualization and Graphical Interface (GUI)

After completing the segmentation and reconstruction process, we need to integrate the two parts and display them on the screen to complete certain interactions with the user. Here we mainly use VTK to realize the visualization of the 3D model, Qt is used as the GUI of the whole software, and then we divide the interface into two parts, the left side is used for 2D image display, processing and interaction, and the right side is used to display the generated 3D model and implement simple interactions. The main functions of the two-dimensional image display processing part on the left are:

1. Read and display two-dimensional medical images, and at the same time display the coordinates and gray value of any point;

2. Obtain the gradient of the medical image and display it;

3. Display the gray histogram of the image and fit it with a curve;

4. Realize simple image editing.

The main function on the right side is to allow the user to select the medical file to be processed, and then display the reconstructed 3D model, and perform interactive operations such as translation and zooming.

3D printing principle

3D printing is a kind of rapid prototyping technology. It is a technology based on digital model files and using bondable materials such as powdered metal or plastic to construct objects by layer-by-layer printing.

3D printing technology appeared in the mid-1990s and is actually the latest rapid prototyping device utilizing technologies such as light curing and paper lamination. The working principle of it is basically the same as that of ordinary printing. The printer is filled with "printing materials" such as liquid or powder. After being connected to the computer, the "printing materials" are superimposed layer by layer through computer control, and finally the blueprint on the computer is turned into a real object. This printing technology is called 3D printing technology.

Nowadays, 3D printers mainly use the following technologies for rapid prototyping: layered solid manufacturing (LOM): cutting foil and paper according to the layered geometric information of the part, and bonding the obtained layers into a three-dimensional solid. Stereolithography (SLA): Use a laser with a specific wavelength and intensity to focus on the surface of the photo-curable material to solidify it sequentially from point to line and from line to surface to complete a layer of drawing work, and then the lifting table is vertically Move the height of one ply before solidifying another. In this way, layers are superimposed to form a three-dimensional entity. Selective Laser Sintering (SLS): Formed using powdered materials. Sprinkle the material powder on the upper surface of the formed part and smooth it. Use a high-intensity CO2 laser to scan the section of the part on the newly laid layer; the material powder is sintered together under high-intensity laser irradiation to obtain the section of the part and connect it with the formed part below. After sintering, a new layer of powder material is laid, and the lower section is selectively sintered until the whole is printed. Fused deposition manufacturing (FDM): The heating head heats the hot-melt material (ABS resin, nylon, wax, etc.) to a critical state, showing semi-fluid properties, and under the control of the computer, moves along the two-dimensional geometric information determined by CAD. The extruder extrudes the material in a semi-fluid state and solidifies to form a thin layer of contour shape. When one layer is completed, the new cambium is lowered by the vertical lifting system for curing. In this way, layers are stacked and bonded to form a three-dimensional solid of a part from bottom to top.

The printing process includes three parts: first design the model, design the model we need through a large number of 3D modeling software such as Ug and auto-cad, and get the corresponding files. Then slice the model file layer by layer, determine the material point of each layer, and plan the routing method, and finally generate a Gcode file, that is, a machine code that the printer can recognize. Finally, the 3D printer accepts the instruction, the nozzle starts to move slowly, the platform starts to descend layer by layer, and finally the entire model is printed.

Validation of the reconstructed model

After introducing the basic principles of 3D printing, we began to print our reconstructed model through 3D printing technology, so as to verify the effectiveness of our model construction through observation and touch.

Because FDM is easy to use and maintain, the operating environment is clean, the material cost is low, and the printed model is relatively hard, so here we mainly use the FDM method for 3D printing. FDM mainly uses two kinds of materials: ABS (acrylonitrile-butadiene-styrene copolymer) and PLA (biodegradable plastic polylactic acid). parameters, and record their respective fits. Through a lot of experiments, we found that the PLA material is better, PLA has no pungent smell of ABS, and generally does not warp corners. In addition, PLA has a low shrinkage rate, and it still performs well even if the size of the model is large. At the same time, the melt strength is low, the printed model is easier to shape, the surface gloss is excellent, and the color is bright. Therefore, in the final bone printing, we also used PLA materials for experiments.

In the experiment, we segmented and reconstructed the sternum and leg bones, and then used a 3D printer to print out the model entities of these two parts. Compared with the model, our printing effect is similar in size, shape and tiny details. It also verified the feasibility of our idea. Although the surface of the printed model is relatively rough, according to the previous experience of printing other samples, this roughness is caused by the accuracy and stability of the printer, so our reconstructed model can basically meet certain accuracy requirements. Of course, 3D printing can not only verify the effect of our reconstruction, but also has many extensive applications in the medical field, and its possible future applications:

By simulating the real environment inside the human body, formulate surgical plans; design human organs, print out bionic ears, bionic eyes and artificial limbs and other personalized appliances through special materials, so as to better help the disabled;

Due to the touchability and detailed representation of the 3D model, it can be used in medical research to help researchers carry out experiments more effectively; the 3D model can be more intuitively displayed to medical students, thereby improving the level of teaching and training.

Claims (10)

1. A segmentation and three-dimensional reconstruction method based on medical image, is characterized in that, comprises the steps: Process the redundant information in the original medical image, and eliminate the information within the upper and lower threshold range of the pixel value; By increasing the contrast of the image, the boundary information between the bone and the background in the above image is extracted to obtain a gradient map; After Gaussian smoothing operation of the gradient map, the canny algorithm is used to detect the edge of the image, and then the contour is determined and the contour limit is added; The image obtained after the above image segmentation is subjected to three-dimensional reconstruction through the MC algorithm. 2. segmentation and three-dimensional reconstruction method according to claim 1, is characterized in that, also comprises the operation step of following visualization and graphical interface: VTK is used to realize the visualization of the 3D model after the above 3D reconstruction, and Qt is used as GUI; And, displaying a two-dimensional image and processing interaction with a user through the first graphical interface; The generated three-dimensional model is displayed and the interaction with the user is processed through the second graphical interface. 3. segmentation and three-dimensional reconstruction method according to claim 1, is characterized in that, when carrying out three-dimensional reconstruction by described MC algorithm, comprises following optimization steps: Use triangular strips to connect scattered triangular patches to reduce memory usage; According to the coordinates of three vertices that need to be stored in a single triangular patch, N triangles need to store 3N points, If these triangular faces are connected into a triangular strip, N+2 points are needed. 4. according to claim 1 or 3 described segmentation and three-dimensional reconstruction method, it is characterized in that, when carrying out three-dimensional reconstruction by described MC algorithm, comprise following optimization steps: Reduce the time and space complexity by reducing the number of triangle faces; First, determine which type each triangle belongs to, and then decide which triangle to delete according to the attributes of each type of triangle; Second, relocate the resulting holes and fill them in until the user-set reduction ratio is reached. 5. segmentation and three-dimensional reconstruction method according to claim 1 or 2, is characterized in that, when carrying out three-dimensional reconstruction by described MC algorithm, comprises following optimization steps: Use connectivity detection to filter discrete noise regions and remove impurities; First, arrange all the triangular faces in sequence, and then use the first triangular face as the root of the tree; Second, gradually find all connected triangles as a whole; Finally, the entire model is divided into several disjoint regions, and the number of triangles in each region is calculated, and the two largest regions are selected. 6. according to claim 1,2 or 3 described segmentation and three-dimensional reconstruction method, it is characterized in that, when carrying out three-dimensional reconstruction by described MC algorithm, comprise following optimization steps: Detect holes and refine the model by filling holes; First, traverse all the edges. If an edge is not used by two or more triangular patches at the same time, after storing the edges, get the set of such single edges, and find out from the set that can be connected end to end The subset of is the hollow contour; Then, triangulate the area to form a new triangular patch and fill the holes; Finally, by calculating the angles between all triangular patches and adjacent patches, finding the abnormal patches among them, forming a set of abnormal patches, and then searching for the closed areas, these depressed areas can be automatically located, and then through The method of adding triangular faces in this area finally completes the repair of the model. 7. A system based on medical image segmentation and three-dimensional reconstruction, characterized in that it includes: an image segmentation module and a three-dimensional reconstruction module, the image segmentation module is used to process redundant information in the original medical image and remove pixels The information within the upper and lower threshold range of the value; by increasing the contrast of the image, extract the boundary information between the bone and the background in the above image to obtain a gradient map; after Gaussian smoothing operation on the gradient map, use the canny algorithm to detect the edge of the image Afterwards, the outline is determined and outline restrictions are added; the 3D reconstruction module is used to receive the processing result of the above-mentioned image segmentation module, and display the generated 3D model through the second graphical interface and process the interaction with the user. 8. The segmentation and three-dimensional reconstruction system according to claim 7, further comprising a 3D printing module for receiving the processing results of the three-dimensional reconstruction module, and constructing the three-dimensional reconstruction by layer-by-layer printing The reconstructed model in the module. 9. A visualization and a system with a graphical interface, characterized in that it includes the segmentation and three-dimensional reconstruction system as claimed in claim 7 or 8, and also includes a display, The display is configured to: adopt VTK to realize the visualization of the three-dimensional model after the above-mentioned three-dimensional reconstruction, and use Qt as a GUI; and, display two-dimensional images and process interaction with users through the first graphical interface; display the generated results through the second graphical interface 3D model and handle interaction with the user. 10. A 3D printing system, adopting the segmentation and three-dimensional reconstruction method as claimed in claim 1 to obtain the reconstruction model, characterized in that, the reconstruction model is 3D printed by FDM fused deposition manufacturing method, wherein the FDM method Used in: one of ABS acrylonitrile-butadiene-styrene copolymer and/or PLA biodegradable plastic polylactic acid.