CN107256563A - 基于差异液位图像序列的水下三维重建系统及其方法 - Google Patents

基于差异液位图像序列的水下三维重建系统及其方法 Download PDF

Info

Publication number
CN107256563A
CN107256563A CN201710441431.5A CN201710441431A CN107256563A CN 107256563 A CN107256563 A CN 107256563A CN 201710441431 A CN201710441431 A CN 201710441431A CN 107256563 A CN107256563 A CN 107256563A
Authority
CN
China
Prior art keywords
mtd
mrow
msub
mtr
msubsup
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
CN201710441431.5A
Other languages
English (en)
Other versions
CN107256563B (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.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
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 National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN201710441431.5A priority Critical patent/CN107256563B/zh
Publication of CN107256563A publication Critical patent/CN107256563A/zh
Application granted granted Critical
Publication of CN107256563B publication Critical patent/CN107256563B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • 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
    • G06T7/55—Depth or shape recovery from multiple images
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
    • 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/10028—Range image; Depth image; 3D point clouds

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Computer Graphics (AREA)
  • Geometry (AREA)
  • Software Systems (AREA)
  • Image Processing (AREA)
  • Image Analysis (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)

Abstract

本发明公开了一种基于差异液位图像序列的水下三维重建系统及其方法。该系统包括：摄像机、计算机系统、主水箱、液位差异测量装置、支架，该系统实施方法包括以下步骤：首先，采用单摄像机拍摄差异液位图像序列；其次，进行图像局部特征提取、匹配及误匹配剔除；再次，估计水平面在相机坐标系下的法向量方向，并计算图像特征点转换坐标；然后，对差异液位图像序列中各图像对应的水面距离进行标定；最后，进行稠密三维场景重建和优化。本发明仅需要一个照相机进行数据采集，且无需额外水下标识物或特殊标定设备，相对于一般水下三维重建方法和系统具有操作便捷、成本低廉、重建快速的优势，适用于任何折射系数已知的透明液体环境。

Description

基于差异液位图像序列的水下三维重建系统及其方法
技术领域
本发明属于计算机视觉技术领域，特别涉及一种基于差异液位图像序列的水下三维重建系统及其方法。
背景技术
目前，绝大多数基于图像的三维重建方法和系统均针对陆上环境，高精度的水下三维重建方法直到最近几年才逐渐在计算机视觉领域引起重视。基于图像的水下三维重建面临的主要挑战是：拍摄水下环境的过程中，光线需要穿过水、空气等不同的介质，由于不同介质的折射率(refractive index)不同，光线会在介质交界面产生折射现象，从而导致成像折射变形(refractive distortion)。普通图像三维重建方法采用透视相机模型，由于该模型假设光线沿直线传播，因此透视相机模型在有折射现象发生的成像系统中不再有效。
一种简单的水下三维重建方法是忽略折射现象，或使用近似模型来补偿折射变形，例如：焦距调节模型、镜头径向畸变模型，以及兼具焦距调节和镜头畸变的近似模型。然而，由于折射变形的高度非线性和场景结构依赖性，忽略折射或采用近似模型来补偿折射理论上都无法消除折射对三维重建的影响。由于折射变形依赖于物体的深度，从理论上来说无法使用任何图像空间(image space)变形模型来描述，使用透视相机模型将导致较大的相机定标误差和三维重建误差。
另一种更加精确的方法是利用符合光学原理的正确的折射相机模型(RefractiveCameraModel)来显式的对折射建模。Chari和Sturm(Chari V,Sturm P.Multiple-viewgeometry of the refractive plane[C].In Proceedings of BritishMachine VisionConference(BMVC).2009:1–11.)分析了位于单个折射介质同侧的两个相机之间的几何关系，并从理论上证实了折射基本矩阵、折射单应矩阵等几何量的存在。然而，由于折射相机模型下的几何量估计仍然是一个极具挑战的课题，因此上述文献并没有给出上述理论结果的实际应用。Chang和Chen(Chang Y,Chen T.Multi-view 3D reconstruction for scenesunder the refractive plane withknown vertical direction[C].In Proceedings ofIEEE International Conference on ComputerVision(ICCV).2011:351–358.)提出一个水下三维重建方法，该方法适用于单折射界面(refractive interface)的情形。在该方法中，折射现象被显式的表达为场景深度的函数，因此是一个精确的折射相机模型。该水下三维重建方法方法的局限性在于该方法是建立在已知相机两个旋转参数的基础上的，这两个参数需要诸如惯性测量单元(Inertial Measurement Unit，IMU)的外部硬件获得；此外，该方法还假设折射界面的法向量已知。Agrawal等(Agrawal A,Ramalingam S,Taguchi Y,etal.A theory of multi-layer flat refractive geometry[C].In Proceedings of IEEEConference on Computer Vision and Pattern Recognition(CVPR).2012:3346–3353.)的研究表明折射几何在本质上与轴向相机(Axial Camera)相同，并提出了一个利用平面定标物体对这种系统进行定标的方法。Sedlazeck和Koch(Sedlazeck A,Koch R.Calibrationof housing parameters forunderwater stereo-camera rigs[C].In Proceedings ofBritish Machine Vision Conference(BMVC).2011:1–11.)提出一种水下双目立体相机外壳参数标定方法，该方法的主要思想是最小化三维场景点在外部折射界面的虚拟透视投影误差，该方法的主要缺点是其优化过程非常耗时，在普通PC上的运行时间长达3小时之久。
综上所述，已有基于图像的水下三维重建方法的主要局限是需要特定的标定物或者额外的标定设备才能实现相机的标定，或者需要借助耗时的优化计算实现相机自定标，因此水下三维重建方法在系统成本、操作便捷性和计算效率等方面均有较大的提升空间。
发明内容
本发明的目的在于克服上述现有技术的不足，提出一种基于差异液位图像序列的水下三维重建系统及其方法。
为实现上述目的，一方面本发明提出一种如下技术方案：
本发明提出的一种基于差异液位图像序列的水下三维重建系统包括：摄像机、计算机系统、主水箱、液位差异测量装置；
所述摄像机用于获取水下场景的彩色照片，摄像机的内部参数已经事先经过标定；
所述计算机系统用于控制摄像机拍照以及执行基于差异液位图像序列的水下三维重建相关计算和结果展示；
所述液位差异测量装置用于测量拍摄不同图像时的液位差异；
所述主水箱用于盛放水及待重建物体，主水箱顶部敞开。
作为所述水下三维重建系统技术方案的进一步改进，所述水下三维重建系统还包括水量控制装置和固定支架；所述水量控制装置用于控制主水箱的水量，所述固定支架包括摄像机固定支架和主水箱固定支架，所述摄像机固定支架用于将摄像机固定在主水箱上方，摄像机与计算机系统之间通过数据线连接。
作为所述水下三维重建系统技术方案的进一步改进，所述计算机系统包括PC机、显示器及外接设备；所述水量控制装置包括进水阀、泄水阀、水管；所述液位差异测量装置具体采用液位尺读取液位高度并计算不同液位的差异。
为实现上述目的，另一方面本发明提出一种如下技术方案：
本发明提出的一种基于差异液位图像序列的水下三维重建方法包括前面所述一种基于差异液位图像序列的水下三维重建系统，该水下三维重建方法具体包括以下步骤：
S1.单摄像机拍摄差异液位图像序列；
S2.图像局部特征提取、匹配及误匹配剔除；
S3.估计水平面在相机坐标系下的法向量方向，并计算图像特征点转换坐标；
S4.差异液位图像序列中各图像对应的水面距离标定；
S5.稠密三维场景重建及其优化。
作为所述水下三维重建方法技术方案的进一步改进，所述步骤S1具体包括以下内容：
在不改变摄像机内部参数和外部参数的情况下，通过多次拍摄静态水下场景获取一个水下图像序列；每次拍摄前通过减少或增加容器内的水量调节液位高度，确保拍摄各图像时所对应的水面距离依次递增或递减；
其中，水面距离指的是相机投影中心与水面的垂直距离，记图像的数量为M(M为整数，且M≥3)，各图像所对应的水面距离记为{Dl|l＝1,2,…,M}，假定各水面距离满足如下约束：D1＜D2＜…＜DM。
作为所述水下三维重建方法技术方案的进一步改进，所述步骤S2具体包括以下内容：
S2.1根据相机的内部参数对图像进行变形矫正，移除由于镜头生产工艺限制导致的图像径向变形，获得矫正水下图像序列
S2.2利用仿射不变性图像局部特征提取和匹配方法计算中两两图像的特征点匹配结果，在此基础上获取M幅图像的匹配特征轨迹，每条特征轨迹由匹配成功的多个特征点坐标组成，并删除特征点数量少于M的特征轨迹；
S2.3对于每条特征轨迹，利用最小二乘法从其特征点坐标中拟合一条直线，若该特征轨迹中各特征点到所拟合的直线的距离均小于1个像素，则保留该特征轨迹，否则删除该特征轨迹，记最终保留下来的特征轨迹为{Lp|p＝1,2,…,N}；
作为所述水下三维重建方法技术方案的进一步改进，所述步骤S3具体包括以下内容：
S3.1从N条特征轨迹拟合直线中随机挑选2条直线并计算其交点，上述过程重复执行P次，记所有交点为{ci＝(xi,yi)T|i＝1,2,…,P}，其中，若N＞150，则P取11325，否则取P＝(N2-N)/2；
根据如下公式计算所有拟合直线的理想交汇点坐标c＝(xc,yc)T：
S3.2对各特征轨迹中的每个特征点v＝(x,y)T，根据如下关系计算各特征点转换坐标v′＝(x′,y′)T：
其中，符号[·]i表示取向量的第i个分量，且：
上式中f为照相机的焦距长度(单位为像素)，θx、θy的计算公式如下：
作为所述水下三维重建方法技术方案的进一步改进，所述步骤S4具体包括以下内容：
S4.1初始化参数Niter＝1、Nin＝0，集合
S4.2从{Lp|p＝1,2,…,N}中随机挑选1条特征轨迹，根据该特征轨迹中各特征点的转换坐标，计算各图像所对应的未知水面距离Dl(l＝{1,2,…M})和该条特征轨迹对应的场景点的深度D；
其中，场景点的深度指的是相机投影中心到场景点所在水平面的垂直距离；
水面距离和场景点深度的最优估计通过求解如下优化问题计算：
其中，符号||·||2表示取向量的二范数，液位差异D1-Dj(j＝2,3,…,M)为已知常量，0α和0β分别为(M-2)×1和(M-1)×1维零向量，1为所有元素均为1的(M-2)×1维向量，I为(M-2)×(M-2)维单位矩阵，Mα＝(w(v′1)-w(v′2),w(v′1)-w(v′3),…,w(v′1)-w(v′M))T，Mβ为对角线上第j(j＝{1,2,…M-1})个元素取值为(w(v′j+1)-a(v′j+1))的(M-1)×(M-1)维对角矩阵；
其中，函数a(·)和w(·)的定义为：
上式中n为液体的已知折射系数，v′l＝(x′l,y′l)T为当前所选中的特征轨迹第l(l＝{1,2,…M})个特征点的转换坐标；
S4.3设定参数N′in＝0，根据S4.2计算的为每条特征轨迹计算一组深度估计值计算公式如下：
然后，根据如下公式计算特征轨迹的发散指标τ：
如果τ＜0.05，则更新参数N′in←N′in+1；在将所有特征轨迹处理完成后，如果N′in＞Nin，则令Nin＝N′in，且更新并将集合设置为包含所有满足τ＜0.05的特征轨迹的集合，记
S4.4更新Niter←Niter+1，如果则转至步骤S4.2；否则，执行步骤S4.5；
其中，参数δ的取值范围0.01～0.02，ε的取值范围是0.30～0.45；
S4.5通过求解如下优化问题获得各水面距离和各特征轨迹对应的场景点深度的最终估计
其中，表示特征轨迹在图像中的特征点转换坐标，0γ为(MNin-Nin)×1维零向量， 为对角线上第j(j＝{1,2,…M-1})个元素取值为的(M-1)×(M-1)维对角矩阵；
最后，更新变量
作为所述水下三维重建方法技术方案的进一步改进，所述步骤S5具体包括以下内容：
S5.1利用图像稠密匹配方法计算图像与图像之间的像素对应关系，记图像中任意像素r＝(xr,yr)T在中的对应像素位置为s＝(xs,ys)T；并利用S3.2中的转换坐标计算方法计算r、s的转换坐标r′＝(x′r,y′r)T、s′＝(x′s,y′s)T；
S5.2对于中每个像素位置r，根据如下公式计算其对应的场景点深度：
然后，根据d(r)估计像素r对应的场景点的初始三维坐标Xr＝(Xr,Yr,Zr)T，计算公式如下：
S5.3对于每个场景点，利用局部非线性优化方法对其初始三维坐标求精，优化过程的目标函数为：
其中，表示图像中以Xr在该图像上的投影为中心、大小为m×m个像素的图像窗口，与图像中以Xr在该图像上的投影为中心、大小为m×m个像素的图像窗口的归一化互相关系数(normalized cross correlation,NCC)；
其中，m的取值范围是6～10个像素，Xr在图像上的投影点的计算基于直接前向投影(forward projection)方法，相机的旋转矩阵根据θx和θy确定；
S5.4将所有经过优化的三维点融合为三维点云，并利用统计外点移除(statistical outlier removal,SOR)点云过滤方法剔除孤立点和细小点云片段，获得最终稠密三维重建结果。
与现有技术相比，本发明具有以下有益效果：
本发明公开了一种基于差异液位图像序列的水下三维重建系统及其方法，其属于一种非接触、被动式水下三维重建系统及方法，本发明仅需要一个照相机进行数据采集，且无需额外水下标识物或特殊标定设备，充分利用不同液位图像的折射差异，相对于一般水下三维重建方法和系统具有操作便捷、成本低廉、重建快速的优势，且本发明系统及其方法适用于任何折射系数已知的透明液体环境。
附图说明
图1为基于差异液位图像序列的水下三维重建系统结构示意图。
图2为基于差异液位图像序列的水下三维重建方法的流程图。
图3为差异液位图像成像几何关系示意图。
图4为特征轨迹拟合直线交汇点估计示例和误差分析示意图。
图5为差异液位图像水面距离标定流程图。
图6为差异液位图像水面距离标定的误差分析示意图。
图7为差异液位图像样张和三维重建示例图。
具体实施方式
本发明为一种基于差异液位图像序列的水下三维重建系统及其方法，在不借助额外水下标识物或特殊标定设备的情况下，充分利用不同液位图像的折射差异，从单个摄像机拍摄的水下图像序列中计算水下场景的三维几何结构，实现水下物体的非接触、被动式三维模型重建。
下面将结合本申请说明书附图，对本发明的一种基于差异液位图像序列的水下三维重建系统及其方法的具体实施例做进一步详细说明，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。
实施例1
实施例1为一种基于差异液位图像序列的水下三维重建系统，该水下三维重建系统结构具体如图1所示，该系统主要包括摄像机、计算机系统、主水箱、水量控制装置、液位差异测量装置、支架等部分。其中，摄像机用于获取水下场景的彩色照片，摄像机的内部参数已经事先经过标定；计算机系统用于控制摄像机拍照以及执行基于差异液位图像序列的水下三维重建相关计算和结果展示，包括PC机、显示器及其他外接设备；主水箱用于盛放水及待重建物体，主水箱顶部敞开；水量控制装置用于控制主水箱的水量，包括进水阀、泄水阀、水管等；液位差异测量装置用于测量拍摄不同图像时的液位差异，采用液位尺读取液位高度并计算不同液位的差异；主水箱固定支架设置在主水箱的底部，其用于固定主水箱。
该系统中，摄像机通过摄像机固定支架固定在主水箱上方，摄像机与计算机系统之间通过数据线连接；摄像机镜头朝向待重建物体，摄像机能够观察到整个待重建物体范围；待重建物体放置在主水箱底部，拍摄照片时水面需为静止状态；拍摄差异液位图像序列过程中水深保持足够的深度，以便图像中具有较明显的折射效应，拍摄各图像时对应的液位具备足够的差异。
实施例2
实施例2为一种基于差异液位图像序列的水下三维重建方法，其包括了实施例1中的一种基于差异液位图像序列的水下三维重建系统，该水下三维重建方法的具体步骤流程如图2所示，包括如下步骤：
S1.单摄像机拍摄差异液位图像序列；
利用一个固定在水面上方的摄像机拍摄一个差异液位图像序列。具体方法是：在不改变摄像机内部参数和外部参数的情况下，通过多次拍摄静态水下场景获取一个水下图像序列。每次拍摄前通过减少(或增加)容器内的水量调节液位高度，确保拍摄各图像时所对应的水面距离依次递增(或递减)。其中，水面距离指的是摄相机投影中心与水面的垂直距离。记拍摄的图像的数量为M(M为整数，且M≥3)，各图像所对应的水面距离记为{Dl|l＝1,2,…,M}。不是一般性，假定各水面距离满足如下约束：D1＜D2＜…＜DM。
S2.图像局部特征提取、匹配及误匹配剔除；
S2.1根据相机的内部参数对拍摄的图像进行变形矫正，移除由于镜头生产工艺限制导致的图像径向变形，获得矫正水下图像序列：
图像径向变形矫正采用Hartley R和Zisserman A(Hartley R,ZissermanA.Multiple view geometry in computer vision[M].2nd ed.Cambridge UniversityPress,2004:189-193)提出的方法。
S2.2利用ASIFT特征提取和匹配方法计算中两两图像的特征点匹配结果，在此基础上获取M幅图像的匹配特征轨迹；每条特征轨迹由匹配成功的多个特征点坐标组成，并删除特征点数量少于M的特征轨迹。
S2.3图3为差异液位图像成像几何关系示意图，理想情况下同一个场景点在不同液位条件下的成像点位于同一条直线上，因此可以利用上述特性剔除误匹配。对于每条特征轨迹，利用最小二乘法从其特征点坐标中拟合一条直线，若该特征轨迹中各特征点到所拟合的直线的距离均小于1个像素，则保留该特征轨迹，否则删除该特征轨迹；记最终保留下来的特征轨迹为{Lp|p＝1,2,…,N}。
S3.估计水平面在相机坐标系下的法向量方向，并计算图像特征点转换坐标；
S3.1从N条特征轨迹拟合直线中随机挑选2条直线并计算其交点，上述过程重复执行P次，记所有交点为{ci＝(xi,yi)T|i＝1,2,…,P}。其中，若N＞150，则P取11325；否则取P＝(N2-N)/2。如图3所示，理想情况下所有特征轨迹拟合直线均相交于公共点，例如：由确定的直线与由确定的直线相交于点c。但由于实际图像中存在噪声，上述性质并不成立。根据如下公式计算所有拟合直线的理想交汇点坐标c＝(xc,yc)T：
图4(a)为从10000个候选交点中估计理想交汇点的实例，该示例中摄相机主光轴偏离水面法线方向10度，特征点坐标噪声标准差为0.1个像素。图4(b)为交汇点估计误差与坐标噪声的关系，每个箱形图(box plot)展示相应噪声条件下交汇点估计误差的中值、上、下四分位点以及外点信息，上述统计信息基于两幅差异液位图像的60次统计结果。交汇点估计误差定义为∠cocgt,cgt为交汇点坐标的真实值(ground truth)。
S3.2对各特征轨迹中的每个特征点v＝(x,y)T，根据如下关系计算各特征点转换坐标v′＝(x′,y′)T：
其中，符号[·]i表示取向量的第i个分量，且：
上式中f为照相机的焦距长度(单位为像素)，θx、θy的计算公式如下：
S4.差异液位图像序列中各图像对应的水面距离标定，其中，差异液位图像水面距离标定流程如图5所示，具体包括如下步骤：
S4.1初始化参数Niter＝1、Nin＝0，集合
S4.2从{Lp|p＝1,2,…,N}中随机挑选1条特征轨迹，根据该特征轨迹中各特征点的转换坐标，计算各图像所对应的未知水面距离Dl(l＝{1,2,…M})和该条特征轨迹对应的场景点的深度D。其中，场景点的深度指的是相机投影中心到场景点所在水平面的垂直距离。水面距离和场景点深度的最优估计通过求解如下优化问题计算：
其中，符号||·||2表示取向量的二范数。液位差异D1-Dj(j＝2,3,…,M)为已知常量。0α和0β分别为(M-2)×1和(M-1)×1维零向量，1为所有元素均为1的(M-2)×1维向量，I为(M-2)×(M-2)维单位矩阵，Mα＝(w(v′1)-w(v′2),w(v′1)-w(v′3),…,w(v′1)-w(v′M))T，Mβ为对角线上第j(j＝{1,2,…M-1})个元素取值为(w(v′j+1)-a(v′j+1))的(M-1)×(M-1)维对角矩阵。其中，函数a(·)和w(·)的定义为：
上式中n为液体的已知折射系数，v′l＝(x′l,y′l)T为当前所选中的特征轨迹第l(l＝{1,2,…M})个特征点的转换坐标。
S4.3设定参数N′in＝0。根据S4.2计算的为每条特征轨迹计算一组深度估计值计算公式如下：
然后，根据如下公式计算特征轨迹的发散指标τ：
如果τ＜0.05，则更新参数N′in←N′in+1。按照上述方法将所有特征轨迹处理完成后，如果N′in＞Nin，则令Nin＝N′in，且更新(l＝1,2,…,M)，并将集合设置为包含所有满足τ＜0.05的特征轨迹的集合，记
S4.4更新Niter←Niter+1。如果则转至S4.2；否则，执行S4.5。其中，参数δ的取值是0.015，ε的取值是0.40。
S4.5通过求解如下优化问题获得各水面距离和各特征轨迹对应的场景点深度的最终估计
其中，表示特征轨迹在图像中的特征点转换坐标，0γ为(MNin-Nin)×1维零向量。 为对角线上第j(j＝{1,2,…M-1})个元素取值为的(M-1)×(M-1)维对角矩阵。最后，更新变量
差异液位图像水面距离标定的误差分析结果如图6所示，其中图6(a)为差异液位层数对标定精度的影响，图6(b)为液位差异测量误差对标定精度的影响。上述结果为Nin＝100、特征点坐标噪声标准差为1.0个像素的条件下的60次统计结果，标定误差定义为所有距离的相对误差的均方根(root mean square,RMS)。可见，增加差异液位图像有助于提高标定精度，而液位差异测量误差对标定精度的影响在可接受范围之内。例如：液位差异测量误差在5.5％时，标定误差的中值仅为2％左右。
S5.稠密三维场景重建及其优化；
S5.1利用已有图像稠密匹配方法计算图像与图像之间的像素对应关系，记图像中任意像素r＝(xr,yr)T在中的对应像素位置为s＝(xs,ys)T。利用S3.2中的转换坐标计算方法计算r、s的转换坐标r′＝(x′r,y′r)T、s′＝(x′s,y′s)T。
S5.2对于中每个像素位置r，根据如下公式计算其对应的场景点深度：
然后，根据d(r)估计像素r对应的场景点的初始三维坐标Xr＝(Xr,Yr,Zr)T，计算公式如下：
S5.3对于每个场景点，利用局部非线性优化方法对其初始三维坐标求精，优化过程的目标函数为：
其中，表示图像中以Xr在该图像上的投影为中心、大小为m×m个像素的图像窗口，与图像中以Xr在该图像上的投影为中心、大小为m×m个像素的图像窗口的归一化互相关系数(normalized cross correlation,NCC)。其中，m的取值范围是6个像素，Xr在图像上的投影点的计算基于已有直接前向投影(forward projection)方法，相机的旋转矩阵根据θx和θy确定。
S5.4将所有经过优化的三维点融合为三维点云，并利用已有统计外点移除(statistical outlier removal,SOR)点云过滤方法剔除孤立点和细小点云片段，获得最终稠密三维重建结果。
本具体实施例中差异液位图像样张和三维重建示例如图7所示，其中，图7(a)为不同液位条件下获取的图像样张，图7(b)分别为采用3幅、6幅、10幅差异液位图像的重建结果。根据图7的测试结果表明，增加差异液位图像的数量有助于提高三维重建精度，该结果与图6(a)中的规律是相吻合的。
本领域技术人员将清楚本发明的范围不限制于以上讨论的示例，有可能对其进行若干改变和修改，而不脱离所附权利要求书限定的本发明的范围。尽管己经在附图和说明书中详细图示和描述了本发明，但这样的说明和描述仅是说明或示意性的，而非限制性的。本发明并不限于所公开的实施例。

Claims (9)

1.一种基于差异液位图像序列的水下三维重建系统，其特征在于，所述水下三维重建系统包括：摄像机、计算机系统、主水箱、液位差异测量装置；
所述摄像机用于获取水下场景的彩色照片，摄像机的内部参数已经事先经过标定；
所述计算机系统用于控制摄像机拍照以及执行基于差异液位图像序列的水下三维重建相关计算和结果展示；
所述液位差异测量装置用于测量拍摄不同图像时的液位差异；
所述主水箱用于盛放水及待重建物体，主水箱顶部敞开。
2.根据权利要求1所述的一种基于差异液位图像序列的水下三维重建系统，其特征在于，所述水下三维重建系统还包括水量控制装置和固定支架；所述水量控制装置用于控制主水箱的水量，所述固定支架包括摄像机固定支架和主水箱固定支架，所述摄像机固定支架用于将摄像机固定在主水箱上方，摄像机与计算机系统之间通过数据线连接。
3.根据权利要求1所述的一种基于差异液位图像序列的水下三维重建系统，其特征在于，所述计算机系统包括PC机、显示器及外接设备；所述水量控制装置包括进水阀、泄水阀、水管；所述液位差异测量装置具体采用液位尺读取液位高度并计算不同液位的差异。
4.一种基于差异液位图像序列的水下三维重建方法，其特征在于，该水下三维重建方法包括权利要求1-3任一项所述的一种基于差异液位图像序列的水下三维重建系统，该水下三维重建方法具体包括以下步骤：
S1.单摄像机拍摄差异液位图像序列；
S2.图像局部特征提取、匹配及误匹配剔除；
S3.估计水平面在相机坐标系下的法向量方向，并计算图像特征点转换坐标；
S4.差异液位图像序列中各图像对应的水面距离标定；
S5.稠密三维场景重建及其优化。
5.根据权利要求4所述的一种基于差异液位图像序列的水下三维重建方法，其特征在于，所述步骤S1具体包括以下内容：
在不改变摄像机内部参数和外部参数的情况下，通过多次拍摄静态水下场景获取一个水下图像序列，每次拍摄前通过减少或增加容器内的水量调节液位高度，确保拍摄各图像时所对应的水面距离依次递增或递减；
其中，水面距离指的是摄相机投影中心与水面的垂直距离，记拍摄的图像的数量为M(M为整数，且M≥3)，各图像所对应的水面距离记为{Dl|l＝1,2,…,M}，假定各水面距离满足如下约束：D1＜D2＜…＜DM。
6.根据权利要求5所述的一种基于差异液位图像序列的水下三维重建方法，其特征在于，所述步骤S2具体包括以下内容：
S2.1根据相机的内部参数对图像进行变形矫正，移除由于镜头生产工艺限制导致的图像径向变形，获得矫正水下图像序列
S2.2利用仿射不变性图像局部特征提取和匹配方法计算中两两图像的特征点匹配结果，在此基础上获取M幅图像的匹配特征轨迹，每条特征轨迹由匹配成功的多个特征点坐标组成，并删除特征点数量少于M的特征轨迹；
S2.3对于每条特征轨迹，利用最小二乘法从其特征点坐标中拟合一条直线，若该特征轨迹中各特征点到所拟合的直线的距离均小于1个像素，则保留该特征轨迹，否则删除该特征轨迹，记最终保留下来的特征轨迹为{Lp|p＝1,2,…,N}。
7.根据权利要求6所述的一种基于差异液位图像序列的水下三维重建方法，其特征在于，所述步骤S3具体包括以下内容：
S3.1从N条特征轨迹拟合直线中随机挑选2条直线并计算其交点，上述过程重复执行P次，记所有交点为{ci＝(xi,yi)T|i＝1,2,…,P}，其中，若N＞150，则P取11325，否则取P＝(N2-N)/2；
根据如下公式计算所有拟合直线的理想交汇点坐标c＝(xc,yc)T：
<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>x</mi> <mi>c</mi> </msub> <mo>=</mo> <mi>m</mi> <mi>e</mi> <mi>d</mi> <mo>{</mo> <msub> <mi>x</mi> <mn>1</mn> </msub> <mo>,</mo> <msub> <mi>x</mi> <mn>2</mn> </msub> <mo>,</mo> <mo>...</mo> <mo>,</mo> <msub> <mi>x</mi> <mi>P</mi> </msub> <mo>}</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>y</mi> <mi>c</mi> </msub> <mo>=</mo> <mi>m</mi> <mi>e</mi> <mi>d</mi> <mo>{</mo> <msub> <mi>y</mi> <mn>1</mn> </msub> <mo>,</mo> <msub> <mi>y</mi> <mn>2</mn> </msub> <mo>,</mo> <mo>...</mo> <mo>,</mo> <msub> <mi>y</mi> <mi>P</mi> </msub> <mo>}</mo> </mrow> </mtd> </mtr> </mtable> </mfenced>
S3.2对各特征轨迹中的每个特征点v＝(x,y)T，根据如下关系计算各特征点转换坐标v′＝(x′,y′)T：
<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msup> <mi>x</mi> <mo>&amp;prime;</mo> </msup> <mo>=</mo> <mfrac> <msub> <mrow> <mo>&amp;lsqb;</mo> <mi>H</mi> <mi>v</mi> <mo>&amp;rsqb;</mo> </mrow> <mn>1</mn> </msub> <msub> <mrow> <mo>&amp;lsqb;</mo> <mi>H</mi> <mi>v</mi> <mo>&amp;rsqb;</mo> </mrow> <mn>3</mn> </msub> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msup> <mi>y</mi> <mo>&amp;prime;</mo> </msup> <mo>=</mo> <mfrac> <msub> <mrow> <mo>&amp;lsqb;</mo> <mi>H</mi> <mi>v</mi> <mo>&amp;rsqb;</mo> </mrow> <mn>2</mn> </msub> <msub> <mrow> <mo>&amp;lsqb;</mo> <mi>H</mi> <mi>v</mi> <mo>&amp;rsqb;</mo> </mrow> <mn>3</mn> </msub> </mfrac> </mrow> </mtd> </mtr> </mtable> </mfenced>
其中，符号[·]i表示取向量的第i个分量，且：
<mrow> <mi>H</mi> <mo>=</mo> <mfenced open = "(" close = ")"> <mtable> <mtr> <mtd> <mi>f</mi> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mi>f</mi> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> </mtr> </mtable> </mfenced> <mfenced open = "(" close = ")"> <mtable> <mtr> <mtd> <mrow> <msub> <mi>cos&amp;theta;</mi> <mi>y</mi> </msub> </mrow> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mrow> <msub> <mi>sin&amp;theta;</mi> <mi>y</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>-</mo> <msub> <mi>sin&amp;theta;</mi> <mi>y</mi> </msub> </mrow> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mrow> <msub> <mi>cos&amp;theta;</mi> <mi>y</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mfenced open = "(" close = ")"> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mrow> <msub> <mi>cos&amp;theta;</mi> <mi>x</mi> </msub> </mrow> </mtd> <mtd> <mrow> <mo>-</mo> <msub> <mi>sin&amp;theta;</mi> <mi>x</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mrow> <msub> <mi>sin&amp;theta;</mi> <mi>x</mi> </msub> </mrow> </mtd> <mtd> <mrow> <msub> <mi>cos&amp;theta;</mi> <mi>x</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <msup> <mfenced open = "(" close = ")"> <mtable> <mtr> <mtd> <mi>f</mi> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mi>f</mi> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>1</mn> </mtd> </mtr> </mtable> </mfenced> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> </mrow>
上式中f为照相机的焦距长度(单位为像素)，θx、θy的计算公式如下：
<mrow> <mo>{</mo> <mrow> <mtable> <mtr> <mtd> <mrow> <msub> <mi>&amp;theta;</mi> <mi>x</mi> </msub> <mo>=</mo> <mi>a</mi> <mi>r</mi> <mi>c</mi> <mi>t</mi> <mi>a</mi> <mi>n</mi> <mfrac> <msub> <mi>y</mi> <mi>c</mi> </msub> <mi>f</mi> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>&amp;theta;</mi> <mi>y</mi> </msub> <mo>=</mo> <mo>-</mo> <mi>a</mi> <mi>r</mi> <mi>c</mi> <mi>t</mi> <mi>a</mi> <mi>n</mi> <mfrac> <msub> <mi>x</mi> <mi>c</mi> </msub> <msqrt> <mrow> <msup> <mi>f</mi> <mn>2</mn> </msup> <mo>+</mo> <msubsup> <mi>y</mi> <mi>c</mi> <mn>2</mn> </msubsup> </mrow> </msqrt> </mfrac> </mrow> </mtd> </mtr> </mtable> <mo>.</mo> </mrow> </mrow>
8.根据权利要求7所述的一种基于差异液位图像序列的水下三维重建方法，其特征在于，所述步骤S4具体包括以下内容：
S4.1初始化参数Niter＝1、Nin＝0，集合
S4.2从{Lp|p＝1,2,…,N}中随机挑选1条特征轨迹，根据该特征轨迹中各特征点的转换坐标，计算各图像所对应的未知水面距离Dl(l＝{1,2,…M})和该条特征轨迹对应的场景点的深度D；
其中，场景点的深度D指的是相机投影中心到场景点所在水平面的垂直距离；
水面距离和场景点深度的最优估计通过求解如下优化问题计算：
<mrow> <msup> <mrow> <mo>(</mo> <msup> <mi>D</mi> <mo>*</mo> </msup> <mo>,</mo> <msubsup> <mi>D</mi> <mn>1</mn> <mo>*</mo> </msubsup> <mo>,</mo> <msubsup> <mi>D</mi> <mn>2</mn> <mo>*</mo> </msubsup> <mo>,</mo> <mn>...</mn> <mo>,</mo> <msubsup> <mi>D</mi> <mi>M</mi> <mo>*</mo> </msubsup> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mo>=</mo> <munder> <mrow> <mi>arg</mi> <mi>min</mi> </mrow> <msup> <mrow> <mo>(</mo> <mi>D</mi> <mo>,</mo> <msub> <mi>D</mi> <mn>1</mn> </msub> <mo>,</mo> <msub> <mi>D</mi> <mn>2</mn> </msub> <mo>,</mo> <mn>...</mn> <mo>,</mo> <msub> <mi>D</mi> <mi>M</mi> </msub> <mo>)</mo> </mrow> <mi>T</mi> </msup> </munder> <msub> <mrow> <mo>||</mo> <mrow> <mfenced open = "(" close = ")"> <mtable> <mtr> <mtd> <msub> <mi>M</mi> <mi>&amp;alpha;</mi> </msub> </mtd> <mtd> <mrow> <mo>(</mo> <mi>a</mi> <mo>(</mo> <msubsup> <mi>v</mi> <mn>1</mn> <mo>&amp;prime;</mo> </msubsup> <mo>)</mo> <mo>-</mo> <mi>w</mi> <mo>(</mo> <msubsup> <mi>v</mi> <mn>1</mn> <mo>&amp;prime;</mo> </msubsup> <mo>)</mo> <mo>)</mo> <mfenced open = "(" close = ")"> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> </mtr> <mtr> <mtd> <mn>1</mn> </mtd> </mtr> </mtable> </mfenced> </mrow> </mtd> <mtd> <msub> <mi>M</mi> <mi>&amp;beta;</mi> </msub> </mtd> </mtr> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </mtd> <mtd> <msubsup> <mn>0</mn> <mi>&amp;alpha;</mi> <mi>T</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <msub> <mn>0</mn> <mi>&amp;alpha;</mi> </msub> </mtd> <mtd> <mrow> <mo>-</mo> <mi>I</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> <mfenced open = "(" close = ")"> <mtable> <mtr> <mtd> <mi>D</mi> </mtd> </mtr> <mtr> <mtd> <msub> <mi>D</mi> <mn>1</mn> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>D</mi> <mn>2</mn> </msub> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <msub> <mi>D</mi> <mi>M</mi> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mfenced open = "(" close = ")"> <mtable> <mtr> <mtd> <msub> <mn>0</mn> <mi>&amp;beta;</mi> </msub> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>D</mi> <mn>1</mn> </msub> <mo>-</mo> <msub> <mi>D</mi> <mn>2</mn> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>D</mi> <mn>1</mn> </msub> <mo>-</mo> <msub> <mi>D</mi> <mn>3</mn> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>D</mi> <mn>1</mn> </msub> <mo>-</mo> <msub> <mi>D</mi> <mi>M</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow> <mo>||</mo> </mrow> <mn>2</mn> </msub> </mrow>
其中，符号||·||2表示取向量的二范数，液位差异D1-Dj(j＝2,3,…,M)为已知常量，0α和0β分别为(M-2)×1和(M-1)×1维零向量，1为所有元素均为1的(M-2)×1维向量，I为(M-2)×(M-2)维单位矩阵，Mα＝(w(v′1)-w(v′2),w(v′1)-w(v′3),…,w(v′1)-w(v′M))T，Mβ为对角线上第j(j＝{1,2,…M-1})个元素取值为(w(v′j+1)-a(v′j+1))的(M-1)×(M-1)维对角矩阵；
其中，函数a(·)和w(·)的定义为：
<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>a</mi> <mrow> <mo>(</mo> <msubsup> <mi>v</mi> <mi>l</mi> <mo>&amp;prime;</mo> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <msqrt> <mrow> <msubsup> <mi>x</mi> <mi>l</mi> <mrow> <mo>&amp;prime;</mo> <mn>2</mn> </mrow> </msubsup> <mo>+</mo> <msubsup> <mi>y</mi> <mi>l</mi> <mrow> <mo>&amp;prime;</mo> <mn>2</mn> </mrow> </msubsup> </mrow> </msqrt> <mi>f</mi> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>w</mi> <mrow> <mo>(</mo> <msubsup> <mi>v</mi> <mi>l</mi> <mo>&amp;prime;</mo> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mi>a</mi> <mrow> <mo>(</mo> <msubsup> <mi>v</mi> <mi>l</mi> <mo>&amp;prime;</mo> </msubsup> <mo>)</mo> </mrow> </mrow> <msqrt> <mrow> <mo>(</mo> <msup> <mi>n</mi> <mn>2</mn> </msup> <mo>-</mo> <mn>1</mn> <mo>)</mo> <msup> <mrow> <mo>(</mo> <mi>a</mi> <mo>(</mo> <msubsup> <mi>v</mi> <mi>l</mi> <mo>&amp;prime;</mo> </msubsup> <mo>)</mo> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>n</mi> <mn>2</mn> </msup> </mrow> </msqrt> </mfrac> </mrow> </mtd> </mtr> </mtable> </mfenced>
上式中n为液体的已知折射系数，v′l＝(x′l,y′l)T为当前所选中的特征轨迹第l(l＝{1,2,…M})个特征点的转换坐标；
S4.3设定参数N′in＝0，根据S4.2计算的为每条特征轨迹计算一组深度估计值计算公式如下：
<mrow> <msub> <mover> <mi>d</mi> <mo>^</mo> </mover> <mi>k</mi> </msub> <mo>=</mo> <mfrac> <mrow> <mo>(</mo> <mi>w</mi> <mo>(</mo> <msubsup> <mi>v</mi> <mn>1</mn> <mo>&amp;prime;</mo> </msubsup> <mo>)</mo> <mo>-</mo> <mi>a</mi> <mo>(</mo> <msubsup> <mi>v</mi> <mn>1</mn> <mo>&amp;prime;</mo> </msubsup> <mo>)</mo> <mo>)</mo> <msubsup> <mi>D</mi> <mn>1</mn> <mo>*</mo> </msubsup> <mo>+</mo> <mo>(</mo> <mi>a</mi> <mo>(</mo> <msubsup> <mi>v</mi> <mn>1</mn> <mo>&amp;prime;</mo> </msubsup> <mo>)</mo> <mo>-</mo> <mi>w</mi> <mo>(</mo> <msubsup> <mi>v</mi> <mn>1</mn> <mo>&amp;prime;</mo> </msubsup> <mo>)</mo> <mo>)</mo> <msubsup> <mi>D</mi> <mi>l</mi> <mo>*</mo> </msubsup> </mrow> <mrow> <mi>w</mi> <mrow> <mo>(</mo> <msubsup> <mi>v</mi> <mn>1</mn> <mo>&amp;prime;</mo> </msubsup> <mo>)</mo> </mrow> <mo>-</mo> <mi>w</mi> <mrow> <mo>(</mo> <msubsup> <mi>v</mi> <mn>1</mn> <mo>&amp;prime;</mo> </msubsup> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow>
然后，根据如下公式计算特征轨迹的发散指标τ：
<mrow> <mi>&amp;tau;</mi> <mo>=</mo> <mfrac> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> <mo>{</mo> <msub> <mover> <mi>d</mi> <mo>^</mo> </mover> <mi>k</mi> </msub> <mo>}</mo> <mo>-</mo> <mi>m</mi> <mi>i</mi> <mi>n</mi> <mo>{</mo> <msub> <mover> <mi>d</mi> <mo>^</mo> </mover> <mi>k</mi> </msub> <mo>}</mo> </mrow> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> <mo>{</mo> <msub> <mover> <mi>d</mi> <mo>^</mo> </mover> <mi>k</mi> </msub> <mo>}</mo> </mrow> </mfrac> </mrow>
如果τ＜0.05，则更新参数N′in←N′in+1；在将所有特征轨迹处理完成后，如果N′in＞Nin，则令Nin＝N′in，且更新并将集合设置为包含所有满足τ＜0.05的特征轨迹的集合，记
S4.4更新Niter←Niter+1，如果则转至步骤S4.2；否则，执行步骤S4.5；
其中，参数δ的取值是0.015，ε的取值是0.40；
S4.5通过求解如下优化问题获得各水面距离和各特征轨迹对应的场景点深度的最终估计
<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <msup> <mrow> <mo>(</mo> <msubsup> <mi>D</mi> <msub> <mi>p</mi> <mn>1</mn> </msub> <mo>*</mo> </msubsup> <mo>,</mo> <msubsup> <mi>D</mi> <msub> <mi>p</mi> <mn>2</mn> </msub> <mo>*</mo> </msubsup> <mo>,</mo> <mn>...</mn> <mo>,</mo> <msubsup> <mi>D</mi> <msub> <mi>p</mi> <msub> <mi>N</mi> <mrow> <mi>i</mi> <mi>n</mi> </mrow> </msub> </msub> <mo>*</mo> </msubsup> <mo>,</mo> <msubsup> <mi>D</mi> <mn>1</mn> <mo>*</mo> </msubsup> <mo>,</mo> <msubsup> <mi>D</mi> <mn>2</mn> <mo>*</mo> </msubsup> <mo>,</mo> <mn>...</mn> <mo>,</mo> <msubsup> <mi>D</mi> <mi>M</mi> <mo>*</mo> </msubsup> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mo>=</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <munder> <mrow> <mi>arg</mi> <mi>min</mi> </mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>D</mi> <msub> <mi>p</mi> <mn>1</mn> </msub> </msub> <mo>,</mo> <msub> <mi>D</mi> <msub> <mi>p</mi> <mn>2</mn> </msub> </msub> <mo>,</mo> <mn>...</mn> <mo>,</mo> <msub> <mi>D</mi> <mrow> <msub> <mi>pN</mi> <mrow> <mi>i</mi> <mi>n</mi> </mrow> </msub> </mrow> </msub> <mo>,</mo> <msub> <mi>D</mi> <mn>1</mn> </msub> <mo>,</mo> <msub> <mi>D</mi> <mn>2</mn> </msub> <mo>,</mo> <mn>...</mn> <mo>,</mo> <msub> <mi>D</mi> <mi>M</mi> </msub> <mo>)</mo> </mrow> <mi>T</mi> </msup> </munder> <msub> <mrow> <mo>||</mo> <mrow> <mfenced open = "(" close = ")"> <mtable> <mtr> <mtd> <msubsup> <mi>M</mi> <mi>&amp;alpha;</mi> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </msubsup> </mtd> <mtd> <mrow> <mo>(</mo> <mi>a</mi> <mo>(</mo> <msubsup> <mi>v</mi> <mn>1</mn> <mrow> <mo>&amp;prime;</mo> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> <mo>-</mo> <mi>w</mi> <mo>(</mo> <msubsup> <mi>v</mi> <mn>1</mn> <mrow> <mo>&amp;prime;</mo> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> <mo>)</mo> <mfenced open = "(" close = ")"> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> </mtr> <mtr> <mtd> <mn>1</mn> </mtd> </mtr> </mtable> </mfenced> </mrow> </mtd> <mtd> <msubsup> <mi>M</mi> <mi>&amp;beta;</mi> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </msubsup> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>M</mi> <mi>&amp;alpha;</mi> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </msubsup> </mtd> <mtd> <mrow> <mo>(</mo> <mi>a</mi> <mo>(</mo> <msubsup> <mi>v</mi> <mn>1</mn> <mrow> <mo>&amp;prime;</mo> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> <mo>-</mo> <mi>w</mi> <mo>(</mo> <msubsup> <mi>v</mi> <mn>1</mn> <mrow> <mo>&amp;prime;</mo> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> <mo>)</mo> <mfenced open = "(" close = ")"> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> </mtr> <mtr> <mtd> <mn>1</mn> </mtd> </mtr> </mtable> </mfenced> </mrow> </mtd> <mtd> <msubsup> <mi>M</mi> <mi>&amp;beta;</mi> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </msubsup> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <msubsup> <mi>M</mi> <mi>&amp;alpha;</mi> <mrow> <mo>(</mo> <msub> <mi>p</mi> <msub> <mi>N</mi> <mrow> <mi>i</mi> <mi>n</mi> </mrow> </msub> </msub> <mo>)</mo> </mrow> </msubsup> </mtd> <mtd> <mrow> <mo>(</mo> <mi>a</mi> <mo>(</mo> <msubsup> <mi>v</mi> <mn>1</mn> <mrow> <mo>&amp;prime;</mo> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> <mo>-</mo> <mi>w</mi> <mo>(</mo> <msubsup> <mi>v</mi> <mn>1</mn> <mrow> <mo>&amp;prime;</mo> <mrow> <mo>(</mo> <msub> <mi>p</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> <mo>)</mo> <mfenced open = "(" close = ")"> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> </mtr> <mtr> <mtd> <mn>1</mn> </mtd> </mtr> </mtable> </mfenced> </mrow> </mtd> <mtd> <msubsup> <mi>M</mi> <mi>&amp;beta;</mi> <mrow> <mo>(</mo> <msub> <mi>p</mi> <msub> <mi>N</mi> <mrow> <mi>i</mi> <mi>n</mi> </mrow> </msub> </msub> <mo>)</mo> </mrow> </msubsup> </mtd> </mtr> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </mtd> <mtd> <msubsup> <mn>0</mn> <mi>&amp;alpha;</mi> <mi>T</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <msub> <mn>0</mn> <mi>&amp;alpha;</mi> </msub> </mtd> <mtd> <mrow> <mo>-</mo> <mi>I</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> <mfenced open = "(" close = ")"> <mtable> <mtr> <mtd> <msub> <mi>D</mi> <msub> <mi>p</mi> <mn>1</mn> </msub> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>D</mi> <msub> <mi>p</mi> <mn>2</mn> </msub> </msub> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <msub> <mi>D</mi> <msub> <mi>p</mi> <msub> <mi>N</mi> <mrow> <mi>i</mi> <mi>n</mi> </mrow> </msub> </msub> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>D</mi> <mn>1</mn> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>D</mi> <mn>2</mn> </msub> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <msub> <mi>D</mi> <mi>M</mi> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mfenced open = "(" close = ")"> <mtable> <mtr> <mtd> <msub> <mn>0</mn> <mi>&amp;gamma;</mi> </msub> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>D</mi> <mn>1</mn> </msub> <mo>-</mo> <msub> <mi>D</mi> <mn>2</mn> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>D</mi> <mn>1</mn> </msub> <mo>-</mo> <msub> <mi>D</mi> <mn>3</mn> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>D</mi> <mn>1</mn> </msub> <mo>-</mo> <msub> <mi>D</mi> <mi>M</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow> <mo>||</mo> </mrow> <mn>2</mn> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced>
其中，表示特征轨迹在图像中的特征点转换坐标，0γ为(MNin-Nin)×1维零向量， 为对角线上第j(j＝{1,2,…M-1})个元素取值为的(M-1)×(M-1)维对角矩阵；
最后，更新变量
9.根据权利要求8所述的一种基于差异液位图像序列的水下三维重建方法，其特征在于，所述步骤S5具体包括以下内容：
S5.1利用图像稠密匹配方法计算图像与图像之间的像素对应关系，记图像中任意像素r＝(xr,yr)T在中的对应像素位置为s＝(xs,ys)T；并利用S3.2中的转换坐标计算方法计算r、s的转换坐标r′＝(x′r,y′r)T、s′＝(x′s,y′s)T；
S5.2对于中每个像素位置r，根据如下公式计算其对应的场景点深度：
<mrow> <mi>d</mi> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mo>(</mo> <mi>w</mi> <mo>(</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> <mo>-</mo> <mi>a</mi> <mo>(</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> <mo>)</mo> <msub> <mover> <mi>D</mi> <mo>^</mo> </mover> <mn>1</mn> </msub> <mo>+</mo> <mo>(</mo> <mi>a</mi> <mo>(</mo> <msup> <mi>s</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> <mo>-</mo> <mi>w</mi> <mo>(</mo> <msup> <mi>s</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> <mo>)</mo> <msub> <mover> <mi>D</mi> <mo>^</mo> </mover> <mi>M</mi> </msub> </mrow> <mrow> <mi>w</mi> <mrow> <mo>(</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> </mrow> <mo>-</mo> <mi>w</mi> <mrow> <mo>(</mo> <msup> <mi>s</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow>
然后，根据d(r)估计像素r对应的场景点的初始三维坐标Xr＝(Xr,Yr,Zr)T，计算公式如下：
<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>X</mi> <mi>r</mi> </msub> <mo>=</mo> <mfrac> <msubsup> <mi>x</mi> <mi>r</mi> <mo>&amp;prime;</mo> </msubsup> <msqrt> <mrow> <msubsup> <mi>x</mi> <mi>r</mi> <mrow> <mo>&amp;prime;</mo> <mn>2</mn> </mrow> </msubsup> <mo>+</mo> <msubsup> <mi>y</mi> <mi>r</mi> <mrow> <mo>&amp;prime;</mo> <mn>2</mn> </mrow> </msubsup> </mrow> </msqrt> </mfrac> <mrow> <mo>(</mo> <msub> <mover> <mi>D</mi> <mo>^</mo> </mover> <mn>1</mn> </msub> <mi>a</mi> <mo>(</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> <mo>+</mo> <mo>(</mo> <mrow> <mi>d</mi> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mover> <mi>D</mi> <mo>^</mo> </mover> <mn>1</mn> </msub> </mrow> <mo>)</mo> <mi>w</mi> <mo>(</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>Y</mi> <mi>r</mi> </msub> <mo>=</mo> <mfrac> <msubsup> <mi>y</mi> <mi>r</mi> <mo>&amp;prime;</mo> </msubsup> <msqrt> <mrow> <msubsup> <mi>x</mi> <mi>r</mi> <mrow> <mo>&amp;prime;</mo> <mn>2</mn> </mrow> </msubsup> <mo>+</mo> <msubsup> <mi>y</mi> <mi>r</mi> <mrow> <mo>&amp;prime;</mo> <mn>2</mn> </mrow> </msubsup> </mrow> </msqrt> </mfrac> <mrow> <mo>(</mo> <msub> <mover> <mi>D</mi> <mo>^</mo> </mover> <mn>1</mn> </msub> <mi>a</mi> <mo>(</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> <mo>+</mo> <mo>(</mo> <mrow> <mi>d</mi> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mover> <mi>D</mi> <mo>^</mo> </mover> <mn>1</mn> </msub> </mrow> <mo>)</mo> <mi>w</mi> <mo>(</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>Z</mi> <mi>r</mi> </msub> <mo>=</mo> <mi>d</mi> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced>
S5.3对于每个场景点，利用局部非线性优化方法对其初始三维坐标求精，优化过程的目标函数为：
<mrow> <msubsup> <mi>X</mi> <mi>r</mi> <mo>*</mo> </msubsup> <mo>=</mo> <munder> <mi>argmin</mi> <msub> <mi>X</mi> <mi>r</mi> </msub> </munder> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>2</mn> </mrow> <mi>M</mi> </munderover> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <mi>h</mi> <mo>(</mo> <mrow> <msub> <mi>I</mi> <mn>1</mn> </msub> <mo>,</mo> <msub> <mi>I</mi> <mi>l</mi> </msub> <mo>,</mo> <msub> <mi>X</mi> <mi>r</mi> </msub> <mo>,</mo> <mi>m</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> </mrow> 4
其中，表示图像中以Xr在该图像上的投影为中心、大小为m×m个像素的图像窗口，与图像中以Xr在该图像上的投影为中心、大小为m×m个像素的图像窗口的归一化互相关系数(normalized cross correlation,NCC)；
其中，m的取值范围是6～10个像素，Xr在图像上的投影点的计算基于直接前向投影方法，相机的旋转矩阵根据θx和θy确定；
S5.4将所有经过优化的三维点融合为三维点云，并利用统计外点移除点云过滤方法剔除孤立点和细小点云片段，获得最终稠密三维重建结果。
CN201710441431.5A 2017-06-13 2017-06-13 基于差异液位图像序列的水下三维重建系统及其方法 Active CN107256563B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710441431.5A CN107256563B (zh) 2017-06-13 2017-06-13 基于差异液位图像序列的水下三维重建系统及其方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710441431.5A CN107256563B (zh) 2017-06-13 2017-06-13 基于差异液位图像序列的水下三维重建系统及其方法

Publications (2)

Publication Number Publication Date
CN107256563A true CN107256563A (zh) 2017-10-17
CN107256563B CN107256563B (zh) 2020-04-07

Family

ID=60023158

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710441431.5A Active CN107256563B (zh) 2017-06-13 2017-06-13 基于差异液位图像序列的水下三维重建系统及其方法

Country Status (1)

Country Link
CN (1) CN107256563B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109754428A (zh) * 2018-11-26 2019-05-14 西北工业大学 一种用于水下双目视觉定位误差的测量方法
CN110111413A (zh) * 2019-04-08 2019-08-09 西安电子科技大学 一种基于水陆共存场景的稀疏点云三维模型建立方法
CN112465950A (zh) * 2020-11-26 2021-03-09 江苏国和智能科技有限公司 深海网箱渔网水下距离测量装置、方法、电子设备及介质

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102903101A (zh) * 2012-09-06 2013-01-30 北京航空航天大学 使用多台相机进行水面数据采集与重建的方法
CN103744085A (zh) * 2014-01-17 2014-04-23 哈尔滨工程大学 水下机器人五分量测距声纳斜井三维成像系统及成像方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102903101A (zh) * 2012-09-06 2013-01-30 北京航空航天大学 使用多台相机进行水面数据采集与重建的方法
CN103744085A (zh) * 2014-01-17 2014-04-23 哈尔滨工程大学 水下机器人五分量测距声纳斜井三维成像系统及成像方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
CHRISTIAN FRIIH ETC: "3D Model Generation for Cities Using Aerial Photographs and Ground Level Laser Scans", 《PROCEEDINGS OF THE 2001 IEEE COMPUTER SOCIETY CONFERENCE ON COMPUTER VISION AND PATTERN RECOGNITION》 *
杨宇: "水下多通道真彩色三维重建与颜色还原方法研究", 《中国博士学位论文全文数据库基础科学辑》 *
谢毓湘: "一种基于局部不变特征的图像特定场景检测方法", 《国防科技大学学报》 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109754428A (zh) * 2018-11-26 2019-05-14 西北工业大学 一种用于水下双目视觉定位误差的测量方法
CN109754428B (zh) * 2018-11-26 2022-04-26 西北工业大学 一种用于水下双目视觉定位误差的测量方法
CN110111413A (zh) * 2019-04-08 2019-08-09 西安电子科技大学 一种基于水陆共存场景的稀疏点云三维模型建立方法
CN112465950A (zh) * 2020-11-26 2021-03-09 江苏国和智能科技有限公司 深海网箱渔网水下距离测量装置、方法、电子设备及介质

Also Published As

Publication number Publication date
CN107256563B (zh) 2020-04-07

Similar Documents

Publication Publication Date Title
CN105678742B (zh) 一种水下相机标定方法
CN104376552B (zh) 一种3d模型与二维图像的虚实配准方法
CN106485690A (zh) 基于点特征的点云数据与光学影像的自动配准融合方法
CN105957007B (zh) 基于特征点平面相似度的图像拼接方法
Zhang et al. A UAV-based panoramic oblique photogrammetry (POP) approach using spherical projection
CN108876749A (zh) 一种鲁棒的镜头畸变校正方法
CN108288292A (zh) 一种三维重建方法、装置及设备
CN106651767A (zh) 一种获取全景图像的方法及装置
CN107204010A (zh) 一种单目图像深度估计方法与系统
CN107358631A (zh) 一种虑及三维畸变的双目视觉重建方法
CN107424181A (zh) 一种改进的图像拼接关键帧快速提取方法
CN104036542B (zh) 一种基于空间光线聚集性的像面特征点匹配方法
CN109615663A (zh) 全景视频校正方法及终端
CN104156957B (zh) 一种稳定高效的高分辨率立体匹配方法
CN104616247B (zh) 一种用于基于超像素sift航拍地图拼接的方法
CN105654547B (zh) 三维重建方法
CN108629829B (zh) 一种球幕相机与深度相机结合的三维建模方法和系统
CN107481279A (zh) 一种单目视频深度图计算方法
CN102903101B (zh) 使用多台相机进行水面数据采集与重建的方法
CN107274483A (zh) 一种物体三维模型构建方法
CN109615661A (zh) 光场相机内参数标定装置及方法
CN105931222A (zh) 用低精度二维平面靶标实现高精度相机标定的方法
CN104155765A (zh) 在拼接式集成成像显示器中校正三维图像的方法和设备
CN106910208A (zh) 一种存在运动目标的场景图像拼接方法
CN107103056A (zh) 一种基于局部标识的双目视觉室内定位数据库建立方法及定位方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
CB03 Change of inventor or designer information

Inventor after: Kang Lai

Inventor after: Wu Lingda

Inventor after: Bai Liang

Inventor after: Wei Yingmei

Inventor after: Lao Songyang

Inventor after: Jiang Jie

Inventor after: Xie Yuxiang

Inventor before: Wu Lingda

Inventor before: Bai Liang

Inventor before: Wei Yingmei

Inventor before: Lao Songyang

Inventor before: Jiang Jie

Inventor before: Liu Yuxiang

CB03 Change of inventor or designer information
GR01 Patent grant
GR01 Patent grant