CN109829942B - Automatic quantification method for retinal vessel diameter of fundus image - Google Patents
Automatic quantification method for retinal vessel diameter of fundus image Download PDFInfo
- Publication number
- CN109829942B CN109829942B CN201910128419.8A CN201910128419A CN109829942B CN 109829942 B CN109829942 B CN 109829942B CN 201910128419 A CN201910128419 A CN 201910128419A CN 109829942 B CN109829942 B CN 109829942B
- Authority
- CN
- China
- Prior art keywords
- image
- pixel
- retinal
- scale
- tracking
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 79
- 210000001210 retinal vessel Anatomy 0.000 title claims abstract description 50
- 238000011002 quantification Methods 0.000 title claims abstract description 14
- 210000004204 blood vessel Anatomy 0.000 claims abstract description 65
- 230000011218 segmentation Effects 0.000 claims abstract description 21
- 230000002792 vascular Effects 0.000 claims abstract description 19
- 230000002207 retinal effect Effects 0.000 claims abstract description 16
- 238000012545 processing Methods 0.000 claims abstract description 15
- 238000007781 pre-processing Methods 0.000 claims abstract description 14
- 238000004458 analytical method Methods 0.000 claims abstract description 9
- 238000003745 diagnosis Methods 0.000 claims abstract description 9
- 230000002708 enhancing effect Effects 0.000 claims abstract description 3
- 230000008569 process Effects 0.000 claims description 24
- 238000005259 measurement Methods 0.000 claims description 22
- 239000011159 matrix material Substances 0.000 claims description 12
- 238000001514 detection method Methods 0.000 claims description 11
- 210000001525 retina Anatomy 0.000 claims description 11
- 238000004364 calculation method Methods 0.000 claims description 9
- 238000003708 edge detection Methods 0.000 claims description 8
- 238000000605 extraction Methods 0.000 claims description 6
- 238000003384 imaging method Methods 0.000 claims description 6
- 210000003733 optic disk Anatomy 0.000 claims description 5
- 238000013139 quantization Methods 0.000 claims description 5
- 230000004256 retinal image Effects 0.000 claims description 5
- 230000003902 lesion Effects 0.000 claims description 4
- 230000009466 transformation Effects 0.000 claims description 4
- 238000012937 correction Methods 0.000 claims description 3
- 125000004122 cyclic group Chemical group 0.000 claims description 3
- 238000005286 illumination Methods 0.000 claims description 3
- 238000013507 mapping Methods 0.000 claims description 3
- 239000003550 marker Substances 0.000 claims description 3
- 238000003672 processing method Methods 0.000 claims description 3
- 238000004088 simulation Methods 0.000 claims description 3
- 230000003287 optical effect Effects 0.000 claims 5
- 206010025421 Macule Diseases 0.000 claims 1
- 241001632422 Radiola linoides Species 0.000 claims 1
- 238000010586 diagram Methods 0.000 description 9
- 230000000694 effects Effects 0.000 description 5
- 241001632427 Radiola Species 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 210000002189 macula lutea Anatomy 0.000 description 2
- 230000008719 thickening Effects 0.000 description 2
- 208000024172 Cardiovascular disease Diseases 0.000 description 1
- 208000026106 cerebrovascular disease Diseases 0.000 description 1
- 238000003759 clinical diagnosis Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 230000002526 effect on cardiovascular system Effects 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000008092 positive effect Effects 0.000 description 1
- 208000024891 symptom Diseases 0.000 description 1
Images
Abstract
The invention belongs to the technical field of medical image processing, and discloses an automatic quantification method for retinal vascular tube diameters of fundus images; automatically detecting the video disc by using a two-step positioning method combining thickness positioning and fine positioning; preprocessing the fundus image, including denoising and enhancing the fundus image; separating retinal blood vessels from fundus image background by adopting a multi-scale analysis method, and removing pseudo blood vessels in segmentation by utilizing mathematical morphology reconstruction; the vessel diameter is measured by using a boundary method, and the result is displayed. On the basis of fundus image characteristics, the invention firstly detects the optic disc, then enhances the image contrast, cuts out the retinal blood vessel, finally measures the size of the retinal blood vessel diameter, displays the result on a system interface, provides reference for early prediction and auxiliary diagnosis of clinicians, saves labor and improves efficiency.
Description
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)
wherein ,
(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)
wherein ,
(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.
Claims (2)
1. The automatic quantification method for the retinal vascular tube diameter of the fundus image is characterized by comprising the following steps of:
step one, detecting the video disc by using a two-step positioning method combining thickness positioning and fine positioning;
step two, preprocessing the fundus image, denoising and enhancing, separating retinal blood vessels from the background of the fundus image by adopting a multi-scale method, and removing pseudo blood vessels existing in segmentation;
measuring the vessel diameter by adopting a boundary method, and displaying the result on a system interface;
in the first step, the automatic detection of the disc positioning in the optic disc is the basis of the fundus image retinal vessel tracking and macula and lesion extraction work:
(1) Optic disc coarse positioning
Binarizing the image according to the gray level distribution condition of the fundus image, removing partial interference by adopting mathematical morphology knowledge to obtain a video disc interval, roughly positioning the mass center of the video disc, roughly positioning the video disc into a circular shape; assuming that the mass center of the optic disc is O, finding an upper point A, B, a lower point A, B, a left point C, D and a right point C, D on the edge of the optic disc along the horizontal direction and the vertical direction of the O point respectively, finding the mass center by using a least square fitting circle, and then obtaining an image of the optic disc;
the coordinates of centroid O (x o ，y o ) The calculation is as follows:
wherein xi and yi Respectively the abscissa and the ordinate of the video disc edge points, and N is the total number of the video 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 the coordinate of any 3 points in A, B, C, D, and the center coordinate mean value under 4 conditions is respectively used as the mass center coordinate of the video disc;
(2) Precision positioning of optic disk
The optical disc is coarsely positioned by utilizing the brightness information of the optical disc and the local blood vessel information, and the optimal optical disc boundary curve is conveniently obtained by adopting a dynamic contour tracking method based on a non-reinitialized variation level set, so that the optical disc positioning method has good fault tolerance and robustness, and realizes the accurate positioning of the optical disc:
the image energy function is expressed as:
in the omega table image field,
is a sign distance function; />
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;
in the second step, distortion phenomenon occurs in the individual images, and distortion correction is realized by adopting a linear interpolation method and a normalized image processing method; for the influence of non-ideal acquisition conditions and characteristic differences of an imaging sensor in the fundus image imaging process, the acquired images have the problems of noise, uneven illumination and lower image contrast; the preprocessing mainly comprises image denoising, contrast enhancement and mathematical morphology processing;
(1) Fundus image denoising and contrast enhancement based on Contourlet transformation; noise is suppressed in a Contourlet domain by adopting a soft threshold denoising method; carrying out enhancement treatment on the transformed sub-band coefficients in each 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)
wherein ,
(2) Retinal vessel segmentation based on multi-scale analysis linear tracking
The method comprises the steps of representing an image into a plurality of scales, and dividing a retina blood vessel into blood vessels with different areas and widths under different scales by using a multi-scale analysis method; 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 to obtain an initial blood vessel;
(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, T, of the pixel of the image at point (x, y) 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 The points belonging to the bright areas of the blood vessels of the retina image, the confidence of the pixel points 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 coefficient of a certain pixel in the confidence coefficient matrix is higher, the probability of the pixel belonging to the blood vessel is relatively higher; initially, setting elements under all scales in a confidence matrix to 0;
C w (x,y)：＝0 (6)
(2) initialization of linear tracking
The initialization process of linear tracking of retinal blood vessels is represented by formula (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 Tracking the pixel sequence for a new one;
(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, and does not contain 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 The updating is as follows:
C w (x,y):＝C w (x,y)+1,(x,y):＝(x c ,y c ) (10)
when all parameter values in the formula (9) are smaller than T, searching for the next seed pixel point (T: =t+1);
(4) multi-scale linear tracking
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 larger; through multiple simulation experiments, the dimension w= 3,5,7,9,11 is proper;
(5) initial extraction of retinal blood vessels
Carrying out initial estimation on the blood vessel by adopting a mapping quantization method; the initial vessel is determined by confidence matrix greater than T C Is constructed by pixel points of T C The value is equal to the value of the dimension W;
in the third step, the specific process of measuring the vessel diameter by a boundary method is as follows:
(1) Vessel measurement image region selection
Selecting a region of interest from the segmented retinal blood vessel images to intercept the blood vessel images for edge detection, and then measuring the width of the blood vessels of the sub-images to reduce the processing time;
(2) Calculation of vessel diameter width
The blood vessel direction in the sub-image is parallel to the horizontal direction as much as possible, if not, the sub-image is parallel to the horizontal direction as much as possible by adopting a rotation angle method; the automatic measurement of the retinal vascular caliber 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 expressed as:
wherein ,wi Is the width value of column i, (x) i ,y i) and (xj ,y j ) Coordinates of the upper and lower edges of the blood vessel respectively; the distance method has the advantages of wide application range, good stability and high precision;
(3) Scale positioning and pixel pitch calculation
In order to convert the vessel diameter marked by pixels into the actual size of the vessel diameter, the pixel spacing needs to be measured, so that the positioning of the scale in the image must be solved;
calculating the physical length of the pixels according to the minimum scale value and the number of pixels between the adjacent scale lines; setting a scale in an image with the same size as the intercepted sub-image, wherein the scale is 10mm multiplied by 10mm, performing edge detection processing on the scale image to obtain a binary image, setting the pixel value of the scale to be 1 and setting the pixel value of the background to be 0; 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 is defined by
Determining that the unit of p is mm/pixel;
(4) Converting pixel size of vessel diameter into physical size
The physical size of the retinal vessel diameter width is known by multiplying p by the measured pixel at each marker point.
2. A fundus image retinal vessel diameter automatic quantization system for implementing the fundus image retinal vessel diameter automatic quantization method of claim 1, characterized in that the fundus image retinal vessel diameter automatic quantization system comprises:
the optic disc detection module is used for realizing the optic disc detection of fundus images;
the image preprocessing and segmentation module is used for realizing fundus image denoising, contrast enhancement, retinal vessel segmentation and edge detection;
and the quantification and application module is used for realizing retinal vascular caliber measurement, clinical early prediction and auxiliary diagnosis.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910128419.8A CN109829942B (en) | 2019-02-21 | 2019-02-21 | Automatic quantification method for retinal vessel diameter of fundus image |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910128419.8A CN109829942B (en) | 2019-02-21 | 2019-02-21 | Automatic quantification method for retinal vessel diameter of fundus image |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109829942A CN109829942A (en) | 2019-05-31 |
| CN109829942B true CN109829942B (en) | 2023-04-28 |
Family
ID=66863968
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910128419.8A Active CN109829942B (en) | 2019-02-21 | 2019-02-21 | Automatic quantification method for retinal vessel diameter of fundus image |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109829942B (en) |
Families Citing this family (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110448267B (en) * | 2019-09-06 | 2021-05-25 | 重庆贝奥新视野医疗设备有限公司 | Multimode fundus dynamic imaging analysis system and method |
| CN112668360A (en) * | 2019-10-15 | 2021-04-16 | 深圳市理邦精密仪器股份有限公司 | Pupil measuring method and device, medical equipment and storage medium |
| CN111489353B (en) * | 2020-05-07 | 2023-06-23 | 清华大学深圳国际研究生院 | Fundus image fovea positioning method |
| CN112617789A (en) * | 2020-07-28 | 2021-04-09 | 上海大学 | Laser speckle blood flow imaging method and system |
| CN112288794B (en) * | 2020-09-04 | 2021-09-07 | 深圳硅基智能科技有限公司 | Method and device for measuring blood vessel diameter of fundus image |
| CN111840794A (en) * | 2020-09-07 | 2020-10-30 | 复旦大学附属眼耳鼻喉科医院 | Method and system for rapidly measuring implantation depth of cochlear implant electrode |
| CN113344895A (en) * | 2021-06-23 | 2021-09-03 | 依未科技(北京)有限公司 | High-precision fundus blood vessel diameter measuring method, device, medium and equipment |
| CN113724315B (en) * | 2021-09-03 | 2024-04-02 | 上海海事大学 | Fundus retina blood vessel width measuring method, electronic equipment and computer readable storage medium |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108986106A (en) * | 2017-12-15 | 2018-12-11 | 浙江中医药大学 | Retinal vessel automatic division method towards glaucoma clinical diagnosis |
Family Cites Families (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP1691714A4 (en) * | 2003-11-25 | 2010-01-27 | F D Cardio Ltd | Stent positioning using inflation tube |
| US20050157939A1 (en) * | 2004-01-16 | 2005-07-21 | Mark Arsenault | Processes, products and systems for enhancing images of blood vessels |
| CA2728662C (en) * | 2008-10-14 | 2020-06-16 | Lightlab Imaging, Inc. | Methods for stent strut detection and related measurement and display using optical coherence tomography |
| CN101984916B (en) * | 2010-11-17 | 2012-10-31 | 哈尔滨工程大学 | Blood vessel diameter measuring method based on digital image processing technology |
| JP6358590B2 (en) * | 2013-08-09 | 2018-07-18 | 富士通株式会社 | Blood vessel data generation device, blood vessel data generation method, and blood vessel data generation program |
| US9757023B2 (en) * | 2015-05-27 | 2017-09-12 | The Regents Of The University Of Michigan | Optic disc detection in retinal autofluorescence images |
| CN107292835B (en) * | 2017-05-31 | 2020-03-13 | 瑞达昇医疗科技(大连)有限公司 | Method and device for automatically vectorizing retinal blood vessels of fundus image |
| CN109166124B (en) * | 2018-11-20 | 2021-12-14 | 中南大学 | Retinal blood vessel morphology quantification method based on connected region |
-
2019
- 2019-02-21 CN CN201910128419.8A patent/CN109829942B/en active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108986106A (en) * | 2017-12-15 | 2018-12-11 | 浙江中医药大学 | Retinal vessel automatic division method towards glaucoma clinical diagnosis |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109829942A (en) | 2019-05-31 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109829942B (en) | Automatic quantification method for retinal vessel diameter of fundus image | |
| CN109493954B (en) | SD-OCT image retinopathy detection system based on category distinguishing and positioning | |
| US8605979B2 (en) | Automatic detection and quantification of plaque in the coronary arteries of subjects from CT scans | |
| CN105279759B (en) | The abdominal cavity aortic aneurysm outline dividing method constrained with reference to context information arrowband | |
| CN110717888B (en) | Automatic identification method for intravascular Optical Coherence Tomography (OCT) vessel wall inner contour | |
| CN108186051B (en) | Image processing method and system for automatically measuring double-apical-diameter length of fetus from ultrasonic image | |
| CN108378869B (en) | Image processing method and processing system for automatically measuring head circumference length of fetus from ultrasonic image | |
| CN109816655B (en) | Pulmonary nodule image feature detection method based on CT image | |
| WO2012126070A1 (en) | Automatic volumetric analysis and 3d registration of cross sectional oct images of a stent in a body vessel | |
| CN109035227A (en) | The system that lung tumors detection and diagnosis is carried out to CT image | |
| CN108846827B (en) | Method for rapidly segmenting fundus optic disk based on multiple circles | |
| CN108961278B (en) | Method and system for abdominal wall muscle segmentation based on image data | |
| CN110706218B (en) | Breast tumor positioning analysis method based on dynamic enhanced magnetic resonance imaging | |
| CN110060261B (en) | Blood vessel segmentation method based on optical coherence tomography system | |
| CN113012127A (en) | Cardiothoracic ratio measuring method based on chest medical image | |
| CN106372593B (en) | Optic disk area positioning method based on vascular convergence | |
| CN111145155B (en) | Meibomian gland identification method | |
| Purnama et al. | Follicle detection on the usg images to support determination of polycystic ovary syndrome | |
| CN109492653B (en) | Method and device for measuring breast lesion volume, computer equipment and storage medium | |
| CN110930346B (en) | Automatic detection method and storage device for eyeground image microangioma | |
| CN112529918B (en) | Method, device and equipment for segmenting brain room area in brain CT image | |
| CN115841472A (en) | Method, device, equipment and storage medium for identifying high-density characteristics of middle cerebral artery | |
| Rashid et al. | Segmenting melanoma lesion using single shot detector (SSD) and level set segmentation technique | |
| CN109949243B (en) | Calcification artifact eliminating method and device and computer storage medium | |
| Jagoe et al. | Quantification of retinal damage during cardiopulmonary bypass |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |