Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] Annotation of /branches/vis15/src/lib/include/diderot/kdtree-inst.hxx
ViewVC logotype

Annotation of /branches/vis15/src/lib/include/diderot/kdtree-inst.hxx

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4372 - (view) (download) (as text)

1 : jhr 4365 /*! \file kdtree-inst.hxx
2 :     *
3 :     * \author John Reppy
4 :     */
5 :    
6 :     /*
7 :     * This code is part of the Diderot Project (http://diderot-language.cs.uchicago.edu)
8 :     *
9 :     * COPYRIGHT (c) 2016 The University of Chicago
10 :     * All rights reserved.
11 :     */
12 :    
13 :     #ifndef _DIDEROT_KDTREE_INST_HXX_
14 :     #define _DIDEROT_KDTREE_INST_HXX_
15 :    
16 :     #ifndef _DIDEROT_KDTREE_HXX_
17 :     # error kdtree-inst.hxx should not be directly included
18 :     #endif
19 :    
20 :     #include <stack>
21 :    
22 :     namespace diderot {
23 :    
24 :     namespace __details {
25 :    
26 :     // generic test for two points being within a given radius^2
27 :     template <const uint32_t D, typename REAL>
28 : jhr 4372 bool within_sphere (const REAL pos1[D], const REAL pos2[D], REAL radius2)
29 : jhr 4365 {
30 : jhr 4372 float sum = REAL(0);
31 : jhr 4365 for (uint32_t i = 0; i < D; i++) {
32 :     REAL d = pos1[i] - pos2[i];
33 : jhr 4372 sum += (d * d);
34 : jhr 4365 }
35 : jhr 4372 return sum <= radius2;
36 : jhr 4365 }
37 :    
38 :     // generic test for two points being within a specified retangular radius?
39 :     template <const uint32_t D, typename REAL>
40 : jhr 4372 bool within_box (const REAL pos1[D], const REAL pos2[D], REAL radius)
41 : jhr 4365 {
42 :     for (uint32_t i = 0; i < D; i++) {
43 :     if (std::abs(pos1 - pos2) > radius) {
44 :     return false;
45 :     }
46 :     }
47 :     return true;
48 :     }
49 :    
50 :     } // namespace __details
51 :    
52 :     template <const uint32_t D, typename REAL, typename S>
53 : jhr 4372 kdtree<D,REAL,S>::kdtree (uint32_t nStrands)
54 :     : _nStrands(nStrands), _poolSz(nStrands),
55 : jhr 4365 _parts(nullptr), _pool(nullptr), _strands(nullptr)
56 :     {
57 : jhr 4372 // FIXME: for "grid" programs, we don't need the extra space!
58 :     this->_partsSz = nStrands + (nStrands >> 2);
59 :     delete[] this->_parts;
60 :     this->_parts = new uint32_t[nStrands + (nStrands >> 2)];
61 :     for (int i = 0; i < nStrands; i++) {
62 :     this->_parts[i] = i;
63 :     }
64 :    
65 :     uint32_t reqNumNodes =
66 :     2 * ((nStrands + kdtree<D,REAL,S>::STRANDS_PER_LEAF - 1) / STRANDS_PER_LEAF) - 1;
67 :     if (reqNumNodes < this->_poolSz) {
68 :     delete[] this->_pool;
69 :     this->_pool = new node[reqNumNodes];
70 :     this->_poolSz = reqNumNodes;
71 :     }
72 :    
73 : jhr 4365 }
74 :    
75 :     template <const uint32_t D, typename REAL, typename S>
76 :     kdtree<D,REAL,S>::~kdtree ()
77 :     {
78 :     delete[] _parts;
79 :     delete[] _pool;
80 :     }
81 :    
82 :     template <const uint32_t D, typename REAL, typename S>
83 :     void kdtree<D,REAL,S>::swap_parts (uint32_t i, uint32_t j)
84 :     {
85 :     uint32_t tmp = this->_parts[i];
86 :     this->_parts[i] = this->_parts[j];
87 :     this->_parts[j] = tmp;
88 :     }
89 :    
90 :     // partition _parts[lo..hi] such that _parts[lo..pivotIx] <= X < _parts[pivotIx+1..hi], where
91 :     // X is the initial value of _strands[_parts[pivotIx]].pos()[axis].
92 :     //
93 :     template <const uint32_t D, typename REAL, typename S>
94 : jhr 4367 uint32_t kdtree<D,REAL,S>::partition (uint32_t axis, uint32_t lo, uint32_t hi, uint32_t pivotIx)
95 : jhr 4365 {
96 : jhr 4372 REAL X = this->strand(pivotIx)->pos()[axis];
97 : jhr 4365
98 :     // move pivot element to end
99 :     this->swap_parts(pivotIx, hi);
100 :    
101 :     uint32_t ix = lo;
102 :     for (uint32_t jx = lo; jx < hi-1; jx++) {
103 : jhr 4372 if (this->strand(jx)->pos()[axis] <= X) {
104 : jhr 4365 this->swap_parts (ix, jx);
105 :     ix++;
106 :     }
107 :     }
108 :    
109 :     this->swap_parts (ix, pivotIx);
110 : jhr 4367 return ix;
111 : jhr 4365 }
112 : jhr 4372
113 : jhr 4365 // partition _parts[lo..hi] into _parts[lo..m] and _parts[m+1..hi] such that the strand
114 :     // with id _parts[m] has the median position on the specified axis.
115 : jhr 4367 // We use the "Quick Select" method (https://en.wikipedia.org/wiki/Quickselect)
116 :     //
117 : jhr 4365 template <const uint32_t D, typename REAL, typename S>
118 :     uint32_t kdtree<D,REAL,S>::median (uint32_t axis, uint32_t lo, uint32_t hi)
119 :     {
120 : jhr 4367 assert (hi - lo >= STRANDS_PER_LEAF);
121 :     // the mid-point of the interval [lo..hi]
122 :     uint32_t mid = (hi + lo) >> 1;
123 :    
124 :     // partition until we are within +/- 25% of the mid-point
125 :     uint32_t tol = ((mid - lo) >> 2);
126 :    
127 :     while (true) {
128 :     uint32_t pivotIx = this->partition (axis, lo, hi, mid);
129 : jhr 4372 if (std::abs(static_cast<int>(pivotIx) - static_cast<int>(mid)) <= tol) {
130 : jhr 4367 return pivotIx;
131 :     }
132 :     else if (pivotIx < mid) {
133 : jhr 4372 lo = pivotIx + 1;
134 : jhr 4367 }
135 :     else {
136 : jhr 4372 hi = pivotIx - 1;
137 : jhr 4367 }
138 :     }
139 :    
140 : jhr 4365 }
141 :    
142 :     template <const uint32_t D, typename REAL, typename S>
143 : jhr 4367 uint32_t kdtree<D,REAL,S>::builder (uint32_t axis, uint32_t lo, uint32_t hi)
144 :     {
145 : jhr 4372 assert (lo <= hi);
146 : jhr 4367
147 :     // allcate the node
148 :     uint32_t nd = this->_nextNode++;
149 :    
150 :     uint32_t n = hi - lo + 1;
151 :     if (n <= kdtree<D,REAL,S>::STRANDS_PER_LEAF){
152 :     // allocate a leaf
153 :     this->_pool[nd]._lc = 0;
154 :     this->_pool[nd]._u._leaf._first = lo;
155 :     this->_pool[nd]._u._leaf._last = hi;
156 :     }
157 :     else {
158 :     uint32_t mid = this->median (axis, lo, hi);
159 :     // INV: strands indexed by _parts[lo..mid-1] are <= strand[_parts[mid]]
160 :     // and strand[_parts[mid]] < strands indexed by _parts[mid+1..hi]
161 :     this->_pool[nd]._u._nd._id = this->_parts[mid];
162 :     this->_pool[nd]._u._nd._axis = axis;
163 :     axis = (axis + 1) % D;
164 :     this->_pool[nd]._lc = builder (axis, lo, mid-1);
165 :     this->_pool[nd]._rc = builder (axis, mid+1, hi);
166 :     }
167 :    
168 :     return nd;
169 :     }
170 :    
171 :     template <const uint32_t D, typename REAL, typename S>
172 : jhr 4365 void kdtree<D,REAL,S>::rebuild (uint32_t nStrands, const S *strands)
173 :     {
174 :     this->_strands = strands;
175 : jhr 4367
176 :     if (this->_nStrands > nStrands) {
177 :     // # of strands has shrunk from last call to rebuild
178 :     for (uint32_t i = 0; i < nStrands; i++) {
179 :     this->_parts[i] = i;
180 :     }
181 :     }
182 :     else if (this->_nStrands < nStrands) {
183 :     // # of strands has grown from last call to rebuild
184 :     if (this->_partsSz < nStrands) {
185 :     // need to reallocate the _parts array
186 :     this->_partsSz = nStrands + (nStrands >> 2);
187 :     delete[] this->_parts;
188 :     this->_parts = new uint32_t[nStrands + (nStrands >> 2)];
189 :     }
190 :     for (int i = 0; i < nStrands; i++) {
191 :     this->_parts[i] = i;
192 :     }
193 :    
194 :     // a conservative bound on the number of nodes is 2*(ceil(nStrands / STRANDS_PER_LEAF)) - 1
195 :     uint32_t reqNumNodes =
196 : jhr 4372 2 * ((nStrands + kdtree<D,REAL,S>::STRANDS_PER_LEAF - 1) / STRANDS_PER_LEAF) - 1;
197 : jhr 4367 if (reqNumNodes < this->_poolSz) {
198 :     delete[] this->_pool;
199 :     this->_pool = new node[reqNumNodes];
200 :     this->_poolSz = reqNumNodes;
201 :     }
202 :     }
203 :    
204 :     this->_nextNode = 0;
205 :     uint32_t root = this->builder (0, 0, nStrands-1);
206 :     assert (root == 0);
207 :    
208 : jhr 4365 }
209 :    
210 :     // FIXME: need to filter out the strand that is doing the query!!
211 :     template <const uint32_t D, typename REAL, typename S>
212 :     dynseq<uint32_t> kdtree<D,REAL,S>::sphere_query (const REAL pos[D], REAL radius)
213 :     {
214 :     dynseq<uint32_t> result;
215 :    
216 : jhr 4372 // return empty sequence on
217 : jhr 4365 if (radius < 0.0) {
218 :     return result;
219 :     }
220 :    
221 :     // stack of nodes for which we must still visit
222 :     std::stack<const node *> stk;
223 :    
224 :     REAL radius2 = radius * radius;
225 :     const node *nd = this->root();
226 :     do {
227 :     if (nd->isLeaf()) {
228 :     // check strands in leaf to see if they are within the sphere
229 :     for (uint32_t i = nd->_u._leaf._first; i <= nd->_u._leaf._last; i++) {
230 :     uint32_t id = this->_parts[i];
231 : jhr 4372 if (__details::within_sphere<D,REAL>(this->_strands[id].pos(), pos, radius2)) {
232 : jhr 4365 // add the strand to the result list
233 :     result.append (id);
234 :     }
235 :     }
236 :     if (! stk.empty()) {
237 :     // continue searching
238 :     nd = stk.top ();
239 :     stk.pop ();
240 :     }
241 :     else {
242 :     nd = nullptr;
243 :     }
244 :     }
245 :     else {
246 :     uint32_t axis = nd->axis();
247 :     REAL sPos = this->strand(nd)->pos()[axis];
248 :     if (pos[axis] < sPos - radius) {
249 :     nd = this->left(nd);
250 :     }
251 :     else if (sPos + radius < pos[axis]) {
252 :     nd = this->right(nd);
253 :     }
254 :     else {
255 : jhr 4372 if (__details::within_sphere<D,REAL>(this->strand(nd)->pos(), pos, radius2)) {
256 : jhr 4365 // add the strand to the result list
257 :     result.append (nd->_u._nd._id);
258 : jhr 4372 std::cout << "result.append(" << nd->_u._nd._id << "); length = " << result.length() << std::endl;
259 : jhr 4365 }
260 :     stk.push(this->right(nd)); // visit right child later
261 :     nd = this->left(nd);
262 :     }
263 :     }
264 :     } while (nd != nullptr);
265 :    
266 :     return result;
267 :     }
268 :    
269 :     } // namespace diderot
270 :    
271 :     #endif // !_DIDEROT_KDTREE_INST_HXX_

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0