Bug Summary

File:src/ten/grads.c
Location:line 411, column 3
Description:Value stored to 'done' 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#include "ten.h"
25#include "privateTen.h"
26
27tenGradientParm *
28tenGradientParmNew(void) {
29 tenGradientParm *ret;
30
31 ret = (tenGradientParm *)calloc(1, sizeof(tenGradientParm));
32 if (ret) {
33 ret->initStep = 1.0;
34 ret->jitter = 0.2;
35 ret->minVelocity = 0.000000001;
36 ret->minPotentialChange = 0.000000001;
37 ret->minMean = 0.0001;
38 ret->minMeanImprovement = 0.00005;
39 ret->single = AIR_FALSE0;
40 ret->insertZeroVec = AIR_FALSE0;
41 ret->verbose = 1;
42 ret->snap = 0;
43 ret->report = 400;
44 ret->expo = 1;
45 ret->expo_d = 0;
46 ret->seed = 42;
47 ret->maxEdgeShrink = 20;
48 ret->minIteration = 0;
49 ret->maxIteration = 1000000;
50 ret->step = 0;
51 ret->nudge = 0;
52 ret->itersUsed = 0;
53 ret->potential = 0;
54 ret->potentialNorm = 0;
55 ret->angle = 0;
56 ret->edge = 0;
57 }
58 return ret;
59}
60
61tenGradientParm *
62tenGradientParmNix(tenGradientParm *tgparm) {
63
64 airFree(tgparm);
65 return NULL((void*)0);
66}
67
68int
69tenGradientCheck(const Nrrd *ngrad, int type, unsigned int minnum) {
70 static const char me[]="tenGradientCheck";
71
72 if (nrrdCheck(ngrad)) {
73 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: basic validity check failed", me);
74 return 1;
75 }
76 if (!( 3 == ngrad->axis[0].size && 2 == ngrad->dim )) {
77 char stmp[AIR_STRLEN_SMALL(128+1)];
78 biffAddf(TENtenBiffKey, "%s: need a 3xN 2-D array (not a %sx? %u-D array)", me,
79 airSprintSize_t(stmp, ngrad->axis[0].size), ngrad->dim);
80 return 1;
81 }
82 if (nrrdTypeDefault != type && type != ngrad->type) {
83 biffAddf(TENtenBiffKey, "%s: requested type %s but got type %s", me,
84 airEnumStr(nrrdType, type), airEnumStr(nrrdType, ngrad->type));
85 return 1;
86 }
87 if (nrrdTypeBlock == ngrad->type) {
88 biffAddf(TENtenBiffKey, "%s: sorry, can't use %s type", me,
89 airEnumStr(nrrdType, nrrdTypeBlock));
90 return 1;
91 }
92 if (!( minnum <= ngrad->axis[1].size )) {
93 char stmp[AIR_STRLEN_SMALL(128+1)];
94 biffAddf(TENtenBiffKey, "%s: have only %s gradients, need at least %d", me,
95 airSprintSize_t(stmp, ngrad->axis[1].size), minnum);
96 return 1;
97 }
98
99 return 0;
100}
101
102/*
103******** tenGradientRandom
104**
105** generates num random unit vectors of type double
106*/
107int
108tenGradientRandom(Nrrd *ngrad, unsigned int num, unsigned int seed) {
109 static const char me[]="tenGradientRandom";
110 double *grad, len;
111 unsigned int gi;
112
113 if (nrrdMaybeAlloc_va(ngrad, nrrdTypeDouble, 2,
114 AIR_CAST(size_t, 3)((size_t)(3)), AIR_CAST(size_t, num)((size_t)(num)))) {
115 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate output", me);
116 return 1;
117 }
118 airSrandMT(seed);
119 grad = AIR_CAST(double*, ngrad->data)((double*)(ngrad->data));
120 for (gi=0; gi<num; gi++) {
121 do {
122 grad[0] = AIR_AFFINE(0, airDrandMT(), 1, -1, 1)( ((double)(1)-(-1))*((double)(airDrandMT())-(0)) / ((double)
(1)-(0)) + (-1));
123 grad[1] = AIR_AFFINE(0, airDrandMT(), 1, -1, 1)( ((double)(1)-(-1))*((double)(airDrandMT())-(0)) / ((double)
(1)-(0)) + (-1));
124 grad[2] = AIR_AFFINE(0, airDrandMT(), 1, -1, 1)( ((double)(1)-(-1))*((double)(airDrandMT())-(0)) / ((double)
(1)-(0)) + (-1));
125 len = ELL_3V_LEN(grad)(sqrt((((grad))[0]*((grad))[0] + ((grad))[1]*((grad))[1] + ((
grad))[2]*((grad))[2])));
126 } while (len > 1 || !len);
127 ELL_3V_SCALE(grad, 1.0/len, grad)((grad)[0] = (1.0/len)*(grad)[0], (grad)[1] = (1.0/len)*(grad
)[1], (grad)[2] = (1.0/len)*(grad)[2]);
128 grad += 3;
129 }
130 return 0;
131}
132
133/*
134******** tenGradientIdealEdge
135**
136** edge length of delauney triangulation of idealized distribution of
137** N gradients (2*N points), but also allowing a boolean "single" flag
138** saying that we actually care about N points
139*/
140double
141tenGradientIdealEdge(unsigned int N, int single) {
142
143 return sqrt((!single ? 4 : 8)*AIR_PI3.14159265358979323846/(sqrt(3)*N));
144}
145
146/*
147******** tenGradientJitter
148**
149** moves all gradients by amount dist on tangent plane, in a random
150** direction, and then renormalizes. The distance is a fraction
151** of the ideal edge length (via tenGradientIdealEdge)
152*/
153int
154tenGradientJitter(Nrrd *nout, const Nrrd *nin, double dist) {
155 static const char me[]="tenGradientJitter";
156 double *grad, perp0[3], perp1[3], len, theta, cc, ss, edge;
157 unsigned int gi, num;
158
159 if (nrrdConvert(nout, nin, nrrdTypeDouble)) {
160 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: trouble converting input to double", me);
161 return 1;
162 }
163 if (tenGradientCheck(nout, nrrdTypeDouble, 3)) {
164 biffAddf(TENtenBiffKey, "%s: didn't get valid gradients", me);
165 return 1;
166 }
167 grad = AIR_CAST(double*, nout->data)((double*)(nout->data));
168 num = AIR_UINT(nout->axis[1].size)((unsigned int)(nout->axis[1].size));
169 /* HEY: possible confusion between single and not */
170 edge = tenGradientIdealEdge(num, AIR_FALSE0);
171 for (gi=0; gi<num; gi++) {
172 ELL_3V_NORM(grad, grad, len)(len = (sqrt((((grad))[0]*((grad))[0] + ((grad))[1]*((grad))[
1] + ((grad))[2]*((grad))[2]))), ((grad)[0] = (1.0/len)*(grad
)[0], (grad)[1] = (1.0/len)*(grad)[1], (grad)[2] = (1.0/len)*
(grad)[2]));
173 ell_3v_perp_d(perp0, grad);
174 ELL_3V_CROSS(perp1, perp0, grad)((perp1)[0] = (perp0)[1]*(grad)[2] - (perp0)[2]*(grad)[1], (perp1
)[1] = (perp0)[2]*(grad)[0] - (perp0)[0]*(grad)[2], (perp1)[2
] = (perp0)[0]*(grad)[1] - (perp0)[1]*(grad)[0]);
175 theta = AIR_AFFINE(0, airDrandMT(), 1, 0, 2*AIR_PI)( ((double)(2*3.14159265358979323846)-(0))*((double)(airDrandMT
())-(0)) / ((double)(1)-(0)) + (0));
176 cc = dist*edge*cos(theta);
177 ss = dist*edge*sin(theta);
178 ELL_3V_SCALE_ADD3(grad, 1.0, grad, cc, perp0, ss, perp1)((grad)[0] = (1.0)*(grad)[0] + (cc)*(perp0)[0] + (ss)*(perp1)
[0], (grad)[1] = (1.0)*(grad)[1] + (cc)*(perp0)[1] + (ss)*(perp1
)[1], (grad)[2] = (1.0)*(grad)[2] + (cc)*(perp0)[2] + (ss)*(perp1
)[2]);
179 ELL_3V_NORM(grad, grad, len)(len = (sqrt((((grad))[0]*((grad))[0] + ((grad))[1]*((grad))[
1] + ((grad))[2]*((grad))[2]))), ((grad)[0] = (1.0/len)*(grad
)[0], (grad)[1] = (1.0/len)*(grad)[1], (grad)[2] = (1.0/len)*
(grad)[2]));
180 grad += 3;
181 }
182
183 return 0;
184}
185
186void
187tenGradientMeasure(double *pot, double *minAngle, double *minEdge,
188 const Nrrd *npos, tenGradientParm *tgparm,
189 int edgeNormalize) {
190 /* static const char me[]="tenGradientMeasure"; */
191 double diff[3], *pos, atmp=0, ptmp, edge, len;
192 unsigned int ii, jj, num;
193
194 /* allow minAngle NULL */
195 if (!(pot && npos && tgparm )) {
196 return;
197 }
198
199 num = AIR_UINT(npos->axis[1].size)((unsigned int)(npos->axis[1].size));
200 pos = AIR_CAST(double *, npos->data)((double *)(npos->data));
201 edge = (edgeNormalize
202 ? tenGradientIdealEdge(num, tgparm->single)
203 : 1.0);
204 *pot = 0;
205 if (minAngle) {
206 *minAngle = AIR_PI3.14159265358979323846;
207 }
208 if (minEdge) {
209 *minEdge = 2;
210 }
211 for (ii=0; ii<num; ii++) {
212 for (jj=0; jj<ii; jj++) {
213 ELL_3V_SUB(diff, pos + 3*ii, pos + 3*jj)((diff)[0] = (pos + 3*ii)[0] - (pos + 3*jj)[0], (diff)[1] = (
pos + 3*ii)[1] - (pos + 3*jj)[1], (diff)[2] = (pos + 3*ii)[2]
- (pos + 3*jj)[2]);
214 len = ELL_3V_LEN(diff)(sqrt((((diff))[0]*((diff))[0] + ((diff))[1]*((diff))[1] + ((
diff))[2]*((diff))[2])));
215 if (minEdge) {
216 *minEdge = AIR_MIN(*minEdge, len)((*minEdge) < (len) ? (*minEdge) : (len));
217 }
218 if (tgparm->expo) {
219 ptmp = airIntPow(edge/len, tgparm->expo);
220 } else {
221 ptmp = pow(edge/len, tgparm->expo_d);
222 }
223 *pot += ptmp;
224 if (minAngle) {
225 atmp = ell_3v_angle_d(pos + 3*ii, pos + 3*jj);
226 *minAngle = AIR_MIN(atmp, *minAngle)((atmp) < (*minAngle) ? (atmp) : (*minAngle));
227 }
228 if (!tgparm->single) {
229 *pot += ptmp;
230 ELL_3V_ADD2(diff, pos + 3*ii, pos + 3*jj)((diff)[0] = (pos + 3*ii)[0] + (pos + 3*jj)[0], (diff)[1] = (
pos + 3*ii)[1] + (pos + 3*jj)[1], (diff)[2] = (pos + 3*ii)[2]
+ (pos + 3*jj)[2]);
231 len = ELL_3V_LEN(diff)(sqrt((((diff))[0]*((diff))[0] + ((diff))[1]*((diff))[1] + ((
diff))[2]*((diff))[2])));
232 if (minEdge) {
233 *minEdge = AIR_MIN(*minEdge, len)((*minEdge) < (len) ? (*minEdge) : (len));
234 }
235 if (tgparm->expo) {
236 *pot += 2*airIntPow(edge/len, tgparm->expo);
237 } else {
238 *pot += 2*pow(edge/len, tgparm->expo_d);
239 }
240 if (minAngle) {
241 *minAngle = AIR_MIN(AIR_PI-atmp, *minAngle)((3.14159265358979323846 -atmp) < (*minAngle) ? (3.14159265358979323846
-atmp) : (*minAngle));
242 }
243 }
244 }
245 }
246 return;
247}
248
249/*
250** Do asynchronous update of positions in "npos', based on force
251** calculations wherein the distances are normalized "edge". Using a
252** small "edge" allows forces to either underflow to zero, or be
253** finite, instead of exploding to infinity, for high exponents.
254**
255** The smallest seen edge length is recorded in "*edgeMin", which is
256** initialized to the given "edge". This allows, for example, the
257** caller to try again with a smaller edge normalization.
258**
259** The mean velocity of the points through the update is recorded in
260** "*meanVel".
261**
262** Based on the observation that when using large exponents, numerical
263** difficulties arise from the (force-based) update of the positions
264** of the two (or few) closest particles, this function puts a speed
265** limit (variable "limit") on the distance a particle may move during
266** update, expressed as a fraction of the normalizing edge length.
267** "limit" has been set heuristically, according to the exponent (we
268** have to clamp speeds more aggresively with higher exponents), as
269** well as (even more heuristically) according to the number of times
270** the step size has been decreased. This latter factor has to be
271** bounded, so that the update is not unnecessarily bounded when the
272** step size gets very small at the last stages of computation.
273** Without the step-size-based speed limit, the step size would
274** sometimes (e.g. num=200, expo=300) have to reduced to a miniscule
275** value, which slows subsequent convergence terribly.
276**
277** this function is not static, though it could be, so that mac's
278** "Sampler" app can profile this
279*/
280int
281_tenGradientUpdate(double *meanVel, double *edgeMin,
282 Nrrd *npos, double edge, tenGradientParm *tgparm) {
283 /* static const char me[]="_tenGradientUpdate"; */
284 double *pos, newpos[3], grad[3], ngrad[3],
285 dir[3], len, rep, step, diff[3], limit, expo;
286 int num, ii, jj, E;
287
288 E = 0;
289 pos = AIR_CAST(double *, npos->data)((double *)(npos->data));
290 num = AIR_UINT(npos->axis[1].size)((unsigned int)(npos->axis[1].size));
291 *meanVel = 0;
292 *edgeMin = edge;
293 expo = tgparm->expo ? tgparm->expo : tgparm->expo_d;
294 limit = expo*AIR_MIN(sqrt(expo),((sqrt(expo)) < (log(1 + tgparm->initStep/tgparm->step
)) ? (sqrt(expo)) : (log(1 + tgparm->initStep/tgparm->step
)))
295 log(1 + tgparm->initStep/tgparm->step))((sqrt(expo)) < (log(1 + tgparm->initStep/tgparm->step
)) ? (sqrt(expo)) : (log(1 + tgparm->initStep/tgparm->step
)));
296 for (ii=0; ii<num; ii++) {
297 ELL_3V_SET(grad, 0, 0, 0)((grad)[0] = (0), (grad)[1] = (0), (grad)[2] = (0));
298 for (jj=0; jj<num; jj++) {
299 if (ii == jj) {
300 continue;
301 }
302 ELL_3V_SUB(dir, pos + 3*ii, pos + 3*jj)((dir)[0] = (pos + 3*ii)[0] - (pos + 3*jj)[0], (dir)[1] = (pos
+ 3*ii)[1] - (pos + 3*jj)[1], (dir)[2] = (pos + 3*ii)[2] - (
pos + 3*jj)[2]);
303 ELL_3V_NORM(dir, dir, len)(len = (sqrt((((dir))[0]*((dir))[0] + ((dir))[1]*((dir))[1] +
((dir))[2]*((dir))[2]))), ((dir)[0] = (1.0/len)*(dir)[0], (dir
)[1] = (1.0/len)*(dir)[1], (dir)[2] = (1.0/len)*(dir)[2]));
304 *edgeMin = AIR_MIN(*edgeMin, len)((*edgeMin) < (len) ? (*edgeMin) : (len));
305 if (tgparm->expo) {
306 rep = airIntPow(edge/len, tgparm->expo+1);
307 } else {
308 rep = pow(edge/len, tgparm->expo_d+1);
309 }
310 ELL_3V_SCALE_INCR(grad, rep/num, dir)((grad)[0] += (rep/num)*(dir)[0], (grad)[1] += (rep/num)*(dir
)[1], (grad)[2] += (rep/num)*(dir)[2]);
311 if (!tgparm->single) {
312 ELL_3V_ADD2(dir, pos + 3*ii, pos + 3*jj)((dir)[0] = (pos + 3*ii)[0] + (pos + 3*jj)[0], (dir)[1] = (pos
+ 3*ii)[1] + (pos + 3*jj)[1], (dir)[2] = (pos + 3*ii)[2] + (
pos + 3*jj)[2]);
313 ELL_3V_NORM(dir, dir, len)(len = (sqrt((((dir))[0]*((dir))[0] + ((dir))[1]*((dir))[1] +
((dir))[2]*((dir))[2]))), ((dir)[0] = (1.0/len)*(dir)[0], (dir
)[1] = (1.0/len)*(dir)[1], (dir)[2] = (1.0/len)*(dir)[2]));
314 *edgeMin = AIR_MIN(*edgeMin, len)((*edgeMin) < (len) ? (*edgeMin) : (len));
315 if (tgparm->expo) {
316 rep = airIntPow(edge/len, tgparm->expo+1);
317 } else {
318 rep = pow(edge/len, tgparm->expo_d+1);
319 }
320 ELL_3V_SCALE_INCR(grad, rep/num, dir)((grad)[0] += (rep/num)*(dir)[0], (grad)[1] += (rep/num)*(dir
)[1], (grad)[2] += (rep/num)*(dir)[2]);
321 }
322 }
323 ELL_3V_NORM(ngrad, grad, len)(len = (sqrt((((grad))[0]*((grad))[0] + ((grad))[1]*((grad))[
1] + ((grad))[2]*((grad))[2]))), ((ngrad)[0] = (1.0/len)*(grad
)[0], (ngrad)[1] = (1.0/len)*(grad)[1], (ngrad)[2] = (1.0/len
)*(grad)[2]));
324 if (!( AIR_EXISTS(len)(((int)(!((len) - (len))))) )) {
325 /* things blew up, either in incremental force
326 additions, or in the attempt at normalization */
327 E = 1;
328 *meanVel = AIR_NAN(airFloatQNaN.f);
329 break;
330 }
331 if (0 == len) {
332 /* if the length of grad[] underflowed to zero, we can
333 legitimately zero out ngrad[] */
334 ELL_3V_SET(ngrad, 0, 0, 0)((ngrad)[0] = (0), (ngrad)[1] = (0), (ngrad)[2] = (0));
335 }
336 step = AIR_MIN(len*tgparm->step, edge/limit)((len*tgparm->step) < (edge/limit) ? (len*tgparm->step
) : (edge/limit));
337 ELL_3V_SCALE_ADD2(newpos,((newpos)[0] = (1.0)*(pos + 3*ii)[0] + (step)*(ngrad)[0], (newpos
)[1] = (1.0)*(pos + 3*ii)[1] + (step)*(ngrad)[1], (newpos)[2]
= (1.0)*(pos + 3*ii)[2] + (step)*(ngrad)[2])
338 1.0, pos + 3*ii,((newpos)[0] = (1.0)*(pos + 3*ii)[0] + (step)*(ngrad)[0], (newpos
)[1] = (1.0)*(pos + 3*ii)[1] + (step)*(ngrad)[1], (newpos)[2]
= (1.0)*(pos + 3*ii)[2] + (step)*(ngrad)[2])
339 step, ngrad)((newpos)[0] = (1.0)*(pos + 3*ii)[0] + (step)*(ngrad)[0], (newpos
)[1] = (1.0)*(pos + 3*ii)[1] + (step)*(ngrad)[1], (newpos)[2]
= (1.0)*(pos + 3*ii)[2] + (step)*(ngrad)[2]);
340 ELL_3V_NORM(newpos, newpos, len)(len = (sqrt((((newpos))[0]*((newpos))[0] + ((newpos))[1]*((newpos
))[1] + ((newpos))[2]*((newpos))[2]))), ((newpos)[0] = (1.0/len
)*(newpos)[0], (newpos)[1] = (1.0/len)*(newpos)[1], (newpos)[
2] = (1.0/len)*(newpos)[2]));
341 ELL_3V_SUB(diff, pos + 3*ii, newpos)((diff)[0] = (pos + 3*ii)[0] - (newpos)[0], (diff)[1] = (pos +
3*ii)[1] - (newpos)[1], (diff)[2] = (pos + 3*ii)[2] - (newpos
)[2]);
342 *meanVel += ELL_3V_LEN(diff)(sqrt((((diff))[0]*((diff))[0] + ((diff))[1]*((diff))[1] + ((
diff))[2]*((diff))[2])));
343 ELL_3V_COPY(pos + 3*ii, newpos)((pos + 3*ii)[0] = (newpos)[0], (pos + 3*ii)[1] = (newpos)[1]
, (pos + 3*ii)[2] = (newpos)[2]);
344 }
345 *meanVel /= num;
346
347 return E;
348}
349
350/*
351** assign random signs to the vectors and measures the length of their
352** mean, as quickly as possible
353*/
354static double
355party(Nrrd *npos, airRandMTState *rstate) {
356 double *pos, mean[3];
357 unsigned int ii, num, rnd, rndBit;
358
359 pos = (double *)(npos->data);
360 num = AIR_UINT(npos->axis[1].size)((unsigned int)(npos->axis[1].size));
361 rnd = airUIrandMT_r(rstate);
362 rndBit = 0;
363 ELL_3V_SET(mean, 0, 0, 0)((mean)[0] = (0), (mean)[1] = (0), (mean)[2] = (0));
364 for (ii=0; ii<num; ii++) {
365 if (32 == rndBit) {
366 rnd = airUIrandMT_r(rstate);
367 rndBit = 0;
368 }
369 if (rnd & (1 << rndBit++)) {
370 ELL_3V_SCALE(pos + 3*ii, -1, pos + 3*ii)((pos + 3*ii)[0] = (-1)*(pos + 3*ii)[0], (pos + 3*ii)[1] = (-
1)*(pos + 3*ii)[1], (pos + 3*ii)[2] = (-1)*(pos + 3*ii)[2]);
371 }
372 ELL_3V_INCR(mean, pos + 3*ii)((mean)[0] += (pos + 3*ii)[0], (mean)[1] += (pos + 3*ii)[1], (
mean)[2] += (pos + 3*ii)[2]);
373 }
374 ELL_3V_SCALE(mean, 1.0/num, mean)((mean)[0] = (1.0/num)*(mean)[0], (mean)[1] = (1.0/num)*(mean
)[1], (mean)[2] = (1.0/num)*(mean)[2]);
375 return ELL_3V_LEN(mean)(sqrt((((mean))[0]*((mean))[0] + ((mean))[1]*((mean))[1] + ((
mean))[2]*((mean))[2])));
376}
377
378/*
379** parties until the gradients settle down
380*/
381int
382tenGradientBalance(Nrrd *nout, const Nrrd *nin,
383 tenGradientParm *tgparm) {
384 static const char me[]="tenGradientBalance";
385 double len, lastLen, improv;
386 airRandMTState *rstate;
387 Nrrd *ncopy;
388 unsigned int iter, maxIter;
389 int done;
390 airArray *mop;
391
392 if (!nout || tenGradientCheck(nin, nrrdTypeUnknown, 2) || !tgparm) {
393 biffAddf(TENtenBiffKey, "%s: got NULL pointer (%p,%p) or invalid nin", me,
394 AIR_VOIDP(nout)((void *)(nout)), AIR_VOIDP(tgparm)((void *)(tgparm)));
395 return 1;
396 }
397 if (nrrdConvert(nout, nin, nrrdTypeDouble)) {
398 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: can't initialize output with input", me);
399 return 1;
400 }
401
402 mop = airMopNew();
403 ncopy = nrrdNew();
404 airMopAdd(mop, ncopy, (airMopper)nrrdNuke, airMopAlways);
405 rstate = airRandMTStateNew(tgparm->seed);
406 airMopAdd(mop, rstate, (airMopper)airRandMTStateNix, airMopAlways);
407 /* HEY: factor of 100 is an approximate hack */
408 maxIter = 100*tgparm->maxIteration;
409
410 lastLen = 1.0;
411 done = AIR_FALSE0;
Value stored to 'done' is never read
412 do {
413 iter = 0;
414 do {
415 iter++;
416 len = party(nout, rstate);
417 } while (len > lastLen && iter < maxIter);
418 if (iter >= maxIter) {
419 if (tgparm->verbose) {
420 fprintf(stderr__stderrp, "%s: stopping at max iter %u\n", me, maxIter);
421 }
422 if (nrrdCopy(nout, ncopy)) {
423 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: trouble copying", me);
424 airMopError(mop); return 1;
425 }
426 done = AIR_TRUE1;
427 } else {
428 if (nrrdCopy(ncopy, nout)) {
429 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: trouble copying", me);
430 airMopError(mop); return 1;
431 }
432 improv = lastLen - len;
433 lastLen = len;
434 if (tgparm->verbose) {
435 fprintf(stderr__stderrp, "%s: (iter %u) improvement: %g (mean length = %g)\n",
436 me, iter, improv, len);
437 }
438 done = (improv <= tgparm->minMeanImprovement
439 || len < tgparm->minMean);
440 }
441 } while (!done);
442
443 airMopOkay(mop);
444 return 0;
445}
446
447/*
448******** tenGradientDistribute
449**
450** Takes the given list of gradients, normalizes their lengths,
451** optionally jitters their positions, does point repulsion, and then
452** (optionally) selects a combination of directions with minimum vector sum.
453**
454** The complicated part of this is the point repulsion, which uses a
455** gradient descent with variable set size. The progress of the system
456** is measured by decrease in potential (when its measurement doesn't
457** overflow to infinity) or an increase in the minimum angle. When a
458** step results in negative progress, the step size is halved, and the
459** iteration is attempted again. Based on the observation that at
460** some points the step size must be made very small to get progress,
461** the step size is cautiously increased ("nudged") at every
462** iteration, to try to avoid using an overly small step. The amount
463** by which the step is nudged is halved everytime the step is halved,
464** to avoid endless cycling through step sizes.
465*/
466int
467tenGradientDistribute(Nrrd *nout, const Nrrd *nin,
468 tenGradientParm *tgparm) {
469 static const char me[]="tenGradientDistribute";
470 char filename[AIR_STRLEN_SMALL(128+1)];
471 unsigned int ii, num, iter, oldIdx, newIdx, edgeShrink;
472 airArray *mop;
473 Nrrd *npos[2];
474 double *pos, len, meanVelocity, pot, potNew, potD,
475 edge, edgeMin, angle, angleNew;
476 int E;
477
478 if (!nout || tenGradientCheck(nin, nrrdTypeUnknown, 2) || !tgparm) {
479 biffAddf(TENtenBiffKey, "%s: got NULL pointer or invalid input", me);
480 return 1;
481 }
482
483 num = AIR_UINT(nin->axis[1].size)((unsigned int)(nin->axis[1].size));
484 mop = airMopNew();
485 npos[0] = nrrdNew();
486 npos[1] = nrrdNew();
487 airMopAdd(mop, npos[0], (airMopper)nrrdNuke, airMopAlways);
488 airMopAdd(mop, npos[1], (airMopper)nrrdNuke, airMopAlways);
489 if (nrrdConvert(npos[0], nin, nrrdTypeDouble)
490 || nrrdConvert(npos[1], nin, nrrdTypeDouble)) {
491 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: trouble allocating temp buffers", me);
492 airMopError(mop); return 1;
493 }
494
495 pos = (double*)(npos[0]->data);
496 for (ii=0; ii<num; ii++) {
497 ELL_3V_NORM(pos, pos, len)(len = (sqrt((((pos))[0]*((pos))[0] + ((pos))[1]*((pos))[1] +
((pos))[2]*((pos))[2]))), ((pos)[0] = (1.0/len)*(pos)[0], (pos
)[1] = (1.0/len)*(pos)[1], (pos)[2] = (1.0/len)*(pos)[2]));
498 pos += 3;
499 }
500 if (tgparm->jitter) {
501 if (tenGradientJitter(npos[0], npos[0], tgparm->jitter)) {
502 biffAddf(TENtenBiffKey, "%s: problem jittering input", me);
503 airMopError(mop); return 1;
504 }
505 }
506
507 /* initialize things prior to first iteration; have to
508 make sure that loop body tests pass 1st time around */
509 meanVelocity = 2*tgparm->minVelocity;
510 potD = -2*tgparm->minPotentialChange;
511 oldIdx = 0;
512 newIdx = 1;
513 tgparm->step = tgparm->initStep;
514 tgparm->nudge = 0.1;
515 tenGradientMeasure(&pot, &angle, NULL((void*)0),
516 npos[oldIdx], tgparm, AIR_TRUE1);
517 for (iter = 0;
518 ((!!tgparm->minIteration && iter < tgparm->minIteration)
519 ||
520 (iter < tgparm->maxIteration
521 && (!tgparm->minPotentialChange
522 || !AIR_EXISTS(potD)(((int)(!((potD) - (potD)))))
523 || -potD > tgparm->minPotentialChange)
524 && (!tgparm->minVelocity
525 || meanVelocity > tgparm->minVelocity)
526 && tgparm->step > FLT_MIN1.17549435e-38F));
527 iter++) {
528 /* copy positions from old to new */
529 memcpy(npos[newIdx]->data, npos[oldIdx]->data, 3*num*sizeof(double))__builtin___memcpy_chk (npos[newIdx]->data, npos[oldIdx]->
data, 3*num*sizeof(double), __builtin_object_size (npos[newIdx
]->data, 0));
530 edge = tenGradientIdealEdge(num, tgparm->single);
531 edgeShrink = 0;
532 /* try to do a position update, which will fail if repulsion values
533 explode, from having an insufficiently small edge normalization,
534 so retry with smaller edge next time */
535 do {
536 E = _tenGradientUpdate(&meanVelocity, &edgeMin,
537 npos[newIdx], edge, tgparm);
538 if (E) {
539 if (edgeShrink > tgparm->maxEdgeShrink) {
540 biffAddf(TENtenBiffKey, "%s: %u > %u edge shrinks (%g), update still failed",
541 me, edgeShrink, tgparm->maxEdgeShrink, edge);
542 airMopError(mop); return 1;
543 }
544 edgeShrink++;
545 /* re-initialize positions (HEY ugly code logic) */
546 memcpy(npos[newIdx]->data, npos[oldIdx]->data, 3*num*sizeof(double))__builtin___memcpy_chk (npos[newIdx]->data, npos[oldIdx]->
data, 3*num*sizeof(double), __builtin_object_size (npos[newIdx
]->data, 0));
547 edge = edgeMin;
548 }
549 } while (E);
550 tenGradientMeasure(&potNew, &angleNew, NULL((void*)0),
551 npos[newIdx], tgparm, AIR_TRUE1);
552 if ((AIR_EXISTS(pot)(((int)(!((pot) - (pot))))) && AIR_EXISTS(potNew)(((int)(!((potNew) - (potNew))))) && potNew <= pot)
553 || angleNew >= angle) {
554 /* there was progress of some kind, either through potential
555 decrease, or angle increase */
556 potD = 2*(potNew - pot)/(potNew + pot);
557 if (!(iter % tgparm->report) && tgparm->verbose) {
558 fprintf(stderr__stderrp, "%s(%d): . . . . . . step = %g, edgeShrink = %u\n"
559 " velo = %g<>%g, phi = %g ~ %g<>%g, angle = %g ~ %g\n",
560 me, iter, tgparm->step, edgeShrink,
561 meanVelocity, tgparm->minVelocity,
562 pot, potD, tgparm->minPotentialChange,
563 angle, angleNew - angle);
564 }
565 if (tgparm->snap && !(iter % tgparm->snap)) {
566 sprintf(filename, "%05d.nrrd", iter/tgparm->snap)__builtin___sprintf_chk (filename, 0, __builtin_object_size (
filename, 2 > 1 ? 1 : 0), "%05d.nrrd", iter/tgparm->snap
);
567 if (tgparm->verbose) {
568 fprintf(stderr__stderrp, "%s(%d): . . . . . . saving %s\n",
569 me, iter, filename);
570 }
571 if (nrrdSave(filename, npos[newIdx], NULL((void*)0))) {
572 char *serr;
573 serr = biffGetDone(NRRDnrrdBiffKey);
574 if (tgparm->verbose) { /* perhaps shouldn't have this check */
575 fprintf(stderr__stderrp, "%s: iter=%d, couldn't save snapshot:\n%s"
576 "continuing ...\n", me, iter, serr);
577 }
578 free(serr);
579 }
580 }
581 tgparm->step *= 1 + tgparm->nudge;
582 tgparm->step = AIR_MIN(tgparm->initStep, tgparm->step)((tgparm->initStep) < (tgparm->step) ? (tgparm->initStep
) : (tgparm->step));
583 pot = potNew;
584 angle = angleNew;
585 /* swap buffers */
586 newIdx = 1 - newIdx;
587 oldIdx = 1 - oldIdx;
588 } else {
589 /* oops, did not make progress; back off and try again */
590 if (tgparm->verbose) {
591 fprintf(stderr__stderrp, "%s(%d): ######## step %g --> %g\n"
592 " phi = %g --> %g ~ %g, angle = %g --> %g\n",
593 me, iter, tgparm->step, tgparm->step/2,
594 pot, potNew, potD, angle, angleNew);
595 }
596 tgparm->step /= 2;
597 tgparm->nudge /= 2;
598 }
599 }
600
601 /* when the for-loop test fails, we stop before computing the next
602 iteration (which starts with copying from npos[oldIdx] to
603 npos[newIdx]) ==> the final results are in npos[oldIdx] */
604
605 if (tgparm->verbose) {
606 fprintf(stderr__stderrp, "%s: .......................... done distribution:\n"
607 " (%d && %d) || (%d \n"
608 " && (%d || %d || %d) \n"
609 " && (%d || %d) \n"
610 " && %d) is false\n", me,
611 !!tgparm->minIteration, iter < tgparm->minIteration,
612 iter < tgparm->maxIteration,
613 !tgparm->minPotentialChange,
614 !AIR_EXISTS(potD)(((int)(!((potD) - (potD))))), AIR_ABS(potD)((potD) > 0.0f ? (potD) : -(potD)) > tgparm->minPotentialChange,
615 !tgparm->minVelocity, meanVelocity > tgparm->minVelocity,
616 tgparm->step > FLT_MIN1.17549435e-38F);
617 fprintf(stderr__stderrp, " iter=%d, velo = %g<>%g, phi = %g ~ %g<>%g;\n",
618 iter, meanVelocity, tgparm->minVelocity, pot,
619 potD, tgparm->minPotentialChange);
620 fprintf(stderr__stderrp, " minEdge = %g; idealEdge = %g\n",
621 2*sin(angle/2), tenGradientIdealEdge(num, tgparm->single));
622 }
623
624 tenGradientMeasure(&pot, NULL((void*)0), NULL((void*)0), npos[oldIdx], tgparm, AIR_FALSE0);
625 tgparm->potential = pot;
626 tenGradientMeasure(&pot, &angle, &edge, npos[oldIdx], tgparm, AIR_TRUE1);
627 tgparm->potentialNorm = pot;
628 tgparm->angle = angle;
629 tgparm->edge = edge;
630 tgparm->itersUsed = iter;
631
632 if ((tgparm->minMeanImprovement || tgparm->minMean)
633 && !tgparm->single) {
634 if (tgparm->verbose) {
635 fprintf(stderr__stderrp, "%s: optimizing balance:\n", me);
636 }
637 if (tenGradientBalance(nout, npos[oldIdx], tgparm)) {
638 biffAddf(TENtenBiffKey, "%s: failed to minimize vector sum of gradients", me);
639 airMopError(mop); return 1;
640 }
641 if (tgparm->verbose) {
642 fprintf(stderr__stderrp, "%s: .......................... done balancing.\n", me);
643 }
644 } else {
645 if (tgparm->verbose) {
646 fprintf(stderr__stderrp, "%s: .......................... (no balancing)\n", me);
647 }
648 if (nrrdConvert(nout, npos[oldIdx], nrrdTypeDouble)) {
649 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: couldn't set output", me);
650 airMopError(mop); return 1;
651 }
652 }
653
654 airMopOkay(mop);
655 return 0;
656}
657
658/*
659** note that if tgparm->insertZeroVec, there will be one sample more
660** along axis 1 of nout than the requested #gradients "num"
661*/
662int
663tenGradientGenerate(Nrrd *nout, unsigned int num, tenGradientParm *tgparm) {
664 static const char me[]="tenGradientGenerate";
665 Nrrd *nin;
666 airArray *mop;
667
668 if (!(nout && tgparm)) {
669 biffAddf(TENtenBiffKey, "%s: got NULL pointer", me);
670 return 1;
671 }
672 if (!( num >= 3 )) {
673 biffAddf(TENtenBiffKey, "%s: can generate minimum of 3 gradient directions "
674 "(not %d)", me, num);
675 return 1;
676 }
677 mop = airMopNew();
678 nin = nrrdNew();
679 airMopAdd(mop, nin, (airMopper)nrrdNuke, airMopAlways);
680
681 if (tenGradientRandom(nin, num, tgparm->seed)
682 || tenGradientDistribute(nout, nin, tgparm)) {
683 biffAddf(TENtenBiffKey, "%s: trouble", me);
684 airMopError(mop); return 1;
685 }
686 if (tgparm->insertZeroVec) {
687 /* this is potentially confusing: the second axis (axis 1)
688 is going to come back one longer than the requested
689 number of gradients! */
690 Nrrd *ntmp;
691 ptrdiff_t padMin[2] = {0, -1}, padMax[2];
692 padMax[0] = AIR_CAST(ptrdiff_t, nout->axis[0].size-1)((ptrdiff_t)(nout->axis[0].size-1));
693 padMax[1] = AIR_CAST(ptrdiff_t, num-1)((ptrdiff_t)(num-1));
694 ntmp = nrrdNew();
695 airMopAdd(mop, ntmp, (airMopper)nrrdNuke, airMopAlways);
696 if (nrrdPad_nva(ntmp, nout, padMin, padMax,
697 nrrdBoundaryPad, 0.0)
698 || nrrdCopy(nout, ntmp)) {
699 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: trouble adding zero vector", me);
700 airMopError(mop); return 1;
701 }
702 }
703
704 airMopOkay(mop);
705 return 0;
706}