graph_geometric.hh 5.26 KB
Newer Older
1
2
// graph-tool -- a general graph modification and manipulation thingy
//
Tiago Peixoto's avatar
Tiago Peixoto committed
3
// Copyright (C) 2006-2015 Tiago de Paula Peixoto <tiago@skewed.de>
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 3
// of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.

#ifndef GRAPH_GEOMETRIC_HH
#define GRAPH_GEOMETRIC_HH

#include <iostream>
Tiago Peixoto's avatar
Tiago Peixoto committed
22
#include <unordered_map>
23

24
25
26
#include <boost/functional/hash.hpp>
#include "graph_util.hh"

Tiago Peixoto's avatar
Tiago Peixoto committed
27
#ifndef __clang__
28
29
#include <ext/numeric>
using __gnu_cxx::power;
Tiago Peixoto's avatar
Tiago Peixoto committed
30
31
32
33
34
35
36
#else
template <class Value>
Value power(Value value, int n)
{
    return pow(value, n);
}
#endif
37

38
39
40
41
42
43
namespace graph_tool
{
using namespace std;
using namespace boost;


44
45
46
template <class Point, class Range>
void get_box(const Point& p, double w, vector<int>& box,
             const Range& ranges, bool periodic)
47
48
49
50
{
    if (box.size() != p.size())
        box.resize(p.size());
    for (size_t i = 0; i < p.size(); ++i)
51
52
53
54
55
    {
        box[i] = int(floor(p[i] / w));
        if (periodic && p[i] == ranges[i].second)
            box[i] -= 1;
    }
56
57
}

58
59
60
template <class Point, class Range>
double get_dist(const Point& p1, const Point& p2,
                const Range& ranges, bool periodic)
61
{
62
    double r = 0, diff, size;
63
    for (size_t i = 0; i < p1.size(); ++i)
64
65
66
67
68
69
70
71
72
    {
        diff = abs(p1[i] - p2[i]);
        if (periodic)
        {
            size = abs(ranges[i].second - ranges[i].first);
            diff = min(diff, abs(diff - size));
        }
        r += pow(diff, 2);
    }
73
74
75
    return sqrt(r);
}

76
void periodic(vector<int>& box, const vector<pair<int,int> >& ranges)
77
{
78
    for (size_t i = 0; i < box.size(); ++i)
79
    {
80
81
82
83
        if (box[i] >= ranges[i].second)
            box[i] = ranges[i].first;
        if (box[i] < ranges[i].first)
            box[i] = ranges[i].second - 1;
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
template <class Graph>
bool is_adjacent(typename graph_traits<Graph>::vertex_descriptor v1,
                 typename graph_traits<Graph>::vertex_descriptor v2,
                 Graph& g)
{
    typename graph_traits<Graph>::out_edge_iterator e, e_end;
    for (tie(e, e_end) = out_edges(v1, g); e != e_end; ++e)
        if (target(*e, g) == v2)
            return true;
    return false;
}

struct get_geometric
{
    template <class Graph, class Pos>
    void operator()(Graph& g, Pos upos, vector<vector<double> >& points,
                    vector<pair<double, double> >& ranges,
                    double r, bool periodic_boundary) const
    {
        int N = points.size();
        typename Pos::checked_t pos = upos.get_checked();
        vector<int> box;
110
111
112
113
114
115
116
117
118
119
120
121
        vector<pair<int, int> > box_ranges;

        double w = 2 * r;
        if (periodic_boundary)
        {
            box_ranges.resize(ranges.size());
            for (size_t i = 0; i < ranges.size(); ++i)
            {
                box_ranges[i].first = floor(ranges[i].first / w);
                box_ranges[i].second = ceil(ranges[i].second / w);
            }
        }
122

Tiago Peixoto's avatar
Tiago Peixoto committed
123
        std::unordered_multimap<vector<int>,
124
                                typename graph_traits<Graph>::vertex_descriptor> boxes;
125
126
127
128
129
130

        for (int i = 0; i < N; ++i)
        {
            typename graph_traits<Graph>::vertex_descriptor v = add_vertex(g);
            pos[v].resize(points[i].size());
            copy(points[i].begin(), points[i].end(), pos[v].begin());
131
            get_box(points[i], w, box, ranges, periodic_boundary);
132
133
134
            boxes.insert(make_pair(box, v));
        }

Tiago Peixoto's avatar
Tiago Peixoto committed
135
        int i;
136
        #pragma omp parallel for default(shared) private(i, box) \
137
            schedule(runtime) if (N > 100)
Tiago Peixoto's avatar
Tiago Peixoto committed
138
        for (i = 0; i < N; ++i)
139
140
141
        {
            typename graph_traits<Graph>::vertex_descriptor v = vertex(i, g);

142
            get_box(points[i], w, box, ranges, periodic_boundary);
143
            for (int k = 0; k < power(3, int(box.size())); ++k)
144
            {
145
                for (int j = 0; j < int(box.size()); ++j)
146
                  box[j] += ((k / power(3, j)) % 3) - 1;
147

148
149
150
                if (periodic_boundary)
                    periodic(box, box_ranges);

151
152
153
154
155
156
                typeof(boxes.begin()) iter, end;
                for (tie(iter, end) = boxes.equal_range(box);
                     iter != end; ++iter)
                {
                    typename graph_traits<Graph>::vertex_descriptor w =
                        iter->second;
157
158
                    double d = get_dist(pos[v], pos[w], ranges,
                                        periodic_boundary);
159

160
161
                    if (w > v && d <= r &&
                        (!periodic_boundary || !is_adjacent(v, w, g)))
162
163
164
165
166
                    {
                        #pragma omp critical
                        add_edge(v, w, g);
                    }
                }
167
                get_box(points[i], w, box, ranges, periodic_boundary);
168
169
170
171
172
173
174
175
            }
        }
    }
};

} // namespace graph_tool

#endif // GRAPH_GEOMETRIC_HH