CN113066016B - 一种基于图信号处理的三维动态点云修复方法 - Google Patents

一种基于图信号处理的三维动态点云修复方法 Download PDF

Info

Publication number
CN113066016B
CN113066016B CN202110190191.2A CN202110190191A CN113066016B CN 113066016 B CN113066016 B CN 113066016B CN 202110190191 A CN202110190191 A CN 202110190191A CN 113066016 B CN113066016 B CN 113066016B
Authority
CN
China
Prior art keywords
block
point
point cloud
frame
points
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
CN202110190191.2A
Other languages
English (en)
Other versions
CN113066016A (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.)
Peking University
Original Assignee
Peking University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Peking University filed Critical Peking University
Priority to CN202110190191.2A priority Critical patent/CN113066016B/zh
Publication of CN113066016A publication Critical patent/CN113066016A/zh
Application granted granted Critical
Publication of CN113066016B publication Critical patent/CN113066016B/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
    • G06T5/00—Image enhancement or restoration
    • G06T5/77—Retouching; Inpainting; Scratch removal
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/90—Details of database functions independent of the retrieved data types
    • G06F16/901—Indexing; Data structures therefor; Storage structures
    • G06F16/9024—Graphs; Linked lists
    • 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/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/22—Matching criteria, e.g. proximity measures
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T15/00—3D [Three Dimensional] image rendering
    • G06T15/08—Volume rendering
    • 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)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Databases & Information Systems (AREA)
  • Mathematical Physics (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Software Systems (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Computation (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Algebra (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Evolutionary Biology (AREA)
  • Computer Graphics (AREA)
  • Image Analysis (AREA)

Abstract

本发明公开了一种基于图信号处理的三维动态点云修复方法，其步骤包括：1)将点云序列S有缺失区域Pf分割成n个立方块并从中选择包含缺失数据的立方块作为目标块ct；2)选取与ct相似度最高的候选块作为帧内源块cs，并对几何结构进行位置匹配，得到最终的帧内源块
Figure DDA0002943758280000015
3)在Pf‑1、Pf+1中找到与ct相对位置相同的块c′t、c″t；4)在c′t、c″t的包围盒中搜索ct中每个点的最近邻点；5)寻找Pf‑1、Pf+1中包含ct的最近邻点最多的立方块，作为帧间源块
Figure DDA0002943758280000011
并对几何结构进行位置匹配，得到帧间源块
Figure DDA0002943758280000012
6)利用
Figure DDA0002943758280000013
和
Figure DDA0002943758280000014
修复ct。

Description

一种基于图信号处理的三维动态点云修复方法
技术领域
本发明属于图信号处理领域，涉及一种点云修复方法，尤其涉及基于图信号处理的三维动态点云修复方法。
背景技术
使用三维扫描仪和深度感测设备采集到的动态点云数据会不可避免地产生数据误差：对于一些扫描不到的区域，在扫描结果中会以缺失数据的孔洞的形式呈现。缺失数据的原因包括由于扫描角度的不完全，激光扫描仪本身的局限性，物体几何结构过于复杂繁复，以及物体运动变化太快等。
目前还没有关于动态点云修复的算法，但有一些静态点云修复的算法，可以通过对动态点云中的每一帧含有缺失区域的静态点云应用这些算法，实现动态点云修复。这些静态点云修复算法主要分为两类：修复物体本身的缺失区域，修复在采集过程中产生的缺失区域，和修复特殊数据。
修复物体本身的缺失区域的方法的缺失区域一般比较大，主要依赖线上数据库作为填充内容的参考。Sahay等人先尝试通过参考缺失区域周围的几何信息，从已有的点云库中寻找相似数据进行填充；后又提出另一种将点云投影到深度图上的方法，使用字典学习从在线的深度图数据库中搜索相似的深度图，将未知部分对应的信息粘贴到空缺区域，最后将补好的深度图投影回点云。此类方法对被搜索的数据库要求较高，需要足够的与待修补点云相似的序列，且直接填充难免在衔接处产生几何结构的不对齐，这是在最后平滑处理时难以解决的问题。此外，在深度图方法中的两次投影过程将不可避免地产生几何结构的损失。另外，Shankar等人通过提取孔洞周围的边界点来提取其几何特征，拟合缺失部分的曲面后采样插值得到填充的点。Dinesh等人先确定缺失区域的修复顺序，然后基于最小的旋转差异搜索最佳匹配区域来填充缺失区域。这类方法由于参考的数据来源单一，结果不可避免地有几何失真。
修复在采集过程中产生的缺失区域方法关注的是由于扫描设备的限制而产生的孔洞。这种孔洞通常要比上一类的孔洞小，因此数据本身的信息通常足以用于修复，一般通过对缺失区域周围的点的分析来预测缺失区域内部的几何结构。例如Quinsat等人和Wang等人基于输入点云创建三角形网格，然后基于孔附近的网格结构在孔洞上构建网格，最后对缺失区域进行插值得到点。但这类方法的结果依赖于网格构造的质量，并且构建网格过程的计算复杂度很高。另外，Lozes等人提出了通过偏微分方法求解孔洞周围点云的优化方程来进行修补的方法；Lin等人则将孔洞分解为几个子洞后，分别利用张量投票方法提取几何特征单独填充；Muraki等人在孔洞附近拟合了一个曲面并进行插值来进行修复。以上的方法由于仅参考了局部的邻域信息，所以其修补结果比实际内容要更趋于平面。而且当几何结构复杂时，在边界附近可能会出现伪影。
此外，还有一些针对静态点云的工作处理具有特殊数据形式的点云数据。例如，Friedman的方法中的建筑物的几何规则点云，Li的方法中的扫描线数据，Wu的方法中的人体数据中扁平的条状孔洞。这类方法具有较强的数据针对性，不适用于其他形式的点云。
发明内容
针对上述问题和相关方法的缺陷，本发明提出了一种基于图信号处理的三维动态点云修复方法。
在介绍本发明方法的主要步骤之前，首先介绍一下谱图理论中的基本概念，谱图理论将是本发明中的方法的核心。谱图理论包括以下几个内容：
(1)图和图拉普拉斯算子：
我们定义一个无向图
Figure BDA0002943758260000021
其中
Figure BDA0002943758260000022
为图
Figure BDA0002943758260000023
上顶点的集合，顶点的数量
Figure BDA0002943758260000024
ε为边的集合；W为带权邻接矩阵。其中W是一个N×N的实对称矩阵，wi，j是连接顶点i与顶点j的边的权重，通常使用非负值作为权重。
图拉普拉斯算子通常被定义为
Figure BDA0002943758260000025
其中D为对角矩阵，对角矩阵D中对角线上的第i个元素
Figure BDA0002943758260000026
根据带权邻接矩阵W得到对角矩阵D，进而得到
Figure BDA0002943758260000027
(2)图信号平滑算子：
图信号是指驻留在图的顶点上的数据，例如社交网络，交通网络，传感器网络和神经元网络。例如，如果我们在点云上构造一个K近邻图(K-NN图)，那么每个点的坐标可以被看作是在K-NN图上定义的图信号。
对于在图G上定义的信号z，如果满足下述条件，那么称信号z相对于G的拓扑结构是平滑的：
Figure BDA0002943758260000028
其中，∈是一个参数，i～j表示图中相连的两个顶点。为了满足上述公式，zi和zj相对于权重较大的边时必须相似，对于权重较小的边时zi和zj可能完全不同。因此，上述公式强制z适应G的拓扑结构，从而创造了图信号的平滑性。由于zTLz＝∑i～jwi，j(zi-zj)2，上述公式还可以写为zTLz＜∈。
这个算子将在后续的点云修复步骤中作为一个正则化项使用。
在介绍完基本的谱图理论之后，开始讨论本发明中的方法。对输入的点云序列S＝{P1，P2，…，Pq}中每一帧有缺失区域的点云Pf，本发明将实施五个主要步骤进行修复：
1)点云的预处理：
(1)将当前帧点云Pf分割成互相有重叠部分的立方块{c1，c2，…，cn}作为衡量点云的几何信息的单位。
(2)选择包含缺失数据的立方块作为目标块ct。
(3)其他具有足够点数的块及其镜面对称的块作为候选块cc，我们将从中寻找与目标块最相似的块。
2)搜索帧内相似区域：
(1)计算目标块ct和所有候选块的法向量的直流分量(Direct Component，DC)。
(2)在目标块ct和所有候选块上应用K-NN算法对块中的点建立K-NN图。
(3)计算目标块ct和每一候选块的法向量的各向异性图全变分(AnisotropicGraph Total Variation，AGTV)。
(4)计算每一候选块与目标块ct的DC差距和AGTV差距，得到所有候选块与目标块间的相似度。
(5)选取相似度最高的候选块作为帧内源块cs。利用专利ZL201811610195.6中的方法进行cs中的几何结构与ct中几何结构的相对位置的匹配，得到最终的帧内源块
Figure BDA0002943758260000031
用于最后的修复步骤。
3)搜索帧间相关区域：
(1)在与点云Pf相邻的前一帧Pf-1中找到与ct相对位置相同的块c′t。
(2)围绕c′t在Pf-1中创建一个包围盒
Figure BDA0002943758260000032
(3)在
Figure BDA0002943758260000033
中搜索目标块ct中每个点的最近邻点。
(4)在包围盒
Figure BDA0002943758260000034
中创建一个滑动立方窗口
Figure BDA0002943758260000035
寻找帧Pf-1中包含ct的最近邻点最多的立方块，作为帧间源块
Figure BDA0002943758260000036
(5)利用专利ZL201811610195.6中的方法进行帧间源块
Figure BDA0002943758260000041
中的几何结构与ct中几何结构的相对位置的匹配，得到最终的前一帧的帧间源块
Figure BDA0002943758260000042
同理，通过后续步骤(6)～(10)可以在后一帧中得到帧间源块
Figure BDA0002943758260000043
(6)在与点云Pf相邻的后一帧Pf+1中找到与ct相对位置相同的块c″t。
(7)围绕c″t在Pf+1中创建一个包围盒
Figure BDA0002943758260000044
(8)在
Figure BDA0002943758260000045
中搜索目标块ct中每个点的最近邻点。
(9)在包围盒
Figure BDA0002943758260000046
中创建一个滑动立方窗口
Figure BDA0002943758260000047
寻找帧Pf+1中包含ct的最近邻点最多的立方块，作为帧间源块
Figure BDA0002943758260000048
(10)在帧间源块
Figure BDA0002943758260000049
上进行与对
Figure BDA00029437582600000410
相同的块内几何结构匹配，得到最终的前一帧的帧间源块
Figure BDA00029437582600000411
4)空间-时间图初始化：
(1)建立三重立方块g。
(2)在g上为修复结果块cr中的点建图
Figure BDA00029437582600000412
对于cr内的点之间，所建的子图是K-NN图。
(3)对于cr和
Figure BDA00029437582600000413
的点之间，在cr中的每个点qr，k和其在
Figure BDA00029437582600000414
中的最近邻点
Figure BDA00029437582600000415
之间连接一条无权边；qr，k为cr中的第k个点，
Figure BDA00029437582600000416
为
Figure BDA00029437582600000417
中的第1个点，是qr，k在
Figure BDA00029437582600000418
中对应的最邻近点。在cr和
Figure BDA00029437582600000419
之间建的图同理。
(4)对于cr和
Figure BDA00029437582600000420
的点之间，在cr中的每个点qr，k和其在
Figure BDA00029437582600000421
中的最近邻点
Figure BDA00029437582600000422
之间连接一条无权边。
5)修复缺失区域：
(1)利用最终的帧内源块
Figure BDA00029437582600000423
和
Figure BDA00029437582600000424
中的信息修复目标块ct中的缺失区域，将这个修复工作转变为一个正则项为图平滑算子和二阶时间连贯项的优化问题，列出优化方程。
(2)将优化方程进行简化和重整。
(3)通过对重整后的优化方程中的每个变量进行循环迭代求解得到最终的修复结果块cr，用结果块在点云中替换目标块，得到最终的修复后的点云。
与现有技术相比，本发明的积极效果为：
本发明是第一个解决三维动态点云修复问题的技术，提出了搜索帧间相关区域的方法，以及将三维动态点云修复操作公式化并求解的技术。与现有技术相比，本发明的修复结果具有较高水平。
使用公知的针对点云的几何失真性的度量GPSNR和NSHD作为评估指标。其中，GPSNR是基于传统的PSNR设计出的衡量点云失真或其噪声水平的客观标准，两个点云之间的GPSNR值越大，则越相似。NSHD是基于单边Hausdorff距离设计的归一化的Hausdorff距离，衡量两个点云之间的距离差，两个点云之间的NSHD值越小，则越相似。
表1和表2分别展示了本发明方法针对人为孔洞的修复结果的GPSNR和NSHD客观评测结果。可以看到本发明的修复结果从GPSNR和NSHD两个指标来看，都明显优于其他所有竞争方法。具体而言，本发明方法相对于基于网格的方法的GPSNR结果平均增益28.07dB，相对于基于微分方程的方法平均增益17.18dB，相对于专利ZL201811610195.6平均增益7.76dB。在表2中，本发明方法结果的NSHD比其他方法低得多，比竞争方法中的最佳方法——专利ZL201811610195.6低至少三倍。
表1：与ground truth之间的GPSNR评估对比
基于网格 基于微分方程 专利ZL201811610195.6 本发明
Longdress 11.7926 30.3883 41.5686 44.5519
Loot 16.4451 27.3715 40.1546 48.7913
Maria 12.2943 28.4910 34.9747 41.1820
Phili 9.9593 19.4861 35.4972 41.6200
Queen 14.0248 22.4650 29.6086 39.1679
Redandblack 13.1772 24.4810 33.9921 41.8005
Sarah 14.5581 26.4108 30.2719 39.4186
Skiing 20.1247 29.9233 36.0379 45.9341
Soldier 17.4697 23.1571 34.5062 44.9148
UlliWegner 24.9424 31.5455 41.3037 48.1035
表2：与ground truth之间的NSHD评估对比
Figure BDA0002943758260000051
Figure BDA0002943758260000061
附图说明
图1为本发明提出的三维动态点云修复算法框架。
具体实施方式
下面将针对搜索帧内相似区域步骤、搜索帧间相关区域步骤、空间-时间图初始化步骤和修复缺失区域步骤，对本发明的详细方法流程作进一步地描述：
搜索帧内相似区域：
(1)将块中点数少于目标块ct中点数80％的候选块去除，并将剩余的块关于x，y平面进行镜面对称得到新的一部分候选块。
(2)计算目标块的法向量的DC值d(ct)和所有候选块的法向量的DC值d(cc)。一个立方块ci包含一组点{ci，1，ci，2，…，ci，m}(m是块内点的数量)，每一个点都对应一个法向量表示为{ni，1，ni，2，…，ni，m}。计算一个立方块ci的DC值d(ci)如下：
Figure BDA0002943758260000062
(3)在目标块和所有候选块上应用K-NN算法对块中的点建立图。以点的坐标为图信号，基于两点间的距离远近建立谱图理论中提到的K-NN图，其权重为：
Figure BDA0002943758260000063
其中k～l表示点pk和点pl之间有一条边。
(4)计算目标块的法向量的AGTV值v(ct)和所有候选块的法向量的AGTV值v(cc)。计算一个立方块ci的AGTV值v(ci)如下：
Figure BDA0002943758260000064
(5)计算所有候选块与目标块的DC差距δD(ct，cc)和AGTV差距δV(ct，cc)如下：
δD(ct，cc)＝|<d(ct)，d(cc)>|，
δV(ct，cc)＝|v(ct)-v(cc)|，
由此计算出所有候选块与目标块间的相似度如下：
δ(ct，cc)＝exp{-[δD(ct，cc)+δV(ct，cc)]}，
(6)选取相似度最高的候选块作为帧内源块cs。利用我们之前申请的专利(ZL201811610195.6)中的方法进行cs中的几何结构与ct中几何结构的相对位置的匹配，得到最终的帧内源块
Figure BDA0002943758260000071
用于最后的修复步骤。
搜索帧间相关区域：
(1)在Pf-1中找到与ct相对位置相同的块c′t：
s(c′t)＝s(ct)，
其中s(c′t)和s(ct)分别表示块c′t和ct的中心点的坐标，也就是说c′t和ct在各自帧内的位置是相同的。
(2)围绕c′t在Pf-1中创建一个包围盒
Figure BDA0002943758260000072
(即比切割出的立方块ci大一圈的立方块)：
Figure BDA0002943758260000073
其中
Figure BDA0002943758260000074
和s(c′t)分别表示块
Figure BDA0002943758260000075
和c′t的中心点的坐标，也就是说
Figure BDA0002943758260000076
和c′t在各自帧内的位置是相同的。
(3)在
Figure BDA0002943758260000077
中搜索目标块ct中每个点的最近邻点。
(4)在包围盒
Figure BDA0002943758260000078
中创建一个滑动立方窗口
Figure BDA0002943758260000079
则帧Pf-1中的帧间源立方块
Figure BDA00029437582600000710
为：
Figure BDA00029437582600000711
其中
Figure BDA00029437582600000712
是
Figure BDA00029437582600000713
中包含的ct的最近邻点数量，即
Figure BDA00029437582600000714
包含最多的ct在Pf-1中的最近邻点。
(5)对帧内源块
Figure BDA00029437582600000715
中的几何结构与c′t中的几何结构进行位置匹配，得到帧间源块
Figure BDA00029437582600000716
(6)在Pf+1中找到与ct相对位置相同的块c″t：
s(c″t)＝s(ct)，
其中s(c″t)和s(ct)分别表示块c′t和ct的中心点的坐标，也就是说c′t和ct在各自帧内的位置是相同的。
(7)围绕c″t在Pf+1中创建一个包围盒
Figure BDA00029437582600000717
(即比切割出的立方块ci大一圈的立方块)：
Figure BDA0002943758260000081
其中
Figure BDA0002943758260000082
和s(c″t)分别表示块
Figure BDA0002943758260000083
和c″t的中心点的坐标，也就是说
Figure BDA0002943758260000084
和c″t在各自帧内的位置是相同的。
(8)在
Figure BDA0002943758260000085
中搜索目标块ct中每个点的最近邻点。
(9)在包围盒
Figure BDA0002943758260000086
中创建一个滑动立方窗口
Figure BDA0002943758260000087
则帧Pf-1中的帧间源立方块
Figure BDA0002943758260000088
为：
Figure BDA0002943758260000089
其中
Figure BDA00029437582600000810
是
Figure BDA00029437582600000811
中包含的ct的最近邻点数量，即
Figure BDA00029437582600000812
包含最多的ct在Pf-1中的最近邻点。
(10)对帧内源块
Figure BDA00029437582600000813
中的几何结构与c″t中的几何结构进行位置匹配，得到帧间源块
Figure BDA00029437582600000814
空间-时间图初始化：
(1)建立三重立方块
Figure BDA00029437582600000815
cr为目标块ct对应的修复结果块。
(2)在g上为cr中的点建图
Figure BDA00029437582600000816
对于cr内的点之间建立空间图，这个子图是K-NN图：
Figure BDA00029437582600000817
其中k～l表示点pk和点pl之间有一条边。根据K-NN图，当点pl属于点pk的最近的前K个邻居时，两者之间连一条边。
(3)对于cr和
Figure BDA00029437582600000818
的点之间建立时间图，在cr中的每个点qr，k和其在
Figure BDA00029437582600000819
中的最近邻点
Figure BDA00029437582600000820
之间连接一条无权边(当点
Figure BDA00029437582600000821
是点qr，k在源块
Figure BDA00029437582600000822
中的最近邻点时，两点之间连接一条无权边)：
Figure BDA00029437582600000823
(4)在cr和
Figure BDA0002943758260000091
之间建的时间图同理。对于cr和
Figure BDA0002943758260000092
的点之间建立时间图，在cr中的每个点qr，k和其在
Figure BDA0002943758260000093
中的最近邻点
Figure BDA0002943758260000094
之间连接一条无权边：
Figure BDA0002943758260000095
由于cr是未知的，所以在建立第一种边的时候用
Figure BDA0002943758260000096
替代cr；在建立第二种边时在未知区域用
Figure BDA0002943758260000097
替代cr，而在已知区域用ct替代cr。
修复缺失区域：
(1)利用最终的源块
Figure BDA0002943758260000098
和
Figure BDA0002943758260000099
中的信息修复目标块中的缺失区域，将这个修复工作转变为一个正则项为图平滑算子和二阶时间连贯项的优化问题，列出优化方程如下：
Figure BDA00029437582600000910
其中：
·
Figure BDA00029437582600000911
是想要得到的结果块；
·Ω是一个M3×M3的对角矩阵，其对角元素Ωi.i∈{0，1}，其中0表示已知点，1表示未知点，所以Ωcr和
Figure BDA00029437582600000912
分别提取出了cr和
Figure BDA00029437582600000913
内的缺失区域，而
Figure BDA00029437582600000914
与Ω互补，是表示块内的已知区域。
·Wf-1，f是在
Figure BDA00029437582600000915
与cr之间建立的时间图连接的权重矩阵，Wf+1，f则是在
Figure BDA00029437582600000916
与cr之间建立的时间图连接的权重矩阵。
·
Figure BDA00029437582600000917
是在g上建立的空间-时间图的拉普拉斯矩阵。
·α，β和γ是三个权重参数。其中α会影响修复区域和帧内源块的对应区域的参考度；β会影响修复结果块中几何结构在空间上的平滑度和时间上的连贯性；而γ则会影响修复后的动态点云动作上的连贯性，即二阶时间连贯性。α越大就会使结果块与帧内源块间的差距越小，即会使结果块在修复区域与帧内源块越相似；β越大就会使结果块与帧间源块间的差距越小，即会使结果块在修复区域与帧间源块越相似；γ越大就会使相邻两帧的图间的差距越小，即会使相邻帧间更连贯。
·tr()表示矩阵的迹，tr(A)表示方阵A的迹。
(2)将优化方程进行简化和重整。
g是由cr，
Figure BDA00029437582600000918
和
Figure BDA00029437582600000919
组成的，所以图
Figure BDA00029437582600000920
的权重矩阵Wg可以写作如下形式：
Figure BDA00029437582600000921
其中Wf-1，f-1表示
Figure BDA0002943758260000101
中的点之间的权重，Wf-1，f表示
Figure BDA0002943758260000102
中的点和cr中的点之间的权重，以此类推。由于Wg是对称的，所以有
Figure BDA0002943758260000103
和
Figure BDA0002943758260000104
在谱图理论中提到，图
Figure BDA0002943758260000105
的度矩阵Dg是对角矩阵，所以可以写作：
Figure BDA0002943758260000106
度矩阵Dg通过Wg计算得到的，计算结果也是一个对角矩阵，这里与Wg对应，分为九块表示；所以图
Figure BDA0002943758260000107
的拉普拉斯矩阵
Figure BDA0002943758260000108
可以计算得到：
Figure BDA0002943758260000109
所以帧间连贯项
Figure BDA00029437582600001010
可以展开得到：
Figure BDA00029437582600001011
在展开式中可以观察到展开后的一些项与修复任务无关，如一、三、七、九项都仅与
Figure BDA00029437582600001012
和
Figure BDA00029437582600001013
相关，而与要求解的cr无关，因此可以将它们视为优化目标中的常量。此外，一些成对的项是相等的。例如，第二项和第四项互为转置，所以它们的值相等。同理，第六项和第八项也相等。因此，它们可以合并为一个项。所以可将上式简化为：
Figure BDA00029437582600001014
其中C是一个常数，
Figure BDA00029437582600001015
Figure BDA00029437582600001016
是展开式中的第一、三、七、九项。
Figure BDA00029437582600001017
是一个广义拉普拉斯矩阵，即一个在非对角项上具有非正值的对称矩阵。
则优化方程可以简化为：
Figure BDA0002943758260000111
Figure BDA0002943758260000112
Figure BDA0002943758260000113
Figure BDA0002943758260000114
Figure BDA0002943758260000115
Figure BDA0002943758260000116
其中
Figure BDA0002943758260000117
是全0向量，
Figure BDA0002943758260000118
是全1向量。
(3)通过对重整后的优化方程中的每个变量进行循环迭代求解：
cr的优化求解：按照空间-时间图初始化Wf-1，f，Wf+1，f，
Figure BDA0002943758260000119
后得到关于cr的优化方程为
Figure BDA00029437582600001110
这个问题方程实际上是一个二次规划问题，对上式用cr进行求导，并将导数置为0后，可得到cr的最优的闭式解为
Figure BDA00029437582600001111
Figure BDA00029437582600001112
的优化求解：按照空间-时间图初始化Wf-1，f，Wf+1，f，并通过cr的优化求解更新cr后，得到关于
Figure BDA00029437582600001113
的优化方程为
Figure BDA00029437582600001114
Figure BDA00029437582600001115
Figure BDA00029437582600001116
Figure BDA00029437582600001117
这是一个凸优化问题，可以通过现有的凸优化工具包，例如MATLAB中的CVX工具解决。
Wf-1，f和Wf+1，f的优化求解：按照上述优化求解更新cr和
Figure BDA00029437582600001221
后，固定Wf+1，f得到关于Wf-1，f的优化方程为
Figure BDA0002943758260000121
s.t.0＜Wf-1，f·1≤1
将矫正有限内存-拟牛顿法(Orthant-Wise Limited-memory Quasi-Newton，OWL-QN)拓展到具有矩阵变量的l1优化问题，通过限制目标函数的象限来保证l1范数的连续性和可微性。具体来说，上式的目标函数
Figure BDA0002943758260000122
的牛顿迭代公式为：
Figure BDA0002943758260000123
其中
Figure BDA0002943758260000124
和
Figure BDA0002943758260000125
分别是第i和i+1次迭代后的Wf-1，f。
Figure BDA0002943758260000126
是目标函数第i次迭代时关于
Figure BDA0002943758260000127
的一阶导数。
F(Wf-1，f)不连续可微，所以在迭代中通过限制Wf-1，f的象限来计算伪梯度作为F’(Wf-1，f)。首先，对于Wf-1，f中的元素wk，l，限制其迭代的象限方向如下：
Figure BDA0002943758260000128
其中符号函数σ根据实变量分别为{负数，零，正数}，对应函数值分别为{-1，0，1}。然后如OWL-QN的方法，函数
Figure BDA0002943758260000129
在
Figure BDA00029437582600001210
的伪梯度F’(Wf-1，f)定义为：
Figure BDA00029437582600001211
其中
Figure BDA00029437582600001212
和
Figure BDA00029437582600001213
分别为
Figure BDA00029437582600001214
的左偏导和右偏导，定义为：
Figure BDA00029437582600001215
其中
Figure BDA00029437582600001216
是目标函数
Figure BDA00029437582600001217
中的线性项，
Figure BDA00029437582600001218
是目标函数中l1范数的左偏导和右偏导，其中的元素
Figure BDA00029437582600001219
是关于
Figure BDA00029437582600001220
的函数：
Figure BDA0002943758260000131
这样就通过计算伪梯度的改进牛顿迭代法OWL-QN实现了Wf-1，f的求解，同理可以固定Wf-1，f来求解Wf+1，f。
通过以上交替迭代求解各个变量cr，Wf-1，f，Wf+1，f和
Figure BDA0002943758260000132
可以得到优化方程的最优解。最后，在目标帧Pf中用结果块cr替换目标块，得到修复后的点云。
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (4)

1.一种基于图信号处理的三维动态点云修复方法，其步骤包括：
1)对于输入的点云序列S＝{P1，P2，…，Pq}，Pq为第q帧对应的点云，如果点云序列S中第f帧对应的点云Pf中有缺失区域，则将点云Pf分割成n个立方块{c1，c2，…，cn}并从中选择包含缺失数据的立方块作为目标块ct；以及从{c1，c2，…，cn}中选择点数超过设定阈值的立方块及其镜面对称块作为候选块；
2)从所述候选块中选取与该目标块ct相似度最高的候选块作为帧内源块cs；对cs中的几何结构与ct中的几何结构进行位置匹配，得到最终的帧内源块
Figure FDA0003608426980000011
3)在点云Pf-1中找到与ct相对位置相同的块c′t、在点云Pf+1中找到与ct相对位置相同的块c″t；其中第f帧的前一帧第f-1帧对应的点云为Pf-1，第f帧的后一帧第f+1帧对应的点云为Pf+1；
4)围绕c′t在Pf-1中创建一个包围盒
Figure FDA0003608426980000012
并在
Figure FDA0003608426980000013
中搜索目标块ct中每个点的最近邻点；围绕c″t在Pf+1中创建一个包围盒
Figure FDA0003608426980000014
并在
Figure FDA0003608426980000015
中搜索目标块ct中每个点的最近邻点；
5)寻找点云Pf-1中包含ct的最近邻点最多的立方块，作为帧间源块
Figure FDA0003608426980000016
寻找点云Pf+1中包含ct的最近邻点最多的立方块，作为帧间源块
Figure FDA0003608426980000017
6)对帧内源块
Figure FDA0003608426980000018
中的几何结构与ct中的几何结构进行位置匹配，得到帧间源块
Figure FDA0003608426980000019
对帧内源块
Figure FDA00036084269800000110
中的几何结构与ct中的几何结构进行位置匹配，得到帧间源块
Figure FDA00036084269800000111
7)利用
Figure FDA00036084269800000112
和
Figure FDA00036084269800000113
中的信息修复目标块ct中的缺失区域，其方法为：71)建立三重立方块
Figure FDA00036084269800000114
在该三重立方块g上为修复结果块cr中的点建立一无向图
Figure FDA00036084269800000115
其中，cr为目标块ct对应的修复结果块，
Figure FDA00036084269800000116
为无向图
Figure FDA00036084269800000117
上顶点的集合，顶点的数量
Figure FDA00036084269800000118
εg为边的集合；Wg为带权邻接矩阵；基于无向图
Figure FDA00036084269800000119
对cr内的点建立若干空间图，在cr中的每个点qr，k和其在
Figure FDA00036084269800000120
中的最近邻点
Figure FDA00036084269800000121
之间连接一条无权边，建立第一时间图；在cr中的每个点qr，k和其在
Figure FDA00036084269800000122
中的最近邻点
Figure FDA00036084269800000123
之间连接一条无权边，建立第二时间图；72)将修复目标块ct中的缺失区域的处理转变为一个正则项为图平滑算子和二阶时间连贯项的优化问题；其中优化方程为
Figure FDA0003608426980000021
Figure FDA0003608426980000022
其中，α、β和γ是三个权重参数，Ωcr为cr内的缺失区域、
Figure FDA0003608426980000023
为cr内的已知区域、
Figure FDA0003608426980000024
为
Figure FDA0003608426980000025
内的缺失区域、
Figure FDA0003608426980000026
为ct内的已知区域；Wf-1，f是在
Figure FDA0003608426980000027
与cr之间建立的第一时间图连接的权重矩阵，Wf+1，f是在
Figure FDA0003608426980000028
与cr之间建立的第二时间图连接的权重矩阵；
Figure FDA0003608426980000029
是在g上建立的空间-时间图的拉普拉斯矩阵；73)通过对该优化方程中的每个变量进行循环迭代求解得到最终的修复结果块cr，用最终得到的修复结果块cr替换目标块ct。
2.如权利要求1所述的方法，其特征在于，获取帧内源块cs的方法为：
21)计算目标块ct和所有候选块的法向量的直流分量DC；
22)在目标块ct和所有候选块上应用K-NN算法对块中的点建立K-NN图；
23)计算目标块ct和每一候选块的法向量的各向异性图全变分AGTV；
24)根据每一候选块与目标块ct的DC差距、AGTV差距，确定各候选块与目标块间的相似度；
25)选取相似度最高的候选块作为帧内源块cs。
3.如权利要求1所述的方法，其特征在于，对所述优化方程进行简化和重整，然后进行步骤73)；其中对所述优化方程进行简化和重整的方法为：
a)获取图
Figure FDA00036084269800000210
的权重矩阵
Figure FDA00036084269800000211
其中Wf-1，f-1表示
Figure FDA00036084269800000212
中的点之间的权重，Wf-1，f表示
Figure FDA00036084269800000213
中的点和cr中的点之间的权重，Wf-1，f+1表示
Figure FDA00036084269800000214
中的点和
Figure FDA00036084269800000215
中的点之间的权重，Wf，f-1表示cr中的点和
Figure FDA00036084269800000216
中的点之间的权重，Wf，f表示cr中的点之间的权重，Wf，f+1表示cr中的点和
Figure FDA00036084269800000217
中的点之间的权重，Wf+1，f-1表示
Figure FDA00036084269800000218
中的点和
Figure FDA00036084269800000219
中的点之间的权重，Wf+1，f表示
Figure FDA00036084269800000220
中的点和cr中的点之间的权重，Wf+1，f+1表示
Figure FDA00036084269800000221
中的点之间的权重；
b)计算图
Figure FDA00036084269800000222
的度矩阵
Figure FDA00036084269800000223
c)计算图
Figure FDA00036084269800000224
的拉普拉斯矩阵
Figure FDA0003608426980000031
d)根据a)～c)，得到
Figure FDA0003608426980000032
其中C是一个常数，
Figure FDA0003608426980000033
是一个在非对角项上具有非正值的对称矩阵；
e)所述优化方程进行简化为：
Figure FDA0003608426980000034
其中，
Figure FDA0003608426980000035
Figure FDA0003608426980000036
Figure FDA0003608426980000037
是全0向量，
Figure FDA0003608426980000038
是全1向量。
4.如权利要求1所述的方法，其特征在于，n个立方块{c1，c2，…，cn}中相邻两立方块之间互相有重叠部分。
CN202110190191.2A 2021-02-18 2021-02-18 一种基于图信号处理的三维动态点云修复方法 Active CN113066016B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110190191.2A CN113066016B (zh) 2021-02-18 2021-02-18 一种基于图信号处理的三维动态点云修复方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110190191.2A CN113066016B (zh) 2021-02-18 2021-02-18 一种基于图信号处理的三维动态点云修复方法

Publications (2)

Publication Number Publication Date
CN113066016A CN113066016A (zh) 2021-07-02
CN113066016B true CN113066016B (zh) 2022-08-05

Family

ID=76558912

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110190191.2A Active CN113066016B (zh) 2021-02-18 2021-02-18 一种基于图信号处理的三维动态点云修复方法

Country Status (1)

Country Link
CN (1) CN113066016B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114310872B (zh) * 2021-11-29 2023-08-22 杭州电子科技大学 一种基于dgg点云分割网络的机械臂自动打菜方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106204473A (zh) * 2016-06-30 2016-12-07 扬州大学 基于Laplacian优化的非完备模型局部变形的恢复方法
CN109859114A (zh) * 2018-12-27 2019-06-07 北京大学 基于局域平滑性和非局域相似性的三维点云修复方法
CN110349091A (zh) * 2018-04-08 2019-10-18 北京大学 基于图信号处理的点云修复方法、装置及终端
CN111145094A (zh) * 2019-12-26 2020-05-12 北京工业大学 一种基于表面法向引导与图拉普拉斯先验约束的深度图增强方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10474161B2 (en) * 2017-07-03 2019-11-12 Baidu Usa Llc High resolution 3D point clouds generation from upsampled low resolution lidar 3D point clouds and camera images
TWI676153B (zh) * 2018-07-25 2019-11-01 國立中央大學 利用2d影像資訊修補不完整3d深度影像之方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106204473A (zh) * 2016-06-30 2016-12-07 扬州大学 基于Laplacian优化的非完备模型局部变形的恢复方法
CN110349091A (zh) * 2018-04-08 2019-10-18 北京大学 基于图信号处理的点云修复方法、装置及终端
CN109859114A (zh) * 2018-12-27 2019-06-07 北京大学 基于局域平滑性和非局域相似性的三维点云修复方法
CN111145094A (zh) * 2019-12-26 2020-05-12 北京工业大学 一种基于表面法向引导与图拉普拉斯先验约束的深度图增强方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
【点云系列】基于图结构的点云快速重采样;^_^ 晅菲;《CSDN》;20200617;第1-16页 *
3D DYNAMIC POINT CLOUD INPAINTING VIA TEMPORAL CONSISTENCY ON GRAPHS;Zeqing Fu等;《2020 IEEE 多媒体与博览会国际会议（ICME）》;20200609;第1-6页 *
基于类别-实例分割的室内点云场景修复补全;缪永伟等;《计算机学报》;20201231;第1-15页 *

Also Published As

Publication number Publication date
CN113066016A (zh) 2021-07-02

Similar Documents

Publication Publication Date Title
CN108492281B (zh) 一种基于生成式对抗网络的桥梁裂缝图像障碍物检测与去除的方法
CN111063021B (zh) 一种空间运动目标的三维重建模型建立方法及装置
CN108038906B (zh) 一种基于图像的三维四边形网格模型重建方法
Lozes et al. Partial difference operators on weighted graphs for image processing on surfaces and point clouds
CN108460760B (zh) 一种基于生成式对抗网络的桥梁裂缝图像判别修复方法
CN109859114B (zh) 基于局域平滑性和非局域相似性的三维点云修复方法
WO2017201751A1 (zh) 虚拟视点视频、图像的空洞填充方法、装置和终端
Zhang et al. Critical regularizations for neural surface reconstruction in the wild
Lee et al. Perceptual organization of 3D surface points
CN113077553B (zh) 一种基于表面属性的三维模型分割方法
CN111242864B (zh) 一种基于Gabor纹理约束的手指静脉图像修复方法
CN110189339A (zh) 深度图辅助的主动轮廓抠图方法及系统
US20080062171A1 (en) Method for simplifying maintenance of feature of three-dimensional mesh data
Li et al. Refinement of LiDAR point clouds using a super voxel based approach
Schertler et al. Field-aligned online surface reconstruction
CN113870128A (zh) 一种基于深度卷积对抗式网络的数字壁画图像修复方法
CN113393577B (zh) 一种倾斜摄影地形重建方法
Dinesh et al. Local 3D point cloud denoising via bipartite graph approximation & total variation
Laycock et al. Aligning archive maps and extracting footprints for analysis of historic urban environments
Lhuillier Surface reconstruction from a sparse point cloud by enforcing visibility consistency and topology constraints
CN108171790B (zh) 一种基于字典学习的目标重建方法
CN113066016B (zh) 一种基于图信号处理的三维动态点云修复方法
CN113177592A (zh) 一种图像分割方法、装置、计算机设备及存储介质
Zhang et al. Towards unbiased volume rendering of neural implicit surfaces with geometry priors
CN108109205B (zh) 一种基于孔洞填充的残缺模型检索与重建方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant