Line data Source code
1 : /*
2 : * Copyright (c) 2009-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Sebastian Reiter
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__LIB_GRID__KD_TREE_IMPL__
34 : #define __H__LIB_GRID__KD_TREE_IMPL__
35 :
36 : #include <list>
37 : #include <vector>
38 : #include "lib_grid/grid/grid.h"
39 : #include "common/math/ugmath.h"
40 :
41 : namespace ug
42 : {
43 :
44 : typedef std::list<KDVertexDistance> KDVertexDistanceList;
45 :
46 : template<class TPositionAttachment, int numDimensions, class TVector>
47 : void
48 0 : KDTreeStatic<TPositionAttachment, numDimensions, TVector>
49 : ::Node
50 : ::clear()
51 : {
52 0 : if(m_pChild[0])
53 0 : delete m_pChild[0];
54 0 : if(m_pChild[1])
55 0 : delete m_pChild[1];
56 0 : m_pChild[0] = m_pChild[1] = NULL;
57 0 : if(m_pvVertices)
58 0 : delete m_pvVertices;
59 0 : m_pvVertices = NULL;
60 0 : }
61 :
62 : template<class TPositionAttachment, int numDimensions, class TVector>
63 : void
64 : KDTreeStatic<TPositionAttachment, numDimensions, TVector>::
65 : clear()
66 : {
67 0 : m_pGrid = NULL;
68 0 : m_parentNode.clear();
69 : }
70 :
71 : template<class TPositionAttachment, int numDimensions, class TVector>
72 : template <class TVrtIterator>
73 : bool
74 : KDTreeStatic<TPositionAttachment, numDimensions, TVector>::
75 : create_from_grid(Grid& grid, TVrtIterator vrtsBegin, TVrtIterator vrtsEnd,
76 : TPositionAttachment& aPos, int maxTreeDepth,
77 : int splitThreshold, KDSplitDimension splitDimension)
78 : {
79 : Grid::VertexAttachmentAccessor<TPositionAttachment> aaPos(grid, aPos);
80 : return create_from_grid(grid, vrtsBegin, vrtsEnd, aaPos, maxTreeDepth,
81 : splitThreshold, splitDimension);
82 : }
83 :
84 : template<class TPositionAttachment, int numDimensions, class TVector>
85 : template <class TVrtIterator>
86 : bool
87 0 : KDTreeStatic<TPositionAttachment, numDimensions, TVector>::
88 : create_from_grid(Grid& grid, TVrtIterator vrtsBegin, TVrtIterator vrtsEnd,
89 : Grid::VertexAttachmentAccessor<TPositionAttachment> aaPos,
90 : int maxTreeDepth, int splitThreshold,
91 : KDSplitDimension splitDimension)
92 : {
93 : clear();
94 0 : m_pGrid = &grid;
95 0 : m_aaPos = aaPos;
96 0 : m_iSplitThreshold = splitThreshold;
97 0 : m_splitDimension = splitDimension; // how the split dimensions are chosen
98 0 : return create_barycentric(vrtsBegin, vrtsEnd, grid.num_vertices(), &m_parentNode, 0, maxTreeDepth);
99 : }
100 :
101 : template<class TPositionAttachment, int numDimensions, class TVector>
102 : bool
103 0 : KDTreeStatic<TPositionAttachment, numDimensions, TVector>::
104 : get_neighbourhood(std::vector<Vertex*>& vrtsOut, typename TPositionAttachment::ValueType& pos, int numClosest)
105 : {
106 : vrtsOut.clear();
107 0 : m_numNeighboursFound = 0;
108 : KDVertexDistanceList pdList;
109 :
110 0 : neighbourhood(pdList, &m_parentNode, pos, numClosest);
111 0 : for(KDVertexDistanceList::iterator iter = pdList.begin(); iter != pdList.end(); iter++)
112 0 : vrtsOut.push_back((*iter).vertex);
113 0 : return true;
114 : }
115 :
116 : template<class TPositionAttachment, int numDimensions, class TVector>
117 : bool
118 : KDTreeStatic<TPositionAttachment, numDimensions, TVector>::
119 : get_points_in_box(std::vector<Vertex*>& vrtsOut, const TVector& boxMin, const TVector& boxMax)
120 : {
121 : vrtsOut.clear();
122 : return get_vertices_in_box(vrtsOut, &m_parentNode, boxMin, boxMax);
123 : }
124 :
125 : template<class TPositionAttachment, int numDimensions, class TVector>
126 : void
127 : KDTreeStatic<TPositionAttachment, numDimensions, TVector>::
128 : get_leafs(std::vector<Node*>& vLeafsOut)
129 : {
130 : vLeafsOut.clear();
131 : get_leafs_recursive(vLeafsOut, &m_parentNode);
132 : }
133 :
134 :
135 : ////////////////////////////////////////////////////////////////////////
136 : // protected
137 : template<class TPositionAttachment, int numDimensions, class TVector>
138 : bool
139 : KDTreeStatic<TPositionAttachment, numDimensions, TVector>::
140 : get_points_in_box(std::vector<Vertex*>& vrtsOut, Node* pNode, const TVector& boxMin, const TVector& boxMax)
141 : {
142 : // check if we have reached a leaf
143 : if(pNode->m_pvVertices)
144 : {
145 : // we're in a leaf -> add the points that are inseide the box to pointsOut
146 : for(VertexVec::iterator iter = pNode->m_pvVertices->begin(); iter != pNode->m_pvVertices->end(); iter++)
147 : {
148 : bool bAdd = true;
149 : for(int i = 0; i < numDimensions; ++i)
150 : {
151 : if((m_aaPos[*iter].coord(i) < boxMin[i]) || (m_aaPos[*iter].coord(i) > boxMax[i]))
152 : {
153 : bAdd = false;
154 : break;
155 : }
156 : }
157 : if(bAdd)
158 : vrtsOut.push_back(*iter);
159 : }
160 : // done
161 : return true;
162 : }
163 :
164 :
165 : // we haven't reached a leaf yet
166 : // check in wich subnodes we have to descent
167 : if(pNode->m_pChild[0] || pNode->m_pChild[1])
168 : {
169 : bool bSuccess = true;
170 :
171 : if(boxMin[pNode->m_iSplitDimension] < pNode->m_fSplitValue)
172 : bSuccess &= get_points_in_box(vrtsOut, pNode->m_pChild[1], boxMin, boxMax);
173 : if(boxMax[pNode->m_iSplitDimension] >= pNode->m_fSplitValue)
174 : bSuccess &= get_points_in_box(vrtsOut, pNode->m_pChild[0], boxMin, boxMax);
175 :
176 : return bSuccess;
177 : }
178 : return false;
179 : }
180 :
181 : template<class TPositionAttachment, int numDimensions, class TVector>
182 : void
183 0 : KDTreeStatic<TPositionAttachment, numDimensions, TVector>::
184 : neighbourhood(KDVertexDistanceList& vrtsOut, Node* pNode, TVector& pos, int numClosest)
185 : {
186 : // if there are points in this node add them to the list (if they match certain criterias)
187 0 : if(pNode->m_pvVertices)
188 : {
189 : // loop through the points associated with the current node
190 0 : for(VertexVec::iterator iter = pNode->m_pvVertices->begin(); iter != pNode->m_pvVertices->end(); iter++)
191 : {
192 : // calculate the distance of the current point and the test point
193 : float distSQ = 0;
194 0 : for(int i = 0; i < numDimensions; ++i)
195 : {
196 0 : float t = pos[i] - m_aaPos[*iter].coord(i);
197 0 : distSQ += (t*t);
198 : }
199 :
200 : int bInsert = 0;
201 0 : if(m_numNeighboursFound == numClosest)
202 : {
203 0 : if(distSQ < m_maxDistSQ)
204 : bInsert = 2;
205 : }
206 : else
207 : bInsert = 1;
208 :
209 : if(bInsert)
210 : {
211 : KDVertexDistanceList::iterator pdIter;
212 : // first we will insert the point into pointlist
213 0 : for(pdIter = vrtsOut.begin(); pdIter != vrtsOut.end(); pdIter++)
214 : {
215 0 : if(distSQ < (*pdIter).distSQ)
216 : {
217 0 : pdIter = vrtsOut.insert(pdIter, KDVertexDistance(*iter, distSQ));
218 0 : m_numNeighboursFound++;
219 0 : break;
220 : }
221 : }
222 :
223 0 : if(pdIter == vrtsOut.end())
224 : {
225 0 : vrtsOut.push_back(KDVertexDistance(*iter, distSQ));
226 0 : m_numNeighboursFound++;
227 : }
228 0 : m_maxDistSQ = vrtsOut.back().distSQ;
229 :
230 0 : if(m_numNeighboursFound > numClosest)
231 : {
232 : // erase the first one
233 : KDVertexDistanceList::iterator tmpIter = vrtsOut.end();
234 : tmpIter--;
235 : vrtsOut.erase(tmpIter);
236 0 : m_numNeighboursFound--;
237 : // find the new maxDistSQ
238 0 : m_maxDistSQ = vrtsOut.back().distSQ;
239 :
240 : }
241 : }
242 : }
243 : }
244 : // if the node has children visit them
245 0 : if(pNode->m_pChild[0] || pNode->m_pChild[1]) //either both or none are NULL
246 : {
247 : // check in wich subnode the specified point lies.
248 : int bestNodeIndex;
249 0 : if(pos[pNode->m_iSplitDimension] >= pNode->m_fSplitValue)
250 : bestNodeIndex = 0;
251 : else
252 : bestNodeIndex = 1;
253 :
254 : // first we will visit the better ranked node (just a guess)
255 0 : if(pNode->m_pChild[bestNodeIndex])
256 0 : neighbourhood(vrtsOut, pNode->m_pChild[bestNodeIndex], pos, numClosest);
257 :
258 : // if we found less than numClosest we got to continue searching anyway
259 : // else if we may find points on the other side of the splitting plane that are closer than the points we found till now we'll have to visit the worse ranked node too.
260 0 : if(pNode->m_pChild[1 - bestNodeIndex])
261 : {
262 0 : if(m_numNeighboursFound < numClosest)
263 : neighbourhood(vrtsOut, pNode->m_pChild[1 - bestNodeIndex], pos, numClosest);
264 : else
265 : {
266 0 : float t = pos[pNode->m_iSplitDimension] - pNode->m_fSplitValue;
267 0 : if(t*t < m_maxDistSQ)
268 : neighbourhood(vrtsOut, pNode->m_pChild[1 - bestNodeIndex], pos, numClosest);
269 : }
270 : }
271 : }
272 0 : }
273 :
274 : template<class TPositionAttachment, int numDimensions, class TVector>
275 : template <class TVertexIterator>
276 : bool
277 0 : KDTreeStatic<TPositionAttachment, numDimensions, TVector>::
278 : create_barycentric(TVertexIterator vrts_begin, TVertexIterator vrts_end, int numVertices,
279 : Node* pNode, int actDimension, int maxTreeDepth)
280 : {
281 : // check if we are in a leaf
282 : {
283 0 : if((maxTreeDepth < 1) || (numVertices <= m_iSplitThreshold))
284 : {
285 : // we are
286 0 : pNode->m_pvVertices = new VertexVec;
287 : // loop through the points and add them to the node
288 0 : for(TVertexIterator iter = vrts_begin; iter != vrts_end; iter++)
289 0 : pNode->m_pvVertices->push_back(*iter);
290 : // we're done
291 0 : return true;
292 : }
293 : }
294 :
295 : // loop through the points and calculate the barycentre
296 : float barycentre = 0;
297 : {
298 0 : for(TVertexIterator iter = vrts_begin; iter != vrts_end; iter++)
299 0 : barycentre += m_aaPos[*iter].coord(actDimension);
300 0 : barycentre /= (float)numVertices;
301 : }
302 :
303 : // fill the lists for the poitive and negative subnodes
304 : std::list<Vertex*> lstPos;
305 : std::list<Vertex*> lstNeg;
306 : int numPos = 0;
307 : int numNeg = 0;
308 : {
309 0 : for(TVertexIterator iter = vrts_begin; iter != vrts_end; iter++, numPos++, numNeg++)
310 : {
311 0 : if(m_aaPos[*iter].coord(actDimension) >= barycentre)
312 0 : lstPos.push_back(*iter);
313 : else
314 0 : lstNeg.push_back(*iter);
315 : }
316 : }
317 : // create the subnodes
318 : bool bSuccess = true;
319 : {
320 0 : pNode->m_iSplitDimension = actDimension;
321 0 : pNode->m_fSplitValue = barycentre;
322 :
323 0 : for(int i = 0; i < 2; ++i)
324 0 : pNode->m_pChild[i] = new Node;
325 : // the positive one
326 0 : if(!lstPos.empty())
327 0 : bSuccess &= create_barycentric(lstPos.begin(), lstPos.end(), numPos, pNode->m_pChild[0], get_next_split_dimension(actDimension, lstPos.begin(), lstPos.end()), maxTreeDepth - 1);
328 : // the negative one
329 0 : if(!lstNeg.empty())
330 0 : bSuccess &= create_barycentric(lstNeg.begin(), lstNeg.end(), numNeg, pNode->m_pChild[1], get_next_split_dimension(actDimension, lstNeg.begin(), lstNeg.end()), maxTreeDepth - 1);
331 : }
332 : // done
333 : return bSuccess;
334 : }
335 :
336 : template<class TPositionAttachment, int numDimensions, class TVector>
337 : template <class TVertexIterator>
338 : int
339 0 : KDTreeStatic<TPositionAttachment, numDimensions, TVector>::
340 : get_largest_dimension(TVertexIterator vrts_begin, TVertexIterator vrts_end)
341 : {
342 : using namespace std;
343 : TVector boxMin;
344 : TVector boxMax;
345 :
346 0 : TVertexIterator iter = vrts_begin;
347 : // assign initial values
348 : {
349 0 : for(int i = 0; i < numDimensions; i++)
350 0 : boxMin[i] = boxMax[i] = m_aaPos[*iter].coord(i);
351 : }
352 : iter++;
353 :
354 : // loop through the points and calculate the bounding box
355 0 : for(; iter != vrts_end; iter++)
356 : {
357 0 : for(int i = 0; i < numDimensions; ++i)
358 : {
359 0 : boxMin[i] = min(boxMin[i], m_aaPos[*iter].coord(i));
360 0 : boxMax[i] = max(boxMax[i], m_aaPos[*iter].coord(i));
361 : }
362 : }
363 : // calculate extension in each dimension
364 : TVector extension;
365 : {
366 0 : for(int i = 0; i < numDimensions; ++i)
367 0 : extension[i] = boxMax[i] - boxMin[i];
368 : }
369 : // return the index of the largest extension
370 : int bCI = 0;
371 : {
372 0 : for(uint i = 1; i < TVector::Size; ++i)
373 : {
374 0 : if(extension[i] > extension[bCI])
375 0 : bCI = (int)i;
376 : }
377 : }
378 0 : return bCI;
379 : }
380 :
381 : template<class TPositionAttachment, int numDimensions, class TVector>
382 : template <class TVertexIterator>
383 : int
384 0 : KDTreeStatic<TPositionAttachment, numDimensions, TVector>::
385 : get_next_split_dimension(int actSplitDimension, TVertexIterator vrts_begin, TVertexIterator vrts_end)
386 : {
387 0 : switch(m_splitDimension)
388 : {
389 0 : case KDSD_LARGEST:
390 0 : return get_largest_dimension(vrts_begin, vrts_end);
391 : break;
392 0 : case KDSD_CIRCULAR:
393 0 : return (actSplitDimension+1) % numDimensions;
394 : break;
395 :
396 : }
397 : // default: SD_CIRCULAR
398 0 : return (actSplitDimension+1) % numDimensions;
399 : }
400 :
401 : template<class TPositionAttachment, int numDimensions, class TVector>
402 : void
403 : KDTreeStatic<TPositionAttachment, numDimensions, TVector>::
404 : get_leafs_recursive(std::vector<Node*>& vLeafsOut, Node* pNode)
405 : {
406 : // check whether the node is a leaf
407 : if(!(pNode->m_pChild[0] || pNode->m_pChild[1]))
408 : {
409 : // it is a leaf. push it to the list
410 : vLeafsOut.push_back(pNode);
411 : }
412 : else
413 : {
414 : for(int i = 0; i < 2; ++i)
415 : {
416 : if(pNode->m_pChild[i] != NULL)
417 : get_leafs_recursive(vLeafsOut, pNode->m_pChild[i]);
418 : }
419 : }
420 : }
421 :
422 : }// end of namespace
423 :
424 : #endif
|
