A kind of medical image Organizational Intelligence dividing method based on three-dimensional reconstruction
Technical field
The present invention relates to the technical fields of medical image segmentation, refer in particular to a kind of medical image group based on three-dimensional reconstruction
Knit intelligent scissor method.
Background technique
Medical image segmentation refers to according to a certain specific character, will have the area of different meanings in medical image
Domain separates, and identical region is made to meet region consistency.The dividing method of common medical image includes: based on threshold value
Dividing method, the dividing method based on region, the dividing method based on edge and the dividing method based on cluster.
Dividing method based on threshold value is using one or a set of gray threshold, by the image data mark in threshold range
It is denoted as area-of-interest, the image data outside range is labeled as background data, achievees the purpose that image segmentation with this.Based on threshold value
Dividing method, have it is simple, quickly, the low advantage of computation complexity;But the selection of threshold value is affected by noise big, it is difficult to choose
Degree is high;Simultaneously as merely with the grayscale information of image, it is difficult to divide the high area-of-interest of gray scale similarity.And medicine
Image generally has the characteristic that intensity profile is uneven, each tissue gray scale similarity is high, and therefore, which is being used for medical image
Segmentation of Multi-target when, segmentation effect is poor.
Dividing method based on region, the methods of main inclusion region growth, split degree method, level set;Wherein, due to
It calculates simply, region-growing method is widely used in medical image segmentation.But traditional algorithm of region growing high degree
Dependent on initial seed point selection and grow criterion setting, it is difficult to realize the automatic segmentation of multizone.
Dividing method based on edge, be using between different zones, often exist more violent grey scale change this
Characteristic comes the edge between detection zone and carries out image segmentation.This method is for edge gray value transition than more significant and noise
The segmentation effect of lesser simple image is ideal;For structure is complicated and presence is compared with the medical image of very noisy, then can produce
The case where raw pseudo-edge or missing edges.
Dividing method based on cluster is to indicate the pixel in image space with corresponding feature space point, according to it
Feature space is split in the aggregation of feature space, they are then mapped back into original image image space, obtains segmentation result.
Dividing method based on cluster mainly includes following difficult point: the determination of cluster class number;Determine the validity criterion of cluster;It is poly-
The setting of class initial value;Computational complexity is larger.Traditional clustering algorithm does not consider the spatial information of image data, makes an uproar to image
Sound is very sensitive.
Currently, traditional Medical image segmentation algorithm still remains following problem:
1. being split using priori knowledge, for complicated medical image, any single partitioning algorithm is all difficult
To obtain satisfactory segmentation result.
2. most partitioning algorithm is all based on the segmentation that two-dimensional image data carries out, image is not accounted in three-dimensional
Location information in space.
3. current Algorithm on Automatic Segmentation is difficult to meet the needs of segmentation precision, therefore cutting procedure usually introduces manually
Interactive information divides low efficiency.
Summary of the invention
It is an object of the invention to overcome the shortcomings of the prior art, a kind of medicine shadow based on three-dimensional reconstruction is proposed
As Organizational Intelligence dividing method, this method can automatically extract region growing institute by statistic histogram using three-dimensional data
Threshold value needed for the initial seed point and growth criterion that need carries out region growing then in conjunction with the edge contour information of image, and
Initial seed point and growing threshold are continued to optimize by iteration growth, so as to improve the result of segmentation.Meanwhile this method is by small
Wave variation, image data was reduced to the 1/4 of original image, so as to shorten operation time.
To achieve the above object, a kind of technical solution provided by the present invention are as follows: medical image group based on three-dimensional reconstruction
Knit intelligent scissor method, comprising the following steps:
1) DICOM medical image sequences are inputted, window is pre-processed and adjusted to image;
2) Wavelet transformation is carried out to image, image data is reduced to the 1/4 of original image;
3) grey level histogram is counted, initial seed point and growing threshold are extracted;
4) three-dimensional edges detection is carried out to image, obtains edge contour figure;
5) gray scale and marginal information are combined, three-dimensional spatial area growth is carried out;
6) area attribute after segmentation is calculated, seed point and growing threshold are optimized.
In step 1), one group of DICOM format medical image sequences is inputted, image sequence is adjusted in pretreatment stage
Window and Gaussian smoothing.So-called tune window turns the gray value being located in window in image data that is, according to the window width of precognition and window position
Be changed to value when display in most bright and most dark range, be set to higher than the part of window tonal range it is most bright, be lower than window gray scale model
Part in enclosing is set to most dark.Adjust window formula as follows:
Wherein, F'xyzFor gray value of the pixel after adjusting window, F at three dimensional space coordinate (x, y, z)xyzFor the pixel
Original gray value, C is window position, and W is window width, and Max is the maximum gradation value that display is shown.
For the image sequence after tune window, Gaussian smoothing is carried out, i.e., image is weighted and averaged, each pixel
Value, the especially value of itself and other interior pixels of neighborhood obtain after being weighted averagely.
Three-dimensional Gaussian function is defined in above formula, wherein G (x, y, z) indicates the Gaussian template being calculated, (x, y, z)
The relative coordinate of other pixels and central pixel point in expression field, σ indicate the standard deviation of Gaussian Profile.
In step 2), two layers of three-dimensional discrete wavelet variation are carried out to image using " db1 " wavelet basis, that is, to image
X, Y, the filtering that Z-direction carries out low-and high-frequency respectively decomposes, and every layer of Wavelet transformation will all generate 8 components, wherein low comprising 1
Frequency component and 7 high fdrequency components continue the low frequency component that first layer generates to do filtering decomposition, finally obtain 15 components.It takes
The low frequency component generated after second layer Wavelet transformation carries out subsequent image segmentation, image data can be reduced to original image
1/4.
In step 3), the pixel grey scale of statistical picture generates grey level histogram, scans histogram, obtains each sense
The gray value and its growing threshold of interest region initial seed point, detail are as follows:
3.1) grey level histogram is generated according to the following formula, defines peak set P, valley set V;
In formula, H (i) indicates that gray scale is the percentage of the total pixel number of pixel number Zhan of i, wherein i ∈ 0,1,2,
3 ..., Max-1 }, X, Y, Z indicate image in three-dimensional space X, and Y, the pixel sum in Z-direction, x, y, z indicates that pixel exists
The coordinate of three-dimensional space, FxyzIndicate gray value of the image pixel at (x, y, z), Equalxyz(i) for judging at (x, y, z)
Pixel grey scale whether be equal to i, if being equal to, return to 1, otherwise, return 0;
3.2) the corresponding point of maximum value of histogram is found as first peak point, and recording the corresponding gray value of point is
p0, and be added among set P;
3.3) find next peak point, corresponding grey scale value is p', Ying Man: centered on p', t is the range of step-length
It is interior, that is, tonal range is in the range of [p'- t, p'] and [p', p'+t], monotonic increase and list is presented in grayscale image respectively
Adjust decrement states;
3.4) step 3.3) is repeated until finding all peak points and their corresponding gray scale value set P, from small arrival
Order rearrangement column set P in element;
3.5) set P is traversed, gray value corresponding to the smallest point of histogram value between every two peak point is sequentially found
As valley, set V is formed;
3.6) set V is traversed from small to large, in pairs as growing threshold, is found between two valleies according to set P
Peak value p, randomly choosing three gray values is the pixel of p as initial seed point.
In step 4), three-dimensional Sobel edge detection operator is designed, calculates the pixel gradient of image, threshold value is set, it is high
In the edge pixel that the pixel of threshold value is considered as image, it is set to 1, is otherwise set to 0, the edge contour of medical image is obtained with this
Figure.
In step 5), using the initial seed point and its corresponding growing threshold for counting acquisition in histogram, in conjunction with figure
The edge contour information of picture carries out the region growing of three-dimensional space, that is, is determining whether point to be grown meets appraisal of growth criterion
When, gray scale bound needed for not only meeting growth, also needing to meet rim condition, specific step is as follows:
5.1) initial seed point extracted in step 3) is pressed into seed stack, the corresponding growing threshold of record seed point;
5.2) judge whether seed stack is empty, is that sky then terminates to grow;It is not sky, then the seed point of stack top is taken out, with this
Point is initial point, searches for its six fields pixel (x, y, z), centered on (x, y, z), generates the super voxel of 3*3*3, using super body
The average gray value of element is as the gray value at (x, y, z) pixel;
5.3) judge whether the gray value of (x, y, z) meets the corresponding growing threshold of seed point, if satisfied, judging that the point is
No is marginal point, and marginal point is divided into three kinds of situations: 1. the point is not marginal point, then with seed point labeled as the same area, and at
For new seed point, it is pressed into seed stack;2. the point is marginal point, but next point in the direction of growth is still marginal point, then
It is labeled as the same area with seed point, and becomes new seed point and is pressed into storehouse;3. be marginal point, and in the direction of growth under
Some non-edge points, the direction stop growing;
5.4) step 5.2) is repeated, until growth terminates.
In step 6), after each seed point growth, the flat of the area-of-interest of its generation is calculated according to the following formula
Equal gray scale and its inter-class variance size with background area;
In formula, μROIIndicate the average gray of area-of-interest, NROIIndicate the number of pixels that area-of-interest includes, μBTable
Show the average gray of background area, NBIndicate the number of pixels that background area includes, FxyzIndicate the ash at pixel (x, y, z)
Angle value, θ indicate that the inter-class variance in right interest region and background area, μ indicate image ensemble average gray scale;
Merge the initial seed point and growing threshold of area-of-interest similar in average gray;Set minimum inter-class variance threshold
Value judges whether there is the area-of-interest that inter-class variance is less than the threshold value, and if it exists, then assert that the cut zone is meaningless,
Merge the region initial seed point and growing threshold to adjacent average gray area-of-interest;According to new seed point and life
Long threshold value regrows, until there is no meaningless cut zone or reaching highest the number of iterations, finally returns that comprising more
The segmented image of class label.
Compared with prior art, the present invention have the following advantages that with the utility model has the advantages that
1. the present invention is based on 3 d image datas to be split, the spatial positional information of image is effectively utilized, it can be with
Obtain better segmentation effect.
2. the present invention solves the disadvantage that algorithm of region growing relies on artificial selected seed point and design growth criterion, can
The gray threshold that automatic selected seed point and its growth need to meet, realizes the automation of segmentation.
3. the present invention combines the edge contour information of image on the basis of region growing, it ensure that and chosen in automation
The inaccurate situation of growing threshold under the accuracy rate that grows, improve the accuracy of segmentation.
4. the present invention carries out region growing by iteration, calculates average gray and inter-class variance improves segmentation result, into one
Step reduces automation and chooses the invalid influence of seed point and growing threshold to segmentation result in obtained seed point, improves segmentation
Accuracy.
5. the present invention reduces image data using Wavelet transformation, reduce calculation amount and runing time, so that algorithm has
There is better real-time.
Detailed description of the invention
Fig. 1 is logical flow diagram of the present invention.
Fig. 2 is the Gaussian template figure that the present invention uses.
Fig. 3 is that two layers of 3 D wavelet that the present invention uses change schematic diagram.
Fig. 4 is the Sobel edge detection operator structure chart that the present invention uses.
Fig. 5 is the region growing flow chart that the present invention designs.
Specific embodiment
The present invention is further explained in the light of specific embodiments.
As shown in Figure 1, the medical image Organizational Intelligence dividing method provided by the present invention based on three-dimensional reconstruction, including such as
Lower step:
The first step, input one group of DICOM format medical image sequences, pretreatment stage to image sequence carry out adjust window and
Gaussian smoothing.So-called tune window is by the grayvalue transition in window is located in image data that is, according to the window width of precognition and window position
Value when display in most bright and most dark range, be set to higher than the part of window tonal range it is most bright, lower than in window tonal range
Part be set to it is most dark.Adjust window formula as follows:
Wherein, F'xyzFor gray value of the pixel after adjusting window, F at three dimensional space coordinate (x, y, z)xyzFor the pixel
Original gray value, C is window position, and W is window width, and Max is the maximum gradation value that display is shown.
In this example, the brain medical image sequences of DICOM format are inputted, C=50 is arranged in size 512*512*200,
W=90, Max=256.For the image sequence after tune window, Gaussian smoothing is carried out, i.e., image is weighted and averaged, each
The value of pixel, the especially value of itself and other interior pixels of neighborhood obtain after being weighted averagely.
Three-dimensional Gaussian function is defined in above formula, wherein G (x, y, z) indicates the coefficient for the Gaussian template being calculated,
(x, y, z) indicates the relative coordinate of other pixels and central pixel point in field, and σ indicates the standard deviation of Gaussian Profile.According to
The size of the σ value of setting determines smooth effect, and value is bigger, then smooth effect is more obvious.As shown in Figure 2, it is shown that a 3*3*
σ=0.8 is arranged in the coordinate of each point in 3 Gaussian template, this example, and coordinate value each in Fig. 2 is substituted into Gaussian function, meter
Obtained G (x, y, z) value is the weight coefficient of Gaussian template.
Second step carries out the variation of two layers of three-dimensional discrete wavelet to image using " db1 " wavelet basis, that is, to the X of image,
Y, the filtering that Z-direction carries out low-and high-frequency respectively are decomposed, and every layer of Wavelet transformation will all generate 8 components, wherein including 1 low frequency point
Amount and 7 high fdrequency components continue to do filtering to the low frequency component that first layer generates and decompose, finally obtain it is shown in Fig. 3 [LLL2,
(LLH2, LHL2, HLL2, LHH2, HLH2, HHL2, HHH2), (LLH1, LHL1, HLL1, LHH1, HLH1, HHL1, HHH1)] altogether
15 components.Take the low frequency component generated after second layer Wavelet transformation, i.e., " LLL2 " component carries out subsequent image segmentation, can be with
Image data is reduced to the 1/4 of original image.In this example, the image of 512*512*200 size by two layers of three-dimensional from
After dissipating Wavelet transformation, output image size is 128*128*40.
Third step, the pixel grey scale of statistical picture generate grey level histogram, scan histogram, it is interested to obtain each
The gray value and its growing threshold of region initial seed point, detail are as follows:
3.1, grey level histogram is generated according to the following formula, defines peak set P, valley set V;
In formula, H (i) indicates that gray scale is the percentage of the total pixel number of pixel number Zhan of i, wherein i ∈ 0,1,2,
3 ..., 255 }, X=128, Y=128, Z=40, x, y, z indicate coordinate of the pixel in three-dimensional space, FxyzIndicate image pixel
Gray value at (x, y, z), if gray value is i, Equalxyz(i) it is equal to 1, the pixel that gray scale is i in corresponding H (i)
Number plus 1.
3.2, the corresponding point of maximum value of histogram is found as first peak point, and recording the corresponding gray value of point is
p0, and be added among set P;
3.3, next peak point is found, corresponding grey scale value is p', Ying Man: centered on p', step-length t=10 is set,
In the range of namely tonal range is [p'- 10, p'] and [p', p'+10], monotonic increase and dullness is presented in grayscale image respectively
Decrement states;
3.4, step 3.3 is repeated until finding all peak points and their corresponding gray scale value set P, from small arrival
Order rearrangement arranges the element in set P;
3.5, set P is traversed, gray value corresponding to the smallest point of histogram value between every two peak point is sequentially found
As valley, set V is formed.
3.6, set V is traversed from small to large, in pairs as growing threshold, is found between two valleies according to set P
Peak value p, randomly choosing three gray values is the pixel of p as initial seed point.
4th step designs each of the Sobel edge detection operator of three-dimensional shown in Fig. 4, traversal image three-dimensional space
Pixel calculates the gradient of the point using Sobel edge detection operator.Threshold value is set, and the pixel higher than threshold value is considered as figure
The edge pixel of picture is set to 1, is otherwise set to 0, obtains the edge contour figure of medical image with this.
5th step, using the initial seed point and its corresponding growing threshold for counting acquisition in histogram, in conjunction with image
Edge contour information carries out the region growing of three-dimensional space shown in fig. 5, that is, is determining whether point to be grown meets appraisal of growth
When criterion, gray scale bound needed for not only meeting growth also needs to meet rim condition.Specific step is as follows:
5.1, the initial seed point extracted in third step is pressed into seed stack, the corresponding growing threshold of record seed point;
5.2, judge whether seed stack is empty, is that sky then terminates to grow;It is not sky, then takes out the seed point (x of stack top0,y0,
z0), using the point as initial point, search for its six fields pixel (x, y, z), (x, y, z) ∈ { (x0-1,y0,z0), (x0+1,y0,z0),
(x0,y0-1,z0),(x0,y0+1,z0),(x0,y0,z0-1),(x0,y0-1,z0+ 1) }, centered on (x, y, z), 3*3*3 is generated
Super voxel, using the average gray value of super voxel, i.e., centered on (x, y, z) add 26 pixels of surrounding average gray
Value is as the gray value at (x, y, z) pixel;
5.3, judge whether the gray value of (x, y, z) meets the corresponding growing threshold of seed point, if satisfied, judging that the point is
No is marginal point, and marginal point is divided into three kinds of situations: 1. the point is not marginal point, then with seed point labeled as the same area, and at
For new seed point, it is pressed into seed stack;2. the point is marginal point, but next point in the direction of growth is still marginal point, then
It is labeled as the same area with seed point, and becomes new seed point and is pressed into storehouse;3. be marginal point, and in the direction of growth under
Some non-edge points, the direction stop growing;
5.4, step 5.2 is repeated, until growth terminates.
6th step calculates the average ash of the area-of-interest of its generation after each seed point growth according to the following formula
Degree and its inter-class variance size with background area.
In formula, μROIIndicate the average gray of area-of-interest, NROIIndicate the number of pixels that area-of-interest includes, μBTable
Show the average gray of background area, NBIndicate the number of pixels that background area includes, FxyzIndicate the ash at pixel (x, y, z)
Angle value, θ indicate that the inter-class variance in right interest region and background area, μ indicate image ensemble average gray scale.
Merge the initial seed point and growing threshold of area-of-interest similar in average gray;This example is set between infima species
Variance threshold values are 5, judge whether there is area-of-interest and background area inter-class variance less than 5, and if it exists, then assert that the sense is emerging
The segmentation in interesting region is meaningless, merge the region initial seed point and growing threshold to adjacent average gray region of interest
Domain;It is regrowed according to new seed point and growing threshold, is changed until meaningless cut zone is not present or reaches highest
Generation number finally returns that the segmented image comprising multiclass label.
Embodiment described above is only the preferred embodiments of the invention, is not limited thereto and makes practical range of the invention,
Therefore all shapes according to the present invention, change made by principle, it should all be included within the scope of protection of the present invention.