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 "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
tend_satinSphereEigen(float * eval,float * evec,float x,float y,float z,float parm,float mina,float maxa,float thick,float bnd,float evsc)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
tend_satinTorusEigen(float * eval,float * evec,float x,float y,float z,float parm,float mina,float maxa,float thick,float bnd,float evsc)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
tend_satinGen(Nrrd * nout,float parm,float mina,float maxa,int wsize,float thick,float bnd,float bndRm,float evsc,int torus)115 tend_satinGen(Nrrd *nout, float parm, float mina, float maxa, int wsize,
116 float thick, float bnd, float bndRm, float evsc, int torus) {
117 static const char me[]="tend_satinGen";
118 char buff[AIR_STRLEN_SMALL];
119 Nrrd *nconf, *neval, *nevec;
120 float *conf, *eval, *evec;
121 size_t xi, yi, zi, size[3];
122 float x, y, z, min[3], max[3];
123
124 if (torus) {
125 ELL_3V_SET(size, 2*wsize, 2*wsize, wsize);
126 ELL_3V_SET(min, -2, -2, -1);
127 ELL_3V_SET(max, 2, 2, 1);
128 } else {
129 ELL_3V_SET(size, wsize, wsize, wsize);
130 ELL_3V_SET(min, -1, -1, -1);
131 ELL_3V_SET(max, 1, 1, 1);
132 }
133 if (nrrdMaybeAlloc_va(nconf=nrrdNew(), nrrdTypeFloat, 3,
134 size[0], size[1], size[2]) ||
135 nrrdMaybeAlloc_va(neval=nrrdNew(), nrrdTypeFloat, 4,
136 AIR_CAST(size_t, 3), size[0], size[1], size[2]) ||
137 nrrdMaybeAlloc_va(nevec=nrrdNew(), nrrdTypeFloat, 4,
138 AIR_CAST(size_t, 9), size[0], size[1], size[2])) {
139 biffMovef(TEN, NRRD, "%s: trouble allocating temp nrrds", me);
140 return 1;
141 }
142
143 conf = (float *)nconf->data;
144 eval = (float *)neval->data;
145 evec = (float *)nevec->data;
146 for (zi=0; zi<size[2]; zi++) {
147 z = AIR_CAST(float, AIR_AFFINE(0, zi, size[2]-1, min[2], max[2]));
148 for (yi=0; yi<size[1]; yi++) {
149 y = AIR_CAST(float, AIR_AFFINE(0, yi, size[1]-1, min[1], max[1]));
150 for (xi=0; xi<size[0]; xi++) {
151 x = AIR_CAST(float, AIR_AFFINE(0, xi, size[0]-1, min[0], max[0]));
152 *conf = 1.0;
153 if (torus) {
154 float aff;
155 aff = AIR_CAST(float, AIR_AFFINE(0, yi, size[1]-1, 0, 1));
156 tend_satinTorusEigen(eval, evec, x, y, z, parm,
157 mina, maxa, thick, bnd + bndRm*aff, evsc);
158 } else {
159 float aff;
160 aff = AIR_CAST(float, AIR_AFFINE(0,yi, size[1]-1, 0, 1));
161 tend_satinSphereEigen(eval, evec, x, y, z, parm,
162 mina, maxa, thick, bnd + bndRm*aff, evsc);
163 }
164 conf += 1;
165 eval += 3;
166 evec += 9;
167 }
168 }
169 }
170
171 if (tenMake(nout, nconf, neval, nevec)) {
172 biffAddf(TEN, "%s: trouble generating output", me);
173 return 1;
174 }
175
176 nrrdNuke(nconf);
177 nrrdNuke(neval);
178 nrrdNuke(nevec);
179 nrrdAxisInfoSet_va(nout, nrrdAxisInfoLabel, "tensor", "x", "y", "z");
180 sprintf(buff, "satin(%g,%g,%g)", parm, mina, maxa);
181 nout->content = airStrdup(buff);
182 return 0;
183 }
184
185 int
tend_satinMain(int argc,const char ** argv,const char * me,hestParm * hparm)186 tend_satinMain(int argc, const char **argv, const char *me,
187 hestParm *hparm) {
188 int pret;
189 hestOpt *hopt = NULL;
190 char *perr, *err;
191 airArray *mop;
192
193 int wsize, torus;
194 float parm, maxa, mina, thick, bnd, bndRm, evsc;
195 Nrrd *nout;
196 char *outS;
197 gageShape *shape;
198 double spd[4][4], orig[4];
199
200 hestOptAdd(&hopt, "t", "do torus", airTypeInt, 0, 0, &torus, NULL,
201 "generate a torus dataset, instead of the default spherical");
202 hestOptAdd(&hopt, "p", "aniso parm", airTypeFloat, 1, 1, &parm, NULL,
203 "anisotropy parameter. 0.0 for one direction of linear (along "
204 "the equator for spheres, or along the larger circumference for "
205 "toruses), 1.0 for planar, 2.0 for the other direction of linear "
206 "(from pole to pole for spheres, or along the smaller "
207 "circumference for toruses)");
208 hestOptAdd(&hopt, "max", "max ca1", airTypeFloat, 1, 1, &maxa, "1.0",
209 "maximum anisotropy in dataset, according to the \"ca1\" "
210 "anisotropy metric. \"1.0\" means "
211 "completely linear or completely planar anisotropy");
212 hestOptAdd(&hopt, "min", "min ca1", airTypeFloat, 1, 1, &mina, "0.0",
213 "minimum anisotropy in dataset");
214 hestOptAdd(&hopt, "b", "boundary", airTypeFloat, 1, 1, &bnd, "0.05",
215 "parameter governing how fuzzy the boundary between high and "
216 "low anisotropy is. Use \"-b 0\" for no fuzziness");
217 hestOptAdd(&hopt, "br", "ramp", airTypeFloat, 1, 1, &bndRm, "0.0",
218 "how much to ramp upeffective \"b\" along Y axis. "
219 "Use \"-b 0\" for no such ramping.");
220 hestOptAdd(&hopt, "th", "thickness", airTypeFloat, 1, 1, &thick, "0.3",
221 "parameter governing how thick region of high anisotropy is");
222 hestOptAdd(&hopt, "evsc", "eval scale", airTypeFloat, 1, 1, &evsc, "1.0",
223 "scaling of eigenvalues");
224 hestOptAdd(&hopt, "s", "size", airTypeInt, 1, 1, &wsize, "32",
225 "dimensions of output volume. For size N, the output is "
226 "N\tx\tN\tx\tN for spheres, and 2N\tx\t2N\tx\tN for toruses");
227 hestOptAdd(&hopt, "o", "nout", airTypeString, 1, 1, &outS, "-",
228 "output filename");
229
230 mop = airMopNew();
231 airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
232 USAGE(_tend_satinInfoL);
233 JUSTPARSE();
234 airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
235
236 nout = nrrdNew();
237 airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);
238
239 if (tend_satinGen(nout, parm, mina, maxa, wsize, thick,
240 bnd, bndRm, evsc, torus)) {
241 airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways);
242 fprintf(stderr, "%s: trouble making volume:\n%s\n", me, err);
243 airMopError(mop); return 1;
244 }
245
246 /* use gageShape to determine orientation info */
247 nrrdAxisInfoSet_va(nout, nrrdAxisInfoCenter,
248 nrrdCenterUnknown, nrrdCenterCell,
249 nrrdCenterCell, nrrdCenterCell);
250 shape = gageShapeNew();
251 airMopAdd(mop, shape, (airMopper)gageShapeNix, airMopAlways);
252 /* this is a weird mix of new and legacy code. At some point
253 prior to Wed May 27 19:23:55 CDT 2009, it was okay to pass
254 in a volume to gageShapeSet that had absolutely no notion
255 of spacing or orientation. Then gageShapeSet was used to
256 get a plausible set of space directions and space origin.
257 Now, we're setting some spacings, so that gageShapeSet can
258 do its thing, then (below) nan-ing out those spacings so
259 that the nrrd is self-consistent */
260 nout->axis[1].spacing = 1.0;
261 nout->axis[2].spacing = 1.0;
262 nout->axis[3].spacing = 1.0;
263 if (gageShapeSet(shape, nout, tenGageKind->baseDim)) {
264 airMopAdd(mop, err=biffGetDone(GAGE), airFree, airMopAlways);
265 fprintf(stderr, "%s: trouble doing shape:\n%s\n", me, err);
266 airMopError(mop); return 1;
267 }
268 /* the ItoW is a 4x4 matrix, but
269 we really only care about the first three rows */
270 ELL_4V_SET(spd[0], AIR_NAN, AIR_NAN, AIR_NAN, AIR_NAN);
271 ELL_4MV_COL0_GET(spd[1], shape->ItoW); ELL_4V_SCALE(spd[1], 32, spd[1]);
272 ELL_4MV_COL1_GET(spd[2], shape->ItoW); ELL_4V_SCALE(spd[2], 32, spd[2]);
273 ELL_4MV_COL2_GET(spd[3], shape->ItoW); ELL_4V_SCALE(spd[3], 32, spd[3]);
274 ELL_4MV_COL3_GET(orig, shape->ItoW); ELL_4V_SCALE(orig, 32, orig);
275 nrrdSpaceSet(nout, nrrdSpaceRightAnteriorSuperior);
276 nrrdSpaceOriginSet(nout, orig);
277 nrrdAxisInfoSet_va(nout, nrrdAxisInfoSpaceDirection,
278 spd[0], spd[1], spd[2], spd[3]);
279 nout->axis[1].spacing = AIR_NAN;
280 nout->axis[2].spacing = AIR_NAN;
281 nout->axis[3].spacing = AIR_NAN;
282 nrrdAxisInfoSet_va(nout, nrrdAxisInfoKind,
283 nrrdKind3DMaskedSymMatrix, nrrdKindSpace,
284 nrrdKindSpace, nrrdKindSpace);
285 nout->measurementFrame[0][0] = 1;
286 nout->measurementFrame[1][0] = 0;
287 nout->measurementFrame[2][0] = 0;
288 nout->measurementFrame[0][1] = 0;
289 nout->measurementFrame[1][1] = 1;
290 nout->measurementFrame[2][1] = 0;
291 nout->measurementFrame[0][2] = 0;
292 nout->measurementFrame[1][2] = 0;
293 nout->measurementFrame[2][2] = 1;
294
295 if (nrrdSave(outS, nout, NULL)) {
296 airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways);
297 fprintf(stderr, "%s: trouble writing:\n%s\n", me, err);
298 airMopError(mop); return 1;
299 }
300
301 airMopOkay(mop);
302 return 0;
303 }
304 /* TEND_CMD(satin, INFO); */
305 unrrduCmd tend_satinCmd = { "satin", INFO, tend_satinMain };
306