1 /*
2 Teem: Tools to process and visualize scientific data and images .
3 Copyright (C) 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 "push.h"
25 #include "privatePush.h"
26
27 unsigned int
_pushPointTotal(pushContext * pctx)28 _pushPointTotal(pushContext *pctx) {
29 unsigned int binIdx, pointNum;
30 pushBin *bin;
31
32 pointNum = 0;
33 for (binIdx=0; binIdx<pctx->binNum; binIdx++) {
34 bin = pctx->bin + binIdx;
35 pointNum += bin->pointNum;
36 }
37 return pointNum;
38 }
39
40 int
_pushProbe(pushTask * task,pushPoint * point)41 _pushProbe(pushTask *task, pushPoint *point) {
42 static const char me[]="_pushProbe";
43 double posWorld[4], posIdx[4];
44
45 ELL_3V_COPY(posWorld, point->pos); posWorld[3] = 1.0;
46 ELL_4MV_MUL(posIdx, task->gctx->shape->WtoI, posWorld);
47 ELL_4V_HOMOG(posIdx, posIdx);
48 posIdx[0] = AIR_CLAMP(-0.5, posIdx[0], task->gctx->shape->size[0]-0.5);
49 posIdx[1] = AIR_CLAMP(-0.5, posIdx[1], task->gctx->shape->size[1]-0.5);
50 posIdx[2] = AIR_CLAMP(-0.5, posIdx[2], task->gctx->shape->size[2]-0.5);
51
52 if (gageProbe(task->gctx, posIdx[0], posIdx[1], posIdx[2])) {
53 biffAddf(PUSH, "%s: gageProbe failed:\n (%d) %s\n", me,
54 task->gctx->errNum, task->gctx->errStr);
55 return 1;
56 }
57
58 TEN_T_COPY(point->ten, task->tenAns);
59 TEN_T_COPY(point->inv, task->invAns);
60 ELL_3V_COPY(point->cnt, task->cntAns);
61 if (tenGageUnknown != task->pctx->gravItem) {
62 point->grav = task->gravAns[0];
63 ELL_3V_COPY(point->gravGrad, task->gravGradAns);
64 }
65 if (tenGageUnknown != task->pctx->seedThreshItem) {
66 point->seedThresh = task->seedThreshAns[0];
67 }
68 return 0;
69 }
70
71 int
pushOutputGet(Nrrd * nPosOut,Nrrd * nTenOut,Nrrd * nEnrOut,pushContext * pctx)72 pushOutputGet(Nrrd *nPosOut, Nrrd *nTenOut, Nrrd *nEnrOut,
73 pushContext *pctx) {
74 static const char me[]="pushOutputGet";
75 unsigned int binIdx, pointRun, pointNum, pointIdx;
76 int E;
77 float *posOut, *tenOut, *enrOut;
78 pushBin *bin;
79 pushPoint *point;
80
81 pointNum = _pushPointTotal(pctx);
82 E = AIR_FALSE;
83 if (nPosOut) {
84 E |= nrrdMaybeAlloc_va(nPosOut, nrrdTypeFloat, 2,
85 AIR_CAST(size_t, 3),
86 AIR_CAST(size_t, pointNum));
87 }
88 if (nTenOut) {
89 E |= nrrdMaybeAlloc_va(nTenOut, nrrdTypeFloat, 2,
90 AIR_CAST(size_t, 7),
91 AIR_CAST(size_t, pointNum));
92 }
93 if (nEnrOut) {
94 E |= nrrdMaybeAlloc_va(nEnrOut, nrrdTypeFloat, 1,
95 AIR_CAST(size_t, pointNum));
96 }
97 if (E) {
98 biffMovef(PUSH, NRRD, "%s: trouble allocating outputs", me);
99 return 1;
100 }
101 posOut = nPosOut ? (float*)(nPosOut->data) : NULL;
102 tenOut = nTenOut ? (float*)(nTenOut->data) : NULL;
103 enrOut = nEnrOut ? (float*)(nEnrOut->data) : NULL;
104
105 pointRun = 0;
106 for (binIdx=0; binIdx<pctx->binNum; binIdx++) {
107 bin = pctx->bin + binIdx;
108 for (pointIdx=0; pointIdx<bin->pointNum; pointIdx++) {
109 point = bin->point[pointIdx];
110 if (posOut) {
111 ELL_3V_SET_TT(posOut + 3*pointRun, float,
112 point->pos[0], point->pos[1], point->pos[2]);
113 }
114 if (tenOut) {
115 TEN_T_COPY_TT(tenOut + 7*pointRun, float, point->ten);
116 }
117 if (enrOut) {
118 enrOut[pointRun] = AIR_CAST(float, point->enr);
119 }
120 pointRun++;
121 }
122 }
123
124 return 0;
125 }
126
127 int
_pushPairwiseEnergy(pushTask * task,double * enrP,double frc[3],pushEnergySpec * ensp,pushPoint * myPoint,pushPoint * herPoint,double YY[3],double iscl)128 _pushPairwiseEnergy(pushTask *task,
129 double *enrP,
130 double frc[3],
131 pushEnergySpec *ensp,
132 pushPoint *myPoint, pushPoint *herPoint,
133 double YY[3], double iscl) {
134 static const char me[]="_pushPairwiseEnergy";
135 double inv[7], XX[3], nXX[3], rr, mag, WW[3];
136
137 if (task->pctx->midPntSmp) {
138 pushPoint _tmpPoint;
139 double det;
140 ELL_3V_SCALE_ADD2(_tmpPoint.pos,
141 0.5, myPoint->pos,
142 0.5, herPoint->pos);
143 if (_pushProbe(task, &_tmpPoint)) {
144 biffAddf(PUSH, "%s: at midpoint of %u and %u", me,
145 myPoint->ttaagg, herPoint->ttaagg);
146 *enrP = AIR_NAN; return 1;
147 }
148 TEN_T_INV(inv, _tmpPoint.ten, det);
149 } else {
150 TEN_T_SCALE_ADD2(inv,
151 0.5, myPoint->inv,
152 0.5, herPoint->inv);
153 }
154 TEN_TV_MUL(XX, inv, YY);
155 ELL_3V_NORM(nXX, XX, rr);
156
157 ensp->energy->eval(enrP, &mag, rr*iscl, ensp->parm);
158 if (mag) {
159 mag *= iscl;
160 TEN_TV_MUL(WW, inv, nXX);
161 ELL_3V_SCALE(frc, mag, WW);
162 } else {
163 ELL_3V_SET(frc, 0, 0, 0);
164 }
165
166 return 0;
167 }
168
169 #define EPS_PER_MAX_DIST 200
170 #define SEEK_MAX_ITER 30
171
172 int
pushBinProcess(pushTask * task,unsigned int myBinIdx)173 pushBinProcess(pushTask *task, unsigned int myBinIdx) {
174 static const char me[]="pushBinProcess";
175 pushBin *myBin, *herBin, **neighbor;
176 unsigned int myPointIdx, herPointIdx;
177 pushPoint *myPoint, *herPoint;
178 double enr, frc[3], delta[3], deltaLen, deltaNorm[3], warp[3], limit,
179 maxDiffLenSqrd, iscl, diff[3], diffLenSqrd;
180
181 if (task->pctx->verbose > 2) {
182 fprintf(stderr, "%s(%u): doing bin %u\n", me, task->threadIdx, myBinIdx);
183 }
184 maxDiffLenSqrd = (task->pctx->maxDist)*(task->pctx->maxDist);
185 myBin = task->pctx->bin + myBinIdx;
186 iscl = 1.0/(2*task->pctx->scale);
187 for (myPointIdx=0; myPointIdx<myBin->pointNum; myPointIdx++) {
188 myPoint = myBin->point[myPointIdx];
189 myPoint->enr = 0;
190 ELL_3V_SET(myPoint->frc, 0, 0, 0);
191
192 if (1.0 <= task->pctx->neighborTrueProb
193 || airDrandMT_r(task->rng) <= task->pctx->neighborTrueProb
194 || !myPoint->neighArr->len) {
195 neighbor = myBin->neighbor;
196 if (1.0 > task->pctx->neighborTrueProb) {
197 airArrayLenSet(myPoint->neighArr, 0);
198 }
199 while ((herBin = *neighbor)) {
200 for (herPointIdx=0; herPointIdx<herBin->pointNum; herPointIdx++) {
201 herPoint = herBin->point[herPointIdx];
202 if (myPoint == herPoint) {
203 /* can't interact with myself */
204 continue;
205 }
206 ELL_3V_SUB(diff, herPoint->pos, myPoint->pos);
207 diffLenSqrd = ELL_3V_DOT(diff, diff);
208 if (diffLenSqrd > maxDiffLenSqrd) {
209 /* too far away to interact */
210 continue;
211 }
212 if (_pushPairwiseEnergy(task, &enr, frc, task->pctx->ensp,
213 myPoint, herPoint, diff, iscl)) {
214 biffAddf(PUSH, "%s: between points %u and %u, A", me,
215 myPoint->ttaagg, herPoint->ttaagg);
216 return 1;
217 }
218 myPoint->enr += enr/2;
219 if (ELL_3V_DOT(frc, frc)) {
220 ELL_3V_INCR(myPoint->frc, frc);
221 if (1.0 > task->pctx->neighborTrueProb) {
222 unsigned int idx;
223 idx = airArrayLenIncr(myPoint->neighArr, 1);
224 myPoint->neigh[idx] = herPoint;
225 }
226 }
227 if (!ELL_3V_EXISTS(myPoint->frc)) {
228 biffAddf(PUSH, "%s: bad myPoint->frc (%g,%g,%g) @ bin %p end", me,
229 myPoint->frc[0], myPoint->frc[1], myPoint->frc[2],
230 AIR_VOIDP(herBin));
231 return 1;
232 }
233 }
234 neighbor++;
235 }
236 } else {
237 /* we are doing neighborhood list optimization, and this is an
238 iteration where we use the list. So the body of this loop
239 has to be the same as the meat of the above loop */
240 unsigned int neighIdx;
241 for (neighIdx=0; neighIdx<myPoint->neighArr->len; neighIdx++) {
242 herPoint = myPoint->neigh[neighIdx];
243 ELL_3V_SUB(diff, herPoint->pos, myPoint->pos);
244 if (_pushPairwiseEnergy(task, &enr, frc, task->pctx->ensp,
245 myPoint, herPoint, diff, iscl)) {
246 biffAddf(PUSH, "%s: between points %u and %u, B", me,
247 myPoint->ttaagg, herPoint->ttaagg);
248 return 1;
249 }
250 myPoint->enr += enr/2;
251 ELL_3V_INCR(myPoint->frc, frc);
252 }
253 }
254 if (!ELL_3V_EXISTS(myPoint->frc)) {
255 biffAddf(PUSH, "%s: post-nei myPoint->frc (%g,%g,%g) doesn't exist", me,
256 myPoint->frc[0], myPoint->frc[1], myPoint->frc[2]);
257 return 1;
258 }
259
260 /* each point sees containment forces */
261 ELL_3V_SCALE(frc, task->pctx->cntScl, myPoint->cnt);
262 ELL_3V_INCR(myPoint->frc, frc);
263 myPoint->enr += task->pctx->cntScl*(1 - myPoint->ten[0]);
264
265 /* each point also maybe experiences gravity */
266 if (tenGageUnknown != task->pctx->gravItem) {
267 ELL_3V_SCALE(frc, -task->pctx->gravScl, myPoint->gravGrad);
268 myPoint->enr +=
269 task->pctx->gravScl*(myPoint->grav - task->pctx->gravZero);
270 ELL_3V_INCR(myPoint->frc, frc);
271 }
272 if (!ELL_3V_EXISTS(myPoint->frc)) {
273 biffAddf(PUSH, "%s: post-grav myPoint->frc (%g,%g,%g) doesn't exist", me,
274 myPoint->frc[0], myPoint->frc[1], myPoint->frc[2]);
275 return 1;
276 }
277
278 /* each point in this thing also maybe experiences wall forces */
279 if (task->pctx->wall) {
280 /* there's an effort here to get the forces and energies, which
281 are actually computed in index space, to be correctly scaled
282 into world space, but no promises that its right ... */
283 double enrIdx[4]={0,0,0,0}, enrWorld[4];
284 unsigned int ci;
285 double posWorld[4], posIdx[4], len, frcIdx[4], frcWorld[4];
286 ELL_3V_COPY(posWorld, myPoint->pos); posWorld[3] = 1.0;
287 ELL_4MV_MUL(posIdx, task->pctx->gctx->shape->WtoI, posWorld);
288 ELL_4V_HOMOG(posIdx, posIdx);
289 for (ci=0; ci<3; ci++) {
290 if (1 == task->pctx->gctx->shape->size[ci]) {
291 frcIdx[ci] = 0;
292 } else {
293 len = posIdx[ci] - -0.5;
294 if (len < 0) {
295 len *= -1;
296 frcIdx[ci] = task->pctx->wall*len;
297 enrIdx[ci] = task->pctx->wall*len*len/2;
298 } else {
299 len = posIdx[ci] - (task->pctx->gctx->shape->size[ci] - 0.5);
300 if (len > 0) {
301 frcIdx[ci] = -task->pctx->wall*len;
302 enrIdx[ci] = task->pctx->wall*len*len/2;
303 } else {
304 frcIdx[ci] = 0;
305 enrIdx[ci] = 0;
306 }
307 }
308 }
309 }
310 frcIdx[3] = 0.0;
311 enrIdx[3] = 0.0;
312 ELL_4MV_MUL(frcWorld, task->pctx->gctx->shape->ItoW, frcIdx);
313 ELL_4MV_MUL(enrWorld, task->pctx->gctx->shape->ItoW, enrIdx);
314 ELL_3V_INCR(myPoint->frc, frcWorld);
315 myPoint->enr += ELL_3V_LEN(enrWorld);
316 } /* wall */
317 if (!ELL_3V_EXISTS(myPoint->frc)) {
318 biffAddf(PUSH, "%s: post-wall myPoint->frc (%g,%g,%g) doesn't exist", me,
319 myPoint->frc[0], myPoint->frc[1], myPoint->frc[2]);
320 return 1;
321 }
322
323 task->energySum += myPoint->enr;
324
325 /* -------------------------------------------- */
326 /* force calculation done, now update positions */
327 /* -------------------------------------------- */
328
329 ELL_3V_SCALE(delta, task->pctx->step, myPoint->frc);
330 ELL_3V_NORM(deltaNorm, delta, deltaLen);
331 if (0 == deltaLen) {
332 /* an unforced point, but this isn't an error */
333 return 0;
334 }
335 if (!(AIR_EXISTS(deltaLen) && ELL_3V_EXISTS(deltaNorm))) {
336 biffAddf(PUSH, "%s: deltaLen %g or deltaNorm (%g,%g,%g) doesn't exist",
337 me, deltaLen, deltaNorm[0], deltaNorm[1], deltaNorm[2]);
338 return 1;
339 }
340 if (deltaLen) {
341 double newDelta;
342 TEN_TV_MUL(warp, myPoint->inv, delta);
343 /* limit is some fraction of glyph radius along direction of delta */
344 limit = (task->pctx->deltaLimit
345 *task->pctx->scale*deltaLen/(FLT_MIN + ELL_3V_LEN(warp)));
346 newDelta = limit*deltaLen/(limit + deltaLen);
347 /* by definition newDelta <= deltaLen */
348 task->deltaFracSum += newDelta/deltaLen;
349 ELL_3V_SCALE_INCR(myPoint->pos, newDelta, deltaNorm);
350 if (!ELL_3V_EXISTS(myPoint->pos)) {
351 biffAddf(PUSH, "%s: myPoint->pos %g*(%g,%g,%g) --> (%g,%g,%g) "
352 "doesn't exist", me,
353 newDelta, deltaNorm[0], deltaNorm[1], deltaNorm[2],
354 myPoint->pos[0], myPoint->pos[1], myPoint->pos[2]);
355 return 1;
356 }
357 }
358 if (2 == task->pctx->dimIn) {
359 double posIdx[4], posWorld[4], posOrig[4];
360 ELL_3V_COPY(posOrig, myPoint->pos); posOrig[3] = 1.0;
361 ELL_4MV_MUL(posIdx, task->pctx->gctx->shape->WtoI, posOrig);
362 ELL_4V_HOMOG(posIdx, posIdx);
363 posIdx[task->pctx->sliceAxis] = 0.0;
364 ELL_4MV_MUL(posWorld, task->pctx->gctx->shape->ItoW, posIdx);
365 ELL_34V_HOMOG(myPoint->pos, posWorld);
366 if (!ELL_3V_EXISTS(myPoint->pos)) {
367 biffAddf(PUSH, "%s: myPoint->pos (%g,%g,%g) -> (%g,%g,%g) "
368 "doesn't exist", me,
369 posOrig[0], posOrig[1], posOrig[2],
370 myPoint->pos[0], myPoint->pos[1], myPoint->pos[2]);
371 return 1;
372 }
373 }
374 if (1.0 <= task->pctx->probeProb
375 || airDrandMT_r(task->rng) <= task->pctx->probeProb) {
376 if (_pushProbe(task, myPoint)) {
377 biffAddf(PUSH, "%s: probing at new field pos", me);
378 return 1;
379 }
380 }
381
382 /* the point lived, count it */
383 task->pointNum += 1;
384 } /* for myPointIdx */
385
386 return 0;
387 }
388