CN1826080A - 用于激光外科手术和其它光学应用的迭代傅立叶重建 - Google Patents

用于激光外科手术和其它光学应用的迭代傅立叶重建 Download PDF

Info

Publication number
CN1826080A
CN1826080A CNA2004800209971A CN200480020997A CN1826080A CN 1826080 A CN1826080 A CN 1826080A CN A2004800209971 A CNA2004800209971 A CN A2004800209971A CN 200480020997 A CN200480020997 A CN 200480020997A CN 1826080 A CN1826080 A CN 1826080A
Authority
CN
China
Prior art keywords
eyes
optical
gradient
optical texture
gradient fields
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
CNA2004800209971A
Other languages
English (en)
Other versions
CN100446718C (zh
Inventor
迪米特利·切恩亚克
查尔斯·E.·坎贝尔
艾里克·格罗斯
西玛·萨曼尼
戴光明
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.)
AMO Manufacturing USA LLC
Original Assignee
Visx Inc
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 Visx Inc filed Critical Visx Inc
Publication of CN1826080A publication Critical patent/CN1826080A/zh
Application granted granted Critical
Publication of CN100446718C publication Critical patent/CN100446718C/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B3/00Apparatus for testing the eyes; Instruments for examining the eyes
    • A61B3/10Objective types, i.e. instruments for examining the eyes independent of the patients' perceptions or reactions
    • A61B3/1015Objective types, i.e. instruments for examining the eyes independent of the patients' perceptions or reactions for wavefront analysis
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01JMEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
    • G01J9/00Measuring optical phase difference; Determining degree of coherence; Measuring optical wavelength

Abstract

提供使用傅立叶变换算法为光学组织系统确定光学表面模型的方法、系统和软件。重建眼睛的光学组织的方法包括通过眼睛的光学组织传输图像的步骤。穿过眼睛的光学组织从被传输图像测量表面梯度。应用傅立叶变换算法于表面梯度,以重建与眼睛的光学组织相对应的光学表面模型。

Description

用于激光外科手术和其它 光学应用的迭代傅立叶重建
相关申请的交叉引用
本申请是在2003年6月20日提交的美国专利申请No.10/601048(律师备审案件No.018158-021800US)的部分继续申请,并涉及在2004年3月3日提交的美国临时专利申请No.60/550514(律师备审案件No.018158-024300US),其公开的全部内容在此包含引作参考。
技术领域
本发明通常涉及测量光学系统的光学误差。更具体而言,本发明涉及为眼睛的光学组织系统确定光学表面模型的改进的方法和系统、用于重建眼睛的光学组织的波阵面表面/立视图(elevation map)的改进的方法和系统、和用于计算烧蚀模式的改进的系统。
背景技术
已知的激光眼科手术过程通常使用紫外或红外激光,以从眼睛的角膜上去除基质组织的显微层。激光典型去除选择形状的角膜组织,以经常用来校正眼睛的折射误差。紫外激光烧蚀导致角膜组织的光解作用,但是通常不会对眼睛的附近组织和下部组织造成重大的热损伤。受照射的分子以光化学的方式侵入到较小的挥发性碎片中,从而直接破坏分子间的键。
激光烧蚀过程可以去除角膜的目标基质,以为校正近视、远视、散光等各种目的改变角膜的外形。可以由各种系统和方法提供对于角膜上的烧蚀能量的分布的控制,包括使用可烧蚀掩模、固定或可移动的孔径、受控扫描系统、眼睛移动跟踪机构等。在已知的系统中,激光束常常包含一系列离散的激光能量脉冲,被去除组织的总形状和总量由在角膜上照射的激光能量脉冲的形状、大小、位置和/或数量决定。各种算法可以用来计算用于对角膜进行整形以校正眼睛折射误差的激光脉冲的模式。已知的系统利用各种形式的激光器和/或激光能量影响校正,包含红外激光器、紫外激光器、飞秒激光器、波长倍增固态激光器等。替代的视力校正技术利用角膜的径向切开、眼内透镜、可移动的角膜支撑结构等。
已知的角膜校正处理方法通常已经在校正诸如近视、远视、散光等的标准视力误差方面取得成功。但是,如所有的成功一样,仍希望得到进一步的改进。为此,现在可获得波阵面测量系统,以精确地测量特定患者的眼睛的折射特性。一种示意性波阵面技术系统是VISXWaveScanSystem,它使用Hartmann-Shack波阵面透镜组阵列,该波阵面透镜组阵列可量化患者眼睛的整个光学系统的像差,包括一阶和二阶球柱(sphero-cylindrical)误差、慧形像差、与慧形像差有关的三阶和四阶像差、散光、以及球面像差。
可以使用眼睛的波阵面测量,产生可允许评估眼睛的整个光程的像差例如内部像差和角膜表面上像差的眼睛像差图、高阶像差图或波阵面立视图。然后,使用该像差图计算允许外科激光系统校正患者的眼睛中或眼睛上的复杂像差的定制烧蚀模式。使用波阵面传感器数据计算定制烧蚀图像的已知方法通常包含用展开级数技术在数学上使眼睛的光学表面模型化。
人眼像差的波阵面或光程差(OPD)的重建对各种使用情况有益。例如,波阵面图、波阵面折射、点展开函数和处理表(treatment table)均可取决于重建的波阵面。
公知的波阵面重建可分类为两种方法:区域(zonal)重建和模态(modal)重建。区域重建用于早期的自适应光学系统。最近,由于Zernike多项式的使用,模态重建开始流行。可以通过已知的拟合技术导出Zernike多项式的系数,并且,如数学级数展开模型所指示,可通过使用眼睛的光学表面的形状确定折射校正过程。
表面重建的Zernike函数方法及其用于一般眼睛的精度具有限度。例如,6阶Zernike多项式不能精确地表示在所有情况下的实际波阵面。对于具有圆锥角膜的眼睛,这种差异可能会非常大。已知的Zernike多项式模型化方法还可以导致误差或“噪声”,该误差可导致小于理想的折射校正。并且,已知的表面模型化技术有些不直接,并且可能导致不必要的计算误差,以及导致缺乏对要执行的物理校正的理解。
鉴于此,所以希望提供对眼睛的光学组织进行数学模型化的改进的方法和系统。
发明内容
本发明提供新颖的能够用于解决丢失、错误或其它不充分的数据点的迭代傅立叶变换方法和系统。本发明还提供可用于确定最佳或合理的迭代次数的确定的目标或度量(metrics)。而且,本发明提供使用傅立叶变换算法在光学系统中测量误差并重建波阵面立视图的系统、软件和方法。
在第一方面中,本发明提供一种确定眼睛的光学组织系统的光学表面模型的方法,该方法包括以下步骤:a)从眼睛的光学组织系统输入光学数据,该光学数据包括一组局部梯度;b)基于该组梯度建立第一组合梯度场;c)通过应用傅立叶变换于第一组合梯度场导出第一重建波阵面;d)基于该第一重建波阵面提供第一修正梯度场;e)基于该第一修正梯度场建立第二组合梯度场;f)通过应用傅立叶变换于第二组合梯度场导出第二重建波阵面;和g)基于该第二重建波阵面确定光学表面模型。
在某方面中,所述光学数据包含波阵面数据。在某方面中,所述波阵面数据包含孔径内的一组局部或表面梯度。在某方面中,所述孔径与眼睛的瞳孔相对应。
在其他方面中,所述第一组合梯度场包含第一外部梯度场和测量梯度场,从而使该测量梯度场设置在第一外部梯度场的内部,并且该测量梯度场与所述组局部梯度相对应;所述第二组合梯度场包含第二外部梯度场和测量梯度场,从而使该第二外部梯度场对应于第一修正梯度场,并且使该测量梯度场设置在第二外部梯度场的内部。
在相关的方面中,确定眼睛光学组织系统的光学表面模型的方法还包括:选择第二重建波阵面中的至少一部分,提供最终的重建波阵面,并基于该最终的重建波阵面确定光学表面模型。
在相关的方面中,该方法还包括反复执行上述步骤(d)到(f),以通过应用傅立叶变换于第n个组合梯度场导出第n个重建波阵面;选择第n个重建波阵面中的至少一部分,提供最终的重建波阵面;和基于该最终的重建波阵面确定光学表面模型。
在相关的方面中,本发明提供一种绘制眼睛的波阵面误差图的方法。该方法包括如上所述确定光学表面模型,和基于光学表面模型绘制眼睛的波阵面误差图。在其他的相关方面中,本发明提供一种计算眼睛的光学组织系统的校正烧蚀模式的方法。该方法包括如上所述确定光学表面模型,和基于光学表面模型计算校正烧蚀模式。
在其他方面中,所述第一组合梯度场包含第一外部梯度场和测量梯度场,从而使得该测量梯度场设置在第一外部梯度场的内部,并且使该测量梯度场与所述组局部梯度相对应;以及所述第二组合梯度场包含第二外部梯度场和测量梯度场,从而使该第二外部梯度场对应于第一修正梯度场,并且使测量梯度场设置在第二外部梯度场的内部。
在相关的方面中,确定眼睛光学组织系统的光学表面模型的方法还包括选择第二重建波阵面中的至少一部分,提供最终的重建波阵面,并基于该最终的重建波阵面确定光学表面模型。
在相关的方面中,该方法还包括反复执行上述步骤(d)到(f),以通过应用傅立叶变换于第n个组合梯度场导出第n个重建波阵面;选择第n个重建波阵面中的至少一部分,提供最终的重建波阵面;和基于该最终的重建波阵面确定光学表面模型。
在相关的方面中,本发明提供一种改变眼睛的光学组织表面的方法。该方法包括如上所述计算校正烧蚀模式;和通过将激光烧蚀施加到眼睛上,根据校正烧蚀模式改变光学组织表面。
在第二方面中,本发明提供一种用于确定眼睛的光学组织系统的光学表面模型的系统。该系统包括a)用于将图像传输穿过光学组织系统的光源;b)被定向成通过检测被传输的图像为光学组织系统确定一组局部梯度的传感器;c)与传感器耦合的处理器;和d)与处理器耦合的存储器,该存储器被配置为存储多个由处理器执行的代码模块。所述多个代码模块包括:i)用于基于一组局部梯度建立第一组合梯度场的模块;ii)用于通过应用傅立叶变换于第一组合梯度场导出第一重建波阵面的模块;iii)用于基于第一重建波阵面提供第一修正梯度场的模块;iv)用于基于第一修正梯度场建立第二组合梯度场的模块;v)用于通过应用傅立叶变换于第二组合梯度场导出第二重建波阵面的模块;和vi)用于基于第二重建波阵面为眼睛的光学组织系统确定光学表面模型的模块。在本方面以及整个申请中描述的模块可以包括数据处理软件和/或硬件,并可与其它数据处理结构集成。
在相关的方面中,本发明提供一种绘制眼睛的光学组织系统的波阵面误差的系统。该系统包括用于为眼睛的光学组织系统确定光学表面模型的上述系统、以及基于该光学表面模型绘制眼睛的波阵面误差图的模块。在进一步的相关方面中,本发明提供一种为眼睛的光学组织系统计算校正烧蚀模式的系统。该系统包括为眼睛的光学组织系统确定光学表面模型的上述系统、以及基于光学表面模型为眼睛计算校正烧蚀模式的模块。在再一个方面中,本发明提供一种用于改变眼睛的光学组织表面的系统。该系统可包括为眼睛的光学组织系统计算校正烧蚀模式的上述系统,以及基于该校正烧蚀模式改变眼睛的光学组织表面的激光器。
在第三方面中,本发明提供一种用于确定对应于眼睛的光学组织系统的光学表面模型的系统。该系统包括a)用于将图像传输穿过光学组织系统的的装置;b)在来自图像的光程中、基于被传输图像模式为光学组织系统确定一组局部梯度的装置;c)与局部梯度确定装置耦合的、基于一组局部梯度建立第一组合梯度场的装置;d)与第一组合梯度场建立装置耦合的、基于应用傅立叶变换于第一组合梯度场导出第一重建波阵面的装置;e)与第一重建波阵面导出装置耦合的、基于第一重建波阵面提供第一修正梯度场的装置;f)与第一修正梯度场提供装置耦合的、基于第一修正梯度场建立第二组合梯度场的装置;g)与第二组合梯度场建立装置耦合的、基于应用傅立叶变换于第二组合梯度场导出第二重建波阵面的装置;和h)与第二重建波阵面导出装置耦合的、基于第二重建波阵面为眼睛的光学组织系统确定光学表面模型的装置。
在第四方面中,本发明提供一种存储在计算机可读存储介质上的计算机程序。该计算机程序包括a)用于将图像传输穿过眼睛的光学组织系统的代码;b)用于基于被传输图像为眼睛的光学组织系统确定一组局部梯度的代码;c)用于基于一组局部梯度建立第一组合梯度场的代码;d)用于基于应用傅立叶变换给第一组合梯度场导出第一重建波阵面的代码;e)用于基于第一重建波阵面提供第一修正梯度场的代码;f)用于基于第一修正梯度场建立第二组合梯度场的代码;g)用于基于应用傅立叶变换给第二组合梯度场导出第二重建波阵面的代码;和h)用于基于第二重建波阵面为眼睛的光学组织系统确定光学表面模型的代码。在相关的方面中,该计算机程序还包括:用于基于光学表面模型计算校正烧蚀模式的代码。在另一相关方面中,计算机程序还可包括用于基于校正烧蚀模式向眼睛传送激光能量的代码。
在一个方面中,本发明提供一种重建眼睛的光学组织的方法。该方法包括将图像传输穿过眼睛的光学组织。在眼睛的光学组织上测量来自被传输图像的表面梯度。对表面梯度应用傅立叶变换算法,重建与眼睛的光学组织相对应的表面。
在一些实施例中,由Hartmann-Shack传感器组件实施表面梯度数据的测量。图像作为多个光束组(beamlet)被光学组织传输,表面梯度将是梯度的阵列的形式。每一梯度与眼睛的光学组织的相关部分相对应,每一光束组根据相应的梯度传输穿过光学组织。在这些实施例中,被测量的表面将是波阵面表面或波阵面立视图的形式。
在一个实施例中,傅立叶变换算法可以应用快速傅立叶变换或离散傅立叶分解和逆离散傅立叶变换。一些傅立叶变换算法可以包括平均梯度场,以便从重建的表面去除斜角。与基于单位圆的Zernike多项式重建不同,傅立叶变换在重建中使用所有可获得的信息。
如傅立叶重建的表面所指示,可以计算基于眼睛的重建光学组织计算的校正烧蚀模式并与眼睛对准。该校正烧蚀模式典型包括在光学组织的高度中提出的变化,以实现所希望的眼睛光学特性的变化。在导出校正烧蚀模式之后,可以使用激光烧蚀改变光学组织表面。
在另一方面中,本发明提供一种测量眼睛的光学组织的方法。该方法包括穿过光学组织传输图像。从该传输的图像确定穿过光学组织的局部梯度。通过将傅立叶变换算法应用于穿过眼睛的光学组织的表面梯度绘制眼睛的波阵面误差图。
通过Hartmann-Shack传感器组件可以实施局部梯度的测量。可以向波阵面误差添加平均梯度场以校正斜角。在绘制波阵面误差图之后,可以基于绘制的眼睛光学组织的波阵面误差图产生激光烧蚀处理表,并且根据校正烧蚀模式通过激光烧蚀模式改变光学组织表面。
在另一方面中,本发明还提供一种用于测量眼睛组织的波阵面误差的系统。该系统包括与处理器耦合的存储器。存储器可被配置为存储多个用于由处理器执行的代码模块。该多个代码模块包括:用于穿过光学组织传输图像的模块;用于从该传输的图像确定穿过光学组织的局部梯度的模块;和用于通过将傅立叶变换算法应用于穿过眼睛的光学组织的表面梯度绘制眼睛的波阵面误差图的模块。
该系统还包括用于穿过眼睛的光学组织传输源图像的、与处理器耦合的图像源。该图像可以为光点或小斑点、或任何其它适当的图像。该系统还可以包括与处理器耦合的波阵面传感器系统诸如Hartmann-Shack传感器组件。
该系统可以包括一个或更多个跟踪光学组织的位置的照相机。在这些实施例中,多个代码模块还包括用于记录关于由照相机得到的光学组织图像的波阵面误差的模块。
在一些实施例中,该系统可以包括基于光学组织的重建表面计算定制的激光烧蚀程序或烧蚀表的模块。激光系统可与上述测量系统通信。该激光系统可以包括可编程用于根据校正烧蚀模式向光学组织传送激光能量的激光器。
在另一方面中,本发明提供一种存储在计算机可读存储介质上的计算机程序。该计算机程序包括用于将图像传输穿过眼睛的光学组织的代码,用于从穿过眼睛的光学组织的传输图像测量表面梯度的代码,和用于通过将傅立叶变换算法应用于穿过眼睛的光学组织的表面梯度绘制眼睛的波阵面误差图的代码。
在一些实施例中,计算机程序可以包括用于基于由傅立叶重建表面表示的眼睛的光学组织计算校正烧蚀模式的代码,用于根据校正烧蚀模式向光学组织传送激光能量的代码,和/或用于将绘制的波阵面误差图与眼睛的光学组织的图像对准的代码。
本发明的方法和设备可以以用于这些用途的一个或更多工具箱来提供。例如,这些工具箱可包括用于确定与眼睛的光学组织系统相对应的光学表面模型的系统。可选地,这些工具箱还可以包括与本发明相关描述的任何其它系统组件以及与本发明相关的任何其它资料或项目。使用说明能阐述上述方法中的任何之一。还应理解的是,根据本发明的系统可以配置成执行上述方法步骤中的任何之一。
在以下附图、说明书和权利要求书中,这些和其它方面将会显而易见。
附图说明
图1示例了根据本发明实施例的激光烧蚀系统。
图2示例了根据本发明实施例的简化的计算机系统。
图3示例了根据本发明实施例的波阵面测量系统。
图3A示例了根据本发明实施例的另一个波阵面测量系统。
图4图示了执行本发明的一种方法的一组简化模块。
图5是图示使用傅立叶变换算法确定角膜烧蚀处理程序的方法的流程图。
图6图示了对在6mm瞳孔上具有+2烧蚀的表面的直接积分重建、6阶Zernike多项式重建、10阶Zernike多项式重建和傅立叶变换重建的比较。
图7图示了对远视眼表面的直接积分重建、6阶Zernike多项式重建、10阶Zernike多项式重建和傅立叶变换重建的比较。
图8图示了对另一远视眼表面的直接积分重建、6阶Zernike多项式重建、10阶Zernike多项式重建和傅立叶变换重建的比较。
图9示例了从由傅立叶变换重建算法(F梯度)、Zernike多项式重建算法(Z梯度)、直接积分重建算法(D梯度)重建的波阵面和直接测量的波阵面计算的梯度场的差异。
图10示例了傅立叶、10阶Nernike和直接积分重建的重建波阵面的强度图。
图11示例了傅立叶、6阶Nernike和直接积分重建的重建波阵面的强度图。
图12示例了根据本发明的一个实施例的算法流程图。
图13示例了根据本发明的一个实施例的波阵面重建的表面图。
图14示例了根据本发明的一个实施例的波阵面重建的表面图。
图15示例了具有和不具有丢失数据的波阵面的比较。
具体实施方式
本发明提供了使用高速和精确的傅立叶或迭代傅立叶变换算法以数学地为眼睛的光学组织系统确定光学表面模型或另外数学地重建眼睛的光学组织的系统、软件、和方法。
本发明通常对提高诸如光折射角膜切除术(PRK)、光疗角膜切除术(PTK)和激光原位角膜上皮磨镶术(keratomileusis)(LASIK)的激光眼科手术过程的精度和效率有用。本发明通过改善用于测量眼睛的光学误差的方法能够提供增强光学精度的折射过程,并由此计算更精确的折射烧蚀程序。在一个特定的实施例中,本发明涉及对有病的眼睛进行基于处理波阵面的烧蚀。
本发明可容易地适用于现有的激光系统、波阵面测量系统和其它光学测量设备。本发明通过提供更直接的测量并校正光学系统误差的方法(因此,倾向于产生很少的噪声和其它误差),可以方便对角膜的整形,从而使被处理的眼睛有规律地超过希望视力的正常20/20阈值。虽然主要以激光眼科手术系统为上下文描述本发明的系统、软件和方法,但是应该理解,本发明也可适用于诸如眼镜、眼内透镜、隐形眼镜、角膜环移植、胶原角膜组织热改型等的眼睛处理过程和系统。
现在参照图1,本发明的激光眼科手术系统10包括产生激光束14的激光器12。激光器12光耦合至将激光束14引向患者P眼睛的激光传输光学器件16。传输光学器件支撑结构(这里为了清楚起见未示出)从支撑激光器12的框架18延伸。显微镜20安装在传输光学器件支撑结构上,显微镜常常用于对眼睛的角膜成像。
激光器12通常包括准分子激光器,该准分子激光器在理想情况下包括产生具有约193nm波长的激光脉冲的氩氟激光器。激光器12将优选设计用于在患者的眼睛上提供经由激光传输光学器件16传输的反馈稳定的注量(fluence)。本发明也可使用替代的紫外或红外辐射源,特别是可以在不对眼睛的附近和/或下层组织造成重大损伤的情况下适于可控地烧蚀角膜组织的辐射源。在替换的实施例中,激光束源使用如授权给Lin的美国专利No.5520679和No.5144630和授权给Mead的No.5742626所述具有193~215nm波长的固态激光源,其公开的全部内容在此包含引作参考。在另一个实施例中,激光源是授权给Telfair的No.5782822和No.6090102所述红外激光器,其公开的全部内容在此包含引作参考。因此,虽然准分子激光器是示例性的烧蚀束的来源,但本发明可使用其它的激光器。
激光器12和激光传输光学器件16通常在计算机系统22的指引下将激光束14引向患者P的眼睛。计算机系统22将常常选择性地调节激光束14,以将角膜的部分暴露于激光能量的脉冲下,以便实现对角膜的预定整形,并改变眼睛的折射特性。在许多的实施例中,激光器12和激光传输光学系统16将均在计算机系统22的控制下,影响希望的激光整形处理,同时该计算系统影响(并可选地修改)激光脉冲的模式。可以以处理表的形式在实体介质29的机器可读数据中汇总脉冲的模式,并且可根据响应于从烧蚀监视系统反馈系统提供的实时反馈数据,从自动图像分析系统输入到计算机系统22的反馈(或由系统操作员手动输入处理器)调整处理表。激光处理系统10和计算机系统22可响应该反馈继续和/或终止整形处理,并且也可以有选择地至少部分地基于反馈修改计划的整形。
如本领域的普通技术人员应该理解,激光系统10可以包括其它组件和子系统。例如,如美国专利No.5646791所描述的,可以包括空间和/或时间积分器,以控制激光束内的能量分布,其公开的全部内容在此包含引作参考。烧蚀流出排放器/过滤器、吸气器和激光外科系统的其它辅助部件是本领域中公知的。在共同受让的美国专利No.4665913、4669466、4732148、4770172、4773414、5207668、5108388、5219343、5646791和5163934中可能够找到用于执行激光烧蚀过程的适当系统的进一步细节,其公开的全部内容在此包含引作参考。适当的系统还包括诸如由Alcon、Bausch&Lomb、Nidek、WaveLight、LaserSight、Schwind、和Zeiss-Meditect等制造和/或销售的市场上可获得的折射激光系统。
图2是可由本发明的激光外科系统10使用的示意性计算机系统22的简化方框图。计算机系统22典型地包括可以经由总线子系统54与大量的外设通信的至少一个处理器52。这些外设可以包括包含内存子系统58和文件存储子系统60的存储子系统56、用户接口输入设备62、用户接口输出设备64、和网络接口子系统66。网络接口子系统66提供与外部网络68和/或其它设备诸如波阵面测量系统30的接口。
用户接口输入设备62可以包括键盘、诸如鼠标、跟踪球、接触板或图形输入板的指针设备、扫描仪、脚踏板、操纵杆、包含在显示器中的接触屏、诸如语音识别系统的音频输入设备、麦克风和其它类型的输入设备。用户输入设备62将常常用于从实体存储介质29下载实施本发明任何方法的计算机可执行代码。一般而言,使用术语“输入设备”目的在于包括各种常规和专有设备和向计算机系统22输入信息的方式。
用户接口输出设备64可以包括显示子系统、打印机、传真机、或诸如音频输出设备的非可视化显示器。显示子系统可以是阴极射线管(CRT)、诸如液晶显示器(LCD)的平板显示设备、投影设备等。该显示子系统可以诸如经由音频输出设备提供非可视化显示。一般而言,使用术语“输出设备”目的在于包括各种常规和专有设备和从计算机系统22向用户输出信息的方式。
存储子系统56存储提供本发明各个实施例的功能的基本程序设计和数据结构。例如,可以在存储子系统56中存储这里所述的实现本发明方法的功能的数据库和模块。这些软件模块一般由处理器52执行。在分布式环境中,软件模块可存储在多个计算机系统上并由多个计算机系统的处理器执行。存储子系统56典型包括内存子系统58和文件存储子系统60。
内存子系统58典型包括大量的存储器,这些存储器包括用于在程序执行期间存储指令和数据的主随机存取存储器(RAM)70和存储固定指令的只读存储器(ROM)72。文件存储子系统60提供程序和数据文件的持久性(非易失性)存储,并包括可选的实体存储介质29(图1),这些介质可以包括有波阵面传感器数据、波阵面梯度、波阵面立视图、处理图、和/或烧蚀表。文件存储子系统60可以包括硬盘驱动器、软盘驱动器以及相关的可移动介质、光盘只读存储器(CD-ROM)驱动器、光驱动器、DVD、CD-R、CD-RW、固态移动存储器和/或其它移动介质盒式磁盘或盘。一个或更多的驱动器可位于与计算机系统22耦合的其它地点上的其它相连接计算机上的远程位置上。实现本发明功能的模块可由文件存储子系统60存储。
总线子系统54提供让计算机系统22的各个组件和子系统按照目的相互通信的机制。计算机系统22的各个子系统和组件不需要位于相同的物理位置,而是可以分布在分布式网络内的各个位置上。虽然总线子系统54被示意性示出为单条总线,但总线子系统的替换实施例可使用多条总线。
计算机系统22本身可以为各种的类型,包括个人计算机、便携式计算机、工作站、计算机终端、网络计算机、波阵面测量系统或激光外科系统中的控制系统、主计算机或任何其它数据处理系统。由于计算机和网络的不断变化的本性,图2中示出的计算机系统22的描述仅仅是为了示例本发明的一个实施例的具体实例。计算机系统22的许多其他配置有可能具有比图2中所示的计算机系统或多或或少的组件。
现在参照图3,该图以简化的形式示意性地示出了波阵面测量系统30的一个实施例。一般而言,波阵面测量系统30被配置为检测离开患者眼睛的梯度图的局部斜率。基于Hartmann-Shack原理的设备通常包含在孔径上均匀地对梯度图取样的透镜组阵列,该孔径典型为眼睛的出射光瞳。之后,分析梯度图的局部斜率,以便重建波阵面表面或图。
更具体而言,一个波阵面测量系统30包括诸如激光器的图像源32,该图像源32通过眼睛E的光学组织34投射源图像,以在视网膜R的表面上形成图像44。来自视网膜R的图像通过眼睛的光学系统(例如光学组织34)传输,并通过系统光学器件37成像到波阵面传感器36上。波阵面传感器36将信号传送到用于测量光学组织34中光学误差和/或确定光学组织烧蚀处理程序的计算机系统22′。计算机系统22′可以包括与图1和图2中所示例的计算机系统22相同或相似的硬件。计算机系统22′可以与控制激光外科系统10的计算机系统22进行通信,或者,波阵面测量系统30和激光外科系统10的计算机系统22、22′的部分或全部组件可以是组合的或者是单独的。如果需要,可以经由实体介质29、I/O端口、诸如企业内部网或因特网等网络连接66从波阵面传感器36发送数据给激光计算机系统22。
波阵面传感器36通常包括透镜组阵列38和图像器传感器40。当来自视网膜40的图像通过光学组织34传输并在图像传感器40的表面上成像以及眼睛瞳孔P的图像类似地在透镜组阵列38的表面上成像时,透镜组阵列将该传输的图像分开为光束组的阵列42,并且(与系统的其它光学组件组合)使在传感器40的表面上成像分开的光束组。传感器40典型包括电荷耦合器件或“CCD”,并检测这些各个光束组的特性,使用该特性可确定光学组织34的相关区域的特性。特别地,在图像44包括光点或小斑点时,由光束组所成像的传输斑点的位置可直接指示光学组织的相关区域的局部梯度。
眼睛E通常定义前取向ANT和后取向POS。如图3所示,图像源32通常通过光学组织34将后取向的图像投影到视网膜R上。光学组织34再从视网膜向前向波阵面传感器36传输图像44。当图像源通过光学组织34原始地传输时,在视网膜R上实际形成的图像44可被眼睛的光学系统中的任何缺陷畸变。可选地,图像源投影光学器件46可配置为或适于减少图像44的任何畸变。
在一些实施例中,图像源光学器件46可通过补偿光学组织34的球面和/或柱面误差减少较低阶光学误差。光学组织的较更高阶的光学误差也可以通过使用诸如可变形镜(以下对此描述)的自适应光学元件得到补偿。使用选择的用于在视网膜R上定义图像44上的点或小斑点的图像源32可便于分析由波阵面传感器36提供的数据。由于瞳孔的中心部分与周边部分相比较少倾向于产生光学误差,因此通过穿过小于瞳孔50的光学组织34的中心区域48传输源图像,可以限制图像44的畸变。不管特定的图像源结构如何,通常有益的是在视网膜R上具有完美定义和精确形成的图像44。
在一个实施例中,波阵面数据可以以两个单独的阵列存储在波阵面传感器系统30的计算机可读介质29或存储器中,该阵列包含由瞳孔照相机51(图3)图像测量的、从Hartmann-Shack传感器图像的图像点分析得到的x和y波阵面梯度值、加上与Hartmann-Shack透镜组阵列的标称中心的x和y瞳孔中心偏离量。这种信息包含所有可获得的有关眼睛的波阵面误差的信息,并足以重建波阵面或其任何部分。在这些实施例中,不需要多于一次地重新处理Hartmann-Shack图像,并且存储梯度阵列所需要的数据空间不大。例如,要容纳具有8mm直径的瞳孔的图像,20×20大小的阵列(即,400个单元)常常足够。可以理解,在其它实施例中,波阵面数据可以以单个阵列或多个阵列存储在波阵面传感器系统的存储器中。
虽然通常参照检测图像44描述本发明的方法,但应该理解,可以进行一系列波阵面传感器数据读取。例如,波阵面数据读取的时间序列可帮助提供眼睛组织像差的更精确的全面确定。由于眼睛组织的形状可在很短的时间周期内变化,因此多个时间分开的波阵面传感器测量能够避免依赖光学特性的单个快照作为折射校正过程的基础。还可获得其他的替换,包括在使眼睛处于不同配置、位置和/或取向时获取眼睛的波阵面传感器数据。例如,如美国专利No.6004313中所描述的那样,患者将经常通过在固定目标上聚焦帮助保持眼睛与波阵面测量系统30对准,其公开的全部内容在此包含引作参考。如该参考文献所述通过改变固定目标的位置,可在眼睛调节或自适应成像不同距离和/或角度的视场的同时确定眼睛的光学特性。
可以通过参照由瞳孔照相机52提供的数据校验眼睛的光轴位置。在示意性的实施例中,瞳孔照相机52对瞳孔50进行成像,以便为记录与光学组织有关的波阵面传感器数据确定瞳孔的位置。
在图3A中示例波阵面测量系统的替换实施例。图3A的系统的主要组件与图3的主要组件相似。另外,图3A包括可变形镜形式的自适应光学元件53。源图像在向视网膜R传输的过程中被可变形镜98反射,并且可变形镜也沿着用于在视网膜R和成像传感器40之间形成传输图像的光程。可变形镜98可在计算机系统22的控制下变形,以限制在视网膜上形成的图像的畸变或由在视网膜上形成的图像形成的随后图像的畸变,并且可以提高所产生的波阵面数据的精度。在美国专利No.6095651中更全面地描述了图3A的系统的结构和使用,其公开的全部内容在此包含引作参考。
用于测量眼睛和烧蚀的波阵面测量系统的实施例的组件包括可从加利福尼亚州圣克拉拉市的VISX,INCORPORATED可获得的VISX WaveScan的元件。一个实施例包括具有上述可变形镜的WaveScan。在美国专利No.6271915中描述了波阵面测量系统的替换实施例,其公开的全部内容在此包含引作参考。
W.H.Southwell在“Wave-front estimation from wave-frontslope measurements,”J.Opt.Soc.Am.70:998-1006(1980)中详细地讨论了具有Zernike多项式的模态重建的使用以及模态和区域重建的比较。相关地,G.Dai在“Modal wave-front reconstruction withZernike polynomials and Karhunen-Loeve functions,”J.Opt.Soc.Am.13:1218-1225(1996)中提供了模态重建具有Zernike多项式的各种波阵面重建误差的详细分析。如Liang等人在“Objective Measurementof Wave Aberrations of the Human Eye with the Use of aHarman-Shack front Sensor,”J.Opt.Soc.Am.11(7):1949-1957(1994)中提出的那样,使用Zernike多项式使光学表面模型化。这些参考中每一篇的全部内容在此包含引作参考。
Schweigerling,J等人在“Using corneal height maps andpolynomial decomposition to determine corneal aberrations,”Opt.Vis.Sci.,第74卷,第11期(1997)中和Guirao,A等人在“Corneal waveaberration from videokeratography:Accuracy and limitation of theprocedure,”J.Opt.Soc.Am.第17卷第6期(2000)中为规则的角膜形状广泛地研究了正常眼睛的表面重建及其精读的Zernike函数方法。D.R.Ishkander等人在“An Alternative PolynomialRepresentation of the Wavefront Error Function,”IEEE Transactionson Biomedical Engineering,第49卷,第4期(2002)中报道了,当与Bhatia-Wolf多项式相比较时,6阶Zernike多项式重建方法提供了差的拟合。
模态波阵面重建典型包含将波阵面展开成一组基函数。由于Zernike多项式对环形瞳孔是一组完全且正交的函数,因此使用Zernike多项式作为波阵面展开基函数已在波阵面技术领域中接受。另外,诸如散焦、散光、慧形像差和球面像差的一些较低阶Zernike模代表经典的像差。不幸的是,使用Zernike多项式可能存在缺点。由于Zernike基函数在孔径附近具有快速的波动,特别是对于更高的阶次,因此Zernike系数的轻微的变化就会大大影响波阵列表面。并且,由于在低阶和高阶Zernike模之间平衡的像差,Zernike级数的断点常导致不一致的Zernike系数。
为了解决Zernike重建的上述问题,我们寻求其它基函数。由于具有稳定快速的傅立叶变换(FFT)算法,因此傅立叶级数看起来似乎是有利的基函数集。并且,傅立叶级数的导数仍是傅立叶级数。对于无界函数(即,没有边界条件),可以使用傅立叶直接从一组梯度数据估计函数。但是,由于波阵面重建一般涉及有界函数或具有瞳孔孔径的函数,因此难以直接将该技术应用于波阵面技术。
迭代傅立叶重建技术可应用于具有无限制的孔径函数的有界函数。也就是说,函数的孔径可以是圆形、环形、卵形线、正方形、矩形或任何其它形状。在Roddier等人的“Wavefront reconstructionusing iterative Fourier transforms,”Appl.Opt.30,1325-1327(1991)中讨论了这种方法,其公开的全部内容在此包含引作参考。但是,由于角膜反射、坏的CCD像素等等产生的丢失数据点,这种解决方法大大得到改善。
I.为眼睛的光学组织系统确定光学表面模型
本发明提供使用高速和精确的迭代傅立叶变换算法以为眼睛的光学组织系统数学地确定光学表面模型的系统、软件和方法。
A.从眼睛的光学组织系统输入光学数据
有各种用于从光学组织系统产生光学数据的设备和方法。像差镜或像差计的类别包括经典的综合屈光检查仪和波阵面方法。也可使用基于形态学的测量设备和方法产生光学数据。波阵面设备常用于同时测量光学组织系统的低阶和高阶像差。特别地,波阵面分析典型包括以下步骤:穿过眼睛的光学系统传输图像,并基于该传输的图像确定光学组织系统的一组表面梯度。这些表面梯度可用于确定光学数据。
B.通过将迭代傅立叶变换应用于光学数据确定光学表面模型
图4图示了用于执行根据本发明一个实施例的方法的一组简化模块。这些模块可以是处理器52(图2)所处理的计算机可读介质上的软件模块、硬件模块或其组合。波阵面像差模块80典型从波阵面传感器接收数据,并测量成像的整个光学组织系统的像差和其它光学特性。如上所述,典型通过经由光学组织传输图像(诸如光的小斑点或光点)产生来自波阵面传感器的数据。波阵面像差模块80产生光学梯度的阵列或梯度图。来自波阵面像差模块80的光学梯度数据可被传输到傅立叶变换模块82,在该傅立叶变换模块82中,光学表面或其模型或波阵面立面图可以根据光学梯度数据数学地重建。
应当理解,由于梯度数据将表示实际位于整个眼睛组织系统中的像差的影响,因此光学表面或其模型不需要精确地与实际的组织表面匹配。但是,施加于光学组织表面上校正从梯度导出的像差的校正应该校正光学组织系统。这里使用的术语诸如“光学组织表面”或“光学表面模型”可包含理论上的组织表面(例如从波阵面传感器数据导出)、实际的组织表面和/或为处理目的形成的组织表面(例如,通过在LASIK过程中切割角膜组织以替换角膜上皮层和基质的片并暴露下面的基质)。
一旦波阵面立面图由傅立叶变换模块82产生,波阵面梯度图就可以传输到激光处理模块84,用于产生处理处理或改进光学组织中的光学误差的激光烧蚀处理方案。
图5是示例一个基于傅立叶产生激光烧蚀处理方案的方法的数据流和方法步骤的详细流程图。图示的方法典型由包括处理器和与处理器相耦合的存储器的系统执行。该存储器可配置为存储具有用于执行该方法步骤的指令和算法的多个模块。
可以理解,本发明不应限于所示例出步骤顺序或特定的步骤,并且,在不背离本发明的范围的情况下,可以对该方法进行各种修改诸如具有更多或更少的步骤。为了进行比较,产生波阵面立视图的级数展开方法以虚线的方式示出,并且它们是可选的步骤。
包括波阵面传感器(诸如Hartmann-Shack传感器)的波阵面测量系统可用于获得眼睛光学组织的一个或更多个位移图90(例如,Hartmann-Shack位移图)。可以通过穿过眼睛的光学组织传输图像并检测正在离开的波阵面表面得到该位移图。
从位移图90,有可能计算在眼睛光学组织上的表面梯度或梯度图92(例如Hartmann-Shack梯度图)。梯度图92可以包括从Hartmann-Shack传感器的每个透镜组的每一位置计算的局部化梯度的阵列。
可以将傅立叶变换应用于梯度图,以数学重建光学组织或确定光学表面模型。傅立叶变换将典型地以波阵面立视图的形式输出重建的光学组织或光学表面模型。对本发明来说,术语傅立叶变换还包含迭代傅立叶变换。
已经发现诸如快速傅立叶变换(FFT)的傅立叶变换重建方法比当前使用的6阶Zernike或多项式重建方法快很多倍,并产生更精确的实际波阵面的重建。有利的是,傅立叶重建对于可获得的数据密度来说将重建中使用的特定频率限制到Nyquist极限,并给出没有混淆现象的较高分辨率。出于一些先验的原因,如果需要限制使用的空间频率,这可通过舍弃计算的傅立叶变换空间后半段中的梯度的变换完成。如果希望对可获得的波阵面的小部分取样或使其偏心,这可以通过在傅立叶变换操作前对梯度数据简单的掩模(mask)操作完成。与需要瞳孔尺寸和瞳孔集中的Zernike重建方法不同,这些考虑不影响快速傅立叶变换。
而且,由于波阵面传感器在有规律间隔的网格上测量梯度图的x和y分量,因此,数据受到频带限制,并且数据不包含比与装置中的透镜组的间隔(典型地,透镜组的间隔不大于约0.8mm和约0.1mm,并且典型为约0.4mm)对应的Nyquist率大的空间频率。由于数据在有规律间隔的笛卡尔网格上,因此诸如傅立叶变换的非径向(non-radial)重建方法良好地适用于带限的数据。
与傅立叶变换不同,当使用级数展开技术以根据梯度图92产生波阵面立视图100时,使用梯度图92和选择的展开级数96导出适当的展开级数系数98。用于将组织表面模型化的数学级数展开的特别有益的形式是Zernike多项式。使用包含0至6阶或0至10阶项的典型Zernike多项式组。例如可以通过使用标准的最小二乘法拟合技术确定用于每个Zernike多项式Zn的系数an。Zernike多项式系数an的数量可以受到限制(例如,限制为约28个系数)
当通常考虑方便光学表面的模型化以便产生立视图时,Zernike多项式(以及也许所有的级数展开)可引入误差。但是,将Zernike多项式与它们的系数组合并对Zernike系数99进行求和使得可以计算波阵面立视图100,并且在一些情况下,可以非常精确地重建波阵面立视图100。
已经发现在一些情况下,特别是在眼睛光学组织中的误差为球形的情况下,Zernike重建可以比傅立叶变换重建更精确。因此,在一些实施例中,本发明的模块可包括傅立叶变换模块94和Zernike模块96、98、99。在这些实施例中,可以通过比较模块(未示出)来比较由两个模块得到的重建表面,以确定两个模块中的哪一个模块提供更精确的波阵面立视图。然后可以由100、102使用更精确的波阵面立视图,以分别计算处理图和烧蚀表。
在一个实施例中,波阵面立视图模块100可从每一个模块计算波阵面立视图,并且可以从每一个波阵面立视图计算梯度场。在一个配置中,比较模块可应用价值函数以确定每一个梯度图和原始测量的梯度图之间的差。价值函数的一个例子是从以下等式发现的均方根梯度误差RMSgrad
RMS grad = Σ alldatapoints { ( ∂ W ( x , y ) / ∂ x - Dx ( x , y ) 2 ) + ( ∂ W ( x , y ) / ∂ y - Dy ( x , y ) 2 ) } N
这里,N是取样位置的数量
(x,y)是取样位置
W(x,y)/x是重建的波阵面梯度的x分量
W(x,y)/y是重建的波阵面梯度的y分量
Dx(x,y)是梯度数据的x分量
Dy(x,y)是梯度数据的y分量
如果根据Zernike重建的梯度图更精确,则使用Zernike重建。如果傅立叶重建更精确,则使用傅立叶重建。
在计算波阵面立视图后,可以从波阵面立视图100计算处理图102,以去除光学组织的规则(球面和/或柱面)和不规则误差。通过将处理图102与特定激光系统的激光烧蚀脉冲特性104组合,可以产生烧蚀脉冲位置、尺寸、形状和/或数量的烧蚀表106。
激光处理烧蚀表106可包括用于一系列脉冲中的每一激光束脉冲的、眼睛上的激光束的水平和垂直位置。在处理的过程中光束的直径可在约0.65mm至6.5mm之间变化。处理烧蚀表106典型包括,在几百个脉冲到五千或更多个脉冲之间,激光束脉冲的数量随被去除的材料的量和激光处理表使用的激光束直径变化。可通过对各个脉冲排序可选地优化烧蚀表106,以避免局部化加热,并在中断处理程序等的情况下最小化不规则的烧蚀。然后可通过激光烧蚀根据处理表106对眼睛进行烧蚀。
在一个实施例中,激光烧蚀表106可通过使用各种替换的机制调节激光束14,以产生想要的整形。可通过使用一个或多个可变孔径选择性地限制激光束14。在美国专利No.5713892中描述了具有可变光圈和可变宽度缝隙的示意性可变孔径系统,其公开的全部内容在此包含引作参考。如美国专利No.5683379和6203539所描述的,以及如在1999年3月22日提交的美国申请No.09/274999所描述的,可以通过改变与眼轴的激光点的尺寸和偏移量设计激光束,其公开的全部内容在此包含引作参考。
还有可能存在其他的替换方案,包括:例如在美国专利No.4665913(其公开的全部内容在此包含引作参考)中所描述的,在眼睛的表面上扫描激光束,并控制每一位置上的脉冲数量和/或停留时间;如在1995年6月6日提交的美国专利申请No.08/468898(其公开的全部内容在此包含引作参考)中所描述的,在烧蚀的激光束14的光程中使用掩模,以改变入射到角膜上的光束的轮廓;在角膜上扫描可变尺寸光束(一般由可变宽度缝隙和/或可变直径虹彩光阑(irisdiaphragm)控制)的混合轮廓扫描系统等。
在美国专利N0.6673062中描述了用于准备这种烧蚀表的一个示意性方法和系统,其公开的全部内容在此包含引作参考。
现在将描述根据本发明的一个实施例使用傅立叶变换算法从表面梯度数据进行表面重建的数学推导。这种数学算法典型地包含在傅立叶变换模块82(图4)、傅立叶变换步骤94(图5)或其它相当的软件或硬件模块中,以重建波阵面表面。可以理解,下述的傅立叶变换算法仅是例子,本发明不应限于该特定的实现方式。
首先,假设存在可由函数s(x,y)表示的表面,并假设存在给出该表面梯度
Figure A20048002099700312
的数据。目的是从该梯度数据求得表面s(x,y)。
假设该表面在所有空间上都局部可积分,从而它可由傅立叶变换表示。然后由下式给出该表面的傅立叶变换:
F ( s ( x , y ) ) = 1 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ s ( x , y ) e - i ( ux + vy ) dxdy = S ( u , v ) - - - ( 1 )
然后可利用下式根据变换系数S(u,v)重建表面:
s ( x , y ) = 1 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ S ( u , y ) e i ( ux + vy ) dudv - - - ( 2 )
现在可以利用等式(2)关于表面的傅立叶系数给出梯度的x分量 的表示:
∂ s ( x , y ) ∂ x = ∂ ( 1 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ S ( u , y ) e i ( ux + vy ) dudy ) ∂ x
然后,积分下的微分给出:
∂ s ( x , y ) ∂ x = 1 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ iuS ( u , v ) e i ( ux + vy ) dudv - - - ( 3 )
与(3)相类似的等式关于傅立叶系数给出梯度的y分量的表示:
∂ s ( x , y ) ∂ y = 1 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ ivS ( u , v ) e i ( ux + vy ) dudv - - - ( 4 )
梯度的x分量是x和y的函数,因此它也可以通过从傅立叶变换得到的系数表示。使 dx ( x , y ) = ∂ s ( x , y ) ∂ x , 从而,按照导出(2)的逻辑
dx ( x , y ) = 1 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ Dx ( u , v ) e i ( ux + vy ) dudv = ∂ s ( x , y ) ∂ x - - - ( 5 )
这里,
F ( dx ( x , y ) ) = ∫ - ∞ ∞ ∫ - ∞ ∞ dx ( x , y ) e - i ( wx + vy ) dxdy = Dx ( u , v ) - - - ( 6 )
等式(3)必须等于等式(5),并且通过检查相似的项可以看出,如果
Dx(u,v)=iuS(u,v)
或者
S ( u , v ) = - iDx ( u , v ) u - - - ( 7 )
通常这才是成立的。
使用(4)进行的y梯度分量的类似推导得出
S ( u , v ) = - iDy ( u , v ) v - - - ( 8 )
注意,(7)和(8)表示Dx(u,v)和Dy(u,v)在函数上取决于关系:
vDx(u,v)=uDy(u,v)
现在可以通过第一次执行两个梯度场dx和dy的离散傅立叶分解根据梯度数据重建表面,以产生离散傅立叶梯度系数Dx(u,v)和Dy(u,v)。使用这些分量(7)和(8)以求得表面S(u,v)的傅立叶系数。这些反过来与逆离散傅立叶变换一起使用以重建表面s(x,y)。
上述处理进行离散傅立叶梯度系数的非对称使用,其中,使用一个或另一个离散傅立叶梯度系数以求得表面的傅立叶系数。该方法利用由下式给出的拉普拉斯算子、多项式、2阶微分算子:
L ≡ ∂ 2 ∂ 2 x + ∂ 2 ∂ 2 y L ≡ ∂ ∂ x ( ∂ ∂ x ) + ∂ ∂ y ( ∂ ∂ y )
因此当拉普拉斯作用于表面函数s(x,y)时,其结果是
Ls ( x , y ) = ∂ 2 s ( x , y ) ∂ 2 x + ∂ 2 s ( x , y ) ∂ 2 y Ls ( x , y ) = ∂ ∂ x ( ∂ s ( x , y ) ∂ x ) + ∂ ∂ y ( ∂ s ( x , y ) ∂ y )
使用上面给出的第二形式,并用(3)替换求和的第一积分中的并用(4)替换第二项中的
Figure A20048002099700336
求得表面的拉普拉斯算子为:
Ls ( x , y ) = ∂ ∂ x ( 1 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ iuS ( u , v ) e i ( ux + vy ) dudv ) + ∂ ∂ y ( 1 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ ivS ( u , v ) e i ( ux + vy ) dudv )
Ls ( x , y ) = 1 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ - u 2 S ( u , v ) e i ( ux + vy ) dudv + 1 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ - v 2 S ( u , v ) e i ( ux + vy ) dudv
Ls ( x , y ) = 1 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ - ( u 2 + v 2 ) S ( u , v ) e i ( ux + vy ) dudv - - - ( 9 )
等式(9)表示,二维函数的拉普拉斯算符的傅立叶系数等于函数本身的傅立叶系数的-(u2+v2)倍,使得
S ( u , v ) = - F ( Ls ( x , y ) ) ( u 2 + v 2 )
现在使拉普拉斯算符用上面限定的dx(x,y)和dy(x,y)表达,使得
Ls ( x , y ) = ∂ ∂ x ( dx ( x , y ) ) + ∂ ∂ y ( dy ( x , y ) )
并且,通过使用(5)和用于dy(x,y)的类似表达式
Ls ( x , y ) = ∂ ∂ x ( 1 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ Dx ( u , v ) e i ( ux + vy ) dudv ) + ∂ ∂ y ( 1 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ Dy ( u , v ) e i ( ux + vy ) dudv )
Ls ( x , y ) = 1 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ iuDx ( u , v ) e i ( ux + vy ) dudv + 1 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ ivDy ( u , v ) e i ( ux + vy ) dudv
Ls ( x , y ) = 1 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ i ( uDx ( u , v ) + vDy ( u , v ) ) e i ( ux + vy ) dudv - - - ( 10 )
(9)和(10)一定相等,并且将它们比较,可以看到要求:
-(u2+v2)S(u,v)=i(uDx(u,v)+vD(u,v))
S ( u , v ) = - i ( uDx ( u , v ) + vDy ( u , v ) ) u 2 + v 2 - - - ( 11 )
如上所述,通过对测量的梯度场分量进行傅立叶变换,求得Dx(u,v)和Dy(u,v)。然后在(11)中使用它们,以求得表面本身的傅立叶系数,反过来根据这些系数重建该表面。这种方法具有在重建中使用所有可获得的信息的效果,但Zernike多项式方法不能使用所有的可获得的信息。
应该注意,s(x,y)可被表达为新函数s(x,y)′和穿过原点的斜面表面的和。该和由下式给出:
s(x,y)=s(x,y)’+ax+by
然后f(x,y)关于x和y的第一偏导数由下式给出:
∂ s ( x , y ) ∂ x = ∂ s ( x , y ) ′ ∂ x + a
∂ s ( x , y ) ∂ y = ∂ s ( x , y ) ′ ∂ y + b
现在按照与导出(6)相同的过程,求得这些偏导数Dx(u,v)和Dy(u,v)的傅立叶变换为:
Dx ( u , v ) = 1 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ ∂ s ( x , y ) ∂ x e - i ( wx + vy ) dxdy = 1 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ ∂ s ( x , y ) ′ ∂ x e - i ( wx + vy ) dxdy + a 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ e - i ( wx + vy ) dxdy
Dx ( u , v ) = Dx ( u , v ) ′ + aδ ( u , v ) ) - - - ( 12 )
Dy ( u , v ) = 1 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ ∂ s ( x , y ) ∂ y e - i ( wx + vy ) dxdy = 1 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ ∂ s ( x , y ) ′ ∂ y e - i ( wx + vy ) dxdy + b 2 π ∫ - ∞ ∞ ∫ - ∞ ∞ e - i ( wx + vy ) dxdy
Dy ( u , v ) = Dy ( u , v ) ′ + bδ ( u , v ) ) - - - ( 13 )
在(12)和(13)中,δ(u,v)是迪拉克δ函数,该函数当u=v=0时取值为1且在其它情况下取值为0,在(11)中使用(12)和(13),表面的傅立叶变换的表达式可写为:
S ( u , v ) = - i ( uDx ( u , v ) ′ + vDy ( u , v ) ′ + ( ua + vb ) δ ( u , v ) ) u 2 + v 2
但应认识到,如果u和v不同时为零,则δ函数为零使得该项变为零,因此,上述等式中的项不会对S(u,v)的值产生任何影响。但是,在仅有的另一种情况下,u和v同时为零,这也会导致该项为零。这意味着重建的表面不是唯一的,而是一族表面中的一个,每一个表面具有不同的斜面(或直线)部分。因此,要重建由给定的梯度场表示的唯一表面,必须恢复正确的斜率。可以以几种方式进行斜率校正。
由于斜面梯度的分量a和b在表面上的每一点上相同,因此,如果能够在一点求得校正值,那么通过使初始重建的表面s(x,y)′增加表面(ax+by),可以立即将它们应用于任何位置。这可通过求得一点上的s(x,y)′的梯度并从给定的分量中减去这些分量而完成。差值是恒定的梯度值a和b。但是,当使用真实数据(real data)时,常存在噪声,因此最好用平均值求得a和b。完成这一点的一种有用的方法是,求得给定梯度场的平均梯度分量并使用这些平均值作为a和b。
<s/x>=a        <s/y>=b
重建的表面由下式给出:
s(x,y)=s(x,y)’+<s/x>x+<s/y>y         (14)
这里,通过使用上面推导的傅立叶重建方法求得s(x,y)′。
现在讨论使用离散快速傅立叶变换方法的实现本方法。二维离散快速傅立叶变换使用的算法由下式给出:
F ( k , l ) = &Sigma; m = 1 M &Sigma; n = 1 N f ( n , m ) e - i 2 &pi; ( ( k - 1 ) ( n - 1 ) N + ( l - 1 ) ( m - 1 ) M ) - - - ( 12 )
逆变换由下式给出:
F ( n , m ) = 1 MN &Sigma; k = 1 M &Sigma; i = 1 N F ( k , 1 ) e i 2 &pi; ( ( k - 1 ) ( n - 1 ) N + ( l - 1 ) ( m - 1 ) M ) - - - ( 13 )
当实现这些等式时,N和M一般被选择为相等。出于对计算速度的考虑,它们一般取2的幂。(12)和(13)假定在被间隔dx和dy分开的位置上对函数取样。为了算法简化起见,如下所述,dx和dy一般被设为相等。
在等式(12)中,使n为数组f(n,m)中的x数据的指数,k为变换数组F(k,l)中的变量u的指数。
我们从假定在离散的情况下x值相距dx,因此,当使用等式(12)时,每次增加n时,将x改变dx的量。因此,通过适当地选择坐标系,我们可以通过下式表示瞳孔数据的位置x:
x=(n-1)dx
使得,
(n-1)=x/dx
类似地,(m-1)可被设置为等于y/dy。通过利用这些关系式,(12)可被写为:
F ( k , 1 ) = &Sigma; m = 1 M &Sigma; n = 1 N f ( n , m ) e - i 2 &pi; ( ( k - 1 ) x Ndx + ( l - 1 ) y Mdy )
将该离散等式中的指数项与其积分形式(1)中的指数项相比较,可以看出,
u = 2 &pi; ( k - 1 ) Ndx 并且 v = 2 &pi; ( l - 1 ) Mdy
在这些等式中,注意到Ndx是取样区的x宽度,Mdy是取样区的y宽度。使得
(14)Ndx =X,总的x宽度
Mdy =Y,总的y宽度
以上等式变为
u ( k ) = 2 &pi; ( k - 1 ) X - - - ( 15 ) 并且 v ( 1 ) = 2 &pi; ( l - 1 ) Y
等式(15)允许从梯度分量dx(n,m)和dy(n,m)的离散快速傅立叶变换求得的傅立叶系数Dx(k,l)和Dy(k,l)转换为如下所示的表面的离散傅立叶系数S(k,l)。
S ( k , l ) = - i 2 &pi; ( k - 1 ) X Dx ( k , l ) - i 2 &pi; ( l - 1 ) Y Dy ( k , l ) ( 2 &pi; ( k - 1 ) X ) 2 + ( 2 &pi; ( l - 1 ) Y ) 2
如果上述N被选择等于M,dx被选择等于dy,使得X=Y,那么该等式可大大简化。从而,该等式变为
S ( k , l ) = ( - iX 2 &pi; ) ( k - 1 ) Dx ( k , l ) + ( l - 1 ) Dy ( k , l ) ( k - 1 ) 2 + ( l - 1 ) 2 - - - ( 16 )
我们现在考虑当k在1到N变化并且l在1到M变化时u(k)和v(l)在(13)中的取值。当k=l=1时,u=v=0并且指数取值1。当u和v增加1使得k=l=2时,u和v增加单位增量du和dv,使得
u ( 2 ) = 2 &pi; X = du 并且 v ( 2 ) = 2 &pi; Y = dv
k或l的各个增量分别由量du和dv表示,使得对于k或l的任何值,u(k)和v(l)可被写为:
u(k)=(k-1)du,v(1)=(1-1)dv
该过程可继续,直到k=N且l=M,此时
u ( N ) = 2 &pi; ( N - 1 ) X = 2 &pi; ( N ) X - 2 &pi; X = 2 &pi; ( N ) X - du
v ( M ) = 2 &pi; ( M - 1 ) X = 2 &pi; ( M ) X - 2 &pi; X = 2 &pi; ( M ) X - dv
现在考虑当保持这些条件时指数在(13)中的取值。以下,指数被表示为乘积
e i 2 &pi; ( N - 1 ) ( n - 1 ) N e i 2 &pi; ( M - 1 ) ( m - 1 ) M = e ( i 2 &pi; ( n - 1 c ) - i 2 &pi; ( n - 1 ) N ) e ( i 2 &pi; ( m - 1 c ) - i 2 &pi; ( m - 1 ) M ) = e - i 2 &pi; ( n - 1 ) N e - i 2 &pi; ( m - 1 ) M
使用与用于得到等式(15)所相同的逻辑,u(N)和v(M)的值为
u(N)=-du并且v(M)=-dv
利用该事实,对于k>1和l>1,现在可以在u(k+1)和u(N-k)之间以及v(l+1)和v(M-1)之间建立以下相关性
u(k)=-u(N-k+2)   v(l)=-v(M-1+2)
鉴于等式(15)
u ( N - k + 2 ) = - 2 &pi; ( k - 1 ) X - - - ( 17 ) v ( M - 1 + 2 ) = - 2 &pi; ( l - 1 ) Y
为了实现(15),首先注意,Dx(k,l)和Dy(k,l)形成为矩阵数组,因此最好形成作为矩阵数组的系数(k-1)和(l-1),使得可利用矩阵相乘方法以形成作为矩阵数组的S(k,l)。
假定Dx和Dy是N×N个单元的方阵列,使(k-1)阵列、K(k,l)形成为N×N阵列,其所有的行都相同,均包含从0(k=1)开始并以整数间隔1前进到顶盖值(value ceil)(N/2)-1的整数。“顶盖”算子将值N/2舍入到下一个更高的整数,用于考虑到N是奇数的情况。然后,鉴于(17)中给出的相互关系,下一行单元的值被赋予底值(value-floor)(N/2)。“底”算子将值N/2舍入到下一个更低的整数,也用于N是奇数的情况。具有底值(N/2)的行单元后面的行单元增加1,并继续下去直到到达最后的单元(k=N)并取值-1。这样,当矩阵|Dx|是倍数矩阵|K|项的被乘项时,具有相同的值k的|Dx|的各项乘以校正整数并由此乘以校正u值。
类似地,使矩阵|L(k,l)|形成为N×N阵列,其所有的列都相同,均包含从0(k=1)开始并以整数间隔1前进到顶盖值(N/2)-1的整数。然后,鉴于(17)中给出的相互关系,下一列单元的值被赋予底值(N/2)。具有底值(N/2)的列单元后面的列单元增加1,并继续下去直到到达最后的单元(l=N)并取值-1。这样,当矩阵|Dy|是倍数矩阵|L|项的被乘项时,具有相同的值l的|Dy|的各项乘以校正整数并由此乘以校正v值。
通过从通过逐项自乘|K|倍和自乘|L|倍形成的矩阵和产生矩阵|D|的(15)的分母。|D|的(1,1)单元总是为零,并且,为了避免被零除的问题,当|D|初始地形成后将其设为1。由于此时|K|和|L|的(1,1)单元也为零,因此这具有将|S|的(1,1)单元设为等于零的效果。这反过来意味着重建表面的平均高度为零,这一点可这样理解:当k=l=1时(12)的值是函数f(x,y)的所有值的求和。如果该求和为零,那么函数的平均值也为零。
使两个矩阵|A|和|B|的逐项乘由|A|.*|B|表示,使|A|被|B|的逐项除由|A|./|B|表示。那么,以矩阵的形式(16)可被写为:
| S | = ( - iX 2 &pi; ) ( | K | . * | Dx | + | L | . * Dy ( k , l ) ) / ( | K | . * | K | + | L | . * | L | ) - - - ( 18 )
公因子
Figure A20048002099700392
既不是位置的函数,也不是“频率”(傅立叶变换空间的变量)的函数。因此它是球形比例因子。
作为当将(18)编码时的实际问题,如果通过使用将u=v=0单元放置在位置(floor(N/2)+1,floor(N/2)+1)上的标准离散傅立叶变换象限移位技术第一次“移位”变换矩阵Dx和Dy,那么形成K和L很简单。从而,K的各行和L的各列可由下式形成:
row=[1,2,3,...N-2,N-1,N]-(floor(N/2)+1)
                  column=rowT
在通过使用移位矩阵用(18)求得矩阵|S|后,通过使用逆离散反傅立叶变换(13)在求得s(x,y)的各值之前反向移位|S|。
最后的步骤是求得梯度场dx(n,m)和dy(n,m)的平均值。这些平均值乘以评估的每一个表面点的各个x和y值,并被加到在以上的步骤中求得的s(x,y)的值上,以给出完全重建的表面。
实验结果
现在描述比较展开级数(例如,Zernike多项式)重建方法、直接积分重建方法和傅立叶变换重建方法的表面重建的一些试验方法的详细描述。
虽然这里没有详述,但应理解,本发明还包含使用用于重建波阵面立视图的直接积分算法和模块。傅立叶变换模块、直接积分模块和Zernike模块的使用不矛盾或不相互排斥,如果需要,它们可以被组合。例如,除了包含示例的模块或作为示例模块的替换,图5的模块还可以包含直接积分模块。在发明名称同为“Direct Wavefront-BasedCorneal Ablation Treatment Program”的在2001年12月6日提交的共同受让的美国专利申请No.10/006992和2001年12月6日提交的PCT申请No.PCT/US01/46573中更完整地描述了直接积分模块和方法,在此引入它们的全部公开内容作为参考。
要比较各种方法,在塑料上烧蚀表面,并将各种重建方法与直接表面测量相比,以确定各种方法的精确度。试验产生下述三种不同的试验表面:
(1)6mm瞳孔上的+2烧蚀,其中烧蚀中心相对于瞳孔中心偏离约1mm;
(2)具有1.5μm高、偏心1.0mm的直径为2.5mm的“凸点”的远视眼形状I。
(3)具有1.0μm高、偏心0.5mm的直径为2.0mm的“凸点”的远视眼形状II。
烧蚀的表面由波阵面传感器系统30成像(参见图3和图3A),并且处理Hartmann-Shack点图以得到波阵面梯度。烧蚀的表面还通过Phase Shift Technologies制造的表面测绘干涉仪Micro XCAM扫描,以产生高精度的表面立视图。将由Micro XCAM直接测量的立视图与由各种不同的算法重建的立视图相比。具有最小的均方根(RMS)误差的算法被认为是在重建表面中最为有效。
在直接测量和数学重建中,都存在需要校正的系统“倾斜”。对于直接测量,通过减去适合表面的平面从数据中消除表面的倾斜(由于保持样品的样品台的倾斜而被引入)。
对于数学重建,相对于波阵面测量系统中透镜组阵列的表面的角位置和空间位置在重建表面中引入了倾斜和偏心。通过识别诸如顶部的顶点的主要特征并解释整个表面数据以与重建中的该特征的位置匹配,完成校正“偏心”对准。
为了消除倾斜,在一个实施例中,将沿x轴和y轴的重建表面的线剖面与测量的表面的相应剖面相比较。估计重建表面相对于测量表面的斜率。并且,还确定用于中心对准的相同的主要特征(例如,顶部)的高度差。从重建的表面减去由该斜率和高度差限定的表面。在另一个实施例中,确定了傅立叶变换算法中的倾斜可来自在重建过程中设置为零的x和y梯度的傅立叶变换的DC分量。结果,整个波阵面的净梯度缺失。在平均梯度场“复原(untip)”中加入重建表面。可以理解,这些方法可以包含在本发明的模块中以从重建中消除倾斜。
图6中示例了对于偏心+2透镜的重建表面和直接测量表面的比较。从图6所示例,所有重建方法都与表面匹配得很好。重建中的RMS误差如下:
  傅立叶   0.2113e-3
  直接积分   0.2912e-3
  Zernike(6阶)   0.2264e-3
图7示出了远视眼形状I重建的断面。可以看出,Zernike 6阶重建过度地加宽凸点特征。其它重建提供与表面的更好的匹配。四种重建方法的RMS误差如下:
  傅立叶   0.1921e-3
  直接积分   0.1849e-3
  Zernike(6阶)   0.3566e-3
  Zernike(10阶)   0.3046e-3
图8示出了远视眼形状II重建的断面。数据定性类似于图7。四种重建方法的RMS误差如下:
  傅立叶   0.1079e-3
  直接积分   0.1428e-3
  Zernike(6阶)   0.1836e-3
  Zernike(10阶)   0.1413e-3
从以上结果可以看出,对于具有大于约1-2毫米的特征的平滑表面,6阶Zernike重建足够精确。但是,对于诸如远视眼形状的凸点的锐变特征,与其它重建方法相比,6阶Zernike重建给出与实际表面的较差的匹配。
在一些病变的眼睛和手术处理的眼睛中,可能存在角膜表面的锐变特征或曲率的局部快速变化。另外,具有小的和锐利特征的处理可被应用于远视眼和一些高度畸形的眼睛。
申请人相信,傅立叶变换提供更好的结果的部分原因在于,与Zernike重建算法(被限定在圆上,并将瞳孔近似为圆)不同,傅立叶变换算法(以及直接积分算法)利用所有可获得的数据并允许基于瞳孔的实际形状(典型为略椭圆状)进行计算。离散傅立叶分析的带宽是波阵面测量设备的取样频率的一半。因此,傅立叶方法可利用所有的梯度场数据点。并且,由于傅立叶变换算法固有地具有截止频率,因此傅立叶算法滤除(即,设置为零)所有比由数据取样间隔表示的频率更高的频率,从而防止诸如混叠的人工因素引入重建。最后,由于许多波阵面测量系统在正方形网格上对波阵面表面进行取样,并且在正方形网格上执行傅立叶方法,因此,傅立叶方法更好地适于来自波阵面设备的输入数据。
相反,Zernike方法使用径向和角度项(例如极坐标),因此,Zernike方法不平等地加重中心点和周边点。当使用更高阶的多项式再现波阵面中的一些细节时,作为半径函数的幅度的振荡不均匀。另外,对于任何给定的多项式,非零的纵向指数值的纵向项是正弦函数。由该Zernike项引入的峰值和谷值更大,较远的一个远离瞳孔的中心。并且,它还引入波阵面的非均匀空间频率取样。因此,瞳孔的中心与周边相比,相同的多项式项可适应波阵面中非常较小的变化。为了在瞳孔边缘得到局部变化的良好取样,必须使用较大数量的Zernike项。不幸的是,这些较大数量的Zernike项会导致在瞳孔中心过取样以及引入诸如混叠的人工因素。由于傅立叶方法提供均匀的空间取样,因此可以避免引入这种人工因素。
图9到11中示出了临床数据的其它试验结果。比较重建波阵面的傅立叶方法与根据临床数据重建波阵面的6阶Zernike方法和直接积分方法。然后,对重建的波阵面求微分,以计算梯度场。计算和测量的梯度场的均方根(RMS)被用作重建质量的度量。
重建的试验方法如下:通过使用三种算法(例如,Zernike、傅立叶和直接积分)重建与具有大量像差的眼睛相对应的波阵面。计算中使用的瞳孔尺寸是3mm的半径。比较重建的波阵面的梯度场与测量的梯度场。每一个取样点的梯度的x和y分量被平方并相加起来。总和的平方根提供有关表面的曲率的信息。这样的数值等于乘以取样点总数的梯度的平均幅度。例如,0的量对应于平整或平面波阵面。具有该量的梯度场的RMS偏离的比值给出重建质量的度量。例如,比值越小,重建的波阵面越接近于直接测量的波阵面。具有不同重建的量的RMS偏离(上述)的比值如下:
  Zernike(6阶)   1.09
  Zernike(10阶)   0.82
  直接积分   0.74
  傅立叶   0.67
图9示例了计算和测量的梯度场之差的矢量图。Zernike图(由“Z场”表示)是用于使用最多10阶的项重建。图11示例了使用最多6阶的项的Zernike重建算法不能正确地再现波阵面上小的和锐变的特征。如图10所示,使用最多10阶的项的Zernike算法可更好地再现小的和锐变的特征。从结果可以看出,对于傅立叶方法,与测量的梯度之间的RMS偏离最小。
现在描述根据本发明实施例的使用迭代傅立叶变换算法确定光学表面模型的数学推导。
在波阵面技术中,可以测量诸如人眼的光学系统的光程差(OPD)。在波阵面检测方面有不同的技术,并且,Hartmann-Shack波阵面检测已变成最流行的用于测量眼睛像差的技术。Hartmann-Shack设备通常将诸如瞳孔的孔径分成一组子孔径;每一个子孔径对应于从透镜组阵列投影的一个区域。由于Hartmann-Shack设备测量每一个子孔径的局部斜率(或梯度),因此希望使用局部斜率数据用于波阵面重建。
假定W(x,y)是要重建的波阵面,那么W(x,y)沿x轴的局部斜率是
Figure A20048002099700441
沿y轴的局部斜率是
Figure A20048002099700442
进一步假定c(u,z)是W(x,y)的傅立叶变换,那么W(x,y)将是c(u,v)的逆傅立叶变换。因此,我们得到
W(x,y)=∫∫c(u,v)exp[i2π(ux+vy)]dudv,     (19)
这里,c(u,v)是展开系数。在等式(19)中分别取x和y的偏导数,我们得到
&PartialD; W ( x , y ) &PartialD; x = i 2 &pi; &Integral; &Integral; uc ( u , v ) exp [ i 2 &pi; ( ux + vy ) ] dudv &PartialD; W ( x , y ) &PartialD; y = i 2 &pi; &Integral; &Integral; vc ( u , v ) exp [ i 2 &pi; ( ux + vy ) ] dudv - - - ( 20 )
cu表示W(x,y)的x导数的傅立叶变换,cv表示W(x,y)的y导数的傅立叶变换。根据傅立叶变换的定义,我们得到
c u ( u , v ) = &Integral; &Integral; &PartialD; W ( x , y ) &PartialD; x exp [ - i 2 &pi; ( ux + vy ) ] dxdy c v ( u , v ) = &Integral; &Integral; &PartialD; W ( x , y ) &PartialD; y exp [ - i 2 &pi; ( ux + vy ) ] dxdy - - - ( 21 )
可以以逆傅立叶变换的形式将等式(21)写为
&PartialD; W ( x , y ) &PartialD; x = &Integral; &Integral; c u ( u , v ) exp [ i 2 &pi; ( ux + vy ) ] dudv &PartialD; W ( x , y ) &PartialD; y = &Integral; &Integral; c v ( u , v ) exp [ i 2 &pi; ( ux + vy ) ] dudv - - - ( 22 )
比较等式(20)和(22),我们得到
cu(u,v)=i2πuc(u,v)                      (23)
cv(u,v)=i2πvc(u,v)                      (24)
如果我们在等式(23)的两个边同时乘以u,并在等式(24)的两边同时乘以v,并将它们相加,我们得到
ucu(u,v)+vcv(u,v)=i2π(u2+v2)c(u,v).    (25)
从等式(25),我们得到傅立叶展开系数如下:
c ( u , v ) = - i u c u ( u , v ) + v c v ( u , v ) 2 &pi; ( u 2 + v 2 ) - - - ( 26 )
因此,可得到波阵面的傅立叶变换:
c ( u , v ) = - i 2 &pi; ( u 2 + v 2 ) [ u &Integral; &Integral; &PartialD; W ( x , y ) &PartialD; x exp [ - i 2 &pi; ( ux + vy ) ] + v &Integral; &Integral; &PartialD; W ( x , y ) &PartialD; y exp [ - i 2 &pi; ( ux + vy ) ] ] - - - ( 27 )
由此,取等式(27)中的逆傅立叶变换,我们得到波阵面如下:
W(x,y)=∫∫c(u,v)exp[i2π(ux+vy]dudv.    (28)
等式(28)是波阵面重建的最终解。也就是说,如果我们知道波阵面斜率数据,我们可以通过使用等式(27)计算傅立叶级数的系数。使用等式(28),然后可以重建未知的波阵面。在Hartmann-Shack方法中,测量一组局部波阵面斜率,因此,这种方法可应用于等式(27)。
但是,在一些情况下,上面的波阵面重建方法可能限制于无界函数。为了利用边界条件(即孔径边界)得到波阵面估算,申请人已发现,迭代重建方法是有用的。首先,可遵循上述方法提供初始解,该初始解将函数值赋予比函数边界大的方格。这类似于如下面讨论的那样将数据点设为小的非零值。然后可以计算整个方格的估算波阵面的局部斜率。下一步,所有已知的局部斜率值,即,从Hartmann-Shack设备测量的梯度,可以重写计算的斜率。基上更新的斜率,可以重新应用上述方法,并可得到波阵面的新的估算。重复该程序,直到达到预定次数的迭代或满足预定的准则。
在WaveScan软件中,在实现傅立叶重建中使用三种主要的算法。这些算法是用于实现整个迭代傅立叶重建的基础。第一种算法是迭代傅立叶重建本身。第二种算法用于计算WaveScan设备中的显示器的折射。第三种算法用于报告均方根(RMS)误差。
A.波阵面表面重建
图12示例了示意性的迭代方法。该方法从输入来自眼睛的光学组织系统的光学数据开始。通常,光学数据是由波阵面测量设备产生的波阵面数据,并作为测量的梯度场200输入,这里,测量的梯度场对应于孔径内的一组局部梯度。然后将迭代傅立叶变换应用于该光学数据,以确定光学表面模型。该方法建立第一组合的梯度场210,该梯度场210包括设置在第一外部梯度场内的测量的梯度场200。第一外部梯度场可对应于在平面上具有恒定值W(x,y)并可与任何孔径一起被使用的平面波或无界函数。
在一些情况下,测量的梯度场200可包含丢失、错误或不充分的数据。在这些情况下,当建立组合的梯度场210时,有可能忽视这些数据点,并只使用那些认为是好的测量梯度场200的值。测量梯度场200中要忽视的点被指定与第一外部梯度场相对应的值。通过应用傅立叶变换,第一组合梯度场200用于导出第一重建波阵面220,然后,使用该第一重建波阵面200提供第一修正梯度场230。
建立包括设置在第一修正梯度场230内的测量梯度场200的第二组合梯度场240。本质上,第二外部梯度场是未被测量梯度场200所代替的第一修正梯度场230的部分。当建立第二组合梯度场240时,有可能以与上述方式相同的方式,只使用那些认为是对的测量梯度场200的值。通过应用傅立叶变换,第二组合梯度场240用于导出第二重建波阵面250。第二重建波阵面250或至少其一部分,可被用于提供最终的重建波阵面290。然后,可基于该最终的重建波阵面290确定光学表面模型。
可选地,可进一步迭代第二组合梯度场。例如,第二重建波阵面250可用于提供第(n-1)个梯度场260。然后,可建立包括被设置在第(n-1)个修正梯度场260内的测量梯度场200的第(n)个组合梯度场270。本质上,第(n)个外部梯度场是未被测量的梯度场200所代替的第(n-1)个修正梯度场260的部分。通过应用傅立叶变换,第(n)个组合梯度场270用于导出第(n)个重建波阵面280。第(n)个重建波阵面280或至少其一部分,可被用于提供最终的重建波阵面290。然后,可基于该最终的重建波阵面290确定光学表面模型。实际上,每一次迭代都可使每一连续的重建波阵面更接近实际,特别是对于瞳孔或孔径的边界或周边。
假定Hartmann-Shack测量由dZx和dZy表示的局部波阵面斜率,这里dZx代表沿x方向的波阵面斜率,dZy代表沿y方向的波阵面斜率。在计算波阵面估算值时,可以使用两个临时数组cx和cy以存储估算的波阵面w的局部斜率。它还帮助实现诸如FFT、iFFT、FFTShift和iFFTShift的标准函数。
下面描述示意性的算法:
1.将数据点设为小的但非零的值,这里,在测量(从Hartmann-Shack设备)中没有数据表示(标记(mark)=1.2735916e-99)。
2.开始10次迭代的迭代重建
a.对于梯度场不等于标记的原始数据点,将梯度场dZx和dZy复制到梯度场数组cx和cy
b.分别计算cx和cy的快速傅立叶变换(FFT)
c.对在步骤b中得到的数组进行象限交换(FFTShift)
d.根据等式(26)计算c(u,v)
e.对在步骤d中得到的数组进行象限交换(iFFTShift)
f.对在步骤e中得到的数组进行逆傅立叶变换(iFFT)
g.计算更新的表面估算值w(来自步骤e的数组的实部)
h.从w计算更新的梯度(w对x和y的导数)
i.当迭代的次数等于10次时,结束
3.通过使用来自步骤2.h的估算值,计算平均梯度
4.从从步骤2.h得到的梯度场减去平均梯度,以消除斜坡/倾斜分量
5.应用步骤2.b-g以得到波阵面的最终估算
B.波阵面折射度(refraction)计算
当构造波阵面时,波阵面折射度的计算可能比使用Zernike重建时更加困难。原因是,一旦用Zernike重建得到Zernike系数,那么可以用下公式直接计算波阵面折射度:
C = - 4 6 ( c 2 - 2 ) 2 + ( c 2 2 ) 2 R 2 , - - - ( 29 )
S = - 4 3 c 2 0 R 2 - 0.5 C , - - - ( 30 )
&theta; = 1 2 tan - 1 [ c 2 2 c 2 - 2 ] . - - - ( 31 )
这里,c2 -2、c2 0和c2 2代表第二阶中的三个Zernike系数,S代表球面,C代表圆柱面,θ代表圆柱轴。但是,使用傅立叶重建,傅立叶系数的任何一个都与经典像差无关。由此,为了使用等式(29)~(31)计算折射度,需要Zernike分解以得到Zernike系数。
Zernike分解尝试用一组Zernike多项式函数拟合表面,该Zernike多项式函数在拟合被最小化后具有最小平方的含义即均方根(RMS)误差。为了实现最小平方准则,可以使用奇异值分解(SVD),该奇异值分解是基于最小平方准则的迭代算法。
假定波阵面由下式的Zernike展开表示:
W ( r , &theta; ) = &Sigma; i = 0 N c i Z i ( r , &theta; ) - - - ( 32 )
或以矩阵形式数字化表示为:
W=Z·c,                        (33)
这里,W是波阵面图的2维MxM矩阵,Z是矩阵的具有N个层的MxMxN张量,每一个代表具有单位系数的特定Zernike模的一个表面,并且c是包含Zernike系数组的列向量。
假定已知W求c,如果我们得到所谓通用的Z的逆矩阵,那么很简单,
c=Z+·W.                        (34)
奇异值分解(SVD)算法可以以最小平方的含义计算任何矩阵的通用的逆矩阵。因此,如果我们得到
Z=U·w·VT,                    (35)
那么该组系数的最终解将是
c=V·w-1·UT·W.                (36)
SVD中的一个考虑是确定截止本征值。在以上等式中,w是对角矩阵,它具有在对角线中本征值从最大到最小排列的单元。但是,在很多情况下,最小本征值过于接近零,以至于该值的倒数太大,因此它会放大输入表面矩阵中的噪声。为了防止矩阵求逆的问题,希望将条件数r限定为最大本征值与截止本征值的比值。任何小于截止本征值的本征值都将不被用于求逆或简单地将零作为倒数。在一个实施例中,可以使用100到1000的条件数。在另一实施例中,可以使用200的条件数。
一旦实现Zernike分解,那么可以通过使用等式(29)~(31)得到球面、圆柱以及圆柱轴的计算。但是,一般在与测量面不同的顶点距离给出折射度。假定d代表顶点距离,那么有可能使用以下公式计算新的折射度(圆柱轴不变):
S &prime; = S 1 + dS - - - ( 19 )
C &prime; = S + C 1 + d ( S + C ) - S &prime;
该算法可如下描述:
1.将球面和圆柱的预补偿添加到由迭代傅立叶重建算法估算的波阵面
2.分解来自步骤1的表面,以得到第一组五个Zernike系数
3.应用等式(29)~(31)计算折射度
4.使用等式(37)将折射度重新调整到顶点距离
5.根据柱面符号显示折射度
C.波阵面均方根(RMS)计算
最后,可以计算波阵面均方根(RMS)误差。重新使用Zernike重建直接计算RMS误差。但是,使用迭代傅立叶重建,会如上所讨论的那样更加困难。在这种情况下,需要Zernike分解,以计算波阵面折射度,并因此可得到Zernike分解用于计算RMS误差。
要得到RMS误差,可以使用三种不同的类别:低阶RMS、高阶RMS和总RMS。对于低阶RMS,可以使用以下公式:
lo . r . m . s . = c 3 2 + c 4 2 + c 5 2 - - - ( 38 )
这里,c3、c4和c5分别是散光、散焦和散光的Zernike系数。对于高阶RMS,可以用下式使用整个波阵面:
ho . r . m . s . = &Sigma; n ( v i - v &OverBar; ) 2 n - - - ( 39 )
这里,vi代表第i个位置上的波阵面表面值,v代表瞳孔内的平均波阵面表面值,n代表瞳孔内的位置的总数。要计算总RMS,可以使用下式:
r . m . s . = lo . r . m . s . 2 + ho . r . m . s . 2 - - - ( 40 )
该算法是:
1.对于低阶RMS,使用等式(38)
2.对于高阶RMS,使用等式(39)
3.对于总RMS,使用等式(40)
收敛性
可以使用收敛性,以估算在迭代傅立叶变换算法中需要迭代的次数。如上所述,迭代傅立叶重建算法对于无界函数有效。但是,在上述实施例中,由于瞳孔函数被用作边界,因此等式(27)和(28)不能提供适当的解。但使用根据本发明的迭代算法,可以得到相当好的边界函数的解。表1表示一些Zernike模的重建后的均方根(RMS)值,每一个具有一微米RMS输入。
表1从重建的波阵面得到的RMS值
 #迭代   Z3   Z4   Z5   Z6   Z7   Z10   Z12   Z17   Z24   实际
  1   0.211   0.986   0.284   0.247   1.772   0.236   0.969   1.995   0.938   0.828
  2   0.490   0.986   0.595   0.538   1.353   0.518   0.969   1.522   0.938   0.891
  5   0.876   0.986   0.911   0.877   1.030   0.861   0.969   1.069   0.938   0.966
  10   0.967   0.986   0.956   0.943   0.987   0.935   0.969   0.982   0.938   0.979
  20   0.981   0.986   0.962   0.955   0.982   0.951   0.969   0.968   0.938   0.981
  50   0.987   0.986   0.966   0.963   0.980   0.960   0.969   0.963   0.938   0.981
其中实际是用于具有总量为1微米误差的组合Zernike模的波阵面。
作为例子,图13表示分别使用一次、二次、五次和十次迭代的迭代傅立叶技术的、散光项(Z3)的波阵面重建的表面图。为了更逼真,图14表示分别使用一次、二次、五次和十次迭代的迭代傅立叶技术的、实际眼睛的波阵面重建的表面图,它表示它比单个非对称Zernike项更快地收敛。很显然,10次迭代看起来实现了大于90%恢复具有Zernike输入的输入RMS误差,但是,除非在眼睛中存在纯圆柱,否则5次迭代可以是充分的。
外插
迭代傅立叶变换方法和系统可解释丢失、错误或其它不充分的数据点。例如,在一些情况下,测量的梯度场200可能包含有缺陷的数据。在这些情况下,当建立组合的梯度场210时,可以忽视这些数据点,并只使用被认为是好的测量梯度场200的值。
开发了称为WaveTool的研究软件程序供研究使用。该软件以C++写成,其中仔细测试了迭代傅立叶重建的实现并将结果与用Matlab代码得到的结果相比较。在试验过程中,假定顶行、底行或顶行和底行都为缺失数据,使得傅立叶重建必须在波阵面重建的期间估算梯度场。在另一种情况下,假定其中的一个中间模式由于角膜反射而缺失、模拟数据缺失。绘制具有和不具有预补偿的重建波阵面,以表示这种变化。同时,比较均方根(RMS)误差以及折射度。用10次迭代重建每个波阵面。
在计算中仅使用一只眼睛。原始的H-S图像由最大可计算的6mm瞳孔的15×15阵列梯度场构成。当数据缺失时,外插法对于计算存在缺失数据时的6mm瞳孔的波阵面是有用的。表2示出了在一些缺失数据的情况下折射度、总RMS误差以及表面RMS误差(与没有缺失数据的情况相比)的变化。
由于照相机典型为矩形,因此,测量的梯度场沿垂直方向可能具有缺失的边缘。经常是,捕获沿水平方向的所有数据,但沿垂直方向可能存在缺失数据。在这些情况下,测量的梯度场会具有缺失的顶行、底行或顶行和底行。
表2:具有缺失数据的重建波阵面的折射度和RMS的比较
  情况   Rx   总RMS   RMS误差
  没有缺失数据   -2.33DS/-1.02DC×170°@12.5   3.77μm   -
  缺失顶行   -2.33DS/-1.03DC×170°@12.5   3.78μm   0.0271μm
  缺失底行   -2.37DS/-0.97DC×169°@12.5   3.75μm   0.0797μm
  同时缺失顶行和底行   -2.37DS/-0.99DC×170°@12.5   3.76μm   0.0874μm
  缺失一点   -2.33DS/-1.02DC×170°@12.5   3.77μm   0.0027μm
  缺失四点   -2.32DS/-1.03DC×170°@12.5   3.76μm   0.0074μm
图15表示对于不同的缺失数据的情况具有和不具有预补偿的重建波阵面。顶行表示具有预补偿的波阵面,底行表示不具有预补偿的波阵面。以下情况被如下示例:(a)没有缺失数据;(b)缺失顶行;(c)缺失底行;(d)同时缺失顶行和底行;(e)缺失中间点;(f)缺失四个点。结果表明,缺失少量的数据影响很小,并且算法可以重建相当精确的波阵面。
采用10次迭代,迭代傅立叶重建可以提供与输入数据相比大于90%的精度。在由于角膜反射而缺失瞳孔内的测量数据或缺失CCD检测器外的测量数据的情况下,这种方法也是有益的。
在本发明的范围内,可以进行各种修改。例如,本发明基于傅立叶的方法可用于上述烧蚀监视系统反馈系统中,用于在各激光脉冲的过程中和/或各激光脉冲之间对患者的眼睛进行实时外科内(intrasurgical)测量。由于较高的测量精度和较高的速度,因此基于傅立叶的方法尤其适合于这些用途。可以在示意性的方法步骤或系统模块中包含各种参数、变量、因素等。虽然为了便于理解通过例子详细描述了具体实施例,但各种改进、变化和修改对于本领域的普通技术人员而言是显而易见的。虽然具体参照使用透镜组的波阵面系统描述了本发明,但也可使用测量穿过眼睛的光角度的其它适当的波阵面系统。例如,可以在本发明中使用使用光线跟踪像差计、切尔宁(tscherning)像差计和动态测眼膜术的原理的系统。上述各系统可分别从德克萨期州Bellaire市的TRACEY Technologies公司、德国埃尔兰根市的Wavelight公司和加利福尼亚州弗里蒙特市的Nidek公司得到。也可以用在美国专利No.6099125、6000800和5258791中描述的空间分辨折射仪实践本发明,在此引入它们的全部公开内容作为参考。除了折射激光角膜手术外,受益于本发明的处理还包括眼内透镜、隐形眼镜、眼镜和其它外科方法。因此,本发明的范围只由所附的权利要求书限定。

Claims (56)

1.一种确定眼睛的光学组织系统的光学表面模型的方法,该方法包括以下步骤:
a)从眼睛的光学组织系统输入光学数据,所述光学数据包括一组局部梯度;
b)基于该组局部梯度建立第一组合梯度场;
c)通过应用傅立叶变换于所述第一组合梯度场导出第一重建波阵面;
d)基于该第一重建波阵面提供第一修正梯度场;
e)基于该第一修正梯度场建立第二组合梯度场;
f)通过应用傅立叶变换于所述第二组合梯度场导出第二重建波阵面;和
g)基于该第二重建波阵面确定所述光学表面模型。
2.根据权利要求1所述的方法,其中所述光学数据包含波阵面数据。
3.根据权利要求2所述的方法,其中所述波阵面数据包含孔径内的一组局部梯度。
4.根据权利要求3所述的方法,其中所述孔径包含眼睛的瞳孔。
5.根据权利要求1所述的方法,其中:
所述第一组合梯度场包括第一外部梯度场和测量梯度场,使得所述测量梯度场设置在所述第一外部梯度场的内部,并且所述测量梯度场对应于所述组局部梯度;以及
所述第二组合梯度场包括第二外部梯度场和测量梯度场,使得所述第二外部梯度场对应于所述第一修正梯度场,并且所述测量梯度场设置在所述第二外部梯度场的内部。
6.根据权利要求5所述的方法,还包括以下步骤:
选择所述第二重建波阵面的至少一部分,以提供最终的重建波阵面,并基于该最终的重建波阵面确定所述光学表面模型。
7.根据权利要求1所述的方法,还包括以下步骤:
反复执行步骤(d)到(f),以通过应用傅立叶变换于第n个组合梯度场导出第n个重建波阵面;选择所述第n个重建波阵面的至少一部分,以提供最终的重建波阵面;和基于该最终的重建波阵面确定所述光学表面模型。
8.一种绘制眼睛的波阵面误差图的方法,该方法包括以下步骤:
根据权利要求1所述的方法确定光学表面模型;和
基于该光学表面模型绘制所述眼睛的波阵面误差图。
9.一种为眼睛的光学组织系统计算校正烧蚀模式的方法,该方法包括以下步骤:
根据权利要求1所述的方法确定光学表面模型;和
基于该光学表面模型计算所述校正烧蚀模式。
10.根据权利要求9所述的方法,其特征在于,
所述第一组合梯度场包括第一外部梯度场和测量梯度场,使得所述测量梯度场设置在所述第一外部梯度场的内部,并且所述测量梯度场对应于所述组局部梯度;以及
所述第二组合梯度场包括第二外部梯度场和测量梯度场,使得所述第二外部梯度场对应于所述第一修正梯度场,并且所述测量梯度场设置在所述第二外部梯度场的内部。
11.根据权利要求10所述的方法,还包括以下步骤:
选择所述第二重建波阵面的至少一部分,以提供最终的重建波阵面,并基于该最终的重建波阵面确定所述光学表面模型。
12.根据权利要求11所述的方法,还包括以下步骤:
反复执行步骤(d)到(f),以通过应用傅立叶变换于第n个组合梯度场导出第n个重建波阵面;选择所述第n个重建波阵面的至少一部分,以提供最终的重建波阵面;和基于该最终的重建波阵面确定所述光学表面模型。
13.一种修改眼睛的光学组织表面的方法,该方法包括以下步骤:
根据权利要求9所述的方法计算校正烧蚀模式;和
通过将激光烧蚀应用到眼睛上根据所述校正烧蚀模式修改所述光学组织表面。
14.一种用于确定眼睛的光学组织系统的光学表面模型的系统,该系统包括:
a)用于通过所述光学组织系统传输图像的光源;
b)被定向用于通过检测该传输图像为所述光学组织系统确定一组局部梯度的传感器;
c)与所述传感器耦合的处理器;和
d)与所述处理器耦合的存储器,所述存储器被配置为存储多个由所述处理器执行的代码模块,所述多个代码模块包括:
i)用于基于所述组局部梯度建立第一组合梯度场的模块;
ii)用于通过应用傅立叶变换于该第一组合梯度场导出第一重建波阵面的模块;
iii)用于基于该第一重建波阵面提供第一修正梯度场的模块;
iv)用于基于该第一修正梯度场建立第二组合梯度场的模块;
v)用于通过应用傅立叶变换于所述第二组合梯度场导出第二重建波阵面的模块;和
vi)用于基于该第二重建波阵面为所述眼睛的光学组织系统确定光学表面模型的模块。
15.一种绘制眼睛的光学组织系统的波阵面误差图的系统,所述系统包括权利要求14所述的系统,其中所述多个代码模块还包括基于所述光学表面模型绘制所述眼睛的波阵面误差图的模块。
16.一种为眼睛的光学组织系统计算校正烧蚀模式的系统,所述系统包括权利要求14所述的系统,其中所述多个代码模块还包括基于所述光学表面模型为所述眼睛计算校正烧蚀模式的模块。
17.一种用于修改眼睛的光学组织表面的系统,所述系统包括权利要求16所述的系统,其中所述多个代码模块还包括基于所述校正烧蚀模式修改所述眼睛的光学组织表面的激光器。
18.一种用于确定与眼睛的光学组织系统相对应的光学表面模型的系统,该系统包括:
a)用于通过所述光学组织系统传输图像的装置;
b)位于来自所述图像的光程中的、基于所传输的图像为所述光学组织系统确定一组局部梯度的装置;
c)与该局部梯度确定装置耦合的、基于该组局部梯度建立第一组合梯度场的装置;
d)与该第一组合梯度场建立装置耦合的、基于应用傅立叶变换于所述第一组合梯度场导出第一重建波阵面的装置;
e)与该第一重建波阵面导出装置耦合的、基于所述第一重建波阵面提供第一修正梯度场的装置;
f)与该第一修正梯度场提供装置耦合的、基于所述第一修正梯度场建立第二组合梯度场的装置;
g)与该第二组合梯度场建立装置耦合的、基于应用傅立叶变换于所述第二组合梯度场导出第二重建波阵面的装置;和
h)与该第二重建波阵面导出装置耦合的、基于所述第二重建波阵面为所述眼睛的光学组织系统确定光学表面模型的装置。
19.一种存储在计算机可读存储介质上的计算机程序,所述计算机程序包括:
a)用于通过眼睛的光学组织系统传输图像的代码;
b)用于基于被传输的图像为所述眼睛的光学组织系统确定一组局部梯度的代码;
c)用于基于所述一组局部梯度建立第一组合梯度场的代码;
d)用于基于应用傅立叶变换于所述第一组合梯度场导出第一重建波阵面的代码;
e)用于基于该第一重建波阵面提供第一修正梯度场的代码;
f)用于基于该第一修正梯度场建立第二组合梯度场的代码;
g)用于基于应用傅立叶变换于所述第二组合梯度场导出第二重建波阵面的代码;和
h)用于基于该第二重建波阵面为所述眼睛的光学组织系统确定光学表面模型的代码。
20.根据权利要求19所述的计算机程序,还包括:
用于基于所述光学表面模型计算校正烧蚀模式的代码。
21.根据权利要求20所述的计算机程序,还包括:
用于基于所述校正烧蚀模式向所述眼睛传送激光能量的代码。
22.一种重建眼睛的光学组织的方法,所述方法包括以下步骤:
通过所述眼睛的光学组织传输图像;
从通过所述眼睛的光学组织的所传输的图像测量表面梯度;和
对所述表面梯度应用傅立叶变换算法,以重建与所述眼睛的光学组织相对应的表面。
23.根据权利要求22所述的方法,包括以下步骤:
在测量所述表面梯度的过程中,将所述眼睛的光学组织的重建表面与已得到的所述眼睛的图像对准。
24.根据权利要求22或23所述的方法,包括以下步骤:
基于如由所述傅立叶重建表面表示的所述眼睛的光学组织计算校正烧蚀模式。
25.根据权利要求24所述的方法,其中计算校正烧蚀模式的步骤包括导出所述光学组织的高度的所提出的变化,以便实现所述眼睛的光学性能的想要的变化。
26.根据权利要求25所述的方法,还包括根据所述校正烧蚀模式通过激光烧蚀修改所述光学组织表面的步骤。
27.根据权利要求22所述的方法,还包括增加平均梯度场,以从被重建的表面消除倾斜。
28.根据权利要求22所述的方法,其中测量所述表面梯度包括在孔径上均匀地对所传输的图像进行取样。
29.根据权利要求28所述的方法,其中所述孔径是眼睛的瞳孔。
30.根据权利要求22所述的方法,其中用Hartmann-Shack传感器组件进行测量表面梯度数据。
31.根据权利要求22所述的方法,其中所述表面是波阵面表面。
32.根据权利要求22所述的方法,其中应用傅立叶变换包括应用离散傅立叶分解和逆离散傅立叶变换。
33.根据权利要求22所述的方法,其中傅立叶变换在重建中使用所有可得到的信息。
34.根据权利要求22所述的方法,其中应用傅立叶变换计算所述眼睛的光学组织的X线断层摄影波阵面误差图。
35.根据权利要求22所述的方法,其中所述图像作为多个光束组(beamlet)由所述光学组织传输,所述表面梯度包括梯度的阵列,
其中每一梯度与所述眼睛的光学组织的相关部分相对应,每一光束组根据相应的梯度通过所述光学组织传输。
36.一种测量眼睛的光学组织的方法,所述方法包括以下步骤:
通过所述光学组织传输图像;
从该传输图像确定跨过所述光学组织的局部梯度;和
通过将傅立叶变换算法应用到跨过所述眼睛的光学组织的表面梯度绘制所述眼睛的波阵面误差图。
37.根据权利要求36所述的方法,还包括向所述波阵面误差添加平均梯度场以对倾斜进行校正。
38.根据权利要求36所述的方法,其中由Hartmann-Shack传感器组件进行确定跨过所述光学组织的局部梯度。
39.根据权利要求36所述的方法,包括基于绘制的所述眼睛的光学组织的波阵面误差图产生激光烧蚀处理表。
40.根据权利要求39所述的方法,包括根据所述校正烧蚀模式通过激光烧蚀修改所述光学组织表面。
41.一种用于测量眼睛组织的波阵面误差的系统,所述系统包括:
处理器;
与所述处理器耦合的存储器,所述存储器被配置为存储多个由所述处理器执行的代码模块,所述多个代码模块包括:
用于通过所述光学组织传输图像的模块;
用于从该传输的图像确定跨过所述光学组织的局部梯度的模块;和
用于通过应用傅立叶变换算法于跨过所述眼睛的光学组织的表面梯度绘制所述眼睛的波阵面误差图的模块。
42.根据权利要求41所述的系统,还包括与所述处理器耦合的用于通过所述眼睛的光学组织传输源图像的图像源。
43.根据权利要求41所述的系统,还包括与所述处理器耦合的波阵面传感器系统。
44.根据权利要求43所述的系统,其中所述波阵面传感器系统包括Hartmann-Shack传感器组件。
45.根据权利要求41所述的系统,其中所述代码模块还包括用于基于如由所述傅立叶重建表面表示的所述眼睛的光学组织计算校正烧蚀模式的模块。
46.一种与权利要求45所述的系统通信的激光系统,其中该激光系统包括可编程用于根据所述校正烧蚀模式向所述光学组织传送激光能量的激光器。
47.根据权利要求41所述的系统,还包括跟踪所述光学组织的位置的摄像机,
其中所述代码模块还包括用于登记关于所述光学组织的波阵面误差的模块。
48.根据权利要求41所述的系统,还包括与所述处理器耦合的自适应光学元件。
49.根据权利要求48所述的系统,其中所述自适应光学元件是可变形镜。
50.一种存储在计算机可读存储介质上用于测量光学组织的计算机程序,该计算机程序包括:
a)用于通过眼睛的光学组织传输图像的代码;
b)用于从通过所述眼睛的光学组织的被传输图像测量表面梯度的代码;和
c)用于通过应用傅立叶变换算法于跨过所述眼睛的光学组织的所述表面梯度绘制所述眼睛的波阵面误差图的代码。
51.根据权利要求50所述的计算机程序,还包括:
基于如由所述傅立叶重建表面表示的所述眼睛的光学组织计算校正烧蚀模式的代码。
52.根据权利要求51所述的计算机程序,还包括:
根据所述校正烧蚀模式向所述光学组织传送激光能量的代码。
53.根据权利要求50所述的计算机程序,还包括用于将绘制的波阵面误差图与所述眼睛的光学组织的图像对准的代码。
54.一种用于测量眼睛的光学组织的系统,包括:
用于通过所述光学组织传输图像的装置;
用于从该传输图像确定跨过所述光学组织的局部梯度的装置;和
用于通过应用傅立叶变换于跨过所述眼睛的光学组织的表面梯度绘制所述眼睛的波阵面误差图的装置。
55.根据权利要求54所述的系统,还包括基于如由所述傅立叶重建表面表示的所述眼睛的光学组织计算校正烧蚀模式的装置。
56.根据权利要求55所述的系统,还包括根据所述校正烧蚀模式通过激光烧蚀修改所述光学组织表面的装置。
CNB2004800209971A 2003-06-20 2004-06-17 用于激光外科手术和其它光学应用的迭代傅立叶重建 Expired - Fee Related CN100446718C (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
US10/601,048 US7175278B2 (en) 2003-06-20 2003-06-20 Wavefront reconstruction using fourier transformation and direct integration
US10/601,048 2003-06-20

Related Child Applications (1)

Application Number Title Priority Date Filing Date
CN2008100962701A Division CN101444418B (zh) 2003-06-20 2004-06-17 用于激光外科手术和其它光学应用的迭代傅立叶重建

Publications (2)

Publication Number Publication Date
CN1826080A true CN1826080A (zh) 2006-08-30
CN100446718C CN100446718C (zh) 2008-12-31

Family

ID=33517887

Family Applications (2)

Application Number Title Priority Date Filing Date
CN2008100962701A Expired - Fee Related CN101444418B (zh) 2003-06-20 2004-06-17 用于激光外科手术和其它光学应用的迭代傅立叶重建
CNB2004800209971A Expired - Fee Related CN100446718C (zh) 2003-06-20 2004-06-17 用于激光外科手术和其它光学应用的迭代傅立叶重建

Family Applications Before (1)

Application Number Title Priority Date Filing Date
CN2008100962701A Expired - Fee Related CN101444418B (zh) 2003-06-20 2004-06-17 用于激光外科手术和其它光学应用的迭代傅立叶重建

Country Status (9)

Country Link
US (3) US7175278B2 (zh)
EP (2) EP1648290B1 (zh)
JP (1) JP5016307B2 (zh)
CN (2) CN101444418B (zh)
AT (1) ATE476907T1 (zh)
CA (1) CA2529328C (zh)
DE (1) DE602004028627D1 (zh)
MX (1) MXPA05013980A (zh)
WO (1) WO2004112588A2 (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102113000A (zh) * 2008-08-04 2011-06-29 卡尔蔡斯视觉股份有限公司 用于定制镜片的系统和方法
CN102821674A (zh) * 2010-02-05 2012-12-12 爱尔康手术激光股份有限公司 在激光手术系统中整合有局部成像的梯度搜索
CN104271087A (zh) * 2012-04-20 2015-01-07 威孚莱有限公司 用于控制角膜切除激光的技术

Families Citing this family (31)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7460288B2 (en) * 2002-12-06 2008-12-02 Amo Manufacturing Usa, Llc Methods for determining refractive corrections from wavefront measurements
US7175278B2 (en) * 2003-06-20 2007-02-13 Visx, Inc. Wavefront reconstruction using fourier transformation and direct integration
US7296893B2 (en) * 2004-03-03 2007-11-20 Visx, Incorporated Transformation methods of wavefront maps from one vertex distance to another
US7547102B2 (en) 2004-03-03 2009-06-16 Amo Manufacturing Usa, Llc Wavefront propagation from one plane to another
US8597282B2 (en) * 2005-01-13 2013-12-03 Amo Manufacturing Usa, Llc Database system for centralized clinical and research applications with data from wavefront aberrometers
US7331674B2 (en) * 2005-09-02 2008-02-19 Visx, Incorporated Calculating Zernike coefficients from Fourier coefficients
DE102005042436C5 (de) * 2005-09-07 2010-05-27 Carl Zeiss Surgical Gmbh Ophthalmo-Operationsmikroskop mit Messeinrichtung
AU2007227371B2 (en) * 2006-03-23 2012-02-02 Amo Manufacturing Usa, Llc Systems and methods for wavefront reconstruction for aperture with arbitrary shape
US9295584B2 (en) 2007-05-17 2016-03-29 Amo Development, Llc Customized laser epithelial ablation systems and methods
US7832864B2 (en) * 2007-06-15 2010-11-16 The Arizona Board Of Regents On Behalf Of The University Of Arizona Inverse optical design
US8529558B2 (en) * 2008-04-22 2013-09-10 Amo Development Llc. High-order optical correction during corneal laser surgery
DE102008046834A1 (de) * 2008-09-11 2010-03-18 Iroc Ag Steuerprogramm zum Steuern elektromagnetischer Strahlung für eine Quervernetzung von Augengewebe
US9504376B2 (en) * 2009-12-22 2016-11-29 Amo Wavefront Sciences, Llc Optical diagnosis using measurement sequence
CN101968383B (zh) * 2010-09-02 2012-01-25 北京理工大学 一种抗扰动的时频域波前检测方法
US8978660B2 (en) 2011-07-21 2015-03-17 Amo Development, Llc Tilt compensation, measurement, and associated adjustment of refractive prescriptions during surgical and other treatments of the eye
CN102721478B (zh) * 2012-07-10 2014-01-22 中国科学院光电技术研究所 一种应用于曲率波前传感器的波前复原方法
EP2903576A1 (en) * 2012-10-02 2015-08-12 AMO Development, LLC Systems and methods for treatment target deconvolution
US9332899B2 (en) * 2012-11-06 2016-05-10 Clarity Medical Systems, Inc. Electronic eye marking/registration
ITFI20120240A1 (it) * 2012-11-07 2014-05-08 Strumenti Oftalmici C S O S R L Costruzioni Metodo e apparato per la misura di aberrazioni del sistema ottico di un essere vivente
CN103162846B (zh) * 2013-02-07 2015-02-18 中国科学院光电技术研究所 一种构建Zernike多项式像差模式与Walsh函数像差模式之间系数转换矩阵的方法
US9554696B2 (en) 2013-11-26 2017-01-31 Abbott Medical Optics Inc. System and method for measuring dysphotopsia
WO2015114631A1 (en) * 2014-02-01 2015-08-06 Shimon Mizrahi Detecting and treating skin surface features
CN104976962B (zh) * 2014-04-09 2017-05-17 南京理工大学 基于共轭差分法测量平面镜绝对面形的方法
CN105318843B (zh) * 2014-06-04 2017-12-12 南京理工大学 一种应用共轭差分法检测柱面镜绝对面形的方法
CN104976964B (zh) * 2015-06-09 2017-07-14 中国科学院长春光学精密机械与物理研究所 球面镜和平面镜面形误差的绝对检测方法
US20180042771A1 (en) 2016-08-10 2018-02-15 Amo Development, Llc Epithelial ablation systems and methods
EP3522771B1 (en) 2016-10-25 2022-04-06 Amo Groningen B.V. Realistic eye models to design and evaluate intraocular lenses for a large field of view
US10739227B2 (en) 2017-03-23 2020-08-11 Johnson & Johnson Surgical Vision, Inc. Methods and systems for measuring image quality
EP3687447A1 (en) 2017-11-30 2020-08-05 AMO Groningen B.V. Intraocular lenses that improve post-surgical spectacle independent and methods of manufacturing thereof
CN111667548B (zh) * 2020-06-12 2022-03-29 暨南大学 一种多模式显微图像数值重建方法
CN114998161A (zh) * 2022-06-02 2022-09-02 中国科学院西安光学精密机械研究所 一种基于完美傅里叶变换的傅里叶叠层显微术高精度重构图像方法

Family Cites Families (63)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US3936160A (en) 1973-10-15 1976-02-03 Karlheinz Von Bieren Interferometer for the measurement of wavefront sections of general imaging systems including the human eye
US4732148A (en) 1983-11-17 1988-03-22 Lri L.P. Method for performing ophthalmic laser surgery
US4773414A (en) 1983-11-17 1988-09-27 Lri L.P. Method of laser-sculpture of the optically used portion of the cornea
US5207668A (en) 1983-11-17 1993-05-04 Visx Incorporated Method for opthalmological surgery
US4770172A (en) 1983-11-17 1988-09-13 Lri L.P. Method of laser-sculpture of the optically used portion of the cornea
US5219343A (en) 1983-11-17 1993-06-15 Visx Incorporated Apparatus for performing ophthalmogolical surgery
US4665913A (en) 1983-11-17 1987-05-19 Lri L.P. Method for ophthalmological surgery
US5108388B1 (en) 1983-12-15 2000-09-19 Visx Inc Laser surgery method
US4669466A (en) 1985-01-16 1987-06-02 Lri L.P. Method and apparatus for analysis and correction of abnormal refractive errors of the eye
US5163934A (en) 1987-08-05 1992-11-17 Visx, Incorporated Photorefractive keratectomy
US5233517A (en) 1990-04-30 1993-08-03 Jindra Lawrence F Early glaucoma detection by Fourier transform analysis of digitized eye fundus images
WO1992001417A1 (en) 1990-07-19 1992-02-06 Horwitz Larry S Vision measurement and correction
US5258791A (en) 1990-07-24 1993-11-02 General Electric Company Spatially resolved objective autorefractometer
US5220360A (en) 1990-10-24 1993-06-15 Ophthalmic Imaging Systems, Inc. Apparatus and method for topographical analysis of the retina
US5144630A (en) 1991-07-29 1992-09-01 Jtt International, Inc. Multiwavelength solid state laser using frequency conversion techniques
CA2073802C (en) 1991-08-16 2003-04-01 John Shimmick Method and apparatus for combined cylindrical and spherical eye corrections
US5818957A (en) 1991-10-08 1998-10-06 Computed Anatomy Incorporated Processing of keratoscopic images
DE4232915A1 (de) 1992-10-01 1994-04-07 Hohla Kristian Vorrichtung zur Formung der Cornea durch Abtragen von Gewebe
US5520679A (en) 1992-12-03 1996-05-28 Lasersight, Inc. Ophthalmic surgery method using non-contact scanning laser
CO4230054A1 (es) 1993-05-07 1995-10-19 Visx Inc Metodo y sistemas para tratamiento con laser de errores refractivos utilizando formacion de imagenes de desplazamiento
US5745309A (en) * 1994-06-13 1998-04-28 The United States Of America As Represented By The United States Department Of Energy Method for removing tilt control in adaptive optics systems
US5646791A (en) 1995-01-04 1997-07-08 Visx Incorporated Method and apparatus for temporal and spatial beam integration
JPH08299269A (ja) 1995-05-10 1996-11-19 Nikon Corp 検影型他覚眼屈折力測定装置
US5782822A (en) 1995-10-27 1998-07-21 Ir Vision, Inc. Method and apparatus for removing corneal tissue with infrared laser radiation
US5936720A (en) 1996-07-10 1999-08-10 Neal; Daniel R. Beam characterization by wavefront sensor
US6130419A (en) 1996-07-10 2000-10-10 Wavefront Sciences, Inc. Fixed mount wavefront sensor
US5742626A (en) 1996-08-14 1998-04-21 Aculight Corporation Ultraviolet solid state laser, method of using same and laser surgery apparatus
US20010041884A1 (en) * 1996-11-25 2001-11-15 Frey Rudolph W. Method for determining and correcting vision
US6271914B1 (en) 1996-11-25 2001-08-07 Autonomous Technologies Corporation Objective measurement and correction of optical systems using wavefront analysis
US5777719A (en) * 1996-12-23 1998-07-07 University Of Rochester Method and apparatus for improving vision and the resolution of retinal images
US6090102A (en) 1997-05-12 2000-07-18 Irvision, Inc. Short pulse mid-infrared laser source for surgery
US6000800A (en) 1998-06-22 1999-12-14 Schepens Eye Research Institute Coaxial spatially resolved refractometer
US6099125A (en) 1998-12-07 2000-08-08 Schepens Eye Research Foundation Coaxial spatially resolved refractometer
US6004313A (en) 1998-06-26 1999-12-21 Visx, Inc. Patient fixation system and method for laser eye surgery
US6184974B1 (en) 1999-07-01 2001-02-06 Wavefront Sciences, Inc. Apparatus and method for evaluating a target larger than a measuring aperture of a sensor
US6199986B1 (en) 1999-10-21 2001-03-13 University Of Rochester Rapid, automatic measurement of the eye's wave aberration
US6419671B1 (en) 1999-12-23 2002-07-16 Visx, Incorporated Optical feedback system for vision correction
US6673062B2 (en) 2000-03-14 2004-01-06 Visx, Inc. Generating scanning spot locations for laser eye surgery
US6659613B2 (en) 2000-03-27 2003-12-09 Board Of Regents, The University Of Texas System Methods and systems for measuring local scattering and aberration properties of optical media
AU781722B2 (en) 2000-04-25 2005-06-09 Alcon Refractivehorizons, Inc. Spatial filter for enhancing Hartmann-Shack images and associated methods
US6460997B1 (en) 2000-05-08 2002-10-08 Alcon Universal Ltd. Apparatus and method for objective measurements of optical systems using wavefront analysis
PT1210003E (pt) * 2000-05-08 2004-11-30 Alcon Inc Medicao objectiva e correccao de sistemas opticos que utilizam a analise da frente de onda
US6471654B2 (en) * 2000-05-10 2002-10-29 Asahi Kogaku Kogyo Kabushiki Kaisha Ultrasonic endoscope
NL1015935C2 (nl) 2000-08-15 2002-02-18 Stichting Tech Wetenschapp Werkwijze voor het bepalen van de kwaliteit van een oogoptiek.
US6738511B1 (en) * 2000-10-04 2004-05-18 Veeco Instruments, Inc. Reduced noise sensitivity method and apparatus for converting an interferogram phase map to a surface profile map
EP1324689B1 (en) 2000-10-10 2006-08-02 University of Rochester Determination of ocular refraction from wavefront aberration data
US6827444B2 (en) 2000-10-20 2004-12-07 University Of Rochester Rapid, automatic measurement of the eye's wave aberration
WO2002034126A1 (en) 2000-10-20 2002-05-02 Wavefront Sciences, Inc. Method for computing visual performance from objective ocular aberration measurements
AU2002239515A1 (en) 2000-12-08 2002-06-18 Visx Incorporated Direct wavefront-based corneal ablation treatment program
US6736507B2 (en) 2001-01-22 2004-05-18 Alexis Kudryashov High resolution, multispectral, wide field of view retinal imager
US6331059B1 (en) 2001-01-22 2001-12-18 Kestrel Corporation High resolution, multispectral, wide field of view retinal imager
WO2002075367A2 (en) 2001-03-15 2002-09-26 Wavefront Sciences, Inc. Tomographic wavefront analysis system
CN100353907C (zh) * 2001-04-18 2007-12-12 博士伦公司 获得客观式显然验光的装置
US6709108B2 (en) * 2001-08-31 2004-03-23 Adaptive Optics Associates, Inc. Ophthalmic instrument with adaptive optic subsystem that measures aberrations (including higher order aberrations) of a human eye and that provides a view of compensation of such aberrations to the human eye
US6781681B2 (en) * 2001-12-10 2004-08-24 Ophthonix, Inc. System and method for wavefront measurement
US6924899B2 (en) * 2002-05-31 2005-08-02 Optical Physics Company System for measuring wavefront tilt in optical systems and method of calibrating wavefront sensors
JP4185337B2 (ja) * 2002-09-13 2008-11-26 株式会社トプコン 矯正要素判定装置及び方法
US7434936B2 (en) 2002-12-06 2008-10-14 Amo Manufacturing Usa, Llc Residual accommodation threshold for correction of presbyopia and other presbyopia correction using patient data
US7365893B2 (en) 2003-06-20 2008-04-29 Amo Manufacturing Usa, Llc Iterative Fourier reconstruction for laser surgery and other optical applications
US7168807B2 (en) 2003-06-20 2007-01-30 Visx, Incorporated Iterative fourier reconstruction for laser surgery and other optical applications
US7175278B2 (en) * 2003-06-20 2007-02-13 Visx, Inc. Wavefront reconstruction using fourier transformation and direct integration
US20060126019A1 (en) * 2004-12-10 2006-06-15 Junzhong Liang Methods and systems for wavefront analysis
US7331674B2 (en) * 2005-09-02 2008-02-19 Visx, Incorporated Calculating Zernike coefficients from Fourier coefficients

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102113000A (zh) * 2008-08-04 2011-06-29 卡尔蔡斯视觉股份有限公司 用于定制镜片的系统和方法
CN105844103A (zh) * 2008-08-04 2016-08-10 卡尔蔡斯视觉股份有限公司 用于定制镜片的系统和方法
US11521174B2 (en) 2008-08-04 2022-12-06 Carl Zeiss Vision Gmbh Systems and methods for determining a lens prescription
CN102821674A (zh) * 2010-02-05 2012-12-12 爱尔康手术激光股份有限公司 在激光手术系统中整合有局部成像的梯度搜索
CN102821674B (zh) * 2010-02-05 2015-06-17 爱尔康手术激光股份有限公司 在激光手术系统中整合有局部成像的梯度搜索
CN104271087A (zh) * 2012-04-20 2015-01-07 威孚莱有限公司 用于控制角膜切除激光的技术

Also Published As

Publication number Publication date
MXPA05013980A (es) 2006-03-02
US20100179793A1 (en) 2010-07-15
EP2229873A1 (en) 2010-09-22
EP1648290B1 (en) 2010-08-11
US7175278B2 (en) 2007-02-13
US7731363B2 (en) 2010-06-08
CN101444418A (zh) 2009-06-03
US20080212031A1 (en) 2008-09-04
EP1648290A2 (en) 2006-04-26
ATE476907T1 (de) 2010-08-15
DE602004028627D1 (de) 2010-09-23
EP1648290A4 (en) 2007-08-15
WO2004112588A2 (en) 2004-12-29
US20040257530A1 (en) 2004-12-23
CN101444418B (zh) 2012-07-18
EP2229873B1 (en) 2015-12-16
CA2529328A1 (en) 2004-12-29
CA2529328C (en) 2015-11-24
JP2007528745A (ja) 2007-10-18
WO2004112588A3 (en) 2005-03-24
JP5016307B2 (ja) 2012-09-05
CN100446718C (zh) 2008-12-31
US8228586B2 (en) 2012-07-24

Similar Documents

Publication Publication Date Title
CN1826080A (zh) 用于激光外科手术和其它光学应用的迭代傅立叶重建
CN1154431C (zh) 估计用于矫正病人眼睛散光的外科手术参数的方法
US7717562B2 (en) Scaling Zernike coefficients to smaller pupil sizes for refractive treatments
CA2687100C (en) Combined wavefront and topography systems and methods
US7547102B2 (en) Wavefront propagation from one plane to another
US7331674B2 (en) Calculating Zernike coefficients from Fourier coefficients
CN1575146A (zh) 根据波前象差数据确定眼睛折光度
CN1650148A (zh) 跟踪扭转的眼睛的方向和位置
CN101048691A (zh) 使用具有受控球面像差范围和中心遮拦孔径的多焦点透镜的扩展景深
US7168807B2 (en) Iterative fourier reconstruction for laser surgery and other optical applications
US7365893B2 (en) Iterative Fourier reconstruction for laser surgery and other optical applications

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20081231

Termination date: 20200617