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 |
|
|
|
27 |
|
|
#define INFO "Calculate balanced gradient directions for DWI acquisition" |
28 |
|
|
static const char *_tend_gradsInfoL = |
29 |
|
|
(INFO |
30 |
|
|
", based on a simulation of anti-podal point pairs repelling each other " |
31 |
|
|
"on the unit sphere surface. This can either distribute more uniformly " |
32 |
|
|
"a given set of gradients, or it can make a new distribution from scratch. " |
33 |
|
|
"A more clever implementation could decrease drag with time, as the " |
34 |
|
|
"solution converges, to get closer to the minimum energy configuration " |
35 |
|
|
"faster. In the mean time, you can run a second pass on the output of " |
36 |
|
|
"the first pass, using lower drag. A second phase of the algorithm " |
37 |
|
|
"tries sign changes in gradient directions in trying to find an optimally " |
38 |
|
|
"balanced set of directions. This uses a randomized search, so if it " |
39 |
|
|
"doesn't seem to be finishing in a reasonable amount of time, try " |
40 |
|
|
"restarting with a different \"-seed\"."); |
41 |
|
|
|
42 |
|
|
int |
43 |
|
|
tend_gradsMain(int argc, const char **argv, const char *me, |
44 |
|
|
hestParm *hparm) { |
45 |
|
|
int pret; |
46 |
|
2 |
hestOpt *hopt = NULL; |
47 |
|
1 |
char *perr, *err; |
48 |
|
|
airArray *mop; |
49 |
|
|
|
50 |
|
1 |
int num, E; |
51 |
|
1 |
Nrrd *nin, *nout; |
52 |
|
1 |
char *outS; |
53 |
|
|
tenGradientParm *tgparm; |
54 |
|
1 |
unsigned int seed; |
55 |
|
|
|
56 |
|
1 |
mop = airMopNew(); |
57 |
|
1 |
tgparm = tenGradientParmNew(); |
58 |
|
1 |
airMopAdd(mop, tgparm, (airMopper)tenGradientParmNix, airMopAlways); |
59 |
|
|
|
60 |
|
1 |
hestOptAdd(&hopt, "n", "# dir", airTypeInt, 1, 1, &num, "6", |
61 |
|
|
"desired number of diffusion gradient directions"); |
62 |
|
1 |
hestOptAdd(&hopt, "i", "grads", airTypeOther, 1, 1, &nin, "", |
63 |
|
|
"initial gradient directions to start with, instead " |
64 |
|
|
"of default random initial directions (overrides \"-n\")", |
65 |
|
1 |
NULL, NULL, nrrdHestNrrd); |
66 |
|
1 |
hestOptAdd(&hopt, "seed", "value", airTypeUInt, 1, 1, &seed, "42", |
67 |
|
|
"seed value to used with airSrandMT()"); |
68 |
|
1 |
hestOptAdd(&hopt, "step", "step", airTypeDouble, 1, 1, &(tgparm->initStep), |
69 |
|
|
"1.0", |
70 |
|
|
"time increment in solver"); |
71 |
|
1 |
hestOptAdd(&hopt, "single", NULL, airTypeInt, 0, 0, &(tgparm->single), NULL, |
72 |
|
|
"instead of the default behavior of tracking a pair of " |
73 |
|
|
"antipodal points (appropriate for determining DWI gradients), " |
74 |
|
|
"use only single points (appropriate for who knows what)."); |
75 |
|
1 |
hestOptAdd(&hopt, "snap", "interval", airTypeInt, 1, 1, &(tgparm->snap), "0", |
76 |
|
|
"specifies an interval between which snapshots of the point " |
77 |
|
|
"positions should be saved out. By default (not using this " |
78 |
|
|
"option), there is no such snapshot behavior"); |
79 |
|
1 |
hestOptAdd(&hopt, "jitter", "jitter", airTypeDouble, 1, 1, |
80 |
|
1 |
&(tgparm->jitter), "0.1", |
81 |
|
|
"amount by which to perturb points when given an input nrrd"); |
82 |
|
1 |
hestOptAdd(&hopt, "miniter", "# iters", airTypeInt, 1, 1, |
83 |
|
1 |
&(tgparm->minIteration), "0", |
84 |
|
|
"max number of iterations for which to run the simulation"); |
85 |
|
1 |
hestOptAdd(&hopt, "maxiter", "# iters", airTypeInt, 1, 1, |
86 |
|
1 |
&(tgparm->maxIteration), "1000000", |
87 |
|
|
"max number of iterations for which to run the simulation"); |
88 |
|
1 |
hestOptAdd(&hopt, "minvelo", "vel", airTypeDouble, 1, 1, |
89 |
|
1 |
&(tgparm->minVelocity), "0.00001", |
90 |
|
|
"low threshold on mean velocity of repelling points, " |
91 |
|
|
"at which point repulsion phase of algorithm terminates."); |
92 |
|
1 |
hestOptAdd(&hopt, "exp", "exponent", airTypeDouble, 1, 1, |
93 |
|
1 |
&(tgparm->expo_d), "1", |
94 |
|
|
"the exponent n that determines the potential energy 1/r^n."); |
95 |
|
1 |
hestOptAdd(&hopt, "dp", "potential change", airTypeDouble, 1, 1, |
96 |
|
1 |
&(tgparm->minPotentialChange), "0.000000001", |
97 |
|
|
"low threshold on fractional change of potential at " |
98 |
|
|
"which point repulsion phase of algorithm terminates."); |
99 |
|
1 |
hestOptAdd(&hopt, "minimprov", "delta", airTypeDouble, 1, 1, |
100 |
|
1 |
&(tgparm->minMeanImprovement), "0.00005", |
101 |
|
|
"in the second phase of the algorithm, " |
102 |
|
|
"when stochastically balancing the sign of the gradients, " |
103 |
|
|
"the (small) improvement in length of mean gradient " |
104 |
|
|
"which triggers termination (as further improvements " |
105 |
|
|
"are unlikely."); |
106 |
|
1 |
hestOptAdd(&hopt, "minmean", "len", airTypeDouble, 1, 1, |
107 |
|
1 |
&(tgparm->minMean), "0.0001", |
108 |
|
|
"if length of mean gradient falls below this, finish " |
109 |
|
|
"the balancing phase"); |
110 |
|
1 |
hestOptAdd(&hopt, "izv", "insert", airTypeBool, 1, 1, |
111 |
|
1 |
&(tgparm->insertZeroVec), "false", |
112 |
|
|
"adding zero vector at beginning of grads"); |
113 |
|
1 |
hestOptAdd(&hopt, "o", "filename", airTypeString, 1, 1, &outS, "-", |
114 |
|
|
"file to write output nrrd to"); |
115 |
|
1 |
airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways); |
116 |
✓✗ |
2 |
USAGE(_tend_gradsInfoL); |
117 |
|
|
JUSTPARSE(); |
118 |
|
|
airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways); |
119 |
|
|
nout = nrrdNew(); |
120 |
|
|
airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways); |
121 |
|
|
|
122 |
|
|
/* see if it was an integral exponent */ |
123 |
|
|
tgparm->expo = AIR_CAST(unsigned int, tgparm->expo_d); |
124 |
|
|
if (tgparm->expo == tgparm->expo_d) { |
125 |
|
|
/* ooo, it was */ |
126 |
|
|
tgparm->expo_d = 0; |
127 |
|
|
} else { |
128 |
|
|
/* no, its non-integral, indicate this as follows */ |
129 |
|
|
tgparm->expo = 0; |
130 |
|
|
} |
131 |
|
|
tgparm->seed = seed; |
132 |
|
|
if (tgparm->snap) { |
133 |
|
|
tgparm->report = tgparm->snap; |
134 |
|
|
} |
135 |
|
|
E = (nin |
136 |
|
|
? tenGradientDistribute(nout, nin, tgparm) |
137 |
|
|
: tenGradientGenerate(nout, num, tgparm)); |
138 |
|
|
if (E) { |
139 |
|
|
airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways); |
140 |
|
|
fprintf(stderr, "%s: trouble making distribution:\n%s\n", me, err); |
141 |
|
|
airMopError(mop); return 1; |
142 |
|
|
} |
143 |
|
|
|
144 |
|
|
if (nrrdSave(outS, nout, NULL)) { |
145 |
|
|
airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways); |
146 |
|
|
fprintf(stderr, "%s: trouble writing:\n%s\n", me, err); |
147 |
|
|
airMopError(mop); return 1; |
148 |
|
|
} |
149 |
|
|
|
150 |
|
|
airMopOkay(mop); |
151 |
|
|
return 0; |
152 |
|
1 |
} |
153 |
|
|
TEND_CMD(grads, INFO); |