CN107657605A - A kind of Measurement of surface deepth method before sieve plate based on active profile and energy constraint - Google Patents

A kind of Measurement of surface deepth method before sieve plate based on active profile and energy constraint Download PDF

Info

Publication number
CN107657605A
CN107657605A CN201710810756.6A CN201710810756A CN107657605A CN 107657605 A CN107657605 A CN 107657605A CN 201710810756 A CN201710810756 A CN 201710810756A CN 107657605 A CN107657605 A CN 107657605A
Authority
CN
China
Prior art keywords
mrow
sieve plate
msub
pixel
bmo
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.)
Granted
Application number
CN201710810756.6A
Other languages
Chinese (zh)
Other versions
CN107657605B (en
Inventor
陈再良
彭鹏
沈海澜
刘晴
魏浩
段宣初
欧阳平波
廖胜辉
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN201710810756.6A priority Critical patent/CN107657605B/en
Publication of CN107657605A publication Critical patent/CN107657605A/en
Application granted granted Critical
Publication of CN107657605B publication Critical patent/CN107657605B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/0002Inspection of images, e.g. flaw detection
    • G06T7/0012Biomedical image inspection
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/11Region-based segmentation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/13Edge detection
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/50Depth or shape recovery
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10101Optical tomography; Optical coherence tomography [OCT]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30004Biomedical image processing
    • G06T2207/30041Eye; Retina; Ophthalmic

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Quality & Reliability (AREA)
  • Image Analysis (AREA)

Abstract

The invention discloses a kind of Measurement of surface deepth method before sieve plate based on active profile and energy constraint, this method is using surface segmentation method before the Bruch films opening point assay method based on k mean clusters and active profile and the sieve plate based on energy constraint, initial profile first using the dendrogram of k mean clusters as C V active contour models, extract the Bruch films opening point in profile, the area-of-interest of surface segmentation before sieve plate is obtained further according to the position of opening point, the preceding surface of sieve plate is partitioned into using energy constraint method, case depth before sieve plate is finally gone out according to the structure measurement of two steps.This method acquired results are better than existing method, it is and consistent with expert's craft calibration result, it can solve the problem that expert needs the problem of wasting time and energy of case depth before demarcation measurement sieve plate manually in clinical diagnosis, there is positive effect to the early screening and clinical diagnosis of glaucoma.

Description

A kind of Measurement of surface deepth method before sieve plate based on active profile and energy constraint
Technical field
The invention belongs to image identification technical field, is related to a kind of preceding surface of the sieve plate based on active profile and energy constraint Depth measurement method.
Background technology
Glaucoma is second-biggest-in-the-world blinding disease, and it to the peripapillary optic ganglion cell aixs cylinder of people by breaking It is bad, cause the visual field of people to lack.Due to the irreversibility of glaucoma, morning detection, early discovery and early treatment to glaucoma can Slow down the process of disease development.But because the pathogenesis of glaucoma is not yet completely clearly, therefore to the risk factors of glaucoma Research be still current hot issue.
Lamina Cribrosa (lamina cribrosa) is located at sclera, is that retinal ganglion cell axon passes eyes Position, in sieve texture, it is responsible for providing nutrition and structural support to retinal ganglion cell axon, and help to maintain intraocular And outside eye between barometric gradient.The anatomical structure of sieve plate is as follows:Ten thousand retinal ganglion cell axons of about 1.2-1.5 It is gathered in and regards nipple, eyes is partly left by the inside (i.e. Bruch's films opening) and outside (i.e. sclera) of canalis opticus.Depending on Nerve is a part for central nervous system, and outside is surrounded by three layers of envelope (endocranium, the spider web continued by three layers of meninx Film, pia mater), endocranium cavity of resorption and cavum subarachnoidale are also extended to around optic nerve therewith, and cavum subarachnoidale is full of cerebrospinal fluid. At sclera metapore, 1/3 forms thin layer fiber every there is many sieve apertures to allow retinal ganglial cellses in dividing plate before scleral tissue Aixs cylinder is walked, and rear 2/3 connects with endocranium.Research shows the sieve plate risk factors maximum with glaucoma --- intraocular pressure has one Determine relation, would generally be along with the rise of intraocular pressure in the pathogenic process of glaucoma, and the rise of intraocular pressure can cause simultaneously The rearward displacement of sieve plate and flexural deformation etc..Wang Ning profits et al. propose optic ganglion cell axonal injury and across sieve plate pressure differential (difference of intraocular pressure and intracerebral pressure) is relevant, and sieve plate can also show the feature of rearward displacement under this pressure differential.
Clinic of means of optical coherence tomography (Optic Coherence Tomography, the OCT) technology in ophthalmology Using only 20 years, its technological innovation was rapid, now had become one of most important inspection of clinical ophthalmology.It is right that the technology passes through Tissue emissions coherent light, reclaims the reflected light and scattering light of tissue, and is combined with the delay of time, can get tissue The information of two-dimentional fault information or overall retina three-dimensional image.Except in real time monitoring, it is noninvasive the advantages that outside, OCT is most important Feature is its high resolution, and the fine structure observed is cross-sectional structure, meets pathological normal observation custom, is The live body retinal morphology research of ophthalmology researcher provides technical support.The OCT course of work is as follows:Sent by low-coherence light source Low coherence light, two-beam line is divided into by interferometer, it is a branch of to enter detection light path, direct incident intraocular, by intraocular difference group The interface knitted fires back, to provide the thickness and range information of the various tissues of intraocular;Another beam enters with reference to light path, by Know that the reference mirror of space length reflects.Two-beam is integrated into one in optical fiber coupler, due to being transmitted into coming for reference mirror The distance returned distance and be transmitted into the given structure of intraocular accurately matches, therefore produces interference, is detected by light-sensitive detector.Adjust Signal input computer after solution carries out computing, obtains the optical coherence tomography image of testee.Because eye inner tissue has Different depth and space structure, so meeting generation time is poor between two-beam line, this time difference is called the optical delay time. The time difference that detector is calculated is related to tolerance principle using low-coherent light, the information of Tissue reflectance can be obtained.Get After these information, one-dimensional scan image information is calculated by computer, usually the row in two dimensional image.
In recent years, the progress of OCT technology gradually solves the problems, such as that deep tissues can not be imaged.Heidelberg company OCT's Strengthen Depth Imaging (Enhanced Depth Imaging, EDI) pattern can by light collection in the deeper position on eyeground, Such as sclera, choroid and sieve plate.This imaging technique have beneficial to the pathologic of live body sieve plate change such as sieve plate to Backward shift etc. is further studied, and improves the accuracy of glaucoma early diagnosis.We are usually using case depth before sieve plate To represent the rearward displacement of sieve plate, the measuring method of this sieve plate depth is:With Bruch films opening point (Bruch ' s Membrane Opening, BMO) composition BMO reference planes (sieved with the preceding surface of sieve plate as reference plane, measurement datum The distance between interface between being organized before plate tissue and sieve plate).Therefore, the measurement of case depth is divided into two steps before sieve plate: BMO is determined and surface segmentation before sieve plate.After obtaining the position on BMO reference planes and the preceding surface of sieve plate, case depth can before sieve plate To be drawn by measuring the distance between two faces.
Because the research to sieve plate just starts to walk, some other existing methods can not also be advantageously applied to sieve plate figure In handling, so related achievement in research is relatively fewer.
A) BMO is determined
BMO points are one of most important biomarkers in OCT image, and it represents optic disk edge in OCT faultage images Position, therefore it is widely used in optic disk segmentation, disk is the problems such as measurement along.
2013, Fu et al. proposed a kind of BMO assay methods rebuild based on dictionary learning and low-rank matrix, this method Be partitioned into first retina internal limiting membrane (ILM) layer and retinal pigment epithelium (RPE) layer and divide the image into training region and Candidate region, low-rank dictionary then is trained with the training region of image, and goes out whole image with the data reconstruction of candidate region, Obtain reconstruction error curve, including luminance errors, local binary patterns (LBP) error and apart from biasing.Finally by three mistakes Poor curve draws the position of the position i.e. BMO points of optic disk.The shortcomings that this method, is the BMO that can not effectively determine glaucoma lesion Position.
2014, Akram et al. is proposed a kind of determined BMO position based on the model of deconvolution.The model is thought Each layer in OCT image can be modeled by two compositions, be respectively:With a curve come model the skeleton of layer and with one or One group of wave filter models the thickness of layer.As long as the parameter of the two compositions is obtained, it is possible to uniquely determine every on OCT image One layer of position.The model is entered using Monte Carlo Markov Chain (Monte Carlo Markov Chain, MCMC) method Row parametric measurement, reach higher accuracy rate, but also bring the problem of less efficient simultaneously.
2015, Wang et al. proposed a kind of BMO assay method, the method use five kinds of methods and determines respectively BMO position, shape facility of the first method based on ILM layers, texture of remaining four kinds of method based on BMO neighbouring positions are special Sign.After obtaining the result of five kinds of methods, select wherein three immediate values and ask their average value to be surveyed as final BMO Determine result.The problems such as less efficient, robustness is poor be present in this method.
Hussian et al. proposes a kind of BMO assay methods based on graph theory, and this method measures three reference layers first Position, it is retinal nerve fiber (RNFL) layer, external plexiform layer (ONL) layer and RPE layers respectively.Then the position using layer and spy The weight as figure is levied, BMO position is drawn using a kind of graph search algorithm, has reached higher accuracy rate.
B) surface segmentation before sieve plate
The preceding surface of sieve plate is sieve plate tissue and the interface organized before sieve plate, not solid because the sieve plate of people comes in every shape Fixed shape facility, contrast unobvious sometimes can also be traditional based on gradient because the disease such as glaucoma causes to lack Feature, the dividing method effect of color and shape and bad, some existing Hierarchical Segmentation methods are not suitable for table before sieve plate yet In the segmentation of face.
2015, Akram et al. proposed the automatic division method on the preceding surface of the first sieve plate, and this method is based on Ma Erke Husband's random field and MCMC sampling algorithms, classification problem is converted into by segmentation problem, it is only necessary to knows whether each point belongs to sieve The preceding surface of plate.In units of each point, extract the feature of the point and its neighbours' point establish a Markov with Airport, the parameter of markov random file is then determined using MCMC algorithms, and sample by iteration to obtain final surface.
However, above-mentioned existing method generally exist it is less efficient, to glaucoma lesion image not robust the problems such as, cause it Practicality in the examination and clinical diagnosis of glaucoma it is poor.
The content of the invention
The invention provides a kind of Measurement of surface deepth method before sieve plate based on active profile and energy constraint, its purpose It is, overcomes the problem of efficiency is low in the prior art and robustness is not strong.
A kind of Measurement of surface deepth method before sieve plate based on active profile and energy constraint, comprises the following steps:
Step 1:After gray level image to be detected is carried out into inverse processing, all pixels point in image is handled to inverse and adopted Clustered with k means clustering methods, obtain dendrogram;
Step 2:The energy function of C-V active contour models is constructed, C-V active profile dies are used as by the use of the profile of dendrogram The initial curve of type, using the curve during energy function minimums of C-V active contour models as objective contour, with objective contour The terminating point of lower boundary obtains BMO and refers to projection line as BMO points, two BMO points of connection;
Step 3:Handle and chosen in image positioned at BMO with reference to the region below projection line from reflection, it is corresponding with selected areas Area-of-interest of the minimum rectangular area as surface segmentation before sieve plate, wherein, the coboundary of selected areas is BMO with reference to throwing The projection line in the horizontal direction of hachure;
Step 4:The preceding surface candidate point of one sieve plate, and foundation are extracted from area-of-interest based on energy function in each row Integrity constraint, rejecting processing is carried out to all preceding surface candidate points of sieve plate, obtain controlling point set;Control point set is entered Row curve matching obtains floor lever of sieve tray contour line;
Step 5:Projection line and floor lever of sieve tray contour line are referred to based on BMO, calculate case depth before sieve plate;
Using BMO with reference to the midpoint of projection line and its nasal side and two points of each 100 microns of temporo side as measurement point, at three Vertical line is made with reference to projection line to BMO in measurement point and intersected with surface profile line before sieve plate, the average value of three length of perpendicular is Case depth before sieve plate to be measured.
Further, the dendrogram is M (i, j, k), is to handle inverse to belong to the minimum cluster of pixel value in image The pixel of the class of center representative retains, and the gray value of rest of pixels point is set to 0:
Wherein, cmAnd ckM-th and k-th of cluster centre, and c are represented respectivelymRepresent the current cluster of pixel (i, j) Center, I (i, j) be image in pixel (i, j) pixel value, I (cm) and I (ck) it is respectively cluster centre cmAnd ckPlace Pixel value.
Further, the energy function of the C-V active contour models is:
E (C)=μ L (C)+v*Area (inside (C))+λ ∫inside(C)(I(x,y)-Ii)2dxdy
+ω∫outside(C)(I(x,y)-Io)2dxdy
Wherein, I (x, y) is the pixel value of the pixel (x, y) in dendrogram, and L (C) represents length of a curve, Area (inside (C)) represents the area that contour curve is surrounded, IiAnd IoThe respectively pixel of outer and in profile the pixel of profile Value;Inside (C) represents the pixel in profile, and outside (C) represents the pixel outside profile in dendrogram;μ joins for length Number, span is (0,1);V is area parameters, value 0;λ and ω is respectively the weight coefficient of self-energy and outer energy, λ =ω=1.
Further, according to following during the energy function based on surface segmentation before sieve plate respectively arranges from area-of-interest The preceding surface candidate point of one sieve plate of formulas Extraction:
Wherein, SjThe preceding surface candidate point of sieve plate in jth row is represented, n represents total columns of pixel in area-of-interest;Ij (i) the brightness value of the ith pixel point of jth row pixel in area-of-interest, gradient (I are representedj(i) I) is representedj(i) Gradient Features value, α and β represent the brightness of pixel and the weight parameter of Gradient Features respectively, and α span is [0.3,0.4], β span is [0.5,0.6], and arg max () represent to find the position for the point for making energy function maximum.
Total columns of pixel is the width of region of interest area image in area-of-interest;
Arg max () are to cause α * Ij(i)+β*gradient(Ij(i) when) obtaining maximum, corresponding pixel is made For Sj
Further, the integrity constraint is:
Wherein, Control (j) represents control point set, and d represents offset distance control parameter, and span is (10,15);S For whole candidate point ordinate set, YjThe ordinate of the candidate point on j-th of A-scan is represented, average (S) waits to be whole The average value of reconnaissance set ordinate.
Further, using B-spline fitting algorithm point set will be controlled to carry out curve fitting to obtain floor lever of sieve tray contour line.
Beneficial effect
The present invention proposes a kind of Measurement of surface deepth method before sieve plate based on active profile and energy constraint, this method Based on active contour model, extract the BMO points in OCT image and BMOBMO refers to projection line, reuse energy constraint method It is partitioned into the preceding surface of sieve plate.Finally using the result of surface segmentation before BMO measure and sieve plate, case depth before sieve plate is measured. This method has used the active contour model based on region in BMO points measure, avoids making measure using traditional Gradient Features As a result by the abnormal influence such as blood vessel shade, lesion;Before sieve plate in surface segmentation, this method employs the think of of energy function Want solve the problem of preceding surface of sieve plate is without the fixed traditional characteristic such as shape, size, and sieve is fitted using curve matching The lack part of plate;Test result indicates that the accuracy rate peace of this method surface segmentation before the BMO accuracy determined and sieve plate Existing method, the result more demarcated manually close to expert are superior in slip.The present invention can solve the problem that expert in clinical diagnosis When need manually demarcation measurement sieve plate before the problem of wasting time and energy of case depth, early screening and clinical diagnosis to glaucoma With positive effect.
Brief description of the drawings
Fig. 1 is the flow chart of the present invention, wherein, (a) is its logical flow chart, and (b) is that it performs the flow chart of step;
Fig. 2 is the schematic diagram of Measurement of surface deepth before sieve plate;
Fig. 3 is the BMO measurement results of this method, wherein, (a) is measurement result Fig. 1 of this method;(b) to be corresponding with (a) The result demarcated by hand of expert;(c) it is measurement result Fig. 2 of this method;(d) expert corresponding with (c) demarcates by hand As a result;(e) it is measurement result Fig. 3 of this method;(f) result demarcated by hand for expert corresponding with (e);
Fig. 4 is surface segmentation result before the sieve plate of this method;
The contrast schematic diagram that Fig. 5 is demarcated by hand for case depth before the sieve plate of this method measurement with expert, wherein, (a) is Measurement result Fig. 1 of this method;(b) result demarcated by hand for expert corresponding with (a);(c) it is the measurement result of this method Fig. 2;(d) result demarcated by hand for expert corresponding with (c);(e) it is measurement result Fig. 3 of this method;(f) it is right with (e) The result that the expert answered demarcates by hand.
Embodiment
Below in conjunction with accompanying drawing, the present invention will be further described.
The present invention proposes a kind of Measurement of surface deepth method before sieve plate based on active profile and energy constraint, including four Individual step:First with initial profile of the dendrogram of k mean algorithms as C-V active contour models, then contouring lower boundary Terminating point, as BMO points, then extract in image comprising BMO with reference to projection line region below minimum rectangle as sieve The area-of-interest of surface segmentation before plate, candidate point is extracted according to energy function to each row of interesting image regions, will The candidate point for not meeting constraints is left out, and is partitioned into whole surface with curve-fitting method, finally calculates table before sieve plate Face depth.
Shown in idiographic flow such as Fig. 1 (a), shown in implementation procedure such as Fig. 1 (b) of this method, its specific implementation step is as follows:
Step 1:After gray level image to be detected is carried out into inverse processing, all pixels point in image is handled to inverse and adopted Clustered with k means clustering methods, obtain dendrogram;
Image to be detected I (m, n) after given inverse, k mean algorithms are divided into m × n pixel in image predefined K class, it is smaller according to the otherness of pixel in each class class, and the principle that otherness between class and class is larger, Wo Menke The problem of making cost function (representing the otherness between pixel) minimum so that problem is converted into, then obtain every a kind of cluster Center so that the otherness of each point and nearest cluster centre is minimum.Used cost function is as follows:
Wherein c is predefined classification number, d (i, j, k)=| | I (i, j)-I (ck)||2For two points on description image it Between diversity factor metric function because OCT image has level feature, therefore use the pixel value tag of image as difference Property measurement, I (i, j) be image in pixel, ckFor cluster centre.So make the minimum dendrogram M (i, j, k) of cost function It is calculated according to lower formula.
Wherein, cmAnd ckM-th and k-th of cluster centre, and c are represented respectivelymRepresent the current cluster of pixel (i, j) Center, I (i, j) be image in pixel (i, j) pixel value, I (cm) and I (ck) it is respectively cluster centre cmAnd ckPlace Pixel value.The pixel for belonging to the class that the minimum cluster centre of pixel value represents in figure is retained, remaining sets to 0, as to be measured Determine dendrogram.After obtaining dendrogram, the central point per a kind of all pixels is exactly new cluster centre.
Step 2:The energy function of C-V active contour models is constructed, C-V active profile dies are used as by the use of the profile of dendrogram The initial curve of type, using the curve during energy function minimums of C-V active contour models as objective contour, with objective contour The terminating point of lower boundary obtains BMO and refers to projection line as BMO points, two BMO points of connection;
Because the result of active contour model is largely influenceed by given initial curve, therefore by upper one Step obtained by k mean clusters figure as initial curve C, the C-V active contour model of active contour model energy function such as Following formula:
E (C)=μ L (C)+v*Area (inside (C))+Einside+Eoutside
=μ L (C)+v*Area (inside (C))+λ ∫inside(C)(I(x,y)-Ii)2dxdy
+ω∫outside(C)(I(x,y)-Io)2dxdy
Wherein, I (x, y) is the pixel value of the pixel (x, y) in dendrogram, and L (C) represents length of a curve, Area (inside (C)) represents the area that contour curve is surrounded, IiAnd IoThe respectively pixel of outer and in profile the pixel of profile Value;Inside (C) represents the pixel in profile, and outside (C) represents the pixel outside profile in dendrogram;μ joins for length Number, its value is determined by the size of target in image, if the size of target is larger, μ value is also larger, otherwise μ Value will be smaller, span 0-1, generally takes μ=0.5.V is area parameters, and λ and ω are respectively self-energy and outer energy Weight coefficient.Generally take λ=ω=1, v=0.C-V models obtain the distribution of its energy by constructing energy function, with energy Function ever-reduced method promotes curve to approach the edge of target, and most target is split from background at last.It breaks away from Dependence in classical edge partitioning algorithm such as LOG and traditional Snake models to image gradient, have to OCT image good Segmentation ability.
Step 3:Handle and chosen in image positioned at BMO with reference to the region below projection line from reflection, it is corresponding with selected areas Area-of-interest of the minimum rectangular area as surface segmentation before sieve plate, wherein, the coboundary of selected areas is BMO with reference to throwing The projection line in the horizontal direction of hachure;
After obtaining positions of the BMO with reference to projection line, area-of-interest is used as with reference to the floor projection line of projection line using BMO The width of coboundary, i.e. rectangular area;Region lower boundary is image base, and right boundary is respectively the position of two BMO points.
Step 4:The preceding surface candidate point of one sieve plate, and foundation are extracted from area-of-interest based on energy function in each row Integrity constraint, rejecting processing is carried out to all preceding surface candidate points of sieve plate, obtain controlling point set;Control point set is entered Row curve matching obtains floor lever of sieve tray contour line;
Brightness and Gradient Features around the preceding surface of sieve plate build the energy function of surface segmentation before sieve plate, and The preceding surface candidate point of a sieve plate, such as following formula are selected according to energy function in each A-scan (i.e. each row, similarly hereinafter):
Wherein, SjThe preceding surface candidate point of sieve plate in jth row is represented, n represents total columns of pixel in area-of-interest;Ij (i) the brightness value of the ith pixel point of jth row pixel in area-of-interest, gradient (I are representedj(i) I) is representedj(i) Gradient Features value, α and β represent the brightness of pixel and the weight parameter of Gradient Features respectively, and α span is [0.3,0.4], β span is [0.5,0.6], and arg max () represent to find the position for the point for making energy function maximum Put.After obtaining candidate point set, because the morphological feature of sieve plate varies with each individual, and the influence of blood vessel shade, there are some times Reconnaissance is in fact, it is necessary to find out and deleted the common trait of this part candidate point not on the position on the preceding surface of sieve plate Go.Integrity constraint is as follows:
Wherein, Control (j) represents control point set, and d represents offset distance control parameter, is needed with this state modulator The excessive candidate point of the deviant of removal, value is between 10-15.S is whole candidate point ordinate set, YjRepresent j-th The ordinate of candidate point on A-scan, average (S) are the average value of whole candidate point set ordinate.In order to ensure Smooth surface before the sieve plate arrived, we are carried out curve fitting on obtained control point set Control (j), and the fitting used is calculated Method is B-spline fitting algorithm.
Step 5:Projection line and floor lever of sieve tray contour line are referred to based on BMO, calculate case depth before sieve plate;
Obtain before BMO measurement results and sieve plate after surface segmentation result, with BMO with reference to the midpoint of projection line and its nasal side and Two points of each 100 microns of temporo side as measurement point, in three measurement points to BMO with reference to projection line make vertical line and with before sieve plate Intersect on surface.Such as Fig. 2, the average value of three length of perpendicular is case depth before sieve plate to be measured.
To verify the performance of method proposed by the present invention, 30 width images of 18 samples are gathered as data set, and with special Ground Truth figures corresponding to the preceding surface conduct of BMO points and sieve plate that family demarcates by hand.Based on described in step 3, α and β parameters Whether effective method has been largely fixed it, has been found through experiments that, when α increases, curve can be close to sieve plate region;When β increases When big, curve can be close to internal limiting membrane region.For the robustness of ensuring method, α=0.3-0.4, β=0.5-0.6 are set.
The BMO assay methods and NLSC sieve plates that the present invention proposes the computational methods of proposition and Hussian in 2015 et al. Preceding surface segmentation method is compared.In BMO measure, we are with the mean error and variance of measuring point and Ground Truth As evaluation criterion;Before sieve plate in surface segmentation, we are then using error rate as evaluation criterion.
Table 1BMO measurement results
From table 1 it follows that method proposed by the present invention (48.17 microns of mean error) in accuracy is better than The BMO assay methods (54.18 microns of mean error) that Hussian et al. is proposed, and in stability, the standard deviation of this method It is 51.32 microns, is also slightly better than its 53.74 microns.Fig. 3 lists this method and obtains segmentation result and expert's manual segmentation BMO points and BMO reference planes, as can be seen from Figure 3 due to avoiding Gradient Features, the present invention can be very good to avoid by big The problem of image blurring that blood vessel shadow band comes, can also obtain accurately splitting very much on the EDI-OCT images more than blood vessel shade As a result, it is and consistent with the result that expert demarcates manually.
(error be present with Ground Truth in the error rate before sieve plate in surface segmentation using the preceding surface sampling point of sieve plate Point the ratio between number and total sample number) evaluation criterion is used as, error rate grade is defined as from 0 to 2:It is small that grade 0 represents mean error The ratio between total sample number is accounted in the point of 3 pixels, grade 1 represents point of the mean error between 3 pixels and 5 pixels and accounts for sample The ratio between sum, grade 2 represent mean error and account for the ratio between total sample number more than the point of 5 pixels.By point of proposition method of the present invention Cut result and a kind of BFPS edges dividing method and the preceding surface of currently the only sieve plate based on Markov chain Monte-Carlo Automatic division method --- NLSC methods are contrasted, as a result as shown in table 2.
Surface segmentation result before the sieve plate of table 2
100-150 evaluation point is up-sampled in every pictures.As can be seen from Table 2 compared with Ground Truth, In error rate grade 0 and grade 1, method proposed by the invention respectively reaches 76.5% and 17.6%, higher than NLSC methods 73.7% and 16.1% and the 40.9% of BFPS methods and 22.8%, and in the higher error rate grade 2 of error, institute of the present invention Moving party's rule only has 5.9%, is better than the 36.3% of the 10.2% and BFPS methods of NLSC methods.Institute of the present invention as can be seen here The method of proposition is better than NLSC methods and BFPS methods on the whole.Fig. 4 lists pair of experimental result and Ground Truth Than, it can be seen that the result that experimental result is demarcated manually with expert is very close on image, and curve is also very smooth, from Also there is reasonability from the point of view of in histology.
Sieve plate depth measured by last the inventive method is compared with Ground Truth, still using mean error It is as shown in table 3 as evaluation criterion, its result with standard deviation.
The preceding surface automatic measurement result of the sieve plate of table 3
It can be seen that, the sieve plate depth and the resultant error of expert's demarcation that the present invention measures are smaller from table 3, basic one Causing, Fig. 5 gives more intuitively result, it can also be seen that the craft demarcation of this method and expert are basically identical from figure, and Accuracy rate is higher, and it is 22.1 microns averagely to have error in label, and no error in label is 28.7 microns.
Specific embodiment described herein is only to spirit explanation for example of the invention.Technology belonging to the present invention is led The technical staff in domain can be made various modifications or supplement to described specific embodiment or be replaced using similar mode Generation, but without departing from the spiritual of the present invention or surmount scope defined in appended claims.

Claims (6)

1. a kind of Measurement of surface deepth method before sieve plate based on active profile and energy constraint, it is characterised in that including following Step:
Step 1:After gray level image to be detected is carried out into inverse processing, all pixels point in image is handled to inverse and uses k Means clustering method is clustered, and obtains dendrogram;
Step 2:The energy function of C-V active contour models is constructed, C-V active contour models are used as by the use of the profile of dendrogram Initial curve, using the curve during energy function minimums of C-V active contour models as objective contour, with the following of objective contour The terminating point on boundary obtains BMO and refers to projection line as BMO points, two BMO points of connection;
Step 3:Handle and chosen in image positioned at BMO with reference to the region below projection line from reflection, with corresponding to selected areas most Area-of-interest of the small rectangular area as surface segmentation before sieve plate, wherein, the coboundary of selected areas refers to projection line for BMO Projection line in the horizontal direction;
Step 4:The preceding surface candidate point of one sieve plate is extracted in each row from area-of-interest based on energy function, and according to complete Property constraints, rejecting processing is carried out to all preceding surface candidate points of sieve plate, obtain control point set;Point set march will be controlled Line is fitted to obtain floor lever of sieve tray contour line;
Step 5:Projection line and floor lever of sieve tray contour line are referred to based on BMO, calculate case depth before sieve plate;
Using BMO with reference to the midpoint of projection line and its nasal side and two points of each 100 microns of temporo side as measurement point, in three measurements Vertical line is made with reference to projection line to BMO on point and intersected with surface profile line before sieve plate, the average value of three length of perpendicular is to be measured Case depth before the sieve plate of amount.
2. according to the method for claim 1, it is characterised in that the dendrogram is M (i, j, k), is to scheme inverse processing The pixel for belonging to the class that the minimum cluster centre of pixel value represents as in retains, and the gray value of rest of pixels point is set to 0:
<mrow> <mi>M</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>I</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>i</mi> <mi>f</mi> <mo>|</mo> <mo>|</mo> <mi>I</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>-</mo> <mi>I</mi> <mrow> <mo>(</mo> <msub> <mi>c</mi> <mi>m</mi> </msub> <mo>)</mo> </mrow> <mo>|</mo> <msub> <mo>|</mo> <mn>2</mn> </msub> <mo>&amp;le;</mo> <mo>|</mo> <mo>|</mo> <mi>I</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>-</mo> <mi>I</mi> <mrow> <mo>(</mo> <msub> <mi>c</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mo>|</mo> <msub> <mo>|</mo> <mn>2</mn> </msub> <mo>,</mo> <mi>m</mi> <mo>&amp;NotEqual;</mo> <mi>k</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mrow> <mi>o</mi> <mi>t</mi> <mi>h</mi> <mi>e</mi> <mi>r</mi> <mi>s</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow>
Wherein, cmAnd ckM-th and k-th of cluster centre, and c are represented respectivelymThe current cluster centre of pixel (i, j) is represented, I (i, j) be image in pixel (i, j) pixel value, I (cm) and I (ck) it is respectively cluster centre cmAnd ckThe pixel at place Value.
3. according to the method for claim 2, it is characterised in that the energy function of the C-V active contour models is:
E (C)=μ L (C)+v*Area (inside (C))+λ ∫inside(C)(I(x,y)-Ii)2dxdy+ω∫outside(C)(I(x,y)- Io)2dxdy
Wherein, I (x, y) is the pixel value of the pixel (x, y) in dendrogram, and L (C) represents length of a curve, Area (inside (C) area that contour curve is surrounded, I) are representediAnd IoThe respectively pixel value of outer and in profile the pixel of profile;inside (C) pixel in profile is represented, outside (C) represents the pixel outside profile in dendrogram;μ is length parameter, value model Enclose for (0,1);V is area parameters, value 0;λ and ω is respectively the weight coefficient of self-energy and outer energy, λ=ω=1.
4. according to the method for claim 1, it is characterised in that the energy function based on surface segmentation before sieve plate is from sense In interest region the preceding surface candidate point of a sieve plate is extracted in each row according to below equation:
<mrow> <msub> <mi>S</mi> <mi>j</mi> </msub> <mo>=</mo> <mi>arg</mi> <munder> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> </mrow> <mi>i</mi> </munder> <mrow> <mo>(</mo> <mi>&amp;alpha;</mi> <mo>*</mo> <msub> <mi>I</mi> <mi>j</mi> </msub> <mo>(</mo> <mi>i</mi> <mo>)</mo> <mo>+</mo> <mi>&amp;beta;</mi> <mo>*</mo> <mi>g</mi> <mi>r</mi> <mi>a</mi> <mi>d</mi> <mi>i</mi> <mi>e</mi> <mi>n</mi> <mi>t</mi> <mo>(</mo> <mrow> <msub> <mi>I</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mo>,</mo> <mi>j</mi> <mo>=</mo> <mn>1</mn> <mo>...</mo> <mi>n</mi> </mrow>
Wherein, SjThe preceding surface candidate point of sieve plate in jth row is represented, n represents total columns of pixel in area-of-interest;Ij(i) table Show the brightness value of the ith pixel point of jth row pixel in area-of-interest, gradient (Ij(i) I) is representedj(i) ladder Spend characteristic value, α and β represent the brightness of pixel and the weight parameter of Gradient Features respectively, α span for [0.3, 0.4], β span is [0.5,0.6], and argmax () represents to find the position for the point for making energy function maximum.
5. according to the method for claim 4, it is characterised in that the integrity constraint is:
<mrow> <mi>C</mi> <mi>o</mi> <mi>n</mi> <mi>t</mi> <mi>r</mi> <mi>o</mi> <mi>l</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>S</mi> <mi>e</mi> <mi>t</mi> <mo>{</mo> <mi>j</mi> <mo>|</mo> <mfrac> <mrow> <msub> <mi>dY</mi> <mi>j</mi> </msub> </mrow> <mrow> <mi>d</mi> <mi>j</mi> </mrow> </mfrac> <mo>*</mo> <mfrac> <mrow> <msub> <mi>dY</mi> <mrow> <mi>j</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> </mrow> <mrow> <mi>d</mi> <mi>j</mi> </mrow> </mfrac> <mo>&amp;le;</mo> <mn>0</mn> <mo>,</mo> <mo>|</mo> <msub> <mi>Y</mi> <mi>j</mi> </msub> <mo>-</mo> <mi>a</mi> <mi>v</mi> <mi>e</mi> <mi>r</mi> <mi>a</mi> <mi>g</mi> <mi>e</mi> <mrow> <mo>(</mo> <mi>S</mi> <mo>)</mo> </mrow> <mo>|</mo> <mo>&lt;</mo> <mi>d</mi> <mo>}</mo> <mo>,</mo> <mi>j</mi> <mo>=</mo> <mn>1</mn> <mo>...</mo> <mi>n</mi> </mrow>
Wherein, Control (j) represents control point set, and d represents offset distance control parameter, and span is (10,15);S is whole Individual candidate point ordinate set, YjThe ordinate of the candidate point on j-th of A-scan is represented, average (S) is whole candidate point Gather the average value of ordinate.
6. according to the method described in claim any one of 1-5, it is characterised in that will control point set using B-spline fitting algorithm Carry out curve fitting to obtain floor lever of sieve tray contour line.
CN201710810756.6A 2017-09-11 2017-09-11 A kind of sieve plate front surface depth measurement method based on active profile and energy constraint Active CN107657605B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710810756.6A CN107657605B (en) 2017-09-11 2017-09-11 A kind of sieve plate front surface depth measurement method based on active profile and energy constraint

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710810756.6A CN107657605B (en) 2017-09-11 2017-09-11 A kind of sieve plate front surface depth measurement method based on active profile and energy constraint

Publications (2)

Publication Number Publication Date
CN107657605A true CN107657605A (en) 2018-02-02
CN107657605B CN107657605B (en) 2019-12-03

Family

ID=61128062

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710810756.6A Active CN107657605B (en) 2017-09-11 2017-09-11 A kind of sieve plate front surface depth measurement method based on active profile and energy constraint

Country Status (1)

Country Link
CN (1) CN107657605B (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109528155A (en) * 2018-11-19 2019-03-29 复旦大学附属眼耳鼻喉科医院 A kind of intelligent screening system and its method for building up suitable for the concurrent open-angle glaucoma of high myopia
CN110288540A (en) * 2019-06-04 2019-09-27 东南大学 A kind of online imaging standards method of carbon-fibre wire radioscopic image

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103236065A (en) * 2013-05-09 2013-08-07 中南大学 Biochip analysis method based on active contour model and cell neural network
EP2967317A2 (en) * 2013-03-14 2016-01-20 Carl Zeiss Meditec AG Multimodal integration of ocular data acquisition and analysis
CN106373168A (en) * 2016-11-24 2017-02-01 北京三体高创科技有限公司 Medical image based segmentation and 3D reconstruction method and 3D printing system

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2967317A2 (en) * 2013-03-14 2016-01-20 Carl Zeiss Meditec AG Multimodal integration of ocular data acquisition and analysis
CN103236065A (en) * 2013-05-09 2013-08-07 中南大学 Biochip analysis method based on active contour model and cell neural network
CN106373168A (en) * 2016-11-24 2017-02-01 北京三体高创科技有限公司 Medical image based segmentation and 3D reconstruction method and 3D printing system

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109528155A (en) * 2018-11-19 2019-03-29 复旦大学附属眼耳鼻喉科医院 A kind of intelligent screening system and its method for building up suitable for the concurrent open-angle glaucoma of high myopia
CN110288540A (en) * 2019-06-04 2019-09-27 东南大学 A kind of online imaging standards method of carbon-fibre wire radioscopic image
CN110288540B (en) * 2019-06-04 2021-07-06 东南大学 Carbon fiber wire X-ray image online imaging standardization method

Also Published As

Publication number Publication date
CN107657605B (en) 2019-12-03

Similar Documents

Publication Publication Date Title
Bogunović et al. RETOUCH: The retinal OCT fluid detection and segmentation benchmark and challenge
Kafieh et al. A review of algorithms for segmentation of optical coherence tomography from retina
US10194866B2 (en) Methods and apparatus for reducing artifacts in OCT angiography using machine learning techniques
Baroni et al. Towards quantitative analysis of retinal features in optical coherence tomography
Sigal et al. Eye-specific IOP-induced displacements and deformations of human lamina cribrosa
WO2022007352A1 (en) Three-dimensional choroidal vessel imaging and quantitative analysis method and apparatus based on optical coherence tomography system
Norman et al. Dimensions of the human sclera: thickness measurement and regional changes with axial length
Esmaeelpour et al. Choroidal Haller's and Sattler's layer thickness measurement using 3-dimensional 1060-nm optical coherence tomography
US7992999B2 (en) Automated assessment of optic nerve head with spectral domain optical coherence tomography
Knott et al. Spatial correlation of mouse photoreceptor-RPE thickness between SD-OCT and histology
Xu et al. Automated geographic atrophy segmentation for SD-OCT images based on two-stage learning model
Guo et al. A framework for classification and segmentation of branch retinal artery occlusion in SD-OCT
Zhang et al. A novel technique for robust and fast segmentation of corneal layer interfaces based on spectral-domain optical coherence tomography imaging
Hussain et al. An automated method for choroidal thickness measurement from Enhanced Depth Imaging Optical Coherence Tomography images
Dhaini et al. Automated detection and measurement of corneal haze and demarcation line in spectral-domain optical coherence tomography images
CN107657605B (en) A kind of sieve plate front surface depth measurement method based on active profile and energy constraint
CN109744996B (en) OCT image BMO position positioning method
Eladawi et al. Optical coherence tomography: A review
Monemian et al. Mathematical analysis of texture indicators for the segmentation of optical coherence tomography images
EP3417401B1 (en) Method for reducing artifacts in oct using machine learning techniques
Piorkowski et al. Selected aspects of corneal endothelial segmentation quality
US20180315185A1 (en) System and method to analyze various retinal layers
Mari et al. A digital staining algorithm for optical coherence tomography images of the optic nerve head
Cheung et al. Region of interest densitometry analysis of Descemet membrane endothelial keratoplasty dehiscence on anterior segment optical coherence tomography
Ibrahim et al. Volumetric quantification of choroid and Haller's sublayer using OCT scans: An accurate and unified approach based on stratified smoothing

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