CN112330638B - 一种视网膜oct图像水平配准和图像增强方法 - Google Patents

一种视网膜oct图像水平配准和图像增强方法 Download PDF

Info

Publication number
CN112330638B
CN112330638B CN202011237202.XA CN202011237202A CN112330638B CN 112330638 B CN112330638 B CN 112330638B CN 202011237202 A CN202011237202 A CN 202011237202A CN 112330638 B CN112330638 B CN 112330638B
Authority
CN
China
Prior art keywords
image
registered
registration
correlation coefficient
pixels
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
Application number
CN202011237202.XA
Other languages
English (en)
Other versions
CN112330638A (zh
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.)
Suzhou University
Original Assignee
Suzhou 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 Suzhou University filed Critical Suzhou University
Priority to CN202011237202.XA priority Critical patent/CN112330638B/zh
Publication of CN112330638A publication Critical patent/CN112330638A/zh
Application granted granted Critical
Publication of CN112330638B publication Critical patent/CN112330638B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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
    • G06T5/00Image enhancement or restoration
    • G06T5/70Denoising; Smoothing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • 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/30Determination of transform parameters for the alignment of images, i.e. image registration
    • 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/20Special algorithmic details
    • G06T2207/20024Filtering details
    • G06T2207/20032Median filtering
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20092Interactive image processing based on input by user
    • G06T2207/20104Interactive definition of region of interest [ROI]
    • 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)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Medical Informatics (AREA)
  • Quality & Reliability (AREA)
  • Radiology & Medical Imaging (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Image Processing (AREA)
  • Eye Examination Apparatus (AREA)

Abstract

本申请实施例公开一种视网膜OCT图像水平配准和图像增强方法,该方法包括以下步骤:S1、图像配准:将待配准图与基准图一一进行配准,配准方向包括垂直方向与水平方向;S2、图像处理与结果获取:将配准后的图像和基准图进行叠加平均,对平均图分出外核层ONL,对图像中外核层ONL以上的像素进行加权自适应伽马校正AGCWD,以进一步增强对比度,得到最终的结果。本申请的视网膜OCT图像水平配准和图像增强方法,在使用任意一张待配准图与基准图配准时,优化配准流程,提升水平方向上的配准精度与速度。同时,还对平均后图像中外核层ONL以上部分进行局部对比度增强,进一步提升图像质量,使得层次结构与细节特征更加鲜明与突出。

Description

一种视网膜OCT图像水平配准和图像增强方法
技术领域
本申请涉及视网膜图像处理方法技术领域,具体是一种视网膜OCT图像水平配准和图像增强方法。
背景技术
光学相干断层扫描成像(Optical Coherence Tomography,OCT)因具有高分辨率、无损伤等特点,在眼科领域特别是眼底视网膜和脉络膜的组织成像观测中广泛应用。散斑噪声由光源热噪声和生物组织随机散射等引起,是OCT***的固有噪声,难以使用基于硬件的设计方法去除,从而降低了OCT图像的对比度并且可能会遮挡视网膜的细节特征,影响对相关组织结构异常的辨识和判断。OCT设备可以对同一位置进行扫描多次以获取多幅眼底图像,通过配准平均以得到高清图像。但因人眼在拍摄过程中不自觉的微小位移而导致OCT图像出现几十个像素以上的偏移,不配准而直接平均将导致得到的图像模糊,视网膜分层结构与细节特征变得难以分辨,因此需要精确的配准以得到清晰的图像。
现有技术[1]-《“Enhancing the signal-to-noise ratio in ophthalmicoptical coherence tomography by image registration—method and clinicalexamples”,Journal of Biomedical Optics,12(4),041208(2007)》中公开了一种方法是在轴向上运用正则动态规划技术校正轴向偏移,在水平方向上通过计算一系列位置的互相关度来获取最佳配准位置。该方法首次提出了基于水平与垂直方向的配准技术并有较佳的临床实例,但在水平方向上需要一个个位置进行计算,难以在短时间内进行较大范围的横向配准,计算效率较低。
现有技术[2]-申请号为CN106504228A的中国专利《“一种眼科OCT图像的大范围高清快速配准方法”》尝试解决这一问题,该方法在横向上使用二分法来进行水平运动伪差校正,通过比较配准范围两端的图像互相关度值以确定配准方向并缩小范围,逐次迭代直至找到最佳配准位置。这样做极大地减少了运算量,提高了配准算法的运算效率。但实际情况中越偏离图像最佳配准位置处计算相似度越容易产生误差,可能导致第一次迭代后配准方向选择错误,从而进行了无效配准。另外,通过配准平均得到的图像中,光反射率较大的脉络膜至外核层ONL(outer nuclear layer,外核层)亮度较大,受平均的影响较小。但在反射率较低的结构中(即视网膜外核层ONL以上部分,称为内视网膜区域),经过平均,其亮度值会轻微降低,导致局部对比度较低,层次结构分辨不清。
综上,现有技术中存在下列问题:传统方法需要进行较为复杂的互相关度计算,计算时间成本高;传统方法需要一个个位置进行计算,难以在短时间内进行较大范围的横向配准,另外所提出的二分法会出现首次迭代时配准方向选择错误的情况,从而导致了无效配准,配准准确性受到影响;传统方法在平均后视网膜层次较为模糊,特别是反射率较低的内视网膜区域经过平均后,其亮度值会轻微降低,导致局部对比度较低。
发明内容
本申请旨在解决上述技术问题,提供一种视网膜OCT图像水平配准和图像增强方法,在使用任意一张待配准图与基准图配准时,优化配准流程,提升水平方向上的配准精度与速度。同时,还对平均后图像中外核层ONL以上部分进行局部对比度增强,进一步提升图像质量,使得层次结构与细节特征更加鲜明与突出。
为实现上述目的,本申请公开了一种视网膜OCT图像水平配准和图像增强方法,该方法包括以下步骤:
S1、图像配准:将待配准图与基准图一一进行配准,配准方向包括垂直方向与水平方向;
S2、图像处理与结果获取:将配准后的图像和基准图进行叠加平均,对平均图分出外核层ONL,对图像中外核层ONL以上的像素进行加权自适应伽马校正AGCWD,以进一步增强对比度,得到最终的结果。
作为优选,所述S1图像配准具体包括:
S1-1、将基准图与第i幅待配准图每隔s_col列进行采样,认定采样列的轴向偏移与采样列所在s_col列的轴向偏移是一致的;
S1-2、将第i幅采样后的待配准图向左偏移R个像素,R为所设定的待配准图最大的水平偏移距离;将当前爬山位置记为p_cur=-R,进行轴向配准,根据评价函数计算与采样后的基准图的相关系数并给予保留;然后进行步骤S1-3;
S1-3、将当前爬山位置处的图像继续向右偏移step个像素,再进行轴向配准,根据评价函数得到对应位置的相关系数并给予保留,当前位置变成p_cur=p_cur+step,其中step为所设定的初始爬山步长;重复操作直至当前位置的相关系数小于前一位置的相关系数;当当前位置的相关系数小于前一位置的相关系数时,改变移动方向,将步长缩短至一半,即step’=step/2,进入步骤S1-4;
S1-4、将当前爬山位置处的图像向左偏移step’个像素,当前位置更新为p_cur=p_cur-step;判断是否存在当前位置的相关系数,若不存在,继续进行轴向配准,根据评价函数计算得到该位置的相关系数并给予保留;重复操作直至当前位置的相关系数小于前一位置的相关系数;当当前位置的相关系数小于前一位置的相关系数时,改变移动方向,将步长缩短至一半,step”=step/2,进入步骤S1-5;
S1-5、重复步骤S1-3与步骤S1-4步骤,直至step”<1;若当step”<1,则前一爬山位置即为此幅图像的最佳水平偏移位置,经前一爬山位置偏移再做轴向配准后,得到的图像即为第i幅配准后图像;
S1-6、参照步骤S1-1~步骤S1-5,将剩余待配准图像与基准图一一进行配准,得到多幅配准后的图像。
作为优选,所述评价函数的计算公式为:
Figure BDA0002767081780000031
其中,I、J为尺寸大小是M*N(M为图像行数,N为图像列数)的两幅图像,I(i,j)与J(i,j)为各自图像中的像素,
Figure BDA0002767081780000032
与/>
Figure BDA0002767081780000033
分别为各自图像像素的平均值。
作为优选,所述轴向配准具体包括:
T1、在采样后的基准图中分割出视网膜色素上皮层RPE:对合成图进行中值滤波以消除噪声的干扰,将图像像素强度作为代价函数,运用正则动态规划算法,选取出一条灰度值最大的路径,使用savitzky-golay滤波器对路径进行平滑,平滑后的曲线即为视网膜色素上皮层RPE;
T2、确定视网膜色素上皮层RPE,在合成图中分割出内界膜层ILM:使用sobel算子对滤波后的合成图进行卷积,以检测图像垂直方向上的边缘,选取正边缘梯度作为评价函数,将内界膜层ILM的搜索范围限定在视网膜色素上皮层RPE上方30像素以上;运用正则动态规划算法,选取出正边缘梯度最大的一条路径,使用savitzky-golay滤波器对路径进行平滑,平滑后的曲线即为内界膜层ILM;
T3、在采样后的第i幅待配准图中按照步骤T1、步骤T2分割出视网膜色素上皮层RPE与内界膜层ILM;采样后基准图内界膜层ILM与采样后第i幅待配准图内界膜层ILM各列的差值即为待配准图各列的轴向偏移距离,原第i幅待配准图的对应s_col列的轴向偏移距离与采样列轴向偏移距离一致,将原第i幅待配准图的各列进行偏移校正,从而完成轴向配准。
作为优选,所述S2图像处理与结果获取具体包括:
S2-1、将配准后的图像和基准图进行叠加平均,得到合成图;
S2-2、在合成图中,分割出外核层ONL,对外核层ONL以上部分进行加权自适应伽马校正AGCWD。
作为优选,所述S2-2中外核层ONL的分割步骤具体包括:
在合成图中分割出视网膜色素上皮层RPE;分割出视网膜色素上皮层RPE后,使用sobel算子对滤波后的合成图进行卷积检测图像垂直方向上的边缘;将边缘梯度大于0的值置为0,只选取负边缘梯度的绝对值作为评价函数,将外核层ONL的搜索范围限定与视网膜色素上皮层RPE上方5像素至30像素之间,运用正则动态规划算法找到一条使得评价函数最大的路径,该路径即为粗略的ONL曲线,并对该曲线进行savitzky-golay滤波,得到最终的外核层ONL。
作为优选,所述视网膜色素上皮层RPE的步骤具体包括:
对合成图进行中值滤波以消除噪声的干扰,将图像像素强度作为代价函数,运用正则动态规划算法,选取出一条灰度值最大的路径,使用savitzky-golay滤波器对路径进行平滑,平滑后的曲线即为视网膜色素上皮层RPE。
作为优选,所述S2-2具体包括:
S2-2-1、计算合成图的灰度值直方图pdf,然后计算每个灰度值的加权直方图参数pdfw(l),通过加权直方图参数可计算出各灰度值的加权分布函数cdfw(l),进入步骤S2-2-2;
S2-2-2、找到合成图中最大的灰度值lmax,选取合成图中外核层ONL以上的部分作为感兴趣区域ROI,使用加权分布函数cdfw(l)作为自适应参数,对感兴趣区域ROI中各灰度值的像素进行自适应伽马校正AGCWD,得到经局部增强后的结果图。
作为优选,所述的加权直方图参数pdfw(l)的计算公式为:
Figure BDA0002767081780000041
所述加权分布函数cdfw(l)的计算公式为:
Figure BDA0002767081780000042
式中l为图像的灰度值,pdfmax与pdfmin分别为直方图统计结果的最大值与最小值,α为调整参数。
1作为优选,所述自适应伽马校正AGCWD的计算公式为:
Figure BDA0002767081780000043
式中T(l)即为每个灰度值经对比度增强后的结果。
有益效果:本发明能在较短时间内实现高精度的视网膜OCT图像配准,在多幅平均后,保证了视网膜OCT图像边缘的锐利性。同时克服了多幅平均后合成图像中视网膜局部对比度不高的问题,将视网膜的细节特征与分层结构更清晰的展现出来,显著提高了视网膜OCT图像质量,方便对细小组织结构异常的辨识。
附图说明
为了更清楚地说明本申请实施例或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅是本申请的一些实施例,对于本领域技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。
图1是所设定水平配准范围中每一位置偏移后的图像与基准图计算出的互相关度折线图。
图2是采样后的基准图。
图3为配准后的图像及其局部放大图,其中(a)为仅垂直配准后的图像、(b)为经过水平垂直配准后的图像。
图4是对平均后图像分层后的得到的结果。
图5是区域对比度提升前后的结果图及局部放大图。
具体实施方式
为了使本技术领域的人员更好地理解本申请中的技术方案,下面将结合本申请实施例中的附图,对本申请实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅是本申请一部分实施例,而不是全部的实施例。基于本申请中的实施例,本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例,都应当属于本申请保护的范围。
实施例:一种视网膜OCT图像水平配准和图像增强方法,该方法包括以下步骤:
S1、图像配准:将待配准图与基准图一一进行配准,配准方向包括垂直方向与水平方向;
S2、图像处理与结果获取:将配准后的图像和基准图进行叠加平均,对平均图分出外核层ONL,对图像中外核层ONL以上的像素进行加权自适应伽马校正AGCWD,以进一步增强对比度,得到最终的结果。
具体地,步骤S1展开为以下步骤:
S1-1、首先从待配准图像中拿出一幅图进行下采样,每隔s_col(此例s_col=2)列进行采样,认定采样列的轴向偏移与采样列所在s_col列的轴向偏移是一致的。同时按上述方法对基准图进行下采样。
S1-2、将采样后的待配准图向左偏移32个位置,即所设定的最大水平偏移距离R=32。然后进行轴向配准,根据评价函数计算出互相关度并保存在数组中,当前位置p_cur=-R;
S1-3、设定初始爬山步长step=16,为避免陷入局部最大值的情况,初始步长应设定为较大值。开始爬山,将原位置的图像向右移动step个像素,当前位置随之更新,p_cur=p_cur+step。然后进行轴向配准,并根据评价函数算出相关系数,保存在数组中。继续爬山,重复上述操作,直至爬过了山顶,即当前位置的相关系数小于前一位置的相关系数,改变移动方向,缩短步长step’=step/2,进入下一动作。
S1-4、将当前爬山位置处的图像向左偏移step’个像素,当前位置更新为p_cur=p_cur-step。判断先前是否经过此位置,是否存在当前位置的相关系数,若不存在,继续进行轴向配准,根据评价函数计算得到该位置的相关系数并给予保留,若存在则无需进行重复计算,可减少运算时间。重复这一操作直至当前位置的相关系数小于前一位置的相关系数,即再次爬过了山顶,改变移动方向,将步长缩短至一半,step”=step’/2,进入下一步动作(即步骤S1-5)。
S1-5、重复S1-3与S1-4步骤,直至step”<1。若当step”<1,则前一爬山位置即为此幅图像的最佳水平偏移位置,经前一爬山位置偏移再做轴向配准后,得到的图像即为第i幅配准后图像。进入下一动作(即步骤S1-6)。
评价函数的计算公式为:
Figure BDA0002767081780000061
其中,I、J为尺寸大小是M*N(M为图像行数,N为图像列数)的两幅图像,I(i,j)与J(i,j)为各自图像中的像素,/>
Figure BDA0002767081780000062
与/>
Figure BDA0002767081780000063
分别为各自图像像素的平均值。
在本实施例中,正则动态规划方法参见现有技术[1]。
进一步地,所述的轴向配准方法见下:
T1、在采样后的基准图中分割出视网膜色素上皮层RPE(retinal pigmentedepithelium,视网膜色素上皮层):对合成图进行中值滤波(此例选择尺寸为11的卷积核)以消除噪声的干扰,将图像像素强度作为代价函数,运用正则动态规划算法,选取出一条灰度值最大的路径,使用savitzky-golay滤波器(设定滑动窗口尺寸为5,下同)对路径进行平滑,平滑后的曲线即为视网膜色素上皮层RPE;
T2、确定视网膜色素上皮层RPE,在合成图中分割出内界膜层ILM(internallimiting membrane,内界膜):使用sobel算子(尺寸为1)对滤波后的合成图进行卷积,以检测图像垂直方向上的边缘,选取正边缘梯度作为评价函数,将内界膜层ILM的搜索范围限定在视网膜色素上皮层RPE上方30像素以上;运用正则动态规划算法,选取出正边缘梯度最大的一条路径,使用savitzky-golay滤波器对路径进行平滑,平滑后的曲线即为内界膜层ILM;
T3、在采样后的第i幅待配准图中按照步骤T1、步骤T2分割出视网膜色素上皮层RPE与内界膜层ILM;采样后基准图内界膜层ILM与采样后第i幅待配准图内界膜层ILM各列的差值即为待配准图各列的轴向偏移距离,原第i幅待配准图的对应s_col列(此例为2)的轴向偏移距离与采样列轴向偏移距离一致,将原第i幅待配准图的各列进行偏移校正,从而完成轴向配准。本专利采用分层方法以进行轴向配准,不需要进行复杂的互相关计算,提升了算法的运行速度。
请参阅图1,我们在一组眼底OCT数据中选取了一幅具有水平偏移的图像,所设定的水平配准范围是-32~32像素,我们算出经每一位置偏移后的图像与基准图的互相关系数,绘制成了一幅折线图。从图中可以看到,数据大体上呈现了一条抛物线的形状,所述的最佳水平偏移位置即为所求曲线的峰值位置。另外,越远离最佳偏移位置的互相关系数越容易出现误差(理想的形状是具有单峰性的一条抛物线),这就导致现有技术[2]中所提到的二分法会对此类图像造成无效配准。该方法首先计算了-32,0,32位置处的相关系数,在此幅图像中-32位置处的相关系数大于32位置处的相关系数,因此其会在-32位置与0位置之间继续进行二分,从而确定最佳偏移位置,而正确的偏移位置在0位置与32位置之间。本发明的爬山法则可以避免这一问题,使配准更为精确。
参阅图2,图2是该组数据经下采样(每隔两列采样一列)后的基准图,以采样后的基准图内界膜层ILM为基准,计算采样后的待配准图内界膜层ILM与基准的偏移距离,认定原待配准图相邻列与采样列的偏移距离一致,从而完成轴向偏移校正。
参阅图3,(a)图为图1所选取的图像仅经过轴向配准所得到的配准图,(b)图为该幅图像经爬山法进行水平配准与轴向配准所得到的图像。可见本发明的爬山法可精确校准待配准图的水平偏移。框图部分为局部放大图,水平配准后的黄斑边缘相较于未水平配准的黄斑边缘更为锐利。
S1-6、重复S11~S15步骤,将待配准图与基准图一一进行配准,得到若干张配准后的图像,进入下一步骤S2。
具体地,S2步骤具体包括以下步骤:
S2-1、将若干张配准后的图像与基准图进行叠加平均,得到合成图,进入步骤S2-2;
S2-2、在合成图中分割出外核层ONL,在外核层ONL以上部分进行局部对比度增强。
进一步地,S2-2中外核层ONL的分割步骤具体包括:
首先分割出视网膜色素上皮层RPE以确定外核层ONL的搜索区域。对合成图进行中值滤波平滑边缘并滤除部分噪声,选取图像强度作为代价函数,运用正则动态规划算法找到一条灰度值最大的路径,经savitzky-golay滤波器平滑后得到视网膜色素上皮层RPE,使用sobel算子(尺寸为1)对滤波后的合成图进行卷积,以检测图像垂直方向上的边缘。将边缘梯度大于0的值置为0,只选取负边缘梯度的绝对值作为评价函数,将外核层ONL的搜索范围限定与视网膜色素上皮层RPE上方5像素至30像素之间,运用正则动态规划算法找到负梯度最大的路径,并对该路径进行savitzky-golay滤波,得到最终的外核层ONL。
请参阅图4,图3即为该组眼底OCT数据进行配准叠加后得到的合成图,图中分层符合视网膜特征结构。
进一步地,S2-2具体包括:
S2-2-1、计算合成图的灰度值直方图pdf,然后计算每个灰度值的加权直方图参数pdfw(l),通过加权直方图参数可计算出各灰度值的加权分布函数cdfw(l),进入步骤S2-2-2;
S2-2-2、找到合成图中最大的灰度值lmax,选取合成图中外核层ONL以上的部分作为感兴趣区域ROI,使用加权分布函数cdfw(l)作为自适应参数,对感兴趣区域ROI中各灰度值的像素进行自适应伽马校正AGCWD,得到经局部增强后的结果图。
加权直方图参数pdfw(l)的计算公式为:
Figure BDA0002767081780000081
加权分布函数cdfw(l)的计算公式为:
Figure BDA0002767081780000082
式中l为图像的灰度值,pdfmax与pdfmin分别为直方图统计结果的最大值与最小值,α为调整参数。
计算出合成图的灰度直方图,并通过直方图计算各灰度值的加权直方图,其目的在于轻微修改统计直方图使数据分布更为均匀,其中α设为0.01。通过加权直方图参数计算加权分布函数,加权分布函数作为一个自适应补偿值,能够逐渐增强较低亮度的像素并避免过度增强较高亮度的像素。
自适应伽马校正AGCWD的计算公式为:
Figure BDA0002767081780000083
式中T(l)即为每个灰度值经对比度增强后的结果。
请参阅图5,图(a)为未进行局部对比度增强的结果图,图(b)为进行局部对比度增强后的结果图,框图内为局部放大图。可见图像中层次结构更为清晰,信号强度较弱的外界膜ELM(external limiting membrane,外界膜)也能较为明显地分辨。
以上描述是为了进行图示说明而不是为了进行限制。通过阅读上述描述,在所提供的示例之外的许多实施方式和许多应用对本领域技术人员来说都将是显而易见的。因此,本教导的范围不应该参照上述描述来确定,而是应该参照所附权利要求以及这些权利要求所拥有的等价物的全部范围来确定。出于全面之目的,所有文章和参考包括专利申请和公告的公开都通过参考结合在本文中。在前述权利要求中省略这里公开的主题的任何方面并不是为了放弃该主体内容,也不应该认为申请人没有将该主题考虑为所公开的申请主题的一部分。

Claims (9)

1.一种视网膜OCT图像水平配准和图像增强方法,其特征在于,该方法包括以下步骤:
S1、图像配准:将待配准图与基准图一一进行配准,配准方向包括垂直方向与水平方向;其具体包括:
S1-1、将基准图与第i幅待配准图每隔s_col列进行采样,认定采样列的轴向偏移与采样列所在s_col列的轴向偏移是一致的;
S1-2、将第i幅采样后的待配准图向左偏移R个像素,R为所设定的待配准图最大的水平偏移距离;将当前爬山位置记为p_cur=-R,进行轴向配准,根据评价函数计算与采样后的基准图的相关系数并给予保留;然后进行步骤S1-3;
S1-3、将当前爬山位置处的图像继续向右偏移step个像素,再进行轴向配准,根据评价函数得到对应位置的相关系数并给予保留,当前位置变成p_cur=p_cur+step,其中step为所设定的初始爬山步长;重复操作直至当前位置的相关系数小于前一位置的相关系数;当当前位置的相关系数小于前一位置的相关系数时,改变移动方向,将步长缩短至一半,即step’=step/2,进入步骤S1-4;
S1-4、将当前爬山位置处的图像向左偏移step’个像素,当前位置更新为p_cur=p_cur-step;判断是否存在当前位置的相关系数,若不存在,继续进行轴向配准,根据评价函数计算得到该位置的相关系数并给予保留;重复操作直至当前位置的相关系数小于前一位置的相关系数;当当前位置的相关系数小于前一位置的相关系数时,改变移动方向,将步长缩短至一半,step”=step/2,进入步骤S1-5;
S1-5、重复步骤S1-3与步骤S1-4步骤,直至step”<1;若当step”<1,则前一爬山位置即为此幅图像的最佳水平偏移位置,经前一爬山位置偏移再做轴向配准后,得到的图像即为第i幅配准后图像;
S1-6、参照步骤S1-1~步骤S1-5,将剩余待配准图像与基准图一一进行配准,得到多幅配准后的图像;
S2、图像处理与结果获取:将配准后的图像和基准图进行叠加平均,对平均图分出外核层ONL,对图像中外核层ONL以上的像素进行加权自适应伽马校正AGCWD,以进一步增强对比度,得到最终的结果。
2.根据权利要求1所述的视网膜OCT图像水平配准和图像增强方法,其特征在于,所述评价函数的计算公式为:
Figure FDA0004227484180000011
其中,I、J为尺寸大小是M*N(M为图像行数,N为图像列数)的两幅图像,I(i,j)与J(i,j)为各自图像中的像素,
Figure FDA0004227484180000012
与/>
Figure FDA0004227484180000013
分别为各自图像像素的平均值。
3.根据权利要求1所述的视网膜OCT图像水平配准和图像增强方法,其特征在于,所述轴向配准具体包括:
T1、在采样后的基准图中分割出视网膜色素上皮层RPE:对合成图进行中值滤波以消除噪声的干扰,将图像像素强度作为代价函数,运用正则动态规划算法,选取出一条灰度值最大的路径,使用savitzky-golay滤波器对路径进行平滑,平滑后的曲线即为视网膜色素上皮层RPE;
T2、确定视网膜色素上皮层RPE,在合成图中分割出内界膜层ILM:使用sobel算子对滤波后的合成图进行卷积,以检测图像垂直方向上的边缘,选取正边缘梯度作为评价函数,将内界膜层ILM的搜索范围限定在视网膜色素上皮层RPE上方30像素以上;运用正则动态规划算法,选取出正边缘梯度最大的一条路径,使用savitzky-golay滤波器对路径进行平滑,平滑后的曲线即为内界膜层ILM;
T3、在采样后的第i幅待配准图中按照步骤T1、步骤T2分割出视网膜色素上皮层RPE与内界膜层ILM;采样后基准图内界膜层ILM与采样后第i幅待配准图内界膜层ILM各列的差值即为待配准图各列的轴向偏移距离,原第i幅待配准图的对应s_col列的轴向偏移距离与采样列轴向偏移距离一致,将原第i幅待配准图的各列进行偏移校正,从而完成轴向配准。
4.根据权利要求1所述的视网膜OCT图像水平配准和图像增强方法,其特征在于,所述S2图像处理与结果获取具体包括:
S2-1、将配准后的图像和基准图进行叠加平均,得到合成图;
S2-2、在合成图中,分割出外核层ONL,对外核层ONL以上部分进行加权自适应伽马校正AGCWD。
5.根据权利要求4所述的视网膜OCT图像水平配准和图像增强方法,其特征在于,所述S2-2中外核层ONL的分割步骤具体包括:
在合成图中分割出视网膜色素上皮层RPE;分割出视网膜色素上皮层RPE后,使用sobel算子对滤波后的合成图进行卷积检测图像垂直方向上的边缘;将边缘梯度大于0的值置为0,只选取负边缘梯度的绝对值作为评价函数,将外核层ONL的搜索范围限定与视网膜色素上皮层RPE上方5像素至30像素之间,运用正则动态规划算法找到一条使得评价函数最大的路径,该路径即为粗略的ONL曲线,并对该曲线进行savitzky-golay滤波,得到最终的外核层ONL。
6.根据权利要求5所述的视网膜OCT图像水平配准和图像增强方法,其特征在于,所述视网膜色素上皮层RPE的步骤具体包括:
对合成图进行中值滤波以消除噪声的干扰,将图像像素强度作为代价函数,运用正则动态规划算法,选取出一条灰度值最大的路径,使用savitzky-golay滤波器对路径进行平滑,平滑后的曲线即为视网膜色素上皮层RPE。
7.根据权利要求4所述的视网膜OCT图像水平配准和图像增强方法,其特征在于,所述S2-2具体包括:
S2-2-1、计算合成图的灰度值直方图pdf,然后计算每个灰度值的加权直方图参数pdfw(l),通过加权直方图参数可计算出各灰度值的加权分布函数cdfw(l),进入步骤S2-2-2;
S2-2-2、找到合成图中最大的灰度值lmax,选取合成图中外核层ONL以上的部分作为感兴趣区域ROI,使用加权分布函数cdfw(l)作为自适应参数,对感兴趣区域ROI中各灰度值的像素进行自适应伽马校正AGCWD,得到经局部增强后的结果图。
8.根据权利要求7所述的视网膜OCT图像水平配准和图像增强方法,其特征在于,所述的加权直方图参数pdfw(l)的计算公式为:
Figure FDA0004227484180000031
所述加权分布函数cdfw(l)的计算公式为:
Figure FDA0004227484180000032
式中l为图像的灰度值,pdfmax与pdfmin分别为直方图统计结果的最大值与最小值,α为调整参数。
9.根据权利要求8所述的视网膜OCT图像水平配准和图像增强方法,其特征在于,所述自适应伽马校正AGCWD的计算公式为:
Figure FDA0004227484180000033
式中T(l)即为每个灰度值经对比度增强后的结果。
CN202011237202.XA 2020-11-09 2020-11-09 一种视网膜oct图像水平配准和图像增强方法 Active CN112330638B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011237202.XA CN112330638B (zh) 2020-11-09 2020-11-09 一种视网膜oct图像水平配准和图像增强方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011237202.XA CN112330638B (zh) 2020-11-09 2020-11-09 一种视网膜oct图像水平配准和图像增强方法

Publications (2)

Publication Number Publication Date
CN112330638A CN112330638A (zh) 2021-02-05
CN112330638B true CN112330638B (zh) 2023-06-16

Family

ID=74316867

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011237202.XA Active CN112330638B (zh) 2020-11-09 2020-11-09 一种视网膜oct图像水平配准和图像增强方法

Country Status (1)

Country Link
CN (1) CN112330638B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115409689B (zh) * 2021-05-28 2023-09-29 南京博视医疗科技有限公司 一种多模态视网膜眼底图像的配准方法及装置
CN114092464B (zh) * 2021-11-29 2024-06-07 唯智医疗科技(佛山)有限公司 Oct图像的处理方法及装置
CN115670370B (zh) * 2022-12-29 2023-04-07 汕头大学·香港中文大学联合汕头国际眼科中心 一种去除眼底图像玻璃体混浊斑的视网膜成像方法和装置

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103810709A (zh) * 2014-02-25 2014-05-21 南京理工大学 基于血管的眼底图像与sd-oct投影图像配准方法
CN106489152A (zh) * 2014-04-10 2017-03-08 Sync-Rx有限公司 在存在医学设备的情况下的图像分析
CN106504228A (zh) * 2016-09-30 2017-03-15 深圳市莫廷影像技术有限公司 一种眼科oct图像的大范围高清快速配准方法及装置
CN107240128A (zh) * 2017-05-09 2017-10-10 北京理工大学 一种基于轮廓特征的x线片和彩色照片配准方法
CN109272507A (zh) * 2018-07-11 2019-01-25 武汉科技大学 基于结构随机森林模型的相干光断层图像的层分割方法
CN109934887A (zh) * 2019-03-11 2019-06-25 吉林大学 一种基于改进的脉冲耦合神经网络的医学图像融合方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103810709A (zh) * 2014-02-25 2014-05-21 南京理工大学 基于血管的眼底图像与sd-oct投影图像配准方法
CN106489152A (zh) * 2014-04-10 2017-03-08 Sync-Rx有限公司 在存在医学设备的情况下的图像分析
CN106504228A (zh) * 2016-09-30 2017-03-15 深圳市莫廷影像技术有限公司 一种眼科oct图像的大范围高清快速配准方法及装置
CN107240128A (zh) * 2017-05-09 2017-10-10 北京理工大学 一种基于轮廓特征的x线片和彩色照片配准方法
CN109272507A (zh) * 2018-07-11 2019-01-25 武汉科技大学 基于结构随机森林模型的相干光断层图像的层分割方法
CN109934887A (zh) * 2019-03-11 2019-06-25 吉林大学 一种基于改进的脉冲耦合神经网络的医学图像融合方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
"眼底OCT 图像中糖尿病性黄斑水肿的分割";何*** 等;《光电工程》;第45卷(第7期);第170605-1至170605-8页 *
"胸主动脉CTA与X线透视图配准研究";郝永达;《信息科技辑》;第1-5章 *
pan lingjiao etal.."OCTRexpert:a feature-based 3d registration method for retinal OCT images".《IEEE》.2020,第3885-3897页. *

Also Published As

Publication number Publication date
CN112330638A (zh) 2021-02-05

Similar Documents

Publication Publication Date Title
CN112330638B (zh) 一种视网膜oct图像水平配准和图像增强方法
CN106558030B (zh) 三维大视野扫频光学相干断层成像中脉络膜的分割方法
CN109345469A (zh) 一种基于条件生成对抗网络的oct成像中散斑去噪方法
CN104203081A (zh) 将复数的眼睛影像结合成为多聚焦全光影像的方法
WO2003053228A2 (en) Method and apparatus for eye registration
US20110134393A1 (en) Image processing apparatus, image processing method, and program storage medium
CN109829942A (zh) 一种眼底图像视网膜血管管径自动量化方法
CN111710012B (zh) 一种基于两维复合配准的octa成像方法与装置
CN108836257A (zh) 一种眼底oct图像中视网膜分层方法
EP3591614A1 (en) Method and computer program for segmentation of optical coherence tomography images of the retina
CN110648302B (zh) 基于边缘增强引导滤波的光场全聚焦图像融合方法
CN109325955B (zh) 一种基于oct图像的视网膜分层方法
CN117391955A (zh) 基于多帧光学相干层析图像的凸集投影超分辨率重建方法
CN114092405A (zh) 一种针对黄斑水肿oct图像的视网膜层自动分割方法
CN115393239A (zh) 一种多模态眼底图像配准和融合方法及***
CN110033496B (zh) 时间序列三维视网膜sd-oct图像的运动伪差校正方法
Shi et al. Automated choroid segmentation in three-dimensional 1-μ m wide-view OCT images with gradient and regional costs
CN111861977A (zh) 一种基于机器视觉的眼前节断层图像的特征提取方法
CN115690183A (zh) 一种图像配准数据处理方法及***
EP2693397B1 (en) Method and apparatus for noise reduction in an imaging system
CN109003284A (zh) 基于层厚统计信息模型的相干光断层扫描图像的层分割方法
Nourrit et al. Blind deconvolution for high-resolution confocal scanning laser ophthalmoscopy
CN110852977B (zh) 融合边缘灰度直方图与人眼视觉感知特性的图像增强方法
Li et al. MCFDFusion: Multi-focus image fusion based on multi-scale cross-difference and focus detection
CN111652847B (zh) 一种基于机器视觉的睑板腺图像预处理的方法

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