1 |
|
|
/* |
2 |
|
|
Teem: Tools to process and visualize scientific data and images . |
3 |
|
|
Copyright (C) 2013, 2012, 2011, 2010, 2009 University of Chicago |
4 |
|
|
Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann |
5 |
|
|
Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah |
6 |
|
|
|
7 |
|
|
This library is free software; you can redistribute it and/or |
8 |
|
|
modify it under the terms of the GNU Lesser General Public License |
9 |
|
|
(LGPL) as published by the Free Software Foundation; either |
10 |
|
|
version 2.1 of the License, or (at your option) any later version. |
11 |
|
|
The terms of redistributing and/or modifying this software also |
12 |
|
|
include exceptions to the LGPL that facilitate static linking. |
13 |
|
|
|
14 |
|
|
This library is distributed in the hope that it will be useful, |
15 |
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of |
16 |
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
17 |
|
|
Lesser General Public License for more details. |
18 |
|
|
|
19 |
|
|
You should have received a copy of the GNU Lesser General Public License |
20 |
|
|
along with this library; if not, write to Free Software Foundation, Inc., |
21 |
|
|
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA |
22 |
|
|
*/ |
23 |
|
|
|
24 |
|
|
#include "pull.h" |
25 |
|
|
#include "privatePull.h" |
26 |
|
|
|
27 |
|
|
/* |
28 |
|
|
** because the pullContext keeps an array of bins (not pointers to them) |
29 |
|
|
** we have Init and Done functions (not New and Nix) |
30 |
|
|
*/ |
31 |
|
|
void |
32 |
|
|
_pullBinInit(pullBin *bin) { |
33 |
|
|
|
34 |
|
|
bin->point = NULL; |
35 |
|
|
bin->pointNum = 0; |
36 |
|
|
bin->pointArr = NULL; |
37 |
|
|
bin->neighBin = NULL; |
38 |
|
|
return; |
39 |
|
|
} |
40 |
|
|
|
41 |
|
|
/* |
42 |
|
|
** bins own the points they contain- so this frees them |
43 |
|
|
*/ |
44 |
|
|
void |
45 |
|
|
_pullBinDone(pullBin *bin) { |
46 |
|
|
unsigned int idx; |
47 |
|
|
|
48 |
|
|
for (idx=0; idx<bin->pointNum; idx++) { |
49 |
|
|
bin->point[idx] = pullPointNix(bin->point[idx]); |
50 |
|
|
} |
51 |
|
|
bin->pointArr = airArrayNuke(bin->pointArr); |
52 |
|
|
bin->neighBin = (pullBin **)airFree(bin->neighBin); |
53 |
|
|
return; |
54 |
|
|
} |
55 |
|
|
|
56 |
|
|
int |
57 |
|
|
_pullBinNeighborSet(pullContext *pctx, pullBin *bin) { |
58 |
|
|
static const char me[]="_pullBinNeighborSet"; |
59 |
|
|
unsigned int neiIdx, neiNum, be[4], binIdx; |
60 |
|
|
unsigned int xi, yi, zi, si, xx, yy, zz, ss, |
61 |
|
|
xmax, ymax, zmax, smax; |
62 |
|
|
pullBin *nei[3*3*3*3]; |
63 |
|
|
int xmin, ymin, zmin, smin; |
64 |
|
|
|
65 |
|
|
binIdx = AIR_UINT(bin - pctx->bin); |
66 |
|
|
/* annoyingly, have to recover the bin coordinates */ |
67 |
|
|
ELL_4V_COPY(be, pctx->binsEdge); |
68 |
|
|
xi = binIdx % be[0]; |
69 |
|
|
binIdx = (binIdx - xi)/be[0]; |
70 |
|
|
yi = binIdx % be[1]; |
71 |
|
|
binIdx = (binIdx - yi)/be[1]; |
72 |
|
|
zi = binIdx % be[2]; |
73 |
|
|
si = (binIdx - zi)/be[2]; |
74 |
|
|
neiNum = 0; |
75 |
|
|
bin->neighBin = (pullBin **)airFree(bin->neighBin); |
76 |
|
|
smin = AIR_MAX(0, (int)si-1); |
77 |
|
|
smax = AIR_MIN(si+1, be[3]-1); |
78 |
|
|
zmin = AIR_MAX(0, (int)zi-1); |
79 |
|
|
zmax = AIR_MIN(zi+1, be[2]-1); |
80 |
|
|
ymin = AIR_MAX(0, (int)yi-1); |
81 |
|
|
ymax = AIR_MIN(yi+1, be[1]-1); |
82 |
|
|
xmin = AIR_MAX(0, (int)xi-1); |
83 |
|
|
xmax = AIR_MIN(xi+1, be[0]-1); |
84 |
|
|
for (ss=smin; ss<=smax; ss++) { |
85 |
|
|
for (zz=zmin; zz<=zmax; zz++) { |
86 |
|
|
for (yy=ymin; yy<=ymax; yy++) { |
87 |
|
|
for (xx=xmin; xx<=xmax; xx++) { |
88 |
|
|
nei[neiNum++] = pctx->bin + xx + be[0]*(yy + be[1]*(zz + be[2]*ss)); |
89 |
|
|
} |
90 |
|
|
} |
91 |
|
|
} |
92 |
|
|
} |
93 |
|
|
if (!( bin->neighBin = AIR_CALLOC(1+neiNum, pullBin*) )) { |
94 |
|
|
biffAddf(PULL, "%s: couldn't calloc array of %u neighbor pointers", me, 1+neiNum); |
95 |
|
|
return 1; |
96 |
|
|
} |
97 |
|
|
for (neiIdx=0; neiIdx<neiNum; neiIdx++) { |
98 |
|
|
bin->neighBin[neiIdx] = nei[neiIdx]; |
99 |
|
|
} |
100 |
|
|
/* NULL-terminate the bin array */ |
101 |
|
|
bin->neighBin[neiIdx] = NULL; |
102 |
|
|
return 0; |
103 |
|
|
} |
104 |
|
|
|
105 |
|
|
/* |
106 |
|
|
** bins on boundary now extend to infinity; so the only time this |
107 |
|
|
** returns NULL (indicating error) is for non-existent positions |
108 |
|
|
*/ |
109 |
|
|
pullBin * |
110 |
|
|
_pullBinLocate(pullContext *pctx, double *posWorld) { |
111 |
|
|
static const char me[]="_pullBinLocate"; |
112 |
|
|
unsigned int axi, eidx[4], binIdx; |
113 |
|
|
|
114 |
|
|
if (!ELL_4V_EXISTS(posWorld)) { |
115 |
|
|
biffAddf(PULL, "%s: non-existent position (%g,%g,%g,%g)", me, |
116 |
|
|
posWorld[0], posWorld[1], posWorld[2], posWorld[3]); |
117 |
|
|
return NULL; |
118 |
|
|
} |
119 |
|
|
|
120 |
|
|
if (pctx->flag.binSingle) { |
121 |
|
|
binIdx = 0; |
122 |
|
|
} else { |
123 |
|
|
for (axi=0; axi<4; axi++) { |
124 |
|
|
eidx[axi] = airIndexClamp(pctx->bboxMin[axi], |
125 |
|
|
posWorld[axi], |
126 |
|
|
pctx->bboxMax[axi], |
127 |
|
|
pctx->binsEdge[axi]); |
128 |
|
|
} |
129 |
|
|
binIdx = (eidx[0] |
130 |
|
|
+ pctx->binsEdge[0]*( |
131 |
|
|
eidx[1] + pctx->binsEdge[1]*( |
132 |
|
|
eidx[2] + pctx->binsEdge[2] * eidx[3]))); |
133 |
|
|
} |
134 |
|
|
|
135 |
|
|
return pctx->bin + binIdx; |
136 |
|
|
} |
137 |
|
|
|
138 |
|
|
/* |
139 |
|
|
** this makes the bin the owner of the point |
140 |
|
|
*/ |
141 |
|
|
int |
142 |
|
|
_pullBinPointAdd(pullContext *pctx, pullBin *bin, pullPoint *point) { |
143 |
|
|
static const char me[]="_pullBinPointAdd"; |
144 |
|
|
int pntI; |
145 |
|
|
pullPtrPtrUnion pppu; |
146 |
|
|
|
147 |
|
|
AIR_UNUSED(pctx); |
148 |
|
|
if (!(bin->pointArr)) { |
149 |
|
|
pppu.points = &(bin->point); |
150 |
|
|
bin->pointArr = airArrayNew(pppu.v, &(bin->pointNum), |
151 |
|
|
sizeof(pullPoint *), _PULL_BIN_INCR); |
152 |
|
|
if (!( bin->pointArr )) { |
153 |
|
|
biffAddf(PULL, "%s: couldn't create point array", me); |
154 |
|
|
return 1; |
155 |
|
|
} |
156 |
|
|
} |
157 |
|
|
if (!( bin->neighBin )) { |
158 |
|
|
/* set up neighbor bin vector if not done so already */ |
159 |
|
|
if (_pullBinNeighborSet(pctx, bin)) { |
160 |
|
|
biffAddf(PULL, "%s: couldn't initialize neighbor bins", me); |
161 |
|
|
return 1; |
162 |
|
|
} |
163 |
|
|
} |
164 |
|
|
pntI = airArrayLenIncr(bin->pointArr, 1); |
165 |
|
|
bin->point[pntI] = point; |
166 |
|
|
return 0; |
167 |
|
|
} |
168 |
|
|
|
169 |
|
|
/* |
170 |
|
|
** the bin loses track of the point, caller responsible for ownership, |
171 |
|
|
** even though caller never identifies it by pointer, which is weird |
172 |
|
|
*/ |
173 |
|
|
void |
174 |
|
|
_pullBinPointRemove(pullContext *pctx, pullBin *bin, int loseIdx) { |
175 |
|
|
|
176 |
|
|
AIR_UNUSED(pctx); |
177 |
|
|
bin->point[loseIdx] = bin->point[bin->pointNum-1]; |
178 |
|
|
airArrayLenIncr(bin->pointArr, -1); |
179 |
|
|
return; |
180 |
|
|
} |
181 |
|
|
|
182 |
|
|
/* |
183 |
|
|
** adds point to context |
184 |
|
|
*/ |
185 |
|
|
int |
186 |
|
|
pullBinsPointAdd(pullContext *pctx, pullPoint *point, pullBin **binP) { |
187 |
|
|
static const char me[]="pullBinsPointAdd"; |
188 |
|
|
pullBin *bin; |
189 |
|
|
|
190 |
|
|
if (binP) { |
191 |
|
|
*binP = NULL; |
192 |
|
|
} |
193 |
|
|
if (!( bin = _pullBinLocate(pctx, point->pos) )) { |
194 |
|
|
biffAddf(PULL, "%s: can't locate point %p %u", |
195 |
|
|
me, AIR_CAST(void*, point), point->idtag); |
196 |
|
|
return 1; |
197 |
|
|
} |
198 |
|
|
if (binP) { |
199 |
|
|
*binP = bin; |
200 |
|
|
} |
201 |
|
|
if (_pullBinPointAdd(pctx, bin, point)) { |
202 |
|
|
biffAddf(PULL, "%s: trouble adding point %p %u", |
203 |
|
|
me, AIR_CAST(void*, point), point->idtag); |
204 |
|
|
return 1; |
205 |
|
|
} |
206 |
|
|
return 0; |
207 |
|
|
} |
208 |
|
|
|
209 |
|
|
int |
210 |
|
|
pullBinsPointMaybeAdd(pullContext *pctx, pullPoint *point, |
211 |
|
|
/* output */ |
212 |
|
|
pullBin **binP, int *added) { |
213 |
|
|
static const char me[]="pullBinsPointMaybeAdd"; |
214 |
|
|
pullBin *bin; |
215 |
|
|
unsigned int idx; |
216 |
|
|
int okay; |
217 |
|
|
|
218 |
|
|
if (binP) { |
219 |
|
|
*binP = NULL; |
220 |
|
|
} |
221 |
|
|
if (!(pctx && point && added)) { |
222 |
|
|
biffAddf(PULL, "%s: got NULL pointer", me); |
223 |
|
|
return 1; |
224 |
|
|
} |
225 |
|
|
if (!( bin = _pullBinLocate(pctx, point->pos) )) { |
226 |
|
|
biffAddf(PULL, "%s: can't locate point %p %u", |
227 |
|
|
me, AIR_CAST(void*, point), point->idtag); |
228 |
|
|
return 1; |
229 |
|
|
} |
230 |
|
|
if (binP) { |
231 |
|
|
*binP = bin; |
232 |
|
|
} |
233 |
|
|
if (pctx->flag.restrictiveAddToBins) { |
234 |
|
|
okay = AIR_TRUE; |
235 |
|
|
for (idx=0; idx<bin->pointNum; idx++) { |
236 |
|
|
double diff[4], len; |
237 |
|
|
ELL_4V_SUB(diff, point->pos, bin->point[idx]->pos); |
238 |
|
|
ELL_3V_SCALE(diff, 1/pctx->sysParm.radiusSpace, diff); |
239 |
|
|
diff[3] /= pctx->sysParm.radiusScale; |
240 |
|
|
len = ELL_4V_LEN(diff); |
241 |
|
|
if (len < _PULL_BINNING_MAYBE_ADD_THRESH) { |
242 |
|
|
okay = AIR_FALSE; |
243 |
|
|
break; |
244 |
|
|
} |
245 |
|
|
} |
246 |
|
|
if (okay) { |
247 |
|
|
if (_pullBinPointAdd(pctx, bin, point)) { |
248 |
|
|
biffAddf(PULL, "%s: trouble adding point %p %u", |
249 |
|
|
me, AIR_CAST(void*, point), point->idtag); |
250 |
|
|
return 1; |
251 |
|
|
} |
252 |
|
|
*added = AIR_TRUE; |
253 |
|
|
} else { |
254 |
|
|
*added = AIR_FALSE; |
255 |
|
|
} |
256 |
|
|
} else { |
257 |
|
|
if (_pullBinPointAdd(pctx, bin, point)) { |
258 |
|
|
biffAddf(PULL, "%s: trouble adding point %p %u", |
259 |
|
|
me, AIR_CAST(void*, point), point->idtag); |
260 |
|
|
return 1; |
261 |
|
|
} |
262 |
|
|
*added = AIR_TRUE; |
263 |
|
|
} |
264 |
|
|
return 0; |
265 |
|
|
} |
266 |
|
|
|
267 |
|
|
int |
268 |
|
|
_pullBinSetup(pullContext *pctx) { |
269 |
|
|
static const char me[]="_pullBinSetup"; |
270 |
|
|
unsigned ii; |
271 |
|
|
double volEdge[4], width; |
272 |
|
|
|
273 |
|
|
/* the maximum distance of interaction is when one particle is sitting |
274 |
|
|
on the edge of another particle's sphere of influence, *NOT*, when |
275 |
|
|
the spheres of influence of two particles are tangent: particles |
276 |
|
|
interact with potential fields of other particles, but there is |
277 |
|
|
no interaction between potential fields. */ |
278 |
|
|
width = (pctx->sysParm.radiusSpace ? pctx->sysParm.radiusSpace : 0.1); |
279 |
|
|
pctx->maxDistSpace = pctx->sysParm.binWidthSpace*width; |
280 |
|
|
width = (pctx->sysParm.radiusScale ? pctx->sysParm.radiusScale : 0.1); |
281 |
|
|
pctx->maxDistScale = 1*width; |
282 |
|
|
|
283 |
|
|
if (pctx->verbose) { |
284 |
|
|
printf("%s: radiusSpace = %g -(%g)-> maxDistSpace = %g\n", me, |
285 |
|
|
pctx->sysParm.radiusSpace, pctx->sysParm.binWidthSpace, |
286 |
|
|
pctx->maxDistSpace); |
287 |
|
|
printf("%s: radiusScale = %g ----> maxDistScale = %g\n", me, |
288 |
|
|
pctx->sysParm.radiusScale, pctx->maxDistScale); |
289 |
|
|
} |
290 |
|
|
|
291 |
|
|
if (pctx->flag.binSingle) { |
292 |
|
|
pctx->binsEdge[0] = 1; |
293 |
|
|
pctx->binsEdge[1] = 1; |
294 |
|
|
pctx->binsEdge[2] = 1; |
295 |
|
|
pctx->binsEdge[3] = 1; |
296 |
|
|
pctx->binNum = 1; |
297 |
|
|
} else { |
298 |
|
|
volEdge[0] = pctx->bboxMax[0] - pctx->bboxMin[0]; |
299 |
|
|
volEdge[1] = pctx->bboxMax[1] - pctx->bboxMin[1]; |
300 |
|
|
volEdge[2] = pctx->bboxMax[2] - pctx->bboxMin[2]; |
301 |
|
|
volEdge[3] = pctx->bboxMax[3] - pctx->bboxMin[3]; |
302 |
|
|
if (pctx->verbose) { |
303 |
|
|
printf("%s: volEdge = %g %g %g %g\n", me, |
304 |
|
|
volEdge[0], volEdge[1], volEdge[2], volEdge[3]); |
305 |
|
|
} |
306 |
|
|
pctx->binsEdge[0] = AIR_CAST(unsigned int, |
307 |
|
|
floor(volEdge[0]/pctx->maxDistSpace)); |
308 |
|
|
pctx->binsEdge[0] = pctx->binsEdge[0] ? pctx->binsEdge[0] : 1; |
309 |
|
|
pctx->binsEdge[1] = AIR_CAST(unsigned int, |
310 |
|
|
floor(volEdge[1]/pctx->maxDistSpace)); |
311 |
|
|
pctx->binsEdge[1] = pctx->binsEdge[1] ? pctx->binsEdge[1] : 1; |
312 |
|
|
pctx->binsEdge[2] = AIR_CAST(unsigned int, |
313 |
|
|
floor(volEdge[2]/pctx->maxDistSpace)); |
314 |
|
|
pctx->binsEdge[2] = pctx->binsEdge[2] ? pctx->binsEdge[2] : 1; |
315 |
|
|
pctx->binsEdge[3] = AIR_CAST(unsigned int, |
316 |
|
|
floor(volEdge[3]/pctx->maxDistScale)); |
317 |
|
|
pctx->binsEdge[3] = pctx->binsEdge[3] ? pctx->binsEdge[3] : 1; |
318 |
|
|
/* hack to observe things at bin boundaries |
319 |
|
|
ELL_3V_SET(pctx->binsEdge, 3, 3, 3); |
320 |
|
|
*/ |
321 |
|
|
if (pctx->verbose) { |
322 |
|
|
printf("%s: binsEdge=(%u,%u,%u,%u)\n", me, |
323 |
|
|
pctx->binsEdge[0], pctx->binsEdge[1], |
324 |
|
|
pctx->binsEdge[2], pctx->binsEdge[3]); |
325 |
|
|
} |
326 |
|
|
pctx->binNum = (pctx->binsEdge[0]*pctx->binsEdge[1] |
327 |
|
|
*pctx->binsEdge[2]*pctx->binsEdge[3]); |
328 |
|
|
} |
329 |
|
|
if (pctx->binNum > PULL_BIN_MAXNUM) { |
330 |
|
|
biffAddf(PULL, "%s: sorry, #bins %u > PULL_BIN_MAXNUM %u. Try " |
331 |
|
|
"increasing pctx->sysParm.binWidthSpace (%g)", me, |
332 |
|
|
pctx->binNum, PULL_BIN_MAXNUM, pctx->sysParm.binWidthSpace); |
333 |
|
|
return 1; |
334 |
|
|
} |
335 |
|
|
if (pctx->verbose) { |
336 |
|
|
printf("%s: trying to allocate %u bins . . . \n", me, pctx->binNum); |
337 |
|
|
} |
338 |
|
|
pctx->bin = AIR_CALLOC(pctx->binNum, pullBin); |
339 |
|
|
if (!( pctx->bin )) { |
340 |
|
|
biffAddf(PULL, "%s: couln't allocate %u bins", me, pctx->binNum); |
341 |
|
|
return 1; |
342 |
|
|
} |
343 |
|
|
if (pctx->verbose) { |
344 |
|
|
printf("%s: done allocating. Initializing . . . \n", me); |
345 |
|
|
} |
346 |
|
|
for (ii=0; ii<pctx->binNum; ii++) { |
347 |
|
|
_pullBinInit(pctx->bin + ii); |
348 |
|
|
} |
349 |
|
|
if (pctx->verbose) { |
350 |
|
|
printf("%s: done initializing.\n", me); |
351 |
|
|
} |
352 |
|
|
if (pctx->flag.binSingle) { |
353 |
|
|
if (!( pctx->bin[0].neighBin = AIR_CALLOC(2, pullBin*) )) { |
354 |
|
|
biffAddf(PULL, "%s: trouble allocating for single bin?", me); |
355 |
|
|
return 1; |
356 |
|
|
} |
357 |
|
|
pctx->bin[0].neighBin[0] = pctx->bin + 0; |
358 |
|
|
pctx->bin[0].neighBin[1] = NULL; |
359 |
|
|
} |
360 |
|
|
return 0; |
361 |
|
|
} |
362 |
|
|
|
363 |
|
|
void |
364 |
|
|
_pullBinFinish(pullContext *pctx) { |
365 |
|
|
unsigned int ii; |
366 |
|
|
|
367 |
|
|
for (ii=0; ii<pctx->binNum; ii++) { |
368 |
|
|
_pullBinDone(pctx->bin + ii); |
369 |
|
|
} |
370 |
|
|
pctx->bin = (pullBin *)airFree(pctx->bin); |
371 |
|
|
ELL_4V_SET(pctx->binsEdge, 0, 0, 0, 0); |
372 |
|
|
pctx->binNum = 0; |
373 |
|
|
} |
374 |
|
|
|
375 |
|
|
/* |
376 |
|
|
** sets pctx->stuckNum |
377 |
|
|
** resets all task[]->stuckNum |
378 |
|
|
** reallocates pctx->tmpPointPerm and pctx->tmpPointPtr |
379 |
|
|
** the point of this is to do rebinning |
380 |
|
|
** |
381 |
|
|
** This function is only called by the master thread, this |
382 |
|
|
** does *not* have to be thread-safe in any way |
383 |
|
|
*/ |
384 |
|
|
int |
385 |
|
|
_pullIterFinishDescent(pullContext *pctx) { |
386 |
|
|
static const char me[]="_pullIterFinishDescent"; |
387 |
|
|
unsigned int oldBinIdx, pointIdx, taskIdx, runIdx, pointNum; |
388 |
|
|
pullBin *oldBin; |
389 |
|
|
pullPoint *point; |
390 |
|
|
int added; |
391 |
|
|
|
392 |
|
|
_pullNixTheNixed(pctx); |
393 |
|
|
|
394 |
|
|
pctx->stuckNum = 0; |
395 |
|
|
for (taskIdx=0; taskIdx<pctx->threadNum; taskIdx++) { |
396 |
|
|
pctx->stuckNum += pctx->task[taskIdx]->stuckNum; |
397 |
|
|
pctx->task[taskIdx]->stuckNum = 0; |
398 |
|
|
} |
399 |
|
|
|
400 |
|
|
/* even w/ a single bin, we may still have to permute the points */ |
401 |
|
|
pointNum = pullPointNumber(pctx); |
402 |
|
|
if (pointNum != pctx->tmpPointNum) { |
403 |
|
|
if (pctx->verbose) { |
404 |
|
|
printf("!%s: changing total point # %u --> %u\n", me, |
405 |
|
|
pctx->tmpPointNum, pointNum); |
406 |
|
|
} |
407 |
|
|
airFree(pctx->tmpPointPerm); |
408 |
|
|
airFree(pctx->tmpPointPtr); |
409 |
|
|
pctx->tmpPointPtr = AIR_CAST(pullPoint **, |
410 |
|
|
calloc(pointNum, sizeof(pullPoint*))); |
411 |
|
|
pctx->tmpPointPerm = AIR_CAST(unsigned int *, |
412 |
|
|
calloc(pointNum, sizeof(unsigned int))); |
413 |
|
|
if (!( pctx->tmpPointPtr && pctx->tmpPointPerm )) { |
414 |
|
|
biffAddf(PULL, "%s: couldn't allocate tmp buffers %p %p", me, |
415 |
|
|
AIR_VOIDP(pctx->tmpPointPtr), AIR_VOIDP(pctx->tmpPointPerm)); |
416 |
|
|
return 1; |
417 |
|
|
} |
418 |
|
|
pctx->tmpPointNum = pointNum; |
419 |
|
|
} |
420 |
|
|
runIdx = 0; |
421 |
|
|
for (oldBinIdx=0; oldBinIdx<pctx->binNum; oldBinIdx++) { |
422 |
|
|
oldBin = pctx->bin + oldBinIdx; |
423 |
|
|
while (oldBin->pointNum) { |
424 |
|
|
/* tricky: we can't traverse bin->point[], because of how it is |
425 |
|
|
re-ordered on point removal, so we always grab point[0] */ |
426 |
|
|
pctx->tmpPointPtr[runIdx++] = oldBin->point[0]; |
427 |
|
|
_pullBinPointRemove(pctx, oldBin, 0); |
428 |
|
|
} |
429 |
|
|
} |
430 |
|
|
airShuffle_r(pctx->task[0]->rng, |
431 |
|
|
pctx->tmpPointPerm, pointNum, pctx->flag.permuteOnRebin); |
432 |
|
|
if (pctx->flag.permuteOnRebin && pctx->verbose) { |
433 |
|
|
printf("%s: permuting %u points\n", me, pointNum); |
434 |
|
|
} |
435 |
|
|
for (pointIdx=0; pointIdx<pointNum; pointIdx++) { |
436 |
|
|
point = pctx->tmpPointPtr[pctx->tmpPointPerm[pointIdx]]; |
437 |
|
|
/* |
438 |
|
|
if (1 == pctx->iter && 102 == point->idtag) { |
439 |
|
|
_pullDebugSanity(pctx->task[0], point); |
440 |
|
|
} |
441 |
|
|
*/ |
442 |
|
|
/* Sun Jul 14 01:30:41 CDT 2013: the previous code was basically |
443 |
|
|
just pullBinsPointAdd(). But when working with codim-3 |
444 |
|
|
constraints (like finding spatial maxima with maximal strength |
445 |
|
|
along scale), its easy for many points to start piling on top |
446 |
|
|
of each other; that is the problem that |
447 |
|
|
pullFlagRestrictiveAddToBins was designed to solve. */ |
448 |
|
|
if (pctx->constraint && 0 == pctx->constraintDim) { |
449 |
|
|
if (pullBinsPointMaybeAdd(pctx, point, NULL, &added)) { |
450 |
|
|
biffAddf(PULL, "%s: trouble binning? point %u", me, point->idtag); |
451 |
|
|
return 1; |
452 |
|
|
} |
453 |
|
|
if (!added) { |
454 |
|
|
/* the point wasn't owned by any bin, and now it turns out |
455 |
|
|
no bin wanted to own it, so we have to remove it */ |
456 |
|
|
/* in the case of point maxima searching; this mainly happened |
457 |
|
|
because two points at the same spatial positition slowly |
458 |
|
|
descended along scale towards each other at the scale of |
459 |
|
|
maximal strength */ |
460 |
|
|
point = pullPointNix(point); |
461 |
|
|
} |
462 |
|
|
} else { |
463 |
|
|
if (pullBinsPointAdd(pctx, point, NULL)) { |
464 |
|
|
biffAddf(PULL, "%s: trouble binning point %u", me, point->idtag); |
465 |
|
|
return 1; |
466 |
|
|
} |
467 |
|
|
point = NULL; |
468 |
|
|
} |
469 |
|
|
pctx->tmpPointPtr[pctx->tmpPointPerm[pointIdx]] = NULL; |
470 |
|
|
} |
471 |
|
|
|
472 |
|
|
return 0; |
473 |
|
|
} |
474 |
|
|
|