Bug Summary

File:src/pull/pointPull.c
Location:line 1320, column 5
Description:Value stored to 'point' is never read

Annotated Source Code

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
25#include "pull.h"
26#include "privatePull.h"
27
28#define DEBUG1 1
29
30/*
31** HEY: this has to be threadsafe, at least threadsafe when there
32** are no errors, because this can now be called from multiple
33** tasks during population control
34*/
35pullPoint *
36pullPointNew(pullContext *pctx) {
37 static const char me[]="pullPointNew";
38 pullPoint *pnt;
39 unsigned int ii;
40 size_t pntSize;
41 pullPtrPtrUnion pppu;
42
43 if (!pctx) {
44 biffAddf(PULLpullBiffKey, "%s: got NULL pointer", me);
45 return NULL((void*)0);
46 }
47 if (!pctx->infoTotalLen) {
48 biffAddf(PULLpullBiffKey, "%s: can't allocate points w/out infoTotalLen set\n", me);
49 return NULL((void*)0);
50 }
51 /* Allocate the pullPoint so that it has pctx->infoTotalLen doubles.
52 The pullPoint declaration has info[1], hence the "- 1" below */
53 pntSize = sizeof(pullPoint) + sizeof(double)*(pctx->infoTotalLen - 1);
54 pnt = AIR_CAST(pullPoint *, calloc(1, pntSize))((pullPoint *)(calloc(1, pntSize)));
55 if (!pnt) {
56 biffAddf(PULLpullBiffKey, "%s: couldn't allocate point (info len %u)\n", me,
57 pctx->infoTotalLen - 1);
58 return NULL((void*)0);
59 }
60
61 pnt->idtag = pctx->idtagNext++;
62 pnt->idCC = 0;
63 pnt->neighPoint = NULL((void*)0);
64 pnt->neighPointNum = 0;
65 pppu.points = &(pnt->neighPoint);
66 pnt->neighPointArr = airArrayNew(pppu.v, &(pnt->neighPointNum),
67 sizeof(pullPoint *),
68 PULL_POINT_NEIGH_INCR16);
69 pnt->neighPointArr->noReallocWhenSmaller = AIR_TRUE1;
70 pnt->neighDistMean = 0;
71 ELL_10V_ZERO_SET(pnt->neighCovar)((pnt->neighCovar)[0]=0, (pnt->neighCovar)[1]=0, (pnt->
neighCovar)[2]=0, (pnt->neighCovar)[3]=0, (pnt->neighCovar
)[4]=0, (pnt->neighCovar)[5]=0, (pnt->neighCovar)[6]=0,
(pnt->neighCovar)[7]=0, (pnt->neighCovar)[8]=0, (pnt->
neighCovar)[9]=0);
72 pnt->stability = 0.0;
73#if PULL_TANCOVAR1
74 ELL_6V_ZERO_SET(pnt->neighTanCovar)((pnt->neighTanCovar)[0]=(0), (pnt->neighTanCovar)[1]=(
0), (pnt->neighTanCovar)[2]=(0), (pnt->neighTanCovar)[3
]=(0), (pnt->neighTanCovar)[4]=(0), (pnt->neighTanCovar
)[5]=(0));
75#endif
76 pnt->neighInterNum = 0;
77 pnt->stuckIterNum = 0;
78#if PULL_PHIST0
79 pnt->phist = NULL((void*)0);
80 pnt->phistNum = 0;
81 pnt->phistArr = airArrayNew(AIR_CAST(void**, &(pnt->phist))((void**)(&(pnt->phist))),
82 &(pnt->phistNum),
83 _PHN6*sizeof(double), 32);
84#endif
85 pnt->status = 0;
86 ELL_4V_SET(pnt->pos, AIR_NAN, AIR_NAN, AIR_NAN, AIR_NAN)((pnt->pos)[0] = ((airFloatQNaN.f)), (pnt->pos)[1] = ((
airFloatQNaN.f)), (pnt->pos)[2] = ((airFloatQNaN.f)), (pnt
->pos)[3] = ((airFloatQNaN.f)));
87 pnt->energy = AIR_NAN(airFloatQNaN.f);
88 ELL_4V_SET(pnt->force, AIR_NAN, AIR_NAN, AIR_NAN, AIR_NAN)((pnt->force)[0] = ((airFloatQNaN.f)), (pnt->force)[1] =
((airFloatQNaN.f)), (pnt->force)[2] = ((airFloatQNaN.f)),
(pnt->force)[3] = ((airFloatQNaN.f)));
89 pnt->stepEnergy = pctx->sysParm.stepInitial;
90 pnt->stepConstr = pctx->sysParm.stepInitial;
91 for (ii=0; ii<pctx->infoTotalLen; ii++) {
92 pnt->info[ii] = AIR_NAN(airFloatQNaN.f);
93 }
94 return pnt;
95}
96
97pullPoint *
98pullPointNix(pullPoint *pnt) {
99
100 pnt->neighPointArr = airArrayNuke(pnt->neighPointArr);
101#if PULL_PHIST0
102 pnt->phistArr = airArrayNuke(pnt->phistArr);
103#endif
104 airFree(pnt);
105 return NULL((void*)0);
106}
107
108#if PULL_PHIST0
109void
110_pullPointHistInit(pullPoint *point) {
111
112 airArrayLenSet(point->phistArr, 0);
113 return;
114}
115
116void
117_pullPointHistAdd(pullPoint *point, int cond, double val) {
118 static const char me[]="_pullPointHistAdd";
119 unsigned int phistIdx;
120
121 phistIdx = airArrayLenIncr(point->phistArr, 1);
122 ELL_4V_COPY(point->phist + _PHN*phistIdx, point->pos)((point->phist + 6*phistIdx)[0] = (point->pos)[0], (point
->phist + 6*phistIdx)[1] = (point->pos)[1], (point->
phist + 6*phistIdx)[2] = (point->pos)[2], (point->phist
+ 6*phistIdx)[3] = (point->pos)[3]);
123
124 fprintf(stderr__stderrp, "!%s: point %p pos = %.17g %.17g %.17g %.17g (%g)\n", me,
125 point, point->pos[0], point->pos[1], point->pos[2], point->pos[3],
126 val);
127
128 (point->phist + _PHN6*phistIdx)[4] = cond;
129 (point->phist + _PHN6*phistIdx)[5] = val;
130 return;
131}
132#endif
133
134/*
135** HEY: there should be something like a "map" over all the points,
136** which could implement all these redundant functions
137*/
138
139unsigned int
140pullPointNumberFilter(const pullContext *pctx,
141 unsigned int idtagMin,
142 unsigned int idtagMax) {
143 unsigned int binIdx, pointNum;
144 const pullBin *bin;
145 const pullPoint *point;
146
147 pointNum = 0;
148 for (binIdx=0; binIdx<pctx->binNum; binIdx++) {
149 unsigned int pointIdx;
150 bin = pctx->bin + binIdx;
151 if (0 == idtagMin && 0 == idtagMax) {
152 pointNum += bin->pointNum;
153 } else {
154 for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) {
155 point = bin->point[pointIdx];
156 pointNum += (idtagMin <= point->idtag
157 && (0 == idtagMax
158 || point->idtag <= idtagMax));
159 }
160 }
161 }
162 return pointNum;
163}
164
165unsigned int
166pullPointNumber(const pullContext *pctx) {
167
168 return pullPointNumberFilter(pctx, 0, 0);
169}
170
171double
172_pullEnergyTotal(const pullContext *pctx) {
173 unsigned int binIdx, pointIdx;
174 const pullBin *bin;
175 const pullPoint *point;
176 double sum;
177
178 sum = 0;
179 for (binIdx=0; binIdx<pctx->binNum; binIdx++) {
180 bin = pctx->bin + binIdx;
181 for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) {
182 point = bin->point[pointIdx];
183 sum += point->energy;
184 }
185 }
186 return sum;
187}
188
189void
190_pullPointStepEnergyScale(pullContext *pctx, double scale) {
191 unsigned int binIdx, pointIdx;
192 const pullBin *bin;
193 pullPoint *point;
194
195 for (binIdx=0; binIdx<pctx->binNum; binIdx++) {
196 bin = pctx->bin + binIdx;
197 for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) {
198 point = bin->point[pointIdx];
199 point->stepEnergy *= scale;
200 point->stepEnergy = AIR_MIN(point->stepEnergy,((point->stepEnergy) < (3.40282347e+38F) ? (point->stepEnergy
) : (3.40282347e+38F))
201 _PULL_STEP_ENERGY_MAX)((point->stepEnergy) < (3.40282347e+38F) ? (point->stepEnergy
) : (3.40282347e+38F));
202 }
203 }
204 return;
205}
206
207double
208_pullStepInterAverage(const pullContext *pctx) {
209 unsigned int binIdx, pointIdx, pointNum;
210 const pullBin *bin;
211 const pullPoint *point;
212 double sum, avg;
213
214 sum = 0;
215 pointNum = 0;
216 for (binIdx=0; binIdx<pctx->binNum; binIdx++) {
217 bin = pctx->bin + binIdx;
218 pointNum += bin->pointNum;
219 for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) {
220 point = bin->point[pointIdx];
221 sum += point->stepEnergy;
222 }
223 }
224 avg = (!pointNum ? AIR_NAN(airFloatQNaN.f) : sum/pointNum);
225 return avg;
226}
227/* ^^^ vvv HEY HEY HEY: COPY + PASTE COPY + PASTE COPY + PASTE */
228double
229_pullStepConstrAverage(const pullContext *pctx) {
230 unsigned int binIdx, pointIdx, pointNum;
231 const pullBin *bin;
232 const pullPoint *point;
233 double sum, avg;
234
235 sum = 0;
236 pointNum = 0;
237 for (binIdx=0; binIdx<pctx->binNum; binIdx++) {
238 bin = pctx->bin + binIdx;
239 pointNum += bin->pointNum;
240 for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) {
241 point = bin->point[pointIdx];
242 sum += point->stepConstr;
243 }
244 }
245 avg = (!pointNum ? AIR_NAN(airFloatQNaN.f) : sum/pointNum);
246 return avg;
247}
248
249/*
250** convenience function for learning a scalar AND its gradient or hessian
251**
252*/
253double
254pullPointScalar(const pullContext *pctx, const pullPoint *point, int sclInfo,
255 /* output */
256 double grad[3], double hess[9]) {
257 static const char me[]="pullPointScalar";
258 double scl;
259 const pullInfoSpec *ispec;
260 int gradInfo[1+PULL_INFO_MAX23] = {
261 0, /* pullInfoUnknown */
262 0, /* pullInfoTensor */
263 0, /* pullInfoTensorInverse */
264 0, /* pullInfoHessian */
265 pullInfoInsideGradient, /* pullInfoInside */
266 0, /* pullInfoInsideGradient */
267 pullInfoHeightGradient, /* pullInfoHeight */
268 0, /* pullInfoHeightGradient */
269 0, /* pullInfoHeightHessian */
270 0, /* pullInfoHeightLaplacian */
271 0, /* pullInfoSeedPreThresh */
272 0, /* pullInfoSeedThresh */
273 0, /* pullInfoLiveThresh */
274 0, /* pullInfoLiveThresh2 */
275 0, /* pullInfoLiveThresh3 */
276 0, /* pullInfoTangent1 */
277 0, /* pullInfoTangent2 */
278 0, /* pullInfoNegativeTangent1 */
279 0, /* pullInfoNegativeTangent2 */
280 pullInfoIsovalueGradient, /* pullInfoIsovalue */
281 0, /* pullInfoIsovalueGradient */
282 0, /* pullInfoIsovalueHessian */
283 0, /* pullInfoStrength */
284 };
285 int hessInfo[1+PULL_INFO_MAX23] = {
286 0, /* pullInfoUnknown */
287 0, /* pullInfoTensor */
288 0, /* pullInfoTensorInverse */
289 0, /* pullInfoHessian */
290 0, /* pullInfoInside */
291 0, /* pullInfoInsideGradient */
292 pullInfoHeightHessian, /* pullInfoHeight */
293 0, /* pullInfoHeightGradient */
294 0, /* pullInfoHeightHessian */
295 0, /* pullInfoHeightLaplacian */
296 0, /* pullInfoSeedPreThresh */
297 0, /* pullInfoSeedThresh */
298 0, /* pullInfoLiveThresh */
299 0, /* pullInfoLiveThresh2 */
300 0, /* pullInfoLiveThresh3 */
301 0, /* pullInfoTangent1 */
302 0, /* pullInfoTangent2 */
303 0, /* pullInfoNegativeTangent1 */
304 0, /* pullInfoNegativeTangent2 */
305 pullInfoIsovalueHessian, /* pullInfoIsovalue */
306 0, /* pullInfoIsovalueGradient */
307 0, /* pullInfoIsovalueHessian */
308 0, /* pullInfoStrength */
309 };
310 const unsigned int *infoIdx;
311
312 infoIdx = pctx->infoIdx;
313 ispec = pctx->ispec[sclInfo];
314 /* NB: this "scl" is not scale-space scale; its the scaling
315 of the scalar. this is getting confusing ... */
316 scl = point->info[infoIdx[sclInfo]];
317 scl = (scl - ispec->zero)*ispec->scale;
318 /* if (289 == pctx->iter) {
319 fprintf(stderr, "!%s(%04u,%s)@(%g,%g): (%g - %g)*%g == %g\n",
320 me, point->idtag, airEnumStr(pullInfo, sclInfo),
321 point->pos[0], point->pos[1],
322 point->info[infoIdx[sclInfo]], ispec->zero, ispec->scale, scl);
323 } */
324 if (0 && _pullVerbose) {
325 if (pullInfoSeedThresh == sclInfo) {
326 printf("!%s: seed thresh (%g - %g)*%g == %g\n", me,
327 point->info[infoIdx[sclInfo]], ispec->zero, ispec->scale, scl);
328 }
329 }
330 /*
331 learned: this wasn't thought through: the idea was that the height
332 *laplacian* answer should be transformed by the *height* zero and
333 scale. scale might make sense, but not zero. This cost a few
334 hours of tracking down the fact that the first zero-crossing
335 detection phase of the lapl constraint was failing because the
336 laplacian was vacillating around hspec->zero, not 0.0 . . .
337 if (pullInfoHeightLaplacian == sclInfo) {
338 const pullInfoSpec *hspec;
339 hspec = pctx->ispec[pullInfoHeight];
340 scl = (scl - hspec->zero)*hspec->scale;
341 } else {
342 scl = (scl - ispec->zero)*ispec->scale;
343 }
344 */
345 /*
346 fprintf(stderr, "!%s: %s = (%g - %g)*%g = %g*%g = %g = %g\n", me,
347 airEnumStr(pullInfo, sclInfo),
348 point->info[infoIdx[sclInfo]],
349 ispec->zero, ispec->scale,
350 point->info[infoIdx[sclInfo]] - ispec->zero, ispec->scale,
351 (point->info[infoIdx[sclInfo]] - ispec->zero)*ispec->scale,
352 scl);
353 */
354 if (grad && gradInfo[sclInfo]) {
355 const double *ptr = point->info + infoIdx[gradInfo[sclInfo]];
356 ELL_3V_SCALE(grad, ispec->scale, ptr)((grad)[0] = (ispec->scale)*(ptr)[0], (grad)[1] = (ispec->
scale)*(ptr)[1], (grad)[2] = (ispec->scale)*(ptr)[2]);
357 }
358 if (hess && hessInfo[sclInfo]) {
359 const double *ptr = point->info + infoIdx[hessInfo[sclInfo]];
360 ELL_3M_SCALE(hess, ispec->scale, ptr)((((hess)+0)[0] = ((ispec->scale))*((ptr)+0)[0], ((hess)+0
)[1] = ((ispec->scale))*((ptr)+0)[1], ((hess)+0)[2] = ((ispec
->scale))*((ptr)+0)[2]), (((hess)+3)[0] = ((ispec->scale
))*((ptr)+3)[0], ((hess)+3)[1] = ((ispec->scale))*((ptr)+3
)[1], ((hess)+3)[2] = ((ispec->scale))*((ptr)+3)[2]), (((hess
)+6)[0] = ((ispec->scale))*((ptr)+6)[0], ((hess)+6)[1] = (
(ispec->scale))*((ptr)+6)[1], ((hess)+6)[2] = ((ispec->
scale))*((ptr)+6)[2]));
361 }
362 return scl;
363}
364
365int
366pullProbe(pullTask *task, pullPoint *point) {
367 static const char me[]="pullProbe";
368 unsigned int ii, gret=0;
369 int edge;
370 /*
371 fprintf(stderr, "!%s: task->probeSeedPreThreshOnly = %d\n", me,
372 task->probeSeedPreThreshOnly);
373 */
374#if 0
375 static int opened=AIR_FALSE0;
376 static FILE *flog;
377#endif
378
379 /*
380 fprintf(stderr, "%s(%u,%u): A volNum = %u\n", me, task->pctx->iter, point->idtag,task->pctx->volNum);
381 */
382#if 0
383 static int logIdx=0, logDone=AIR_FALSE0, logStarted=AIR_FALSE0;
384 static Nrrd *nlog;
385 static double *log=NULL((void*)0);
386 if (!logStarted) {
387 if (81 == point->idtag) {
388 printf("\n\n%s: ###### HELLO begin logging . . .\n\n\n", me);
389 /* knowing the logIdx at the end of logging . . . */
390 nlog = nrrdNew();
391 nrrdMaybeAlloc_va(nlog, nrrdTypeDouble, 2,
392 AIR_CAST(size_t, 25)((size_t)(25)),
393 AIR_CAST(size_t, 2754)((size_t)(2754)));
394 log = AIR_CAST(double*, nlog->data)((double*)(nlog->data));
395 logStarted = AIR_TRUE1;
396 }
397 }
398#endif
399
400 if (!ELL_4V_EXISTS(point->pos)((((int)(!(((point->pos)[0]) - ((point->pos)[0]))))) &&
(((int)(!(((point->pos)[1]) - ((point->pos)[1]))))) &&
(((int)(!(((point->pos)[2]) - ((point->pos)[2]))))) &&
(((int)(!(((point->pos)[3]) - ((point->pos)[3]))))))) {
401 fprintf(stderr__stderrp, "%s: pnt %u non-exist pos (%g,%g,%g,%g)\n\n!!!\n\n\n",
402 me, point->idtag, point->pos[0], point->pos[1],
403 point->pos[2], point->pos[3]);
404 biffAddf(PULLpullBiffKey, "%s: pnt %u non-exist pos (%g,%g,%g,%g)",
405 me, point->idtag, point->pos[0], point->pos[1],
406 point->pos[2], point->pos[3]);
407 return 1;
408 /* can't probe, but make it go away as quickly as possible */
409 /*
410 ELL_4V_SET(point->pos, 0, 0, 0, 0);
411 point->status |= PULL_STATUS_NIXME_BIT;
412 return 0;
413 */
414 }
415 if (task->pctx->verbose > 3) {
416 printf("%s: hello; probing %u volumes\n", me, task->pctx->volNum);
417 }
418 edge = AIR_FALSE0;
419 task->pctx->count[pullCountProbe] += 1;
420 /*
421 fprintf(stderr, "%s(%u,%u): B volNum = %u\n", me, task->pctx->iter, point->idtag,task->pctx->volNum);
422 */
423 for (ii=0; ii<task->pctx->volNum; ii++) {
424 pullVolume *vol;
425 vol = task->vol[ii];
426 if (task->probeSeedPreThreshOnly
427 && !(vol->forSeedPreThresh)) {
428 /* we're here *only* to probe SeedPreThresh,
429 and this volume isn't used for that */
430 continue;
431 }
432 if (task->pctx->iter && vol->seedOnly) {
433 /* its after the 1st iteration (#0), and this vol is only for seeding */
434 continue;
435 }
436 /* HEY should task->vol[ii]->scaleNum be the using-scale-space test? */
437 if (!task->vol[ii]->ninScale) {
438 gret = gageProbeSpace(task->vol[ii]->gctx,
439 point->pos[0], point->pos[1], point->pos[2],
440 AIR_FALSE0 /* index-space */,
441 AIR_TRUE1 /* clamp */);
442 } else {
443 if (task->pctx->verbose > 3) {
444 printf("%s: vol[%u] has scale (%u)-> "
445 "gageStackProbeSpace(%p) (v %d)\n",
446 me, ii, task->vol[ii]->scaleNum,
447 AIR_VOIDP(task->vol[ii]->gctx)((void *)(task->vol[ii]->gctx)),
448 task->vol[ii]->gctx->verbose);
449 }
450 /*
451 if (81 == point->idtag) {
452 printf("%s: probing vol[%u] @ %g %g %g %g\n", me, ii,
453 point->pos[0], point->pos[1], point->pos[2], point->pos[3]);
454 }
455 */
456 gret = gageStackProbeSpace(task->vol[ii]->gctx,
457 point->pos[0], point->pos[1], point->pos[2],
458 (task->pctx->flag.scaleIsTau
459 ? gageSigOfTau(point->pos[3])airSigmaOfTau(point->pos[3])
460 : point->pos[3]),
461 AIR_FALSE0 /* index-space */,
462 AIR_TRUE1 /* clamp */);
463 }
464 if (gret) {
465 biffAddf(PULLpullBiffKey, "%s: probe failed on vol %u/%u: (%d) %s", me,
466 ii, task->pctx->volNum,
467 task->vol[ii]->gctx->errNum, task->vol[ii]->gctx->errStr);
468 return 1;
469 }
470 /*
471 if (!edge
472 && AIR_ABS(point->pos[1] - 67) < 1
473 && AIR_ABS(point->pos[2] - 67) < 1
474 && point->pos[3] > 3.13
475 && !!task->vol[ii]->gctx->edgeFrac) {
476 fprintf(stderr, "!%s(%u @ %g,%g,%g,%g): "
477 "vol[%u]->gctx->edgeFrac %g => edge bit on\n",
478 me, point->idtag,
479 point->pos[0], point->pos[1],
480 point->pos[2], point->pos[3],
481 ii, task->vol[ii]->gctx->edgeFrac);
482 }
483 */
484 edge |= !!task->vol[ii]->gctx->edgeFrac;
485 }
486 if (edge) {
487 point->status |= PULL_STATUS_EDGE_BIT(1<< 4);
488 } else {
489 point->status &= ~PULL_STATUS_EDGE_BIT(1<< 4);
490 }
491
492 /* maybe is a little stupid to have the infos indexed this way,
493 since it means that we always have to loop through all indices,
494 but at least the compiler can unroll it . . . */
495 for (ii=0; ii<=PULL_INFO_MAX23; ii++) {
496 unsigned int alen, aidx;
497 const pullInfoSpec *ispec;
498 ispec = task->pctx->ispec[ii];
499 if (ispec) {
500 alen = _pullInfoLen[ii];
501 aidx = task->pctx->infoIdx[ii];
502 if (pullSourceGage == ispec->source) {
503 _pullInfoCopy[alen](point->info + aidx, task->ans[ii]);
504 /* if (289 == task->pctx->iter) {
505 fprintf(stderr, "!%s(%u): copied info %u (%s) len %u\n", me, point->idtag,
506 ii, airEnumStr(pullInfo, ii), alen);
507 if (1 == alen) {
508 fprintf(stderr, "!%s(%u): (point->info + %u)[0] = %g\n", me, point->idtag,
509 aidx, (point->info + aidx)[0]);
510 }
511 } */
512 /*
513 if (81 == point->idtag) {
514 pullVolume *vol;
515 pullInfoSpec *isp;
516 isp = task->pctx->ispec[ii];
517 vol = task->pctx->vol[isp->volIdx];
518 if (1 == alen) {
519 printf("!%s: info[%u] %s: %s(\"%s\") = %g\n", me,
520 ii, airEnumStr(pullInfo, ii),
521 airEnumStr(vol->kind->enm, isp->item),
522 vol->name, task->ans[ii][0]);
523 } else {
524 unsigned int vali;
525 printf("!%s: info[%u] %s: %s(\"%s\") =\n", me,
526 ii, airEnumStr(pullInfo, ii),
527 airEnumStr(vol->kind->enm, isp->item), vol->name);
528 for (vali=0; vali<alen; vali++) {
529 printf("!%s: [%u] %g\n", me, vali,
530 task->ans[ii][vali]);
531 }
532 }
533 }
534 */
535 } else if (pullSourceProp == ispec->source) {
536 switch (ispec->prop) {
537 case pullPropIdtag:
538 point->info[aidx] = point->idtag;
539 break;
540 case pullPropIdCC:
541 point->info[aidx] = point->idCC;
542 break;
543 case pullPropEnergy:
544 point->info[aidx] = point->energy;
545 break;
546 case pullPropStepEnergy:
547 point->info[aidx] = point->stepEnergy;
548 break;
549 case pullPropStepConstr:
550 point->info[aidx] = point->stepConstr;
551 break;
552 case pullPropStuck:
553 point->info[aidx] = ((point->status & PULL_STATUS_STUCK_BIT(1<< 1))
554 ? point->stuckIterNum
555 : 0);
556 break;
557 case pullPropPosition:
558 ELL_4V_COPY(point->info + aidx, point->pos)((point->info + aidx)[0] = (point->pos)[0], (point->
info + aidx)[1] = (point->pos)[1], (point->info + aidx)
[2] = (point->pos)[2], (point->info + aidx)[3] = (point
->pos)[3]);
559 break;
560 case pullPropForce:
561 ELL_4V_COPY(point->info + aidx, point->force)((point->info + aidx)[0] = (point->force)[0], (point->
info + aidx)[1] = (point->force)[1], (point->info + aidx
)[2] = (point->force)[2], (point->info + aidx)[3] = (point
->force)[3]);
562 break;
563 case pullPropNeighDistMean:
564 point->info[aidx] = point->neighDistMean;
565 break;
566 case pullPropScale:
567 point->info[aidx] = point->pos[3];
568 break;
569 case pullPropNeighCovar:
570 ELL_10V_COPY(point->info + aidx, point->neighCovar)((point->info + aidx)[0] = (point->neighCovar)[0], (point
->info + aidx)[1] = (point->neighCovar)[1], (point->
info + aidx)[2] = (point->neighCovar)[2], (point->info +
aidx)[3] = (point->neighCovar)[3], (point->info + aidx
)[4] = (point->neighCovar)[4], (point->info + aidx)[5] =
(point->neighCovar)[5], (point->info + aidx)[6] = (point
->neighCovar)[6], (point->info + aidx)[7] = (point->
neighCovar)[7], (point->info + aidx)[8] = (point->neighCovar
)[8], (point->info + aidx)[9] = (point->neighCovar)[9]);
571 break;
572 case pullPropNeighCovar7Ten:
573 TEN_T_SET(point->info + aidx, 1.0f,( (point->info + aidx)[0] = (1.0f), (point->info + aidx
)[1] = (point->neighCovar[0]), (point->info + aidx)[2] =
(point->neighCovar[1]), (point->info + aidx)[3] = (point
->neighCovar[2]), (point->info + aidx)[4] = (point->
neighCovar[4]), (point->info + aidx)[5] = (point->neighCovar
[5]), (point->info + aidx)[6] = (point->neighCovar[7]) )
574 point->neighCovar[0],( (point->info + aidx)[0] = (1.0f), (point->info + aidx
)[1] = (point->neighCovar[0]), (point->info + aidx)[2] =
(point->neighCovar[1]), (point->info + aidx)[3] = (point
->neighCovar[2]), (point->info + aidx)[4] = (point->
neighCovar[4]), (point->info + aidx)[5] = (point->neighCovar
[5]), (point->info + aidx)[6] = (point->neighCovar[7]) )
575 point->neighCovar[1],( (point->info + aidx)[0] = (1.0f), (point->info + aidx
)[1] = (point->neighCovar[0]), (point->info + aidx)[2] =
(point->neighCovar[1]), (point->info + aidx)[3] = (point
->neighCovar[2]), (point->info + aidx)[4] = (point->
neighCovar[4]), (point->info + aidx)[5] = (point->neighCovar
[5]), (point->info + aidx)[6] = (point->neighCovar[7]) )
576 point->neighCovar[2],( (point->info + aidx)[0] = (1.0f), (point->info + aidx
)[1] = (point->neighCovar[0]), (point->info + aidx)[2] =
(point->neighCovar[1]), (point->info + aidx)[3] = (point
->neighCovar[2]), (point->info + aidx)[4] = (point->
neighCovar[4]), (point->info + aidx)[5] = (point->neighCovar
[5]), (point->info + aidx)[6] = (point->neighCovar[7]) )
577 point->neighCovar[4],( (point->info + aidx)[0] = (1.0f), (point->info + aidx
)[1] = (point->neighCovar[0]), (point->info + aidx)[2] =
(point->neighCovar[1]), (point->info + aidx)[3] = (point
->neighCovar[2]), (point->info + aidx)[4] = (point->
neighCovar[4]), (point->info + aidx)[5] = (point->neighCovar
[5]), (point->info + aidx)[6] = (point->neighCovar[7]) )
578 point->neighCovar[5],( (point->info + aidx)[0] = (1.0f), (point->info + aidx
)[1] = (point->neighCovar[0]), (point->info + aidx)[2] =
(point->neighCovar[1]), (point->info + aidx)[3] = (point
->neighCovar[2]), (point->info + aidx)[4] = (point->
neighCovar[4]), (point->info + aidx)[5] = (point->neighCovar
[5]), (point->info + aidx)[6] = (point->neighCovar[7]) )
579 point->neighCovar[7])( (point->info + aidx)[0] = (1.0f), (point->info + aidx
)[1] = (point->neighCovar[0]), (point->info + aidx)[2] =
(point->neighCovar[1]), (point->info + aidx)[3] = (point
->neighCovar[2]), (point->info + aidx)[4] = (point->
neighCovar[4]), (point->info + aidx)[5] = (point->neighCovar
[5]), (point->info + aidx)[6] = (point->neighCovar[7]) );
580 break;
581#if PULL_TANCOVAR1
582 case pullPropNeighTanCovar:
583 TEN_T_SET(point->info + aidx, 1.0f,( (point->info + aidx)[0] = (1.0f), (point->info + aidx
)[1] = (point->neighTanCovar[0]), (point->info + aidx)[
2] = (point->neighTanCovar[1]), (point->info + aidx)[3]
= (point->neighTanCovar[2]), (point->info + aidx)[4] =
(point->neighTanCovar[3]), (point->info + aidx)[5] = (
point->neighTanCovar[4]), (point->info + aidx)[6] = (point
->neighTanCovar[5]) )
584 point->neighTanCovar[0],( (point->info + aidx)[0] = (1.0f), (point->info + aidx
)[1] = (point->neighTanCovar[0]), (point->info + aidx)[
2] = (point->neighTanCovar[1]), (point->info + aidx)[3]
= (point->neighTanCovar[2]), (point->info + aidx)[4] =
(point->neighTanCovar[3]), (point->info + aidx)[5] = (
point->neighTanCovar[4]), (point->info + aidx)[6] = (point
->neighTanCovar[5]) )
585 point->neighTanCovar[1],( (point->info + aidx)[0] = (1.0f), (point->info + aidx
)[1] = (point->neighTanCovar[0]), (point->info + aidx)[
2] = (point->neighTanCovar[1]), (point->info + aidx)[3]
= (point->neighTanCovar[2]), (point->info + aidx)[4] =
(point->neighTanCovar[3]), (point->info + aidx)[5] = (
point->neighTanCovar[4]), (point->info + aidx)[6] = (point
->neighTanCovar[5]) )
586 point->neighTanCovar[2],( (point->info + aidx)[0] = (1.0f), (point->info + aidx
)[1] = (point->neighTanCovar[0]), (point->info + aidx)[
2] = (point->neighTanCovar[1]), (point->info + aidx)[3]
= (point->neighTanCovar[2]), (point->info + aidx)[4] =
(point->neighTanCovar[3]), (point->info + aidx)[5] = (
point->neighTanCovar[4]), (point->info + aidx)[6] = (point
->neighTanCovar[5]) )
587 point->neighTanCovar[3],( (point->info + aidx)[0] = (1.0f), (point->info + aidx
)[1] = (point->neighTanCovar[0]), (point->info + aidx)[
2] = (point->neighTanCovar[1]), (point->info + aidx)[3]
= (point->neighTanCovar[2]), (point->info + aidx)[4] =
(point->neighTanCovar[3]), (point->info + aidx)[5] = (
point->neighTanCovar[4]), (point->info + aidx)[6] = (point
->neighTanCovar[5]) )
588 point->neighTanCovar[4],( (point->info + aidx)[0] = (1.0f), (point->info + aidx
)[1] = (point->neighTanCovar[0]), (point->info + aidx)[
2] = (point->neighTanCovar[1]), (point->info + aidx)[3]
= (point->neighTanCovar[2]), (point->info + aidx)[4] =
(point->neighTanCovar[3]), (point->info + aidx)[5] = (
point->neighTanCovar[4]), (point->info + aidx)[6] = (point
->neighTanCovar[5]) )
589 point->neighTanCovar[5])( (point->info + aidx)[0] = (1.0f), (point->info + aidx
)[1] = (point->neighTanCovar[0]), (point->info + aidx)[
2] = (point->neighTanCovar[1]), (point->info + aidx)[3]
= (point->neighTanCovar[2]), (point->info + aidx)[4] =
(point->neighTanCovar[3]), (point->info + aidx)[5] = (
point->neighTanCovar[4]), (point->info + aidx)[6] = (point
->neighTanCovar[5]) );
590 break;
591#endif
592 case pullPropStability:
593 point->info[aidx] = point->stability;
594 break;
595 default:
596 break;
597 }
598 }
599 }
600 }
601
602#if 0
603 if (logStarted && !logDone) {
604 unsigned int ai;
605 /* the actual logging */
606 log[0] = point->idtag;
607 ELL_4V_COPY(log + 1, point->pos)((log + 1)[0] = (point->pos)[0], (log + 1)[1] = (point->
pos)[1], (log + 1)[2] = (point->pos)[2], (log + 1)[3] = (point
->pos)[3]);
608 for (ai=0; ai<20; ai++) {
609 log[5 + ai] = point->info[ai];
610 }
611 log += nlog->axis[0].size;
612 logIdx++;
613 if (1 == task->pctx->iter && 81 == point->idtag) {
614 printf("\n\n%s: ###### OKAY done logging (%u). . .\n\n\n", me, logIdx);
615 nrrdSave("probelog.nrrd", nlog, NULL((void*)0));
616 nlog = nrrdNuke(nlog);
617 logDone = AIR_TRUE1;
618 }
619 }
620#endif
621
622
623#if 0
624 if (!opened) {
625 flog = fopen("flog.txt", "w");
626 opened = AIR_TRUE1;
627 }
628 if (opened) {
629 fprintf(flog, "%s(%u): spthr(%g,%g,%g,%g) = %g\n", me, point->idtag,
630 point->pos[0], point->pos[1], point->pos[2], point->pos[3],
631 point->info[task->pctx->infoIdx[pullInfoSeedPreThresh]]);
632 }
633#endif
634
635 return 0;
636}
637
638static int
639_threshFail(const pullContext *pctx, const pullPoint *point, int info) {
640 /* static const char me[]="_threshFail"; */
641 double val;
642 int ret;
643
644 if (pctx->ispec[info]) {
645 val = pullPointScalar(pctx, point, info, NULL((void*)0), NULL((void*)0));
646 ret = (val < 0);
647 /*
648 fprintf(stderr, "%s(%s): val=%g -> ret=%d\n", me,
649 airEnumStr(pullInfo, info), val, ret);
650 */
651 } else {
652 ret = AIR_FALSE0;
653 }
654 return ret;
655}
656
657int
658pullPointInitializePerVoxel(const pullContext *pctx,
659 const unsigned int pointIdx,
660 pullPoint *point, pullVolume *scaleVol,
661 /* output */
662 int *createFailP) {
663 static const char me[]="pullPointInitializePerVoxel";
664 unsigned int vidx[3], pix;
665 double iPos[3];
666 airRandMTState *rng;
667 pullVolume *seedVol;
668 gageShape *seedShape;
669 int reject, rejectEdge, constrFail;
670 unsigned int k;
671
672 seedVol = pctx->vol[pctx->ispec[pullInfoSeedThresh]->volIdx];
673 seedShape = seedVol->gctx->shape;
674 rng = pctx->task[0]->rng;
675
676 /* Obtain voxel and indices from pointIdx */
677 /* axis ordering for this is x, y, z, scale */
678 pix = pointIdx;
679 if (pctx->initParm.pointPerVoxel > 0) {
680 pix /= pctx->initParm.pointPerVoxel;
681 } else {
682 pix *= -pctx->initParm.pointPerVoxel;
683 }
684 vidx[0] = pix % seedShape->size[0];
685 pix = (pix - vidx[0])/seedShape->size[0];
686 vidx[1] = pix % seedShape->size[1];
687 pix = (pix - vidx[1])/seedShape->size[1];
688 if (pctx->initParm.ppvZRange[0] <= pctx->initParm.ppvZRange[1]) {
689 unsigned int zrn;
690 zrn = pctx->initParm.ppvZRange[1] - pctx->initParm.ppvZRange[0] + 1;
691 vidx[2] = (pix % zrn) + pctx->initParm.ppvZRange[0];
692 pix = (pix - (pix % zrn))/zrn;
693 } else {
694 vidx[2] = pix % seedShape->size[2];
695 pix = (pix - vidx[2])/seedShape->size[2];
696 }
697 for (k=0; k<=2; k++) {
698 iPos[k] = vidx[k] + pctx->initParm.jitter*(airDrandMT_r(rng)-0.5);
699 }
700 gageShapeItoW(seedShape, point->pos, iPos);
701 if (pctx->flag.zeroZ) {
702 point->pos[2] = 0.0;
703 }
704
705 if (0 && _pullVerbose) {
706 printf("!%s: pointIdx %u -> vidx %u %u %u (%u)\n"
707 " -> iPos %g %g %g -> wPos %g %g %g\n",
708 me, pointIdx, vidx[0], vidx[1], vidx[2], pix,
709 iPos[0], iPos[1], iPos[2],
710 point->pos[0], point->pos[1], point->pos[2]);
711 }
712
713 /* Compute sigma coordinate from pix */
714 if (pctx->haveScale) {
715 int outside;
716 double aidx, bidx;
717 /* pix should already be integer in [0, pctx->samplesAlongScaleNum-1)]. */
718 aidx = pix + pctx->initParm.jitter*(airDrandMT_r(rng)-0.5);
719 bidx = AIR_AFFINE(-0.5, aidx, pctx->initParm.samplesAlongScaleNum-0.5,( ((double)(scaleVol->scaleNum-1)-(0.0))*((double)(aidx)-(
-0.5)) / ((double)(pctx->initParm.samplesAlongScaleNum-0.5
)-(-0.5)) + (0.0))
720 0.0, scaleVol->scaleNum-1)( ((double)(scaleVol->scaleNum-1)-(0.0))*((double)(aidx)-(
-0.5)) / ((double)(pctx->initParm.samplesAlongScaleNum-0.5
)-(-0.5)) + (0.0));
721 point->pos[3] = gageStackItoW(scaleVol->gctx, bidx, &outside);
722 if (pctx->flag.scaleIsTau) {
723 point->pos[3] = gageTauOfSig(point->pos[3])airTauOfSigma(point->pos[3]);
724 }
725 if (0 && _pullVerbose) {
726 printf("!%s(%u): pix %u -> a %g b %g -> wpos %g\n", me, point->idtag,
727 pix, aidx, bidx, point->pos[3]);
728 }
729 } else {
730 point->pos[3] = 0;
731 }
732
733 if (pctx->ispec[pullInfoSeedPreThresh]) {
734 /* we first do a special-purpose probe just for SeedPreThresh */
735 pctx->task[0]->probeSeedPreThreshOnly = AIR_TRUE1;
736 if (pullProbe(pctx->task[0], point)) {
737 biffAddf(PULLpullBiffKey, "%s: pre-probing pointIdx %u of world", me, pointIdx);
738 return 1;
739 }
740 pctx->task[0]->probeSeedPreThreshOnly = AIR_FALSE0;
741 if (_threshFail(pctx, point, pullInfoSeedPreThresh)) {
742 reject = AIR_TRUE1;
743 /* HEY! this obviously need to be re-written */
744 goto finish;
745 }
746 }
747 /* else, we didn't have a SeedPreThresh, or we did, and we passed
748 it. Now we REDO the probe, including possibly re-learning
749 SeedPreThresh, which is silly, but the idea is that this is a
750 small price compared to what has been saved by avoiding all the
751 gageProbe()s on volumes unrelated to SeedPreThresh */
752 if (pullProbe(pctx->task[0], point)) {
753 biffAddf(PULLpullBiffKey, "%s: probing pointIdx %u of world", me, pointIdx);
754 return 1;
755 }
756
757 constrFail = AIR_FALSE0;
758 reject = AIR_FALSE0;
759
760 /* Check we pass pre-threshold */
761 if (!reject) reject |= _threshFail(pctx, point, pullInfoSeedPreThresh);
762
763 if (!pctx->flag.constraintBeforeSeedThresh) {
764 /* we should be guaranteed to have a seed thresh info */
765 if (!reject) reject |= _threshFail(pctx, point, pullInfoSeedThresh);
766 if (pctx->initParm.liveThreshUse) {
767 if (!reject) reject |= _threshFail(pctx, point, pullInfoLiveThresh);
768 if (!reject) reject |= _threshFail(pctx, point, pullInfoLiveThresh2);
769 if (!reject) reject |= _threshFail(pctx, point, pullInfoLiveThresh3);
770 }
771 }
772
773 if (!reject && pctx->constraint) {
774 if (_pullConstraintSatisfy(pctx->task[0], point,
775 10*_PULL_CONSTRAINT_TRAVEL_MAX2,
776 &constrFail)) {
777 biffAddf(PULLpullBiffKey, "%s: on pnt %u",
778 me, pointIdx);
779 return 1;
780 }
781 reject |= constrFail;
782 /* post constraint-satisfaction, we certainly have to assert thresholds */
783 if (!reject) reject |= _threshFail(pctx, point, pullInfoSeedThresh);
784 if (pctx->initParm.liveThreshUse) {
785 if (!reject) reject |= _threshFail(pctx, point, pullInfoLiveThresh);
786 if (!reject) reject |= _threshFail(pctx, point, pullInfoLiveThresh2);
787 if (!reject) reject |= _threshFail(pctx, point, pullInfoLiveThresh3);
788 }
789 if (pctx->flag.nixAtVolumeEdgeSpace
790 && (point->status & PULL_STATUS_EDGE_BIT(1<< 4))) {
791 rejectEdge = AIR_TRUE1;
792 } else {
793 rejectEdge = AIR_FALSE0;
794 }
795 reject |= rejectEdge;
796 if (pctx->verbose > 1) {
797 fprintf(stderr__stderrp, "%s(%u): constr %d, seed %d, thresh %d %d %d, edge %d\n",
798 me, point->idtag, constrFail,
799 _threshFail(pctx, point, pullInfoSeedThresh),
800 _threshFail(pctx, point, pullInfoLiveThresh),
801 _threshFail(pctx, point, pullInfoLiveThresh2),
802 _threshFail(pctx, point, pullInfoLiveThresh3), rejectEdge);
803 }
804 } else {
805 constrFail = AIR_FALSE0;
806 }
807
808 finish:
809 /* Gather consensus */
810 if (reject) {
811 *createFailP = AIR_TRUE1;
812 } else {
813 *createFailP = AIR_FALSE0;
814 }
815
816 return 0;
817}
818
819static void
820_pullUnitToWorld(const pullContext *pctx, const pullVolume *scaleVol,
821 double wrld[4], const double unit[4]) {
822 /* static const char me[]="_pullUnitToWorld"; */
823
824 wrld[0] = AIR_AFFINE(0.0, unit[0], 1.0, pctx->bboxMin[0], pctx->bboxMax[0])( ((double)(pctx->bboxMax[0])-(pctx->bboxMin[0]))*((double
)(unit[0])-(0.0)) / ((double)(1.0)-(0.0)) + (pctx->bboxMin
[0]));
825 wrld[1] = AIR_AFFINE(0.0, unit[1], 1.0, pctx->bboxMin[1], pctx->bboxMax[1])( ((double)(pctx->bboxMax[1])-(pctx->bboxMin[1]))*((double
)(unit[1])-(0.0)) / ((double)(1.0)-(0.0)) + (pctx->bboxMin
[1]));
826 wrld[2] = AIR_AFFINE(0.0, unit[2], 1.0, pctx->bboxMin[2], pctx->bboxMax[2])( ((double)(pctx->bboxMax[2])-(pctx->bboxMin[2]))*((double
)(unit[2])-(0.0)) / ((double)(1.0)-(0.0)) + (pctx->bboxMin
[2]));
827 if (pctx->haveScale) {
828 double sridx;
829 int outside;
830 sridx = AIR_AFFINE(0.0, unit[3], 1.0, 0, scaleVol->scaleNum-1)( ((double)(scaleVol->scaleNum-1)-(0))*((double)(unit[3])-
(0.0)) / ((double)(1.0)-(0.0)) + (0));
831 wrld[3] = gageStackItoW(scaleVol->gctx, sridx, &outside);
832 if (pctx->flag.scaleIsTau) {
833 wrld[3] = gageTauOfSig(wrld[3])airTauOfSigma(wrld[3]);
834 }
835 } else {
836 wrld[3] = 0.0;
837 }
838 /*
839 fprintf(stderr, "!%s: (%g,%g,%g,%g) -- [%g,%g]x[%g,%g]x[%g,%g]--> (%g,%g,%g,%g)\n", me,
840 unit[0], unit[1], unit[2], unit[3],
841 pctx->bboxMin[0], pctx->bboxMin[1], pctx->bboxMin[2],
842 pctx->bboxMax[0], pctx->bboxMax[1], pctx->bboxMax[2],
843 wrld[0], wrld[1], wrld[2], wrld[3]);
844 */
845 return;
846}
847
848int
849pullPointInitializeRandomOrHalton(pullContext *pctx,
850 const unsigned int pointIdx,
851 pullPoint *point, pullVolume *scaleVol) {
852 static const char me[]="pullPointInitializeRandomOrHalton";
853 int reject, verbo;
854 airRandMTState *rng;
855 unsigned int tryCount = 0, threshFailCount = 0, spthreshFailCount = 0,
856 constrFailCount = 0;
857 rng = pctx->task[0]->rng;
858
859 do {
860 double rpos[4];
861 /* fprintf(stderr, "!%s: starting doooo (tryCount %u)!\n", me, tryCount); */
862 tryCount++;
863 reject = AIR_FALSE0;
864 _pullPointHistInit(point);
865 /* Populate tentative random point */
866 if (pullInitMethodHalton == pctx->initParm.method) {
867 /* we generate all 4 coordinates, even if we don't need them all */
868 airHalton(rpos, (pointIdx + threshFailCount + constrFailCount
869 + pctx->haltonOffset + pctx->initParm.haltonStartIndex),
870 airPrimeList, 4);
871 /*
872 fprintf(stderr, "!%s(%u/%u): halton(%u=%u+%u+%u+%u+%u) => %g %g %g %g\n",
873 me, pointIdx, pctx->idtagNext,
874 (pointIdx + threshFailCount + constrFailCount +
875 pctx->haltonOffset + pctx->initParm.haltonStartIndex),
876 pointIdx, threshFailCount, constrFailCount,
877 pctx->haltonOffset, pctx->initParm.haltonStartIndex,
878 rpos[0], rpos[1], rpos[2], rpos[3]);
879 */
880 /*
881 fprintf(stderr, "%g %g %g %g ",
882 rpos[0], rpos[1], rpos[2], rpos[3]);
883 */
884 if (!pctx->haveScale) {
885 rpos[3] = 0;
886 }
887 } else {
888 ELL_3V_SET(rpos, airDrandMT_r(rng), airDrandMT_r(rng), airDrandMT_r(rng))((rpos)[0] = (airDrandMT_r(rng)), (rpos)[1] = (airDrandMT_r(rng
)), (rpos)[2] = (airDrandMT_r(rng)));
889 if (pctx->haveScale) {
890 rpos[3] = airDrandMT_r(rng);
891 } else {
892 rpos[3] = 0;
893 }
894 }
895 _pullUnitToWorld(pctx, scaleVol, point->pos, rpos);
896 if (pctx->flag.zeroZ) {
897 point->pos[2] = 0.0;
898 }
899 /*
900 verbo = (AIR_ABS(-0.246015 - point->pos[0]) < 0.1 &&
901 AIR_ABS(-144.78 - point->pos[0]) < 0.1 &&
902 AIR_ABS(-85.3813 - point->pos[0]) < 0.1);
903 */
904 verbo = AIR_FALSE0;
905 if (verbo) {
906 fprintf(stderr__stderrp, "%s: verbo on for point %u at %g %g %g %g\n", me,
907 point->idtag, point->pos[0], point->pos[1],
908 point->pos[2], point->pos[3]);
909 }
910 _pullPointHistAdd(point, pullCondOld, AIR_NAN);
911
912 /* HEY copy and paste */
913 if (pctx->ispec[pullInfoSeedPreThresh]) {
914 /* we first do a special-purpose probe just for SeedPreThresh */
915 pctx->task[0]->probeSeedPreThreshOnly = AIR_TRUE1;
916 if (pullProbe(pctx->task[0], point)) {
917 biffAddf(PULLpullBiffKey, "%s: pre-probing pointIdx %u of world", me, pointIdx);
918 return 1;
919 }
920 pctx->task[0]->probeSeedPreThreshOnly = AIR_FALSE0;
921 if (_threshFail(pctx, point, pullInfoSeedPreThresh)) {
922 threshFailCount++;
923 spthreshFailCount++;
924 reject = AIR_TRUE1;
925 goto reckoning;
926 }
927 }
928 if (pullProbe(pctx->task[0], point)) {
929 biffAddf(PULLpullBiffKey, "%s: probing pointIdx %u of world", me, pointIdx);
930 return 1;
931 }
932 /* Check we pass pre-threshold */
933#define THRESH_TEST(INFO)if (pctx->ispec[INFO] && _threshFail(pctx, point, INFO
)) { threshFailCount++; reject = 1; goto reckoning; } \
934 if (pctx->ispec[INFO] && _threshFail(pctx, point, INFO)) { \
935 threshFailCount++; \
936 reject = AIR_TRUE1; \
937 goto reckoning; \
938 }
939 /* fprintf(stderr, "!%s: bi ngo 0 (%d) %d %p\n", me,
940 !pctx->flag.constraintBeforeSeedThresh,
941 pctx->initParm.liveThreshUse, pctx->ispec[pullInfoLiveThresh]); */
942 if (!pctx->flag.constraintBeforeSeedThresh) {
943 THRESH_TEST(pullInfoSeedThresh)if (pctx->ispec[pullInfoSeedThresh] && _threshFail
(pctx, point, pullInfoSeedThresh)) { threshFailCount++; reject
= 1; goto reckoning; };
944 if (pctx->initParm.liveThreshUse) {
945 THRESH_TEST(pullInfoLiveThresh)if (pctx->ispec[pullInfoLiveThresh] && _threshFail
(pctx, point, pullInfoLiveThresh)) { threshFailCount++; reject
= 1; goto reckoning; };
946 THRESH_TEST(pullInfoLiveThresh2)if (pctx->ispec[pullInfoLiveThresh2] && _threshFail
(pctx, point, pullInfoLiveThresh2)) { threshFailCount++; reject
= 1; goto reckoning; };
947 THRESH_TEST(pullInfoLiveThresh3)if (pctx->ispec[pullInfoLiveThresh3] && _threshFail
(pctx, point, pullInfoLiveThresh3)) { threshFailCount++; reject
= 1; goto reckoning; };
948 }
949 }
950 /* fprintf(stderr, "!%s: bi ngo 1\n", me); */
951
952 if (pctx->constraint) {
953 int constrFail;
954 /* fprintf(stderr, "!%s: calling _pullConstraintSatisfy(%u)\n", me, pointIdx); */
955 if (_pullConstraintSatisfy(pctx->task[0], point,
956 _PULL_CONSTRAINT_TRAVEL_MAX2,
957 &constrFail)) {
958 biffAddf(PULLpullBiffKey, "%s: trying constraint on point %u", me, pointIdx);
959 return 1;
960 }
961#if PULL_PHIST0
962 if (DEBUG1) {
963 Nrrd *nhist;
964 FILE *fhist;
965 char fname[AIR_STRLEN_SMALL(128+1)];
966 nhist = nrrdNew();
967 if (pullPositionHistoryNrrdGet(nhist, pctx, point)) {
968 biffAddf(PULLpullBiffKey, "%s: trouble", me);
969 return 1;
970 }
971 sprintf(fname, "%04u-%04u-%04u-phist.nrrd", pctx->iter,__builtin___sprintf_chk (fname, 0, __builtin_object_size (fname
, 2 > 1 ? 1 : 0), "%04u-%04u-%04u-phist.nrrd", pctx->iter
, point->idtag, tryCount)
972 point->idtag, tryCount)__builtin___sprintf_chk (fname, 0, __builtin_object_size (fname
, 2 > 1 ? 1 : 0), "%04u-%04u-%04u-phist.nrrd", pctx->iter
, point->idtag, tryCount);
973 if ((fhist = fopen(fname, "w"))) {
974 if (nrrdSave(fname, nhist, NULL((void*)0))) {
975 biffMovef(PULLpullBiffKey, NRRDnrrdBiffKey, "%s: trouble", me);
976 return 1;
977 }
978 fclose(fhist);
979 }
980 nrrdNuke(nhist);
981 }
982#endif
983 if (constrFail) {
984 constrFailCount++;
985 reject = AIR_TRUE1;
986 goto reckoning;
987 }
988 /* fprintf(stderr, "!%s: bi ngo 2\n", me); */
989 /* post constraint-satisfaction, we certainly have to assert thresholds */
990 THRESH_TEST(pullInfoSeedThresh)if (pctx->ispec[pullInfoSeedThresh] && _threshFail
(pctx, point, pullInfoSeedThresh)) { threshFailCount++; reject
= 1; goto reckoning; };
991 if (pctx->initParm.liveThreshUse) {
992 THRESH_TEST(pullInfoLiveThresh)if (pctx->ispec[pullInfoLiveThresh] && _threshFail
(pctx, point, pullInfoLiveThresh)) { threshFailCount++; reject
= 1; goto reckoning; };
993 THRESH_TEST(pullInfoLiveThresh2)if (pctx->ispec[pullInfoLiveThresh2] && _threshFail
(pctx, point, pullInfoLiveThresh2)) { threshFailCount++; reject
= 1; goto reckoning; };
994 THRESH_TEST(pullInfoLiveThresh3)if (pctx->ispec[pullInfoLiveThresh3] && _threshFail
(pctx, point, pullInfoLiveThresh3)) { threshFailCount++; reject
= 1; goto reckoning; };
995 }
996 /* fprintf(stderr, "!%s: bi ngo 3 (reject=%d)\n", me, reject); */
997 }
998
999 reckoning:
1000 if (reject) {
1001 if (threshFailCount + constrFailCount >= _PULL_RANDOM_SEED_TRY_MAX1000000) {
1002 /* Very bad luck; we've too many times */
1003 biffAddf(PULLpullBiffKey, "%s: failed too often (%u times) placing point %u: "
1004 "%u fails on thresh (%u on pre-thresh), %u on constr",
1005 me, _PULL_RANDOM_SEED_TRY_MAX1000000, pointIdx,
1006 threshFailCount, spthreshFailCount, constrFailCount);
1007 return 1;
1008 }
1009 }
1010 } while (reject);
1011
1012 if (pullInitMethodHalton == pctx->initParm.method) {
1013 pctx->haltonOffset += threshFailCount + constrFailCount;
1014 }
1015
1016 return 0;
1017}
1018
1019int
1020pullPointInitializeGivenPos(pullContext *pctx,
1021 const double *posData,
1022 const unsigned int pointIdx,
1023 pullPoint *point,
1024 /* output */
1025 int *createFailP) {
1026 static const char me[]="pullPointInitializeGivenPos";
1027 int reject, rejectEdge;
1028
1029 /* Copy nrrd point into pullPoint */
1030 ELL_4V_COPY(point->pos, posData + 4*pointIdx)((point->pos)[0] = (posData + 4*pointIdx)[0], (point->pos
)[1] = (posData + 4*pointIdx)[1], (point->pos)[2] = (posData
+ 4*pointIdx)[2], (point->pos)[3] = (posData + 4*pointIdx
)[3]);
1031 if (pctx->flag.zeroZ) {
1032 point->pos[2] = 0.0;
1033 }
1034
1035 /*
1036 if (AIR_ABS(247.828 - point->pos[0]) < 0.1 &&
1037 AIR_ABS(66.8817 - point->pos[1]) < 0.1 &&
1038 AIR_ABS(67.0031 - point->pos[2]) < 0.1) {
1039 fprintf(stderr, "%s: --------- point %u at %g %g %g %g\n", me,
1040 point->idtag,
1041 point->pos[0], point->pos[1],
1042 point->pos[2], point->pos[3]);
1043 }
1044 */
1045
1046 /* we're dictating positions, but still have to do initial probe,
1047 and possibly liveThresholding */
1048 if (pullProbe(pctx->task[0], point)) {
1049 biffAddf(PULLpullBiffKey, "%s: probing pointIdx %u of npos", me, pointIdx);
1050 return 1;
1051 }
1052 reject = AIR_FALSE0;
1053 if (pctx->flag.nixAtVolumeEdgeSpace
1054 && (point->status & PULL_STATUS_EDGE_BIT(1<< 4))) {
1055 rejectEdge = AIR_TRUE1;
1056 } else {
1057 rejectEdge = AIR_FALSE0;
1058 }
1059 reject |= rejectEdge;
1060 if (pctx->initParm.liveThreshUse) {
1061 if (!reject) reject |= _threshFail(pctx, point, pullInfoLiveThresh);
1062 if (!reject) reject |= _threshFail(pctx, point, pullInfoLiveThresh2);
1063 if (!reject) reject |= _threshFail(pctx, point, pullInfoLiveThresh3);
1064 }
1065 /*
1066 if (reject) {
1067 fprintf(stderr, "!%s(%u): edge %d thresh %d %d %d\n",
1068 me, point->idtag, rejectEdge,
1069 _threshFail(pctx, point, pullInfoLiveThresh),
1070 _threshFail(pctx, point, pullInfoLiveThresh2),
1071 _threshFail(pctx, point, pullInfoLiveThresh3));
1072 }
1073 */
1074 if (reject) {
1075 *createFailP = AIR_TRUE1;
1076 } else {
1077 *createFailP = AIR_FALSE0;
1078 }
1079
1080 return 0;
1081}
1082
1083/*
1084** _pullPointSetup sets:
1085**
1086** This is only called by the master thread
1087**
1088** this should set stuff to be like after an update stage and
1089** just before the rebinning
1090*/
1091int
1092_pullPointSetup(pullContext *pctx) {
1093 static const char me[]="_pullPointSetup";
1094 char doneStr[AIR_STRLEN_SMALL(128+1)];
1095 unsigned int pointIdx, binIdx, tick, pn;
1096 pullPoint *point;
1097 pullBin *bin;
1098 int createFail,added;
1099 airArray *mop;
1100 Nrrd *npos;
1101 pullVolume *seedVol, *scaleVol;
1102 gageShape *seedShape;
1103 double factor, *posData;
1104 unsigned int totalNumPoints, voxNum, ii;
1105
1106 /* on using pullBinsPointMaybeAdd: This is used only in the context
1107 of constraints; it makes sense to be a little more selective
1108 about adding points, so that they aren't completely piled on top
1109 of each other. This relies on _PULL_BINNING_MAYBE_ADD_THRESH: its
1110 tempting to set this value high, to more aggressively limit the
1111 number of points added, but that's really the job of population
1112 control, and we can't always guarantee that constraint manifolds
1113 will be well-sampled (with respect to pctx->radiusSpace and
1114 pctx->radiusScale) to start with */
1115
1116 if (pctx->verbose) {
1117 printf("%s: beginning (%s). . . ", me,
1118 airEnumStr(pullInitMethod, pctx->initParm.method));
1119 fflush(stdout__stdoutp);
1120 }
1121 mop = airMopNew();
1122 switch (pctx->initParm.method) {
1123 case pullInitMethodGivenPos:
1124 npos = nrrdNew();
1125 airMopAdd(mop, npos, (airMopper)nrrdNuke, airMopAlways);
1126 /* even if npos came in as double, we have to copy it */
1127 if (nrrdConvert(npos, pctx->initParm.npos, nrrdTypeDouble)) {
1128 biffMovef(PULLpullBiffKey, NRRDnrrdBiffKey, "%s: trouble converting npos", me);
1129 airMopError(mop); return 1;
1130 }
1131 posData = AIR_CAST(double *, npos->data)((double *)(npos->data));
1132 if (pctx->initParm.numInitial || pctx->initParm.pointPerVoxel) {
1133 printf("%s: with npos, overriding both numInitial (%u) "
1134 "and pointPerVoxel (%d)\n", me, pctx->initParm.numInitial,
1135 pctx->initParm.pointPerVoxel);
1136 }
1137 totalNumPoints = AIR_CAST(unsigned int, npos->axis[1].size)((unsigned int)(npos->axis[1].size));
1138 break;
1139 case pullInitMethodPointPerVoxel:
1140 npos = NULL((void*)0);
1141 posData = NULL((void*)0);
1142 if (pctx->initParm.numInitial && pctx->verbose) {
1143 printf("%s: pointPerVoxel %d overrides numInitial (%u)\n", me,
1144 pctx->initParm.pointPerVoxel, pctx->initParm.numInitial);
1145 }
1146 /* Obtain number of voxels */
1147 seedVol = pctx->vol[pctx->ispec[pullInfoSeedThresh]->volIdx];
1148 seedShape = seedVol->gctx->shape;
1149 if (pctx->initParm.ppvZRange[0] <= pctx->initParm.ppvZRange[1]) {
1150 unsigned int zrn;
1151 if (!( pctx->initParm.ppvZRange[0] < seedShape->size[2]
1152 && pctx->initParm.ppvZRange[1] < seedShape->size[2] )) {
1153 biffAddf(PULLpullBiffKey, "%s: ppvZRange[%u,%u] outside volume [0,%u]", me,
1154 pctx->initParm.ppvZRange[0], pctx->initParm.ppvZRange[1],
1155 seedShape->size[2]-1);
1156 airMopError(mop); return 1;
1157 }
1158 zrn = pctx->initParm.ppvZRange[1] - pctx->initParm.ppvZRange[0] + 1;
1159 voxNum = seedShape->size[0]*seedShape->size[1]*zrn;
1160 if (pctx->verbose) {
1161 printf("%s: vol size %u %u [%u,%u] -> voxNum %u\n", me,
1162 seedShape->size[0], seedShape->size[1],
1163 pctx->initParm.ppvZRange[0], pctx->initParm.ppvZRange[1],
1164 voxNum);
1165 }
1166 } else {
1167 voxNum = seedShape->size[0]*seedShape->size[1]*seedShape->size[2];
1168 if (pctx->verbose) {
1169 printf("%s: vol size %u %u %u -> voxNum %u\n", me,
1170 seedShape->size[0], seedShape->size[1], seedShape->size[2],
1171 voxNum);
1172 }
1173 }
1174
1175 /* Compute total number of points */
1176 if (pctx->initParm.pointPerVoxel > 0) {
1177 factor = pctx->initParm.pointPerVoxel;
1178 } else {
1179 factor = -1.0/pctx->initParm.pointPerVoxel;
1180 }
1181 if (pctx->haveScale) {
1182 unsigned int sasn;
1183 sasn = pctx->initParm.samplesAlongScaleNum;
1184 totalNumPoints = AIR_CAST(unsigned int, voxNum * factor * sasn)((unsigned int)(voxNum * factor * sasn));
1185 } else {
1186 totalNumPoints = AIR_CAST(unsigned int, voxNum * factor)((unsigned int)(voxNum * factor));
1187 }
1188 break;
1189 case pullInitMethodRandom:
1190 case pullInitMethodHalton:
1191 npos = NULL((void*)0);
1192 posData = NULL((void*)0);
1193 totalNumPoints = pctx->initParm.numInitial;
1194 break;
1195 default:
1196 biffAddf(PULLpullBiffKey, "%s: pullInitMethod %d not handled!", me,
1197 pctx->initParm.method);
1198 airMopError(mop); return 1;
1199 break;
1200 }
1201 if (pctx->verbose) {
1202 printf("%s: initializing/seeding %u pts. . . ", me, totalNumPoints);
1203 fflush(stdout__stdoutp);
1204 }
1205
1206 /* find first scale volume, if there is one; this is used by some
1207 seeders to determine placement along the scale axis */
1208 scaleVol = NULL((void*)0);
1209 for (ii=0; ii<pctx->volNum; ii++) {
1210 if (pctx->vol[ii]->ninScale) {
1211 scaleVol = pctx->vol[ii];
1212 break;
1213 }
1214 }
1215
1216 /* Start adding points */
1217 tick = totalNumPoints/1000;
1218 point = NULL((void*)0);
1219 unsigned int initRorHack = 0;
1220 /* This loop would normally be:
1221 for (pointIdx = 0; pointIdx < totalNumPoints; pointIdx++) {
1222 but because of pctx->flag.nixAtVolumeEdgeSpaceInitRorH we
1223 need to be able to decrement pointIdx to force another redo,
1224 even if pointIdx==0. So loop index is actually one greater
1225 than the index value actually used */
1226 for (pointIdx = 1; pointIdx <= totalNumPoints; pointIdx++) {
1227 int E;
1228 if (pctx->verbose) {
1229 if (tick < 100 || 0 == (pointIdx-1) % tick) {
1230 printf("%s", airDoneStr(0, pointIdx-1, totalNumPoints, doneStr));
1231 fflush(stdout__stdoutp);
1232 }
1233 }
1234 if (pctx->verbose > 5) {
1235 printf("\n%s: setting up point = %u/%u\n", me,
1236 pointIdx-1, totalNumPoints);
1237 }
1238 /* Create point */
1239 if (!point) {
1240 point = pullPointNew(pctx);
1241 }
1242 /* Filling array according to initialization method */
1243 E = 0;
1244 switch(pctx->initParm.method) {
1245 case pullInitMethodRandom:
1246 case pullInitMethodHalton:
1247 /* point index passed here always increments, even if
1248 we failed last time because of being at edge */
1249 E = pullPointInitializeRandomOrHalton(pctx, pointIdx-1 + initRorHack,
1250 point, scaleVol);
1251 if (pctx->flag.nixAtVolumeEdgeSpaceInitRorH) {
1252 createFail = !!(point->status & PULL_STATUS_EDGE_BIT(1<< 4));
1253 if (createFail) {
1254 initRorHack++;
1255 if (initRorHack > 10*totalNumPoints) {
1256 biffAddf(PULLpullBiffKey, "%s: (point %u) handling nixAtVolumeEdgeSpaceInitRorH, "
1257 "have failed so often (%u > 10*#pnts = %u); something "
1258 "must be wrong", me, point->idtag, initRorHack,
1259 totalNumPoints);
1260 airMopError(mop); return 1;
1261 }
1262 /* this is a for-loop, post-body increment will always happen;
1263 we need to subvert it. */
1264 pointIdx--;
1265 }
1266 } else {
1267 createFail = AIR_FALSE0;
1268 }
1269 break;
1270 case pullInitMethodPointPerVoxel:
1271 E = pullPointInitializePerVoxel(pctx, pointIdx-1, point, scaleVol,
1272 &createFail);
1273 break;
1274 case pullInitMethodGivenPos:
1275 E = pullPointInitializeGivenPos(pctx, posData, pointIdx-1, point,
1276 &createFail);
1277 break;
1278 }
1279 if (E) {
1280 biffAddf(PULLpullBiffKey, "%s: trouble trying point %u (id %u)", me,
1281 pointIdx-1, point->idtag);
1282 airMopError(mop); return 1;
1283 }
1284 if (createFail) {
1285 /* We were not successful in creating a point; not an error */
1286 continue;
1287 }
1288
1289 /* else, the point is ready for binning */
1290 if (pctx->constraint) {
1291 if (pullBinsPointMaybeAdd(pctx, point, NULL((void*)0), &added)) {
1292 biffAddf(PULLpullBiffKey, "%s: trouble binning point %u", me, point->idtag);
1293 airMopError(mop); return 1;
1294 }
1295 /*
1296 if (4523 == point->idtag) {
1297 fprintf(stderr, "!%s(%u): ----- added=%d at %g %g %g %g\n",
1298 me, point->idtag, added,
1299 point->pos[0], point->pos[1], point->pos[2], point->pos[3]);
1300 }
1301 */
1302 if (added) {
1303 point = NULL((void*)0);
1304 }
1305 } else {
1306 if (pullBinsPointAdd(pctx, point, NULL((void*)0))) {
1307 biffAddf(PULLpullBiffKey, "%s: trouble binning point %u", me, point->idtag);
1308 airMopError(mop); return 1;
1309 }
1310 point = NULL((void*)0);
1311 }
1312 } /* Done looping through total number of points */
1313 if (pctx->verbose) {
1314 printf("%s\n", airDoneStr(0, pointIdx-1, totalNumPoints,
1315 doneStr));
1316 }
1317 if (point) {
1318 /* we created a new test point, but it was never placed in the volume */
1319 /* so, HACK: undo pullPointNew . . . */
1320 point = pullPointNix(point);
Value stored to 'point' is never read
1321 /* pctx->idtagNext -= 1; */
1322 }
1323
1324 /* Final check: do we have any points? */
1325 pn = pullPointNumber(pctx);
1326 if (!pn) {
1327 char stmp1[AIR_STRLEN_MED(256+1)], stmp2[AIR_STRLEN_MED(256+1)];
1328 int guess=AIR_FALSE0;
1329 sprintf(stmp1, "%s: seeded 0 points", me)__builtin___sprintf_chk (stmp1, 0, __builtin_object_size (stmp1
, 2 > 1 ? 1 : 0), "%s: seeded 0 points", me);
1330 if (pctx->ispec[pullInfoSeedThresh]) {
1331 guess=AIR_TRUE1;
1332 sprintf(stmp2, " (? bad seedthresh %g ?)",__builtin___sprintf_chk (stmp2, 0, __builtin_object_size (stmp2
, 2 > 1 ? 1 : 0), " (? bad seedthresh %g ?)", pctx->ispec
[pullInfoSeedThresh]->zero)
1333 pctx->ispec[pullInfoSeedThresh]->zero)__builtin___sprintf_chk (stmp2, 0, __builtin_object_size (stmp2
, 2 > 1 ? 1 : 0), " (? bad seedthresh %g ?)", pctx->ispec
[pullInfoSeedThresh]->zero);
1334 strcat(stmp1, stmp2)__builtin___strcat_chk (stmp1, stmp2, __builtin_object_size (
stmp1, 2 > 1 ? 1 : 0));
1335 }
1336 if (pctx->flag.nixAtVolumeEdgeSpace) {
1337 guess=AIR_TRUE1;
1338 sprintf(stmp2, " (? flag.nixAtVolumeEdgeSpace true ?)")__builtin___sprintf_chk (stmp2, 0, __builtin_object_size (stmp2
, 2 > 1 ? 1 : 0), " (? flag.nixAtVolumeEdgeSpace true ?)");
1339 strcat(stmp1, stmp2)__builtin___strcat_chk (stmp1, stmp2, __builtin_object_size (
stmp1, 2 > 1 ? 1 : 0));
1340 }
1341 if (!guess) {
1342 sprintf(stmp2, " (no guess as to why)")__builtin___sprintf_chk (stmp2, 0, __builtin_object_size (stmp2
, 2 > 1 ? 1 : 0), " (no guess as to why)");
1343 strcat(stmp1, stmp2)__builtin___strcat_chk (stmp1, stmp2, __builtin_object_size (
stmp1, 2 > 1 ? 1 : 0));
1344 }
1345 biffAddf(PULLpullBiffKey, "%s", stmp1);
1346 airMopError(mop); return 1;
1347 }
1348 if (pctx->verbose) {
1349 fprintf(stderr__stderrp, "%s: initialized to %u points (idtagNext = %u)\n",
1350 me, pn, pctx->idtagNext);
1351 }
1352 /* */
1353 if (1) {
1354 Nrrd *ntmp;
1355 ntmp = nrrdNew();
1356 pullOutputGet(ntmp, NULL((void*)0), NULL((void*)0), NULL((void*)0), 0.0, pctx);
1357 nrrdSave("pos-init.nrrd", ntmp, NULL((void*)0));
1358 nrrdNuke(ntmp);
1359 }
1360 /* */
1361 pctx->tmpPointPtr = AIR_CAST(pullPoint **,((pullPoint **)(calloc(pn, sizeof(pullPoint*))))
1362 calloc(pn, sizeof(pullPoint*)))((pullPoint **)(calloc(pn, sizeof(pullPoint*))));
1363 pctx->tmpPointPerm = AIR_CAST(unsigned int *,((unsigned int *)(calloc(pn, sizeof(unsigned int))))
1364 calloc(pn, sizeof(unsigned int)))((unsigned int *)(calloc(pn, sizeof(unsigned int))));
1365 if (!( pctx->tmpPointPtr && pctx->tmpPointPerm )) {
1366 biffAddf(PULLpullBiffKey, "%s: couldn't allocate tmp buffers %p %p", me,
1367 AIR_VOIDP(pctx->tmpPointPtr)((void *)(pctx->tmpPointPtr)), AIR_VOIDP(pctx->tmpPointPerm)((void *)(pctx->tmpPointPerm)));
1368 airMopError(mop); return 1;
1369 }
1370 pctx->tmpPointNum = pn;
1371
1372 /* now that all points have been added, set their energy to help
1373 inspect initial state */
1374 for (binIdx=0; binIdx<pctx->binNum; binIdx++) {
1375 bin = pctx->bin + binIdx;
1376 for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) {
1377 point = bin->point[pointIdx];
1378 point->energy = _pullPointEnergyTotal(pctx->task[0], bin, point,
1379 /* ignoreImage */
1380 !pctx->haveScale,
1381 point->force);
1382 }
1383 }
1384
1385 if (pctx->verbose) {
1386 printf("%s: all done. ", me);
1387 fflush(stdout__stdoutp);
1388 }
1389 airMopOkay(mop);
1390 return 0;
1391}
1392
1393void
1394_pullPointFinish(pullContext *pctx) {
1395
1396 airFree(pctx->tmpPointPtr);
1397 airFree(pctx->tmpPointPerm);
1398 return ;
1399}