CN102842136A - Optic disc projection location method synthesizing vascular distribution with video disc appearance characteristics - Google Patents

Abstract

The invention discloses an optic disc projection location method synthesizing vascular distribution with video disc appearance characteristics, which comprises the following steps of: (1) extracting an interest retina fundus image region by use of a mask operation; (2) carrying out normalization enhancement on a fundus image based on an image observation model; (3) achieving extraction and segmentation of the fundus image by using a non-vascular structure inhibition operator and combining with a hysteresis multi-threshold processing technology; (4) setting a vertical window of double main blood vessel width to slide on a blood vessel segmentation image along the horizontal direction and calculating a vascular distribution degree D(x) at each horizontal position x to obtain a distribution degree curve of a horizontal projection, wherein the minimum value point of the curve is determined to be optic disc horizontal coordinate xod; (5) setting a rectangular window side length of which is equal to the optic disc diameter to slide up and down at the horizontal coordinate xod, estimating local region brightness IN(xod, y) and edge gradient information gN(xod, y) respectively, drawing a vertical projection curve reflecting a change of a characteristic value f(y)=IN(xod, y)*gN(xod, y), wherein a maximum value point of the curve is optic disc vertical coordinate yod. The method has a simple algorithm, a high success rate and excellent robustness.

Description

The optic disk projection localization method of a kind of comprehensive vascular distribution and optic disk appearance characteristics

Technical field

The present invention relates to the method that the optic disk in a kind of retina eye fundus image is located automatically, the optic disk projection localization method of particularly a kind of comprehensive vascular distribution and optic disk appearance characteristics.

Background technology

Optic disk is one of amphiblestroid major physiological structure, and optic nerve and blood vessel get into eye and extend to ambient radiation from this zone, in the retina eye fundus image, show as circular glassy yellow zone, comprise a large amount of thicker blood vessels simultaneously.The researchist pays special attention to the automatic location of optic disk always, because the accurate location of optic disk helps the diagnosis of eyeground pathological changes, for example: a lot of optic disk partitioning algorithms need provide initial seed point in advance; Because the distance between optic disk and the macula lutea (the sharpest zone of vision) is approximately constant, after the optic disk position is confirmed, confirms that macular region is also just easy; Optic disk and large-sized glassy yellow pathology (overflow pathology) are mixed easily and are disturbed, and the location of optic disk can be got rid of it from the lesion candidates object.

The optic disk localization method roughly is divided into following several types at present:

1, based on the method for optic disk appearance characteristics.These class methods adopt outward appearance such as brightness, contrast, the shape information of optic disk as location feature usually, and the center of generally that brightness maximum or contrast is the strongest border circular areas is as the reference position of optic disk.Utilize appearance characteristics merely, the success ratio of the optic disk location in the normal preferably eye fundus image of quality is higher, but in the pathology image, because pathology causes the optic disk outward appearance to change, or other large tracts of land light tone lesion region is disturbed optic disk location easy error.

2, based on the method for vessel properties.This mainly is to be initiation region and the vessel branch point that blood vessel gets into eye according to optic disk, and should the zone blood vessel the thickest, the vessel density maximum.The optic disk location algorithm that detects based on vessel properties is more complicated all generally, and consuming time, and in inferior quality or pathology image, the blood vessel feature detection is still a relatively problem of difficulty simultaneously.

3, based on the regional manifold method of optic disk.Such method synthesis utilizes the appearance characteristics and the vessel properties of optic disk, and bearing accuracy is higher relatively.But these class methods need be set up than complex model, and its complex algorithm is not easy to real-time application.

In recent years, lot of domestic and foreign scholar has carried out a large amount of research to the optic disk location, and however, the overwhelming majority is studied and is being difficult to take into account aspect the accuracy rate of optic disk location and the speed two, causes the optic disk localization method to be difficult to adapt to the real-time processing requirements of real system.

Summary of the invention

In order to solve the above-mentioned technical matters that existing optic disk location exists, the invention provides the optic disk projection localization method of a kind of comprehensive vascular distribution and appearance characteristics.The present invention has fully utilized the vascular distribution and the appearance characteristics information in optic disk zone; Simultaneously owing to the two-dimensional localization problem with complicacy converts simple relatively one dimension orientation problem into; Algorithm is realized simple, and success ratio is high, has good robustness; Comprehensive accuracy rate and processing speed, the present invention meets the real-time processing requirements of real system.

The technical scheme that the present invention solves the problems of the technologies described above may further comprise the steps:

(1) get the brightest pixel intensity value of red channel component in the original color eye fundus image 5% as threshold value; And go out the bianry image of foreground area according to this Threshold Segmentation; Bianry image is carried out the morphological erosion operation; Obtain the mask template, obtain the area-of-interest of eye fundus image through mask process;

(2), the area-of-interest after the mask process is carried out the normalization of illumination and contrast and handle the retina eye fundus image that is enhanced according to the image viewing model;

(3) the enhancing retina eye fundus image that step (2) is obtained adopts the non-blood vessel structure of Gabor filtering to suppress operator, and many threshold process technology that combine simultaneously to lag behind realize the eye fundus image vessel extraction and cut apart, and obtain blood vessel segmentation figure;

(4) vertical window is set, the height of window is a picture altitude, and width is 2 times a main blood vessel width, from left to right slides along the horizontal direction of blood vessel segmentation figure; In window, calculate the blood vessel distribution degree value D (x) at each horizontal level x place, draw the horizontal projection curve, find the minimum point of this dimension curve to confirm as the horizontal coordinate x of optic disk according to blood vessel distribution degree value _Od

(5) rectangular window is set, the length of window and width are the size of disc diameter, and the horizontal coordinate at its center is fixed optic disk horizontal coordinate x _Od, respectively in former gray-scale map and initial phase

Gabor filter response figure on moving window from top to bottom vertically; And interior mean flow rate of estimation window and average filter response value; Average brightness value and average filter response value are multiplied each other as longitudinal projection's value at y place, corresponding upright position; According to longitudinal projection's value paint projection curve, the coordinate of maximizing point is the vertical coordinate of optic disk.

In the optic disk projection localization method of above-mentioned comprehensive vascular distribution and optic disk appearance characteristics, the morphological erosion operation in the said step (1) is to adopt 9 * 9 square structure elements that bianry image is carried out erosion operation.

In the optic disk projection localization method of above-mentioned comprehensive vascular distribution and optic disk appearance characteristics, the concrete steps of said step (2) are:

1. image is divided into size and is the piecemeal S of s _i, s is M * N of 1/6～1/8, M * N is an original image size, to each piecemeal S _iCalculate the gray average in this piecemeal

And standard deviation Adopt the bicubic interpolation method to try to achieve in the entire image every average again

With standard deviation

2. pass through Manhabolios apart from d _M(x y), judges (whether x y) belongs to the background area to certain pixel, if it is apart from the neighborhood average

Less than certain threshold value t, threshold value t=1, promptly

$d_M(x,y)=\left|\frac{I(x,y)-{\widehat{\mathrm{\&mu;}}}_N(x,y)}{{\widehat{\mathrm{\&sigma;}}}_N(x,y)}\right|<t---\left(11\right)$

Then this pixel belongs to the background area, otherwise belongs to foreground area;

3. in the background area that splits illumination and contrast are estimated sampled point (s in each piecemeal _x, s _y) illumination

And contrast

Average and standard deviation by background pixel point in the piecemeal are estimated to obtain, and the illumination of all the other non-sampled points and non-background pixel point and the contrast drift factor then obtain through the bicubic interpolation method;

4. according to formula (7), realize that the normalization of image strengthens,

$\widehat{I^0}(x,y)=\frac{I(x,y)-\hat{L}(x,y)}{\hat{C}(x,y)}---\left(7\right)$

I (x; Y) be that original observed image is at (x; Y) gray values of pixel points;

and is respectively the illumination and the contrast drift factor of this some estimation, and

is the gray-scale value after the normalization.

In the optic disk projection localization method of above-mentioned comprehensive vascular distribution and optic disk appearance characteristics, said step (3) concrete steps are:

1. select for use 12 Gabor filtering templates on the different directions respectively image to be carried out filtering; Each adjacent direction phase difference of pi/12, filtered image is expression with

:

Wherein,

is the Gabor function; The f representing input images; Keep the peak response of each pixel at 12 different directions

$H\left(z\right)=\left\{\begin{array}{cc}0,& z<0\\ z,& z\&GreaterEqual;0\end{array}\right.---\left(15\right)$

2. be directed against the Gabor function when out of phase

; Detection characteristic to the different characteristic pattern; Employing formula (16) suppresses to handle to bright fringe and non-vessel boundary

$v_{\mathrm{\&lambda;},\mathrm{\&delta;}}(x,y)=H(r_{\mathrm{\&lambda;},\mathrm{\&delta;},\mathrm{\&pi;}}^{\mathrm{AH}}(x,y)-{\mathrm{\&alpha;}}_1{r}_{\mathrm{\&lambda;},\mathrm{\&delta;},0}^{\mathrm{AH}}(x,y)-{\mathrm{\&alpha;}}_2{r}_{\mathrm{\&lambda;},\mathrm{\&delta;},\mathrm{\&pi;}/2}^{\mathrm{AH}}(x,y))---\left(16\right)$

Wherein

Expression respectively

0, the filter response during pi/2, α ₁, α ₂Represent inhibition component respectively to bright fringe and non-vessel boundary

The inhibition degree;

3. choose many group threshold values

and

and consider the connectedness of blood vessel; The two adjacent groups threshold interval should be overlapping; Promptly

selected threshold value suppresses to handle afterwards image to

and

to filtering and carries out list threshold binarization is cut apart, and obtains corresponding binary image

$b_{\mathrm{\&lambda;},\mathrm{\&delta;}}^i(x,y)=H{T}^i\left(v_{\mathrm{\&lambda;},\mathrm{\&delta;}}(x,y)\right),i\&Element;\{1,...,N_i\}---\left(17\right)$

HT representes the processing of standard hysteresis threshold;

4. each bianry image is expressed as the union of some connected components

$b_{\mathrm{\&lambda;},\mathrm{\&delta;}}^i=\underset{k}{\&cup;}{C}_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{i,k}---\left(18\right)$

Wherein, k connected component in the binary image of

expression threshold value when being chosen for and

.Utilize the morphology expansive working, select 2 * 2 square structure element q for use ₂, right

Each branch

Carry out expansion process, obtain:

$D_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{i,k}=C_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{i,k}\&CirclePlus;q_2---\left(19\right)$

5. final blood vessel structure V _{λ, δ}Can handle acquisition in regular turn through following formula,

$B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{(N_t,N_t-1)}=\underset{(b_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t}\&CirclePlus;q_2)\&cap;D_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t-1,k}\&NotEqual;\mathrm{\&phi;}}{\&cup;}(b_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t}\&cup;C_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t-1,k})---\left(20\right)$

$B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{(N_t-1,N_t-2)}=\underset{(B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{(N_t,N_t-1)}\&CirclePlus;q_2)\&cap;D_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t-2,k}\&NotEqual;\mathrm{\&phi;}}{\&cup;}(B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t,N_t-1}\&cup;C_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t-2,k})---\left(21\right)$

.

.

.

$V_{\mathrm{\&lambda;},\mathrm{\&delta;}}=B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{\left(\mathrm{2,1}\right)}=\underset{(B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{\left(\mathrm{3,2}\right)}\&CirclePlus;q_2)\&cap;D_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{1,k}\&NotEqual;\mathrm{\&phi;}}{\&cup;}(B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{\left(\mathrm{3,2}\right)}\&cup;C_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{1,k})---\left(22\right)$

6. because blood vessel presents the network structure of connection in image, consider to remove the branch of all length less than 10 pixels.

The optic disk projection localization method of above-mentioned comprehensive vascular distribution and optic disk appearance characteristics, the computing formula of described step (4) medium vessels distribution degree value D (x) is:

$D\left(x\right)=\frac{(-1)\underset{i=1}{\overset{n_x}{\mathrm{\&Sigma;}}}{p}_i{\log }_2{p}_i}{\underset{i}{\max }\left\{p_i\right\}}---\left(23\right)$

In the above-mentioned formula

Be i section blood vessel connected component proportion in this position vertical window, m _iRepresent that the i section is communicated with the blood vessel pixel count, M representes the total number of pixels of blood vessel in the horizontal slip window,

The number of pixels of representing largest connected vessel segment, n _xFor blood vessel total in this position level moving window is communicated with number.

In the optic disk projection localization method of above-mentioned comprehensive vascular distribution and optic disk appearance characteristics, the concrete steps of said step (5) are:

1. define window rectangular window W _v, W _vLength and width be the disc diameter size, and the fixed optic disk horizontal coordinate of the horizontal coordinate at its center x _OdThe place, ordinate y changes with the upper and lower slip of window;

2. moving window W from top to bottom vertically _vAnd on the response diagram of θ=90 ° phase place, estimate mean flow rate and average filter response value in the window at former figure and Gabor filtering respectively; According to formula (24), average brightness value and average filter response value are multiplied each other as longitudinal projection's value at y place, corresponding upright position

f(y)＝I _N(x _od,y)*g _N(x _od,y) （24）

I _N(x _Od, be y) with (x _Od, y) be the average brightness value of the rectangular window neighborhood at center, g _N(x _Od, y) be average Gabor filter response estimated value in this neighborhood;

3. draw longitudinal projection's curve according to f (y), the coordinate of maximizing point is the vertical coordinate y of optic disk _Od

Compared with prior art, technique effect of the present invention is:

1, after existing vessel extraction method is seldom considered PVR,, is difficult to the correct extraction blood vessel structure owing to receive the lesion region boundary interference.Optic disk, border hemorrhage or other lesion region have the characteristic of similar blood vessel, and these borders exist in pathology retina eye fundus image in a large number.How to suppress non-angiosomes edge effect in the blood vessel feature extraction phases is that the present invention focuses on one of problem that solves.The present invention utilizes under the out of phase; The Gabor filter function is to the response characteristic at different characteristic pattern (dark fringe, bright fringes and edge); Propose a kind of non-blood vessel structure and suppress operator based on Gabor filtering; Many threshold process that combine simultaneously to lag behind are technological, realize that finally the eye fundus image blood vessel structure extracts.The used vessel extraction method of the present invention can suppress the non-blood vessel pixel in the pathology image to a great extent, and more complete reservation blood vessel structure, greatly improved the accuracy of vessel extraction.

2, the present invention has fully utilized vascular distribution and optic disk appearance characteristics.At first (optic disk zone blood vessel vertically stretches according to the distribution character of blood vessel; Connectivity is good, distributes to concentrate) confirm the horizontal coordinate of optic disk, and then utilize optic disk appearance characteristics (brightness, edge gradient) location vertical coordinate; Owing to convert the two-dimensional localization problem of complicacy into simple relatively two one dimension orientation problems; Algorithm is realized simple, and success ratio is high, has good robustness.

Description of drawings

Fig. 1 is a process flow diagram of the present invention.

Fig. 2 is a colored eye fundus image mask template synoptic diagram among the present invention.

Fig. 3 is the eye fundus image normalization synoptic diagram among the present invention.

Fig. 4 is that out of phase is a Gabor function space template synoptic diagram.

Fig. 5 is certain eyeground coloured image and and corresponding vascular distribution figure thereof.

Fig. 6 is the rectangular window synoptic diagram on eyeground gray level image and

phase place Gabor filtering image thereof.

Fig. 7 is the coordinate projection location synoptic diagram of certain eye fundus image.

Embodiment

Below will combine accompanying drawing and embodiment that the present invention is explained further details.

As shown in Figure 1, optic disk projection localization method of the present invention, its idiographic flow is:

1, eye fundus image mask process.Shown in Fig. 2 (a); The retina eye fundus image comprises dark-coloured background and retina fundus imaging zone usually; Wherein retina fundus imaging zone is ophthalmologist, oculist's a area-of-interest, in order to prevent the interference on background and border, is necessary to obtain the mask template of area-of-interest.Because the red component of colored eyeground figure is near saturated, therefore lighting condition in the time of reflecting imaging gets red channel component I in the original color eye fundus image _RThe brightest pixel intensity value 5% as segmentation threshold t _b, and according to the bianry image of this Threshold Segmentation foreground area,

t _b＝0.05*max(I _R) （1）

$I_{\mathrm{bin}}(x,y)=\left\{\begin{array}{cc}1,& I_R(x,y)>t_b\\ 0,& I_R(x,y)\&le;t_b\end{array}\right.---\left(2\right)$

For getting rid of the interference on border as far as possible, be necessary binary map I _BinCorrode operation, the present invention adopts 9 * 9 square structure q ₉Element is to I _BinCarry out erosion operation, obtain final mask template I _Mask, referring to Fig. 2 (b).

representes erosion operation.Adopt the mask operation to be not difficult to obtain the area-of-interest of image.

2, eye fundus image normalization strengthens.Since the influence of illumination in the image acquisition procedures, camber reflection, visual angle, eye motion, non-homogeneous illumination of eye fundus image ubiquity and contrast phenomenon, so the normalization of eye fundus image is the necessary link of subsequent treatment.

Original retina eye fundus image can be regarded the stack combinations of desirable background image and foreground image as, and in image observation (obtaining) process, image illumination and contrast can produce deformation on the original image basis, and this distortion can use following observation model to describe:

I(x,y)＝C(x,y)*I ⁰(x,y)+L(x,y) （4）

$I^0(x,y)=I_b^0(x,y)+I_f^0(x,y)---\left(5\right)$

I ⁰Being original eye fundus image, is the desirable eye fundus image that does not comprise non-homogeneous illumination and contrast influence.

is original background image, and

is original foreground image.L (x, y) and C (x is respectively that (I is the image that finally observes to pixel in the acquisition process for x, y) the illumination drift factor and the contrast factor of drifting about y).According to formula (4), can obtain original eye fundus image,

$I^0(x,y)=\frac{I(x,y)-L(x,y)}{C(x,y)}---\left(6\right)$

The normally the unknown of the real illumination drift factor and contrast factor can only be estimated its value from observed image.Therefore formula (6) becomes

$\widehat{I^0}(x,y)=\frac{I(x,y)-\hat{L}(x,y)}{\hat{C}(x,y)}---\left(7\right)$

Be to original eye fundus image I ⁰Estimation.Must estimate illumination drift L and contrast drift C for obtaining

.Because foreground portion divides characteristic widely different; Therefore the characteristic of prospect part

is beyond expression of words; On the other hand; Background

changes mild; Can describe by normal distribution law

$I_b^0(x,y)~N({\mathrm{\&mu;}}_b,{\mathrm{\&sigma;}}_b)---\left(8\right)$

μ wherein _bThe desirable uniform illumination value of expression background, standard deviation sigma _bReflected that then retina changes in the natural characteristic of spatial domain.

Because the irregular acute variation that can't cause image of illumination means that illumination drift L and contrast drift C mainly concentrate on the low-frequency component of frequency spectrum, so background image is a main foundation of estimating L and C.For realizing the estimation of the drift factor; Might as well original image be divided into foreground image and background image; Can notice be placed on background parts like this; When only considering the background area, can make relevant foreground area

so formula (4) be reduced to:

$I(x,y)=C(x,y)*I^0(x,y)+L(x,y)$

$=C(x,y)*[I_b^0(x,y)+I_f^0(x,y)]+L(x,y)---\left(9\right)$

$=C(x,y)*I_b^0(x,y)+C(x,y)*I_f^0(x,y)+L(x,y)$

is the set that background area pixels constitutes.Obviously, background accord with normal distribution characteristic:

$I_b^0(x,y)~N({\mathrm{\&mu;}}_b,{\mathrm{\&sigma;}}_b)~N(L(x,y),C(x,y))---\left(10\right)$

Can know that from the analysis of last surface model the average and the standard deviation of background area

that can be through observed image are come

and

in the estimated background and finally recovered original image according to formula (7).Whole process can be divided into for two steps: the first step at first is a background extracting; Second step estimated through the background area that the illumination and the contrast drift factor

and

were because the green channel contrast in the RGB image is the strongest then; Other two passage can better reflect the characteristic on eyeground relatively, all on green channel, carries out in this follow-up processing.Its concrete implementation procedure is following:

1. image is divided into size and is the piecemeal S of s (be about M * N of 1/6～1/8, M * N is an original image size) _i, then to each piecemeal S _iCalculate the average in this piecemeal

And standard deviation

Adopt the bicubic interpolation method to try to achieve in the entire image at last at every With

Value;

2. pass through Manhabolios apart from d _M(x y), judges (whether x y) belongs to the background area to certain pixel, if it is apart from the neighborhood average

Less than certain threshold value t (threshold value t=1), promptly

$d_M(x,y)=\left|\frac{I(x,y)-{\widehat{\mathrm{\&mu;}}}_N(x,y)}{{\widehat{\mathrm{\&sigma;}}}_N(x,y)}\right|<t---\left(11\right)$

Then this pixel belongs to the background area, otherwise belongs to foreground area;

3. in the background area that splits illumination and contrast are estimated.Because background is the part of view picture eye fundus image, calculate for avoiding pointwise, adopt the similar 1. sample point sampling estimation in step.Sampled point (s in each piecemeal _x, s _y) illumination

And contrast Average and standard deviation by background pixel point in the piecemeal are estimated to obtain, and the illumination of all the other non-sampled points and non-background pixel point and the contrast drift factor then obtain through the bicubic interpolation method.Finally, obtain

and

of each pixel of view picture eye fundus image

4. according to formula (7), realize that the normalization of image strengthens.

Fig. 3 (a-e) is the normalization processing procedure synoptic diagram of certain exemplary plot; Wherein Fig. 3 (a) is an original image; Fig. 3 (b) is the background estimating figure of image, and Fig. 3 (c) and Fig. 3 (d) are respectively the estimation synoptic diagram of the brightness drift factor and the contrast drift factor, and Fig. 3 (e) is final normalization coloured picture.

3, vessel extraction and cutting apart.Though existing a lot of documents have proposed effective retinal vessel method for distilling; But seldom document is mentioned the vessel extraction problem in the pathology image; the edge of lesion region and optic disk is when having strong contrast; its matched filtering response (for example Gabor wave filter) is stronger than weak, the thin blood vessel of great majority, and therefore when adopting simple threshold values to cut apart, relevant edge usually is mistaken for blood vessel.For obtaining as far as possible accurately complete blood vessel structure, the present invention proposes a kind of non-blood vessel and suppresses operator, combines many hysteresis threshold technology to realize the extraction of blood vessel simultaneously.Idiographic flow is following:

1. because blood vessel is long and narrow dark-coloured strip or linear structure, and its xsect is similar to Gaussian distribution, we adopt has the Gabor wave filter realization blood vessel Feature Extraction that can adjust visual field and direction,

$\tilde{x}=x\cos \mathrm{\&theta;}+y\sin \mathrm{\&theta;}$

$\tilde{y}=-x\sin \mathrm{\&theta;}+y\cos \mathrm{\&theta;}$

Wherein, γ is the space length breadth ratio, has determined to accept the ovality of visual field, is made as constant 0.5 among the present invention, and standard deviation δ confirms to accept the size of visual field.Parameter lambda is the wavelength of cosine factor.Ratio delta/λ determines the bandwidth of spatial frequency domain, has influence on the fringe area quantity of excitatoty and inhibition parallel in the visual field, is taken as definite value δ/λ=0.56 in the literary composition.Angle parameter theta, and θ ∈ [0, π), be predetermined filtering direction.Phase shift parameter

OK

function symmetry.

Gabor filter response

by a Gabor function

convolution of the input image f calculated.Select ({ θ on 12 different directions for use _i=0, π/12,2 π/12 ... 11 π/12}|, each adjacent direction phase difference of pi/12) Gabor filtering template respectively image is carried out filtering, filtered image is used

Expression:

Keep the peak response of each pixel at 12 different directions,

$H\left(z\right)=\left\{\begin{array}{cc}0,& z<0\\ z,& z\&GreaterEqual;0\end{array}\right.---\left(15\right)$

Wherein, H is the half-wave function of adjusting.Under ideal conditions; The blood vessel center pixel has peak response when phase deviation parameter

, response approaches 0 when

.And there is the peak response value at the edge when or pi/2.

2. blood vessel presents the dark fringe structure in the retina gray level image.When Gabor filter function

; The kernel function space structure is shown in Fig. 4 (a), and this moment, the Gabor filter function had peak response to dark fringe; During as

; Kernel function space structure such as Fig. 4 (b), this moment, the Gabor filter function was maximum to the bright fringes response; During as

; Kernel function space structure such as Fig. 4 (c), this moment, there was peak response at the edge.

To the Gabor function when out of phase

; Detection characteristic to the different characteristic pattern; Employing formula (16) suppresses to handle to bright fringe and non-vessel boundary

$v_{\mathrm{\&lambda;},\mathrm{\&delta;}}(x,y)=H(r_{\mathrm{\&lambda;},\mathrm{\&delta;},\mathrm{\&pi;}}^{\mathrm{AH}}(x,y)-{\mathrm{\&alpha;}}_1{r}_{\mathrm{\&lambda;},\mathrm{\&delta;},0}^{\mathrm{AH}}(x,y)-{\mathrm{\&alpha;}}_2{r}_{\mathrm{\&lambda;},\mathrm{\&delta;},\mathrm{\&pi;}/2}^{\mathrm{AH}}(x,y))---\left(16\right)$

Wherein

Expression respectively (dark fringe), 0 (bright fringes), the filter response during pi/2 (edge), α ₁, α ₂Represent inhibition component respectively to bright fringe and non-vessel boundary The inhibition degree.In the actual extracting process, select α ₁=0.5, α ₂=1, suppress algorithm this moment all has zones such as bright fringes and non-vessel boundaries and comprehensively suppresses effect preferably.

3. cut apart through simple threshold values and can extract main blood vessel structure.Although most of non-blood vessel pixel is effectively suppressed after the inhibition operator is handled; But still there are the point-like and the residual pixel of short-term shape of some distributions; For removing these residual risks, on standard hysteresis threshold method basis, the present invention proposes a kind of many hysteresis threshold blood vessel segmentation method.

The standard hysteresis threshold comprises a pair of threshold value t _lAnd t _h, t _l＜t _hThe Gabor response is greater than t _hPixel when output two-value vessel graph, be appointed as 1, and the Gabor response is less than t _lPixel be designated as 0, response is at t _lAnd t _hBetween pixel, if can find a linkage path can with a certain response greater than t _hPoint be communicated with, the pixel response value on the passage is greater than t _l, then this pixel also is designated as 1.

Based on standard hysteresis threshold technology, choose many group threshold values

With

Consider the connectedness of blood vessel, the two adjacent groups threshold interval should be overlapping, promptly For example, in realizing example, adopt 4 pairs of threshold values

With

I=1 ..., 4, threshold parameter is to being chosen for: [0.02,0.04] * max (v _{λ, δ}), [0.03,0.05] * max (v _{λ, δ}), [0.04,0.06] * max (v _{λ, δ}), [0.05,0.08] * max (v _{λ, δ}), max (v _{λ, δ}) be v in the entire image _{λ, δ}Maximal value.Selected threshold value suppresses to handle afterwards image to

and

to filtering carries out list threshold binarization is cut apart, and obtains corresponding binary image

$b_{\mathrm{\&lambda;},\mathrm{\&delta;}}^i(x,y)=H{T}^i\left(v_{\mathrm{\&lambda;},\mathrm{\&delta;}}(x,y)\right),i\&Element;\{1,...,N_i\}---\left(17\right)$

HT representes the processing of standard hysteresis threshold.

4. each binary map is expressed as the union of some connected components

$b_{\mathrm{\&lambda;},\mathrm{\&delta;}}^i=\underset{k}{\&cup;}{C}_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{i,k}---\left(18\right)$

Wherein, k connected component in the binary image of

expression threshold value when being chosen for

and

.Utilize the morphology expansive working, select 2 * 2 square structure element q for use ₂, right

Each branch

Carry out expansion process, obtain:

$D_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{i,k}=C_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{i,k}\&CirclePlus;q_2---\left(19\right)$

5. final blood vessel structure V _{λ, δ}Can handle acquisition in regular turn through following formula,

$B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{(N_t,N_t-1)}=\underset{(b_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t}\&CirclePlus;q_2)\&cap;D_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t-1,k}\&NotEqual;\mathrm{\&phi;}}{\&cup;}(b_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t}\&cup;C_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t-1,k})---\left(20\right)$

$B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{(N_t-1,N_t-2)}=\underset{(B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{(N_t,N_t-1)}\&CirclePlus;q_2)\&cap;D_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t-2,k}\&NotEqual;\mathrm{\&phi;}}{\&cup;}(B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t,N_t-1}\&cup;C_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t-2,k})---\left(21\right)$

.

.

.

$V_{\mathrm{\&lambda;},\mathrm{\&delta;}}=B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{\left(\mathrm{2,1}\right)}=\underset{(B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{\left(\mathrm{3,2}\right)}\&CirclePlus;q_2)\&cap;D_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{1,k}\&NotEqual;\mathrm{\&phi;}}{\&cup;}(B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{\left(\mathrm{3,2}\right)}\&cup;C_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{1,k})---\left(22\right)$

6. because blood vessel presents the network structure of connection in image, consider to remove the branch of all length less than 10 pixels.

The blood vessel trunk structure that the vessel extraction method that the present invention adopts is extracted is complete, clean, can get rid of most of non-blood vessel pixel.

4, locate based on the horizontal coordinate of vascular distribution characteristic.Fig. 5 (b) is the handmarking's of certain eye fundus image (Fig. 5 (a)) a blood vessel binary map, in the width that limits is about 2 times the vertical window of main blood vessel width, is not difficult to observe in the varying level position, and there were significant differences for vascular distribution.In the optic disk position (Fig. 5 (b)), blood vessel UNICOM number is few, and it is very concentrated to distribute, and local vascular density is big, and in other position, blood vessel connection number is many and distribute scattered relatively.According to this distribution character, the idiographic flow of optic disk horizontal coordinate location is following:

1. define window W _h(highly be picture altitude, width is approximately 2 times main blood vessel width), and it is centrally located in horizontal coordinate point x place;

Because the connectedness of blood vessel; The blood vessel trunk structure is even more important for the calculating of blood vessel distribution degree; Disturb for getting rid of weak, thin blood vessel as far as possible, can earlier some short little vessel segments in the vertical window of definition be removed before the calculating blood vessel distribution degree, rule of thumb; Threshold value th=30 promptly neglects and is communicated with the blood vessel pixel less than 30 segment blood vessel.

2. at the horizontal direction of blood vessel segmentation figure moving window W from left to right _h, and calculate the blood vessel distribution degree value D (x) in each window, as the projection value at horizontal level x place;

$D\left(x\right)=\frac{(-1)\underset{i=1}{\overset{n_x}{\mathrm{\&Sigma;}}}{p}_i{\log }_2{p}_i}{\underset{i}{\max }\left\{p_i\right\}}---\left(23\right)$

Wherein D (x) representes the blood vessel distribution degree at horizontal level x place,

Be i section blood vessel connected component proportion in this position vertical window, m _iRepresent that the i section is communicated with the blood vessel pixel count, M representes the number of pixels that the vertical window blood vessel is total.

The number of pixels of representing largest connected vessel segment can to a certain degree be weighed the local density of blood vessel, n _xFor blood vessel total in this position vertical window is communicated with number.Generally speaking, the corresponding distribution degree in vertical window zone that comprises optic disk is less than the distribution degree of other horizontal level.The vertical window along continuous straight runs slides; Can obtain the blood vessel distribution degree value of each horizontal level; Thereby draw its one dimension drop shadow curve such as Fig. 7 (c) in the horizontal direction, no blood vessel or blood vessel pixel account for the blood vessel distribution degree that the vertical window area is less than 0.01% position and are preset as infinitely large quantity.

3. the minimum point of the one dimension drop shadow curve of horizontal direction can be confirmed as the horizontal coordinate x of optic disk _Od

5, locate based on the vertical coordinate of optic disk appearance characteristics.Because the brightness and the edge gradient value in optic disk zone are greater than other zones relatively, after precisely locating the horizontal coordinate of optic disk, can utilize this appearance characteristics further to locate the optic disk vertical coordinate.Specific practice is:

1. select the rectangular window W of optic disk size _v(length and width are the disc diameter size), and the fixed optic disk horizontal coordinate of the horizontal coordinate at its center x _OdThe place, ordinate y changes with the upper and lower slip of window;

2. rectangular window W _vRespectively former gray-scale map (Fig. 6 (a)) and

Vertically slide up and down on the phase place Gabor filter response figure (Fig. 6 (b)).And on the response diagram of

phase place, estimate mean flow rate and average filter response value in the window at former figure and Gabor filtering respectively; According to formula (24); Average brightness value and average filter response value are multiplied each other as the characteristic value f(y of longitudinal projection at y place, corresponding upright position)

f(y)＝I _N(x _od,y)*g _N(x _od,y) （24）

I wherein _N(x _Od, be y) with (x _Od, y) be the average brightness value of the rectangular window neighborhood at center, g _N(x _Od, be that Gabor filtering exists in this neighborhood y)

Average filter response estimation value during phase place has reflected the edge gradient information in the field.Because the brightness and the vessel boundary information in optic disk zone all are greater than other zone; Therefore the big position of f (y) value is that the possibility of optic disk is big; Equally when rectangular window slides up and down, can projection obtain a dimension curve of reflection f (y) situation of change of relevant position, like Fig. 7 (b).

3. the coordinate of maximizing point is the vertical coordinate y of optic disk in one dimension f (y) drop shadow curve _Od

The above only is a preferred implementation of the present invention, and protection scope of the present invention also not only is confined to the foregoing description, and all technical schemes that belongs under the thinking of the present invention all belong to protection scope of the present invention.Should be pointed out that for those skilled in the art some improvement and retouching not breaking away under the principle of the invention prerequisite should be regarded as protection scope of the present invention.

Claims (6)

1. the optic disk projection localization method of comprehensive vascular distribution and optic disk appearance characteristics may further comprise the steps:

(1) get the brightest pixel intensity value of red channel component in the original color eye fundus image 5% as threshold value; And go out the bianry image of foreground area according to this Threshold Segmentation; Bianry image is carried out the morphological erosion operation; Obtain the mask template, obtain the area-of-interest of eye fundus image through mask process;

(2), the area-of-interest after the mask process is carried out the normalization of illumination and contrast and handle the retina eye fundus image that is enhanced according to the image viewing model;

(3) the enhancing retina eye fundus image that step (2) is obtained adopts the non-blood vessel structure of Gabor filtering to suppress operator, and many threshold process technology that combine simultaneously to lag behind realize the eye fundus image vessel extraction and cut apart, and obtain blood vessel segmentation figure;

(4) vertical window is set, the height of window is a picture altitude, and width is 2 times a main blood vessel width, from left to right slides along the horizontal direction of blood vessel segmentation figure; In window, calculate the blood vessel distribution degree value D (x) at each horizontal level x place, draw the horizontal projection curve, find the minimum point of this dimension curve to confirm as the horizontal coordinate x of optic disk according to blood vessel distribution degree value _Od

(5) rectangular window is set, the length of window and width are the size of disc diameter, and the horizontal coordinate at its center is fixed optic disk horizontal coordinate x _Od, respectively in former gray-scale map and initial phase

Gabor filter response figure on moving window from top to bottom vertically; And interior mean flow rate of estimation window and average filter response value; Average brightness value and average filter response value are multiplied each other as longitudinal projection's value at y place, corresponding upright position; According to longitudinal projection's value paint projection curve, the coordinate of maximizing point is the vertical coordinate of optic disk.

2. the optic disk projection localization method of comprehensive vascular distribution according to claim 1 and optic disk appearance characteristics, the morphological erosion operation in the said step (1) is to adopt 9 * 9 square structure elements that bianry image is carried out erosion operation.

3. the optic disk projection localization method of comprehensive vascular distribution according to claim 1 and optic disk appearance characteristics, the concrete steps of said step (2) are:

1. image is divided into size and is the piecemeal S of s _i, s is M * N of 1/6～1/8, M * N is an original image size, to each piecemeal S _iCalculate the gray average in this piecemeal And standard deviation

Adopt the bicubic interpolation method to try to achieve in the entire image every average again

With standard deviation

2. pass through Manhabolios apart from d _M(x y), judges (whether x y) belongs to the background area to certain pixel, if it is apart from the neighborhood average

Less than certain threshold value t, threshold value t=1, promptly

$d_M(x,y)=\left|\frac{I(x,y)-{\widehat{\mathrm{\&mu;}}}_N(x,y)}{{\widehat{\mathrm{\&sigma;}}}_N(x,y)}\right|<t---\left(11\right)$

Then this pixel belongs to the background area, otherwise belongs to foreground area;

3. in the background area that splits illumination and contrast are estimated sampled point (s in each piecemeal _x, s _y) illumination And contrast

Average and standard deviation by background pixel point in the piecemeal are estimated to obtain, and the illumination of all the other non-sampled points and non-background pixel point and the contrast drift factor then obtain through the bicubic interpolation method;

4. according to formula (7), realize that the normalization of image strengthens,

$\widehat{I^0}(x,y)=\frac{I(x,y)-\hat{L}(x,y)}{\hat{C}(x,y)}---\left(7\right)$

I (x; Y) be that original observed image is at (x; Y) gray values of pixel points;

and

is respectively the illumination and the contrast drift factor of this some estimation, and

is the gray-scale value after the normalization.

4. the optic disk projection localization method of comprehensive vascular distribution according to claim 1 and optic disk appearance characteristics, said step (3) concrete steps are:

1. select for use 12 Gabor filtering templates on the different directions respectively image to be carried out filtering; Each adjacent direction phase difference of pi/12, filtered image is expression with

:

Wherein,

is the Gabor function; The f representing input images; Keep the peak response of each pixel at 12 different directions

$H\left(z\right)=\left\{\begin{array}{cc}0,& z<0\\ z,& z\&GreaterEqual;0\end{array}\right.---\left(15\right)$

2. be directed against the Gabor function when out of phase

; Detection characteristic to the different characteristic pattern; Employing formula (16) suppresses to handle to bright fringe and non-vessel boundary

$v_{\mathrm{\&lambda;},\mathrm{\&delta;}}(x,y)=H(r_{\mathrm{\&lambda;},\mathrm{\&delta;},\mathrm{\&pi;}}^{\mathrm{AH}}(x,y)-{\mathrm{\&alpha;}}_1{r}_{\mathrm{\&lambda;},\mathrm{\&delta;},0}^{\mathrm{AH}}(x,y)-{\mathrm{\&alpha;}}_2{r}_{\mathrm{\&lambda;},\mathrm{\&delta;},\mathrm{\&pi;}/2}^{\mathrm{AH}}(x,y))---\left(16\right)$

Wherein Expression respectively

0, the filter response during pi/2, α ₁, α ₂Represent inhibition component respectively to bright fringe and non-vessel boundary

The inhibition degree;

3. choose many group threshold values

and

and consider the connectedness of blood vessel; The two adjacent groups threshold interval should be overlapping; Promptly

selected threshold value suppresses to handle afterwards image to

and

to filtering and carries out list threshold binarization is cut apart, and obtains corresponding binary image

$b_{\mathrm{\&lambda;},\mathrm{\&delta;}}^i(x,y)=H{T}^i\left(v_{\mathrm{\&lambda;},\mathrm{\&delta;}}(x,y)\right),i\&Element;\{1,...,N_i\}---\left(17\right)$

HT representes the processing of standard hysteresis threshold;

4. each bianry image

is expressed as the union of some connected components

$b_{\mathrm{\&lambda;},\mathrm{\&delta;}}^i=\underset{k}{\&cup;}{C}_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{i,k}---\left(18\right)$

Wherein, k connected component in the binary image of

expression threshold value when being chosen for

and

.Utilize the morphology expansive working, select 2 * 2 square structure element q for use ₂, right

Each branch

Carry out expansion process, obtain:

$D_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{i,k}=C_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{i,k}\&CirclePlus;q_2---\left(19\right)$

5. final blood vessel structure V _{λ, δ}Can handle acquisition in regular turn through following formula,

$B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{(N_t,N_t-1)}=\underset{(b_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t}\&CirclePlus;q_2)\&cap;D_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t-1,k}\&NotEqual;\mathrm{\&phi;}}{\&cup;}(b_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t}\&cup;C_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t-1,k})---\left(20\right)$

$B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{(N_t-1,N_t-2)}=\underset{(B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{(N_t,N_t-1)}\&CirclePlus;q_2)\&cap;D_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t-2,k}\&NotEqual;\mathrm{\&phi;}}{\&cup;}(B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t,N_t-1}\&cup;C_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{N_t-2,k})---\left(21\right)$

.

.

.

$V_{\mathrm{\&lambda;},\mathrm{\&delta;}}=B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{\left(\mathrm{2,1}\right)}=\underset{(B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{\left(\mathrm{3,2}\right)}\&CirclePlus;q_2)\&cap;D_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{1,k}\&NotEqual;\mathrm{\&phi;}}{\&cup;}(B_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{\left(\mathrm{3,2}\right)}\&cup;C_{\mathrm{\&lambda;},\mathrm{\&delta;}}^{1,k})---\left(22\right)$

6. because blood vessel presents the network structure of connection in image, consider to remove the branch of all length less than 10 pixels.

5. the optic disk projection localization method of comprehensive vascular distribution according to claim 1 and optic disk appearance characteristics, the computing formula of described step (4) medium vessels distribution degree value D (x) is:

$D\left(x\right)=\frac{(-1)\underset{i=1}{\overset{n_x}{\mathrm{\&Sigma;}}}{p}_i{\log }_2{p}_i}{\underset{i}{\max }\left\{p_i\right\}}---\left(23\right)$

In the above-mentioned formula

Be i section blood vessel connected component proportion in this position vertical window, m _iRepresent that the i section is communicated with the blood vessel pixel count, M representes the total number of pixels of blood vessel in the horizontal slip window,

The number of pixels of representing largest connected vessel segment, n _xFor blood vessel total in this position level moving window is communicated with number.

6. the optic disk projection localization method of comprehensive vascular distribution according to claim 1 and optic disk appearance characteristics, the concrete steps of said step (5) are:

1. define window rectangular window W _v, W _vLength and width be the disc diameter size, and the fixed optic disk horizontal coordinate of the horizontal coordinate at its center x _OdThe place, ordinate y changes with the upper and lower slip of window;

2. moving window W from top to bottom vertically _vAnd on the response diagram of θ=90 ° phase place, estimate mean flow rate and average filter response value in the window at former figure and Gabor filtering respectively; According to formula (24), average brightness value and average filter response value are multiplied each other as longitudinal projection's value at y place, corresponding upright position

f(y)＝I _N(x _od,y)*g _N(x _od,y) （24）

I _N(x _Od, be y) with (x _Od, y) be the average brightness value of the rectangular window neighborhood at center, g _N(x _Od, y) be average Gabor filter response estimated value in this neighborhood;

3. draw longitudinal projection's curve according to f (y), the coordinate of maximizing point is the vertical coordinate y of optic disk _Od

Images