48 typedef TCoordType coord_type;
51 typedef std::pair<point_type, point_index> value_type;
52 typedef std::vector<point_index> point_data_vec;
53 typedef point_data_vec::const_iterator pd_cit;
55 typedef CDT::array<node_index, 2> children_type;
66 data.reserve(NumVerticesInLeaf);
85 -std::numeric_limits<coord_type>::max(),
86 -std::numeric_limits<coord_type>::max()))
88 std::numeric_limits<coord_type>::max(),
89 std::numeric_limits<coord_type>::max()))
91 , m_isRootBoxInitialized(false)
92 , m_tasksStack(InitialStackDepth, NearestTask())
94 m_root = addNewNode();
103 , m_isRootBoxInitialized(true)
104 , m_tasksStack(InitialStackDepth, NearestTask())
106 m_root = addNewNode();
119 insert(
const point_index& iPoint,
const std::vector<point_type>& points)
124 while(!isInsideBox(pos, m_min, m_max))
129 node_index node = m_root;
132 NodeSplitDirection::Enum dir = m_rootDir;
135 NodeSplitDirection::Enum newDir(NodeSplitDirection::X);
140 if(m_nodes[node].isLeaf())
143 point_data_vec& pd = m_nodes[node].data;
144 if(pd.size() < NumVerticesInLeaf)
146 pd.push_back(iPoint);
150 if(!m_isRootBoxInitialized)
152 initializeRootBox(points);
157 calcSplitInfo(min, max, dir, mid, newDir, newMin, newMax);
158 const node_index c1 = addNewNode(), c2 = addNewNode();
159 Node& n = m_nodes[node];
161 point_data_vec& c1data = m_nodes[c1].data;
162 point_data_vec& c2data = m_nodes[c2].data;
164 for(pd_cit it = n.
data.begin(); it != n.
data.end(); ++it)
166 whichChild(points[*it], mid, dir) == 0
167 ? c1data.push_back(*it)
168 : c2data.push_back(*it);
170 n.
data = point_data_vec();
174 calcSplitInfo(min, max, dir, mid, newDir, newMin, newMax);
177 const std::size_t iChild = whichChild(points[iPoint], mid, dir);
178 iChild == 0 ? max = newMax : min = newMin;
179 node = m_nodes[node].children[iChild];
190 const std::vector<point_type>& points)
const
194 coord_type minDistSq = std::numeric_limits<coord_type>::max();
195 m_tasksStack[++iTask] =
196 NearestTask(m_root, m_min, m_max, m_rootDir, minDistSq);
199 const NearestTask t = m_tasksStack[iTask--];
200 if(t.distSq > minDistSq)
202 const Node& n = m_nodes[t.node];
205 for(pd_cit it = n.
data.begin(); it != n.
data.end(); ++it)
209 if(distSq < minDistSq)
220 NodeSplitDirection::Enum newDir;
222 calcSplitInfo(t.min, t.max, t.dir, mid, newDir, newMin, newMax);
224 const coord_type distToMid = t.dir == NodeSplitDirection::X
227 const coord_type toMidSq = distToMid * distToMid;
229 const std::size_t iChild = whichChild(point, mid, t.dir);
230 if(iTask + 2 >=
static_cast<int>(m_tasksStack.size()))
233 m_tasksStack.size() + StackDepthIncrement);
238 m_tasksStack[++iTask] = NearestTask(
239 n.
children[1], newMin, t.max, newDir, toMidSq);
240 m_tasksStack[++iTask] = NearestTask(
241 n.
children[0], t.min, newMax, newDir, toMidSq);
245 m_tasksStack[++iTask] = NearestTask(
246 n.
children[0], t.min, newMax, newDir, toMidSq);
247 m_tasksStack[++iTask] = NearestTask(
248 n.
children[1], newMin, t.max, newDir, toMidSq);
257 node_index addNewNode()
259 const node_index newNodeIndex =
static_cast<node_index
>(m_nodes.size());
260 m_nodes.push_back(Node());
266 std::size_t whichChild(
267 const point_type& point,
268 const coord_type& split,
269 const NodeSplitDirection::Enum dir)
const
271 return static_cast<size_t>(
272 dir == NodeSplitDirection::X ? point.x > split : point.y > split);
276 static void calcSplitInfo(
277 const point_type& min,
278 const point_type& max,
279 const NodeSplitDirection::Enum dir,
281 NodeSplitDirection::Enum& newDirOut,
282 point_type& newMinOut,
283 point_type& newMaxOut)
289 case NodeSplitDirection::X:
290 midOut = (min.x + max.x) / coord_type(2);
291 newDirOut = NodeSplitDirection::Y;
292 newMinOut.x = midOut;
293 newMaxOut.x = midOut;
295 case NodeSplitDirection::Y:
296 midOut = (min.y + max.y) / coord_type(2);
297 newDirOut = NodeSplitDirection::X;
298 newMinOut.y = midOut;
299 newMaxOut.y = midOut;
305 static bool isInsideBox(
307 const point_type& min,
308 const point_type& max)
310 return p.x >= min.x && p.x <= max.x && p.y >= min.y && p.y <= max.y;
315 void extendTree(
const point_type& point)
317 const node_index newRoot = addNewNode();
318 const node_index newLeaf = addNewNode();
321 case NodeSplitDirection::X:
322 m_rootDir = NodeSplitDirection::Y;
323 point.y < m_min.
y ? m_nodes[newRoot].setChildren(newLeaf, m_root)
324 : m_nodes[newRoot].setChildren(m_root, newLeaf);
325 if(point.y < m_min.
y)
326 m_min.
y -= m_max.
y - m_min.
y;
327 else if(point.y > m_max.
y)
328 m_max.
y += m_max.
y - m_min.
y;
330 case NodeSplitDirection::Y:
331 m_rootDir = NodeSplitDirection::X;
332 point.x < m_min.
x ? m_nodes[newRoot].setChildren(newLeaf, m_root)
333 : m_nodes[newRoot].setChildren(m_root, newLeaf);
334 if(point.x < m_min.
x)
335 m_min.
x -= m_max.
x - m_min.
x;
336 else if(point.x > m_max.
x)
337 m_max.
x += m_max.
x - m_min.
x;
344 void initializeRootBox(
const std::vector<point_type>& points)
346 const point_data_vec& data = m_nodes[m_root].data;
347 m_min = points[data.front()];
349 for(pd_cit it = data.begin(); it != data.end(); ++it)
351 const point_type& p = points[*it];
353 std::min(m_min.
x, p.x), std::min(m_min.
y, p.y));
355 std::max(m_max.
x, p.x), std::max(m_max.
y, p.y));
359 const TCoordType padding(1);
360 if(m_min.
x == m_max.
x)
365 if(m_min.
y == m_max.
y)
370 m_isRootBoxInitialized =
true;
375 std::vector<Node> m_nodes;
376 NodeSplitDirection::Enum m_rootDir;
381 bool m_isRootBoxInitialized;
388 NodeSplitDirection::Enum dir;
393 const node_index node_,
394 const point_type& min_,
395 const point_type& max_,
396 const NodeSplitDirection::Enum dir_,
397 const coord_type distSq_)
406 mutable std::vector<NearestTask> m_tasksStack;