Shape inhibition based fundus microstructure detection method

Abstract

The application provides a fundus small target detection method based on shape inhibition. Firstly, a spot enhancement and elongation inhibition filter and an additive bias correction level set are used to extract preliminary candidate targets. The spot enhancement and elongation inhibition filter is used to enhance the spot structure and inhibit the elongation structure. The additive bias correction level set is used to converge to small targets with less pixels by minimizing an energy function after additive bias correction. Then, threshold constraint is performed on the candidate targets to complete the preliminary screening. Secondly, twelve GLCM features are extracted, including features that have not been used in existing methods and have high sensitivity. Classification is performed according to the extracted features and the position information of the candidate targets. Finally, a random oversampling enhancement (RUSBoosting) classifier is used for binary classification in the classification stage.

Description

TECHNICAL FIELD

[0001] The application relates to a fundus small target detection method based on shape inhibition, and belongs to the fields of computer vision and medical image processing. BACKGROUND

[0002] Firstly, small targets are difficult to detect due to their small structure, thereby reducing the precision of small target detection. Secondly, the contrast of small targets is low, and morphological methods and filter-based methods can cause the targets to be missed or even unable to be detected. In addition, a retinal image has excessive noise and a complex source. The noise has similar characteristics of small targets, and it is difficult to eliminate the noise, and the small targets are also affected. SUMMARY

[0003] The application aims to provide a fundus small target detection method based on shape inhibition, so as to realize initial lesion judgment, reduce the burden and pressure of doctors, and reduce the loss of the government.

[0004] In order to achieve the above purpose, the application comprises the following steps:

[0005] Step 1: candidate target extraction: (1) taking the green channel G of a color fundus image as an input image, removing fundus structures with a strip structure by using a spot enhancement and elongation suppression filter to obtain a response graph F1, and selecting a threshold T1 to binarize F1 to obtain a binary graph B1;

[0006] (2) performing additive bias correction on the green channel G to obtain a result graph P, and then finding a circular candidate target region by using an additive bias correction level set, and the specific steps are as follows:

[0007] 1) first, define a circular neighborhood (or rectangular neighborhood) with a radius of p, which takes O y = [x: |x-y| <= p} as the center. The division of the entire domain Omega results in the division of the neighborhood O y , that is the division of O y . For a slowly changing bias field b, the intensity of each partition O y ∩ Omega i is replaced by b(x), where x belongs to O y ∩ Omega i . If we divide the entire domain Omega into two sub-regions: Omega1 and Omega2, the boundary line of Omega1 and Omega2 can be represented by a curve C. The local intensity can be classified by using a K-means clustering algorithm, which can be represented in a continuous form as

[0008]

[0009] where I b is the image after additive bias correction, I b b i (x) represents the bias field of the image O y after additive bias correction, the intensity of the neighborhood O y can be divided into N clusters, with C i ≈b i , i = 1,..., N, u i is a membership function related to Ω i , with the constraint: u i (x) = 1, x ∈ Ω i and u i (x) = 0, F y can be rewritten as:

[0010]

[0011] Due to the characteristics of σ in the calculation, only the data of the small area close to the narrow band is used to calculate the curve C. In order to solve this problem, a truncation function G(y-x) is defined to introduce sub-regions Ω1 and Ω2, and the clustering function can be rewritten as:

[0012]

[0013]

[0014] The purpose of adding G σ (y-x) is to merge O y into the region Ω i , and the minimization E y is realized by integrating y over the image domain Ω, and the function is as follows:

[0015]

[0016] The disjoint regions Ω1 and Ω2 can be represented by and , where represents the Heaviside function, and the level set formula is as follows:

[0017]

[0018] The variational level set is obtained by minimizing the energy function E using the standard gradient descent method, and then the de-parameterization regularization function and the neighborhood filtering are added to replace the regularization term and the length term in the traditional level set method, and then the response graph F2 is obtained. Selecting a threshold T2 to binarize F2 to obtain a binary graph B2;

[0019] (3) performing logical operation on the two binary images in (1) and (2) to obtain a candidate target;

[0020] (4) performing threshold segmentation on the candidate target by using area (Area) and metric;

[0021] Step 2: feature extraction, extracting a gray level co-occurrence matrix (GLCM) feature vector with a dimension of 1x6, i.e. maximum probability, autocorrelation, cluster prominence, cluster shadow, differential entropy, differential variance and energy, from each candidate target region R in each image i , i = 1, 2..., N, where N is the total number of candidate targets in the image;

[0022] Step 3: classification, inputting the candidate target and the feature vector into an AdaBoost classifier for classification.

[0023] The blobness enhancement and elongation suppression filter described in (1) in Step 1, introducing a suppression factor (λ2-λ1), calculating a three-dimensional Hessian matrix, and constructing an enhancement function VR, wherein λ1, λ2 represent eigenvalues calculated from the Hessian matrix, and λ3 = λ2, the formula is as follows:

[0024]

[0025] According to the Hessian eigenvalue relationship satisfied by the blood vessels: λ2≈λ3Λ|λ 2，3 |>>|λ1|, wherein Λ represents a union, and formula (7) is obtained according to λ2>>|λ1|:

[0026]

[0027] Then, according to λ2≈λ3, the following is derived:

[0028]

[0029] The Hessian eigenvalue relationship based on the blob shape structure is λ1≈λ2≈λ3Λ|λ 1，2，3 |, wherein Λ represents a union, and formula (7) can be simplified as:

[0030]

[0031] Regularizing the eigenvalue λ3, the blobness function B is obtained:

[0032]

[0033] Then, the response map thereof is obtained.

[0034] The beneficial effects of this invention are: the method is very effective for detecting small targets in the fundus, and it can detect most small targets and has high sensitivity on five public datasets. Brief description of the attached figures

[0035] Figure 1 This is a flowchart of the fundus small target detection method based on shape suppression of the present invention. Detailed Implementation

[0036] Candidate target extraction: The green channel image of the original image is used as input and processed by a speckle enhancement and elongation suppression filter. The eigenvalues ​​are calculated based on the Hessian matrix of the image, as shown in the following formula:

[0037]

[0038] Where G(x, s) = (2πs) 2 ) -3 / 2 exp(-x T x / 2s 2 ) is a 3D Gaussian, and * indicates convolution.

[0039] For each x, the eigenvalues ​​eigH(x, s) → λ are obtained through eigenvalue decomposition. i , i = 1, 2, 3.

[0040] The blobness function represents the degree of isotropy of the blob shape, taking into account some measures used to describe the structural isotropy and anisotropy of the diffusion tensor, such as volume ratio, which is used to detect approximately spherical diffusion tensors and is not proportional to the eigenvalues.

[0041]

[0042] First, replace λ1 in the denominator of the existing spot enhancement function with (λ2-λ1), then keep λ1 in the numerator unchanged. This enhances the spot structure while suppressing elongated structures. The formula is as follows:

[0043]

[0044] According to the Hessian eigenvalue relationship satisfied by blood vessels: λ2≈λ3Λ|λ 2，3 |>>|λ1|,Formula (14) is first rewritten based on λ2>>|λ1| as:

[0045]

[0046] Then, based on λ2≈λ3, equation (15) can be rewritten as:

[0047]

[0048] simplification,

[0049] According to |λ3|»|λ1|, the formula (14) can be simplified as: approximately equal to 0. The blobness function exhibits eigenvalues with high isotropic magnitude (|λ1|≈|λ2|≈|λ3|»|λ 1，2，3 |). The formula (14) can be simplified as:

[0050]

[0051] approximately equal to 0. The blobness function exhibits eigenvalues with high isotropic magnitude (|λ1|≈|λ2|≈|λ3|»|λ 1，2，3 |). The formula (14) can be simplified as:

[0052] The uniform intensity image region The VR in the formula (14) is too sensitive to noise, and the magnitude of λ3 is low. In order to improve the performance of λ3, it is suggested to regularize the eigenvalue λ3 in the formula (14) as:

[0053]

[0054] where s denotes at each scale, x denotes at each voxel, and τ is a cutoff threshold between 0 and 1, which can be chosen by the user. The final blobness function is:

[0055]

[0056] For structures with low contrast and affected by image noise, the output value of B in the formula (19) will decrease, and the blobness function cannot enhance the structure, so B only enhances structures with contrast higher than the image noise level. When regularizing λ3, a higher τ value will increase the difference between λ1 and λ ρ 3 of low-contrast structures, thereby enhancing low-contrast structures. In the fundus image, small targets have low contrast compared to the surrounding background and are affected by image noise, and the magnitudes of the eigenvalues λ1 and λ3 are both low, so a higher τ value can enhance small targets. Therefore, we set the τ value to 1 to sufficiently enhance low-contrast structures. After the above operations, the binary image is obtained.

[0057] Then the steps of the level set with additive bias correction are as follows: first, define a circular neighborhood (or rectangular neighborhood) with a radius of ρ, which is centered at each point y∈Ω defined by y = {x: |x-y|≤ρ}. The division of the entire domain Ω leads to the division of the neighborhood O y , i.e. forming the division of O y . For a slowly varying bias field b, the intensity of each partition O y ∩Ω i is replaced by b(x), where x∈O y ∩Ωi If we divide the entire domain Ω into two sub-regions: Ω1 and Ω2, the boundary between Ω1 and Ω2 can be represented by curve C. Neighborhood O y The intensity can be divided into N clusters, with C as the center. i ≈b i Let i = 1, ..., N. The K-means clustering algorithm is applied to classify local intensities. The K-means algorithm is an iterative process and can be expressed in continuous form as follows:

[0058]

[0059] Among them I b For the image after additive bias correction, I b b i (x) represents the bias field of the image after additive bias correction, u i It is a Ω i The related member functions have the following membership degree: u i (x) = 1, x ∈ Ω i and u i (x) = 0, F y It can be rewritten as:

[0060]

[0061] The above truncated Gaussian function G σ Introducing (yx), the clustering function can be rewritten as:

[0062]

[0063] Join G σ The purpose of (yx) is to make O y Merge into Ω i Region. Minimize E by integrating y over the image domain Ω. y Its function is as follows:

[0064]

[0065] The non-intersecting regions Ω1 and Ω2 can be used and The membership function is defined as follows. The level set formula is as follows:

[0066]

[0067] Analyzing the local region problem in equation (24), its data items are extracted as follows:

[0068] I b ei (x) = ∫G σ (y - x) | I b (x) - r(y) - I b b i (y) | 2 dy (25)

[0069] The energy function E is minimized by the standard gradient descent method, while the activation function is introduced to balance the data term, which is an odd function. The variational level set formulation is as follows:

[0070]

[0071] where η is a fixed value of 7, and the image standard deviation β is added to the denominator of the data-driven term. The smaller the value of β, the closer the distance between the target and the background, and the smaller the data-driven difference. Adding the coefficient α can improve the robustness of equation (24), because it does not need to adjust the parameters repeatedly for pictures of different sizes and features.

[0072] Subsequently, the de-parameterization regularization function and neighborhood filtering are added to replace the regularization term and length term in the traditional level set method, which are defined as shown in equations (27) and (28):

[0073]

[0074]

[0075] The setting of the window size k: for images with large targets and clear boundaries, k can be increased to contain large targets. It is known that otherwise it will be reduced. The setting of the coefficient α: for targets with weak intensity, the value of α should be increased. At the same time, increasing α is beneficial to speed up the iteration speed. However, for targets with fuzzy boundaries, the value of α should be reduced to prevent the target boundary from being damaged. The setting of the Gaussian scaling parameter σ: for images with uneven lighting or serious background interference, the value of σ should be moderately increased. After the above processing, a binary image is obtained.

[0076] The two binary images are subjected to an or operation, and then a suitable candidate selection algorithm is proposed, which reduces the number of analysis targets to reduce the false detection rate, algorithm complexity and computation time.

[0077] This step processes the blob-enhanced, elongated-inhibited filtered image and performs connected component analysis on the binary image to identify and select objects with an internal satisfying area less than 150px and a circularity of 0.5. Since the experiments were done on different datasets, their parameters were chosen to vary according to the dataset. At the same time, the image after level set processing with additive bias correction is closed to extract linear structures and then removed. The physical result of the closing operation will smooth the outline of the object and connect narrow gaps to form elongated curves. There are still a lot of noise and small targets after closing. Through analysis, we found that blood vessels are linear, and noise and small targets are not linear, so we use morphological analysis to extract linear blood vessels. Then we remove the extracted containers.

[0078] Feature extraction: Feature extraction is the process of extracting some features from the pre-processed image. Gray level co-occurrence matrix (GLCM) focuses on texture information. GLCM uses the spatial relationship between pixel features in the neighborhood to capture numerical feature values. It is a matrix with row and column equal to the number of gray levels. We extracted a total of twelve texture features according to GLCM. Including maximum probability, autocorrelation, cluster prominence, cluster shadow, contrast, correlation, difference entropy, difference variance, dissimilarity, energy, entropy. The features not used by the existing method and similar to the task of detecting small targets in fundus images are also detecting small targets in complex background.

[0079] Classification: To solve the imbalance between positive and negative samples, a random oversampling enhanced classifier is used for binary classification.

[0080] It should be noted that the above only describes the embodiments of the present application, and only explains the present application, and does not limit the scope of the present application. Changes that are only obvious to those skilled in the art and belong to the technical concept of the present application are also within the scope of protection of the present application.

Claims

1. A shape-based inhibition-based fundus small target detection method, comprising the following four steps: Step 1: candidate target extraction: (1) take the green channel G of the color fundus image as the input image, remove the fundus structure with bar structure by using the bar inhibition-based blob enhancement filter, introduce the inhibition factor (λ2-λ1), calculate the three-dimensional Hessian matrix, construct the enhancement function VR, wherein λ1, λ2 represent the eigenvalues calculated by the Hessian matrix, λ3=λ2, the formula is as follows: According to the Hessian eigenvalue relation satisfied by the blood vessels: λ2≈λ3Λ|λ 2，3 |>|λ1|, where Λ denotes the union, formula (7) first obtains according to λ2>>|λ1|: Then according to λ2≈λ3, it is deduced that: The Hessian eigenvalue relation based on the spot shape structure is λ1≈λ2≈λ3Λ|λ 1，2，3 , wherein Λ represents a union, and formula (7) is simplified as: The regularized eigenvalue λ3: The blobness function B is obtained: The response map F1 is obtained, and the threshold T1 is selected to binarize F1 to obtain the binary image B1; (2) Perform additive bias correction on the green channel G to obtain the result image P, and then use the additive bias correction level set to find the circular candidate target area, the specific steps are as follows: 1) First, define a circular neighborhood of radius p centered at each point y e Ω defined by y = {x: |x - y| < p} and divide the entire domain Ω results in a division of the neighborhood O y , i.e. forms a division of O y ; for a slowly varying bias field b, the intensity of each partition O y ∩ Ω i is replaced by b(x) where x e O y ∩ Ω i ; divide the entire domain Ω into two sub-regions: Ω1 and Ω2, with the boundary between Ω1 and Ω2 represented by a curve C; apply the K-means clustering algorithm to classify the local intensities, expressed in continuous form as where I b is the image after additive bias correction, denotes the O y of the image after additive bias correction, neighborhood O y is divided into N clusters with C i ≈b i , i = 1,..., N, u i is a membership function related to Ω i with constraints: u i (x) = 1, x e Ω i and u i (x) = 0, F y is rewritten as: Due to the characteristics of σ in calculation, only the data of the small area close to the narrow band is used to calculate the curve C, in order to solve this problem, a truncation function G(y-x) is defined to introduce sub-regions Ω1 and Ω2, and the clustering function is rewritten as: Join G σ The purpose of (y - x) is to add O y Merging into Ω i Region, the minimization of E y The function is as follows: The disjoint regions Ω1 and Ω2 are denoted by and where denotes the Heaviside function, and the level set formulation is as follows: The variational level set is obtained by minimizing the energy function E using the standard gradient descent method, then the deparameterization regularization function and neighborhood filtering are added to replace the regularization term and length term in the traditional level set method, and then the response map F2 is obtained, and the threshold T2 is selected to binarize F2 to obtain the binary image B2; (3) Perform logical OR operation on the two binary images in (1) and (2) to obtain the candidate target; (4) Threshold segmentation of the candidate target by area Area and circularity metric; Process the image after blob enhancement and elongation inhibition filtering, and perform connected component analysis on the binary image to identify and select objects with an internal area less than 150px and a circularity of 0.5, perform closing operation on the image after additive bias correction level set processing, extract linear structures, and then remove them, the physical result of the closing operation will smooth the outline of the object and connect the narrow gaps to form elongated curves, there are still many noises and small targets after the closing operation, through analysis, blood vessels are linear, noises and small targets are not linear, so morphological analysis is used to extract linear blood vessels, and then the extracted container is removed; Step 2: Feature extraction, for each candidate target region R in each image i , i = 1, 2..., N, where N is the total number of candidate targets in the image, a 1x6-dimensional gray-level co-occurrence matrix (GLCM) feature vector is extracted, i.e., maximum probability, autocorrelation, cluster prominence, cluster shadow, contrast, correlation, difference entropy, difference variance, dissimilarity, energy, entropy; Step 3: classification, input the candidate target and feature vector into the AdaBoost classifier for classification.

Images