CN111028201B - Fundus blood vessel image segmentation system and method based on multi-scale level set - Google Patents

Abstract

The invention provides a fundus blood vessel image segmentation system and method based on a multi-scale level set, and relates to the technical field of image processing. The device comprises an input module, a preprocessing module, a scale construction module, a level set module and an output module; the present invention improves upon existing level set methods based on local features. The fundus blood vessels are different in thickness, and the blood vessel segmentation granularity is controlled by using a single scale, so that the problem of low blood vessel segmentation accuracy is caused. Based on the characteristic of different thicknesses of blood vessels, the optimal response scale is found for each pixel point, namely, the scales of the blood vessels with different thicknesses are different, and the problems of missing division, wrong division and the like of the blood vessels with different thicknesses are overcome, so that the robustness of a segmentation method is enhanced, and an accurate image segmentation result can be obtained by applying the method.

Description

Fundus blood vessel image segmentation system and method based on multi-scale level set

Technical Field

The invention relates to the technical field of image processing, in particular to a fundus blood vessel image segmentation system and method based on a multi-scale level set.

Background

Fundus blood vessel images are the basis for screening eye diseases such as cataract, glaucoma, age-related macular degeneration and the like. The problems of complex fundus blood vessel structure and uneven image gray scale result in the problems of insufficient segmentation of a plurality of algorithms at adjacent blood vessels. On the other hand, the fundus blood vessels are different in thickness, and the existing segmentation method mostly uses a single scale to control the blood vessel segmentation granularity, so that the problem of low blood vessel segmentation accuracy is caused. Therefore, the invention provides a fundus blood vessel image segmentation method based on a multi-scale level set. And selecting the optimal scale for each pixel point, and overcoming the problem of missed division and wrong division of blood vessels with uneven thickness.

Disclosure of Invention

Aiming at the problems existing in the prior art, the invention provides a fundus blood vessel image segmentation system based on a multi-scale level set.

The technical scheme adopted by the invention is as follows: a fundus blood vessel image segmentation system based on a multi-scale level set comprises an input module, a preprocessing module, a scale construction module, a level set module and an output module;

the input module is used for receiving fundus blood vessel images input by a user and then transmitting the fundus blood vessel images to the preprocessing module;

the preprocessing module is used for expanding the received fundus blood vessel image, and the expanded result is transmitted to the scale construction module;

the scale construction module is used for finding the optimal vascular response scale for the ocular fundus vascular image, constructing a Gaussian curve based on the scale, subtracting the original image from the convolution result of the Gaussian curve and the image to obtain an image M, and transmitting the image M to the level set module;

the level set module is used for carrying out blood vessel segmentation on the image M by adopting a level set algorithm, and transmitting a segmentation result to the output module;

the output module is used for outputting and displaying the images of the division results.

On the other hand, the invention provides a fundus blood vessel image segmentation method based on a multi-scale level set, which is realized by the fundus blood vessel image segmentation system based on the multi-scale level set, and comprises the following steps of:

step 1: the input module receives a fundus blood vessel image A to be segmented, wherein the width of the image A is W, and the height of the image A is H;

step 2: the preprocessing module constructs the received fundus blood vessel image A to be segmented to generate a fundus blood vessel image I without edge background, and the resolution of the fundus blood vessel image I is the same as that of the image A; converting the image A into a gray image J, setting the whole area of the image J as omega, directly obtaining an image P, I=J if the area to be segmented is full of omega, and executing the step 3; if the region to be segmented is not full of omega, mirror image expansion is carried out on the effective boundary region in the image to cover the black filling region in the whole image, the expanded image is P, meanwhile, I=P, and step 2.1 is executed; the method comprises the steps of carrying out a first treatment on the surface of the

Step 2.1: determining a threshold Q to divide the gray image J into a foreground R and a background B, wherein R and B are binary images, r=1 at the foreground region, b=0, r=0 at the background region, b=1;

step 2.2: constructing foreground region boundary image R ₂ Corroding the prospect R to obtain R ₁ ，R ₂ ＝R-R ₁ Wherein the image R ₂ Is a binary image;

step 2.3: an initialization variable p=1, where P is an initialization abscissa value of the pixel point in the outer loop;

step 2.4: if p is less than or equal to W, and the initialized ordinate value q=1 of the pixel point in the outer loop, executing the step 2.5, otherwise, executing the step 3;

step 2.5: if q is less than or equal to H, executing step 2.6, otherwise, making p=p+1, and executing step 2.4;

step 2.6: if the current point (p, q) satisfies b=1, executing step 2.7, otherwise, let q=q+1, executing step 2.5;

step 2.7: initializing variable m=1, m is the initial abscissa value of pixel point in image inner circulation, and initializingWherein r is the distance between any two pixels in the initialized image, and the initial coordinates (a, b) are two random positive integers;

step 2.8: if m is less than or equal to W, the variable n=1, n is the initial value of the ordinate of the pixel point in the inner loop, executing the step 2.9, otherwise, executing the step 2.13;

step 2.9: if n is less than or equal to H, executing the step 2.10, otherwise, making m=m+1, and executing the step 2.8;

step 2.10: if the point (m, n) satisfies R ₂ Step 2.11 is performed, otherwise, let n=n+1, step 2.9 is performed;

step 2.11: calculating a distance dist ((p, q), (m, n)) between the point (p, q) and the point (m, n), wherein dist is a distance metric function, if r < dist ((p, q), (m, n)), r=dist ((p, q), (m, n)), while letting a=m, b=n;

step 2.12: let n=n+1, jump to step 2.9;

step 2.13: the points (p, q) take the points (a, b) as the centers, R is the symmetrical distance, a corresponding pixel point L in R is obtained, and the value of the pixel point L is assigned to the current point Z; wherein L is _x ＝2a-p，L _y =2b_q, where L _x Is the abscissa L of the pixel point L _y Is the ordinate of the pixel point L;

step 2.14: let q=q+1, jump to step 2.5;

step 3: for the parameter minimum scale sigma _min Maximum dimension sigma _max Setting constant coefficient gamma, constant coefficient c and scale step delta sigma, wherein sigma _min <σ _max ；

Step 4: for the initial dimension sigma _i Setting sigma _i ＝σ _min ；

Step 5: constructing a Hessian matrix of the image I, constructing a Gaussian function through an optimal response scale sigma, and subtracting the image I from the result of convolution operation of the image I and the Gaussian function to obtain an image M, wherein the method comprises the following specific steps of:

step 5.1: initializing a variable p=1, a scale variable k=1, and a variable v _max＝0, wherein v_max Is the maximum response value variable of the pixel;

step 5.2: if p is less than or equal to W, the variable q=1, the step 5.3 is executed, otherwise, the step 6 is executed;

step 5.3: if q is less than or equal to H, executing step 5.4, otherwise, making p=p+1, and executing step 5.2;

step 5.4: judging the current scale sigma _i Whether or not it is greater than the maximum dimension sigma _max If greater than sigma _max Let sigma _i ＝σ _min Executing the step 5.12, otherwise, executing the step 5.5;

step 5.5: constructing a Hessian matrix of points (p, q):

wherein ,

D＝∫∫I(p-u,q-v)×A(u,v)dudv

E＝∫∫I(p-u,q-v)×B(u,v)dudv

F＝∫∫I(p-u,q-v)×C(u,v)dudv

a (u, v), B (u, v) and C (u, v) are based on the scale sigma _i Construction of the second partial derivative of Gaussian

Wherein D, E, F is the convolution of image I with the corresponding gaussian second partial derivative A, B, C;

step 5.6: calculating two eigenvalues lambda of the Hessian matrix obtained in step 5.5 ₁ 、λ ₂, wherein

wherein

Step 5.7: calculating at the current scale sigma _i The response value V of the pixel (p, q) _σi ；

wherein

in the formula ,R_B Is the degree of anisotropy of the linear shape, S is the norm;

step 5.8: if it isMake->And κ=σ _i ；

Step 5.9: constructing a Gaussian kernel based on the scale kappa to obtain a Gaussian curve g _κ ；

Step 5.10: the Gaussian curve g obtained in the step 5.9 _κ Performing convolution operation with the image I, and subtracting the image I from the convolution result to obtain an image M;

M(p,q)＝(g _κ *I)(p,q)-I(p,q)

step 5.11: updating the scale sigma _i ＝σ _i And (5) jumping to the step 5.4, wherein the delta sigma is the scale step length;

step 5.12: updating the coordinates q=q+1, and jumping to step 5.3;

step 6, setting an initial zero level set segmentation curve on the image M, and generating a high-one-dimensional level set function phi according to the set zero level set;

the set initial zero level set segmentation curve is a closed curve with a specified shape and is a closed curve generated randomly, and the level set function phi is a high one-dimensional symbol distance function generated according to the set initial level set segmentation curve;

step 7, establishing energy functional E (phi) of level set segmentation and setting a weighting coefficient alpha in the level set method ₁ Weighting coefficient alpha ₂ Curve smoothing term coefficient v and area smoothing term coefficient mu;

the energy functional is as follows:

E(φ)＝E _D +vL(φ)+μP(φ)

wherein ,E_D For image data items, L (φ) is a length constraint, P (φ) is an area constraint, and H (φ (P, q)) is a Heaviside function of φ (P, q);

E _D (φ)＝α ₁ ∫ _Ω |M(p,q)-c ₁ | ² H(φ(p,q))dpdq+α ₂ ∫ _Ω |M(p,q)-c ₂ | ² (1-H(φ(p,q)))dpdq

P(φ)＝∫H(φ(p,q))dpdq

wherein c₁ Is the gray average value in the zero level set segmentation curve, c ₂ Is the gray average value outside the zero level set segmentation curve;

step 8: setting a convergence threshold epsilon of the maximum iteration number J and the energy functional E (phi) of the image level set;

step 9: updating the image level set function set H (phi) separately ^j )、1-H(φ ^j )、And->φ ^j Obtaining the energy functional E of the jth iteration ^j ；

Specifically:

where Δt is the time step size,a time partial derivative of the level set function that is iteratively updated, the time partial derivative formula being as follows:

step 10: judging whether the iteration termination condition is reached: if the current iteration number J reaches the maximum iteration number J or the difference value between the iteration results of two adjacent times of the energy functional E (phi) is smaller than a convergence threshold epsilon set in advance, executing the step 11, otherwise, adding 1 to the iteration number, and returning to the step 9;

step 11: from a current updated image level set function phi ^j An image, i.e. a segmentation result of the image, is constructed.

The beneficial effects of adopting above-mentioned technical scheme to produce lie in:

the invention provides a segmentation system and a segmentation method of fundus blood vessel images based on a multi-scale level set, which improve the existing level set method based on local characteristics. The fundus blood vessels are different in thickness, and the blood vessel segmentation granularity is controlled by using a single scale, so that the problem of low blood vessel segmentation accuracy is caused. Based on the characteristic of different thicknesses of the blood vessels, an optimal response scale is found for each pixel point, namely, the scales of the blood vessels with different thicknesses are different, the problems of missing division, wrong division and the like of the blood vessels with different thicknesses are overcome, and therefore the robustness of the segmentation method is enhanced. By the method, an accurate image segmentation result can be obtained.

Drawings

FIG. 1 is a block diagram of a fundus blood vessel image segmentation system;

fig. 2 is a flowchart of a fundus blood vessel image segmentation method

Fig. 3 shows an original image of fundus blood vessel image input;

fig. 4 shows the conversion of fundus blood vessel image artwork into a grey scale image;

fig. 5 is a view showing an expanded fundus blood vessel image;

FIG. 6 shows an enhanced image;

FIG. 7 shows the image after convolution minus the original;

FIG. 8 shows an initial level set;

FIG. 9 shows the final segmentation result of an image;

fig. 10 shows the gold standard of the image.

Detailed Description

The following describes in further detail the embodiments of the present invention with reference to the drawings and examples. The following examples are illustrative of the invention and are not intended to limit the scope of the invention.

A fundus blood vessel image segmentation system based on a multi-scale level set is shown in figure 1, and comprises an input module, a preprocessing module, a scale construction module, a level set module and an output module;

the input module is used for receiving fundus blood vessel images input by a user and then transmitting the fundus blood vessel images to the preprocessing module;

the preprocessing module is used for expanding the received fundus blood vessel image, and the expanded result is transmitted to the scale construction module;

the scale construction module is used for finding the optimal vascular response scale for the ocular fundus vascular image, constructing a Gaussian curve based on the scale, subtracting the original image from the convolution result of the Gaussian curve and the image to obtain an image M, and transmitting the image M to the level set module;

the level set module is used for carrying out blood vessel segmentation on the image M by adopting a level set algorithm, and transmitting a segmentation result to the output module;

the output module is used for outputting and displaying the images of the division results.

In another aspect, the present invention provides a method for segmenting fundus blood vessel images based on a multi-scale level set, which is implemented by the aforementioned fundus blood vessel image segmentation system based on a multi-scale level set, as shown in fig. 2, and includes the following steps:

step 1: the input module receives a fundus blood vessel image A to be segmented, wherein the width of the image A is W, and the height of the image A is H as shown in fig. 3;

step 2: the preprocessing module constructs the received fundus blood vessel image A to be segmented to generate a fundus blood vessel image I without edge background, and the resolution of the fundus blood vessel image I is the same as that of the image A; converting the image A into a gray image J, setting the whole area of the image J as omega as shown in fig. 4, directly obtaining an image P, I=J if the area to be segmented is full of omega, and executing the step 3; if the region to be segmented is not full of omega, mirror image expansion is carried out on the effective boundary region in the image to cover the black filling region in the whole image, the expanded image is P, meanwhile, I=P, and step 2.1 is executed; as shown in fig. 5;

step 2.1: determining a threshold Q to divide the gray image J into a foreground R and a background B, wherein R and B are binary images, r=1 at the foreground region, b=0, r=0 at the background region, b=1;

step 2.2: constructing foreground region boundary image R ₂ Corroding the prospect R to obtain R ₁ ，R ₂ ＝R-R ₁ Wherein the image R ₂ Is a binary image;

step 2.3: an initialization variable p=1, where P is an initialization abscissa value of the pixel point in the outer loop;

step 2.4: if p is less than or equal to W, and the initialized ordinate value q=1 of the pixel point in the outer loop, executing the step 2.5, otherwise, executing the step 3;

step 2.5: if q is less than or equal to H, executing step 2.6, otherwise, making p=p+1, and executing step 2.4;

step 2.6: if the current point (p, q) satisfies b=1, executing step 2.7, otherwise, let q=q+1, executing step 2.5;

step 2.7: initializing variable m=1, m is the initial abscissa value of pixel point in image inner circulation, and initializingWherein r is the distance between any two pixels in the initialized image, and the initial coordinates (a, b) are two random positive integers;

step 2.8: if m is less than or equal to W, the variable n=1, n is the initial value of the ordinate of the pixel point in the inner loop, executing the step 2.9, otherwise, executing the step 2.13;

step 2.9: if n is less than or equal to H, executing the step 2.10, otherwise, making m=m+1, and executing the step 2.8;

step 2.10: if the point (m, n) satisfies R ₂ Step 2.11 is performed, otherwise, let n=n+1, step 2.9 is performed;

step 2.11: calculating a distance dist ((p, q), (m, n)) between the point (p, q) and the point (m, n), wherein dist is a distance metric function, if r < dist ((p, q), (m, n)), r=dist ((p, q), (m, n)), while letting a=m, b=n;

step 2.12: let n=n+1, jump to step 2.9;

step 2.13: the points (p, q) take the points (a, b) as the centers, R is the symmetrical distance, a corresponding pixel point L in R is obtained, and the value of the pixel point L is assigned to the current point Z; wherein L is _x ＝2a-p，L _y =2b_q, where L _x Is the abscissa L of the pixel point L _y Is the ordinate of the pixel point L;

step 2.14: let q=q+1, jump to step 2.5;

step 3: for the parameter minimum scale sigma _min Maximum dimension sigma _max Setting constant coefficient gamma, constant coefficient c and scale step delta sigma, wherein sigma _min <σ _max ；

Step 4: for the initial dimension sigma _i Setting sigma _i ＝σ _min ；

Step 5: constructing a Hessian matrix of the image I, constructing a Gaussian function through an optimal response scale sigma, and subtracting the image I from the result of convolution operation of the image I and the Gaussian function to obtain an image M, wherein the specific steps are as follows:

step 5.1: initializing a variable p=1, a scale variable k=1, and a variable v _max＝0, wherein v_max Is the maximum response value variable of the pixel;

step 5.2: if p is less than or equal to W, the variable q=1, the step 5.3 is executed, otherwise, the step 6 is executed;

step 5.3: if q is less than or equal to H, executing step 5.4, otherwise, making p=p+1, and executing step 5.2;

step 5.4: judging the current scale sigma _i Whether or not it is greater than the maximum dimension sigma _max If greater than sigma _max Let sigma _i ＝σ _min Executing the step 5.12, otherwise, executing the step 5.5;

step 5.5: constructing a Hessian matrix of points (p, q):

wherein ,

D＝∫∫I(p-u,q-v)×A(u,v)dudv

E＝∫∫I(p-u,q-v)×B(u,v)dudv

F＝∫∫I(p-u,q-v)×C(u,v)dudv

a (u, v), B (u, v) and C (u, v) are based on the scale sigma _i Constructing a Gaussian second partial derivative;

wherein D, E, F is the convolution of image I with the corresponding gaussian second partial derivative A, B, C;

step 5.6: calculating two eigenvalues lambda of the Hessian matrix obtained in step 5.5 ₁ 、λ ₂, wherein

wherein

Step 5.7: calculating at the current scale sigma _i The response value V of the pixel (p, q) _σi ；

wherein

in the formula ,R_B Is the degree of anisotropy of the linear shape, S is the norm;

step 5.8: if it isMake->And κ=σ _i ；

Step 5.9: constructing a Gaussian kernel based on the scale kappa to obtain a Gaussian curve g _κ ；

Step 5.10: the Gaussian curve g obtained in the step 5.9 _κ Performing convolution operation with the image I, and subtracting the image I from the convolution result to obtain an image M;

M(p,q)＝(g _κ *I)(p,q)-I(p,q)

step 5.11: updating the scale sigma _i ＝σ _i And (2) jumping to the step size, wherein the delta sigma is the scale step sizeStep 5.4;

step 5.12: updating the coordinates q=q+1, and jumping to step 5.3;

step 6, setting an initial zero level set segmentation curve on the image M, and generating a high-one-dimensional level set function phi according to the set zero level set as shown in fig. 8;

the set initial zero level set segmentation curve is a closed curve with a specified shape and is a closed curve generated randomly, and the level set function phi is a high one-dimensional symbol distance function generated according to the set initial level set segmentation curve;

step 7, establishing energy functional E (phi) of level set segmentation and setting a weighting coefficient alpha in the level set method ₁ Weighting coefficient alpha ₂ Curve smoothing term coefficient v and area smoothing term coefficient mu;

the energy functional is as follows:

E(φ)＝E _D +vL(φ)+μP(φ)

wherein ,E_D For image data items, L (φ) is a length constraint, P (φ) is an area constraint, and H (φ (P, q)) is a Heaviside function of φ (P, q);

E _D (φ)＝α ₁ ∫ _Ω |M(p,q)-c ₁ | ² H(φ(p,q))dpdq+α ₂ ∫ _Ω |M(p,q)-c ₂ | ² (1-H(φ(p,q)))dpdq

P(φ)＝∫H(φ(p,q))dpdq

wherein c₁ Is the gray average value in the zero level set segmentation curve, c ₂ Is the gray average value outside the zero level set segmentation curve;

in this embodiment, the larger the image data item coefficient is, the closer the zero level set segmentation curve is to the boundary of the object to be segmented in the image, the smoother the segmentation curve obtained by the larger the level set segmentation curve length item coefficient is, so that the control parameters of image level set segmentation are set according to the image to be segmented: weighting coefficient alpha ₁ ＝1.0，Weighting coefficient alpha ₂ The curve smoothing term coefficient v=0.001×255×255, and the area smoothing term coefficient μ=2.0=1.0.

Step 8: setting a convergence threshold epsilon of the maximum iteration number J and the energy functional E (phi) of the image level set;

in the present embodiment, the maximum iteration number J is set to 50, and the energy functional E (c ₁ ,c ₂ Phi) is 0.01.

Step 9: updating the image level set function set H (phi) separately ^j )、1-H(φ ^j )、And->φ ^j Obtaining the energy functional E of the jth iteration ^j ；

Specifically:

where Δt is the time step size,a time partial derivative of the level set function that is iteratively updated, the time partial derivative formula being as follows:

step 10: judging whether the iteration termination condition is reached: if the current iteration number J reaches the maximum iteration number J or the difference value between the iteration results of two adjacent times of the energy functional E (phi) is smaller than a convergence threshold epsilon set in advance, executing the step 11, otherwise, adding 1 to the iteration number, and returning to the step 9;

step 11: from a current updated image level set function phi ^j The image, i.e., the segmentation result of the image, is constructed as shown in fig. 9, and fig. 10 is a gold standard of the image.

Finally, it should be noted that: the above embodiments are only for illustrating the technical solution of the present invention, and not for limiting the same; although the invention has been described in detail with reference to the foregoing embodiments, it will be understood by those of ordinary skill in the art that: the technical scheme described in the foregoing embodiments can be modified or some or all of the technical features thereof can be replaced by equivalents; such modifications and substitutions do not depart from the spirit of the corresponding technical solutions, which are defined by the scope of the appended claims.

Claims (1)

1. Fundus blood vessel image segmentation system based on multiscale level set, which is characterized in that: the system comprises an input module, a preprocessing module, a scale construction module, a level set module and an output module;

the input module is used for receiving fundus blood vessel images input by a user and then transmitting the fundus blood vessel images to the preprocessing module;

the preprocessing module is used for expanding the received fundus blood vessel image, and the expanded result is transmitted to the scale construction module;

the scale construction module is used for finding the optimal vascular response scale for the ocular fundus vascular image, constructing a Gaussian curve based on the scale, subtracting the original image from the convolution result of the Gaussian curve and the image to obtain an image M, and transmitting the image M to the level set module;

the level set module is used for carrying out blood vessel segmentation on the image M by adopting a level set algorithm, and transmitting a segmentation result to the output module;

the output module is used for outputting and displaying the images of the division results;

the fundus blood vessel image segmentation system based on the multi-scale level set is used for realizing a fundus blood vessel image segmentation method based on the multi-scale level set, and comprises the following steps of:

step 1: the input module receives a fundus blood vessel image A to be segmented, wherein the width of the image A is W, and the height of the image A is H;

step 2: the preprocessing module constructs the received fundus blood vessel image A to be segmented to generate a fundus blood vessel image I without edge background, and the resolution of the fundus blood vessel image I is the same as that of the image A; converting the image A into a gray image J, setting the whole area of the image J as omega, directly obtaining an image P, I=J if the area to be segmented is full of omega, and executing the step 3; if the region to be segmented is not full of omega, mirror image expansion is carried out on the effective boundary region in the image to cover the black filling region in the whole image, the expanded image is P, meanwhile, I=P, and step 2.1 is executed;

step 2.1: determining a threshold Q to divide the gray image J into a foreground R and a background B, wherein R and B are binary images, r=1 at the foreground region, b=0, r=0 at the background region, b=1;

step 2.2: constructing foreground region boundary image R ₂ Corroding the prospect R to obtain R ₁ ，R ₂ ＝R-R ₁ Wherein the image R ₂ Is a binary image;

step 2.3: an initialization variable p=1, where P is an initialization abscissa value of the pixel point in the outer loop;

step 2.4: if p is less than or equal to W, and the initialized ordinate value q=1 of the pixel point in the outer loop, executing the step 2.5, otherwise, executing the step 3;

step 2.5: if q is less than or equal to H, executing step 2.6, otherwise, making p=p+1, and executing step 2.4;

step 2.6: if the current point (p, q) satisfies b=1, executing step 2.7, otherwise, let q=q+1, executing step 2.5;

step 2.7: initializing variable m=1, m is the initial abscissa value of pixel point in image inner circulation, and initializingWherein r is the distance between any two pixels in the initialized image, and the initial coordinates (a, b) are two random positive integers;

step 2.8: if m is less than or equal to W, the variable n=1, n is the initial value of the ordinate of the pixel point in the inner loop, executing the step 2.9, otherwise, executing the step 2.13;

step 2.9: if n is less than or equal to H, executing the step 2.10, otherwise, making m=m+1, and executing the step 2.8;

step 2.10: if the point (m, n) satisfies R ₂ Step 2.11 is performed, otherwise, let n=n+1, step 2.9 is performed;

step 2.11: calculating a distance dist ((p, q), (m, n)) between the point (p, q) and the point (m, n), wherein dist is a distance metric function, if r < dist ((p, q), (m, n)), r=dist ((p, q), (m, n)), while letting a=m, b=n;

step 2.12: let n=n+1, jump to step 2.9;

step 2.13: the points (p, q) take the points (a, b) as the centers, R is the symmetrical distance, a corresponding pixel point L in R is obtained, and the value of the pixel point L is assigned to the current point Z; wherein L is _x ＝2a-p，L _y =2b_q, where L _x Is the abscissa L of the pixel point L _y Is the ordinate of the pixel point L;

step 2.14: let q=q+1, jump to step 2.5;

step 3: for the parameter minimum scale sigma _min Maximum dimension sigma _max Setting constant coefficient gamma, constant coefficient c and scale step delta sigma, wherein sigma _min <σ _max ；

Step 4: for the initial dimension sigma _i Setting sigma _i ＝σ _min ；

Step 5: constructing a Hessian matrix of the image I, constructing a Gaussian function through an optimal response scale sigma, and subtracting the image I from the result of convolution operation of the image I and the Gaussian function to obtain an image M, wherein the method comprises the following specific steps of:

step 5.1: initializing a variable p=1, a scale variable k=1, and a variable v _max＝0, wherein v_max Is the maximum response value variable of the pixel;

step 5.2: if p is less than or equal to W, the variable q=1, the step 5.3 is executed, otherwise, the step 6 is executed;

step 5.3: if q is less than or equal to H, executing step 5.4, otherwise, making p=p+1, and executing step 5.2;

step 5.4: judging the current scale sigma _i Whether or not it is greater than the maximum dimension sigma _max If greater than sigma _max Let sigma _i ＝σ _min Executing the step 5.12, otherwise, executing the step 5.5;

step 5.5: constructing a Hessian matrix of points (p, q):

wherein ,

D＝∫∫I(p-u,q-v)×A(u,v)dudv

E＝∫∫I(p-u,q-v)×B(u,v)dudv

F＝∫∫I(p-u,q-v)×C(u,v)dudv

a (u, v), B (u, v) and C (u, v) are based on the scale sigma _i Construction of the second partial derivative of Gaussian

Wherein D, E, F is the convolution of image I with the corresponding gaussian second partial derivative A, B, C;

step 5.6: calculating two eigenvalues lambda of the Hessian matrix obtained in step 5.5 ₁ 、λ ₂, wherein

wherein

Step 5.7: calculating at the current scale sigma _i The response value of the pixel (p, q)

wherein

in the formula ,R_B Is the degree of anisotropy of the linear shape, S is the norm;

step 5.8: if it isMake->And κ=σ _i ；

Step 5.9: constructing a Gaussian kernel based on the scale kappa to obtain a Gaussian curve g _κ ；

Step 5.10: the Gaussian curve g obtained in the step 5.9 _κ Performing convolution operation with the image I, and subtracting the image I from the convolution result to obtain an image M;

M(p,q)＝(g _κ *I)(p,q)-I(p,q)

step 5.11: updating the scale sigma _i ＝σ _i And (5) jumping to the step 5.4, wherein the delta sigma is the scale step length;

step 5.12: updating the coordinates q=q+1, and jumping to step 5.3;

step 6, setting an initial zero level set segmentation curve on the image M, and generating a high-one-dimensional level set function phi according to the set zero level set;

the set initial zero level set segmentation curve is a randomly generated closed curve with a specified shape, and the level set function phi is a high one-dimensional symbol distance function generated according to the set initial level set segmentation curve;

step 7, establishing energy functional E (phi) of level set segmentation and setting a weighting coefficient alpha in the level set method ₁ Weighting coefficient alpha ₂ Curve smoothing term coefficient v and area smoothing term coefficient mu;

the energy functional is as follows:

E(φ)＝E _D +vL(φ)+μP(φ)

wherein ,E_D For image data items, L (φ) is a length constraint, P (φ) is an area constraint, and H (φ (P, q)) is a Heaviside function of φ (P, q);

E _D (φ)＝α ₁ ∫ _Ω |M(p,q)-c ₁ | ² H(φ(p,q))dpdq+α ₂ ∫ _Ω |M(p,q)-c ₂ | ² (1-H(φ(p,q)))dpdq

P(φ)＝∫H(φ(p,q))dpdq

wherein c₁ Is the gray average value in the zero level set segmentation curve, c ₂ Is the gray average value outside the zero level set segmentation curve;

step 8: setting a convergence threshold epsilon of the maximum iteration number J and the energy functional E (phi) of the image level set;

step 9: updating images separatelyLevel set function set H (phi) ^j )、1-H(φ ^j )、And->φ ^j Obtaining the energy functional E of the jth iteration ^j ；

Specifically:

where Δt is the time step size,a time partial derivative of the level set function that is iteratively updated, the time partial derivative formula being as follows:

step 10: judging whether the iteration termination condition is reached: if the current iteration number J reaches the maximum iteration number J or the difference value between the iteration results of two adjacent times of the energy functional E (phi) is smaller than a convergence threshold epsilon set in advance, executing the step 11, otherwise, adding 1 to the iteration number, and returning to the step 9;

step 11: from a current updated image level set function phi ^j An image, i.e. a segmentation result of the image, is constructed.