-
Notifications
You must be signed in to change notification settings - Fork 71
/
2014-12-27-Rtree-examples.cpp
183 lines (168 loc) · 6.57 KB
/
2014-12-27-Rtree-examples.cpp
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
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
/**
* @title Create an R-tree data structure using Rcpp and Boost::Geometry
* @author Nagi Teramo and Dirk Eddelbuettel
* @license GPL (>= 2)
* @mathjax false
* @tags boost
* @summary We show the use of Boost::Geometry and R-tree objects.
*
* ### Introduction
*
* The purpose of this post is to show how to use
* [Boost::Geometry library](http://www.boost.org/doc/libs/1_57_0/libs/geometry/doc/html/index.html)
* which was introduced
* [recently](http://dirk.eddelbuettel.com/blog/2014/12/22/#bh_1.55.0-1) in
* Rcpp. Especially, we focus on
* [R-tree data structure](http://www.boost.org/doc/libs/1_57_0/libs/geometry/doc/html/geometry/reference/spatial_indexes/boost__geometry__index__rtree.html)
* for searching objects in space because only one spatial index is
* implemented - R-tree Currently in this library.
*
* Boost.Geometry which is part of the Boost C++ Libraries gives us algorithms
* for solving geometry problems. In this library, the Boost.Geometry.Index
* which is one of components is intended to gather data structures called
* spatial indexes which are often used to searching for objects in space
* quickly. Generally speaking, spatial indexes stores *geometric objects'
* representations* and allows searching for objects occupying some space or
* close to some point in space.
*
* ### R-tree
*
* R-tree is a tree data structure used for spatial searching, i.e., for
* indexing multi-dimensional information such as geographical coordinates,
* rectangles or polygons. R-tree was proposed by Antonin Guttman in 1984 as an
* expansion of B-tree for multi-dimensional data and plays significant role in
* both theoretical and applied contexts. It is only one spatial index
* implemented in Boost::Geometry.
*
* As a real application of this, It is often used to store spatial objects such
* as restaurant locations or the polygons which typical maps are made of:
* streets, buildings, outlines of lakes, coastlines, etc in order to perform a
* spatial query like "Find all stations within 1 km of my current location",
* "Let me know all road segments in 2 km of my location" or "find the nearest
* gas station" which we often ask google seach by your voice recenlty. In this
* way, the R-tree can be used (nearest neighbor) search for some places.
*
* You can find more explanations about R-tree in
* [Wikipedia](http://en.wikipedia.org/wiki/R-tree).
*
* ### Write a wrapper class of rtree in Boost::Geometry using Rcpp
*
* Now, we write a simple C++ wrapper class of rtree class in
* Boost::Geometry::Index that we can use in R.
*
* The most important feature to mention here is the use of Rcpp module to
* expose your own class to R. Although almost all classes in Boost library have
* a lot of functions, , you do not use all in many cases. In that case, you
* should write your wrapper class for making your code simple.
*
* ### Rcpp code
*/
// [[Rcpp::depends(BH)]]
// Enable C++11 via this plugin to suppress 'long long' errors
// [[Rcpp::plugins("cpp11")]]
#include <vector>
#include <Rcpp.h>
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point.hpp>
#include <boost/geometry/geometries/box.hpp>
#include <boost/geometry/index/rtree.hpp>
using namespace Rcpp;
// Mnemonics
namespace bg = boost::geometry;
namespace bgi = boost::geometry::index;
typedef bg::model::point<float, 2, bg::cs::cartesian> point_t;
typedef bg::model::box<point_t> box;
typedef std::pair<box, unsigned int> value_t;
class RTreeCpp {
public:
// Constructor.
// You have to give spatial data as a data frame.
RTreeCpp(DataFrame df) {
int size = df.nrows();
NumericVector id = df[0];
NumericVector bl_x = df[1];
NumericVector bl_y = df[2];
NumericVector ur_x = df[3];
NumericVector ur_y = df[4];
for(int i = 0; i < size; ++i) {
// create a box
box b(point_t(bl_x[i], bl_y[i]), point_t(ur_x[i], ur_y[i]));
// insert new value
rtree_.insert(std::make_pair(b, static_cast<unsigned int>(id[i])));
}
}
// This method(query) is k-nearest neighbor search.
// It returns some number of values nearest to some point(point argument) in space.
std::vector<int> knn(NumericVector point, unsigned int n) {
std::vector<value_t> result_n;
rtree_.query(bgi::nearest(point_t(point[0], point[1]), n), std::back_inserter(result_n));
std::vector<int> indexes;
std::vector<value_t>::iterator itr;
for ( itr = result_n.begin(); itr != result_n.end(); ++itr ) {
value_t value = *itr;
indexes.push_back( value.second );
}
return indexes;
}
private:
// R-tree can be created using various algorithm and parameters
// You can change the algorithm by template parameter.
// In this example we use quadratic algorithm.
// Maximum number of elements in nodes in R-tree is set to 16.
bgi::rtree<value_t, bgi::quadratic<16> > rtree_;
};
// [[Rcpp::export]]
std::vector<int> showKNN(Rcpp::DataFrame df, NumericVector point, unsigned int n) {
RTreeCpp tree(df); // recreate tree each time
return tree.knn(point, n);
}
/**
* ### R code using RTreeCpp
*
* First, we create a sample data set of spatial data.
*/
/*** R
# Sample spatial data(boxes)
points <- data.frame(
id=0:2,
bl_x=c(0, 2, 4),
bl_y=c(0, 2, 4),
ur_x=c(1, 3, 5),
ur_y=c(1, 3, 5))
*/
/*
* To visually the data, we use the following code:
*/
/*** R
size <- nrow(points)
#colors for rectangle area
colors <- rainbow(size)
#Plot these points
plot(c(0, 5), c(0, 5), type= "n", xlab="", ylab="")
for(i in 1:size){
rect(points[i, "bl_x"], points[i, "bl_y"], points[i, "ur_x"], points[i, "ur_y"], col=colors[i])
}
legend(4, 2, legend=points$id, fill=colors)
*/
/**
* One can use the RTreeCpp class as follows:
*/
/*** R
# new RTreeCpp object
# Search nearest neighbor points(return value : id of points data)
showKNN(points, c(0, 0), 1)
showKNN(points, c(0, 0), 2)
showKNN(points, c(0, 0), 3)
*/
/**
* Note the re-creation of the `RTreeCpp` object is of course
* inefficient, but the Rcpp Gallery imposes some constraints on how we
* present code. For actual application a stateful and persistent
* object would be created. This could be done via Rcpp Modules as
* well a number of different ways. Here, however, we need to
* recreate the object for each call as `knitr` (which is used behind
* the scenes) cannot persist objects between code chunks. This is
* simply a technical limitation of the Rcpp Gallery---but not of Rcpp
* itself.
*/