CN108335328B - 摄像机姿态估计方法和摄像机姿态估计装置 - Google Patents

摄像机姿态估计方法和摄像机姿态估计装置 Download PDF

Info

Publication number
CN108335328B
CN108335328B CN201710044456.1A CN201710044456A CN108335328B CN 108335328 B CN108335328 B CN 108335328B CN 201710044456 A CN201710044456 A CN 201710044456A CN 108335328 B CN108335328 B CN 108335328B
Authority
CN
China
Prior art keywords
depth image
feature point
current frame
image
color image
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
CN201710044456.1A
Other languages
English (en)
Other versions
CN108335328A (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.)
Fujitsu Ltd
Original Assignee
Fujitsu Ltd
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 Fujitsu Ltd filed Critical Fujitsu Ltd
Priority to CN201710044456.1A priority Critical patent/CN108335328B/zh
Publication of CN108335328A publication Critical patent/CN108335328A/zh
Application granted granted Critical
Publication of CN108335328B publication Critical patent/CN108335328B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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/10016—Video; Image sequence
    • 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/10024—Color image
    • 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/30241—Trajectory
    • 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/30244—Camera pose

Abstract

公开了一种摄像机姿态估计方法和摄像机姿态估计装置，摄像机姿态估计方法包括：获得彩色图像和深度图像；获取当前帧彩色图像的特征点和相邻帧彩色图像的相应匹配特征点；获取与当前帧彩色图像中的特征点的位置分别对应的当前帧深度图像中的特征点，获取与当前帧深度图像的特征点分别匹配的相邻帧深度图像中的相应匹配特征点；基于当前摄像机姿态和用于校正图像中的特征点的非刚性畸变的参数表示当前帧深度图像中的特征点的世界坐标；在世界坐标系下计算当前帧深度图像和相邻帧深度图像的每对匹配特征点之间的距离，以针对每对匹配的特征点所计算出的距离之和最小为目标，估计当前摄像机姿态和所述参数。可以消除摄像机参数估计中的非刚性畸变。

Description

摄像机姿态估计方法和摄像机姿态估计装置
技术领域
本公开涉及摄像机的跟踪和定位领域，具体涉及一种消除摄像机参数估计中的非刚性畸变的摄像机姿态估计方法和摄像机姿态估计装置。
背景技术
摄像机的跟踪及定位有着非常重要的应用，比如增强现实、三维建模以及机器人导航。例如，对于机器人来说，它需要知道自己在真实世界中的实时位置才能实现导航。尽管这个问题已经被研究了很多年，但是如何提高摄像机姿态估计的准确性和鲁棒性，仍然是一个非常有挑战的问题。
传统的实现摄像机姿态估计的方法可以分为两大类：基于视觉的方法和基于深度的方法。在基于视觉的方法中，主要通过建立稀疏匹配点之间的几何关系来求解摄像机参数。这种方法的缺点在于，参数估计对于彩色图像的质量非常敏感。当没有足够丰富的匹配点时，结果会变得很差。基于深度的方法只利用稠密的深度信息来估计摄像机的姿态，常用的算法为迭代的最近点算法。但是这种方法对参数初始值的精度要求较高，估计较大角度姿态的能力不足。更为重要地是，由于摄像机姿态估计模型中可能存在多种畸变，这会导致姿态估计的错误。因此，如何校正这些畸变也是一个非常关键的问题。
发明内容
在下文中给出了关于本公开的简要概述，以便提供关于本公开的某些方面的基本理解。但是，应当理解，这个概述并不是关于本公开的穷举性概述。它并不是意图用来确定本公开的关键性部分或重要部分，也不是意图用来限定本公开的范围。其目的仅仅是以简化的形式给出关于本公开的某些概念，以此作为稍后给出的更详细描述的前序。
鉴于以上问题，本公开的目的是提供一种消除摄像机参数估计中的非刚性畸变的摄像机姿态估计方法和摄像机姿态估计装置。
根据本公开的一方面，提供了一种摄像机姿态估计方法，包括：通过摄像机获得彩色图像和与彩色图像对应的深度图像；获取当前帧彩色图像中的至少一个特征点和与所述至少一个特征点中的每个分别匹配的相邻帧彩色图像中的至少一个相应匹配特征点；获取与当前帧彩色图像中的所述至少一个特征点中的每个的位置分别对应的当前帧深度图像中的至少一个特征点，并且基于相邻帧彩色图像中的所述至少一个相应匹配特征点，获取与当前帧深度图像中的所述至少一个特征点中的每个分别匹配的相邻帧深度图像中的至少一个相应匹配特征点；基于当前摄像机姿态以及用于校正图像中的特征点的非刚性畸变的参数来表示当前帧深度图像中的每个特征点的世界坐标；以及基于当前帧深度图像中的每个特征点的世界坐标和相邻帧深度图像中的相应匹配特征点的世界坐标来计算该对匹配的特征点之间的距离，以及以针对每对匹配的特征点所计算出的距离之和最小为目标，估计当前摄像机姿态和所述参数。
根据本公开的另一方面，提供了一种摄像机姿态估计装置，包括：彩色图像和深度图像获取单元，被配置成通过摄像机获得彩色图像和与彩色图像对应的深度图像；彩色图像特征点获取单元，被配置成获取当前帧彩色图像中的至少一个特征点和与所述至少一个特征点中的每个分别匹配的相邻帧彩色图像中的至少一个相应匹配特征点；深度图像特征点获取单元，被配置成获取与当前帧彩色图像中的所述至少一个特征点中的每个的位置分别对应的当前帧深度图像中的至少一个特征点，并且基于相邻帧彩色图像中的所述至少一个相应匹配特征点，获取与当前帧深度图像中的所述至少一个特征点中的每个分别匹配的相邻帧深度图像中的至少一个相应匹配特征点；世界坐标表示单元，被配置成基于当前摄像机姿态以及用于校正图像中的特征点的非刚性畸变的参数来表示当前帧深度图像中的每个特征点的世界坐标；以及摄像机姿态和参数估计单元，被配置成基于当前帧深度图像中的每个特征点的世界坐标和相邻帧深度图像中的相应匹配特征点的世界坐标来计算该对匹配的特征点之间的距离，以及以针对每对匹配的特征点所计算出的距离之和最小为目标，估计当前摄像机姿态和所述参数。
根据本公开的其它方面，还提供了用于实现上述根据本公开的方法的计算机程序代码和计算机程序产品以及其上记录有该用于实现上述根据本公开的方法的计算机程序代码的计算机可读存储介质。
在下面的说明书部分中给出本公开实施例的其它方面，其中，详细说明用于充分地公开本公开实施例的优选实施例，而不对其施加限定。
附图说明
本公开可以通过参考下文中结合附图所给出的详细描述而得到更好的理解，其中在所有附图中使用了相同或相似的附图标记来表示相同或者相似的部件。所述附图连同下面的详细说明一起包含在本说明书中并形成说明书的一部分，用来进一步举例说明本公开的优选实施例和解释本公开的原理和优点。其中：
图1示出示出根据本公开的实施例的摄像机姿态估计方法的流程示例的流程图；
图2是示出根据本公开实施例的非刚性畸变的示例的图；
图3是示出根据本公开实施例的偏移向量的示例的图；
图4是示出根据本公开实施例的深度图像中的特征点的三维空间投影的示意图；
图5是示出根据本公开实施例的组中的深度图像中的特征点的三维空间投影的示意图；
图6是示出根据本公开的实施例的摄像机姿态估计装置的功能配置示例的框图；以及
图7是示出作为本公开的实施例中可采用的信息处理装置的个人计算机的示例结构的框图。
具体实施方式
在下文中将结合附图对本公开的示范性实施例进行描述。为了清楚和简明起见，在说明书中并未描述实际实施方式的所有特征。然而，应该了解，在开发任何这种实际实施例的过程中必须做出很多特定于实施方式的决定，以便实现开发人员的具体目标，例如，符合与系统及业务相关的那些限制条件，并且这些限制条件可能会随着实施方式的不同而有所改变。此外，还应该了解，虽然开发工作有可能是非常复杂和费时的，但对得益于本公开内容的本领域技术人员来说，这种开发工作仅仅是例行的任务。
在此，还需要说明的一点是，为了避免因不必要的细节而模糊了本公开，在附图中仅仅示出了与根据本公开的方案密切相关的设备结构和/或处理步骤，而省略了与本公开关系不大的其它细节。
根据本公开的一方面，提出了一种摄像机姿态估计方法，该摄像机姿态估计方法从彩色图像中得到稀疏的对应关系，然后利用非刚性校正方法消除可能存在的畸变，最后通过利用稀疏对应的几何约束以及非刚性校正实现摄像机姿态的估计。
下面结合附图详细说明根据本公开的实施例。
首先，将参照图1描述根据本公开实施例的摄像机姿态估计方法100的流程示例。图1是示出根据本公开的实施例的摄像机姿态估计方法100的流程示例的流程图。如图1所示，根据本公开的实施例的摄像机姿态估计方法100包括彩色图像和深度图像获取步骤S102、彩色图像特征点获取步骤S104、深度图像特征点获取步骤S106、世界坐标表示步骤S108以及摄像机姿态和参数估计步骤S110。
在彩色图像和深度图像获取步骤S102中，可以通过摄像机获得彩色图像和与彩色图像对应的深度图像。这样，可以分别获得彩色图像序列以及包括与彩色图像序列中的每个彩色图像分别对应的深度图像的深度图像序列。
在彩色图像特征点获取步骤S104中，可以获取当前帧彩色图像中的至少一个特征点和与所述至少一个特征点中的每个分别匹配的相邻帧彩色图像中的至少一个相应匹配特征点。
优选地，当前帧彩色图像的所述至少一个特征点和相邻帧彩色图像的所述至少一个相应匹配特征点是通过如下方式获取的：在当前帧彩色图像中检测所述至少一个特征点，并且针对所检测到的所述至少一个特征点中的每个特征点，在相邻帧彩色图像中查找其相应匹配特征点，从而可以获取相邻帧彩色图像中的所述至少一个相应匹配特征点。
作为具体示例，从当前帧彩色图像中进行特征点的检测，从而可以获取当前帧彩色图像中的至少一个特征点。常用的检测方法包括SIFT(尺度不变特征变换)算法、SURF(加速鲁棒特征)算法以及KLT(Kanade-Lucas-Tomasi)算法等。然后在相邻的两帧彩色图像之间进行特征点的匹配，即，针对当前帧彩色图像的所检测到的所述至少一个特征点中的每个特征点，在相邻帧彩色图像中查找其相应匹配特征点，从而获取相邻帧彩色图像中的至少一个相应匹配特征点。每对匹配的特征点随着时间反投影为一个特征跟踪轨迹，特征跟踪轨迹中的每个点对应三维空间中的一个点。令
Figure BDA0001214091780000051
表示第i个轨迹，其表示第i对匹配的特征点随着时间反投影所形成的一个特征跟踪轨迹，其中，i＝0,1,…,N-1，N(N≥1)是相邻两个彩色图像之间的匹配的特征点对的总数量；k为图像的序号；
Figure BDA0001214091780000052
为第k帧彩色图像上与轨迹χi对应的一个点(具体地，
Figure BDA0001214091780000053
表示第k帧彩色图像上的第i个特征点，
Figure BDA0001214091780000054
可以表示为矢量形式)，f(χi)表示轨迹χi经过的所有图像帧的集合。这些特征跟踪轨迹代表稀疏的对应关系，将用来估计摄像机的参数。
优选地，在获取相邻帧彩色图像的匹配特征点时，只有在相邻帧彩色图像中的匹配特征点存在于当前帧彩色图像之前的预定数量的帧中的情况下，才将该匹配特征点确定为最终的匹配特征点，以提高特征点匹配的精度。
对于上述轨迹χi，可以设置一个最小的轨迹长度来去除错误的匹配点。作为具体示例，只有轨迹χi的长度为预定数量的情况下，才将该匹配特征点确定为最终的匹配特征点，以提高特征点匹配的精度。其中，可以根据经验确定上述预定数量。
在深度图像特征点获取步骤S106中，可以获取与当前帧彩色图像中的所述至少一个特征点中的每个的位置分别对应的当前帧深度图像中的至少一个特征点，并且可以基于相邻帧彩色图像中的所述至少一个相应匹配特征点，获取与当前帧深度图像中的所述至少一个特征点中的每个分别匹配的相邻帧深度图像中的至少一个相应匹配特征点。
优选地，当前帧深度图像的所述至少一个特征点和相邻帧深度图像的所述至少一个相应匹配特征点是通过如下方式获取的：可以将当前帧彩色图像中的所述至少一个特征点映射到当前帧深度图像中，以获取当前帧深度图像中的所述至少一个特征点；以及可以将相邻帧彩色图像中的所述至少一个特征点映射到相邻帧深度图像中，以获取相邻帧深度图像中的所述至少一个特征点。
作为具体示例，可以将当前帧彩色图像中某个位置处的特征点映射到当前帧深度图像中的对应位置处，从而获取与当前帧彩色图像中的该特征点对应的当前帧深度图像中的一个特征点，并且可以依次将当前帧彩色图像中的每个特征点映射到当前帧深度图像中的对应位置处，从而获取与当前帧彩色图像中的每个特征点分别对应的当前帧深度图像中的特征点；可以将相邻帧彩色图像中某个位置处的特征点映射到相邻帧深度图像中的对应位置处，从而获取与相邻帧彩色图像中的该特征点对应的相邻帧深度图像中的一个特征点，并且可以依次将相邻帧彩色图像中的每个特征点映射到相邻帧深度图像中的对应位置处，从而获取与相邻帧彩色图像中的每个特征点分别对应的相邻帧深度图像中的特征点。
在世界坐标表示步骤S108中，可以基于当前摄像机姿态以及用于校正图像中的特征点的非刚性畸变的参数来表示当前帧深度图像中的每个特征点的世界坐标。
在一个例子中，摄像机姿态包含三个旋转参数以及三个平移参数，假设采用无径向畸变和无切向畸变的针孔摄像机模型，可以将第k帧图像对应的摄像机姿态表示为一个刚性变换矩阵的形式：
Figure BDA0001214091780000061
其中R为3*3的旋转矩阵，t为3*1的平移向量。变换矩阵Tg,k可以将第k帧的摄像机坐标转化为世界坐标，比如在摄像机坐标下的三维点Xc可以通过公式
Figure BDA0001214091780000062
转换为世界坐标下的三维点Xg。其中
Figure BDA0001214091780000063
表示齐次向量：
Figure BDA0001214091780000064
以下为了方便描述，假设当前帧彩色图像为第k帧彩色图像，并假设第k帧彩色图像的相邻帧彩色图像为第k-1帧彩色图像，因此，当前帧深度图像为第k帧深度图像，第k帧深度图像的相邻帧深度图像为第k-1帧深度图像。
为了实现摄像机姿态的估计，需要建立匹配的特征点之间的几何对应关系。首先，假设第k帧彩色图像上的一个特征点为
Figure BDA0001214091780000065
根据与第k帧彩色图像对应的第k帧深度图像，可获取该特征点
Figure BDA0001214091780000066
在第k帧深度图像中的相同位置处的对应特征点及其深度值
Figure BDA0001214091780000071
同时根据摄像机校正矩阵K，可将第k帧深度图像中的对应特征点反投影到当前摄像机坐标下的三维点Xc，k(u)，公式如下：
Figure BDA0001214091780000072
其中，摄像机校正矩阵K通常是常值，表达式如下：
Figure BDA0001214091780000073
在上式中，f为焦距，(cx，cy)为主点坐标。K的作用为将摄像机坐标转化为图像平面坐标。
上文中以彩色图像为例介绍了χi。本领域技术人员可以容易理解，因为彩色图像序列中的每个彩色图像分别与深度图像序列中的每个深度图像对应，因此深度图像中的每对匹配的特征点随着时间也反投影为一个特征跟踪轨迹，特征跟踪轨迹中的每个点为三维空间中的一个点。除非特别说明，在以下描述中，χi指的是深度图像中的每对匹配的特征点随着时间反投影所形成的特征跟踪轨迹。
在世界坐标系下，第k帧深度图像上与特征跟踪轨迹χi对应的特征点的三维坐标等价于第k帧深度图像上的该特征点处的三维坐标，则第k帧深度图像上与特征跟踪轨迹χi对应的特征点(即，第k帧深度图像上的第i个特征点)在世界坐标系下的三维坐标可通过摄像机姿态Tg,k计算得到：
Figure BDA0001214091780000074
在摄像机姿态估计中，可能存在非刚性畸变。图2是示出根据本公开实施例的非刚性畸变的示例的图。在图2中，图像A中存在两个点a1和a2，图像B中的两个点b1和b2分别为点a1和a2经受了非刚性畸变所对应的点。为了消除摄像机姿态估计中可能存在的非刚性畸变，需要引入非刚性校正。此外，图2的图像B中的两个点b1’和b2’分别为点b1和b2经受了非刚性畸变校正后的点。由图2可以看出，经过非刚性畸变校正，点b1’和b2’的位置与点a1和a2的位置相对应。
作为具体示例，可以将用于校正图像中的特征点的非刚性畸变的参数反映在特征点的世界坐标下中。也就是说，可以基于第k帧对应的摄像机姿态以及用于校正图像中的特征点的非刚性畸变的参数来表示第k帧深度图像中的每个特征点的世界坐标。
优选地，所述参数包括用于校正在当前帧彩色图像中获取匹配的特征点时产生的匹配畸变的参数、用于校正当前帧彩色图像和当前帧深度图像不对齐所产生的不对齐畸变的参数以及用于校正在将当前帧深度图像中的特征点在摄像机坐标中投影时产生的投影畸变的参数中的至少一个。
作为具体示例，在摄像机姿态估计系统中，可能存在三种非刚性畸变。
第一种非刚性畸变是彩色图像上的非刚性畸变。由于不同彩色图像帧之间的角度不同会造成成像颜色并不一致，这使得进行彩色图像特征点匹配时会导致匹配上的位置与正确的位置之间存在一定的偏移。这种偏移属于一种非刚性畸变，会对摄像机姿态估计产生影响。因此，用于校正图像中的特征点的非刚性畸变的参数可以包括用于校正在当前帧彩色图像中获取匹配的特征点时产生的匹配畸变的参数。
第二种非刚性畸变是深度图像上的非刚性畸变。由于深度图像和彩色图像采集的时间以及位置很难完全保持一致，使得深度图像与其对应的彩色图像并不完全对齐。为了修正这种畸变，同样采样非刚性校正。因此，用于校正图像中的特征点的非刚性畸变的参数可以包括用于校正当前帧彩色图像和当前帧深度图像不对齐所产生的不对齐畸变的参数。
第三种非刚性畸变是三维空间中的非刚性畸变。由于传感器的噪声，造成三维空间点与成像点和光心并不共线，从而导致由根据公式(2)反投影计算出来的摄像机坐标中的三维点坐标会偏离正确的位置。这种畸变同样可以采用非刚性变换进行校正。因此，用于校正图像中的特征点的非刚性畸变的参数可以包括用于校正在将当前帧深度图像中的特征点在摄像机坐标中投影时产生的投影畸变的参数。
为了消除非刚性畸变，定义非刚性畸变校正函数F。对于发生非刚性畸变的点u(在下文中，u表示该点的位置向量)，该点处的非刚性畸变校正函数F可表示为：
F(u)＝u+v(u) (4)
在公式(4)中，v(u)为校正向量，F(u)为点u经过非刚性畸变校正后的位置向量。
v(u)可以表示为预定位置上的P(P为大于或等于1的整数)个点的偏移向量的线性组合，定义如下：
Figure BDA0001214091780000091
在公式(5)中，vp为第p个点处的偏移向量，σp为与vp对应的权重系数，权重系数σp可以根据第p个点与点u之间的距离被预先确定，例如，第p个点与点u之间的距离越近，则权重系数σp越大，而第p个点与点u之间的距离越远，则权重系数σp越小，p＝1,2,…P。另外，可以用
Figure BDA0001214091780000092
表示P个点的偏移向量的集合，其中，p＝1,2,…P。。
图3是示出根据本公开实施例的偏移向量的示例的图。在图3中，存在12个点并示出了其中4个点的偏移向量v1、v2、v3以及v4。如上所述，公式(5)中的v(u)可以表示为这12个点的偏移向量的线性组合。
优选地，用于校正匹配畸变的参数是当前帧彩色图像中的多个预定位置处的点的偏移向量，其中，通过对表示当前帧彩色图像中的每个特征点的位置的向量叠加所述偏移向量的线性组合，对该特征点的位置进行校正。
作为具体示例，假设u为第k帧彩色图像中的特征点(u表示该特征点的位置向量)，则定义第k帧彩色图像中的特征点u对应的匹配畸变校正函数为
Figure BDA0001214091780000101
可以将其表示为：
Figure BDA0001214091780000102
在公式(6)中，
Figure BDA0001214091780000103
是第k帧彩色图像上的第p个点处的二维偏移向量，
Figure BDA0001214091780000104
为与
Figure BDA0001214091780000105
对应的权重系数，权重系数
Figure BDA0001214091780000106
可以根据第p个点与点u之间的距离被预先确定，例如，第p个点与点u之间的距离越近，则权重系数
Figure BDA0001214091780000107
越大，而第p个点与点u之间的距离越远，则权重系数
Figure BDA0001214091780000108
越小，p＝1,2,…P，P为第k帧彩色图像中的预定位置的点的总数量。可以定义
Figure BDA0001214091780000109
为第k帧彩色图像上的所述预定位置处的点的P个偏移向量的集合。在以上表示中，在
Figure BDA00012140917800001010
以及
Figure BDA00012140917800001011
的上标中采用“1”以与其他非刚性畸变有关的参数相区分。从公式(6)可知，可以通过对表示第k帧彩色图像中的每个特征点的位置的向量u叠加偏移向量的线性组合，对该特征点的位置进行校正。
优选地，用于校正不对齐畸变的参数是当前帧深度图像中的多个预定位置处的点的偏移向量，其中，通过对表示当前帧深度图像中的每个特征点的位置的向量叠加所述偏移向量的线性组合，对该特征点的位置进行校正。
作为具体示例，假设u为第k帧深度图像中的特征点(u表示该特征点的位置向量)，则定义与第k帧深度图像中的特征点u对应的不对齐畸变校正函数为
Figure BDA0001214091780000111
可以将其表示为：
Figure BDA0001214091780000112
在公式(7)中，
Figure BDA0001214091780000113
是第k帧深度图像上的第m个点处的二维偏移向量，
Figure BDA0001214091780000114
为与
Figure BDA0001214091780000115
对应的权重系数，权重系数
Figure BDA0001214091780000116
可以根据第m个点与点u之间的距离被预先确定，例如，第m个点与点u之间的距离越近，则权重系数
Figure BDA0001214091780000117
越大，而第m个点与点u之间的距离越远，则权重系数
Figure BDA0001214091780000118
越小，m＝1,2,…M，M为第k帧深度图像中的预定位置的点的总数量。可以定义
Figure BDA0001214091780000119
为第k帧深度图像上的所述预定位置处的点的M个偏移向量的集合。在以上表示中，在
Figure BDA00012140917800001110
Figure BDA00012140917800001111
以及
Figure BDA00012140917800001112
的上标中采用“2”以与其他非刚性畸变有关的参数相区分。从公式(7)可知，可以通过对表示第k帧深度图像中的每个特征点的位置的向量u叠加偏移向量的线性组合，对该特征点的位置进行校正。
优选地，用于校正投影畸变的参数是当前帧深度图像中的多个预定位置处的点的摄像机坐标的偏移向量，其中，通过对表示当前帧深度图像中的每个特征点的摄像机坐标的向量叠加所述偏移向量的线性组合，对该特征点的摄像机坐标进行校正。
作为具体示例，定义第k帧深度图像中的特征点在摄像机坐标系下的三维点X(X表示特征点在摄像机坐标下的三维向量)对应的投影畸变校正函数为
Figure BDA0001214091780000121
可以将其表示为：
Figure BDA0001214091780000122
在公式(8)中，
Figure BDA0001214091780000123
是第k帧深度图像上的第q个位置处的点在摄像机坐标下的三维偏移向量，
Figure BDA0001214091780000124
为与
Figure BDA0001214091780000125
对应的权重系数，权重系数
Figure BDA0001214091780000126
可以如上所述被预先确定，q＝1,2,…Q，Q为第k帧深度图像中的预定位置的点的总数量。可以定义
Figure BDA0001214091780000127
为第k帧深度图像上的所述预定位置处的点的摄像机坐标的Q个偏移向量的集合。在以上表示中，在
Figure BDA0001214091780000128
以及
Figure BDA0001214091780000129
的上标中采用“3”以与其他非刚性畸变有关的参数相区分。从公式(8)可知，可以通过对表示当前帧深度图像中的每个特征点的摄像机坐标的向量X叠加偏移向量的线性组合，对该特征点的摄像机坐标进行校正。
可以将针对上述三种畸变的非刚性校正引入到摄像机姿态估计中。通过任意组合上述三种非刚性校正，公式(2)可以被重写为公式(9)到(15)的多种形式。
Figure BDA00012140917800001210
公式(9)是通过用
Figure BDA00012140917800001211
代替公式(2)中的
Figure BDA00012140917800001212
得到的，其中
Figure BDA00012140917800001213
表示对第k帧彩色图像中的第i个特征点进行匹配畸变校正。
Figure BDA00012140917800001214
公式(10)是通过用
Figure BDA0001214091780000131
代替公式(2)中的
Figure BDA0001214091780000132
得到的，其中
Figure BDA0001214091780000133
表示对第k帧深度图像中的第i个特征点进行不对齐畸变校正。
Figure BDA0001214091780000134
公式(11)是通过用
Figure BDA0001214091780000135
代替公式(2)中的
Figure BDA0001214091780000136
得到的，其中
Figure BDA0001214091780000137
表示对第k帧深度图像中的第i个特征点的摄像机坐标进行投影畸变校正。
Figure BDA0001214091780000138
公式(12)是通过用
Figure BDA0001214091780000139
代替公式(2)中的
Figure BDA00012140917800001310
并且用
Figure BDA00012140917800001311
代替公式(2)中的
Figure BDA00012140917800001312
得到的，表明对第k帧彩色图像中的第i个特征点进行匹配畸变校正并且对第k帧深度图像中的第i个特征点进行不对齐畸变校正。
Figure BDA00012140917800001313
公式(13)中的
Figure BDA00012140917800001314
表示对第k帧彩色图像中的第i个特征点进行匹配畸变校正，
Figure BDA00012140917800001315
表示对基于
Figure BDA00012140917800001316
计算出来的第k帧深度图像中的第i个特征点的摄像机坐标进行投影畸变校正。
Figure BDA0001214091780000141
公式(14)中的
Figure BDA0001214091780000142
表示对第k帧深度图像中的第i个特征点进行不对齐畸变校正，
Figure BDA0001214091780000143
表示对基于
Figure BDA0001214091780000144
计算出来的第k帧深度图像中的第i个特征点的摄像机坐标进行投影畸变校正。
Figure BDA0001214091780000145
在公式(15)中，
Figure BDA0001214091780000146
表示对第k帧彩色图像中的第i个特征点进行匹配畸变校正并且
Figure BDA0001214091780000147
表示对第k帧深度图像中的第i个特征点进行不对齐畸变校正，此外，
Figure BDA0001214091780000148
表示对基于
Figure BDA0001214091780000149
和
Figure BDA00012140917800001410
计算出来的第k帧深度图像中的第i个特征点的摄像机坐标进行投影畸变校正。
在公式(9)到(15)中，均包括了用于校正图像中的特征点的非刚性畸变的参数。
另外，如公式(3)所示，由于第k帧深度图像上的第i个特征点在世界坐标系下的三维坐标可通过摄像机姿态Tg,k计算得到，即
Figure BDA0001214091780000151
因此，结合公式(9)到(15)可知，可以基于第k帧深度图像对应的摄像机姿态Tg,k以及用于校正图像中的特征点的非刚性畸变的参数来表示第k帧深度图像中的每个特征点的世界坐标。
在摄像机姿态和参数估计步骤S110中，可以基于当前帧深度图像中的每个特征点的世界坐标和相邻帧深度图像中的相应匹配特征点的世界坐标来计算该对匹配的特征点之间的距离，以及以针对每对匹配的特征点所计算出的距离之和最小为目标，估计当前摄像机姿态和所述参数。
优选地，可以基于当前帧深度图像中的每个特征点的世界坐标和相邻帧深度图像中的相应匹配特征点的世界坐标，计算该特征点与其相应匹配特征点之间的距离。即，在世界坐标系下，计算当前深度图像中的每个特征点和相邻帧深度图像中的相应匹配特征点之间的距离。
图4是示出根据本公开实施例的深度图像中的特征点的三维空间投影的示意图。图4中右侧的图像B上的圆点b1和b2为第k帧深度图像中的两个特征点，图4中左侧的图像A上的圆点a1和a2分别为在作为第k帧深度图像的相邻帧深度图像的第k-1帧深度图像上与特征点b1和b2分别匹配的特征点。在理想情况下，相邻两个帧上的一对匹配特征点在世界坐标系中反投影为三维空间中的同一个点，例如，如图4所示，第k帧深度图像中的特征点b1和第k-1帧深度图像中的相应匹配特征点a1在世界坐标系中反投影为三维空间中的同一个点c1，并且第k帧深度图像中的特征点b2和第k-1帧深度图像中的相应匹配特征点a2在世界坐标系中反投影为三维空间中的同一个点c2。相反，图4还示出了第k帧深度图像中的特征点b2在世界坐标系中反投影为三维空间中的点c2’的情况，在这种情况下，第k帧深度图像中的特征点b2反投影形成的三维空间中的点c2’不同于第k-1帧深度图像的相应匹配特征点a2反投影形成的三维空间中的点c2，为了使得该匹配特征点对反投影为三维空间中的同一个点，应该使图4中示出的圆圈内的两个点c2和c2’在世界坐标系中的距离最小。根据这个原则我们可以通过最小化第k帧深度图像上的每个特征点与已求解的第k-1帧深度图像上的匹配特征点的世界坐标系下的距离之和来估计当前摄像机的姿态。即，通过使得在世界坐标系下的、针对当前帧深度图像(第k帧深度图像)的特征点和相邻帧深度图像(第k-1帧深度图像)的每对匹配特征点所计算出的距离之和最小，估计第k帧的摄像机姿态。作为具体示例的公式如下：
Figure BDA0001214091780000161
在公式(16)中，E()表示误差函数，
Figure BDA0001214091780000162
为与第k帧彩色图像中的第i个特征点
Figure BDA0001214091780000163
对应的第k帧深度图像中的第i个特征点在摄像机坐标下的三维坐标，可以由公式(2)求出，
Figure BDA0001214091780000164
表示与第k帧深度图像中的第i个特征点匹配的第k-1帧深度图像中的第i个特征点的世界坐标，其中，i＝0,1,…,N-1，N为第k帧深度图像和第k-1帧深度图像中的匹配的特征点对的总数量，
Figure BDA0001214091780000165
为在描述公式(6)中提到的、表示第k帧彩色图像上的预定位置处的点的P个偏移向量的集合，
Figure BDA0001214091780000166
为在描述公式(7)中提到的、表示第k帧深度图像上的预定位置处的点的M个偏移向量的集合，
Figure BDA0001214091780000167
为在描述公式(8)中提到的、表示第k帧深度图像上的预定位置处的点的摄像机坐标的Q个偏移向量的集合。
优选地，可以采用交替优化的方式来分别估计当前摄像机姿态和所述参数。
作为具体示例，由于在公式(16)中，Tg,k、
Figure BDA0001214091780000168
以及
Figure BDA0001214091780000169
都是未知变量，因此，可以采用交替优化的方式最小化公式(16)，从而估计
Figure BDA00012140917800001610
Figure BDA00012140917800001611
以及第k帧图像对应的摄像机的参数Tg,k。
优选地，可以将经非刚性畸变校正后的当前帧深度图像作为下一帧深度图像的相邻帧深度图像。作为具体示例，可以将经非刚性畸变校正后的第k帧深度图像作为第k+1帧深度图像的相邻帧深度图像，这样，可以进一步提供摄像机姿态参数估计的精度。
在以上描述中，示例而非限制，将第k-1帧图像作为第k帧图像的相邻帧。此外，还可以将第k帧图像之前的任意帧作为第k帧图像的相邻帧。
优选地，对于一个图像序列，将第一帧深度图像的摄像机坐标系设置为世界坐标系；以及摄像机姿态为用于将摄像机坐标系变换到世界坐标系的变换矩阵。
由公式(3)可以看出，摄像机姿态为用于将摄像机坐标系变换到世界坐标系的变换矩阵。在对一个图像序列进行处理时，设置第一帧(假设k＝1)深度图像的摄像机坐标系为世界坐标系，即将第一帧深度图像的摄像机姿态Tg,1设置为单位矩阵。那么，第一帧上的特征点在世界坐标下的三维点坐标可由公式(2)和公式(3)得到。对于第k帧图像，根据该帧上的特征点所在的特征轨迹，选取出那些已经被求解出三维点坐标
Figure BDA0001214091780000171
的第k-1帧上的特征点，其中，i＝0,1,…,N-1。对于第k帧深度图像中的每个特征点，可以根据公式(9)-(15)中的任一个得到该特征点在摄像机坐标下的三维点坐标
Figure BDA0001214091780000172
然后可以根据公式(3)得到该特征点在世界坐标下的三维点坐标
Figure BDA0001214091780000173
利用公式(16)，在世界坐标系下计算第k-1帧深度图像中的每个特征点和第k-1帧深度图像中相应匹配特征点之间的距离，并且通过使得针对每对匹配的特征点所计算出的距离之和最小，可以估计出第k帧对应的摄像机姿态，例如，在每次进行迭代时，首先固定
Figure BDA0001214091780000174
以及
Figure BDA0001214091780000175
优化Tg,k；然后固定Tg,k，分别优化
Figure BDA0001214091780000176
以及
Figure BDA0001214091780000177
优选地，在估计出当前摄像机姿态的情况下，基于所估计的当前摄像机姿态来计算当前帧深度图像中的所有特征点的世界坐标；以及对于当前帧深度图像中的每个特征点：如果该特征点也存在于当前帧深度图像之前的深度图像中，则将该特征点在当前帧深度图像中的世界坐标与该特征点在当前帧深度图像之前的深度图像中的世界坐标进行平均，作为该特征点在当前帧深度图像中的更新后的世界坐标。
作为具体示例，在估计出第k帧图像所对应的摄像机姿态Tg,k的情况下，基于所估计的摄像机姿态Tg,k来计算第k帧深度图像中的所有特征点的世界坐标；以及对于第k帧深度图像中的每个特征点：如果该特征点也存在于第k帧深度图像之前的深度图像中，则将该特征点在第k帧深度图像中的世界坐标与该特征点在第k帧深度图像之前的深度图像中的世界坐标进行平均，作为该特征点在第k帧深度图像中的更新后的世界坐标。换句话说，如果某个特征点所在的特征跟踪轨迹在第k帧之前已经求解出世界坐标系下的三维点坐标，则将它们取平均作为对应轨迹的新的世界坐标系下的三维点坐标。
在对一个图像序列进行处理时，这种渐进式的摄像机姿态求解方式容易造成时序上的误差累积，从而导致摄像机姿态的估计严重偏离正确的取值。为此，可以采取全局优化的方式来降低累积误差。全局优化即同时优化多帧的摄像机姿态。
优选地，根据本公开实施例的摄像机姿态估计方法还可以包括：将深度图像的序列划分成组；以及对于每个组，基于与该组中的每帧深度图像的所有特征点和该组中的所有其他帧深度图像的所有相应匹配特征点有关的组几何约束，优化该组中的每帧深度图像对应的摄像机姿态和所述参数。
优选地，所述组几何约束为：使得在世界坐标系下的、该组中的所有对匹配特征点之间的距离之和最小。
图5是示出根据本公开实施例的组中的深度图像中的特征点的三维空间投影的示意图。图5中图像C上的圆点c1为第k帧深度图像中的一个特征点，图像B上的圆点b1为在作为第k帧深度图像的相邻帧深度图像的第k-1帧深度图像上与特征点c1匹配的特征点，以及图像A上的圆点a1为在作为第k-1帧深度图像的相邻帧深度图像的第k-2帧深度图像上与特征点b1匹配的特征点。如图5所示，特征点a1在世界坐标系中反投影为点d1、特征点b1在世界坐标系中反投影为点d1’以及特征点c1在世界坐标系中反投影为点d1”。为了使得相互匹配的特征点反投影为三维空间中的同一个点，应该使得在世界坐标系中的、点d1与d1’之间的距离以及d1’与d1”之间的距离之和最小。根据这个原则我们可以通过使组中的所有对匹配特征点之间的距离之和最小来优化该组中的每帧深度图像对应的摄像机的姿态。
作为具体示例，下面以将第1帧至第k帧深度图像划分成一组为例，介绍对该组进行全局优化。对于该组，全局优化的损失函数可以表示为：
Figure BDA0001214091780000191
在公式(17)中，E()为误差函数；
Figure BDA0001214091780000192
和
Figure BDA0001214091780000193
表示在第m帧彩色图像和第n帧彩色图像中的第i对匹配特征点，
Figure BDA0001214091780000194
和
Figure BDA0001214091780000195
在同一条轨迹χi上；i＝0,1,…,M-1，M(M≥1)是该组中匹配的特征点对的总数量；
Figure BDA0001214091780000196
和
Figure BDA0001214091780000197
分别是
Figure BDA0001214091780000198
和
Figure BDA0001214091780000199
在第m帧深度图像和第n帧深度图像中的对应特征点反投影到摄像机坐标下的三维点；Tg,k、Tg,p、Tg,1、Tg,m以及Tg,n分别是第k帧、第p帧、第1帧、第m帧以及第n帧对应的摄像机姿态；
Figure BDA00012140917800001910
为第k帧彩色图像上的预定位置处的点的偏移向量的集合，V1 1为第1帧彩色图像上的预定位置处的点的偏移向量的集合；
Figure BDA00012140917800001911
为第k帧深度图像上的预定位置处的点的偏移向量的集合，V1 2为第1帧深度图像上的预定位置处的点的偏移向量的集合；
Figure BDA00012140917800001912
为第k帧深度图像上的预定位置处的点的摄像机坐标的偏移向量的集合，V1 3为第1帧深度图像上的预定位置处的点的摄像机坐标的偏移向量的集合。最小化公式(17)可以使得包括第1帧至第k帧深度图像的组中所有匹配特征点对应的世界坐标系下的三维点坐标尽可能地相近，从而降低了误差的累积，进一步提高了摄像机姿态估计的精度。全局优化可以在任意时刻进行，但该算法复杂度较高，每帧都进行全局优化是不必要的。可以通过监控当前帧上的误差值，当误差值超过某个阈值时，才采用全局优化的策略。
综上所述，根据本公开的实施例的摄像机姿态估计方法100结合视觉信息和深度信息，并且可以校正存在于深度图像、彩色图像以及三维空间中的非刚性畸变，从而能够实现摄像机参数的精确估计。
与上述摄像机姿态估计方法实施例相对应地，本公开还提供了以下摄像机姿态估计装置的实施例。
图6是示出根据本公开的实施例的摄像机姿态估计装置600的功能配置示例的框图。
如图6所示，根据本公开的实施例的摄像机姿态估计装置600可以包括彩色图像和深度图像获取单元602、彩色图像特征点获取单元604、深度图像特征点获取单元606、世界坐标表示单元608以及摄像机姿态和参数估计单元610。接下来将描述各个单元的功能配置示例。
在彩色图像和深度图像获取单元602中，可以通过摄像机获得彩色图像和与彩色图像对应的深度图像。这样，可以分别获得彩色图像序列以及包括与彩色图像序列中的每个彩色图像分别对应的深度图像的深度图像序列。
在彩色图像特征点获取单元604中，可以获取当前帧彩色图像中的至少一个特征点和与所述至少一个特征点中的每个分别匹配的相邻帧彩色图像中的至少一个相应匹配特征点。
优选地，当前帧彩色图像的所述至少一个特征点和相邻帧彩色图像的所述至少一个相应匹配特征点是通过如下方式获取的：在当前帧彩色图像中检测所述至少一个特征点，并且针对所检测到的所述至少一个特征点中的每个特征点，在相邻帧彩色图像中查找其相应匹配特征点，从而可以获取相邻帧彩色图像中的所述至少一个相应匹配特征点。
获取当前帧彩色图像的所述至少一个特征点和相邻帧彩色图像的所述至少一个相应匹配特征点的示例可以参见以上方法实施例中相应位置的描述，在此不再重复。
每对匹配的特征点随着时间反投影为一个特征跟踪轨迹，特征跟踪轨迹中的每个点为三维空间中的一个点。
Figure BDA0001214091780000201
表示第i个轨迹，其表示第i对匹配的特征点随着时间反投影所形成的一个特征跟踪轨迹。有关参数的具体说明可以参见以上方法实施例中相应位置的描述，在此不再重复。这些特征跟踪轨迹代表稀疏的对应关系，将用来估计摄像机的参数。
优选地，在获取相邻帧彩色图像的匹配特征点时，只有在相邻帧彩色图像中的匹配特征点存在于当前帧彩色图像之前的预定数量的帧中的情况下，才将该匹配特征点确定为最终的匹配特征点，以提高特征点匹配的精度。
对于上述轨迹χi，可以设置一个最小的轨迹长度来去除错误的匹配点。作为具体示例，只有轨迹χi的长度为预定数量的情况下，才将该匹配特征点确定为最终的匹配特征点，以提高特征点匹配的精度。其中，可以根据经验确定上述预定数量。
在深度图像特征点获取单元606中，可以获取与当前帧彩色图像中的所述至少一个特征点中的每个的位置分别对应的当前帧深度图像中的至少一个特征点，并且可以基于相邻帧彩色图像中的所述至少一个相应匹配特征点，获取与当前帧深度图像中的所述至少一个特征点中的每个分别匹配的相邻帧深度图像中的至少一个相应匹配特征点。
优选地，当前帧深度图像的所述至少一个特征点和相邻帧深度图像的所述至少一个相应匹配特征点是通过如下方式获取的：可以将当前帧彩色图像中的所述至少一个特征点映射到当前帧深度图像中，以获取当前帧深度图像中的所述至少一个特征点；以及可以将相邻帧彩色图像中的所述至少一个特征点映射到相邻帧深度图像中，以获取相邻帧深度图像中的所述至少一个特征点。
获取当前帧深度图像的所述至少一个特征点和相邻帧深度图像的所述至少一个相应匹配特征点的示例可以参见以上方法实施例中相应位置的描述，在此不再重复。
在世界坐标表示单元608中，可以基于当前摄像机姿态以及用于校正图像中的特征点的非刚性畸变的参数来表示当前帧深度图像中的每个特征点的世界坐标。
可以将用于校正图像中的特征点的非刚性畸变的参数反映在特征点的世界坐标下中。也就是说，可以基于第k帧对应的摄像机姿态以及用于校正图像中的特征点的非刚性畸变的参数来表示第k帧深度图像中的每个特征点的世界坐标。
优选地，所述参数包括用于校正在当前帧彩色图像中获取匹配的特征点时产生的匹配畸变的参数、用于校正当前帧彩色图像和当前帧深度图像不对齐所产生的不对齐畸变的参数以及用于校正在将当前帧深度图像中的特征点在摄像机坐标中投影时产生的投影畸变的参数中的至少一个。
作为具体示例，在摄像机姿态估计系统中，可能存在匹配畸变、不对齐畸变以及投影畸变这三种非刚性畸变。关于匹配畸变、不对齐畸变以及投影畸变可以参见以上方法实施例中相应位置的描述，在此不再重复。
优选地，用于校正匹配畸变的参数是当前帧彩色图像中的多个预定位置处的点的偏移向量，其中，通过对表示当前帧彩色图像中的每个特征点的位置的向量叠加所述偏移向量的线性组合，对该特征点的位置进行校正。
关于对匹配畸变进行校正可以参见以上方法实施例中相应位置的描述，在此不再重复。
优选地，用于校正不对齐畸变的参数是当前帧深度图像中的多个预定位置处的点的偏移向量，其中，通过对表示当前帧深度图像中的每个特征点的位置的向量叠加所述偏移向量的线性组合，对该特征点的位置进行校正。
关于对不对齐畸变进行校正可以参见以上方法实施例中相应位置的描述，在此不再重复。
优选地，用于校正投影畸变的参数是当前帧深度图像中的多个预定位置处的点的摄像机坐标的偏移向量，其中，通过对表示当前帧深度图像中的每个特征点的摄像机坐标的向量叠加所述偏移向量的线性组合，对该特征点的摄像机坐标进行校正。
关于对投影畸变进行校正可以参见以上方法实施例中相应位置的描述，在此不再重复。
关于将上述三种畸变的非刚性校正引入到摄像机姿态估计中的方法可以参见以上方法实施例中相应位置的描述，在此不再重复。
在摄像机姿态和参数估计单元610中，可以基于当前帧深度图像中的每个特征点的世界坐标和相邻帧深度图像中的相应匹配特征点的世界坐标来计算该对匹配的特征点之间的距离，以及以针对每对匹配的特征点所计算出的距离之和最小为目标，估计当前摄像机姿态和所述参数。
优选地，可以基于当前帧深度图像中的每个特征点的世界坐标和相邻帧深度图像中的相应匹配特征点的世界坐标，计算该特征点与其相应匹配特征点之间的距离。即，在世界坐标系下，计算当前深度图像中的每个特征点和相邻帧深度图像中的相应匹配特征点之间的距离。
优选地，可以采用交替优化的方式来分别估计当前摄像机姿态和所述参数。
计算当前深度图像中的每个特征点的世界坐标和相邻帧深度图像中的相应匹配特征点之间的距离以及估计当前摄像机姿和所述参数的示例可以参见以上方法实施例中相应位置的描述，在此不再重复。
优选地，可以将经非刚性畸变校正后的当前帧深度图像作为下一帧深度图像的相邻帧深度图像。作为具体示例，可以将经非刚性畸变校正后的第k帧深度图像作为第k+1帧深度图像的相邻帧深度图像，这样，可以进一步提供摄像机姿态参数估计的精度。
优选地，对于一个图像序列，将第一帧深度图像的摄像机坐标系设置为世界坐标系；以及摄像机姿态为用于将摄像机坐标系变换到世界坐标系的变换矩阵。具体示例可以参见以上方法实施例中相应位置的描述，在此不再重复。
优选地，在估计出当前摄像机姿态的情况下，基于所估计的当前摄像机姿态来计算当前帧深度图像中的所有特征点的世界坐标；以及对于当前帧深度图像中的每个特征点：如果该特征点也存在于当前帧深度图像之前的深度图像中，则将该特征点在当前帧深度图像中的世界坐标与该特征点在当前帧深度图像之前的深度图像中的世界坐标进行平均，作为该特征点在当前帧深度图像中的更新后的世界坐标。
作为具体示例，在估计出第k帧图像所对应的摄像机姿态的情况下，基于所估计的摄像机姿态来计算第k帧深度图像中的所有特征点的世界坐标；以及对于第k帧深度图像中的每个特征点：如果该特征点也存在于第k帧深度图像之前的深度图像中，则将该特征点在第k帧深度图像中的世界坐标与该特征点在第k帧深度图像之前的深度图像中的世界坐标进行平均，作为该特征点在第k帧深度图像中的更新后的世界坐标。换句话说，如果某个特征点所在的特征跟踪轨迹在第k帧之前已经求解出世界坐标系下的三维点坐标，则将它们取平均作为对应轨迹的新的世界坐标系下的三维点坐标。
在对一个图像序列进行处理时，这种渐进式的摄像机姿态求解方式容易造成时序上的误差累积，从而导致摄像机姿态的估计严重偏离正确的取值。为此，可以采取全局优化的方式来降低累积误差。全局优化即同时优化多帧的摄像机姿态。
优选地，根据本公开实施例的摄像机姿态估计装置还可以包括优化单元，其可以被配置成：将深度图像的序列划分成组；以及对于每个组，基于与该组中的每帧深度图像的所有特征点和该组中的所有其他帧深度图像的所有相应匹配特征点有关的组几何约束，优化该组中的每帧深度图像对应的摄像机姿态和所述参数。
优选地，所述组几何约束为：使得在世界坐标系下的、该组中的所有对匹配特征点之间的距离之和最小。
全局优化的示例可以参见以上方法实施例中相应位置的描述，在此不再重复。
全局优化可以在任意时刻进行，但该算法复杂度较高，每帧都进行全局优化是不必要的。可以通过监控当前帧上的误差值，当误差值超过某个阈值时，才采用全局优化的策略。
综上所述，根据本公开的实施例的摄像机姿态估计装置600结合视觉信息和深度信息，并且可以校正存在于深度图像、彩色图像以及三维空间中的非刚性畸变，从而能够实现摄像机参数的精确估计。
应指出，尽管以上描述了根据本公开的实施例的摄像机姿态估计装置的功能配置，但是这仅是示例而非限制，并且本领域技术人员可根据本公开的原理对以上实施例进行修改，例如可对各个实施例中的功能模块进行添加、删除或者组合等，并且这样的修改均落入本公开的范围内。
此外，还应指出，这里的装置实施例是与上述方法实施例相对应的，因此在装置实施例中未详细描述的内容可参见方法实施例中相应位置的描述，在此不再重复描述。
应理解，根据本公开的实施例的存储介质和程序产品中的机器可执行的指令还可以被配置成执行上述摄像机姿态估计方法，因此在此未详细描述的内容可参考先前相应位置的描述，在此不再重复进行描述。
相应地，用于承载上述包括机器可执行的指令的程序产品的存储介质也包括在本发明的公开中。该存储介质包括但不限于软盘、光盘、磁光盘、存储卡、存储棒等等。
另外，还应该指出的是，上述系列处理和装置也可以通过软件和/或固件实现。在通过软件和/或固件实现的情况下，从存储介质或网络向具有专用硬件结构的计算机，例如图7所示的通用个人计算机700安装构成该软件的程序，该计算机在安装有各种程序时，能够执行各种功能等等。
在图7中，中央处理单元(CPU)701根据只读存储器(ROM)702中存储的程序或从存储部分708加载到随机存取存储器(RAM)703的程序执行各种处理。在RAM 703中，也根据需要存储当CPU 701执行各种处理等时所需的数据。
CPU 701、ROM 702和RAM 703经由总线704彼此连接。输入/输出接口705也连接到总线704。
下述部件连接到输入/输出接口705：输入部分706，包括键盘、鼠标等；输出部分707，包括显示器，比如阴极射线管(CRT)、液晶显示器(LCD)等，和扬声器等；存储部分708，包括硬盘等；和通信部分709，包括网络接口卡比如LAN卡、调制解调器等。通信部分709经由网络比如因特网执行通信处理。
根据需要，驱动器710也连接到输入/输出接口705。可拆卸介质711比如磁盘、光盘、磁光盘、半导体存储器等等根据需要被安装在驱动器710上，使得从中读出的计算机程序根据需要被安装到存储部分708中。
在通过软件实现上述系列处理的情况下，从网络比如因特网或存储介质比如可拆卸介质711安装构成软件的程序。
本领域的技术人员应当理解，这种存储介质不局限于图7所示的其中存储有程序、与设备相分离地分发以向用户提供程序的可拆卸介质711。可拆卸介质711的例子包含磁盘(包含软盘(注册商标))、光盘(包含光盘只读存储器(CD-ROM)和数字通用盘(DVD))、磁光盘(包含迷你盘(MD)(注册商标))和半导体存储器。或者，存储介质可以是ROM 702、存储部分708中包含的硬盘等等，其中存有程序，并且与包含它们的设备一起被分发给用户。
以上参照附图描述了本公开的优选实施例，但是本公开当然不限于以上示例。本领域技术人员可在所附权利要求的范围内得到各种变更和修改，并且应理解这些变更和修改自然将落入本公开的技术范围内。
例如，在以上实施例中包括在一个单元中的多个功能可以由分开的装置来实现。替选地，在以上实施例中由多个单元实现的多个功能可分别由分开的装置来实现。另外，以上功能之一可由多个单元来实现。无需说，这样的配置包括在本公开的技术范围内。
在该说明书中，流程图中所描述的步骤不仅包括以所述顺序按时间序列执行的处理，而且包括并行地或单独地而不是必须按时间序列执行的处理。此外，甚至在按时间序列处理的步骤中，无需说，也可以适当地改变该顺序。
另外，根据本公开的技术还可以如下进行配置。
附记1.一种摄像机姿态估计方法，包括：
通过摄像机获得彩色图像和与所述彩色图像对应的深度图像；
获取当前帧彩色图像中的至少一个特征点和与所述至少一个特征点中的每个分别匹配的相邻帧彩色图像中的至少一个相应匹配特征点；
获取与所述当前帧彩色图像中的所述至少一个特征点中的每个的位置分别对应的当前帧深度图像中的至少一个特征点，并且基于所述相邻帧彩色图像中的所述至少一个相应匹配特征点，获取与所述当前帧深度图像中的所述至少一个特征点中的每个分别匹配的相邻帧深度图像中的至少一个相应匹配特征点；
基于当前摄像机姿态以及用于校正图像中的特征点的非刚性畸变的参数来表示所述当前帧深度图像中的每个特征点的世界坐标；以及
基于所述当前帧深度图像中的每个特征点的世界坐标和所述相邻帧深度图像中的相应匹配特征点的世界坐标来计算该对匹配的特征点之间的距离，以及以针对每对匹配的特征点所计算出的距离之和最小为目标，估计所述当前摄像机姿态和所述参数。
附记2.根据附记1所述的摄像机姿态估计方法，其中，所述参数包括用于校正在所述当前帧彩色图像中获取匹配的特征点时产生的匹配畸变的参数、用于校正所述当前帧彩色图像和所述当前帧深度图像不对齐所产生的不对齐畸变的参数以及用于校正在将所述当前帧深度图像中的特征点在摄像机坐标中投影时产生的投影畸变的参数中的至少一个。
附记3.根据附记2所述的摄像机姿态估计方法，其中，所述用于校正匹配畸变的参数是所述当前帧彩色图像中的多个预定位置处的点的偏移向量，其中，通过对表示所述当前帧彩色图像中的每个特征点的位置的向量叠加所述偏移向量的线性组合，对该特征点的位置进行校正。
附记4.根据附记2所述的摄像机姿态估计方法，其中，所述用于校正不对齐畸变的参数是所述当前帧深度图像中的多个预定位置处的点的偏移向量，其中，通过对表示所述当前帧深度图像中的每个特征点的位置的向量叠加所述偏移向量的线性组合，对该特征点的位置进行校正。
附记5.根据附记2所述的摄像机姿态估计方法，其中，所述用于校正投影畸变的参数是所述当前帧深度图像中的多个预定位置处的点的摄像机坐标的偏移向量，其中，通过对表示所述当前帧深度图像中的每个特征点的摄像机坐标的向量叠加所述偏移向量的线性组合，对该特征点的摄像机坐标进行校正。
附记6.根据附记1所述的摄像机姿态估计方法，其中，采用交替优化的方式来分别估计所述当前摄像机姿态和所述参数。
附记7.根据附记1所述的摄像机姿态估计方法，
其中，在估计出所述当前摄像机姿态的情况下，基于所估计的当前摄像机姿态来计算所述当前帧深度图像中的所有特征点的世界坐标；以及
对于所述当前帧深度图像中的每个特征点：
如果该特征点也存在于所述当前帧深度图像之前的深度图像中，则将该特征点在所述当前帧深度图像中的世界坐标与该特征点在所述当前帧深度图像之前的深度图像中的世界坐标进行平均，作为该特征点在所述当前帧深度图像中的更新后的世界坐标。
附记8.根据附记1所述的摄像机姿态估计方法，还包括：
将深度图像的序列划分成组；以及
对于每个组，基于与该组中的每帧深度图像的所有特征点和该组中的所有其他帧深度图像的所有相应匹配特征点有关的组几何约束，优化该组中的每帧深度图像对应的摄像机姿态和所述参数。
附记9.根据附记8所述的摄像机姿态估计方法，其中，所述组几何约束为：
使得在世界坐标系下的、该组中的所有对匹配的特征点之间的距离之和最小。
附记10.根据附记1所述的摄像机姿态估计方法，其中，将经非刚性畸变校正后的所述当前帧深度图像作为下一帧深度图像的相邻帧深度图像。
附记11.一种摄像机姿态估计装置，包括：
彩色图像和深度图像获取单元，被配置成通过摄像机获得彩色图像和与所述彩色图像对应的深度图像；
彩色图像特征点获取单元，被配置成获取当前帧彩色图像中的至少一个特征点和与所述至少一个特征点中的每个分别匹配的相邻帧彩色图像中的至少一个相应匹配特征点；
深度图像特征点获取单元，被配置成获取与所述当前帧彩色图像中的所述至少一个特征点中的每个的位置分别对应的当前帧深度图像中的至少一个特征点，并且基于所述相邻帧彩色图像中的所述至少一个相应匹配特征点，获取与所述当前帧深度图像中的所述至少一个特征点中的每个分别匹配的相邻帧深度图像中的至少一个相应匹配特征点；
世界坐标表示单元，被配置成基于当前摄像机姿态以及用于校正图像中的特征点的非刚性畸变的参数来表示所述当前帧深度图像中的每个特征点的世界坐标；以及
摄像机姿态和参数估计单元，被配置成基于所述当前帧深度图像中的每个特征点的世界坐标和所述相邻帧深度图像中的相应匹配特征点的世界坐标来计算该对匹配的特征点之间的距离，以及以针对每对匹配的特征点所计算出的距离之和最小为目标，估计所述当前摄像机姿态和所述参数。
附记12.根据附记11所述的摄像机姿态估计装置，其中，所述参数包括用于校正在所述当前帧彩色图像中获取匹配的特征点时产生的匹配畸变的参数、用于校正所述当前帧彩色图像和所述当前帧深度图像不对齐所产生的不对齐畸变的参数以及用于校正在将所述当前帧深度图像中的特征点在摄像机坐标中投影时产生的投影畸变的参数中的至少一个。
附记13.根据附记12所述的摄像机姿态估计装置，其中，所述用于校正匹配畸变的参数是所述当前帧彩色图像中的多个预定位置处的点的偏移向量，其中，通过对表示所述当前帧彩色图像中的每个特征点的位置的向量叠加所述偏移向量的线性组合，对该特征点的位置进行校正。
附记14.根据附记12所述的摄像机姿态估计装置，其中，所述用于校正不对齐畸变的参数是所述当前帧深度图像中的多个预定位置处的点的偏移向量，其中，通过对表示所述当前帧深度图像中的每个特征点的位置的向量叠加所述偏移向量的线性组合，对该特征点的位置进行校正。
附记15.根据附记12所述的摄像机姿态估计装置，其中，所述用于校正投影畸变的参数是所述当前帧深度图像中的多个预定位置处的点的摄像机坐标的偏移向量，其中，通过对表示所述当前帧深度图像中的每个特征点的摄像机坐标的向量叠加所述偏移向量的线性组合，对该特征点的摄像机坐标进行校正。
附记16.根据附记11所述的摄像机姿态估计装置，其中，采用交替优化的方式来分别估计所述当前摄像机姿态和所述参数。
附记17.根据附记11所述的摄像机姿态估计装置，
其中，在估计出所述当前摄像机姿态的情况下，基于所估计的当前摄像机姿态来计算所述当前帧深度图像中的所有特征点的世界坐标；以及
对于所述当前帧深度图像中的每个特征点：
如果该特征点也存在于所述当前帧深度图像之前的深度图像中，则将该特征点在所述当前帧深度图像中的世界坐标与该特征点在所述当前帧深度图像之前的深度图像中的世界坐标进行平均，作为该特征点在所述当前帧深度图像中的更新后的世界坐标。
附记18.根据附记11所述的摄像机姿态估计装置，还包括优化单元，所述优化单元被配置成：
将深度图像的序列划分成组；以及
对于每个组，基于与该组中的每帧深度图像的所有特征点和该组中的所有其他帧深度图像的所有相应匹配特征点有关的组几何约束，优化该组中的每帧深度图像对应的摄像机姿态和所述参数。
附记19.根据附记18所述的摄像机姿态估计装置，其中，所述组几何约束为：
使得在世界坐标系下的、该组中的所有对匹配的特征点之间的距离之和最小。
附记20.根据附记11所述的摄像机姿态估计装置，其中，将经非刚性畸变校正后的所述当前帧深度图像作为下一帧深度图像的相邻帧深度图像。

Claims (10)

1.一种摄像机姿态估计方法，包括：
通过摄像机获得彩色图像和与所述彩色图像对应的深度图像；
获取当前帧彩色图像中的至少一个特征点和与所述至少一个特征点中的每个分别匹配的相邻帧彩色图像中的至少一个相应匹配特征点；
获取与所述当前帧彩色图像中的所述至少一个特征点中的每个的位置分别对应的当前帧深度图像中的至少一个特征点，并且基于所述相邻帧彩色图像中的所述至少一个相应匹配特征点，获取与所述当前帧深度图像中的所述至少一个特征点中的每个分别匹配的相邻帧深度图像中的至少一个相应匹配特征点；
基于当前摄像机姿态以及用于校正图像中的特征点的非刚性畸变的参数来表示所述当前帧深度图像中的每个特征点的世界坐标；以及
基于所述当前帧深度图像中的每个特征点的世界坐标和所述相邻帧深度图像中的相应匹配特征点的世界坐标来计算所述当前帧深度图像中的每个特征点与所述相邻帧深度图像中的相应匹配特征点之间的距离，以及以针对每对匹配的特征点所计算出的距离之和最小为目标，估计所述当前摄像机姿态和所述参数。
2.根据权利要求1所述的摄像机姿态估计方法，其中，所述参数包括用于校正在所述当前帧彩色图像中获取匹配的特征点时产生的匹配畸变的参数、用于校正所述当前帧彩色图像和所述当前帧深度图像不对齐所产生的不对齐畸变的参数以及用于校正在将所述当前帧深度图像中的特征点在摄像机坐标中投影时产生的投影畸变的参数中的至少一个。
3.根据权利要求2所述的摄像机姿态估计方法，其中，所述用于校正匹配畸变的参数是所述当前帧彩色图像中的多个预定位置处的点的偏移向量，其中，通过对表示所述当前帧彩色图像中的每个特征点的位置的向量叠加所述偏移向量的线性组合，对该特征点的位置进行校正。
4.根据权利要求2所述的摄像机姿态估计方法，其中，所述用于校正不对齐畸变的参数是所述当前帧深度图像中的多个预定位置处的点的偏移向量，其中，通过对表示所述当前帧深度图像中的每个特征点的位置的向量叠加所述偏移向量的线性组合，对该特征点的位置进行校正。
5.根据权利要求2所述的摄像机姿态估计方法，其中，所述用于校正投影畸变的参数是所述当前帧深度图像中的多个预定位置处的点的摄像机坐标的偏移向量，其中，通过对表示所述当前帧深度图像中的每个特征点的摄像机坐标的向量叠加所述偏移向量的线性组合，对该特征点的摄像机坐标进行校正。
6.根据权利要求1所述的摄像机姿态估计方法，其中，采用交替优化的方式来分别估计所述当前摄像机姿态和所述参数。
7.根据权利要求1所述的摄像机姿态估计方法，
其中，在估计出所述当前摄像机姿态的情况下，基于所估计的当前摄像机姿态来计算所述当前帧深度图像中的所有特征点的世界坐标；以及
对于所述当前帧深度图像中的每个特征点：
如果该特征点也存在于所述当前帧深度图像之前的深度图像中，则将该特征点在所述当前帧深度图像中的世界坐标与该特征点在所述当前帧深度图像之前的深度图像中的世界坐标进行平均，作为该特征点在所述当前帧深度图像中的更新后的世界坐标。
8.根据权利要求1所述的摄像机姿态估计方法，还包括：
将深度图像的序列划分成组；以及
对于每个组，基于与该组中的每帧深度图像的所有特征点和该组中的所有其他帧深度图像的所有相应匹配特征点有关的组几何约束，优化该组中的每帧深度图像对应的摄像机姿态和所述参数。
9.根据权利要求8所述的摄像机姿态估计方法，其中，所述组几何约束为：
使得在世界坐标系下的、该组中的所有对匹配的特征点之间的距离之和最小。
10.一种摄像机姿态估计装置，包括：
彩色图像和深度图像获取单元，被配置成通过摄像机获得彩色图像和与所述彩色图像对应的深度图像；
彩色图像特征点获取单元，被配置成获取当前帧彩色图像中的至少一个特征点和与所述至少一个特征点中的每个分别匹配的相邻帧彩色图像中的至少一个相应匹配特征点；
深度图像特征点获取单元，被配置成获取与所述当前帧彩色图像中的所述至少一个特征点中的每个的位置分别对应的当前帧深度图像中的至少一个特征点，并且基于所述相邻帧彩色图像中的所述至少一个相应匹配特征点，获取与所述当前帧深度图像中的所述至少一个特征点中的每个分别匹配的相邻帧深度图像中的至少一个相应匹配特征点；
世界坐标表示单元，被配置成基于当前摄像机姿态以及用于校正图像中的特征点的非刚性畸变的参数来表示所述当前帧深度图像中的每个特征点的世界坐标；以及
摄像机姿态和参数估计单元，被配置成基于所述当前帧深度图像中的每个特征点的世界坐标和所述相邻帧深度图像中的相应匹配特征点的世界坐标来计算所述当前帧深度图像中的每个特征点与所述相邻帧深度图像中的相应匹配特征点之间的距离，以及以针对每对匹配的特征点所计算出的距离之和最小为目标，估计所述当前摄像机姿态和所述参数。
CN201710044456.1A 2017-01-19 2017-01-19 摄像机姿态估计方法和摄像机姿态估计装置 Active CN108335328B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710044456.1A CN108335328B (zh) 2017-01-19 2017-01-19 摄像机姿态估计方法和摄像机姿态估计装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710044456.1A CN108335328B (zh) 2017-01-19 2017-01-19 摄像机姿态估计方法和摄像机姿态估计装置

Publications (2)

Publication Number Publication Date
CN108335328A CN108335328A (zh) 2018-07-27
CN108335328B true CN108335328B (zh) 2021-09-24

Family

ID=62922785

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710044456.1A Active CN108335328B (zh) 2017-01-19 2017-01-19 摄像机姿态估计方法和摄像机姿态估计装置

Country Status (1)

Country Link
CN (1) CN108335328B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110132987A (zh) * 2019-06-10 2019-08-16 重庆工商大学 基于数字近景摄影测量岩体节理迹长和密度的测量方法
CN111712857A (zh) * 2019-06-25 2020-09-25 深圳市大疆创新科技有限公司 图像处理方法、装置、云台和存储介质
CN110361005B (zh) * 2019-06-26 2021-03-26 达闼机器人有限公司 定位方法、定位装置、可读存储介质及电子设备
CN113409391B (zh) * 2021-06-25 2023-03-03 浙江商汤科技开发有限公司 视觉定位方法及相关装置、设备和存储介质

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101021947A (zh) * 2006-09-22 2007-08-22 东南大学 三维扫描系统中双摄像机标定方法
WO2006099597A3 (en) * 2005-03-17 2007-11-01 Honda Motor Co Ltd Pose estimation based on critical point analysis
CN104346829A (zh) * 2013-07-29 2015-02-11 中国农业机械化科学研究院 基于pmd相机和摄像头的彩色三维重建系统及其方法
CN104835159A (zh) * 2015-05-07 2015-08-12 中国科学院长春光学精密机械与物理研究所 用于连续变焦距光学成像系统的数字图像校正方法
CN106157367A (zh) * 2015-03-23 2016-11-23 联想(北京)有限公司 三维场景重建方法和设备

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101398937B (zh) * 2008-10-29 2011-05-18 北京航空航天大学 基于同一场景散乱照片集的三维重构方法
CN101815174B (zh) * 2010-01-11 2015-03-04 北京中星微电子有限公司 摄像控制方法及控制装置
CN102855470B (zh) * 2012-07-31 2015-04-08 中国科学院自动化研究所 基于深度图像的人体姿态估计方法
US9251590B2 (en) * 2013-01-24 2016-02-02 Microsoft Technology Licensing, Llc Camera pose estimation for 3D reconstruction
KR102209008B1 (ko) * 2014-02-17 2021-01-28 삼성전자주식회사 카메라 포즈 추정 장치 및 카메라 포즈 추정 방법

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2006099597A3 (en) * 2005-03-17 2007-11-01 Honda Motor Co Ltd Pose estimation based on critical point analysis
CN101021947A (zh) * 2006-09-22 2007-08-22 东南大学 三维扫描系统中双摄像机标定方法
CN104346829A (zh) * 2013-07-29 2015-02-11 中国农业机械化科学研究院 基于pmd相机和摄像头的彩色三维重建系统及其方法
CN106157367A (zh) * 2015-03-23 2016-11-23 联想(北京)有限公司 三维场景重建方法和设备
CN104835159A (zh) * 2015-05-07 2015-08-12 中国科学院长春光学精密机械与物理研究所 用于连续变焦距光学成像系统的数字图像校正方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Rough compressed domain camera pose estimation through object motion;Christian K.等;《2009 16th IEEE International Conference on Image Processing (ICIP)》;20100217;第3481-3484页 *
一种摄像机姿态估计方法;周维 等;《系统仿真学报》;20111231;第23卷(第12期);第2582-2587页 *
应用抛物线运动估计摄像机姿态;吕耀文 等;《光学精密工程》;20140430;第22卷(第4期);第1078-1085页 *

Also Published As

Publication number Publication date
CN108335328A (zh) 2018-07-27

Similar Documents

Publication Publication Date Title
CN108335327B (zh) 摄像机姿态估计方法和摄像机姿态估计装置
CN108027878B (zh) 用于面部对齐的方法
US8953847B2 (en) Method and apparatus for solving position and orientation from correlated point features in images
JP6573354B2 (ja) 画像処理装置、画像処理方法、及びプログラム
CN108335328B (zh) 摄像机姿态估计方法和摄像机姿态估计装置
CN101398934B (zh) 对图像中的对象进行定位的方法和系统
Barath et al. Homography from two orientation-and scale-covariant features
KR101926563B1 (ko) 카메라 추적을 위한 방법 및 장치
US9247139B2 (en) Method for video background subtraction using factorized matrix completion
KR100951309B1 (ko) 광학식 모션 캡처 장비를 위한 다중 카메라 보정 방법
CN109376641B (zh) 一种基于无人机航拍视频的运动车辆检测方法
US10701336B2 (en) Rectifying a sequence of stereo images
US20130208009A1 (en) Method and apparatus for optimization and incremental improvement of a fundamental matrix
US11651581B2 (en) System and method for correspondence map determination
Eichhardt et al. Affine correspondences between central cameras for rapid relative pose estimation
JP2017130067A (ja) 衛星映像の位置正確度改善のための自動映像処理システム及びその方法
KR100792172B1 (ko) 강건한 대응점을 이용한 기본행렬 추정 장치 및 그 방법
CN115330876B (zh) 基于孪生网络和中心位置估计的目标模板图匹配定位方法
Wu et al. Direct pose estimation for planar objects
CN114399547B (zh) 一种基于多帧的单目slam鲁棒初始化方法
CN106651950B (zh) 一种基于二次曲线透视投影不变性的单相机位姿估计方法
Zhang et al. Iterative fitting after elastic registration: An efficient strategy for accurate estimation of parametric deformations
EP2245593A2 (en) A method of estimating a motion of a multiple camera system, a multiple camera system and a computer program product
JP5215615B2 (ja) ３次元位置情報復元装置およびその方法
Duan et al. A Simple and Efficient Method for Radial Distortion Estimation by Relative Orientation

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