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 "Generate a pretty synthetic DT volume" |
28 |
|
|
static const char *_tend_satinInfoL = |
29 |
|
|
(INFO |
30 |
|
|
". The surface of a sphere or torus is covered with either linear or " |
31 |
|
|
"planar anisotropic tensors, or somewhere in between."); |
32 |
|
|
|
33 |
|
|
void |
34 |
|
|
tend_satinSphereEigen(float *eval, float *evec, float x, float y, float z, |
35 |
|
|
float parm, float mina, float maxa, |
36 |
|
|
float thick, float bnd, float evsc) { |
37 |
|
|
float aniso, bound1, bound2, r, norm, tmp[3]; |
38 |
|
|
|
39 |
|
|
r = AIR_CAST(float, sqrt(x*x + y*y + z*z)); |
40 |
|
|
/* 1 on inside, 0 on outside */ |
41 |
|
|
bound1 = AIR_CAST(float, 0.5 - 0.5*airErf((r-0.9)/(bnd + 0.0001))); |
42 |
|
|
/* other way around */ |
43 |
|
|
bound2 = AIR_CAST(float, 0.5 - 0.5*airErf((0.9-thick-r)/(bnd + 0.0001))); |
44 |
|
|
aniso = AIR_CAST(float, AIR_AFFINE(0.0, AIR_MIN(bound1, bound2), 1.0, |
45 |
|
|
mina, maxa)); |
46 |
|
|
|
47 |
|
|
ELL_3V_SET_TT(eval, float, |
48 |
|
|
AIR_LERP(aniso, 1.0/3.0, AIR_AFFINE(0.0, parm, 2.0, 1.0, 0.0)), |
49 |
|
|
AIR_LERP(aniso, 1.0/3.0, AIR_AFFINE(0.0, parm, 2.0, 0.0, 1.0)), |
50 |
|
|
AIR_LERP(aniso, 1.0/3.0, 0)); |
51 |
|
|
ELL_3V_SCALE(eval, evsc, eval); |
52 |
|
|
|
53 |
|
|
/* v0: looking down positive Z, points counter clockwise */ |
54 |
|
|
if (x || y) { |
55 |
|
|
ELL_3V_SET(evec + 3*0, y, -x, 0); |
56 |
|
|
ELL_3V_NORM_TT(evec + 3*0, float, evec + 3*0, norm); |
57 |
|
|
|
58 |
|
|
/* v1: points towards pole at positive Z */ |
59 |
|
|
ELL_3V_SET(tmp, -x, -y, -z); |
60 |
|
|
ELL_3V_NORM_TT(tmp, float, tmp, norm); |
61 |
|
|
ELL_3V_CROSS(evec + 3*1, tmp, evec + 3*0); |
62 |
|
|
|
63 |
|
|
/* v2: v0 x v1 */ |
64 |
|
|
ELL_3V_CROSS(evec + 3*2, evec + 3*0, evec + 3*1); |
65 |
|
|
} else { |
66 |
|
|
/* not optimal, but at least it won't show up in glyph visualizations */ |
67 |
|
|
ELL_3M_IDENTITY_SET(evec); |
68 |
|
|
} |
69 |
|
|
return; |
70 |
|
|
} |
71 |
|
|
|
72 |
|
|
void |
73 |
|
|
tend_satinTorusEigen(float *eval, float *evec, float x, float y, float z, |
74 |
|
|
float parm, float mina, float maxa, |
75 |
|
|
float thick, float bnd, float evsc) { |
76 |
|
|
float bound, R, r, norm, out[3], up[3], aniso; |
77 |
|
|
|
78 |
|
|
thick *= 2; |
79 |
|
|
R = AIR_CAST(float, sqrt(x*x + y*y)); |
80 |
|
|
r = AIR_CAST(float, sqrt((R-1)*(R-1) + z*z)); |
81 |
|
|
/* 1 on inside, 0 on outside */ |
82 |
|
|
bound = AIR_CAST(float, 0.5 - 0.5*airErf((r-thick)/(bnd + 0.0001))); |
83 |
|
|
aniso = AIR_CAST(float, AIR_AFFINE(0, bound, 1, mina, maxa)); |
84 |
|
|
|
85 |
|
|
ELL_3V_SET_TT(eval, float, |
86 |
|
|
AIR_LERP(aniso, 1.0/3.0, AIR_AFFINE(0.0, parm, 2.0, 1.0, 0.0)), |
87 |
|
|
AIR_LERP(aniso, 1.0/3.0, AIR_AFFINE(0.0, parm, 2.0, 0.0, 1.0)), |
88 |
|
|
AIR_LERP(aniso, 1.0/3.0, 0)); |
89 |
|
|
ELL_3V_SCALE(eval, evsc, eval); |
90 |
|
|
|
91 |
|
|
ELL_3V_SET(up, 0, 0, 1); |
92 |
|
|
|
93 |
|
|
if (x || y) { |
94 |
|
|
/* v0: looking down positive Z, points counter clockwise */ |
95 |
|
|
ELL_3V_SET(evec + 3*0, y, -x, 0); |
96 |
|
|
ELL_3V_NORM_TT(evec + 3*0, float, evec + 3*0, norm); |
97 |
|
|
|
98 |
|
|
/* v2: points into core of torus */ |
99 |
|
|
/* out: points away from (x,y)=(0,0) */ |
100 |
|
|
ELL_3V_SET(out, x, y, 0); |
101 |
|
|
ELL_3V_NORM_TT(out, float, out, norm); |
102 |
|
|
ELL_3V_SCALE_ADD2(evec + 3*2, -z, up, (1-R), out); |
103 |
|
|
ELL_3V_NORM_TT(evec + 3*2, float, evec + 3*2, norm); |
104 |
|
|
|
105 |
|
|
/* v1: looking at right half of cross-section, points counter clockwise */ |
106 |
|
|
ELL_3V_CROSS(evec + 3*1, evec + 3*0, evec + 3*2); |
107 |
|
|
} else { |
108 |
|
|
/* not optimal, but at least it won't show up in glyph visualizations */ |
109 |
|
|
ELL_3M_IDENTITY_SET(evec); |
110 |
|
|
} |
111 |
|
|
return; |
112 |
|
|
} |
113 |
|
|
|
114 |
|
|
int |
115 |
|
|
tend_satinGen(Nrrd *nout, float parm, float mina, float maxa, int wsize, |
116 |
|
|
float thick, float scaling, |
117 |
|
|
float bnd, float bndRm, float evsc, int torus) { |
118 |
|
|
static const char me[]="tend_satinGen"; |
119 |
|
|
char buff[AIR_STRLEN_SMALL]; |
120 |
|
|
Nrrd *nconf, *neval, *nevec; |
121 |
|
|
float *conf, *eval, *evec; |
122 |
|
|
size_t xi, yi, zi, size[3]; |
123 |
|
|
float x, y, z, min[3], max[3]; |
124 |
|
|
|
125 |
|
|
if (torus) { |
126 |
|
|
ELL_3V_SET(size, 2*wsize, 2*wsize, wsize); |
127 |
|
|
ELL_3V_SET(min, -2, -2, -1); |
128 |
|
|
ELL_3V_SET(max, 2, 2, 1); |
129 |
|
|
} else { |
130 |
|
|
ELL_3V_SET(size, wsize, wsize, wsize); |
131 |
|
|
ELL_3V_SET(min, -1, -1, -1); |
132 |
|
|
ELL_3V_SET(max, 1, 1, 1); |
133 |
|
|
} |
134 |
|
|
if (nrrdMaybeAlloc_va(nconf=nrrdNew(), nrrdTypeFloat, 3, |
135 |
|
|
size[0], size[1], size[2]) || |
136 |
|
|
nrrdMaybeAlloc_va(neval=nrrdNew(), nrrdTypeFloat, 4, |
137 |
|
|
AIR_CAST(size_t, 3), size[0], size[1], size[2]) || |
138 |
|
|
nrrdMaybeAlloc_va(nevec=nrrdNew(), nrrdTypeFloat, 4, |
139 |
|
|
AIR_CAST(size_t, 9), size[0], size[1], size[2])) { |
140 |
|
|
biffMovef(TEN, NRRD, "%s: trouble allocating temp nrrds", me); |
141 |
|
|
return 1; |
142 |
|
|
} |
143 |
|
|
|
144 |
|
|
conf = (float *)nconf->data; |
145 |
|
|
eval = (float *)neval->data; |
146 |
|
|
evec = (float *)nevec->data; |
147 |
|
|
for (zi=0; zi<size[2]; zi++) { |
148 |
|
|
z = AIR_CAST(float, AIR_AFFINE(0, zi, size[2]-1, min[2], max[2])); |
149 |
|
|
z /= scaling; |
150 |
|
|
for (yi=0; yi<size[1]; yi++) { |
151 |
|
|
y = AIR_CAST(float, AIR_AFFINE(0, yi, size[1]-1, min[1], max[1])); |
152 |
|
|
y /= scaling; |
153 |
|
|
for (xi=0; xi<size[0]; xi++) { |
154 |
|
|
x = AIR_CAST(float, AIR_AFFINE(0, xi, size[0]-1, min[0], max[0])); |
155 |
|
|
x /= scaling; |
156 |
|
|
*conf = 1.0; |
157 |
|
|
if (torus) { |
158 |
|
|
float aff; |
159 |
|
|
aff = AIR_CAST(float, AIR_AFFINE(0, yi, size[1]-1, 0, 1)); |
160 |
|
|
tend_satinTorusEigen(eval, evec, x, y, z, parm, |
161 |
|
|
mina, maxa, thick, bnd + bndRm*aff, evsc); |
162 |
|
|
} else { |
163 |
|
|
float aff; |
164 |
|
|
aff = AIR_CAST(float, AIR_AFFINE(0,yi, size[1]-1, 0, 1)); |
165 |
|
|
tend_satinSphereEigen(eval, evec, x, y, z, parm, |
166 |
|
|
mina, maxa, thick, bnd + bndRm*aff, evsc); |
167 |
|
|
} |
168 |
|
|
conf += 1; |
169 |
|
|
eval += 3; |
170 |
|
|
evec += 9; |
171 |
|
|
} |
172 |
|
|
} |
173 |
|
|
} |
174 |
|
|
|
175 |
|
|
if (tenMake(nout, nconf, neval, nevec)) { |
176 |
|
|
biffAddf(TEN, "%s: trouble generating output", me); |
177 |
|
|
return 1; |
178 |
|
|
} |
179 |
|
|
|
180 |
|
|
nrrdNuke(nconf); |
181 |
|
|
nrrdNuke(neval); |
182 |
|
|
nrrdNuke(nevec); |
183 |
|
|
nrrdAxisInfoSet_va(nout, nrrdAxisInfoLabel, "tensor", "x", "y", "z"); |
184 |
|
|
sprintf(buff, "satin(%g,%g,%g)", parm, mina, maxa); |
185 |
|
|
nout->content = airStrdup(buff); |
186 |
|
|
return 0; |
187 |
|
|
} |
188 |
|
|
|
189 |
|
|
int |
190 |
|
|
tend_satinMain(int argc, const char **argv, const char *me, |
191 |
|
|
hestParm *hparm) { |
192 |
|
|
int pret; |
193 |
|
2 |
hestOpt *hopt = NULL; |
194 |
|
1 |
char *perr, *err; |
195 |
|
|
airArray *mop; |
196 |
|
|
|
197 |
|
1 |
int wsize, torus; |
198 |
|
1 |
float parm, maxa, mina, thick, scaling, bnd, bndRm, evsc; |
199 |
|
|
Nrrd *nout; |
200 |
|
1 |
char *outS; |
201 |
|
|
gageShape *shape; |
202 |
|
1 |
double spd[4][4], orig[4]; |
203 |
|
|
|
204 |
|
1 |
hestOptAdd(&hopt, "t", "do torus", airTypeInt, 0, 0, &torus, NULL, |
205 |
|
|
"generate a torus dataset, instead of the default spherical"); |
206 |
|
1 |
hestOptAdd(&hopt, "p", "aniso parm", airTypeFloat, 1, 1, &parm, NULL, |
207 |
|
|
"anisotropy parameter. 0.0 for one direction of linear (along " |
208 |
|
|
"the equator for spheres, or along the larger circumference for " |
209 |
|
|
"toruses), 1.0 for planar, 2.0 for the other direction of linear " |
210 |
|
|
"(from pole to pole for spheres, or along the smaller " |
211 |
|
|
"circumference for toruses)"); |
212 |
|
1 |
hestOptAdd(&hopt, "max", "max ca1", airTypeFloat, 1, 1, &maxa, "1.0", |
213 |
|
|
"maximum anisotropy in dataset, according to the \"ca1\" " |
214 |
|
|
"anisotropy metric. \"1.0\" means " |
215 |
|
|
"completely linear or completely planar anisotropy"); |
216 |
|
1 |
hestOptAdd(&hopt, "min", "min ca1", airTypeFloat, 1, 1, &mina, "0.0", |
217 |
|
|
"minimum anisotropy in dataset"); |
218 |
|
1 |
hestOptAdd(&hopt, "b", "boundary", airTypeFloat, 1, 1, &bnd, "0.05", |
219 |
|
|
"parameter governing how fuzzy the boundary between high and " |
220 |
|
|
"low anisotropy is. Use \"-b 0\" for no fuzziness"); |
221 |
|
1 |
hestOptAdd(&hopt, "br", "ramp", airTypeFloat, 1, 1, &bndRm, "0.0", |
222 |
|
|
"how much to ramp upeffective \"b\" along Y axis. " |
223 |
|
|
"Use \"-b 0\" for no such ramping."); |
224 |
|
1 |
hestOptAdd(&hopt, "th", "thickness", airTypeFloat, 1, 1, &thick, "0.3", |
225 |
|
|
"parameter governing how thick region of high anisotropy is"); |
226 |
|
1 |
hestOptAdd(&hopt, "scl", "scaling", airTypeFloat, 1, 1, &scaling, "1.0", |
227 |
|
|
"scaling on size of sphere or torus within volume; lowering " |
228 |
|
|
"this below default 1.0 produces more background margin"); |
229 |
|
1 |
hestOptAdd(&hopt, "evsc", "eval scale", airTypeFloat, 1, 1, &evsc, "1.0", |
230 |
|
|
"scaling of eigenvalues"); |
231 |
|
1 |
hestOptAdd(&hopt, "s", "size", airTypeInt, 1, 1, &wsize, "32", |
232 |
|
|
"dimensions of output volume. For size N, the output is " |
233 |
|
|
"N\tx\tN\tx\tN for spheres, and 2N\tx\t2N\tx\tN for toruses"); |
234 |
|
1 |
hestOptAdd(&hopt, "o", "nout", airTypeString, 1, 1, &outS, "-", |
235 |
|
|
"output filename"); |
236 |
|
|
|
237 |
|
1 |
mop = airMopNew(); |
238 |
|
1 |
airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways); |
239 |
✓✗ |
2 |
USAGE(_tend_satinInfoL); |
240 |
|
|
JUSTPARSE(); |
241 |
|
|
airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways); |
242 |
|
|
|
243 |
|
|
nout = nrrdNew(); |
244 |
|
|
airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways); |
245 |
|
|
|
246 |
|
|
if (tend_satinGen(nout, parm, mina, maxa, wsize, thick, scaling, |
247 |
|
|
bnd, bndRm, evsc, torus)) { |
248 |
|
|
airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways); |
249 |
|
|
fprintf(stderr, "%s: trouble making volume:\n%s\n", me, err); |
250 |
|
|
airMopError(mop); return 1; |
251 |
|
|
} |
252 |
|
|
|
253 |
|
|
/* use gageShape to determine orientation info */ |
254 |
|
|
nrrdAxisInfoSet_va(nout, nrrdAxisInfoCenter, |
255 |
|
|
nrrdCenterUnknown, nrrdCenterCell, |
256 |
|
|
nrrdCenterCell, nrrdCenterCell); |
257 |
|
|
shape = gageShapeNew(); |
258 |
|
|
airMopAdd(mop, shape, (airMopper)gageShapeNix, airMopAlways); |
259 |
|
|
/* this is a weird mix of new and legacy code. At some point |
260 |
|
|
prior to Wed May 27 19:23:55 CDT 2009, it was okay to pass |
261 |
|
|
in a volume to gageShapeSet that had absolutely no notion |
262 |
|
|
of spacing or orientation. Then gageShapeSet was used to |
263 |
|
|
get a plausible set of space directions and space origin. |
264 |
|
|
Now, we're setting some spacings, so that gageShapeSet can |
265 |
|
|
do its thing, then (below) nan-ing out those spacings so |
266 |
|
|
that the nrrd is self-consistent */ |
267 |
|
|
nout->axis[1].spacing = 1.0; |
268 |
|
|
nout->axis[2].spacing = 1.0; |
269 |
|
|
nout->axis[3].spacing = 1.0; |
270 |
|
|
if (gageShapeSet(shape, nout, tenGageKind->baseDim)) { |
271 |
|
|
airMopAdd(mop, err=biffGetDone(GAGE), airFree, airMopAlways); |
272 |
|
|
fprintf(stderr, "%s: trouble doing shape:\n%s\n", me, err); |
273 |
|
|
airMopError(mop); return 1; |
274 |
|
|
} |
275 |
|
|
/* the ItoW is a 4x4 matrix, but |
276 |
|
|
we really only care about the first three rows */ |
277 |
|
|
ELL_4V_SET(spd[0], AIR_NAN, AIR_NAN, AIR_NAN, AIR_NAN); |
278 |
|
|
ELL_4MV_COL0_GET(spd[1], shape->ItoW); ELL_4V_SCALE(spd[1], 32, spd[1]); |
279 |
|
|
ELL_4MV_COL1_GET(spd[2], shape->ItoW); ELL_4V_SCALE(spd[2], 32, spd[2]); |
280 |
|
|
ELL_4MV_COL2_GET(spd[3], shape->ItoW); ELL_4V_SCALE(spd[3], 32, spd[3]); |
281 |
|
|
ELL_4MV_COL3_GET(orig, shape->ItoW); ELL_4V_SCALE(orig, 32, orig); |
282 |
|
|
nrrdSpaceSet(nout, nrrdSpaceRightAnteriorSuperior); |
283 |
|
|
nrrdSpaceOriginSet(nout, orig); |
284 |
|
|
nrrdAxisInfoSet_va(nout, nrrdAxisInfoSpaceDirection, |
285 |
|
|
spd[0], spd[1], spd[2], spd[3]); |
286 |
|
|
nout->axis[1].spacing = AIR_NAN; |
287 |
|
|
nout->axis[2].spacing = AIR_NAN; |
288 |
|
|
nout->axis[3].spacing = AIR_NAN; |
289 |
|
|
nrrdAxisInfoSet_va(nout, nrrdAxisInfoKind, |
290 |
|
|
nrrdKind3DMaskedSymMatrix, nrrdKindSpace, |
291 |
|
|
nrrdKindSpace, nrrdKindSpace); |
292 |
|
|
nout->measurementFrame[0][0] = 1; |
293 |
|
|
nout->measurementFrame[1][0] = 0; |
294 |
|
|
nout->measurementFrame[2][0] = 0; |
295 |
|
|
nout->measurementFrame[0][1] = 0; |
296 |
|
|
nout->measurementFrame[1][1] = 1; |
297 |
|
|
nout->measurementFrame[2][1] = 0; |
298 |
|
|
nout->measurementFrame[0][2] = 0; |
299 |
|
|
nout->measurementFrame[1][2] = 0; |
300 |
|
|
nout->measurementFrame[2][2] = 1; |
301 |
|
|
|
302 |
|
|
if (nrrdSave(outS, nout, NULL)) { |
303 |
|
|
airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways); |
304 |
|
|
fprintf(stderr, "%s: trouble writing:\n%s\n", me, err); |
305 |
|
|
airMopError(mop); return 1; |
306 |
|
|
} |
307 |
|
|
|
308 |
|
|
airMopOkay(mop); |
309 |
|
|
return 0; |
310 |
|
1 |
} |
311 |
|
|
/* TEND_CMD(satin, INFO); */ |
312 |
|
|
unrrduCmd tend_satinCmd = { "satin", INFO, tend_satinMain, AIR_FALSE }; |