PLaSK library

Loading...
Searching...
No Matches
spatial_index.cpp
Go to the documentation of this file.
1/*
2 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
3 * Copyright (c) 2022 Lodz University of Technology
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, version 3.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 */
14#include "spatial_index.hpp"
15
16 namespace plask {
17
18
20 template <int DIMS>
21 struct GeometryObjectBBox {
22
23 shared_ptr<const Translation<DIMS> > obj;
24
25 typename Primitive<DIMS>::Box boundingBox;
26
27 GeometryObjectBBox() {}
28
29 GeometryObjectBBox(shared_ptr<const Translation<DIMS> > obj, const typename Primitive<DIMS>::Box& boundingBox)
30 : obj(obj), boundingBox(boundingBox) {}
31
32 GeometryObjectBBox(shared_ptr<const Translation<DIMS> > obj)
33 : obj(obj), boundingBox(obj->getBoundingBox()) {}
34
35};
36
37 template <int DIMS>
38 struct EmptyLeafCacheNode: public SpatialIndexNode<DIMS> {
39
40 shared_ptr<Material> getMaterial(const Vec<DIMS>& /*p*/) const override {
41 return shared_ptr<Material>();
42 }
43
44 bool contains(const Vec<DIMS>& /*p*/) const override {
45 return false;
46 }
47
48 GeometryObject::Subtree getPathsAt(shared_ptr<const GeometryObject> /*caller*/, const Vec<DIMS> &/*point*/, bool /*all*/) const override {
49 return GeometryObject::Subtree();
50 }
51};
52
53 template <int DIMS>
54 struct LeafCacheNode: public SpatialIndexNode<DIMS> {
55
57 typedef std::vector< shared_ptr<const Translation<DIMS> > > ChildVectorT;
58
59 ChildVectorT children;
60
61 LeafCacheNode(const std::vector< GeometryObjectBBox<DIMS> >& children_with_bb) {
62 children.reserve(children_with_bb.size());
63 for (const GeometryObjectBBox<DIMS>& c: children_with_bb)
64 children.push_back(c.obj);
65 }
66
67 LeafCacheNode(const std::vector< shared_ptr< Translation<DIMS> > >& childr) {
68 children.reserve(childr.size());
69 for (const shared_ptr< Translation<DIMS> >& c: childr)
70 children.push_back(c);
71 }
72
73 shared_ptr<Material> getMaterial(const Vec<DIMS>& p) const override {
74 for (auto child_it = children.rbegin(); child_it != children.rend(); ++child_it) {
75 shared_ptr<Material> r = (*child_it)->getMaterial(p);
76 if (r != nullptr) return r;
77 }
78 return shared_ptr<Material>();
79 }
80
81 bool contains(const Vec<DIMS>& p) const override {
82 for (auto child: children) if (child->contains(p)) return true;
83 return false;
84 }
85
86 GeometryObject::Subtree getPathsAt(shared_ptr<const GeometryObject> caller, const Vec<DIMS> &point, bool all) const override {
87 GeometryObject::Subtree result;
88 for (auto child = children.rbegin(); child != children.rend(); ++child) {
89 GeometryObject::Subtree child_path = (*child)->getPathsAt(point, all);
90 if (!child_path.empty()) {
91 result.children.push_back(std::move(child_path));
92 if (!all) break;
93 }
94 }
95 if (!result.children.empty())
96 result.object = caller;
97 return result;
98 /*GeometryObject::Subtree result; //general container version:
99 if (all) {
100 for (auto child = children.begin(); child != children.end(); ++child) {
101 GeometryObject::Subtree child_path = (*child)->getPathsAt(point, true);
102 if (!child_path.empty())
103 result.children.push_back(std::move(child_path));
104 }
105 } else {
106 for (auto child = children.rbegin(); child != children.rend(); ++child) {
107 GeometryObject::Subtree child_path = (*child)->getPathsAt(point, false);
108 if (!child_path.empty()) {
109 result.children.push_back(std::move(child_path));
110 break;
111 }
112 }
113 }
114 if (!result.children.empty())
115 result.object = caller;
116 return result;*/
117 }
118};
119
121 template <int DIMS, int dir>
122 struct InternalCacheNode: public SpatialIndexNode<DIMS> {
123
124 double offset;
125 SpatialIndexNode<DIMS>* lo;
126 SpatialIndexNode<DIMS>* hi;
127
128 InternalCacheNode(const double& offset, SpatialIndexNode<DIMS>* lo, SpatialIndexNode<DIMS>* hi)
129 : offset(offset), lo(lo), hi(hi)
130 {}
131
132 shared_ptr<Material> getMaterial(const Vec<DIMS>& p) const override {
133 return p[dir] < offset ? lo->getMaterial(p) : hi->getMaterial(p);
134 }
135
136 bool contains(const Vec<DIMS>& p) const override {
137 return p[dir] < offset ? lo->contains(p) : hi->contains(p);
138 }
139
140 GeometryObject::Subtree getPathsAt(shared_ptr<const GeometryObject> caller, const Vec<DIMS> &p, bool all) const override {
141 return p[dir] < offset ? lo->getPathsAt(caller, p, all) : hi->getPathsAt(caller, p, all);
142 }
143
144 virtual ~InternalCacheNode() {
145 delete lo;
146 delete hi;
147 }
148};
149
150/*template <int DIMS>
151 inline CacheNode<DIMS>* constructInternalNode(int dir, const double& offset, CacheNode<DIMS>* lo, CacheNode<DIMS>* hi) {
152 static_assert(0, "DIMS must be 2 or 3");
153}
154
155 template <>*/
156 inline SpatialIndexNode<2>* constructInternalNode(int dir, const double& offset, SpatialIndexNode<2>* lo, SpatialIndexNode<2>* hi) {
157 if (dir == 0) return new InternalCacheNode<2, 0>(offset, lo, hi);
158 assert(dir == 1);
159 return new InternalCacheNode<2, 1>(offset, lo, hi);
160}
161
162//template <>
163 inline SpatialIndexNode<3>* constructInternalNode(int dir, const double& offset, SpatialIndexNode<3>* lo, SpatialIndexNode<3>* hi) {
164 if (dir == 0) return new InternalCacheNode<3, 0>(offset, lo, hi);
165 if (dir == 1) return new InternalCacheNode<3, 1>(offset, lo, hi);
166 assert(dir == 2);
167 return new InternalCacheNode<3, 2>(offset, lo, hi);
168}
169
170
171 template <int DIMS>
172 void inPlaceSplit(std::vector< GeometryObjectBBox<DIMS> >& inputAndLo, std::vector< GeometryObjectBBox<DIMS> >& hi, int dir, double offset) {
173 std::vector< GeometryObjectBBox<DIMS> > lo;
174 for (GeometryObjectBBox<DIMS>& i: inputAndLo) {
175 if (i.boundingBox.lower[dir] < offset) lo.push_back(i);
176 if (i.boundingBox.upper[dir] >= offset) hi.push_back(i);
177 }
178 std::swap(lo, inputAndLo);
179}
180
187 template <int DIMS>
188 void calcOptimalSplitOffset(const std::vector< GeometryObjectBBox<DIMS> >& inputSortedByLo, const std::vector< GeometryObjectBBox<DIMS> >& inputSortedByHi,
189 int inputDir, int& bestDir, double& bestOffset, int& bestValue)
190{
191 const int max_allowed_size = int(inputSortedByLo.size()) - 4;
192 std::size_t i_hi = 0;
193 for (std::size_t i_lo = 1; i_lo < inputSortedByLo.size(); ++i_lo) {
194 const double& offset = inputSortedByLo[i_lo].boundingBox.lower[inputDir];
195 while (i_lo+1 < inputSortedByLo.size() && inputSortedByLo[i_lo+1].boundingBox.lower[inputDir] == offset)
196 ++i_lo; //can has more obj. with this lo coordinate
197 //now: obj. from [0, i_lo) will be added to lo set
198 if (int(i_lo) > max_allowed_size)
199 return; //too much obj in lo, i_lo will be increased so we can return
200 while (i_hi < inputSortedByHi.size() && inputSortedByHi[i_hi].boundingBox.upper[inputDir] < offset)
201 ++i_hi;
202 //now: obj. from [i_hi, inputSortedByHi.size()) will be added to hi set
203 const int hi_size = int(inputSortedByHi.size()) - int(i_hi);
204 if (hi_size > max_allowed_size)
205 continue; //too much obj in hi, we must wait for higher i_hi
206 //common part is: [i_hi, i_lo)
207 const int value = (int(i_lo) - int(i_hi)) * 3 //this is number of common obj in two sets * 3, we want to minimalize this
208 + std::abs(int(hi_size) - int(i_lo)); //different of set sizes, we also want to minimalize this
209 if (value < bestValue) {
210 bestValue = value;
211 bestOffset = offset;
212 bestDir = inputDir;
213 }
214 }
215}
216
217#define MIN_CHILD_TO_TRY_SPLIT 16
218
229 template <int DIMS>
230 inline SpatialIndexNode<DIMS>* buildCacheR(std::vector< GeometryObjectBBox<DIMS> >* input, int max_depth = 16) {
231 if (input[0].size() < MIN_CHILD_TO_TRY_SPLIT || max_depth == 0) return new LeafCacheNode<DIMS>(input[0]);
232 double bestOffset;
233 int bestDir;
234 int bestValue = std::numeric_limits<int>::max(); //we will minimalize this value
235 for (int dim = DIMS; dim < DIMS; ++dim)
236 calcOptimalSplitOffset(input[dim*2 + 1], input[dim*2 + 2], dim, bestDir, bestOffset, bestValue);
237 if (bestValue == std::numeric_limits<int>::max()) //there are no enough good split point
238 return new LeafCacheNode<DIMS>(input[0]); //so we will not split more
239 SpatialIndexNode<DIMS> *lo, *hi;
240 {
241 std::vector< GeometryObjectBBox<DIMS> > input_over_offset[1 + DIMS * 2];
242 for (int dim = DIMS; dim < 1 + DIMS * 2; ++dim)
243 inPlaceSplit<DIMS>(input[dim], input_over_offset[dim], bestDir, bestOffset);
244 hi = buildCacheR(input_over_offset, max_depth-1);
245 } //here input_over_offset is deleted
246 lo = buildCacheR(input, max_depth-1);
247 return constructInternalNode(bestDir, bestOffset, lo, hi);
248}
249
250 template <int DIMS>
251std::unique_ptr<SpatialIndexNode<DIMS>> buildSpatialIndex(const std::vector< shared_ptr<Translation<DIMS>> >& children) {
252 if (children.empty()) return std::unique_ptr<SpatialIndexNode<DIMS>>(new EmptyLeafCacheNode<DIMS>());
253 if (children.size() < MIN_CHILD_TO_TRY_SPLIT) return std::unique_ptr<SpatialIndexNode<DIMS>>(new LeafCacheNode<DIMS>(children));
254 std::vector< GeometryObjectBBox<DIMS> > input[1 + DIMS * 2];
255 input[0].reserve(children.size());
256 for (auto& c: children) input[0].emplace_back(c);
257 for (int dir = 0; dir < DIMS; ++dir) {
258 input[2*dir + 1] = input[0];
259 std::sort(input[2*dir + 1].begin(), input[2*dir + 1].end(),
260 [dir](const GeometryObjectBBox<DIMS>& a, const GeometryObjectBBox<DIMS>& b) {
261 return a.boundingBox.lower[dir] < b.boundingBox.lower[dir];
262 }
263 );
264 input[2*dir + 2] = input[0];
265 std::sort(input[2*dir + 2].begin(), input[2*dir + 2].end(),
266 [dir](const GeometryObjectBBox<DIMS>& a, const GeometryObjectBBox<DIMS>& b) {
267 return a.boundingBox.upper[dir] < b.boundingBox.upper[dir];
268 }
269 );
270 }
271 return std::unique_ptr<SpatialIndexNode<DIMS>>(buildCacheR<DIMS>(input));
272}
273
274 template struct PLASK_API SpatialIndexNode<2>;
275 template struct PLASK_API SpatialIndexNode<3>;
276
277 template PLASK_API std::unique_ptr<SpatialIndexNode<2>> buildSpatialIndex(const std::vector< shared_ptr<Translation<2>> >& children);
278 template PLASK_API std::unique_ptr<SpatialIndexNode<3>> buildSpatialIndex(const std::vector< shared_ptr<Translation<3>> >& children);
279
280
281} // namespace plask
282