CN104933709B - Random walk CT lung tissue image automatic segmentation methods based on prior information - Google Patents

Random walk CT lung tissue image automatic segmentation methods based on prior information Download PDF

Info

Publication number
CN104933709B
CN104933709B CN201510309811.4A CN201510309811A CN104933709B CN 104933709 B CN104933709 B CN 104933709B CN 201510309811 A CN201510309811 A CN 201510309811A CN 104933709 B CN104933709 B CN 104933709B
Authority
CN
China
Prior art keywords
image
seed point
random walk
formula
lung tissue
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201510309811.4A
Other languages
Chinese (zh)
Other versions
CN104933709A (en
Inventor
石争浩
司春娇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian University of Technology
Original Assignee
Xian University of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xian University of Technology filed Critical Xian University of Technology
Priority to CN201510309811.4A priority Critical patent/CN104933709B/en
Publication of CN104933709A publication Critical patent/CN104933709A/en
Application granted granted Critical
Publication of CN104933709B publication Critical patent/CN104933709B/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/11Region-based segmentation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30004Biomedical image processing
    • G06T2207/30061Lung

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Processing (AREA)

Abstract

The invention discloses a kind of random walk CT lung tissue image automatic segmentation methods based on prior information, include the following steps:Original CT image is inputted, the CT images of input are carried out two into wavelet transform, obtain image M;Chest area extraction is carried out using the super-pixel method based on entropy rate to image M, image T is obtained after image M extractions chest area processing;The prior information guided using anatomical knowledge obtains seed point number and the position of image T, and carrying out random walk by the seed point of acquisition divides to obtain pulmonary parenchyma initial profile;The reparation for missing cut zone in part is carried out to obtained pulmonary parenchyma initial profile by curvature correction algorithm.Solves the technical issues of being affected by seed point number and position using the segmentation lung tissue image caused by man-machine interactively selected seed point in existing random walk technology.

Description

Random walk CT lung tissue image automatic segmentation methods based on prior information
Technical field
The invention belongs to technical field of image processing, are related to a kind of random walk CT lung tissue images based on prior information Automatic division method.
Background technology
Segmentation of lung parenchyma is the key technology that must be primarily carried out using Computer aided decision system detectio pulmonary disease Step directly affects subsequent detection and analysis result.Studies have shown that the pulmonary disease missing inspection of 5%-7% with lung tissue not Accurate Segmentation is related.Computed tomography (Computer Tomography, CT) is that current Clinical Computer pulmonary disease is auxiliary One of technology relatively conventional in detection is helped, however, since the diversity of biological tissue and the inherence of imaging device are uncertain Property, there is certain fuzziness in the sectioning image that CT transverse scans obtain so that the automatic cutting techniques of computer face very big Problem.Therefore, the automatic accurate segmentation of CT lung tissues image is assisted for improving the computer pulmonary disease based on CT scan The performance of diagnosis and therapy system is of great significance.
In recent years, domestic and foreign scholars have done numerous studies to lung tissue image segmentation.Compare classical, is exactly according to image Gray scale is split, and lung tissue image segmentation can be well realized in conventional threshold values method, at the same time, but fails effectively to remove The regions such as the outer tracheae/bronchus of trunk.Using the region growth method of man-machine interactively selected seed point and growing strategy, Neng Gouzheng Really segmentation CT lung tissue images, and can be very good to retain image diffusivity region, but the method is also easy to produce sky at noise Hole and over-segmentation.Lung tissue image partition method based on movable contour model, by closed curve energy minimization, dynamically to Object boundary is approached to complete to divide, although can preferably retain lung tissue image detail, when lung tissue boundary it is weaker or When person has crack, it is possible to missing edges occur, and it is relatively high to initial profile status requirement, it is longer to calculate the time.This A bit according to the single algorithm of gradation of image, region and marginal information, it is difficult to it is automatic high to meet a large amount of CT slices in clinical application The segmentation demand of effect.It comes into being in conjunction with the integrated dividing method of a variety of single dividing method advantages, is based on Threshold segmentation and number It learns the method that form rolling ball method is combined and uses Mathematical Morphology rolling ball method after threshold method obtains an initial lung areas To correct the false edges formed for the first time.
Recent years inputs seed node, in image essence according to the random walk method that figure optimization is realized by interactive mode Performance is excellent in terms of true segmentation, gradually causes people's great interest.Although this method has well in terms of image Accurate Segmentation Performance, but it but needs man-machine interactively to select initial seed point, and it is sensitive with malposition to the number of seed point.Test table It is bright, though seed point number and the position in segmentation have the variation of very little as a result if having very big difference.
Invention content
The object of the present invention is to provide a kind of random walk CT lung tissue image automatic segmentation methods based on prior information, The dividing method solves in existing random walk technology using the segmentation lung tissue figure caused by man-machine interactively selected seed point The technical issues of as being affected by seed point number and position.
The technical solution adopted in the present invention is random walk CT lung tissues Image Automatic Segmentation side based on prior information Method includes the following steps:
Step 1, original CT image is inputted, the CT images of input are carried out two into wavelet transform, obtain image M;
Step 2, chest area extraction, figure are carried out using the super-pixel method based on entropy rate to the image M obtained in step 1 As obtaining image T after M extractions chest area processing;
Step 3, the prior information guided using anatomical knowledge, the seed point number of image T and position in obtaining step 2 It sets, carrying out random walk by the seed point of acquisition divides to obtain the initial profile of pulmonary parenchyma part;
Step 4:Part accidentally cut zone is carried out to the pulmonary parenchyma initial profile obtained in step 3 by curvature correction algorithm Reparation.
The features of the present invention also characterized in that
The detailed process of step 1 is:
Step 1.1, under two-dimensional case, the wavelet function for choosing three different spaces directions carries out two into discrete wavelet transformer It changes, respectively the wavelet function ψ of horizontal direction1The wavelet function ψ of (x, y), vertical direction2The small echo of (x, y), diagonal Function ψ3The wavelet transformation of (x, y), CT image f (x, y) is:
Wherein, n=1,2,3;J indicates scale;ψnIn (x, y), using the first derivative of cubic B-spline as small echo letter Number, then CT images f (x, y) is expressed as
Wherein ForDual Wavelet function,ForAntithesis scale Function, t1、t2Indicate different time domain, u1、u2Indicate different frequency domains;
After CT image f (x, y) are by two-dimensional binary discrete wavelet fast decoupled, a series of son of different resolutions is obtained Image, expression are as follows:
Wherein, 0≤j≤J=log2N, hkAnd gkIt is the low-pass filter and high-pass filter of wavelet function, l, k table respectively Show that the different sampling processes of low-pass filter and high-pass filter, m, n indicate the row and column of CT images, aj+1For low frequency subgraph, dj + 1,1For the high frequency subgraph of horizontal direction, dJ+1,2The high frequency subgraph of 2 vertical direction, dJ+1,3For the high frequency subgraph of diagonal;
Step 1.2, j=2 is enabled, to low frequency subgraph aj+1Histogram equalization is carried out to improve the grey-scale contrast of CT images; To the high frequency subgraph d of horizontal directionJ+1,1, vertical direction high frequency subgraph dJ+1,2, diagonal high frequency subgraph dJ+1,3It carries out Gaussian filtering is smoothly to remove the noise of CT images;
Step 1.3, to by step 1.2 treated CT images according to iteration two channels filter described in step 1.1 Inverse process image reconstruction is carried out using formula (4), obtain image M after image reconstruction, formula (4) expression formula is as follows:
The detailed process of wherein step 2 is:
Step 2.1, figure G is enabled1=(V1, E1) on the object function of entropy rate super-pixel be:
maxH(A)+λB(A) (5);
Wherein H (A) indicates figure G1=(V1, E1) on random walk entropy rate, B (A) indicate balance term, A is the side of selection Collection, E1For total side collection;λ >=0 is coefficient of balance;
Step 2.2, G will be schemed in step 2.11=(V1, E1) in each vertex regard an individual sub-collective drawing as, each Different sub-collective drawings number is marked in the sub-collective drawing, initializes the A=φ;Pass through Europe using the object function in formula (5) Family name's distance calculates the figure G1=(V1, E1) on E1In each side functional value;
Step 2.3, the maximum side of functional value is selected to be added in A from the calculated functional value of step 2.2, and from E1In subtract The maximum side of this functional value is gone to, the sub-collective drawing at the maximum side both ends of the obtained functional value is merged into one;
Step 2.4, it is new pair with maximum functional value while remaining each after during A is added of functional value in step 2.3 As step 2.3 operation being repeated, until E1For empty or image M sub-collective drawing number NAEqual to setting value 2, each image M Sub-collective drawing be exactly a super-pixel, cluster is generated as 2 mask according to obtained super-pixel result;
Step 2.5, the mask that the cluster that step 2.4 obtains is 2 is multiplied with image M and extracts chest area, extracted The chest area image T gone out.
Wherein step 3 the specific steps are:
Step 3.1, if the lung tissue Some seeds point number of random walk is 2 in image T, non-lung Some seeds point Number is 1;
Step 3.2, it uses in step 2 and the mask clustered as 2 is generated to image T based on the super-pixel method of entropy rate, it will The mask that the cluster generated in this step is 2 multiplies with image T-phase, tentatively obtains the background except lung tissue and lung tissue The seed point in region, lung tissue region is respectively left and right lung tissue areas outside profile middle 10 pixel-shifts inwards The position of amount, the seed point of background area is at the center of background area except lung tissue;
Step 3.3, the seed point location that the seed point number and step 3.2 obtained according to step 3.1 obtains carries out random Migration is divided, and the initial profile of pulmonary parenchyma part is obtained.
The detailed process of wherein step 3.3 is:
Step A:By image T regard as in graph theory without phase weighted graph G2=(V2, E2, W), pixel is that eight neighborhood closes in figure It is the node on side,
Wherein, V2It is point set, E2It is side collection;
W is the weight function on side:wij=exp [- β (gi-gj)2] (9);
Wherein gi, gjIntensity value at respectively pixel i and pixel j, β are free parameters;
Step B, the seed point location that the seed point number and step 3.2 obtained according to step 3.1 obtains, passes through minimum The energy function of formula (10) calculates the transition probability that non-seed point in random walk reaches each seed point, by comparing non-kind Son point reaches the numerical value of the transition probability of each seed point, and non-seed point is divided into the larger a kind of seed point of transition probability In, detailed process is as follows,
The energy function that minimum formula is obtained in conjunction with the formula (9) in step A is
Wherein, L is the joint Laplacian Matrix defined in mapping graph, as i=j, Lij=di, as v in ViAnd vjFor phase When neighbors, Lij=-wij, work as viAnd vjWhen not being adjacent node, Lij=0, wherein diIt is the degree of node i;
Step C, the seed point location that the seed point number and step 3.2 obtained according to step 3.1 obtains will be in image T Node be divided into two set:Seed point set VMWith non-seed point set VU, on this basis, Laplce's square in step B Battle array L is decomposed into diagonal hypermatrix form:
Formula (11) is brought into formula (10) and is obtained,
Wherein, XMAnd XUThe probability value of seed point and non-seed node is corresponded to respectively, and B is intermediate variable;
To function D [X in formula (12)U] differentiate,
LUXU=-BTXM(13);
It enablesIndicate non-seed point ViThe probability value of the seed point of classification S is reached for the first time, and defines seed point set Q (Vi)=s,0 < s≤K, K is the number of seed point, therefore defines K dimension groups M and be:
Probability value size, s ∈ S, by formula (16) to all seed points are calculated to each seed point s by formula (15) Probability value size is calculated, formula (15) and formula (16) are as follows:
LUX=-BM (16);
Wherein, X is indicatedThe column vector of composition, x are transition probabilities, and M indicates msThe column vector of composition, it is known that:
Simultaneous formula (14), formula (15), formula (16), formula (17), solve
In two dimensional image,The pixel of the i-th row jth row is indicated to the transition probability of seed point S, to by each non-kind Son point is integrated into the seed point class belonging to transition probability maximum value, realizes lung tissue image segmentation.
The detailed process of wherein step 4 is:
Step 4.1, the random walk dividing wheel profile that scanning step 3 obtains first is searched using the row scanning for being divided into L L is set as 3 pixels by the crosspoint of rope scan line and random walk segmentation contour, and the profile crosspoint of scanning is divided into 3 classes:It rises Initial pointIntermediate pointDestination nodeDue to starting pointAnd destination nodeIt is located at outside or interior On the lung contour line of side, and intermediate pointTypically in vertical diaphragm or comprising the recess of Lung neoplasm, therefore, starting point is chosenWith destination nodeIt is repaired;
Step 4.2, the curvature k of scan line and the intersection of random walk segmentation contour in step 4.1 is calculatedi, calculate Formula is:
Wherein, (xi,yi) be the i-th row scan line and random walk dividing wheel profile crosspoint, (xi-1,yi-1) and (xi+1,yi+1) crosspoint of the (i-1)-th row and i+1 horizontal scanning line and random walk segmentation contour is indicated respectively;
Step 4.3, with the curvature k of calculated scan line in step 4.2 and the intersection of random walk segmentation contouri For foundation, random walk segmentation contour is repaired, eliminates the initial profile and chest for obtaining pulmonary parenchyma part in the step 3 The connected tubercle of film and Pulmonary Vascular over-segmentation situation.
The detailed process wherein repaired to random walk segmentation contour in step 4.3 is:
With the curvature k of calculated scan line in step 4.2 and the intersection of random walk dividing wheel profileiFor according to According to, setting priori threshold value is 0.2, when adjacent horizontal scanning line upper tracer and random walk dividing wheel profile two crosspoints or Curvature difference between same horizontal scanning line upper tracer and two crosspoints of random walk dividing wheel profile is more than priori threshold value When 0.2, linear interpolation, the adjacent beginning and end of connection random walk segmentation contour defect boundary, mending course are used It completes.
The invention has the advantages that the number of seed point is precipitated in the present invention by anatomy priori statistical first With position, segmentation will be significantly affected by overcoming the seed point number that man-machine interactively in the prior art is chosen and being varied slightly with position As a result the problem of, and Algorithms T-cbmplexity is effectively reduced, improve the accuracy of segmentation.In addition, using two into discrete small Pretreatment of the wave conversion as random walk has not only repaired the defects of image segmentation region in conjunction with curvature correct algorithm, and And the diffusivity region in image is contained, improve segmentation robustness.
Description of the drawings
Fig. 1 is the flow chart of the random walk CT lung tissue image automatic segmentation methods the present invention is based on prior information;
Fig. 2 is original CT image;
Fig. 3 is the figure converted into discrete wavelet to original CT image shown in Fig. 2 using proposed by the present invention two As M;
Fig. 4 is that treated for two into discrete picture M extraction chest areas shown in Fig. 3 image T;
Fig. 5 is using the random walk proposed by the present invention guided based on anatomy priori to chest figure shown in Fig. 4 As the auto divide image of T;
Fig. 6 is to be carried out using curvature correct algorithm based on prior information proposed by the present invention to dividing image shown in Fig. 5 The image that boundary is corrected;
Fig. 7 is using classical threshold method to the images that carry out that treated of CT images shown in Fig. 2;
Fig. 8 is to carry out treated image to CT images shown in Fig. 2 using existing region growing methods;
Fig. 9 is to carry out treated image to CT images shown in Fig. 2 using existing active contour method;
Figure 10 is using traditional random walk method to the images that carry out that treated of CT images shown in Fig. 2;
Figure 11 is to carry out treated image to CT images shown in Fig. 2 using existing threshold value morphological method;
Figure 12 is the result that experienced dept. of radiology expert is split CT images shown in Fig. 2 manually.
Specific implementation mode
The following describes the present invention in detail with reference to the accompanying drawings and specific embodiments.
The present invention is based on the random walk CT lung tissue image automatic segmentation methods of prior information, as shown in Figure 1, include with Lower step:
Step 1, the CT images of input are carried out two into wavelet transform, obtained by input original CT image (as shown in Figure 2) To image M;
Step 1.1:Under two-dimensional case, the wavelet function for choosing three different spaces directions carries out two into discrete wavelet transformer It changes, respectively the wavelet function ψ of horizontal direction1The wavelet function ψ of (x, y), vertical direction2The small echo of (x, y), diagonal Function ψ3The wavelet transformation of (x, y), CT image f (x, y) is:
Wherein, n=1,2,3;J indicates scale;ψnIn (x, y), using the first derivative of cubic B-spline as small echo letter Number, then CT images f (x, y) is expressed as
In formula Wherein,ForDual Wavelet function,ForAntithesis Scaling function, t1、t2Indicate different time domain, u1、u2Indicate different frequency domains;
After CT image f (x, y) are by two-dimensional binary discrete wavelet fast decoupled, a series of son of different resolutions is obtained Image, expression are as follows:
Wherein, 0≤j≤J=log2N, J and j indicate that scale, J are out to out, and j is to refer to specific single scale, hkWith gkIt is the low-pass filter and high-pass filter of wavelet function respectively, l, k indicate the difference of low-pass filter and high-pass filter Sampling process, m, n indicate the row and column of CT images, what CT images were obtained after two-dimentional (M × N) discrete dyadic wavelet fast decoupled Four sub-band images are respectively:aj+1For low frequency subgraph, dJ+1,1For the high frequency subgraph of horizontal direction, dJ+1,2For the height of vertical direction Frequency subgraph, dJ+1,3For the high frequency subgraph of diagonal, low frequency subgraph expresses the approximate information of image, and high frequency subgraph is image The expression of detailed information in different directions.When carrying out dyadic wavelet decomposition, noise is broken down into high frequency subgraph;
Step 1.2, select two scales for best scale layer, to low frequency subgraph aj+1Histogram equalization is carried out to improve CT The grey-scale contrast of image reinforces the edge details of image;To the high frequency subgraph d of horizontal directionJ+1,1, vertical direction high frequency Scheme dJ+1,2, diagonal high frequency subgraph dJ+1,3Gaussian filtering is carried out smoothly to remove the noise that the CT is sliced F images; It in experiment, selects two scales for best scale layer, i.e. j=2, can not only obtain lung tissue boundary characteristic, but also can be effective Ground inhibits noise, and experiment shows if more than two scale, since edge smoothing can lead to over-segmentation;
Step 1.3, to by step 1.2 treated CT images according in step 1.1 iteration two channels filter it is inverse (the positive process of iteration two channels filter is after CT images are decomposed by two-dimentional (M × N) discrete dyadic wavelet, to obtain 4 to process The process of sub-band images) use formula (4) to carry out image reconstruction, two are obtained into discrete wavelet to the original CT after image reconstruction The image M that slice F is converted, formula (4) expression formula are as follows:
Handle independent to each scale is that it is possible to using two into the advantages of wavelet transform.
Two images converted into discrete wavelet to original CT slice F are obtained after carrying out image reconstruction by formula (4) M, as shown in Figure 3.
Step 2, chest area extraction, figure are carried out using the super-pixel method based on closely related rate to the image M obtained in step 1 As obtaining image T after M extractions chest area processing, super-pixel refer to have many characteristics, such as similar grain, color, brightness adjacent picture The region unit that element is constituted.Based on the superpixel segmentation method of entropy rate, it can not only make that the segmentation result of super-pixel is compact and region is equal It is even, and regular shape and size it is similar.Detailed process is:
Step 2.1, figure G is enabled1=(V1, E1) on the object function of entropy rate super-pixel be:
maxH(A)+λB(A) (5);
Wherein H (A) indicates figure G1=(V1, E1) on random walk entropy rate, B (A) indicate balance term, A is the side of selection Collection, E1For total side collection;λ >=0 is coefficient of balance;
Scheme G1=(V1, E1) on the definition of entropy rate H (A) of random walk be:
H (A)=ΣiμiΣjpi,j(A)log(pi,j(A)) (6);
Wherein, i and j indicates the vertex of figure, and μ is the Stationary Distribution of random process,wiIt is to be connected with vertex i The sum of side right weight,It is normalization constant, | V | indicate the number of vertices of figure, pi,j(A) it is random walk on figure Transition probability is defined as:
When compact and homogeneous area side being only added into A, it can just make entropy rate increase most fast, wI, jRepresent vertex i and top Weights between point j on side.
Balance term B (A) is defined as follows:
Wherein NAIndicate that the number of connected subgraph in figure G '=(V, A), G '=(V, A) indicate G1=(V1, E1) in connection Figure, V is the vertex of connected graph, and A is the side in UNICOM domain, ZAIndicate the distribution of connected subgraph cluster,Indicate i-th A subgraph number of vertex proportion, | Si| indicate that the number of vertices of i-th of subgraph, balance term B (A) can make subgraph cluster more equal It is even.
Step 2.2, G will be schemed in step 2.11=(V1, E1) in each vertex regard an individual sub-collective drawing as, each Different sub-collective drawings number is marked in the sub-collective drawing, initializes the A=φ;Pass through Europe using the object function in formula (5) Family name's distance calculates figure G1=(V1, E1) on E1In each side functional value;
Step 2.3, the maximum side of functional value is selected to be added in A from the calculated functional value of step 2.2, and from E1In subtract The maximum side of this functional value is gone to, the sub-collective drawing at the maximum side both ends of the obtained functional value is merged into one;
Step 2.4, it is new pair with maximum functional value while remaining each after during A is added of functional value in step 2.3 As (after the maximum side of functional value that first calculated goes out is added in A, then from remaining functional value to select functional value maximum Side recycles successively), step 2.3 operation is repeated, until E1For empty or described image M sub-collective drawing number NAEqual to setting The sub-collective drawing of value 2, each described image M is exactly a super-pixel, and generate cluster according to obtained super-pixel result covers for 2 Film;
Step 2.5, the mask that the cluster that step 2.4 obtains is 2 is multiplied with described image M and extracts chest area, obtained The chest area image T extracted.The image T of chest area extraction is carried out as shown in figure 4, extraction chest using entropy rate super-pixel Region, which not only meets test requirements document, can also reduce calculating cost.
Step 3, the prior information guided using anatomical knowledge, the number of image T seed points and position in obtaining step 2 It sets, carrying out random walk to seed point by the seed point location of acquisition divides to obtain the initial profile of pulmonary parenchyma part.Specifically Steps are as follows:
Step 3.1, it is dug according to chest solution and learns structure knowledge, there is the area except pulmo and lung tissue in human chest region Domain is constituted, and therefore, in Random Walk Algorithm, the lung tissue Some seeds point number of random walk in image T is 2, non-lung Some seeds point number is 1;
Step 3.2, it uses in step 2 and the mask clustered as 2 is generated to image T based on the super-pixel method of closely related rate, it will The mask that the cluster generated in this step is 2 multiplies with image T-phase, tentatively obtains the background except lung tissue and lung tissue The seed point in region, lung tissue region is respectively left lung tissue areas outside profile middle 10 pixel-shift amounts inwards The position of position, right lung tissue regions lateral profile middle 10 pixel-shift amounts inwards, background area except lung tissue Seed point be background area center at;
Step 3.3, the seed point location that the seed point number and step 3.2 obtained according to step 3.1 obtains carries out random Migration is divided, and the initial profile of pulmonary parenchyma part is obtained.The detailed process of random walk segmentation is as follows:
Step A:By image T regard as in graph theory without phase weighted graph G2=(V2, E2, W), pixel is that eight neighborhood closes in figure It is the node on side,
Wherein, V2It is point set, E2It is side collection;
W is the weight function on side:wij=exp [- β (gi-gj)2] (9);
Wherein gi, gjIntensity value at respectively pixel i and pixel j, β are free parameters;
Step B, the seed point location that the seed point number and step 3.2 obtained according to step 3.1 obtains, passes through minimum The energy function of formula (10) calculates the transition probability that non-seed point in random walk reaches each seed point, by comparing non-kind Son point reaches the numerical value of the transition probability of each seed point, and non-seed point is divided into the larger a kind of seed point of transition probability In, detailed process is as follows,
The energy function that minimum formula is obtained in conjunction with the formula (9) in step A is
Wherein, L is the joint Laplacian Matrix defined in mapping graph, as i=j, Lij=di, as v in ViAnd vjFor phase When neighbors, Lij=-wij, work as viAnd vjWhen not being adjacent node, Lij=0, wherein diIt is the degree of node i;
Step C, the seed point location that the seed point number and step 3.2 obtained according to step 3.1 obtains, by described image Node in T is divided into two set:Seed point set VMWith non-seed point set VU, on this basis, the La Pula in step B This matrix L is decomposed into diagonal hypermatrix form:
Formula (11) is brought into formula (10) and is obtained,
Wherein, XMAnd XUThe probability value of seed point and non-seed node is corresponded to respectively, and B is intermediate variable;
To function D [X in formula (12)U] differentiate,
LUXU=-BTXM(13);
It enablesIndicate non-seed point ViThe probability value of the seed point of classification S is reached for the first time, and defines seed point set Q (Vi)=s,0 < s≤K, K is the number of seed point, therefore it is as follows to define K dimension groups M:
Probability value size, s ∈ S, by formula (16) to all seed points are calculated to each seed point s by formula (15) Probability value size is calculated, formula (15) and formula (16) are as follows:
LUX=-BM (16);
Wherein, X is indicatedThe column vector of composition, x are transition probabilities, and M indicates msThe column vector of composition, it is known that:
Simultaneous formula (14), formula (15), formula (16), formula (17), solve
In two dimensional image,The pixel of the i-th row jth row is indicated to the transition probability of seed point S, to by each non-kind Son point is integrated into the seed point class belonging to transition probability maximum value, realizes lung tissue image segmentation.Utilize above-mentioned dissection Xue Zhi The random walk for knowing guiding is split that the results are shown in Figure 5 to image after being pre-processed shown in Fig. 4, and Random Walk Algorithm can return Become the solution for seeking sparse, symmetric positive definite system a linear equation.
Step 4:Part is carried out by curvature correction algorithm to the pulmonary parenchyma part initial profile obtained in step 3 accidentally to divide The reparation in region, specific step of repairing are:
Step 4.1, the random walk dividing wheel profile that scanning step 3 obtains first calculates the time to reduce, improves and divide The robustness for cutting result scans to search for the crosspoint of scan line and random walk segmentation contour, experiment using the row for being divided into L In, L is set as 3 pixels, the profile crosspoint of scanning is divided into 3 classes:Starting pointIntermediate pointDestination node Due to starting pointAnd destination nodeIt is located on the lung contour line of outside or inside, and intermediate pointTypically In vertical diaphragm or comprising the recess of Lung neoplasm, therefore, starting point is chosenWith destination nodeIt is repaired;
Step 4.2, the curvature k of scan line and the intersection of random walk segmentation contour in step 4.1 is calculatedi, calculate Formula is:
Wherein, (xi,yi) be the i-th row scan line and random walk dividing wheel profile crosspoint, (xi-1,yi-1) and (xi+1,yi+1) crosspoint of the (i-1)-th row and i+1 horizontal scanning line and random walk segmentation contour is indicated respectively;
Step 4.3, with the curvature k of calculated scan line in step 4.2 and the intersection of random walk segmentation contouri For foundation, random walk segmentation contour is repaired, eliminates the initial profile and chest for obtaining pulmonary parenchyma part in the step 3 The connected tubercle of film and Pulmonary Vascular over-segmentation situation.The larger point of curvature value is usually happened at some small pulmonary parenchyma profile and border waves At the normal lung tissue such as dynamic and initial profile bottom, it is therefore desirable to go to reject them using different rules.Experiment hair It is existing, two crosspoints of adjacent horizontal scanning line upper tracer and random walk segmentation contour or same horizontal scanning line upper tracer with When curvature difference between two crosspoints of random walk segmentation contour is more than a priori threshold value, typically at lung tissue edge At the pleura that profile is connected with tubercle or at Pulmonary Vascular, therefore it is as follows to repair principle:
With the curvature k of calculated scan line in step 4.2 and the intersection of random walk dividing wheel profileiFor according to According to, setting priori threshold value is 0.2, when adjacent horizontal scanning line upper tracer and random walk dividing wheel profile two crosspoints or Curvature difference between same horizontal scanning line upper tracer and two crosspoints of random walk dividing wheel profile is more than priori threshold value When 0.2, linear interpolation, the adjacent beginning and end of connection random walk segmentation contour defect boundary, mending course are used It completes.
There is no certain changing rules for the curvature of scan line and random walk segmentation contour intersection, but have been found that, When scan line and the curvature difference in two crosspoints of random walk segmentation contour are more than priori threshold value 0.2, random walk segmentation Profile defects are typically at the pleura being connected with tubercle in lung tissue edge or at Pulmonary Vascular.Fig. 6 gives profile repairing Afterwards as a result, as seen from the figure, right lung profile is than more complete.
With the hand labeled results contrast of dept. of radiology expert shown in Figure 12, threshold method shown in Fig. 7 causes to be excluded in lung For main tracheae except portion region there may be the danger that cancer embolus occurs, growth method segmenting edge in region shown in Fig. 8 is relatively rough, And due to the limitation of growing strategy, as a result generate over-segmentation phenomenon;After the evolution segmentation shown in Fig. 9 using active contour, diffuse Property region lose;Traditional random walk processing shown in Fig. 10, can obviously generate inaccurate segmentation result;Shown in Figure 11 Classical threshold value morphological method, boundary is more smooth, and segmentation result is more accurate.As shown in fig. 6, using the present invention is based on elder generations The random walk CT lung tissue image automatic segmentation methods segmentation of information is tested, compared with other five methods, segmentation result is not only The pleura region being connected with tubercle is contained, also contains diffusivity region well.
Table 1 gives the case where above-mentioned several method is split CT slices F, is mainly tied to segmentation using precision index Fruit is analyzed.Accuracy is to weigh the accuracy of image segmentation, emphasizes the segmentation between segmentation image and graphics standard segmentation Difference.Accuracy is bigger, indicates the image after segmentation closer to Standard Segmentation.Accuracy is defined as:
Wherein NTPIt is the percentage in the lung tissue region correctly divided;NTNIt is the background percentage correctly divided;NFPIt is wrong The accidentally target area percentage of segmentation;NFNIt is the background percentage of erroneous segmentation.
The accuracy that 1 distinct methods of table divide original slice images compares
It can be seen that from the data in table 1 using the image after automatic Segmentation of the present invention, accuracy value is higher than other Method as a result, significantly improving the segmentation performance of image.

Claims (5)

1. random walk CT lung tissue image automatic segmentation methods based on prior information, it is characterised in that:Include the following steps:
Step 1, original CT image is inputted, the CT images of input are carried out two into wavelet transform, obtain image M;
Step 2, chest area extraction, image M are carried out using the super-pixel method based on entropy rate to the image M obtained in step 1 Image T is obtained after the processing of extraction chest area;
Step 2 detailed process is:
Step 2.1, figure G is enabled1=(V1, E1) on the object function of entropy rate super-pixel be:
maxH(A)+λB(A) (5);
Wherein H (A) indicates figure G1=(V1, E1) on random walk entropy rate, B (A) indicate balance term, A is the side collection of selection, E1For Total side collection;λ >=0 is coefficient of balance;
Step 2.2, G will be schemed in step 2.11=(V1, E1) in each vertex regard an individual sub-collective drawing as, it is each described Different sub-collective drawings number is marked in sub-collective drawing, initializes the A=φ;Using the object function in formula (5) by Euclidean away from Figure G described from calculating1=(V1, E1) on E1In each side functional value;
Step 2.3, the maximum side of functional value is selected to be added in A from the calculated functional value of step 2.2, and from E1In subtract this The maximum side of functional value, one is merged by the sub-collective drawing at the maximum side both ends of the obtained functional value;
Step 2.4, with maximum functional value while remaining each after during A is added of functional value in step 2.3 for new object, weight Step 2.3 operation is executed again, until E1For empty or described image M sub-collective drawing number NAEqual to setting value 2, each figure As the sub-collective drawing of M is exactly a super-pixel, cluster is generated as 2 mask according to obtained super-pixel result;
Step 2.5, the mask that the cluster that step 2.4 obtains is 2 is multiplied with described image M and extracts chest area, extracted The chest area image T gone out;
Step 3, the prior information guided using anatomical knowledge, the seed point number of image T and position in obtaining step 2, is led to The seed point progress random walk obtained is crossed to divide to obtain pulmonary parenchyma initial profile;
The detailed process of the step 3 is:
Step 3.1, if the lung tissue Some seeds point number of random walk is 2 in described image T, non-lung Some seeds point Number is 1;
Step 3.2, cluster is generated to described image T based on the super-pixel method of entropy rate covered for 2 described in step 2 is used The mask that the cluster generated in this step is 2 is multiplied with described image T-phase, tentatively obtains lung tissue and lung tissue by film Except background area, the seed point in lung tissue region is respectively left and right lung tissue areas outside profile middle inwards 10 The position of a pixel-shift amount, the seed point of background area is at the center of background area except lung tissue;
Step 3.3, the seed point location that the seed point number and step 3.2 obtained according to step 3.1 obtains carries out random walk Segmentation, obtains the initial profile of pulmonary parenchyma part;
Step 4:Part is carried out to the pulmonary parenchyma initial profile obtained in step 3 by curvature correction algorithm and misses repairing for cut zone It is multiple.
2. random walk CT lung tissue image automatic segmentation methods based on prior information according to claim 1, special Sign is:The detailed process of the step 1 is:
Step 1.1:Under two-dimensional case, the wavelet function for choosing three different spaces directions carries out two into wavelet transform, The respectively wavelet function ψ of horizontal direction1The wavelet function ψ of (x, y), vertical direction2The wavelet function of (x, y), diagonal ψ3The wavelet transformation of (x, y), the CT images f (x, y) is:
Wherein, n=1,2,3;J indicates scale;ψnIn (x, y), using the first derivative of cubic B-spline as wavelet function, then institute CT image f (x, y) are stated to be expressed as
In formula,Wherein, For { ψn(x, y) }1≤n≤3Dual Wavelet function,ForDual scaling functions, t1、t2It indicates not Same time domain, u1、u2Indicate different frequency domains;
After the CT images f (x, y) is by two-dimensional binary discrete wavelet fast decoupled,
A series of subgraph of different resolutions is obtained, expression is as follows:
Wherein, 0≤j≤J=log2N, hkAnd gkIt is the low-pass filter and high-pass filter of wavelet function respectively, l, k indicate low The different sampling processes of bandpass filter and high-pass filter, m, n indicate the row and column of CT images, aj+1For low frequency subgraph, dJ+1,1For The high frequency subgraph of horizontal direction, dJ+1,2For the high frequency subgraph of vertical direction, dJ+1,3For the high frequency subgraph of diagonal;
Step 1.2, j=2 is enabled, to the low frequency subgraph aj+1Histogram equalization is carried out to improve the gray scale pair of the CT images Degree of ratio;To the high frequency subgraph d of the horizontal directionJ+1,1, vertical direction high frequency subgraph dJ+1,2, diagonal high frequency Scheme dJ+1,3Gaussian filtering is carried out smoothly to remove the noise that the CT is sliced F images;
Step 1.3, to the low-pass filter h by step 1.2 treated CT images according to wavelet function described in step 1.1k With high-pass filter gkInverse process image reconstruction is carried out using formula (4), image M, formula (4) expression are obtained after image reconstruction Formula is as follows:
3. random walk CT lung tissue image automatic segmentation methods based on prior information according to claim 1, special Sign is:The detailed process that the step 3.3 carries out random walk segmentation is:
Step A:By described image T regard as in graph theory without phase weighted graph G2=(V2, E2, W), pixel is that eight neighborhood closes in figure It is the node on side,
Wherein, V2It is point set, E2It is side collection;
W is the weight function on side:wij=exp [- β (gi-gj)2]
(9);
Wherein gi, gjIntensity value at respectively pixel i and pixel j, β are free parameters;
Step B, the seed point location that the seed point number and step 3.2 obtained according to step 3.1 obtains, by minimizing formula (10) energy function calculates the transition probability that non-seed point in random walk reaches each seed point, by comparing non-seed point Non-seed point is divided into the larger a kind of seed point of transition probability by the numerical value for reaching the transition probability of each seed point, tool Body process is as follows,
The energy function that minimum formula is obtained in conjunction with the formula (9) in step A is
Wherein, L is the joint Laplacian Matrix defined in mapping graph, as i=j, Lij=di, as v in ViAnd vjFor adjacent segments When point, Lij=-wij, work as viAnd vjWhen not being adjacent node, Lij=0, wherein diIt is the degree of node i;
Step C, the seed point location that the seed point number and step 3.2 obtained according to step 3.1 obtains will be in described image T Node be divided into two set:Seed point set VMWith non-seed point set VU, on this basis, Laplce's square in step B Battle array L is decomposed into diagonal hypermatrix form:
Formula (11) is brought into formula (10) and is obtained,
Wherein, XMAnd XUThe probability value of seed point and non-seed node is corresponded to respectively, and B is intermediate variable;
To function D [X in formula (12)U] differentiate,
LUXU=-BTXM(13);
It enablesIndicate non-seed point ViThe probability value of the seed point of classification S is reached for the first time, and defines seed point set Q (Vi)= S,0 < s≤K, K is the number of seed point, therefore defines K dimension groups M and be:
Probability value size is calculated to each seed point s by formula (15), s ∈ S calculate all seed points by formula (16) Probability value size, formula (15) and formula (16) are as follows:
LUX=-BM (16);
Wherein, X is indicatedThe column vector of composition, x are transition probabilities, and M indicates msThe column vector of composition, it is known that:
Simultaneous formula (14), formula (15), formula (16), formula (17), solve
In two dimensional image,The pixel of the i-th row jth row is indicated to the transition probability of seed point S, thus will each non-seed point It is integrated into the seed point class belonging to transition probability maximum value, realizes lung tissue image segmentation.
4. random walk CT lung tissue image automatic segmentation methods based on prior information according to claim 1, special Sign is:The detailed process of the step 4 is:
Step 4.1, the random walk dividing wheel profile that scanning step 3 obtains first scans to search for sweeping using the row for being divided into L L is set as 3 pixels by the crosspoint for retouching line and random walk segmentation contour, and the profile crosspoint of scanning is divided into 3 classes:Starting point Pi first, intermediate point Pi last, destination node Pi middle(j), due to starting point Pi firstWith destination node Pi middle(j)Be located at outside or On the lung contour line of inside, and intermediate point Pi lastTypically in vertical diaphragm or comprising the recess of Lung neoplasm, therefore, starting point is chosen Pi firstWith destination node Pi middle(j)It is repaired;
Step 4.2, the curvature k of scan line and the intersection of random walk segmentation contour in step 4.1 is calculatedi, calculation formula For:
Wherein, (xi,yi) be the i-th row scan line and random walk dividing wheel profile crosspoint, (xi-1,yi-1) and (xi+1, yi+1) crosspoint of the (i-1)-th row and i+1 horizontal scanning line and random walk segmentation contour is indicated respectively;
Step 4.3, with the curvature k of calculated scan line in step 4.2 and the intersection of random walk segmentation contouriFor according to According to being repaired to random walk segmentation contour, eliminate the initial profile and pleura phase for obtaining pulmonary parenchyma part in the step 3 Connection section and Pulmonary Vascular over-segmentation situation.
5. random walk CT lung tissue image automatic segmentation methods based on prior information according to claim 4, special Sign is:The detailed process repaired to random walk segmentation contour in the step 4.3 is:
With the curvature k of calculated scan line in the step 4.2 and the intersection of random walk dividing wheel profileiFor foundation, It is 0.2 that priori threshold value, which is arranged, when two crosspoints of adjacent horizontal scanning line upper tracer and random walk dividing wheel profile or same Curvature difference between one horizontal scanning line upper tracer and two crosspoints of random walk dividing wheel profile is more than priori threshold value 0.2 When, using linear interpolation, the adjacent beginning and end of connection random walk segmentation contour defect boundary, mending course is complete At.
CN201510309811.4A 2015-06-04 2015-06-04 Random walk CT lung tissue image automatic segmentation methods based on prior information Expired - Fee Related CN104933709B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201510309811.4A CN104933709B (en) 2015-06-04 2015-06-04 Random walk CT lung tissue image automatic segmentation methods based on prior information

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201510309811.4A CN104933709B (en) 2015-06-04 2015-06-04 Random walk CT lung tissue image automatic segmentation methods based on prior information

Publications (2)

Publication Number Publication Date
CN104933709A CN104933709A (en) 2015-09-23
CN104933709B true CN104933709B (en) 2018-09-14

Family

ID=54120863

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201510309811.4A Expired - Fee Related CN104933709B (en) 2015-06-04 2015-06-04 Random walk CT lung tissue image automatic segmentation methods based on prior information

Country Status (1)

Country Link
CN (1) CN104933709B (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110503628A (en) * 2019-07-15 2019-11-26 心医国际数字医疗系统(大连)有限公司 The restorative procedure in the U-shaped region of CT image, the judgment method of long short side and pulmonary parenchyma exterior contour

Families Citing this family (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105354846A (en) * 2015-11-05 2016-02-24 沈阳东软医疗系统有限公司 Method and apparatus for segmenting three-dimensional medical image
WO2017092615A1 (en) 2015-11-30 2017-06-08 上海联影医疗科技有限公司 Computer aided diagnosis system and method
CN107203998B (en) * 2016-03-18 2020-04-03 北京大学 Method for carrying out dentition segmentation on cone beam CT image
CN106355023A (en) * 2016-08-31 2017-01-25 北京数字精准医疗科技有限公司 Open quantitative analysis method and system based on medical image
CN108470331B (en) * 2017-02-23 2021-12-21 富士通株式会社 Image processing apparatus, image processing method, and program
CN106971400B (en) * 2017-03-10 2020-11-10 中国航空工业集团公司洛阳电光设备研究所 Method and device for repairing image dividing line
CN106997596B (en) * 2017-04-01 2019-08-20 太原理工大学 A kind of Lung neoplasm dividing method of the LBF movable contour model based on comentropy and joint vector
CN107316289B (en) * 2017-06-08 2020-05-08 华中农业大学 Method for dividing rice ears in field based on deep learning and superpixel division
CN107341812B (en) * 2017-07-04 2019-11-08 太原理工大学 A kind of sequence Lung neoplasm image partition method based on super-pixel and Density Clustering
CN108520525B (en) * 2018-04-12 2021-11-02 重庆理工大学 Spinal cord segmentation method based on convex constraint seed region growth
CN110874170A (en) * 2018-08-31 2020-03-10 杭州海康威视数字技术股份有限公司 Image area correction method, image segmentation method and device
CN109191452B (en) * 2018-09-12 2021-10-08 南京大学 Peritoneal transfer automatic marking method for abdominal cavity CT image based on active learning
CN109741349B (en) * 2019-01-24 2021-12-07 江门市中心医院 Method for segmenting cerebral arterial thrombosis image
CN111145185B (en) * 2019-12-17 2023-12-22 天津市肿瘤医院 Lung substance segmentation method for extracting CT image based on clustering key frame
CN111932506B (en) * 2020-07-22 2023-07-14 四川大学 Method for extracting discontinuous straight line in image
CN116630329B (en) * 2023-07-26 2023-09-29 山东山森数控技术有限公司 Online visual detection method for multi-axis multi-channel numerical control system

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103258324A (en) * 2013-04-02 2013-08-21 西安电子科技大学 Remote sensing image change detection method based on controllable kernel regression and superpixel segmentation
CN103824295A (en) * 2014-03-03 2014-05-28 天津医科大学 Segmentation method of adhesion vascular pulmonary nodules in lung CT (computed tomography) image

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103258324A (en) * 2013-04-02 2013-08-21 西安电子科技大学 Remote sensing image change detection method based on controllable kernel regression and superpixel segmentation
CN103824295A (en) * 2014-03-03 2014-05-28 天津医科大学 Segmentation method of adhesion vascular pulmonary nodules in lung CT (computed tomography) image

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
"一种基于改进随机游走的肺结节分割方法;依玉峰 等;《东北大学学报(自然科学版)》;20120331;第33卷(第3期);第318-322页 *
An Automatic Method for Lung Segmentation in Thin Slice Computed Tomography Based on Random Walks;Zhenghao Shi等;《Journal of Medical imaging and Health informatics》;20150401;第5卷(第2期);第303-308页 *
一种改进的CT胸部图像肺组织分割方法;孟 等;《东北大学学报(自然科学版)》;20090227;第30卷(第2期);第187-190页 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110503628A (en) * 2019-07-15 2019-11-26 心医国际数字医疗系统(大连)有限公司 The restorative procedure in the U-shaped region of CT image, the judgment method of long short side and pulmonary parenchyma exterior contour
CN110503628B (en) * 2019-07-15 2021-06-25 心医国际数字医疗系统(大连)有限公司 Method for judging U-shaped region, long side and short side of CT image and method for repairing external outline of lung parenchyma

Also Published As

Publication number Publication date
CN104933709A (en) 2015-09-23

Similar Documents

Publication Publication Date Title
CN104933709B (en) Random walk CT lung tissue image automatic segmentation methods based on prior information
CN104751178B (en) Lung neoplasm detection means and method based on shape template matching combining classification device
CN108052977B (en) Mammary gland molybdenum target image deep learning classification method based on lightweight neural network
CN110309860B (en) Method for classifying malignancy degree of lung nodule based on convolutional neural network
CN107154043B (en) Pulmonary nodule false positive sample inhibition method based on 3DCNN
CN109934235B (en) Unsupervised abdominal CT sequence image multi-organ simultaneous automatic segmentation method
CN107067402B (en) Medical image processing apparatus and breast image processing method thereof
US20170169276A1 (en) Systems and methods for automated screening and prognosis of cancer from whole-slide biopsy images
CN109840913B (en) Method and system for segmenting tumor in mammary X-ray image
Maitra et al. Automated digital mammogram segmentation for detection of abnormal masses using binary homogeneity enhancement algorithm
CN112396619B (en) Small particle segmentation method based on semantic segmentation and internally complex composition
Celebi et al. Fast and accurate border detection in dermoscopy images using statistical region merging
CN105389821B (en) It is a kind of that the medical image cutting method being combined is cut based on cloud model and figure
Maitra et al. Detection of abnormal masses using divide and conquer algorithmin digital mammogram
Zidan et al. Level set-based CT liver image segmentation with watershed and artificial neural networks
CN104182755A (en) Mammary gland molybdenum target X-ray image block feature extraction method based on tower-shaped principal component analysis (PCA)
CN112598613A (en) Determination method based on depth image segmentation and recognition for intelligent lung cancer diagnosis
CN106504239B (en) A kind of method of liver area in extraction ultrasound image
CN109363697A (en) A kind of method and device of breast image lesion identification
CN111340770A (en) Method for constructing cancer prognosis model by combining global weighted LBP (local binary pattern) and texture analysis
CN112270667A (en) TI-RADS-based integrated deep learning multi-tag identification method
Soleymanifard et al. Segmentation of whole tumor using localized active contour and trained neural network in boundaries
CN116883341A (en) Liver tumor CT image automatic segmentation method based on deep learning
CN117036288A (en) Tumor subtype diagnosis method for full-slice pathological image
TW201726064A (en) Medical image processing apparatus and breast image processing method thereof

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20180914

Termination date: 20210604