系列文章
《计算几何——算法与应用（第三版）》学习笔记1 - 扫描线（第1-4章）
《计算几何——算法与应用（第三版）》学习笔记2 - 空间索引（第5-7章）
番外篇
计算几何 - C++ 库的调研
-> Boost C++ Libraries学习 - Geometry
Boost C++ Libraries学习 - Boost.Polygon (GTL)

参考 Boost 官方手册

Boost.Geometry contains a dimension-agnostic, coordinate-system-agnostic and scalable kernel, based on concepts, meta-function and tag dispatching. On top of that kernel, algorithms are built: area, length, perimeter, centroid, convex hull, intersection (clipping), within (point in polygon), distance, envelope (bounding box), simplify, transform, and much more. The library supports high precision arithmetic numbers, such as ttmath.

设计原理

Boost 库作为 C++ 中非常常用的库，学习其中的设计原理是对自身的设计能力提高非常有用的。

计算距离的例子

看一个使用自定义点结构计算两点之间的距离的例子，最简单的实现方式如下：

1
2
3
4
5
6
7
8
9
struct mypoint {
  double x, y;
};

double distance(mypoint const& a, mypoint const& b) {
  double dx = a.x - b.x;
  double dy = a.y - b.y;
  return sqrt(dx * dx + dy * dy);
}

计算距离这个需求其实是通用的，如果从库开发者的角度来说，这段代码有如下的问题：

  1. 应该支持任意的“点”概念的结构或类。
  2. 应该不局限于二维。
  3. 应该不局限于笛卡尔坐标系（Cartesian coordinate system）。
  4. 距离概念不止存在于两个点之间，传入任意几何图形的组合也应能计算距离。
  5. double的精度有时不够。
  6. 计算平方根代价相对较大，而且如果是距离比较，绝对的距离不是很重要。

详细修改步骤在这里

解决以上问题的方法：

  1. 使用模板，就可以支持其他“点”了。
  2. 使用 Traits 扩展为维度无关以及坐标类型无关。
  3. 使用 Traits 扩展不同坐标系。
  4. 不同几何之间计算距离，使用 tag dispatching 技术。

使用Boost.Geometry

安装

MacOS 下，使用 homebrew 安装：brew install boost

建立 CMake 工程，并在 CMakeLists.txt 文件中加入以下内容：

1
2
3
4
5
6
find_package(Boost COMPONENTS system filesystem REQUIRED)
find_package(Boost COMPONENTS thread REQUIRED)
if(NOT Boost_FOUND)
    message(FATAL_ERROR "Could not find boost!")
endif()
include_directories(${Boost_INCLUDE_DIRS})

使用

1
2
3
4
5
6
7
8
9
10
11
#include <boost/geometry.hpp>
#include <iostream>

int main(int, char **) {
  boost::geometry::model::d2::point_xy<int> p1(1, 1), p2(2, 2);

  std::cout << "Distance p1-p2 is: " << boost::geometry::distance(p1, p2)
            << std::endl;

  return 0;
}

空间索引（R-Tree）

目前来说，Boost.Geometry 库唯一实现的数据结构是 R-Tree。

TODO: R-Tree 可视化

元素操作

R-Tree 类长这样：

1
2
3
4
5
rtree<Value,
      Parameters,
      IndexableGetter = index::indexable<Value>,
      EqualTo = index::equal_to<Value>,
      Allocator = std::allocator<Value> >

其中：

  • Value: 被存储在容器中的数据类型。
  • Parameters: 参数类型，指明使用的插入/划分算法。
  • IndexableGetter: 函数对象，让 R-Tree 能够理解如何将 Value 转换成 Indexable 对象。
  • EqualTo: 用来比较 Value 的函数对象。
  • Allocator: 创建节点时的申请内存方法。

Indexable 对象是符合 Point, Box 或者 Segment 观念的几何对象。
Value 也可以是 pair 或者 tuple，但需要把 Indexable 对象放在这些数据的第一个。

支持三种算法：

  • Linear: index::rtree< Value, index::linear<16> > rt;
  • Quadratic: index::rtree< Value, index::quadratic<16> > rt;
  • R*-tree: index::rtree< Value, index::rstar<16> > rt;

想要将算法参数在运行时传入 R-Tree，使用下面的写法：

1
2
3
index::rtree<Value, index::dynamic_linear> rt(index::dynamic_linear(16));
index::rtree<Value, index::dynamic_quadratic> rt(index::dynamic_quadratic(16));
index::rtree<Value, index::dynamic_rstar> rt(index::dynamic_rstar(16));

可以使用 R-Tree 中的方法来插入、删除，也可以使用 index:: 中的函数插入、删除。可以插入、删除一个元素，也可以范围操作，具体参考这里。

查询

查询结果要符合一些谓词（predicates）描述，有三种谓词：

  • 空间谓词（spatial predicates）
  • 距离谓词（distance predicates）
  • 自定义谓词（user-defined unary predicate）

查询的写法如下：

1
2
3
std::vector<Value> returned_values;
Box box_region(...);
rt.query(bgi::intersects(box_region), std::back_inserter(returned_values));

由于查询结果一般是用来遍历的，有一种方便的写法（使用|）：

1
2
3
Box box_region(...);
BOOST_FOREACH(Value & v, rt | bgi::adaptors::queried(bgi::intersects(box_region)))
  ; // do something with v

空间查询

不同谓词的效果如下：

intersects(Box) covered_by(Box) disjoint(Box) overlaps(Box) within(Box)
intersects(Segment) intersects(Box) disjoint(Box) intersects(Box) disjoint(Box)

使用如下代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
#include <boost/foreach.hpp>
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/box.hpp>
#include <boost/geometry/geometries/point.hpp>
#include <boost/geometry/index/rtree.hpp>
#include <iostream>
#include <string>
#include <vector>

namespace bg = boost::geometry;
namespace bgi = boost::geometry::index;

typedef bg::model::point<float, 2, bg::cs::cartesian> point;
typedef bg::model::box<point> box;
typedef std::pair<box, unsigned> value;

template <typename T>
void outputQueryResult(const std::string &method, const T &obj,
                       const std::vector<value> &result) {
  std::cout << "---" << std::endl;
  std::cout << "query object:" << std::endl;
  std::cout << bg::wkt<T>(obj) << std::endl;
  std::cout << method << " query result:" << std::endl;
  BOOST_FOREACH (value const &v, result)
    std::cout << bg::wkt<box>(v.first) << " - " << v.second << std::endl;
}

int main(int, char **) {
  bgi::rtree<value, bgi::quadratic<16>> rtree;
  std::vector<value> original_boxes;
  box query_box(point(2.4, 2.4), point(5, 5));
  point query_point(0, 0);

  std::vector<value> result;

  for (unsigned i = 0; i < 10; ++i) {
    box b(point(i + 0.0f, i + 0.0f), point(i + 0.5f, i + 0.5f));
    value v = std::make_pair(b, i);
    rtree.insert(v);
    original_boxes.push_back(v);
  }

  outputQueryResult("original", query_point, original_boxes);

  result.clear();
  rtree.query(bgi::covered_by(query_box), std::back_inserter(result));
  outputQueryResult("covered_by", query_box, result);

  result.clear();
  rtree.query(bgi::covers(query_box), std::back_inserter(result));
  outputQueryResult("covers", query_box, result);

  result.clear();
  rtree.query(bgi::disjoint(query_box), std::back_inserter(result));
  outputQueryResult("disjoint", query_box, result);

  result.clear();
  rtree.query(bgi::intersects(query_box), std::back_inserter(result));
  outputQueryResult("intersects", query_box, result);

  result.clear();
  rtree.query(bgi::overlaps(query_box), std::back_inserter(result));
  outputQueryResult("overlaps", query_box, result);

  result.clear();
  rtree.query(bgi::within(query_box), std::back_inserter(result));
  outputQueryResult("within", query_box, result);

  result.clear();
  rtree.query(bgi::nearest(query_point, 3), std::back_inserter(result));
  outputQueryResult("knn", query_point, result);

  return 0;
}

得到输出：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
---
query object:
POINT(0 0)
original query result:
POLYGON((0 0,0 0.5,0.5 0.5,0.5 0,0 0)) - 0
POLYGON((1 1,1 1.5,1.5 1.5,1.5 1,1 1)) - 1
POLYGON((2 2,2 2.5,2.5 2.5,2.5 2,2 2)) - 2
POLYGON((3 3,3 3.5,3.5 3.5,3.5 3,3 3)) - 3
POLYGON((4 4,4 4.5,4.5 4.5,4.5 4,4 4)) - 4
POLYGON((5 5,5 5.5,5.5 5.5,5.5 5,5 5)) - 5
POLYGON((6 6,6 6.5,6.5 6.5,6.5 6,6 6)) - 6
POLYGON((7 7,7 7.5,7.5 7.5,7.5 7,7 7)) - 7
POLYGON((8 8,8 8.5,8.5 8.5,8.5 8,8 8)) - 8
POLYGON((9 9,9 9.5,9.5 9.5,9.5 9,9 9)) - 9
---
query object:
POLYGON((2.4 2.4,2.4 5,5 5,5 2.4,2.4 2.4))
covered_by query result:
POLYGON((3 3,3 3.5,3.5 3.5,3.5 3,3 3)) - 3
POLYGON((4 4,4 4.5,4.5 4.5,4.5 4,4 4)) - 4
---
query object:
POLYGON((2.4 2.4,2.4 5,5 5,5 2.4,2.4 2.4))
covers query result:
---
query object:
POLYGON((2.4 2.4,2.4 5,5 5,5 2.4,2.4 2.4))
disjoint query result:
POLYGON((0 0,0 0.5,0.5 0.5,0.5 0,0 0)) - 0
POLYGON((1 1,1 1.5,1.5 1.5,1.5 1,1 1)) - 1
POLYGON((6 6,6 6.5,6.5 6.5,6.5 6,6 6)) - 6
POLYGON((7 7,7 7.5,7.5 7.5,7.5 7,7 7)) - 7
POLYGON((8 8,8 8.5,8.5 8.5,8.5 8,8 8)) - 8
POLYGON((9 9,9 9.5,9.5 9.5,9.5 9,9 9)) - 9
---
query object:
POLYGON((2.4 2.4,2.4 5,5 5,5 2.4,2.4 2.4))
intersects query result:
POLYGON((2 2,2 2.5,2.5 2.5,2.5 2,2 2)) - 2
POLYGON((3 3,3 3.5,3.5 3.5,3.5 3,3 3)) - 3
POLYGON((4 4,4 4.5,4.5 4.5,4.5 4,4 4)) - 4
POLYGON((5 5,5 5.5,5.5 5.5,5.5 5,5 5)) - 5
---
query object:
POLYGON((2.4 2.4,2.4 5,5 5,5 2.4,2.4 2.4))
overlaps query result:
POLYGON((2 2,2 2.5,2.5 2.5,2.5 2,2 2)) - 2
---
query object:
POLYGON((2.4 2.4,2.4 5,5 5,5 2.4,2.4 2.4))
within query result:
POLYGON((3 3,3 3.5,3.5 3.5,3.5 3,3 3)) - 3
POLYGON((4 4,4 4.5,4.5 4.5,4.5 4,4 4)) - 4
---
query object:
POINT(0 0)
knn query result:
POLYGON((2 2,2 2.5,2.5 2.5,2.5 2,2 2)) - 2
POLYGON((1 1,1 1.5,1.5 1.5,1.5 1,1 1)) - 1
POLYGON((0 0,0 0.5,0.5 0.5,0.5 0,0 0)) - 0

可视化：

调小查询矩形后的可视化结果：

附上可视化使用的 python 代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
# -*- coding: utf-8 -*-
# author: GuoFan
# time: 2023-9-21

import matplotlib.pyplot as plt
import re

names = ["covered_by", "covers", "disjoint", "intersect", "overlaps", "within"]


def drawPolygons(polygons, color, ax):
    for polygon in polygons:
        print(polygon)
        ax.add_patch(plt.Polygon(xy=polygon, edgecolor=color, fill=False))

    plt.autoscale(True)


def extractPolygon(line):
    result = []
    m = re.match(r"POLYGON\(\((.+)\)\)", line)
    if m:
        coords = m.group(1)
        coord_list = coords.split(",")
        for coord in coord_list:
            xy = coord.split(" ")
            result.append([float(xy[0]), float(xy[1])])

    return result


def main():
    with open("query.txt", "r") as f:
        data = f.read()

    tests_polygons = []

    tests = data.split("---")
    for test in tests:
        lines = test.split("\n")
        polygons = []
        for line in lines:
            polygon = extractPolygon(line)
            if len(polygon) > 0:
                polygons.append(polygon)
        if len(polygons) > 0:
            tests_polygons.append(polygons)

    fig = plt.figure()

    for i in range(6):
        ax = fig.add_subplot(231 + i)
        ax.set_title(names[i])
        drawPolygons(tests_polygons[0], "lightblue", ax)
        drawPolygons(tests_polygons[i + 1], "pink", ax)

    plt.savefig("query.png")
    plt.show()


if __name__ == "__main__":
    main()

如果要使用自定义的谓词，用下面的写法：

1
2
rt.query(index::intersects(box) && index::satisfies([](Value const& v) { return v.is_red(); }),
         std::back_inserter(result));

参考手册

目录