GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: src/ten/tendSatin.c Lines: 21 145 14.5 %
Date: 2017-05-26 Branches: 1 36 2.8 %

Line Branch Exec Source
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 };