CGAL笔记之形状重建——最优运送曲线重构
1 引言
该软件包实现了重建和简化 2D 点集的方法。输入是一组具有质量属性的 2D 点,可能受到噪声和异常值的阻碍。输出是一组近似输入点的线段和隔离点,如图所示。质量属性与每个点的近似重要性相关。左:受噪声阻碍的输入点设置。右:由线段组成的相应重建形状。
在内部,该算法从所有输入点构建初始 2D Delaunay 三角测量,然后简化三角测量,以便三角测量的边和顶点的子集与输入点非常近似。此处的近似值是指基于最佳运输的稳健距离是如何工作的?使用半边折叠、边翻转和顶点重定位运算符的组合来简化三角测量。三角测量在简化过程中仍然有效,即既没有重叠也没有折叠。
重建算法的输出是三角测量的边和顶点的子集。如图 描述了一个示例,其中输出由绿色边和一个隔离顶点组成。绿色边被认为是相关的,因为它们与许多输入点近似。以灰色描绘的边(称为虚影边并被丢弃)不近似任何输入点。以红色表示的边(称为不相关且也丢弃)近似于某些输入点,但不足以被视为相关。(a) 输入点。(b) 输入点的德劳奈三角测量。(c) 简化后,重影边缘为灰色,相关边缘为绿色,不相关边缘为红色。(d) 由若干边和一个孤立顶点组成的最后重建。
1.1 最简单的例子
下面的示例首先在正方形上生成一组输入点。然后,没有质量属性的点将传递给Optimal_transportation_reconstruction_2对象。初始化后,将执行 100 次重建过程的迭代。
// Simplest example for Optimal_transportation_reconstruction_2, with no mass
// attributes for the input points and no Wasserstein tolerance
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/point_generators_2.h>
#include <CGAL/Optimal_transportation_reconstruction_2.h>
#include <fstream>
#include <iostream>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_2 Point;
typedef CGAL::Optimal_transportation_reconstruction_2<K> Otr;
int main ()
{
// Generate a set of random points on the boundary of a square.
std::vector<Point> points;
CGAL::Random_points_on_square_2<Point> point_generator(1.);
std::copy_n(point_generator, 100, std::back_inserter(points));
Otr otr(points);
if (otr.run(100)) //100 steps
std::cerr << "All done." << std::endl;
else
std::cerr << "Premature ending." << std::endl;
return 0;
}
输出重建的目标边沿数与近似误差概念之间没有直接关系。因此,简化算法可以通过均匀到一定距离的最大公差误差来停止。更具体地说,公差误差被指定为每单位质量运输成本的最大平方根,它与距离均匀。
公差是在瓦瑟斯坦距离的意义上给出的。它不是豪斯多夫公差:这并不意味着输入样本和输出折线之间的距离保证小于 。这意味着每质量运输成本的平方根(均匀到一定距离)最多是tolerance。
// Simplest example with tolerance for
// Optimal_transportation_reconstruction_2, with no mass attributes
// for the input points
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/point_generators_2.h>
#include <CGAL/Optimal_transportation_reconstruction_2.h>
#include <fstream>
#include <iostream>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_2 Point;
typedef CGAL::Optimal_transportation_reconstruction_2<K> Otr;
int main ()
{
// Generate 100 random points on the boundary of a square.
std::vector<Point> points;
CGAL::Random_points_on_square_2<Point> point_generator(1.);
std::copy_n(point_generator, 100, std::back_inserter(points));
Otr otr(points); // no mass given, one unit mass per point assumed
otr.run_under_wasserstein_tolerance(0.1);
return 0;
}
1.2 输出示例
重建的输出可以通过两种方式获得:作为 2D 点和线段的序列,或者作为索引格式,其中线段的连通性被编码,因此称为顶点和边。索引格式记录点列表,然后记录所述列表中的边点索引对,以及隔离顶点的点索引。
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Optimal_transportation_reconstruction_2.h>
#include <fstream>
#include <iostream>
#include <string>
#include <iterator>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::FT FT;
typedef K::Point_2 Point;
typedef K::Segment_2 Segment;
typedef CGAL::Optimal_transportation_reconstruction_2<K> Otr_2;
void load_xy_file(const std::string& filename, std::vector<Point>& points)
{
std::ifstream ifs(filename.c_str());
Point point;
while (ifs >> point)
points.push_back(point);
ifs.close();
}
void list_output(Otr_2& otr2)
{
std::cout << "(-------------List output---------- )" << std::endl;
std::vector<Point> isolated_points;
std::vector<Segment> segments;
otr2.list_output(
std::back_inserter(isolated_points), std::back_inserter(segments));
std::vector<Point>::iterator pit;
for (pit = isolated_points.begin(); pit != isolated_points.end(); pit++)
std::cout << *pit << std::endl;
std::vector<Segment>::iterator sit;
for (sit = segments.begin(); sit != segments.end(); sit++)
std::cout << *sit << std::endl;
}
int main ()
{
std::vector<Point> points;
load_xy_file("data/stair-noise00.xy", points);
Otr_2 otr2(points);
otr2.run(100); // 100 steps
list_output(otr2);
return 0;
}
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Optimal_transportation_reconstruction_2.h>
#include <fstream>
#include <iostream>
#include <string>
#include <iterator>
#include <vector>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::FT FT;
typedef K::Point_2 Point;
typedef CGAL::Optimal_transportation_reconstruction_2<K> Otr_2;
void load_xy_file(const std::string& filename, std::vector<Point>& points)
{
std::ifstream ifs(filename.c_str());
Point point;
while (ifs >> point)
points.push_back(point);
ifs.close();
}
void indexed_output(Otr_2& otr2)
{
std::cout << "(-------------Off output---------- )" << std::endl;
std::vector<Point> points;
std::vector<std::size_t> isolated_vertices;
std::vector<std::pair<std::size_t,std::size_t> > edges;
otr2.indexed_output(
std::back_inserter(points),
std::back_inserter(isolated_vertices),
std::back_inserter(edges));
std::cout << "OFF " << points.size() << " 0 " << edges.size() << std::endl;
// points
std::vector<Point>::iterator pit;
for (pit = points.begin(); pit != points.end(); pit++)
std::cout << *pit << std::endl;
// isolated vertices
std::vector<std::size_t>::iterator vit;
for (vit = isolated_vertices.begin(); vit != isolated_vertices.end(); vit++)
std::cout << "1 " << *vit << std::endl;
// edges
std::vector<std::pair<std::size_t, std::size_t> >::iterator eit;
for (eit = edges.begin(); eit != edges.end(); eit++)
std::cout << "2 " << eit->first << " " << eit->second << std::endl;
}
int main ()
{
std::vector<Point> points;
load_xy_file("data/stair-noise00.xy", points);
Otr_2 otr2(points);
otr2.run(100); // 100 steps
indexed_output(otr2);
return 0;
}
1.3 质量属性示例
下面的示例首先从 ASCII 文件中读取一组输入点和质量。使用两个属性映射,将点及其初始质量传递给Optimal_transportation_reconstruction_2对象。初始化后执行重建过程的100次迭代,然后提取重建形状的线段和隔离点并将其打印到控制台。
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Optimal_transportation_reconstruction_2.h>
#include <fstream>
#include <iostream>
#include <string>
#include <iterator>
#include <utility> // std::pair
#include <vector>
#include <CGAL/property_map.h>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::FT FT;
typedef K::Point_2 Point;
typedef K::Segment_2 Segment;
typedef std::pair<Point, FT> PointMassPair;
typedef std::vector<PointMassPair> PointMassList;
typedef CGAL::First_of_pair_property_map <PointMassPair> Point_property_map;
typedef CGAL::Second_of_pair_property_map <PointMassPair> Mass_property_map;
typedef CGAL::Optimal_transportation_reconstruction_2<
K, Point_property_map, Mass_property_map> Otr_2;
void load_xym_file(const std::string& filename, PointMassList& points)
{
std::ifstream ifs(filename.c_str());
Point point;
FT mass;
while (ifs >> point && ifs >> mass)
points.push_back(std::make_pair(point, mass));
ifs.close();
}
int main ()
{
PointMassList points;
load_xym_file("data/stair.xym", points);
Point_property_map point_pmap;
Mass_property_map mass_pmap;
Otr_2 otr2(points, point_pmap, mass_pmap);
otr2.run(100); // 100 steps
std::vector<Point> isolated_vertices;
std::vector<Segment> edges;
otr2.list_output(
std::back_inserter(isolated_vertices), std::back_inserter(edges));
std::cout << "Isolated vertices:" << std::endl;
std::vector<Point>::iterator vit;
for (vit = isolated_vertices.begin(); vit != isolated_vertices.end(); vit++)
std::cout << *vit << std::endl;
std::cerr << "Edges:" << std::endl;
std::vector<Segment>::iterator eit;
for (eit = edges.begin(); eit != edges.end(); eit++)
std::cout << *eit << std::endl;
return 0;
}
2 接口
向用户公开的唯一类是Optimal_transportation_reconstruction_2
类。
2.1 示例调用
/*
K : a geometric kernel.
*/
Optimal_transportation_reconstruction_2<K>
otr2(points.begin(), points.end());
otr2.run(100); // perform 100 edge collapse operators
如果输入不仅仅是没有质量的点,则可以提供与此输入匹配的属性映射。
/*
K : a geometric kernel.
Point_property_map : a PropertyMap for accessing the input points.
Mass_property_map : a PropertyMap for accessing the mass attributes of the
input points.
*/
Optimal_transportation_reconstruction_2<K, Point_property_map, Mass_property_map>
otr2(points.begin(), points.end(), point_pmap, mass_pmap);
otr2.run(100); // perform 100 edge collapse operators
除了调用 run()
之外,还可以调用 run_until()
并指定要保留的输出顶点数,如图所示。
otr2.run_until(20); // perform edge collapse operators until 20 vertices remain.
如图分别由 2000、400 和 200 个输入点组成的数据集中的 20 顶点重建示例。这些示例说明了输入点密度降低时算法的行为。
最后,调用 [run_under_wasserstein_tolerance()
允许用户在输出边数未知时基于距离标准运行算法,如图 所示。
otr2.run_under_wasserstein_tolerance(0.1); // perform collapses until no more can be done within tolerance
具有不同瓦瑟斯坦公差阈值的重建示例。顶部:输入点设置和重建,公差为 0.005。底部:公差为 0.05 且公差为 0.1 的重建。
2.2 全局点重定位
由于噪声和缺失数据可能会阻止重建的形状在正确的位置具有尖角,因此该算法提供了一个功能来重新定位重建的所有点:
otr2.relocate_all_points();
请注意,这些点与基础三角测量的顶点一致。此函数可以在一次简化运行后调用,也可以与多次简化运行交错调用。
选择新的点位置,以便改进输出线段和隔离点与输入点的近似值。更具体地说,搬迁过程在计算给定当前重建的最佳运输计划与在保持当前运输计划不变的情况下重新定位三角测量顶点之间进行迭代。顶点被重新定位,以尽量减少当前运输计划引起的运输成本。
3 用户参数
算法的行为通过以下参数进行控制。
3.1 边缘翻转
在简化内部三角测量期间,需要一些递归边翻转运算符来保证三角测量在应用半边折叠运算符时仍然有效。调用 set_use_flip(false)
可防止算法使用边缘翻转,从而以次优结果为代价产生更短的计算时间,因为并非所有边缘都可以被视为可折叠的。边缘翻转。左:蓝色边缘由于以黑色显示的阻挡边缘而产生折叠。中间:运行递归边缘翻转过程后,蓝色边缘是可折叠的。右:边缘塌陷后的三角测量。
3.2 边缘相关性
从近似角度来看,如果 (1) 边很长,(2) 近似点很多(或点具有质量属性时有大量质量),并且 (3) 具有较小的近似误差,则边是相关的。
3.3 随机样本量
默认情况下,简化依赖于抽取期间半边折叠运算符的穷举优先级队列。为了提高效率,严格大于 0 的参数样本数量将切换到多项选择方法,即,在样本大小的边缘折叠运算符的随机样本中选择最佳选择。样本数量的典型值为 15,但在针对非常粗略的简化时,必须放大此值。
3.4 本地点迁移
除了上述全局重定位函数之外, Optimal_transportation_reconstruction_2
类构造函数的可选参数提供了一种在每个边折叠运算符(可能与边翻转结合使用)之后本地重定位点的方法。本文中的局部意味着仅重新定位每个边折叠算子周围的三角测量中的局部模板的顶点,其过程类似于上述全局重定位函数中描述的过程。局部模具被选为折点在折叠边后剩余的单环邻域。重新定位过程是迭代的,一个参数控制重新定位步骤的数量。
3.5 详细输出
介于 0 和 2 之间的详细参数确定算法生成的控制台输出量。0 值不会生成标准输出的输出。大于 0 的值将生成标准输出std::cerr
的输出。
4 稳健性
该算法的优点是它对噪声和异常值的鲁棒性。如图显示,只要异常值的密度小于输入点的密度,该算法的输出是鲁棒的,它几乎不受噪声和/或异常值的影响。对噪声和异常值的鲁棒性。左:无噪声点设置。中间:嘈杂的点集。右:受噪声和异常值阻碍的点集。
5 可变密度
下图说明了算法在具有均匀质量属性与可变密度的点集上的行为。由于该算法更加重视密集采样区域,因此这将转化为密集采样区域上的较小边缘。在稀疏采样的区域中,算法最初通过一个孤立的顶点近似每个点,然后逐渐用边近似这些点。如图可变密度。左:输入点集。其他三张图片显示了在推动简化时近似值是如何演变的。
6 混合维度
下图 描述了一个输入点集,对一组线段和一个实心区域进行采样。根据输出中的目标点数,实体面积由一组均匀采样的孤立顶点近似。如图混合维度。左:输入点集。中间:孤立的顶点为蓝色,相关边为绿色,不相关边为红色。右:最终输出。
7 可变质量
质量属性提供了一种调整每个点的重要性以进行近似的方法。下图描述了阈值化后灰度级图像的重建,其中像素的灰度用作质量属性。如图可变质量。左:输入灰度图像。中间:阈值化后的图像,以减少用作非零质量点的像素数。右:最终重建。