CN108618749B - 基于便携式数字化眼底照相机的视网膜血管三维重建方法 - Google Patents
基于便携式数字化眼底照相机的视网膜血管三维重建方法 Download PDFInfo
- Publication number
- CN108618749B CN108618749B CN201710174625.3A CN201710174625A CN108618749B CN 108618749 B CN108618749 B CN 108618749B CN 201710174625 A CN201710174625 A CN 201710174625A CN 108618749 B CN108618749 B CN 108618749B
- Authority
- CN
- China
- Prior art keywords
- blood vessel
- fundus
- image
- retinal
- dimensional
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 210000004204 blood vessel Anatomy 0.000 title claims abstract description 54
- 238000000034 method Methods 0.000 title claims abstract description 50
- 210000001525 retina Anatomy 0.000 title abstract description 18
- 210000001210 retinal vessel Anatomy 0.000 claims abstract description 27
- 210000001328 optic nerve Anatomy 0.000 claims abstract description 6
- 239000011159 matrix material Substances 0.000 claims description 18
- 238000004422 calculation algorithm Methods 0.000 claims description 10
- 238000003384 imaging method Methods 0.000 claims description 10
- 210000005252 bulbus oculi Anatomy 0.000 claims description 8
- 230000002207 retinal effect Effects 0.000 claims description 8
- 230000009466 transformation Effects 0.000 claims description 7
- 230000004256 retinal image Effects 0.000 claims description 4
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 238000003672 processing method Methods 0.000 claims description 3
- 238000007781 pre-processing Methods 0.000 claims description 2
- 238000012545 processing Methods 0.000 abstract description 4
- 238000004458 analytical method Methods 0.000 abstract description 3
- 238000006073 displacement reaction Methods 0.000 abstract description 2
- 239000013589 supplement Substances 0.000 abstract description 2
- 230000002792 vascular Effects 0.000 description 8
- 230000011218 segmentation Effects 0.000 description 7
- 210000001508 eye Anatomy 0.000 description 6
- 238000010586 diagram Methods 0.000 description 5
- 230000008859 change Effects 0.000 description 4
- 230000008569 process Effects 0.000 description 4
- 210000003462 vein Anatomy 0.000 description 4
- 210000001367 artery Anatomy 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 3
- 201000010099 disease Diseases 0.000 description 3
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 238000002474 experimental method Methods 0.000 description 3
- 210000004220 fundus oculi Anatomy 0.000 description 3
- 230000002123 temporal effect Effects 0.000 description 3
- 238000013519 translation Methods 0.000 description 3
- 206010025421 Macule Diseases 0.000 description 2
- 208000029091 Refraction disease Diseases 0.000 description 2
- 230000004430 ametropia Effects 0.000 description 2
- 210000003484 anatomy Anatomy 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 238000013507 mapping Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 208000014733 refractive error Diseases 0.000 description 2
- 210000001994 temporal artery Anatomy 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 101100397068 Bacillus subtilis (strain 168) iolD gene Proteins 0.000 description 1
- DFOCUWZXJBAUSQ-UHFFFAOYSA-N Berbamine Natural products O1C(C(=CC=2)O)=CC=2CC(C=23)N(C)CCC3=CC(OC)=C(OC)C=2OC(=CC=23)C(OC)=CC=2CCN(C)C3CC2=CC=C1C=C2 DFOCUWZXJBAUSQ-UHFFFAOYSA-N 0.000 description 1
- 208000010412 Glaucoma Diseases 0.000 description 1
- 206010020772 Hypertension Diseases 0.000 description 1
- 208000006550 Mydriasis Diseases 0.000 description 1
- ASCBRLGHWVZBMG-UHFFFAOYSA-N UNPD71813 Natural products COC=1C(=O)C=2C(O)=C(OC)C(OC)=CC=2OC=1C(C=C1)=CC=C1OC1OC(CO)C(O)C(O)C1O ASCBRLGHWVZBMG-UHFFFAOYSA-N 0.000 description 1
- 230000005856 abnormality Effects 0.000 description 1
- 210000000748 cardiovascular system Anatomy 0.000 description 1
- 238000012888 cubic function Methods 0.000 description 1
- 230000001186 cumulative effect Effects 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 206010012601 diabetes mellitus Diseases 0.000 description 1
- 238000002059 diagnostic imaging Methods 0.000 description 1
- 238000009792 diffusion process Methods 0.000 description 1
- 238000013399 early diagnosis Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 230000036244 malformation Effects 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000003909 pattern recognition Methods 0.000 description 1
- XDKNUCQYHGFFFE-UHFFFAOYSA-N pendulin Natural products COc1cc2ccc3c(OC)c(OC)c(OC)cc3c2c(OC)c1O XDKNUCQYHGFFFE-UHFFFAOYSA-N 0.000 description 1
- 108010034040 pendulin Proteins 0.000 description 1
- 238000004393 prognosis Methods 0.000 description 1
- 230000001105 regulatory effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 102200090266 rs751751942 Human genes 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 238000012731 temporal analysis Methods 0.000 description 1
- 230000036962 time dependent Effects 0.000 description 1
- 238000000700 time series analysis Methods 0.000 description 1
- 238000011426 transformation method Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Images
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B3/00—Apparatus for testing the eyes; Instruments for examining the eyes
- A61B3/10—Objective types, i.e. instruments for examining the eyes independent of the patients' perceptions or reactions
- A61B3/12—Objective types, i.e. instruments for examining the eyes independent of the patients' perceptions or reactions for looking at the eye fundus, e.g. ophthalmoscopes
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B3/00—Apparatus for testing the eyes; Instruments for examining the eyes
- A61B3/0016—Operational features thereof
- A61B3/0025—Operational features thereof characterised by electronic signal processing, e.g. eye models
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B3/00—Apparatus for testing the eyes; Instruments for examining the eyes
- A61B3/10—Objective types, i.e. instruments for examining the eyes independent of the patients' perceptions or reactions
- A61B3/14—Arrangements specially adapted for eye photography
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/10—Image enhancement or restoration using non-spatial domain filtering
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/90—Dynamic range modification of images or parts thereof
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
- G06T7/0012—Biomedical image inspection
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B2560/00—Constructional details of operational features of apparatus; Accessories for medical measuring apparatus
- A61B2560/04—Constructional details of apparatus
- A61B2560/0431—Portable apparatus, e.g. comprising a handle or case
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20048—Transform domain processing
- G06T2207/20064—Wavelet transform [DWT]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30041—Eye; Retina; Ophthalmic
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30101—Blood vessel; Artery; Vein; Vascular
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2210/00—Indexing scheme for image generation or computer graphics
- G06T2210/41—Medical
Landscapes
- Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- Medical Informatics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- General Health & Medical Sciences (AREA)
- Biomedical Technology (AREA)
- Biophysics (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- Heart & Thoracic Surgery (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- Molecular Biology (AREA)
- Ophthalmology & Optometry (AREA)
- Radiology & Medical Imaging (AREA)
- Quality & Reliability (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Signal Processing (AREA)
- Computer Graphics (AREA)
- Geometry (AREA)
- Software Systems (AREA)
- Eye Examination Apparatus (AREA)
Abstract
本发明属于医学图像处理及应用领域，提供了一种基于便携式数字化眼底照相机的视网膜血管三维重建方法。本发明首先使用便携式数字化眼底图像，利用照相机位移拍摄两幅眼底图像，要求两幅图像都包含视网膜视神经盘结构，并产生一定的偏移量；然后对视网膜血管进行分割并进行特征点的提取，并用本发明的基于视网膜血管解剖层次特征点匹配与三维血管重建模型对视网膜血管进行三维重建。本发明的方法能够对现有的基于二维的眼底血管观察和分析提供一种新的补充。
Description
技术领域
本发明涉及医学图像处理及应用领域，特别是涉及一种基于便携式数字化眼底照相机的视网膜血管三维重建方法。
背景技术
眼底视网膜微血管是人体的血管循环系统的一部分，由于可以通过眼底镜直接进行观察，对于临床许多影响心血管系统的疾病如糖尿病、高血压等的早期诊断与预后评判具有重要意义。传统的眼底检查模式是由眼科医师在患者散瞳后直接检查眼底进行观察，但是对眼底的改变不能进行有效的记录。近年来随着数字化技术不断发展，数字照相可以对图像进行方便的获取、保存和传输。数字眼底照相机因为可以对眼底图像进行数字化采集，广泛应用于临床和社区的眼底检查，对于眼部已经全身的疾病监测具有重要意义
但是国内外现有的眼底影像仪器一般只能得到眼底的二维图像，眼科读片人员基于该二维的眼底图像上的结构和异常进行观察和分析。而眼球本身是三维的球形结构，形成的投影二维图像难以准确的反映血管的定量改变。因此，如果能在现有的数字眼底镜基础上将眼底重要的血管结构三维重建出来，具有重要的意义。目前，国内外对眼底图像三维重建的研究尚未开展。
庞云阶等【傅晓烨，田景成，庞云阶.一种眼底图象的三维重建方法[J].中国图象图形学报1997,2(7):509-511】，根据二维图像构造中间几何图素，用传统的计算机图形学技术绘制了眼底的三维示意图，开始了基于眼底图像三维重建的研究。这种方法对逆投影重建具有一定的参考价值，但是它仅使用中间几何图素代表血管的粗细和走向，因此不能真实重建正常的眼底，也不能重现眼底的疾病改变，因此实际临床意义不大。
Corona等【Corona,E.,Mitra,S.,Wilson,M.,etal.Digital stereo imageanalyzer forgeneratingautomated 3-D measures ofopticdisc deformationinglaucoma[J].IEEE Transactions on MedicalImaging,2002,21(10):1244-1253.】研究了青光眼视盘畸形的局部三维图像重建，采用眼底相机的偏心角采集得到眼底视盘的立体观测图像，结合功率倒谱和零均值归一化互相关算法得到视盘的深度信息。
Deguchi等【Deguchi,K.,Kawamata,D.,Mizutani,K.,etal.3D fundus shapereconstructionanddisplayfrom stereo fundus images[J].IEICE TransactionsonInformation andSystems,2000,E83D(7):1408-1414.】将接触式放大镜与眼底相机的透镜组用单个透镜模型代替，计算得到透镜模型参数，采用眼底相机转动微小角度拍摄得到立体眼底图像组，三维恢复并显示了眼底图像。
李超等【李超,梁斌,陈武凡,吴德正.由二维眼底正投影图像向三维曲面逆投影成像的重建算法[J].中国生物医学工程学报,2002,21(4):346-350.】认为，如果在底片上清晰成像，那么眼底必位于焦点上，提出了中心扩散逆投影方法，先找出眼底图像上黄斑区的中心，直接定位于三维曲面的中心位置，然后按照对应关系把黄斑中心周围的点定位于曲面上相应位置。
陈骥等【陈骥,彭承琳.眼底图像的三维重建[J].生物医学工程学杂志,2008,25(01):177-181.】将眼球分为正常眼球与屈光异常眼，采用投影逆变换方法，正常眼球模拟为圆球，半径采用经验值，屈光异常眼模拟为椭球体，长短轴采用眼科超声检查仪测量得到。
Espinosa—Romero等【Espinosa-Romero,A.,Martinez-Perez,M.E.Optical 3Dreconstruction of retinalblood vessels from asequenceof views[C].8thInternationalSymposiumonLaser Metrology,Merida,Mexico，2005:605-612.】提出采用相机白校准的方法获得眼底三维坐标，并进行了改进，提出使用改进的相机.眼球的自校准方法获得眼底三维坐标。在拍摄过程中，通过旋转眼底相机获得患者眼底不同部位的图像，为了获得较大的重叠区域，要求相机视轴的旋转角度不能超过眼底相机视场角的1/4，采用2组各10幅图恢复眼底三维模型。
发明内容
本发明针对现有技术存在的问题和不足，提供一种基于便携式数字化眼底照相机的视网膜血管三维重建方法。
本发明是通过下述技术方案来解决上述技术问题的：
本发明提供一种基于便携式数字化眼底照相机的视网膜血管三维重建方法，其特点在于，其包括以下步骤：
步骤101、利用便携式数字化眼底照相机获取同一眼球的两幅偏移眼底图像并记录数字眼底镜成像的相关参数，两幅偏移眼底图像都包含视网膜视神经盘结构，并产生一定的偏移量；
步骤102、通过图像处理方法对所获得的两幅偏移眼底图像进行预处理操作；
步骤103、对预处理后获得的图像的局部熵的纹理进行提取，利用FCM算法对图像中的视网膜血管进行分割；
步骤104、利用基于距离约束的水平集演化方法对图像中眼底视神经盘进行自动分割，将分割后视网膜血管与视神经盘的二值图像进行掩模操作，获得视网膜血管的起始点，追踪出主要的血管树结构；
步骤105、对视网膜血管的起始点进行标定，追踪分割血管骨架，计算血管骨架的拓扑层次，并对血管树上的交叉点进行标定；
步骤106、利用血管解剖定义模型对血管骨架上的分叉点进行标定，对同一层次的血管特征点进行配对，获得计算三维变换矩阵的血管点集；
步骤107、通过标定的血管点匹配，并由便携式眼底照相机的成像参数计算出三维变换矩阵。
较佳地，在步骤102中，对采集到的偏移眼底图像提取绿色通道，再进行灰度转换，然后利用非下采样的离散小波变换方法对获得的视网膜图像进行增强以获得视网膜灰度图像。
在符合本领域常识的基础上，上述各优选条件，可任意组合，即得本发明各较佳实例。
本发明的积极进步效果在于：
本发明基于视网膜血管解剖特点的特征点标记，不会受图像的大小，旋转等影响。本发明的方法重建出的视网膜血管三维结构为进一步的三维坐标系下视网膜血管形态等参数测量提供了补充。
附图说明
图1为本发明方法的流程图；
图2为本发明方法的数字眼底照相机拍摄原理图；
图3为本发明方法的眼底图像的分割结果，其中A为原始眼底图像，B为A对应的分割结果；
图4为本发明方法的视神经盘分割结果，其中A为原始眼底图像；B为A视神经盘分割结果；
图5为本发明方法的血管特征点的提取结果，其中左图为分割血管树；中图为本方法获得的血管起始点；右图为本方法追踪到的血管树特征点；
图6为本发明方法的主要血管树的识别结果，其中左图为追踪得到的主要血管树；右图为主要血管树在原图上的显示结果；
图7为本发明方法的在一个血管树上的血管特征点标定结果，其中A图为分割出的血管；B图为A图血管骨架与血管特征点标定结果；
图8为本发明方法的视网膜血管特征点的匹配示意图；
图9为本发明方法的重建的3D视网膜血管图像。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明结合同一被拍摄对象的视网膜血管具有相同的解剖学结构特点，没有采用的传统的基于图像角点提取与极线约束得到匹配特征点的计算方法，而是提出了一种基于眼底血管解剖学特点的特征点匹配方法并对获得的双幅眼底图像进行了三维重建。
方法的流程图见图1，具体通过以下步骤实现：
步骤101、使用便携式数字化眼底照相机获取同一眼球的两幅偏移眼底图像并记录数字眼底镜成像的相关参数，两幅偏移眼底图像都包含视网膜视神经盘(optic disk，OD)结构，并产生一定的偏移量。
首先，利用便携式数字化眼底照相机(苏州康捷公司，型号KJ25)拍摄要用于重建的视网膜图像。数字化眼底照相机是一个CCD相机，带有一个低倍显微镜摄像头和照明装置，它用来拍摄眼睛的内表面，包括视网膜、视盘、黄斑和其他结构。尽管眼睛和眼底相机是复杂的，但是，在这里描述的分析中，他们作为一个单镜头被结合和简化。拍摄时相机和眼睛的距离通过镜桶拍摄装置保持近似恒定。成对图像的采集原理如图2，平移便携式数字眼底照相机，并拍摄眼底图像，此时眼球的相对位移和旋转变化的图像与数字眼底照相机系统的成像参数一同被记录下来进行后续操作。
步骤102、通过图像处理方法对所获得的两幅偏移眼底图像进行预处理操作，以提高视网膜血管的对比度。
为了提高视网膜血管的对比度，对采集到的RGB眼底图像(即偏移眼底图像)提取绿色通道，然后进行灰度转换，再利用非下采样的离散小波变换(Undecimated discretewavelet transform，UDWT)方法对获得的视网膜图像进行增强以获得视网膜灰度图像。
采用非下采样离散小波变换【Kui Jiang,Zhixing Zhou,Xingyun Geng,XiaofengZhang,Lemin Tang,Huiqun Wu,Jiancheng Dong.Isotropic undecimated wavelettransform fuzzy algorithm for retinal blood vessel segmentation.Journal ofMedical Imaging and Health Informatics.2015,5(7):1524-1527.】对二维图像信号进行分解。小波变换通过平移基本小波(母小波)可获得信号的时间信息，而通过缩放小波的尺度可获得信号的频率特性。对基本小波(母小波)的缩放和平移操作是为了计算小波的系数，这些系数代表小波和局部信号之间的相互关系。本实验中，设定小波母函数:ψ(t)∈L2(R)，其傅里叶变换Ψ(ω)满足如下的允许条件：
根据视网膜血管纹理自身的特点，本实验中主要用二维离散小波，将图像信号分解成沿行和列的一维问题来处理。离散小波变换针对尺度参数a、平移参数b进行离散化，离散小波变换(DWT)的定义为：
任意L2(R)空间中的x(t)DWT为：
其中，
获得信号的时间信息，而通过缩放小波的尺度可获得信号的频率特性。对基本小波(母小波)的缩放和平移操作是为了计算小波的系数，这些系数代表小波和局部信号之间的相互关系。二维离散小波可以提取眼底图像的低频系数、水平高频系数、垂直高频系数和对角高频系数，但在传统的离散小波变换在每次计算小波系数的过程中会丢弃前一步的尺度和小波结果，形成类金字塔样分解。而UWT在计算DWT过程中，使用一个N×N环形平移的矩阵TX来代替原来的X矩阵。如果X矩阵是一维的矢量X＝[X0；X1；...；XN-2；XN-1]，则
TX＝[Xn-1；X0；X1；…；Xn-2],T-1X＝[X1；X2；…；Xn-2；X0]
然后DWT计算一次X，再计算一次TX，并将各自得到的DWT系数进行融合来形成UWT的小波系数【D.B.Percival,A.T.Walden.Wavelet Methods for Time SeriesAnalysis.Cambridge University Press,Cambridge(2000)】。然后，通过对小波系数的重构得到不同尺度下的视网膜图像。
本实验中，通过UWT得到了视网膜图像不同尺度的小波系数，然后根据不同的目的对小波系数进行相应的处理，然后再对处理后的小波系数进行图像重构从而完成图像的处理。重构后的眼底图像就得到了增强。视网膜血管与视网膜图像背景的对比度明显增强。
步骤103、对该视网膜灰度图像的局部熵的纹理进行提取，利用FCM算法对视网膜血管进行分割。
选用第4层的重构图像结构中提取像素点的灰度和局部熵特征，利用FCM算法对视网膜血管进行自动分割，得到如图3所示血管视网膜的分割结果。
采用基于模糊C-均值聚类分类器的视网膜血管像素分类。模糊聚类分析就是把一个模糊的不确定的没有类别标记的样本按照某种准则划分为若干个子群集，让相似的样本尽可能归为一类，而把不相似的样本划分到不同的类中。是一种无监督模式识别。模糊聚类算法是一种基于函数最优方法的聚类算法，需要使用微积分计算技术求最优代价函数，在基于概率算法的聚类方法中使用概率密度函数。
模糊集合最根本的方法是其所包含元素及其相应的隶属度用函数表示出来，模糊集合的数学定义[13]如下:
称A是给定论域U上的一个模糊子集的含义是：对任何u∈U都对应一个唯一的实数μA(u)∈[0，1]，有序对集合{u，μA(u)|u∈U}均可作为U上的模糊集，μA(u)表示u对A的隶属度。这表明构造了一个映射：
μA(u):U→[0,1],u→μA(u)
这个映射称为模糊集A的隶属度函数。一个定义在空间X＝{x}上的隶属度函数就定义了一个模糊集合A，或者叫定义在论域X＝{x}上的模糊子集
对于有限个对象x1，x2，……，xn模糊集合～A可以表示为：
模糊C-均值聚类(FCM)的核心思想是寻求合适的隶属度和聚类中心，使得聚类内代价函数的方差和迭代误差最小，代价函数的值是像素到聚类中心L2测度的加权累积和，表达形式为：
J(U，V)表示区域的像素到聚类中心加权距离的平方和。uik表示xi属于第k类区域的隶属度。m为模糊指数，m∈[1,∞]，表示控制分类矩阵U的模糊程度，是隶属度的加权指数，决定分类结果的模糊程度，其值越大分类越模糊。同时满足约束条件:
求解此类极小值的问题，可以通过拉格朗日乘数法来求解。首先，由于矩阵u中各列独立，则有：
求解极小值问题则转化为求解
的极小值，且在约束条件
假设拉格朗日函数F为:
通过求解偏导数
可得出，在约束条件下，其隶属函数和聚类中心的更新依赖于以下两个函数：
由这两个公式得出最佳模糊分类矩阵和聚类中心，这一过程可以通过迭代方法实现。
根据需要，设定FCM的类别数为2，隶属度的加权指数m取2。最大迭代次数为100次，判别误差为1e-5，提取眼底图像中属于血管的像素点(图3)。经过测试，本血管分割方法具有较高的灵敏度(0.82)，算法运行迭代至函数收敛次数平均为50次左右，实现分割的平均时间13.4s。
步骤104、利用基于距离约束的水平集演化方法对眼底视神经盘进行自动分割。
由于要获得眼底血管的起始点位置信息，本发明需要对眼底的视神经盘结构进行分割，针对视网膜神经盘通常为圆盘状结构，而且亮度以中心向四周渐变的特点，首先构造一个各项同性的傅里叶平滑滤波器如下(1)：
通过获得眼底上的滤波器得到的系数可以获得初步的视神经盘信息，并以此得到视神经盘的初始中心。在此基础上构造一个以初始中心为中心，长度大于视神经盘直径的方形区域作为感兴趣区域。在此区域内首先构造一个时间相关的水平集模型：
然后利用基于距离约束的水平集演化方法对眼底视神经盘进行自动分割(图4)。[C.Li,C.Xu,C.Gui and M.D.Fox,Distance regularized level set evolution and itsapplication to image segmentation,IEEE Trans.Image Process.19(2010),3243–3254.]
将分割后视网膜血管与视神经盘的二值图像进行掩模操作，可以获得视网膜血管的起始点(图5)，并由此可以追踪出主要的血管树结构(图6)。
步骤105、对起始血管点进行标定，追踪分割血管骨架，计算血管骨架的拓扑层次，并对血管树上的交叉点进行标定；
如图2所示，本发明首先按照图像4邻域法则，识别出眼底血管网络上的特征点如分叉点和终末点。然后根据建立的特征点，并根据图2的层次标定方法，对眼底血管骨架网络进行重新标定，标定的结果如图3所示。本方法提出了一种基于解剖学特征的血管树的分支的编号方法，该方法首先依据视网膜的解剖学特点通常由按照部位和血管性质的不同可以分为颞侧上动脉(SRTA)、颞侧下动脉(IRTA)、颞侧上静脉(SRTV)、颞侧下静脉(IRTV)，鼻侧上动脉(NRTA)、鼻侧下动脉(NRIA)、鼻侧上静脉(NRTV)、鼻侧下静脉(NRIV)主要的血管树。以视神经盘作为视网膜树发出的起点作为一支血管树的根节点，所有的血管树上的分叉子节点都是由根节点派生出来的。节点是按由上往下，由根节点到子节点顺序编号，首先将根编号发出的分支编号为1；然后根节点向下发出的子节点的分支编号+1，以此类推，并由此对一个分割出的血管树上的特征点建立索引号，并将对应的索引号命名规则定义为SRTA.n.m,其中SRTA表示解剖学结构，n表示拓扑层次，m表示该层次下第几个特征点。比如SRTA.2.1表示为分割的颞侧上动脉的第2拓扑层次的第1个分叉点，由此完成对分割的血管树进行标定(图7)。
步骤106、利用匹配后的并用血管解剖定义模型对血管骨架上的分叉点进行标定，对同一层次的血管特征点进行配对。
将同一索引号的血管特征点进行配对，获得下一步计算三维变换矩阵的血管点集。
步骤107、通过标定的血管点匹配，并由便携式眼底照相机的成像参数计算出三维变换矩阵(图8)。
依靠从2D图像重建3D的理论，通过3×3的单数矩阵(基本矩阵A)来表现的2张单一场景的图像是由极线几何所连续的。基本矩阵A得到了2张图像上所有的几何信息，并且必须马上为重建做评估。一个标准线性相机配准矩阵K有以下项：
f是像素的焦距，δ是宽高比，(u0,v0)是主点坐标，s是倾斜因素(矩阵像素倾斜为0)。根据部分实验测试结果，为了获得基本矩阵F，最少需要7个相应的基准点。事实上，越多的点越精细。每个点的3D坐标，Xi能被计算出来。为了获得一个平滑的重建，血管分段中的相应点可以用立方函数被平滑和插入。采用一种迭代特征线性方法(Iterative-Eigen)来进行3D点的三角测量，等式可以被写成：
piT和p′iT是P和P’各自的第i行，wi，j和w′i，j是重要的因素在重复的第j步中：
wi，j＝p3TXi，j-1，w′i，j＝p′3TXi，j-1
在wi，0和w′i，0开始重复时被建为1为了找到Xi，0最初的解法。
由此对血管树上的所有点进行三维变换，得到了三维重建后的视网膜血管(图9)。
虽然以上描述了本发明的具体实施方式，但是本领域的技术人员应当理解，这些仅是举例说明，本发明的保护范围是由所附权利要求书限定的。本领域的技术人员在不背离本发明的原理和实质的前提下，可以对这些实施方式做出多种变更或修改，但这些变更和修改均落入本发明的保护范围。
Claims (2)
1.一种基于便携式数字化眼底照相机的视网膜血管三维重建方法，其特征在于，其包括以下步骤：
步骤101、利用便携式数字化眼底照相机获取同一眼球的两幅偏移眼底图像并记录便携式数字化眼底照相机成像的相关参数，两幅偏移眼底图像都包含视网膜视神经盘结构，并产生一定的偏移量；
步骤102、通过图像处理方法对所获得的两幅偏移眼底图像进行预处理操作；
步骤103、对预处理后获得的图像的局部熵的纹理进行提取，利用FCM算法对图像中的视网膜血管进行分割；
步骤104、利用基于距离约束的水平集演化方法对图像中眼底视神经盘进行自动分割，将分割后视网膜血管与视神经盘的二值图像进行掩模操作，获得视网膜血管的起始点，追踪出主要的血管骨架结构；
步骤105、对视网膜血管的起始点进行标定，追踪分割血管骨架，计算血管骨架的拓扑层次，并对血管骨架上的交叉点进行标定；
步骤106、利用血管解剖定义模型对血管骨架上的交叉点进行标定，对同一层次的血管特征点进行配对，获得计算三维变换矩阵的血管点集；
步骤107、通过标定的血管点匹配，并由便携式眼底照相机的成像参数计算出三维变换矩阵。
2.如权利要求1所述的视网膜血管三维重建方法，其特征在于，在步骤102中，对采集到的偏移眼底图像提取绿色通道，再进行灰度转换，然后利用非下采样的离散小波变换方法对获得的视网膜图像进行增强以获得视网膜灰度图像。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710174625.3A CN108618749B (zh) | 2017-03-22 | 2017-03-22 | 基于便携式数字化眼底照相机的视网膜血管三维重建方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710174625.3A CN108618749B (zh) | 2017-03-22 | 2017-03-22 | 基于便携式数字化眼底照相机的视网膜血管三维重建方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108618749A CN108618749A (zh) | 2018-10-09 |
| CN108618749B true CN108618749B (zh) | 2020-06-19 |
Family
ID=63707103
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710174625.3A Active CN108618749B (zh) | 2017-03-22 | 2017-03-22 | 基于便携式数字化眼底照相机的视网膜血管三维重建方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108618749B (zh) |
Families Citing this family (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109816631A (zh) * | 2018-12-25 | 2019-05-28 | 河海大学 | 一种基于新代价函数的图像分割方法 |
| CN109741349B (zh) * | 2019-01-24 | 2021-12-07 | 江门市中心医院 | 一种缺血性脑卒中图像分割的方法 |
| CN110163119A (zh) * | 2019-04-30 | 2019-08-23 | 中国地质大学（武汉） | 一种手指静脉识别方法及系统 |
| CN110414423A (zh) * | 2019-07-25 | 2019-11-05 | 上海鹰瞳医疗科技有限公司 | 身份识别方法及设备 |
| CN110399833B (zh) * | 2019-07-25 | 2023-03-24 | 上海鹰瞳医疗科技有限公司 | 身份识别方法、建模方法及设备 |
| CN110448267B (zh) * | 2019-09-06 | 2021-05-25 | 重庆贝奥新视野医疗设备有限公司 | 一种多模眼底动态成像分析系统及其方法 |
| CN111145240A (zh) * | 2019-11-18 | 2020-05-12 | 西宁市动物疫病预防控制中心(挂西宁市畜牧兽医站牌子) | 一种基于3d相机的活体西门塔尔牛体尺在线测量方法 |
| CN113658243B (zh) * | 2021-10-18 | 2022-02-22 | 广东唯仁医疗科技有限公司 | 眼底三维模型建立方法、眼底照相机、装置和存储介质 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101301207A (zh) * | 2008-05-28 | 2008-11-12 | 华中科技大学 | 动态模型指导下的血管造影三维重建方法 |
| CN103679801A (zh) * | 2013-12-01 | 2014-03-26 | 北京航空航天大学 | 一种基于多视角x光片的心血管三维重建方法 |
| CN103914874A (zh) * | 2014-04-08 | 2014-07-09 | 中山大学 | 一种无特征提取的紧致sfm三维重建方法 |
| CN104881862A (zh) * | 2015-04-03 | 2015-09-02 | 南通大学 | 一种基于眼底镜图像的视网膜血管迂曲度计算方法及其应用 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6047080A (en) * | 1996-06-19 | 2000-04-04 | Arch Development Corporation | Method and apparatus for three-dimensional reconstruction of coronary vessels from angiographic images |
| US20070165917A1 (en) * | 2005-11-26 | 2007-07-19 | Zhujiang Cao | Fully automatic vessel tree segmentation |
-
2017
- 2017-03-22 CN CN201710174625.3A patent/CN108618749B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101301207A (zh) * | 2008-05-28 | 2008-11-12 | 华中科技大学 | 动态模型指导下的血管造影三维重建方法 |
| CN103679801A (zh) * | 2013-12-01 | 2014-03-26 | 北京航空航天大学 | 一种基于多视角x光片的心血管三维重建方法 |
| CN103914874A (zh) * | 2014-04-08 | 2014-07-09 | 中山大学 | 一种无特征提取的紧致sfm三维重建方法 |
| CN104881862A (zh) * | 2015-04-03 | 2015-09-02 | 南通大学 | 一种基于眼底镜图像的视网膜血管迂曲度计算方法及其应用 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108618749A (zh) | 2018-10-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108618749B (zh) | 基于便携式数字化眼底照相机的视网膜血管三维重建方法 | |
| Abràmoff et al. | Retinal imaging and image analysis | |
| Giancardo et al. | Textureless macula swelling detection with multiple retinal fundus images | |
| CN113450294A (zh) | 多模态医学图像配准融合方法、装置及电子设备 | |
| CN113506334A (zh) | 基于深度学习的多模态医学图像融合方法及系统 | |
| JP7413147B2 (ja) | 画像処理装置、画像処理方法、及びプログラム | |
| CN108416793B (zh) | 基于三维相干断层成像图像的脉络膜血管分割方法及系统 | |
| US11490830B2 (en) | Apparatus and method for quantification of the mapping of the sensory areas of the brain | |
| JP2022549669A (ja) | 時空間データに基づき医用画像を解析するシステム及び方法 | |
| Duan et al. | A generative model for OCT retinal layer segmentation by groupwise curve alignment | |
| JP7332463B2 (ja) | 制御装置、光干渉断層撮影装置、光干渉断層撮影装置の制御方法、及びプログラム | |
| CN115393239A (zh) | 一种多模态眼底图像配准和融合方法及系统 | |
| CN111445575A (zh) | 威利斯环的图像重建方法、装置、电子设备、存储介质 | |
| Hao et al. | Magnetic resonance image segmentation based on multi-scale convolutional neural network | |
| CN112562058B (zh) | 一种基于迁移学习的颅内血管模拟三维模型快速建立方法 | |
| CN116407080A (zh) | 近视患者眼底结构的演变识别及3d可视化系统和方法 | |
| CN112669256B (zh) | 一种基于迁移学习的医学图像分割与显示方法 | |
| CN115147273A (zh) | 一种基于域适应的眼科oct术后轮廓图生成方法及系统 | |
| Al-Fahdawi et al. | An automatic corneal subbasal nerve registration system using FFT and phase correlation techniques for an accurate DPN diagnosis | |
| Yu et al. | 3D reconstruction of fundus vascular structure based on binocular stereovision | |
| Wang et al. | Cross-Domain Multi-disease Ocular Disease Recognition via Data Enhancement | |
| Wang et al. | Two-stage Cross-domain Ocular Disease Recognition with Data Augmentation | |
| Xu et al. | Three dimensional optic disc visualisation from stereo images via dual registration and ocular media optical correction | |
| CN112669439B (zh) | 基于迁移学习的颅内血管造影增强三维模型的建立方法 | |
| Gaddipati | Towards developing a multiple modality fusion technique for automatic detection of Glaucoma |
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 |