Automatic quantification method for retinal vessel diameter of fundus image
Technical Field
The invention belongs to the technical field of medical image processing, and particularly relates to an automatic quantification method for retinal vessel diameter of fundus images.
Background
Currently, the current state of the art commonly used in the industry is as follows: fundus images of retinal blood vessels are an economical and effective clinical aid for ophthalmic diseases, and their development has gained widespread attention. Changes in retinal arteriovenous intersections, arterial branches, venous walking paths, and morphology are important symptoms of arteriosclerotic lesions of the retina. When observing fundus images, doctors mainly adopt a shooting method, a projection method, an ophthalmoscope image-measuring method, a scaling observation method and the like to measure the retinal vascular tube diameter. Most of the methods are judged according to manual observation, manual measurement and experience, and retinal blood vessels can not be technically and accurately detected and extracted, and the sizes of the blood vessel diameters are quantized, so that early prediction and auxiliary diagnosis of cardiovascular and cerebrovascular diseases are not facilitated, and patients do not have clear knowledge of own illness state. Therefore, the manual analysis of fundus images has the defects of strong subjectivity, low diagnosis efficiency, poor repeatability, long time consumption, high labor intensity, difficult guarantee of accuracy and consistency, and the like. The invention utilizes the digital image processing technology to measure the retinal vascular diameter, overcomes the defects of large artificial subjective error and poor repeatability, and has the advantages of objective measurement data, repeated measurement and labor cost saving.
Disclosure of Invention
Aiming at the problems existing in the prior art, the invention provides an automatic quantification method for the retinal vessel diameter of fundus images.
The invention discloses an automatic quantification method of retinal vascular tube diameter of fundus images, which comprises the following steps:
step one, detecting the video disc by using a two-step positioning method combining thickness positioning and fine positioning.
The optic disc positioning is the basis of the operations of fundus image retinal vessel tracking, macula lutea and lesion extraction, and the like, and various parameters of the optic disc have great significance for early prediction and clinical auxiliary diagnosis. The shape of the edge of the optic disc is between an ellipse and a circle, and the rough positioning regards the shape as a circle, and the fine positioning detects the real shape.
(1) Optic disc coarse positioning
Binarization is carried out according to the gray level distribution condition of the fundus image, mathematical morphology knowledge is adopted, partial interference is removed, a video disc interval is obtained, the mass center of the video disc is roughly positioned, and the video disc is roughly positioned into a round shape. Assuming that the mass center of the optic disc is O, finding the upper and lower points A, B and the left and right points C, D of the edge of the optic disc along the horizontal and vertical directions of the O point respectively, finding the mass center by using a least square fitting circle, then obtaining an image of the optic disc,
the coordinates of centroid O (x o ,y o ) The calculation is carried out as follows,
wherein xi and yi The abscissa and ordinate of the optic disc edge points, respectively, and N is the total number of optic disc edge points.
Under the principle of minimum error square sum, the least square estimation of the circle parameters is as follows:
X=(U T U) -1 U T V=U -1 V (2)
wherein
(x 1 ,y 1 ),(x 2 ,y 2 ),(x 3 ,y 3 ) Is any of A, B, C, DAnd 3-point coordinates, and taking the center coordinate mean value under 4 conditions as the center of mass coordinates of the video disc respectively.
(2) Precision positioning of optic disk
The sub-image can be cut by roughly locating the video disc by using the video disc brightness information and the local blood vessel information to roughly determine the position of the video disc. The optimal video disc boundary curve can be conveniently obtained by adopting a dynamic contour tracking method based on a non-reinitialized variation level set, and the video disc boundary curve has good fault tolerance and robustness and realizes the accurate positioning of the video disc.
The image energy function may be expressed as
In the omega table image field,
as a function of symbol distance.
Where d is the shortest distance of point (x, y) to the curve,/and>
table penalty term, which ensures that the level set function is a signed distance function, μ is the penalty term weight, constant coefficients λ and v, λ>0, delta (■) is the Dirac function, stop function +.>
H (■) is a Heaviside function.
And step two, preprocessing the fundus image, denoising and enhancing, separating retinal blood vessels from the fundus image background by adopting a multi-scale method, and removing pseudo blood vessels existing in segmentation.
The distortion phenomenon of the individual image occurs, and the distortion correction is realized by adopting a linear interpolation method and a normalized image processing method. In the fundus image imaging process, the problems of noise, uneven illumination, low image contrast and the like often exist in the acquired images due to the influence of non-ideal acquisition conditions, characteristic differences of an imaging sensor and the like. Thus, the fundus image is preprocessed prior to retinal vessel segmentation. The preprocessing mainly comprises image denoising, contrast enhancement, mathematical morphology processing and the like.
(1) Fundus image denoising and contrast enhancement based on Contourlet transform. Noise is suppressed in a Contourlet domain by adopting a soft threshold denoising method; and (3) carrying out enhancement processing on the transformed sub-band coefficients in the band-pass direction by adopting a nonlinear enhancement operator, and obtaining a contrast enhancement image through inverse transformation.
The nonlinear enhancement operator is as follows:
E(x i,j )=x i,j .sign(x i,j ).tanh(b.u).(1+c.exp(-u.u) (4)
(2) Retinal vessel segmentation based on multi-scale analysis linear tracking
The process of dividing the retinal blood vessel by using the multi-scale analysis method is a process of extracting blood vessels with different areas and widths under different scales. Firstly, selecting a part of seed pixels from an image, starting tracking, endowing higher confidence coefficient to pixel points meeting the condition in the tracking process, and quantifying confidence coefficient matrixes of all pixel points after the tracking is finished, so that an initial blood vessel is obtained.
(1) Linear tracking initial seed selection
The vascular path following the initial seed is defined as:
V s ={(x,y:T low <I(x,y)<T high ) (5)
V s for the selected seed sequence, I (x, y) represents the gray level of the pixel of the image at (x, y), T low and Thigh Representing minimum and maximum gray values of an image, estimated by the size of a vascular region in the image, pixels below T low The points of (2) belong to the darker areas of the retinal image blood vessels, and the pixels are higher than T high Is of the genusIn areas where the retinal image vessels are brighter, the confidence of these pixels is higher.
In the linear tracking process, the confidence of the pixels belonging to the blood vessel is stored in the sequence C w If the confidence of a certain pixel in the confidence matrix is higher, the probability that the pixel belongs to a blood vessel is relatively higher. Initially, elements at all scales in the confidence matrix are set to 0.
C w (x,y):=0 (6)
(2) Initialization of linear tracking
The initialization process of linear tracking of retinal blood vessels can be represented by equation (7):
k:=1,V c (k):=V s (t),C c :={} (7)
wherein ,Vc C, tracking the obtained pixel sequence under the current cyclic variable t c For the new tracking pixel sequence.
(3) Estimation of new tracking pixels
The pixel point tracked at present is the last to enter V c Object C of (2) c Is a sequence formed by candidate pixel points, namely the nearest 8 pixel points around the current tracking pixel, excluding the V c As shown in formula (8):
C c =N 8 (V c (k))-V c (8)
adding the current tracking pixel to V c As shown in formula (9):
k:=k+1,V c (k):=(x,y) (9)
confidence matrix C w Can be updated as:
C w (x,y):=C w (x,y)+1,(x,y):=(x c ,y c ) (10)
when all the parameter values in the formula (9) are smaller than T, searching for the next seed pixel point (T: =t+1) is started.
(4) Multi-scale linear tracking
And repeatedly tracking all seed points according to a certain scale, wherein the scale number is determined by the width of a blood vessel to be detected in the retina image, and if the blood vessel in the retina image has large variability, the corresponding scale number is also large. Through multiple simulation experiments, the effect of the dimension w= 3,5,7,9,11 is better.
(5) Initial extraction of retinal blood vessels
And (5) carrying out initial estimation on the blood vessel by adopting a mapping quantification method. The initial vessel is determined by confidence matrix greater than T C Is usually T C The value is equal to the value of the dimension W.
And thirdly, measuring the vessel diameter by adopting a boundary method, and displaying the result on a system interface.
(1) Vessel measurement image region selection
From the segmented retinal blood vessel images, the blood vessel images of the region of interest are selected for edge detection, and then the blood vessels of the sub-images are subjected to width measurement, so that the processing time can be reduced.
(2) Calculation of vessel diameter width
The blood vessel direction in the sub-image is parallel to the horizontal direction as much as possible, and if not horizontal, a rotation angle method can be adopted to make the blood vessel direction parallel to the horizontal direction as much as possible. The automatic measurement of the retinal vessel diameter width is realized by adopting a minimum distance method.
For the coordinates of the upper and lower edge points in the vertical direction, the distance value is calculated by using a distance formula between the two points, and the calculated distance value is taken as a blood vessel width value and can be expressed as follows:
wherein ,wi Is the width value of column i, (x) i ,y i) and (xk ,y k ) The coordinates of the upper and lower edges of the vessel, respectively. The distance method has the advantages of wide application range, good stability and high precision.
(3) Scale positioning and pixel pitch calculation
To convert the vessel diameter, which is calibrated in pixels, to the actual size of the vessel diameter, the pixel spacing is measured, and therefore the scale must be positioned in the image.
And calculating the physical length of the pixels according to the minimum scale value and the number of the pixels between the adjacent scale lines. The scale is set in the image with the same size as the truncated sub-image, for example, the scale is 10mm×10mm, the edge detection processing is performed on the scale image, a binary image is obtained, the pixel value of the scale is set to be 1, and the pixel value of the background is set to be 0. And traversing the binary image from top to bottom for the binarized scale image, and calculating the pixel number of the binary image.
If the fixed scale is 10mm in size and the number of pixels between scales is Num, the physical length p of the pixels can be determined by
The unit of p is determined in mm/pixel.
(4) Converting pixel size of vessel diameter into physical size
The physical size of the retinal vessel diameter width can be known by multiplying p by the pixel measured at each marker point.
Another object of the present invention is to provide an automatic quantification system for implementing the fundus image retinal vessel caliber measurement method, the system comprising:
the eyeground image optic disk detection module is used for realizing the detection of the eyeground image optic disk;
the preprocessing and blood vessel segmentation module is used for realizing the denoising and contrast enhancement of fundus images and the segmentation of retinal blood vessels;
the blood vessel width measurement and application module is used for realizing retinal blood vessel diameter measurement, early prediction of clinicians and auxiliary diagnosis.
In summary, the invention has the advantages and positive effects that: on the basis of detecting the eyeground image optic disc, preprocessing is carried out on the image, so that the contrast ratio of the eyeground image is improved. The retinal blood vessel is segmented by adopting a multiscale analysis method, the size of the pipe diameter is measured by adopting a boundary method, the result is displayed on a system interface, and references are provided for early prediction and auxiliary diagnosis of clinicians, so that labor is saved, and the efficiency is improved.
Drawings
Fig. 1 is a flowchart of a method for automatically quantifying retinal vascular tube diameters in fundus images provided by an embodiment of the invention.
Fig. 2 is a schematic structural diagram of an automatic quantification system for retinal vascular diameter of fundus images provided by an embodiment of the present invention;
in the figure: 1. an image preprocessing module; 2. a detection and segmentation module; 3. and the quantization and application module.
Fig. 3 is a flowchart of an implementation method for automatically quantifying retinal vascular tube diameters of fundus images according to an embodiment of the present invention.
Fig. 4 is a schematic diagram of the rough positioning process according to the embodiment of the present invention.
In the figure: a. the video disc after interference removal; b. centroid positioning schematic diagram; c. a tray image; d. and (5) roughly positioning the result.
Fig. 5 is a schematic diagram of a fine positioning effect according to an embodiment of the present invention.
In the figure: a. dynamically tracking a process diagram; b. tracking a result graph; c. thickening a tracking curve; d. and (5) fine positioning results.
Fig. 6 is a schematic diagram of preprocessing and retinal vessel segmentation effects provided by an embodiment of the present invention.
In the figure: a. an original image; b. a contrast enhanced image; c. initial segmentation of blood vessels; d. the blood vessel after the interference is removed.
Fig. 7 is a schematic diagram of a measurement effect of a retinal vascular width measurement procedure according to an embodiment of the present invention.
In the figure: a. intercepting a blood vessel; b. edge detection; c. removing isolated noise; d. rotating the blood vessel; e. establishing a scale; f. and measuring the result.
Detailed Description
The present invention will be described in further detail with reference to the following examples in order to make the objects, technical solutions and advantages of the present invention more apparent. It should be understood that the specific embodiments described herein are for purposes of illustration only and are not intended to limit the scope of the invention.
The method aims at solving the problems of high labor intensity and long time consumption caused by manual observation and manual measurement of the vessel diameter in the prior art. The invention automatically measures the retinal vascular caliber through an image processing technology, provides reference for early prediction and auxiliary diagnosis for doctors, reduces the labor intensity of the doctors and improves the diagnosis efficiency.
The principle of application of the invention is described in detail below with reference to the accompanying drawings.
As shown in fig. 1, the automatic quantification method for retinal vascular tube diameters of fundus images provided by the embodiment of the invention comprises the following steps:
s101: automatically detecting the video disc by using a two-step positioning method combining thickness positioning and fine positioning;
s102: preprocessing an eye fundus image, separating retinal blood vessels from a fundus image background, and removing pseudo blood vessels existing in segmentation;
s103: and measuring the vessel diameter by adopting a boundary method, and displaying the result on a system interface.
In the step S101, the automatic detection of disc positioning in the optic disc is the basis of the operations of fundus image retinal vessel tracking, macula lutea and lesion extraction, etc., and various parameters of the optic disc have very important significance for clinical diagnosis. The shape of the edge of the optic disc is between an ellipse and a circle, and the rough positioning regards the shape as a circle, and the fine positioning detects the real shape.
(1) Optic disc coarse positioning
Binarizing the image according to the gray level distribution condition of the fundus image, adopting mathematical morphology knowledge to remove partial interference, obtaining a video disc interval, roughly positioning the mass center of the video disc, and roughly positioning the video disc into a circle. Assuming that the mass center of the optic disc is O, finding the upper and lower points A, B and the left and right points C, D of the edge of the optic disc along the horizontal and vertical directions of the O point respectively, finding the mass center by using a least square fitting circle, then obtaining an image of the optic disc,
the coordinates of centroid O (x o ,y o ) The calculation is carried out as follows,
wherein xi and yi Respectively the abscissa and the ordinate of the video disc edge point, and N is the video disc edge pointIs a sum of (3).
Under the principle of minimum error square sum, the least square estimation of the circle parameters is as follows:
X=(U T U) -1 U T V=U -1 V (2)
wherein
(x 1 ,y 1 ),(x 2 ,y 2 ),(x 3 ,y 3 ) Is the coordinate of any 3 points in A, B, C, D, and the center coordinate mean value of the 4 conditions is respectively used as the barycenter coordinate of the video disc.
(2) Precision positioning of optic disk
The sub-image can be cut by roughly locating the video disc by using the video disc brightness information and the local blood vessel information to roughly determine the position of the video disc. The optimal video disc boundary curve can be conveniently obtained by adopting a dynamic contour tracking method based on a non-reinitialized variation level set, and the video disc boundary curve has good fault tolerance and robustness and realizes the accurate positioning of the video disc.
The image energy function may be expressed as
In the omega table image field,
as a function of symbol distance.
Where d is the shortest distance of point (x, y) to the curve,/and>
table penalty term, which ensures that the level set function is a signed distance function, μ is the penalty term weight, constant coefficients λ and v, λ>0, delta (■) is Dirac function, stop function->
H (■) is a Heaviside function.
In the step S102, distortion occurs in the individual image, and distortion correction is implemented by using a linear interpolation method and a normalized image processing method. In the fundus image imaging process, the problems of noise, uneven illumination, low image contrast and the like often exist in the acquired images due to the influence of non-ideal acquisition conditions, characteristic differences of an imaging sensor and the like. Thus, the fundus image is preprocessed prior to retinal vessel segmentation. The preprocessing mainly comprises image denoising, contrast enhancement, mathematical morphology processing and the like.
(1) Fundus image denoising and contrast enhancement based on Contourlet transform. Noise is suppressed in a Contourlet domain by adopting a soft threshold denoising method; and (3) carrying out enhancement processing on the transformed sub-band coefficients in the band-pass direction by adopting a nonlinear enhancement operator, and obtaining a contrast enhancement image through inverse transformation.
The nonlinear enhancement operator is as follows:
E(x i,j )=x i,j .sign(x i,j ).tanh(b.u).(1+c.exp(-u.u) (4)
(2) Retinal vessel segmentation based on multi-scale analysis linear tracking
The process of dividing the retinal blood vessel by using the multi-scale analysis method is a process of extracting blood vessels with different areas and widths under different scales. Firstly, selecting a part of seed pixels from an image, starting tracking, giving higher confidence to pixel points meeting the conditions in the tracking process, and quantifying confidence matrixes of all the pixel points after the tracking is finished, so that an initial blood vessel is obtained.
(1) Linear tracking initial seed selection
The vascular path following the initial seed is defined as:
V s ={(x,y:T low <I(x,y)<T high ) (5)
vs is the selected seed sequence, I (x, y) represents the gray level of the pixel of the image at (x, y), T low and Thigh Is the minimum gray level and the highest gray level of the image, is estimated by the size of the blood vessel region in the image, and the pixel is lower than T low The points of (2) belong to the darker areas of the retinal image blood vessels, and the pixels are higher than T high Which belongs to a bright area of the retinal image blood vessels, the confidence of these pixel points is high.
In the linear tracking process, the confidence of the pixels belonging to the blood vessel is stored in the sequence C w If the confidence of a certain pixel in the confidence matrix is higher, the probability that the pixel belongs to a blood vessel is relatively higher. Initially, elements at all scales in the confidence matrix are set to 0.
C w (x,y):=0 (6)
(2) Initialization of linear tracking
The initialization process of linear tracking of retinal blood vessels can be represented by equation (7):
k:=1,V c (k):=V s (t),C c :={} (7)
wherein ,Vc C, tracking the obtained pixel sequence under the current cyclic variable t c For the new tracking pixel sequence.
(3) Estimation of new tracking pixels
The pixel point tracked at present is the last to enter V c Object C of (2) c Is a sequence formed by candidate pixel points, namely the nearest 8 pixel points around the current tracking pixel, excluding the V c As shown in formula (8):
C c =N 8 (V c (k))-V c (8)
adding the current tracking pixel to V c As shown in formula (9):
k:=k+1,V c (k):=(x,y) (9)
confidence matrix C w Can be updated as:
C w (x,y):=C w (x,y)+1,(x,y):=(x c ,y c ) (10)
when all the parameter values in the formula (9) are smaller than T, searching for the next seed pixel point (T: =t+1) is started.
(4) Multi-scale linear tracking
And repeatedly tracking all seed points according to a certain scale, wherein the scale number is determined by the width of a blood vessel to be detected in the retina image, and if the blood vessel in the retina image has larger variability, the corresponding scale number is also larger. Through multiple simulation experiments, the effect of the dimension w= 3,5,7,9,11 is better.
(5) Initial extraction of retinal blood vessels
And (5) carrying out initial estimation on the blood vessel by adopting a mapping quantification method. The initial vessel is determined by confidence matrix greater than T C Is usually T C Is equal to the value of the dimension W.
In the step S103, the specific process of measuring the vessel diameter by the boundary method is as follows:
(1) Vessel measurement image region selection
From the segmented retinal blood vessel images, the blood vessel images of the region of interest are selected for edge detection, and then the blood vessels of the sub-images are subjected to width measurement, so that the processing time can be reduced.
(2) Calculation of vessel diameter width
The blood vessel direction in the sub-image is parallel to the horizontal direction as much as possible, and if the blood vessel direction is not horizontal, a rotation angle method can be adopted to enable the blood vessel direction to be parallel to the horizontal direction as much as possible. The automatic measurement of the retinal vessel diameter width is realized by adopting a minimum distance method.
For the coordinates of the upper and lower edge points in the vertical direction, the distance value is calculated by using a distance formula between the two points, and the calculated distance value is taken as a blood vessel width value and can be expressed as follows:
wherein ,wi Is the width value of column i, (x) i ,y i) and (xj ,y j ) The coordinates of the upper and lower edges of the vessel, respectively. The distance method has the advantages of wide application range, good stability and high precision.
(3) Scale positioning and pixel pitch calculation
To convert the vessel diameter, which is calibrated in pixels, to the actual size of the vessel diameter, the pixel spacing is measured, and therefore the scale must be positioned in the image.
And calculating the physical length of the pixels according to the minimum scale value and the number of the pixels between the adjacent scale lines. The scale is set in the image with the same size as the truncated sub-image, the scale is 10mm×10mm, the edge detection processing is carried out on the scale image, a binary image is obtained, the pixel value of the scale is set to be 1, and the pixel value of the background is set to be 0. And traversing the binary image from top to bottom for the binarized scale image, and calculating the pixel number of the binary image.
If the fixed scale is 10mm in size and the number of pixels between scales is Num, the physical length p of the pixels can be determined by
The unit of p is determined in mm/pixel.
(4) Converting pixel size of vessel diameter into physical size
The physical size of the retinal vessel diameter width can be known by multiplying p by the pixel measured at each marker point.
As shown in fig. 2, the automatic quantification system for retinal vessel diameter of fundus image provided by the embodiment of the invention comprises: the system comprises a video disc detection module 1, a preprocessing and blood vessel segmentation module 2 and a blood vessel width measurement and application module 3.
The video disc detection module 1 is used for realizing retina image disc detection.
The preprocessing and blood vessel segmentation module 2 is used for improving the image quality and extracting retinal blood vessels.
The blood vessel width measurement and application module 3 is used for realizing the measurement of the retinal blood vessel diameter and providing reference for clinicians.
The principle of application of the invention will be described in detail with reference to specific embodiments.
Retina optic disc detection:
as shown in fig. 4, the coarse positioning procedure provided by the embodiment of the present invention: finding the mass center of the optic disc, coarsely positioning the mass center into a circle, and intercepting a sub-image.
As shown in fig. 5, the fine positioning process provided by the embodiment of the present invention is: dynamic tracking process diagram, tracking result diagram, tracking curve thickening and fine positioning result.
Pretreatment and retinal vessel segmentation:
as shown in fig. 6, the pretreatment and retinal vascular segmentation process provided in the embodiment of the present invention is as follows: original image, contrast enhancement image, vessel initial segmentation, vessel after removal of the disturbance.
Retinal vessel width measurement:
as shown in fig. 7, the retinal vascular width measurement process provided in the embodiment of the present invention is: intercepting blood vessels, detecting edges, removing isolated noise, rotating blood vessels, establishing a scale and measuring results.
The foregoing description of the preferred embodiments of the invention is not intended to be limiting, but rather is intended to cover all modifications, equivalents, and alternatives falling within the spirit and principles of the invention.