CN116309194B - Oct图像畸变校正方法 - Google Patents
Oct图像畸变校正方法 Download PDFInfo
- Publication number
- CN116309194B CN116309194B CN202310590771.XA CN202310590771A CN116309194B CN 116309194 B CN116309194 B CN 116309194B CN 202310590771 A CN202310590771 A CN 202310590771A CN 116309194 B CN116309194 B CN 116309194B
- Authority
- CN
- China
- Prior art keywords
- scanning
- coordinates
- sample
- coordinate system
- oct
- 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
- 238000000034 method Methods 0.000 title claims abstract description 47
- 238000012937 correction Methods 0.000 title claims abstract description 24
- 238000007781 pre-processing Methods 0.000 claims abstract description 19
- 238000012163 sequencing technique Methods 0.000 claims abstract description 7
- 238000009434 installation Methods 0.000 claims abstract description 5
- 230000009466 transformation Effects 0.000 claims description 5
- 238000000605 extraction Methods 0.000 claims description 4
- 238000012014 optical coherence tomography Methods 0.000 description 100
- 238000010586 diagram Methods 0.000 description 12
- 239000004973 liquid crystal related substance Substances 0.000 description 6
- 238000004364 calculation method Methods 0.000 description 2
- 238000003745 diagnosis Methods 0.000 description 2
- 201000010099 disease Diseases 0.000 description 2
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 2
- 239000011521 glass Substances 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000003325 tomography Methods 0.000 description 2
- 230000004075 alteration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000003759 clinical diagnosis Methods 0.000 description 1
- 238000003702 image correction Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000001356 surgical procedure Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/80—Geometric correction
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- 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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/13—Edge detection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/181—Segmentation; Edge detection involving edge growing; involving edge linking
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/50—Depth or shape recovery
-
- 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/10—Image acquisition modality
- G06T2207/10072—Tomographic images
- G06T2207/10101—Optical tomography; Optical coherence tomography [OCT]
-
- 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
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Quality & Reliability (AREA)
- Radiology & Medical Imaging (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Medical Informatics (AREA)
- Operations Research (AREA)
- General Health & Medical Sciences (AREA)
- Algebra (AREA)
- Health & Medical Sciences (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Microscoopes, Condenser (AREA)
Abstract
本发明涉及一种OCT图像畸变校正方法，包括：采集得到样品OCT图像，对所述OCT图像进行预处理，得到所有样品OCT图像的边缘轮廓数据；在新坐标系下获得三维数据；在振镜扫描起点坐标到扫描起点和终点之间的中心点坐标之间等间距选取N个点，选取第q个点的数据，进行多项式拟合，得到第一次拟合曲线，获得第一次拟合曲线的系数；对获得的第一次拟合曲线的系数进行排序，对每一项系数均进行多项式拟合，得到第二次拟合曲线，获得第二次拟合曲线的系数；建立聚焦物镜长距离安装模式下OCT设备的扫描深度畸变模型，实现对任意线扫任意样品OCT图像的畸变校正。本申请能够在聚焦物镜长工作距离模式下，准确获得校正后样品OCT图像。
Description
技术领域
本发明涉及图像处理技术领域，尤其涉及一种OCT图像畸变校正方法。
背景技术
光学层析成像技术（optical coherence tomography，OCT）在眼科领域应用广泛。对OCT样品扫描模式以及OCT图像畸变的研究越来越重要。通过重构的眼前节光学层析二维图像测量其深度参数，不仅有利于进行眼科疾病的临床诊断，而且对提高眼科疾病临床手术治疗中手术导航的精度和实时性具有重要意义：
一张未校正的任意线扫样品OCT图像不能准确描述样品的特征，如深度值、房角等；而经过校正后，基于校正后任意线扫样品OCT图像可以准确地描述样品特征，方便临床医生进行诊疗诊断。
然而，近年来对OCT图像畸变校正的研究中，包括对样品臂内部折射变换引起的畸变校正研究、OCT图像运动伪差校正研究等，上述研究主要应用于眼底OCT图像成像中。但是，对引入聚焦物镜情况下长工作距离造成的任意线扫眼前节OCT成像畸变校正的研究相对较少。
发明内容
有鉴于此，有必要提供一种OCT图像畸变校正方法，其能够在聚焦物镜长工作距离模式下，准确获得校正后样品OCT图像。
本发明提供一种OCT图像畸变校正方法，该方法包括如下步骤：S1，在OCT扫描坐标系中以规则的角度间隔β°对样品进行B扫描，采集得到180/β张样品OCT图像，对所述OCT图像进行预处理，得到所有样品OCT图像的边缘轮廓数据；其中，5°≤β≤15°，180/β为整数；S2，振镜扫描坐标组成新坐标系的x坐标轴、y坐标轴，样品OCT图像深度值组成新坐标系的z坐标轴；在新坐标系下获得三维数据，所述三维数据由振镜扫描坐标和样品OCT图像深度值组成；S3，在振镜扫描起点坐标到扫描起点和终点之间的中心点坐标之间等间距选取N个点，在步骤S2得到的三维数据中，选取第q个点的数据，进行多项式拟合，得到第一次拟合曲线，获得第一次拟合曲线的系数；其中，1≤q≤N；S4，基于q的顺序，对获得的第一次拟合曲线的系数进行排序，对每一项系数均进行多项式拟合，得到第二次拟合曲线，获得第二次拟合曲线的系数；S5，基于获得的第二次拟合曲线的系数，建立聚焦物镜长距离安装模式下OCT设备的扫描深度畸变模型，基于振镜扫描坐标实现对任意线扫任意样品OCT图像的畸变校正。
具体地，所述步骤S1包括：
步骤S1-1、在OCT扫描坐标系中，得到SS-OCT系统扫描振镜对样品进行扫描的起点坐标和终点坐标；
步骤S1-2、输入起点坐标和终点坐标，输出得到大小为的样品OCT图像共180/β张；
步骤S1-3、对180/β张样品OCT图像进行图像预处理操作，得到180/β张样品OCT图像的轮廓拟合数据，即180/β张样品OCT图像的深度值。
具体地，所述图像预处理操作包括：小波变换、高斯差分、轮廓提取、轮廓数据最小二乘多项式拟合。
具体地，所述的步骤S2包括：
步骤S2-1、根据所述步骤S1-1得到的扫描的起点坐标和终点坐标，及所述步骤S1-3得到的样品OCT图像的深度值，建立扫描振镜的新坐标系：新坐标系的x轴、y轴是振镜扫描坐标的x轴、y轴，用表示；新坐标系的z轴表示样品OCT图像的深度值用/>表示；新坐标系下x轴、y轴、z轴的三维数据表示为/>，建立扫描振镜k线扫描样品的扫描起点坐标/>与相对应的样品OCT图像深度值/>之间的关系，得到；建立扫描振镜k线扫描样品的扫描终点坐标/>与相对应的样品OCT图像深度值/>之间的关系，得到/>；扫描起点坐标/>与扫描终点坐标/>之间的任意一点为/>，建立其与相对应的样品OCT图像的扫描深度/>之间的关系，得到，其中/>，/>；
步骤S2-2、根据新的坐标系中所有振镜扫描的起点和终点坐标、180/β张样品OCT图像的深度值，得到新坐标系下1线扫描、2线扫描，……，k线扫描，……，180/β线扫描的三维数据。
具体地，所述的步骤S3包括：
步骤S3-1、在扫描起点坐标到扫描中心点坐标之间等间距选取N个点，在新坐标系下1线扫描、2线扫描、……、k线扫描、……、180/β线扫描的数据中，选取得到扫描起点坐标到中心点坐标之间N个点中的第q个点的数据，进行多项式拟合；
步骤S3-2、选取在扫描起点坐标到扫描中心点坐标之间等间距N个点中的第1个点得到的数据，进行多项式拟合，得到第一次拟合曲线的系数为：；
选取在扫描起点坐标到扫描中心点坐标之间等间距N个点中的第q个点得到的数据，进行多项式拟合，得到第一次拟合曲线的系数为：；
选取在扫描起点坐标到扫描中心点坐标之间等间距N个点中的第N个点得到的数据，进行多项式拟合，得到第一次拟合曲线的系数为。
具体地，所述的多项式拟合采用三维空间三次多项式拟合方法。
具体地，所述的步骤S4包括：
步骤S4-1、按照顺序q对所有第一次拟合曲线的系数进行排序，
得到新坐标系下与x坐标轴有关的数据:；
得到新坐标系下与y坐标轴有关的数据:；
得到新坐标系下与z坐标轴有关的数据:；
步骤S4-2、对步骤S4-1得到的数据进行多项式拟合，包括对新坐标系下与x坐标轴有关的数据：系数，系数/>，系数，系数/>同时进行多项式拟合；
对新坐标系下与y坐标轴有关的数据：系数，系数，系数/>，系数/>同时进行多项式拟合；
对新坐标系下与z坐标轴有关的数据：系数，系数，系数/>，系数/>同时进行多项式拟合。
具体地，所述的步骤S4还包括：步骤S4-3、步骤S4-2的多项式拟合采用四维空间两次多项式拟合方法；经过所述步骤S4-2后得到第二次拟合曲线的系数。
具体地，所述的步骤S5包括：
步骤S5-1、对于一张未校正的任意线扫样品OCT图像，在扫描振镜坐标系中，扫描振镜的扫描坐标是已知的，用表示，其中/>，/>，/>，W表示任意线扫样品OCT图像，/>表示扫描起点坐标，/>表示扫描终点坐标，通过计算坐标/>与坐标/>之间的距离与扫描半径r的比值，即比值；基于所述步骤S4-3获得的第二次拟合曲线的系数并建立一元二次方程，基于/>值计算得到第一次拟合曲线的系数即/>、、/>；
步骤S5-2、基于新坐标系下与x轴相关的系数建立一元三次方程，通过反解方程/>得到自变量t的值，把自变量t的值带入新坐标系下z轴相关线束建立的一元三次方程/>中，得到扫描坐标/>处扫描样品的深度值/>，得到所有扫描坐标/>（/>，/>，/>）的深度值/>；
步骤S5-3、通过公式计算得到扫描坐标/>处的畸变值；通过公式/>得到校正后任意线扫任意样品OCT图像的像素值/>。
本申请能够在聚焦物镜长距离安装模式下，获得校正后样品OCT图像，实现准确描述样品特征如眼前节深度参数、房角等，方便临床医生进行诊疗诊断。
附图说明
图1为本发明OCT图像畸变校正方法的流程图；
图2为本发明实施例提供的OCT图像采集方式示意图；
图3为本发明实施例提供的1、4、8、12、16线扫描得到样品OCT图像示意图；
图4为本发明实施例提供的图像预处理的流程示意图；
图5为本发明实施例提供的1、4、8、12、16线扫描得到样品OCT图像高斯差分预处理结果示意图；
图6为本发明实施例提供的1、4、8、12、16线扫描得到样品OCT图像边缘轮廓拟合结果示意图；
图7为本发明实施例提供的在新坐标系下样品所有OCT图像轮廓数据示意图；
图8为本发明实施例提供的在新坐标系下所有样品OCT图像边缘轮廓数据中第1个点的示意图；
图9为本发明实施例提供的在新坐标系下所有样品OCT图像边缘轮廓数据中第3个点的示意图；
图10为本发明实施例提供的在新坐标系下所有样品OCT图像边缘轮廓数据中第6个点的示意图；
图11（a）、图11（b）分别为本发明实施例提供的任意线扫样品OCT图像和畸变校正结果示意图。
具体实施方式
下面结合附图及具体实施例对本发明作进一步详细的说明。
请参阅图1，是本发明OCT图像畸变校正方法较佳实施例的作业流程图。
步骤S1，样品OCT图像预处理：在OCT扫描坐标系中以规则的角度间隔β°（5°≤β≤15°）对样品进行B扫描，采集得到180/β张样品OCT图像，进行图像预处理，得到所有样品OCT图像的边缘轮廓数据，即样品OCT图像的深度值。
具体包括：
步骤S1-1、在OCT扫描坐标系中，得到SS-OCT系统扫描振镜对样品进行扫描的起点坐标和终点坐标，包括以角度间隔β°继续B扫描的起点坐标和终点坐标，即扫描振镜1线扫描的起点坐标和终点坐标、扫描振镜2线扫描的起点坐标和终点坐标，……，扫描振镜180/β线扫描的起点坐标和终点坐标；需要说明的是180/β是整数。
步骤S1-2、输入一系列起点坐标和终点坐标，输出得到大小为的样品OCT图像共180/β张，样品OCT图像的像素坐标表示为/>，其中：/>表示行数，表示列数，第1，2，...，180/β张样品OCT图像是SS-OCT系统扫描振镜通过对样品进行1线扫描、2线扫描......,180/β线扫描得到的；
需要说明的是样品是表面光滑平行的玻璃片。
步骤S1-3、对180/β张样品OCT图像进行图像预处理操作，包括：小波变换、高斯差分、轮廓提取、轮廓数据最小二乘多项式拟合等，得到180/β张样品OCT图像的轮廓拟合数据，即为180/β张样品OCT图像的深度值，其中扫描振镜k线扫描（1≤k≤180/β）样品OCT图像的深度值用表示。
以下结合具体实施例及附图进行说明：
步骤S1-1、计算一系列以规则的角度间隔对样品进行B扫描的起点坐标和终点坐标，如图2所示，其中黑圆点表示扫描起点，黑箭头表示扫描终点，本实施例中进行了1线扫描、2线扫描，……，16线扫描，其中8线扫描指的是y线扫，16线扫描指的是x线扫。
步骤S1-2、将步骤S1-1计算得到的一系列起点坐标和终点坐标输入编写好的图像采集程序中，得到16张样品OCT图像，本发明实施例样品选用表面光滑、平行的玻璃片；如图3所示，给出了本发明实施例的1线扫描样品OCT图像、4线扫描样品OCT图像、8线扫描样品OCT图像、12线扫描样品OCT图像、16线扫描样品OCT图像，其中样品OCT图像大小为。
步骤S1-3、采用图像预处理算法如图4所示，包括小波变换、高斯差分、轮廓提取、轮廓数据最小二乘多项式拟合等算法，对16张样品OCT图像进行图像预处理；如图5所示给出了本发明实施例的1线扫描样品OCT图像高斯差分预处理结果、4线扫描样品OCT图像高斯差分预处理结果、8线扫描样品OCT图像高斯差分预处理结果、12线扫描样品OCT图像高斯差分预处理结果、16线扫描样品OCT图像高斯差分预处理结果；如图6所示给出了本发明实施例的1线扫描样品OCT图像边缘轮廓拟合结果、4线扫描样品OCT图像边缘轮廓拟合结果、8线扫描样品OCT图像边缘轮廓拟合结果、12线扫描样品OCT图像边缘轮廓拟合结果、16线扫描样品OCT图像边缘轮廓拟合结果。
步骤S2，获得新坐标系下的三维数据：振镜扫描坐标组成新坐标系的x坐标轴、y坐标轴，样品OCT图像深度值组成新坐标系的z坐标轴；在新坐标系下获得三维数据，所述三维数据由振镜扫描坐标和样品OCT图像深度值组成。
具体包括：
步骤S2-1、获得所述步骤S1-1得到的一系列振镜扫描的起点坐标和终点坐标，获得所述步骤S1-3得到的样品OCT图像的深度值；新坐标系的x轴、y轴是振镜扫描坐标的x轴、y轴，用表示；新坐标系的z轴表示样品OCT图像的深度值用/>表示；新坐标系下x轴、y轴、z轴的三维数据表示为/>，建立扫描振镜k线扫描样品的扫描起点坐标与相对应的样品OCT图像深度值/>之间的关系，得到/>；建立扫描振镜k线扫描样品的扫描终点坐标/>与相对应的样品OCT图像深度值/>之间的关系，得到/>；扫描起点坐标/>与扫描终点坐标/>之间的任意一点为/>，建立其与相对应的样品OCT图像的扫描深度/>之间的关系，得到/>，其中，/>。
步骤S2-2、所有振镜扫描的起点和终点坐标、180/β张样品OCT图像的深度值都在新的坐标系中进行步骤S2-2，得到新坐标系下1线扫描、2线扫描，……，k线扫描，……，180/β线扫描的三维数据，表示为：。
以下结合具体实施例及附图进行说明：
步骤S2-1、获得步骤S1-1计算得到的一系列扫描振镜的扫描起点坐标和扫描终点坐标，获得步骤S1-3得到的16张样品OCT图像边缘轮廓拟合结果数据；建立扫描振镜k，线扫描样品的扫描起点坐标/>和扫描终点坐标/>与相应的样品OCT图像的深度值/>之间的关系：
，/>，；
其中，是扫描起点坐标/>与扫描终点坐标/>之间的任意一点，/>，。
步骤S2-2、在新坐标系下，扫描振镜1,2...,16线扫描的振镜扫描坐标与16张样品OCT图像的深度值组成的三维数据如图7所示，公式表示为：。
步骤S3，获得第一次拟合曲线及其系数：在振镜扫描起点坐标到扫描起点和终点之间的中心点坐标之间等间距选取N个点，在所述步骤S2得到的三维数据中，选取第q个点（1≤q≤N)的数据，进行多项式拟合，得到第一次拟合曲线，获得第一次拟合曲线的系数。
具体包括：
步骤S3-1、扫描振镜1线扫描、2线扫描，……，k线扫描，……，180/β线扫描时的扫描中心点是相同的：；扫描振镜1线扫描、2线扫描、……、k线扫描、……，180/β线扫描时的扫描半径是相同的，扫描半径用r表示：；
在扫描起点坐标到扫描中心点坐标之间等间距选取N个点，在新坐标系下1线扫描、2线扫描，……，k线扫描，……，180/β线扫描的数据中，选取得到扫描起点坐标到中心点坐标之间N个点中的第q个点（1≤q≤N)的数据，进行多项式拟合。
步骤S3-2、选取在扫描起点坐标到扫描中心点坐标之间等间距N个点中的第1个点得到的数据，进行多项式拟合，得到第一次拟合曲线的系数为：；
选取在扫描起点坐标到扫描中心点坐标之间等间距N个点中的第q个点得到的数据，进行多项式拟合，得到第一次拟合曲线的系数为：；
选取在扫描起点坐标到扫描中心点坐标之间等间距N个点中的第N个点得到的数据，进行多项式拟合，得到第一次拟合曲线的系数为；
其中，，即对公式/>进行四舍五入得到整数即为的值；/>。
步骤S3-3、步骤S3-2中的多项式拟合方法采用的三维空间三次多项式拟合方法。
以下结合具体实施例及附图进行说明：
步骤S3-1、采用下面公式计算振镜1线扫描、2线扫描，……，k线扫描，……，16线扫描时的扫描中心点坐标：/>；
采用下面公式计算振镜1线扫描、2线扫描，……，k线扫描，……，16线扫描时的扫描半径r：；
在扫描起点坐标到中心点坐标之间等间距选取N=10个点（包括扫描起点），选取步骤S3-2得到的数据中第q个点（1≤q≤10)的数据为：
；
；其中，，即对公式/>进行四舍五入得到整数即为的值；/>；
如图8所示给出了选取步骤S2-2得到的数据中第1个点数据示意图；如图9所示给出了选取步骤S2-2得到的数据中第3个扫描点数据示意图；如图10所示给出了选取步骤S2-2得到的数据中第6个扫描点数据示意图。
步骤S3-2、对步骤S3-1得到的第1个扫描点数据（如图8所示）进行多项式拟合，即对、/>、同时进行多项式拟合，得到第一次拟合曲线的系数为：
；
对步骤S3-1得到的第q个点数据进行多项式拟合，即对系数，系数/>，系数/>同时进行多项式拟合，其中其中/>，即对公式/>进行四舍五入得到整数即为/>的值；/>得到第一次拟合曲线的系数为：
;
对步骤S3-1得到的第10个点数据进行多项式拟合，即对系数、系数/>、系数同时进行多项式拟合，得到第一次拟合曲线的系数为：
。
步骤S3-3、第一次多项式拟合方法采用三维空间三次多项式拟合方法。
步骤S4，获得第二次拟合曲线及其系数：基于q（1≤q≤N)的顺序，对所述步骤S3获得的第一次拟合曲线系数进行排序，对每一项系数都进行多项式拟合，得到第二次拟合曲线，获得第二次拟合曲线的系数。具体包括：
步骤S4-1、按照顺序q（1≤q≤N)，对所有第一次拟合曲线的系数进行排序，得到新坐标系下与x坐标轴有关的数据:
;得到新坐标系下与y坐标轴有关的数据:
;得到新坐标系下与z坐标轴有关的数据:
。
步骤S4-2、对步骤S4-1得到的数据进行多项式拟合，包括对新坐标系下与x坐标轴有关的数据：系数，系数/>，系数，系数/>同时进行多项式拟合；
对新坐标系下与y坐标轴有关的数据：系数，系数，系数/>，系数/>同时进行多项式拟合；
对新坐标系下与z坐标轴有关的数据：系数，系数，系数/>，系数/>同时进行多项式拟合。
步骤S4-3、步骤S4-2多项式拟合采用的方法是四维空间两次多项式拟合方法；经过所述步骤S4-2后得到第二次拟合曲线的系数。
以下结合具体实施例及附图进行说明：
步骤S4-1、将步骤S3-2得到的第一次拟合曲线的系数分类，分为三类：与新坐标系中x坐标轴有关、与新坐标系中y坐标轴有关、与新坐标系中z坐标轴有关；
与新坐标系中x坐标轴有关的数据为：
;其中q是在扫描起点坐标到中心点坐标之间等间距选取10个点中的第q个点，/>；
与新坐标系中y坐标轴有关的数据为：
;其中q是在扫描起点坐标到中心点坐标之间等间距选取10个点中的第q个点，/>；
将与新坐标系中z坐标轴有关的数据为：
;其中q是在扫描起点坐标到中心点坐标之间等间距选取10个点中的第q个点，/>。
步骤S4-2、对步骤S4-1得到的数据进行多项式拟合，得到第二次拟合曲线，包括对与新坐标系中x坐标轴有关的数据：系数，系数/>，系数/>，系数/>同时进行多项式拟合；
对与新坐标系中y坐标轴有关的数据：系数，系数，系数/>，系数/>同时进行多项式拟合；
对与新坐标系中z坐标轴有关的数据：系数，系数，系数/>，系数/>同时进行多项式拟合。
步骤S4-3、多项式拟合采用的方法是四维空间二次多项式拟合方法；将步骤S4-2得到的第二次拟合曲线的系数保存到json文件中。
步骤S5，样品OCT图像的畸变校正：基于所述步骤S4获得的第二次拟合曲线的系数，建立聚焦物镜长距离安装模式下OCT设备的扫描深度畸变模型，基于振镜扫描坐标实现对任意线扫任意样品OCT图像的畸变校正。具体包括：
步骤S5-1、对于一张未校正的任意线扫样品OCT图像，在扫描振镜坐标系中，扫描振镜的扫描坐标是已知的，用表示，其中/>，/>，/>，W表示任意线扫样品OCT图像，/>表示扫描起点坐标，/>表示扫描终点坐标，通过计算坐标/>与坐标/>之间的距离与扫描半径r的比值，即比值，基于所述步骤S4-3获得的第二次拟合曲线的系数并建立一元二次方程，基于/>值计算得到第一次拟合曲线的系数即/>、、/>。
步骤S5-2、基于新坐标系下与x轴相关的系数建立一元三次方程，通过反解方程/>得到自变量t的值，把自变量t的值带入新坐标系下z轴相关线束建立的一元三次方程/>中，得到扫描坐标/>处扫描样品的深度值/>，可以得到所有扫描坐标/>（/>，/>，/>）的深度值/>。
步骤S5-3、通过公式计算得到扫描坐标/>处的畸变值；通过公式/>得到校正后任意线扫任意样品OCT图像的像素值/>。
以下结合具体实施例及附图进行说明：
步骤S5-1、如图11（a）所示给出了一张未校正的任意线扫样品OCT图像，表示扫描振镜坐标系中，第j个扫描坐标，其中/>；通过计算坐标/>与坐标之间的距离与扫描半径r的比值，即比值/>，调用所述步骤S4-3保存的json文件中第二次拟合曲线的系数并建立一元二次方程，基于/>值计算得到第一次拟合曲线的系数即/>、/>、/>。
步骤S5-2、通过反解方程得到自变量t的值，把自变量t的值带入表达式/>中计算得到扫描坐标/>处扫描样品的深度值/>，通过计算可以得到所有扫描坐标/>的深度值/>。
步骤S5-3、通过公式计算得到扫描坐标处的畸变值；通过公式/>得到校正后任意线扫任意样品OCT图像的像素值/>，步骤S5-1中图11（a）所示的样品OCT图像校正后得到的最终图像如图11（b）所示；
由于聚焦物镜长距离安装，使得样品OCT图像出现较大的畸变如图11（a）中所示，经过校正后畸变消失了如图11（b）所示，本申请很好地实现了任意线扫样品OCT图像畸变校正。
虽然本发明参照当前的较佳实施方式进行了描述，但本领域的技术人员应能理解，上述较佳实施方式仅用来说明本发明，并非用来限定本发明的保护范围，任何在本发明的精神和原则范围之内，所做的任何修饰、等效替换、改进等，均应包含在本发明的权利保护范围之内。
Claims (8)
1.一种OCT图像畸变校正方法，其特征在于，该方法包括如下步骤：
S1，在OCT扫描坐标系中以规则的角度间隔β°对样品进行B扫描，采集得到180/β张样品OCT图像，对所述OCT图像进行预处理，得到所有样品OCT图像的边缘轮廓数据；其中，5°≤β≤15°，180/β为整数；
S2，振镜扫描坐标组成新坐标系的x坐标轴、y坐标轴，样品OCT图像深度值组成新坐标系的z坐标轴；在新坐标系下获得三维数据，所述三维数据由振镜扫描坐标和样品OCT图像深度值组成；
S3，在振镜扫描起点坐标到扫描起点和终点之间的中心点坐标之间等间距选取N个点，在步骤S2得到的三维数据中，选取第q个点的数据，进行多项式拟合，得到第一次拟合曲线，获得第一次拟合曲线的系数；其中，1≤q≤N；
S4，基于q的顺序，对获得的第一次拟合曲线的系数进行排序，对每一项系数均进行多项式拟合，得到第二次拟合曲线，获得第二次拟合曲线的系数；
S5，基于获得的第二次拟合曲线的系数，建立聚焦物镜长距离安装模式下OCT设备的扫描深度畸变模型，基于振镜扫描坐标实现对任意线扫任意样品OCT图像的畸变校正；
所述的步骤S5包括：
步骤S5-1、对于一张未校正的任意线扫样品OCT图像，在扫描振镜坐标系中，扫描振镜的扫描坐标是已知的，用表示，其中/>，/>，/>，W表示任意线扫样品OCT图像的列数，/>表示扫描起点坐标，/>表示扫描终点坐标，通过计算坐标/>与坐标/>之间的距离与扫描半径r的比值，即比值，其中，/>为各线扫描时的扫描中心点坐标；基于所述步骤S4获得的第二次拟合曲线的系数并建立一元二次方程，基于/>值计算得到第一次拟合曲线的系数即/>、/>、/>；
步骤S5-2、基于新坐标系下与x轴相关的系数建立一元三次方程，通过反解方程/>得到自变量t的值，把自变量t的值带入新坐标系下z轴相关线束建立的一元三次方程/>中，得到扫描坐标/>处扫描样品的深度值/>，得到所有扫描坐标/>（/>，/>，/>）的深度值；
步骤S5-3、通过公式计算得到扫描坐标/>处的畸变值；通过公式/>得到校正后任意线扫任意样品OCT图像的像素值/>，其中/>表示样品OCT图像的像素坐标。
2.如权利要求1所述的方法，其特征在于：所述步骤S1包括：
步骤S1-1、在OCT扫描坐标系中，得到SS-OCT系统扫描振镜对样品进行扫描的起点坐标和终点坐标；
步骤S1-2、输入起点坐标和终点坐标，输出得到大小为的样品OCT图像共180/β张；
步骤S1-3、对180/β张样品OCT图像进行图像预处理操作，得到180/β张样品OCT图像的轮廓拟合数据，即180/β张样品OCT图像的深度值。
3.如权利要求2所述的方法，其特征在于，所述图像预处理操作包括：小波变换、高斯差分、轮廓提取、轮廓数据最小二乘多项式拟合。
4.如权利要求3所述的方法，其特征在于，所述的步骤S2包括：
步骤S2-1、根据所述步骤S1-1得到的扫描的起点坐标和终点坐标，及所述步骤S1-3得到的样品OCT图像的深度值，建立扫描振镜的新坐标系：新坐标系的x轴、y轴是振镜扫描坐标的x轴、y轴，用表示；新坐标系的z轴表示样品OCT图像的深度值用/>表示；新坐标系下x轴、y轴、z轴的三维数据表示为/>，建立扫描振镜k线扫描样品的扫描起点坐标/>与相对应的样品OCT图像深度值/>之间的关系，得到/>；建立扫描振镜k线扫描样品的扫描终点坐标/>与相对应的样品OCT图像深度值之间的关系，得到/>；扫描起点坐标/>与扫描终点坐标之间的任意一点为/>，建立其与相对应的样品OCT图像的扫描深度/>之间的关系，得到/>，其中/>，/>；
步骤S2-2、根据新的坐标系中所有振镜扫描的起点和终点坐标、180/β张样品OCT图像的深度值，得到新坐标系下1线扫描、2线扫描，……，k线扫描，……，180/β线扫描的三维数据。
5.如权利要求4所述的方法，其特征在于，所述的步骤S3包括：
步骤S3-1、在扫描起点坐标到扫描中心点坐标之间等间距选取N个点，在新坐标系下1线扫描、2线扫描、……、k线扫描、……、180/β线扫描的数据中，选取得到扫描起点坐标到中心点坐标之间N个点中的第q个点的数据，进行多项式拟合；
步骤S3-2、选取在扫描起点坐标到扫描中心点坐标之间等间距N个点中的第1个点得到的数据，进行多项式拟合，得到第一次拟合曲线的系数为：；
选取在扫描起点坐标到扫描中心点坐标之间等间距N个点中的第q个点得到的数据，进行多项式拟合，得到第一次拟合曲线的系数为：；
选取在扫描起点坐标到扫描中心点坐标之间等间距N个点中的第N个点得到的数据，进行多项式拟合，得到第一次拟合曲线的系数为。
6.如权利要求5所述的方法，其特征在于，所述的多项式拟合采用三维空间三次多项式拟合方法。
7.如权利要求6所述的方法，其特征在于，所述的步骤S4包括：
步骤S4-1、按照顺序q对所有第一次拟合曲线的系数进行排序，
得到新坐标系下与x坐标轴有关的数据:；
得到新坐标系下与y坐标轴有关的数据:；
得到新坐标系下与z坐标轴有关的数据:；
步骤S4-2、对步骤S4-1得到的数据进行多项式拟合，包括对新坐标系下与x坐标轴有关的数据：系数，系数/>，系数/>，系数/>同时进行多项式拟合；
对新坐标系下与y坐标轴有关的数据：系数，系数，系数/>，系数/>同时进行多项式拟合；
对新坐标系下与z坐标轴有关的数据：系数，系数，系数/>，系数/>同时进行多项式拟合。
8.如权利要求7所述的方法，其特征在于，所述的步骤S4还包括：步骤S4-3、步骤S4-2的多项式拟合采用四维空间两次多项式拟合方法；经过所述步骤S4-2后得到第二次拟合曲线的系数。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310590771.XA CN116309194B (zh) | 2023-05-24 | 2023-05-24 | Oct图像畸变校正方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310590771.XA CN116309194B (zh) | 2023-05-24 | 2023-05-24 | Oct图像畸变校正方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116309194A CN116309194A (zh) | 2023-06-23 |
| CN116309194B true CN116309194B (zh) | 2023-08-08 |
Family
ID=86813543
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310590771.XA Active CN116309194B (zh) | 2023-05-24 | 2023-05-24 | Oct图像畸变校正方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116309194B (zh) |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101949689A (zh) * | 2010-06-22 | 2011-01-19 | 深圳市斯尔顿科技有限公司 | 一种oct系统校正方法 |
| CN107862661A (zh) * | 2017-11-06 | 2018-03-30 | 郑州轻工业学院 | 一种光学相干层析成像系统图像校正方法 |
| CN113327202A (zh) * | 2021-03-30 | 2021-08-31 | 苏州微清医疗器械有限公司 | 一种图像畸变的矫正方法及其应用 |
| CN114581326A (zh) * | 2022-03-03 | 2022-06-03 | 上海应用技术大学 | 一种oct成像畸变矫正方法及装置 |
| CN115564822A (zh) * | 2021-07-01 | 2023-01-03 | 杭州三坛医疗科技有限公司 | 一种畸变校准方法、装置、电子设备及介质 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP6740177B2 (ja) * | 2017-06-14 | 2020-08-12 | キヤノン株式会社 | 画像処理装置、画像処理方法及びプログラム |
| JP7419042B2 (ja) * | 2019-11-29 | 2024-01-22 | キヤノン株式会社 | 医用画像処理装置、光干渉断層撮影装置、医用画像処理方法、及びプログラム |
-
2023
- 2023-05-24 CN CN202310590771.XA patent/CN116309194B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101949689A (zh) * | 2010-06-22 | 2011-01-19 | 深圳市斯尔顿科技有限公司 | 一种oct系统校正方法 |
| CN107862661A (zh) * | 2017-11-06 | 2018-03-30 | 郑州轻工业学院 | 一种光学相干层析成像系统图像校正方法 |
| CN113327202A (zh) * | 2021-03-30 | 2021-08-31 | 苏州微清医疗器械有限公司 | 一种图像畸变的矫正方法及其应用 |
| CN115564822A (zh) * | 2021-07-01 | 2023-01-03 | 杭州三坛医疗科技有限公司 | 一种畸变校准方法、装置、电子设备及介质 |
| CN114581326A (zh) * | 2022-03-03 | 2022-06-03 | 上海应用技术大学 | 一种oct成像畸变矫正方法及装置 |
Non-Patent Citations (1)
| Title |
|---|
| BINGYAO TAN et.al.Ultrawide field, distortion-corrected ocular shape estimation with MHz optical coherence tomography (OCT).《Biomedical Optics Express》.2021,第12卷(第9期),第5770-5781页. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116309194A (zh) | 2023-06-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107481228B (zh) | 基于计算机视觉的人体背部脊柱侧弯角度测量方法 | |
| JP5832523B2 (ja) | 光コヒーレンストモグラフィのための動き補正および画像改善の方法および装置 | |
| JP5952564B2 (ja) | 画像処理装置、画像処理方法 | |
| CN111242866B (zh) | 观测者动态眼位条件下ar-hud虚像畸变校正的神经网络插值方法 | |
| CN107564091A (zh) | 一种基于快速对应点搜索的三维重建方法及装置 | |
| CN116309194B (zh) | Oct图像畸变校正方法 | |
| CN108510584A (zh) | 脊椎骨旋转角度计算方法 | |
| EP3861524B1 (en) | Method for automatic shape quantification of an optic nerve head | |
| EP3216387B1 (en) | Method and system for motion artefacts removal in optical coherence tomograpy | |
| CN113034608B (zh) | 一种角膜表面形态测量装置及方法 | |
| CN113936069B (zh) | 一种用于光声断层成像的阵元虚拟插值方法 | |
| CN115690389A (zh) | 一种基于深度学习的白内障手术中角膜中心定位系统 | |
| CN112381952B (zh) | 一种基于多相机的面部轮廓点云模型重构方法及装置 | |
| CN114663544A (zh) | 一种基于深度图像先验的电阻抗图像重建方法 | |
| Andrisani et al. | Applications of PDEs inpainting to magnetic particle imaging and corneal topography | |
| CN111428634A (zh) | 一种采用六点法分块模糊加权的人眼视线追踪定位方法 | |
| CN108269258B (zh) | 从oct角膜图像分割角膜结构的方法和系统 | |
| CN114343565B (zh) | 一种光学相干层析视网膜图像校正方法及装置 | |
| JP2021518771A (ja) | Ｏｃｔにおけるｌｓｏベースの追跡を改良するための後処理方法 | |
| CN114184581B (zh) | 基于oct系统的图像优化方法、装置、电子设备及存储介质 | |
| CN116309594B (zh) | 眼前节oct图像处理方法 | |
| JP6598963B2 (ja) | 画像処理装置、画像処理方法およびプログラム | |
| CN116184668A (zh) | 一种hud动态畸变抑制方法和系统 | |
| Du et al. | A Regularization Method for Deconvolution of Optical Coherence Tomography Image | |
| Antal et al. | Virtually and depth sensor generated Moire pictures in screening and treatment of scoliosis |
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 |