Random field variational optimization medical image segmentation method and device and storage medium

Abstract

The invention relates to the field of medical image processing, and discloses a random field variational optimization medical image segmentation method and device and a storage medium, and the method comprises the steps: obtaining a medical image, and carrying out the solving operation of an initialization model; a solving model is constructed, the target of the model is given image gray Y, and the corresponding optimal segmentation X is solved; changing the model solution into the maximum value of the product of the likelihood probability and the prior probability by adopting a Bayesian formula; a simplified Potts model is adopted to solve a prior probability; the likelihood probability is equivalent to # imgabs0 #, a level set function phi (s) is introduced into the prior probability and the likelihood probability, a regular term is added to ensure the stability of the level set function in the evolution process, and finally a numerical calculation method is adopted for calculation to obtain a segmentation result. The method has the beneficial effects that the gray information and the spatial information of the image are considered at the same time, so that the segmentation result is more accurate.

Description

Technical Field

[0001] The present invention relates to the field of medical image processing, and in particular to a random field variational optimization medical image segmentation method, device and storage medium. Background Art

[0002] Due to the high complexity and ambiguity of the target morphology in medical images and the limitations of various imaging technologies, medical images are more seriously affected by a series of external factors such as noise, field offset, tissue movement, etc. Therefore, compared with natural images, the uneven distribution of grayscale in medical images is more obvious, and traditional image segmentation algorithms can hardly achieve accurate segmentation of medical images.

[0003] At present, the common image segmentation algorithms using the level set method can be divided into two categories: the improved algorithm based on the RSF (Region Scalable Fitting Model) model and the improved algorithm based on the CV (Chan and Vese Model) model. The improved algorithm based on the RSF model considers the local information of the image, which leads to the segmentation result containing many small local areas far away from the target area; the improved algorithm based on the CV model considers the global information of the image, so the model can roughly converge to the vicinity of the target contour, but when processing images with uneven grayscale distribution, the segmentation result of the model is not accurate enough. Summary of the invention

[0004] The purpose of the present invention is to propose a random field variational optimization medical image segmentation method, device and storage medium to solve the current technical problem of inaccurate image segmentation caused by uneven grayscale distribution of medical images.

[0005] The present invention provides a random field variational optimization medical image segmentation method.

[0006] The following steps are involved:

[0007] S1. Obtain a medical image S = {{i, j}|1≤i≤H, 1≤j≤W}, and perform an initialization model solution operation; wherein the total number of pixels of the image is N = H×W; (i, j) is any pixel point of the image;

[0008] S2, build a solution model, the goal of the model is to solve the corresponding optimal segmentation X for a given image grayscale Y;

[0009] S3, using the Bayesian formula to transform the model solution into the maximum value of the product of likelihood probability and prior probability;

[0010] S4, using the simplified Potts model to solve the prior probability;

[0011] S5. Equivalent likelihood probability to in is the local fitting energy term in Li's model;

[0012] S6. The level set function Φ(s) is introduced into the prior probability and likelihood probability, and a regular term is added to ensure the stability of the level set function during the evolution process. Finally, the numerical calculation method is used to obtain the segmentation result.

[0013] A storage medium stores instructions and data for implementing a random field variational optimization medical image segmentation method.

[0014] A random field variational optimization medical image segmentation device comprises: a processor and the storage medium; the processor loads and executes instructions and data in the storage medium to implement a random field variational optimization medical image segmentation method.

[0015] The beneficial effects provided by the present invention are: aiming at the characteristics of uneven grayscale of medical images, with the purpose of improving the segmentation accuracy of the model and enhancing the robustness of the model. A new model for medical image segmentation based on a statistical framework suitable for processing uneven grayscale. The innovation of this model lies in: 1) The new model takes into account both the grayscale information and spatial information of the image, which can make the segmentation result more accurate. 2) The category constraint term can better constrain the diffusion of the local fitting term, avoiding the problem of excessive diffusion of the segmentation contour in the traditional region-based model. 3) The corresponding numerical calculation algorithm of the new model is proposed, which enhances the practicability of the model. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 It is a schematic flow chart of the method of the present invention;

[0017] Figure 2 It is a schematic diagram of the distribution of neighboring regions in the original image domain (blue entity) and the transformed domain (red dashed line);

[0018] Figure 3 is a schematic diagram of the comparison of vascular image segmentation results (1);

[0019] Figure 4 is a schematic diagram of the comparison of vascular image segmentation results (2);

[0020] Figure 5 It is a schematic diagram comparing the segmentation results of different initial contours (1);

[0021] Figure 6 is a schematic diagram comparing the segmentation results of different initial contours (2);

[0022] Figure 7 is a schematic diagram comparing the segmentation results of different initial contours (3);

[0023] Figure 8 It is a schematic diagram comparing the segmentation results of different initial contours (4);

[0024] Fig. 9 It is a working schematic diagram of the hardware device of the present invention. DETAILED DESCRIPTION

[0025] To make the objectives, technical solutions and advantages of the present invention more clear, the embodiments of the present invention will be further described below with reference to the accompanying drawings.

[0026] Before formally explaining the present invention, relevant technical terms in the present invention are explained. In addition, the scheme of the present invention is also generally explained first for easy understanding.

[0027] Please refer to Figure 1 , Figure 1 It is a schematic diagram of the process of the present invention.

[0028] The present invention provides a random field variational optimization medical image segmentation method, comprising the following steps:

[0029] S1. Obtain a medical image S = {{i, j}|1≤i≤H, i≤j≤W}, and perform an initialization model solution operation; wherein the total number of pixels of the image is N = H×W; (i, j) is any pixel point of the image;

[0030] It should be noted that the set S={{i, j}|1≤i≤H, 1≤j≤W} is represented as an image in the present invention.

[0031] Each element in the set S represents a specific pixel, and the total number of pixels in the image is N = H × W, which represents the image domain. Each pixel s∈S has two random variables associated with it:

[0032] Categorical random variable X s ∈Ω≡{Ω 1 ,Ω 2 , ...Ω L}, L represents the number of categories;

[0033] Characteristic random variable Represents the dimension of the image.

[0034] The present invention will X s and Y s A specific random process is represented by a lowercase variable x s ∈Ω and express.

[0035] All random variables X for the entire image s and Y s Combine them with X=(X1 , X 2 , ...X N ) and Y=(Y 1 , Y 2 ,..Y N )express.

[0036] Therefore, the space of X and Y states can be expressed in Cartesian product form: N and

[0037] The specific random process of the X and Y random fields is represented by lowercase variables x = (x 1 , x 2 , ..., x N ) and y=(y 1 ,y 2 , ..., y N )express.

[0038] S2, build a solution model, the goal of the model is to solve the corresponding optimal segmentation X for a given image grayscale Y;

[0039] It should be noted that the solution model refers to the entire solution process mentioned later in this application, and is not aimed at a specific or concrete model. That is, the goal of this application is to solve the corresponding optimal segmentation X given the image grayscale Y.

[0040] S3, using the Bayesian formula to transform the model solution into the maximum value of the product of likelihood probability and prior probability;

[0041] As an example, y and x represent the observation results of random processes Y and X. Using the Bayesian formula, the posterior probability of the segmentation can be expressed as:

[0042]

[0043] The product of the likelihood probability P(Y=y|X=x) and the prior probability P(X=x) is proportional to the posterior probability P(X=x|Y=y). Therefore, finding the maximum posterior probability is equivalent to finding the maximum value of the product of the likelihood probability and the prior probability.

[0044] P(X s =x s ) indicates that an event {X s =x s}; and the prior probability P(X=x) means that the event The probability of occurrence, that is, all pixels s have an event {X s =x s}probability.

[0045] In order to avoid ambiguity, in the present invention, the random variables in the probability function will be omitted, for example, P(x) will be used to represent P(X=x).

[0046] S4, using the simplified Potts model to solve the prior probability;

[0047] Specifically, in order to simplify the prior probability, the present invention assumes that the image satisfies the following Markov properties:

[0048]

[0049] In the expression x -s ={x 1 , x 2 , ..., x s-1 , x s+1 ,...x N},

[0050] Expression set η s The i-th element in .

[0051] event The prior probability of can be approximated as follows:

[0052]

[0053] is a normalization constant used to ensure that the sum of the various probabilities is 1.

[0054] The present invention uses a simplified Potts model to calculate each conditional distribution in formula (2):

[0055]

[0056] The present invention regards the conditional distribution as: in the neighborhood of |η s |×|η s |, the energy of two specific states of the graph is calculated. s is the normalization constant, h(x s ) is when x s =Ω i The Hamiltonian when s ,Ω i ) represents the area η s ∩Ω i is the area of ​​​​the , λ is a constant.

[0057] S5. Equivalent likelihood probability to in is the local fitting energy term in Li's model;

[0058] It should be noted that, through the above steps, this application solves the problem of solving the prior probability, and also needs to solve the problem of solving the likelihood probability P(y|x).

[0059] Therefore, the present application will introduce the likelihood probability equivalent to step S5 as follows: The derivation principle behind it. Specifically:

[0060] Assume that for each pixel {X s =x s The value of} s are independent, that is, the values ​​of other pixels in the image will not affect the value Y of pixel s s Therefore, the problem of solving the likelihood probability P(y|x) is transformed into the likelihood probability P(y s |X s =Ω i ) problem. P(y|x) can be expressed as:

[0061]

[0062] Based on the retinex theory, Li et al. wrote the acquired image I in the form:

[0063] I=b·J+n (19)

[0064] I is the observed image; J is the real image; b is the bias field superimposed on the real image; n is Gaussian noise. In the retinex theory, Land et al. proposed the conditions that the real image J and the bias field b should satisfy in order to simplify the solution of the real image J:

[0065] (1) In the entire image domain Ω, the change process of the bias field b is very slow.

[0066] (2) In each independent region in the image domain, the grayscale of the real image J can be approximated by a constant c. The mathematical definition of this condition can be expressed as: for any s∈Ω i ( is an independent region of the image domain Ω), J(s)≈c i .

[0067] Based on the assumptions of Land et al., the present invention further assumes that:

[0068] In each independent region Ω i The noise in is zero mean and has a covariance of ∑ i Gaussian distribution.

[0069] Therefore, it is not difficult to conclude that: in each independent region Ω iSince the bias field b changes very slowly, the observed image I(s) can be approximately equal to μ i (s), the covariance is ∑ i (s) is a Gaussian distribution. Therefore, Xu Jun et al. s |X s =Ω i ) can be written as:

[0070]

[0071] Zhang et al.'s research shows that: the likelihood probability distribution curve of the above formula is as follows Figure 2 As shown by the blue solid curve in s Misclassification occurs when you fall into the tails of the distribution.

[0072] Zhang et al. proposed a method to reduce misclassification: while keeping the mean of the likelihood probability distribution unchanged, the cancellation profile of the above likelihood probability distribution is compressed to reduce the overlapping tails, thereby reducing the classification error rate (see Figure 2 The red dashed curve in Figure 5).

[0073] Pixel s is its neighborhood O s (O s = {c|||cs||≤ρ}), where ρ is the radius of the neighborhood. Since there are L non-overlapping independent regions Ω in the image i (Ω i is the ith independent region).

[0074] The image domain Ω can be expressed as: Ω = Ω 1 ∪Ω 2 ∪...∪Ω L ,and The present invention defines mapping This mapping transforms the image from the original intensity domain D(T) to another domain R(T):

[0075]

[0076] Assumptions To simplify the following derivation.

[0077] Since the intensity of each pixel s is independently distributed, The corresponding probability density function is still a Gaussian distribution

[0078] Reference Figure 2 As shown in the red dashed curve in , the overlapping tails of the distribution are suppressed to some extent.

[0079] Since intensity inhomogeneities appear as smooth intensity variations on the image, the present invention assumes Since the product of Gaussian PDFs is still Gaussian, we can get:

[0080]

[0081] When the components of image I are independent, the covariance matrix therefore

[0082]

[0083] Depend on We can get:

[0084]

[0085] Therefore, the joint probability can be expressed as:

[0086]

[0087] make We can get:

[0088]

[0089] among

[0090] From this we can get:

[0091]

[0092] Therefore, the likelihood probability

[0093] Obviously, when the image is a grayscale image, the likelihood function of the model of the present invention is equivalent to is the local fitting energy term in Li's model.

[0094] S6. The level set function Φ(s) is introduced into the prior probability and likelihood probability, and a regular term is added to ensure the stability of the level set function during the evolution process. Finally, the numerical calculation method is used to obtain the segmentation result.

[0095] It should be noted that after completing steps S4 and S5, this application will proceed to the final solution process.

[0096] Specifically, the solution process is as follows:

[0097] In order to simplify the estimation of formula (1), the following transformation is made:

[0098] -logP(x|y)=-(logP(y|x)+logP(x)) (28)

[0099] The problem of finding the maximum value of equation (1) can be transformed into the problem of finding the minimum value of equation (14). Substituting equations (3) and (13) into equation (14), we can obtain the following energy functional form:

[0100]

[0101] is a constant. Introduce function K ρ (s, c):

[0102]

[0103] Substituting equation (8) and (16) into equation (15), we can obtain:

[0104]

[0105] The present invention selects the level set function Φ to represent the two independent regions of foreground and background: foreground Ω 1 ={Φ(s)≥0} and background Ω 2 ={Φ(s)<0}.

[0106] The foreground and background areas are represented by the Heaviside function:

[0107]

[0108] The present invention introduces the level set function Φ(s) into the likelihood probability and the prior probability, so it is necessary to add a regularization term to ensure the stability of the level set function during the evolution process. According to the literature "Level set evolution without re-initialization: a new variational formulation", the level set regularization term can be written as follows:

[0109]

[0110] make

[0111] This gives the following energy functional:

[0112]

[0113] The first term in the energy functional of formula (20) is called the local fitting energy term, the second term is called the category constraint term, the third term is called the level set regularization term, and the last term is a constant term.

[0114] Formula (20) contains four variables: Φ, b, c, and ∑. The present invention solves the Euler-Lagrange equation that must be satisfied for formula (20) to achieve the minimum value under this condition by fixing any three of the four variables: Φ, b, c, and ∑.

[0115] By fixing Φ, c, ∑, the present invention obtains the equation that the bias field b should satisfy:

[0116]

[0117] It's the convolution symbol, it's the point multiplication symbol.

[0118] By fixing Φ, b, ∑, the present invention obtains the independent region c i The equations that should be satisfied are:

[0119]

[0120] By fixing Φ, b, c, the present invention obtains the covariance ∑ i The equations that should be satisfied are:

[0121]

[0122] By fixing b, c, ∑, the present invention uses the gradient descent method formula to find the equation that the level set function Φ should satisfy:

[0123]

[0124] δ(Φ) is the Dirac function,

[0125] From this we can get the following evolution equation:

[0126]

[0127] The numerical calculation algorithm of this model can be divided into the following five steps:

[0128] (1) Set the parameter value and the number of iterations k = 1, and initialize the level set function Φ:

[0129]

[0130] where c 0 is a positive constant, and R is any given initialization region in the image domain Ω.

[0131] (2) Initialize the bias field b and set the size of the neighborhood ρ.

[0132] (3) Calculate b, c, ∑ according to equations (21), (22), and (23).

[0133] (4) Calculate Φ according to formula (25).

[0134] (5) Determine whether the iteration result converges. If not, set k=k+1 and return to step 3.

[0135] It should be noted that in numerical calculations, the Heaviside function is given by the following H ε Function to approximate:

[0136]

[0137] Function H ε The derivative of the Dirac function is given by δ ε Function to approximate:

[0138]

[0139] The present invention selects function λ=2λ 0 arctan[0.2(t-1)] / π to update the value of parameter λ (t represents the number of iterations, λ 0 The time step Δt is set to 0.4. The weight parameter α of the level set regularization term is set to 0.1 / Δt. The neighborhood η s is a 21×21 square window centered at pixel s but not including pixel s. The kernel size parameter ρ of the Gaussian kernel function is set to 11. The value of ε is set to 1.

[0140] In order to verify the solution proposed by the present invention, the present invention chooses to use the indicators of Jaccard similarity (JS), Dice similarity coefficient (DSC), false positive rate (RFP) and false negative rate (RFN) to quantitatively compare the segmentation results of different models. The definitions of these indicators are as follows:

[0141]

[0142] N(·) represents the number of pixels in the area, S g represents the foreground area of ​​the real image, S m Represents the foreground area obtained by the model. O is S g and S m common part.

[0143] In order to quantitatively compare the segmentation effects of the models, the present invention selects blood vessel images in medical images for experiments to verify the effectiveness of the model of the present invention.

[0144] Please refer to Figure 3 , Figure 4 ,according to Figure 3 , Figure 4The segmentation results and the quantitative evaluation results in Table 1 show that when processing a blood vessel image with uneven grayscale distribution, the model of the present invention can more accurately segment the area near the edge and does not include local areas outside the target area.

[0145] Table 1 Comparison of index values ​​of segmentation results (1)

[0146]

[0147] The present invention selects blood vessel images and brain images in medical images, and conducts experiments under different initial contour conditions to verify the robustness of the model of the present invention.

[0148] according to Figure 5-Figure 8 The experimental results show that the model of the present invention has stronger generalization ability and is insensitive to the selection of the initial contour, no matter whether the model is segmenting a vascular image or a brain image; no matter whether the initial contour is selected in the foreground area, the background area, or the area contained in both the foreground and the background, the model of the present invention can converge well to the vicinity of the foreground edge, and the segmentation results are also better than those of Xu's model and Zhang's model.

[0149] Experiments show that the time step not only affects the evolution rate of the model of the present invention, but also affects the segmentation results of the model. Therefore, the present invention focuses on the influence of the time step on the image segmentation accuracy and the curve evolution rate. In order to be able to quantitatively compare the segmentation results, the present invention still chooses to use blood vessel images for experiments.

[0150] In order to more quantitatively compare the effects of different time steps on segmentation accuracy, the present invention chooses to calculate various indicators of the segmented contour relative to the hand-drawn contour after each iteration to accurately measure the segmentation result.

[0151] According to the final experimental results, the model of the present invention has an optimal time step setting of about 0.5, which can take into account both the accuracy and segmentation rate of image segmentation.

[0152] See also Fig. 9 , Fig. 9 40 is a schematic diagram of the working of the hardware device of an embodiment of the present invention, wherein the hardware device specifically comprises: a random field variational optimization medical image segmentation device 401, a processor 402 and a storage medium 403.

[0153] A random field variational optimization medical image segmentation device 401: The random field variational optimization medical image segmentation device 401 implements the random field variational optimization medical image segmentation method.

[0154] Processor 402: The processor 402 loads and executes the instructions and data in the storage medium 403 to implement the random field variational optimization medical image segmentation method.

[0155] Storage medium 403: The storage medium 403 stores instructions and data; the storage medium 403 is used to implement the random field variational optimization medical image segmentation method.

[0156] The beneficial effects of the present invention are as follows: the designed model takes into account both the grayscale information and the spatial information of the image, which can make the segmentation result more accurate; the category constraint term can better constrain the diffusion of the local fitting term, avoiding the problem of excessive diffusion of the segmentation contour in the traditional region-based model; and a numerical calculation method corresponding to the new model is proposed, which enhances the practicability of the model.

[0157] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A random field variational optimization medical image segmentation method, characterized by: The method comprises the following steps: S1. Obtain a medical image S = {{i, j}|1≤i≤H, 1≤j≤W}, and perform an initialization model solving operation; wherein the total number of pixels of the image is N = H×W; (i,j) is any pixel in the image; S2, build a solution model, the goal of the model is to solve the corresponding optimal segmentation X for a given image grayscale Y; S3, using the Bayesian formula to transform the model solution into the maximum value of the product of likelihood probability and prior probability; S4, using the simplified Potts model to solve the prior probability; S5. Equivalent likelihood probability to in is the local fitting energy term in Li's model; S6. The level set function Φ(s) is introduced into the prior probability and likelihood probability, and a regular term is added to ensure the stability of the level set function during the evolution process. Finally, the numerical calculation method is used to obtain the segmentation result.

2. The random field variational optimization medical image segmentation method according to claim 1, characterized in that: The initialization model solving operation in step S1 is as follows: S11. Use Ω to represent the image domain; each pixel s∈S has two random variables associated with it: the category random variable X s ∈Ω≡{Ω1,Ω2,...Ω L }, L represents the number of categories; Characteristic random variable Represents the dimension of the image; S12. Change the categorical random variable X s and the characteristic random variable Y s A specific random process is represented by a lowercase variable x s ∈Ω and Represents; random variables X of all the entire images s and Y s Combine them with X = (X1, X2, ...X N ) and Y=(Y1,Y2,..Y N )express; S13. The space of X and Y states can be expressed in the form of Cartesian product: Ω N and The specific random process of the X and Y random fields is represented by lowercase variables x=(x1, x2, ..., x N ) and y=(y1,y2,...,y N )express.

3. A random field variational optimization medical image segmentation method as claimed in claim 2, characterized in that: The specific calculation in step S3 is as follows: Among them, the product of the likelihood probability P(Y=y|X=x) and the prior probability P(X=x) is proportional to the posterior probability P(X=x|Y=y); finding the maximum posterior probability is equivalent to finding the maximum value of the product of the likelihood probability and the prior probability.

4. The random field variational optimization medical image segmentation method according to claim 3, characterized in that: Step S4 is specifically as follows: S41, using P(X s =x s ) indicates that an event {X s =x s }Probability; S42, omitting the random variable in the probability function, and using P(x) to represent P(X=x); S43, using the image S to satisfy the Markov property, the prior probability is simplified as follows: is a normalization constant used to ensure that the sum of all probabilities is 1; S44. Use the simplified Potts model to calculate each conditional distribution in formula (2): The conditional distribution is considered as: in the neighborhood of |η s |×|η s |The energy calculation of two specific states of the graph; Z s is the normalization constant, h(x s ) is when x s =Ω i The Hamiltonian when s ,Ω i ) represents the area η s ∩Ω i The area of ​​​​the product of the prior probability is λ, which is a constant.

5. The random field variational optimization medical image segmentation method according to claim 4, characterized in that: Step S6 is specifically as follows: S61, transform formula (1) as follows: -logP(x|y)=-(logP(y|x)+logP(x)) (4) The maximum value problem of formula (1) can be transformed into the minimum value problem of formula (14); S62. Obtain the energy functional form of formula (14): Among them, I = b·J+n, I is the observed image; J is the real image; b is the bias field superimposed on the real image; n is Gaussian noise; S63, the foreground and background regions are represented by the Heaviside function: S64. Introduce the level set function Φ(s) in the likelihood probability and prior probability, and add the level set regularization term to ensure the stability of the level set function during the evolution process. The level set regularization term can be written as follows: make This gives the following energy functional: The first term in the energy functional of formula (20) is called the local fitting energy term, the second term is called the category constraint term, the third term is called the level set regularization term, and the last term is a constant term; S65, by fixing any three of the four variables Φ, b, c, Σ, solve the Euler-Lagrange equation that equation (20) should satisfy to obtain the minimum value under this condition; S66. Use numerical calculation method to solve the Euler-Lagrange equation.

6. The random field variational optimization medical image segmentation method according to claim 5, characterized in that: In step S65, by fixing Φ, c, ∑, the equation that the bias field b should satisfy is obtained: is the convolution symbol, is the point multiplication symbol; By fixing Φ, b, ∑, we can find the independent region c i The equations that should be satisfied are: By fixing Φ, b, c, we can find the covariance ∑ i The equations that should be satisfied are: By fixing b, c, ∑, and applying the gradient descent formula, we can find the equation that the level set function Φ should satisfy: δ(Φ) is the Dirac function, From this we can get the following evolution equation:

7. The random field variational optimization medical image segmentation method according to claim 6, characterized in that: The numerical calculation method in step S66 is specifically as follows: S661, set parameter values, number of iterations k = 1, and initialize the level set function Φ: Where c0 is a positive constant, R is any given initialization region in the image domain Ω; S662, initialize the bias field b and set the size of the neighborhood ρ; S663, calculate b, c, Σ according to equations (21), (22), and (23); S664, Calculate Φ according to formula (25); S665. Determine whether the iteration result converges; if not, set k=k+1 and return to step S663.

8. A storage medium, characterized in that: The storage medium stores instructions and data for implementing a random field variational optimization medical image segmentation method as described in any one of claims 1 to 7.

9. A random field variational optimization medical image segmentation device, characterized by: include: A processor and a storage medium; the processor loads and executes instructions and data in the storage medium to implement a random field variational optimization medical image segmentation method as described in any one of claims 1 to 7.

Images