GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: src/ten/glyph.c Lines: 94 516 18.2 %
Date: 2017-05-26 Branches: 27 285 9.5 %

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
tenGlyphParm *
28
tenGlyphParmNew() {
29
  tenGlyphParm *parm;
30
31
2
  parm = (tenGlyphParm *)calloc(1, sizeof(tenGlyphParm));
32
1
  if (parm) {
33
1
    parm->verbose = 0;
34
1
    parm->nmask = NULL;
35
1
    parm->anisoType = tenAnisoUnknown;
36
1
    parm->onlyPositive = AIR_TRUE;
37
1
    parm->confThresh = AIR_NAN;
38
1
    parm->anisoThresh = AIR_NAN;
39
1
    parm->maskThresh = AIR_NAN;
40
41
1
    parm->glyphType = tenGlyphTypeUnknown;
42
1
    parm->facetRes = 10;
43
1
    parm->glyphScale = 1.0;
44
1
    parm->sqdSharp = 3.0;
45
1
    ELL_5V_SET(parm->edgeWidth, 0.0f, 0.0f, 0.4f, 0.2f, 0.1f);
46
47
1
    parm->colEvec = 0;  /* first */
48
1
    parm->colMaxSat = 1;
49
1
    parm->colGamma = 1;
50
1
    parm->colIsoGray = 1;
51
1
    parm->colAnisoType = tenAnisoUnknown;
52
1
    parm->colAnisoModulate = 0;
53
1
    ELL_4V_SET(parm->ADSP, 0, 1, 0, 30);
54
55
1
    parm->doSlice = AIR_FALSE;
56
1
    parm->sliceAxis = 0;
57
1
    parm->slicePos = 0;
58
1
    parm->sliceAnisoType = tenAnisoUnknown;
59
1
    parm->sliceOffset = 0.0;
60
1
    parm->sliceBias = 0.05f;
61
1
    parm->sliceGamma = 1.0;
62
1
  }
63
1
  return parm;
64
}
65
66
tenGlyphParm *
67
tenGlyphParmNix(tenGlyphParm *parm) {
68
69
2
  airFree(parm);
70
1
  return NULL;
71
}
72
73
int
74
tenGlyphParmCheck(tenGlyphParm *parm,
75
                  const Nrrd *nten, const Nrrd *npos, const Nrrd *nslc) {
76
  static const char me[]="tenGlyphParmCheck";
77
  int duh;
78
  size_t tenSize[3];
79
  char stmp[5][AIR_STRLEN_SMALL];
80
81
  if (!(parm && nten)) {
82
    biffAddf(TEN, "%s: got NULL pointer", me);
83
    return 1;
84
  }
85
  if (airEnumValCheck(tenAniso, parm->anisoType)) {
86
    biffAddf(TEN, "%s: unset (or invalid) anisoType (%d)",
87
             me, parm->anisoType);
88
    return 1;
89
  }
90
  if (airEnumValCheck(tenAniso, parm->colAnisoType)) {
91
    biffAddf(TEN, "%s: unset (or invalid) colAnisoType (%d)",
92
             me, parm->colAnisoType);
93
    return 1;
94
  }
95
  if (!( parm->facetRes >= 3 )) {
96
    biffAddf(TEN, "%s: facet resolution %d not >= 3", me, parm->facetRes);
97
    return 1;
98
  }
99
  if (!( AIR_IN_OP(tenGlyphTypeUnknown, parm->glyphType,
100
                   tenGlyphTypeLast) )) {
101
    biffAddf(TEN, "%s: unset (or invalid) glyphType (%d)",
102
             me, parm->glyphType);
103
    return 1;
104
  }
105
  if (!( parm->glyphScale > 0)) {
106
    biffAddf(TEN, "%s: glyphScale must be > 0 (not %g)", me, parm->glyphScale);
107
    return 1;
108
  }
109
  if (parm->nmask) {
110
    if (npos) {
111
      biffAddf(TEN, "%s: can't do masking with explicit coordinate list", me);
112
      return 1;
113
    }
114
    if (!( 3 == parm->nmask->dim
115
           && parm->nmask->axis[0].size == nten->axis[1].size
116
           && parm->nmask->axis[1].size == nten->axis[2].size
117
           && parm->nmask->axis[2].size == nten->axis[3].size )) {
118
      biffAddf(TEN, "%s: mask isn't 3-D or doesn't have sizes (%s,%s,%s)", me,
119
               airSprintSize_t(stmp[0], nten->axis[1].size),
120
               airSprintSize_t(stmp[1], nten->axis[2].size),
121
               airSprintSize_t(stmp[2], nten->axis[3].size));
122
      return 1;
123
    }
124
    if (!(AIR_EXISTS(parm->maskThresh))) {
125
      biffAddf(TEN, "%s: maskThresh hasn't been set", me);
126
      return 1;
127
    }
128
  }
129
  if (!( AIR_EXISTS(parm->anisoThresh)
130
         && AIR_EXISTS(parm->confThresh) )) {
131
    biffAddf(TEN, "%s: anisoThresh and confThresh haven't both been set", me);
132
    return 1;
133
  }
134
  if (parm->doSlice) {
135
    if (npos) {
136
      biffAddf(TEN, "%s: can't do slice with explicit coordinate list", me);
137
      return 1;
138
    }
139
    if (!( parm->sliceAxis <=2 )) {
140
      biffAddf(TEN, "%s: slice axis %d invalid", me, parm->sliceAxis);
141
      return 1;
142
    }
143
    if (!( parm->slicePos < nten->axis[1+parm->sliceAxis].size )) {
144
      biffAddf(TEN, "%s: slice pos %s not in valid range [0..%s]", me,
145
               airSprintSize_t(stmp[0], parm->slicePos),
146
               airSprintSize_t(stmp[1], nten->axis[1+parm->sliceAxis].size-1));
147
      return 1;
148
    }
149
    if (nslc) {
150
      if (2 != nslc->dim) {
151
        biffAddf(TEN, "%s: explicit slice must be 2-D (not %d)",
152
                 me, nslc->dim);
153
        return 1;
154
      }
155
      tenSize[0] = nten->axis[1].size;
156
      tenSize[1] = nten->axis[2].size;
157
      tenSize[2] = nten->axis[3].size;
158
      for (duh=parm->sliceAxis; duh<2; duh++) {
159
        tenSize[duh] = tenSize[duh+1];
160
      }
161
      if (!( tenSize[0] == nslc->axis[0].size
162
             && tenSize[1] == nslc->axis[1].size )) {
163
        biffAddf(TEN, "%s: axis %u slice of %sx%sx%s volume != %sx%s", me,
164
                 parm->sliceAxis,
165
                 airSprintSize_t(stmp[0], nten->axis[1].size),
166
                 airSprintSize_t(stmp[1], nten->axis[2].size),
167
                 airSprintSize_t(stmp[2], nten->axis[3].size),
168
                 airSprintSize_t(stmp[3], nslc->axis[0].size),
169
                 airSprintSize_t(stmp[4], nslc->axis[1].size));
170
        return 1;
171
      }
172
    } else {
173
      if (airEnumValCheck(tenAniso, parm->sliceAnisoType)) {
174
        biffAddf(TEN, "%s: unset (or invalid) sliceAnisoType (%d)",
175
                 me, parm->sliceAnisoType);
176
        return 1;
177
      }
178
    }
179
  }
180
  return 0;
181
}
182
183
int
184
tenGlyphGen(limnObject *glyphsLimn, echoScene *glyphsEcho,
185
            tenGlyphParm *parm,
186
            const Nrrd *nten, const Nrrd *npos, const Nrrd *nslc) {
187
  static const char me[]="tenGlyphGen";
188
  gageShape *shape;
189
  airArray *mop;
190
  float *tdata, eval[3], evec[9], *cvec, rotEvec[9], mA_f[16],
191
    absEval[3], glyphScl[3];
192
  double pI[3], pW[3], cl, cp, sRot[16], mA[16], mB[16], msFr[9], tmpvec[3],
193
    R, G, B, qA, qB, qC, glyphAniso, sliceGray;
194
  unsigned int duh;
195
  int slcCoord[3], idx, glyphIdx, axis, numGlyphs,
196
    svRGBAfl=AIR_FALSE;
197
  limnLook *look; int lookIdx;
198
  echoObject *eglyph, *inst, *list=NULL, *split, *esquare;
199
  echoPos_t eM[16], originOffset[3], edge0[3], edge1[3];
200
  char stmp[AIR_STRLEN_SMALL];
201
  /*
202
  int eret;
203
  double tmp1[3], tmp2[3];
204
  */
205
206
  if (!( (glyphsLimn || glyphsEcho) && nten && parm)) {
207
    biffAddf(TEN, "%s: got NULL pointer", me);
208
    return 1;
209
  }
210
  mop = airMopNew();
211
  shape = gageShapeNew();
212
  shape->defaultCenter = nrrdCenterCell;
213
  airMopAdd(mop, shape, (airMopper)gageShapeNix, airMopAlways);
214
  if (npos) {
215
    if (!( 2 == nten->dim && 7 == nten->axis[0].size )) {
216
      biffAddf(TEN, "%s: nten isn't 2-D 7-by-N array", me);
217
      airMopError(mop); return 1;
218
    }
219
    if (!( 2 == npos->dim && 3 == npos->axis[0].size
220
           && nten->axis[1].size == npos->axis[1].size )) {
221
      biffAddf(TEN, "%s: npos isn't 2-D 3-by-%s array", me,
222
               airSprintSize_t(stmp, nten->axis[1].size));
223
      airMopError(mop); return 1;
224
    }
225
    if (!( nrrdTypeFloat == nten->type && nrrdTypeFloat == npos->type )) {
226
      biffAddf(TEN, "%s: nten and npos must be %s, not %s and %s", me,
227
               airEnumStr(nrrdType, nrrdTypeFloat),
228
               airEnumStr(nrrdType, nten->type),
229
               airEnumStr(nrrdType, npos->type));
230
      airMopError(mop); return 1;
231
    }
232
  } else {
233
    if (tenTensorCheck(nten, nrrdTypeFloat, AIR_TRUE, AIR_TRUE)) {
234
      biffAddf(TEN, "%s: didn't get a valid DT volume", me);
235
      airMopError(mop); return 1;
236
    }
237
  }
238
  if (tenGlyphParmCheck(parm, nten, npos, nslc)) {
239
    biffAddf(TEN, "%s: trouble", me);
240
    airMopError(mop); return 1;
241
  }
242
  if (!npos) {
243
    if (gageShapeSet(shape, nten, tenGageKind->baseDim)) {
244
      biffMovef(TEN, GAGE, "%s: trouble", me);
245
      airMopError(mop); return 1;
246
    }
247
  }
248
  if (parm->doSlice) {
249
    ELL_3V_COPY(edge0, shape->spacing);
250
    ELL_3V_COPY(edge1, shape->spacing);
251
    edge0[parm->sliceAxis] = edge1[parm->sliceAxis] = 0.0;
252
    switch(parm->sliceAxis) {
253
    case 0:
254
      edge0[1] = edge1[2] = 0;
255
      ELL_4M_ROTATE_Y_SET(sRot, AIR_PI/2);
256
      break;
257
    case 1:
258
      edge0[0] = edge1[2] = 0;
259
      ELL_4M_ROTATE_X_SET(sRot, AIR_PI/2);
260
      break;
261
    case 2: default:
262
      edge0[0] = edge1[1] = 0;
263
      ELL_4M_IDENTITY_SET(sRot);
264
      break;
265
    }
266
    ELL_3V_COPY(originOffset, shape->spacing);
267
    ELL_3V_SCALE(originOffset, -0.5, originOffset);
268
    originOffset[parm->sliceAxis] *= -2*parm->sliceOffset;
269
  }
270
  if (glyphsLimn) {
271
    /* create limnLooks for diffuse and ambient-only shading */
272
    /* ??? */
273
    /* hack: save old value of setVertexRGBAFromLook, and set to true */
274
    svRGBAfl = glyphsLimn->setVertexRGBAFromLook;
275
    glyphsLimn->setVertexRGBAFromLook = AIR_TRUE;
276
  }
277
  if (glyphsEcho) {
278
    list = echoObjectNew(glyphsEcho, echoTypeList);
279
  }
280
  if (npos) {
281
    numGlyphs = AIR_UINT(nten->axis[1].size);
282
  } else {
283
    numGlyphs = shape->size[0] * shape->size[1] * shape->size[2];
284
  }
285
  /* find measurement frame transform */
286
  if (3 == nten->spaceDim
287
      && AIR_EXISTS(nten->measurementFrame[0][0])) {
288
    /*     msFr        nten->measurementFrame
289
    **   0  1  2      [0][0]   [1][0]   [2][0]
290
    **   3  4  5      [0][1]   [1][1]   [2][1]
291
    **   6  7  8      [0][2]   [1][2]   [2][2]
292
    */
293
    msFr[0] = nten->measurementFrame[0][0];
294
    msFr[3] = nten->measurementFrame[0][1];
295
    msFr[6] = nten->measurementFrame[0][2];
296
    msFr[1] = nten->measurementFrame[1][0];
297
    msFr[4] = nten->measurementFrame[1][1];
298
    msFr[7] = nten->measurementFrame[1][2];
299
    msFr[2] = nten->measurementFrame[2][0];
300
    msFr[5] = nten->measurementFrame[2][1];
301
    msFr[8] = nten->measurementFrame[2][2];
302
  } else {
303
    ELL_3M_IDENTITY_SET(msFr);
304
  }
305
  for (idx=0; idx<numGlyphs; idx++) {
306
    tdata = (float*)(nten->data) + 7*idx;
307
    if (parm->verbose >= 2) {
308
      fprintf(stderr, "%s: glyph %d/%d: hello %g    %g %g %g %g %g %g\n",
309
              me, idx, numGlyphs, tdata[0],
310
              tdata[1], tdata[2], tdata[3],
311
              tdata[4], tdata[5], tdata[6]);
312
    }
313
    if (!( TEN_T_EXISTS(tdata) )) {
314
      /* there's nothing we can do here */
315
      if (parm->verbose >= 2) {
316
        fprintf(stderr, "%s: glyph %d/%d: non-existent data\n",
317
                me, idx, numGlyphs);
318
      }
319
      continue;
320
    }
321
    if (npos) {
322
      ELL_3V_COPY(pW, (float*)(npos->data) + 3*idx);
323
      if (!( AIR_EXISTS(pW[0]) && AIR_EXISTS(pW[1]) && AIR_EXISTS(pW[2]) )) {
324
        /* position doesn't exist- perhaps because its from the push
325
           library, which might kill points by setting coords to nan */
326
        continue;
327
      }
328
    } else {
329
      NRRD_COORD_GEN(pI, shape->size, 3, idx);
330
      /* this does take into account full orientation */
331
      gageShapeItoW(shape, pW, pI);
332
      if (parm->nmask) {
333
        if (!( nrrdFLookup[parm->nmask->type](parm->nmask->data, idx)
334
               >= parm->maskThresh )) {
335
          if (parm->verbose >= 2) {
336
            fprintf(stderr, "%s: glyph %d/%d: doesn't meet mask thresh\n",
337
                    me, idx, numGlyphs);
338
          }
339
          continue;
340
        }
341
      }
342
    }
343
    tenEigensolve_f(eval, evec, tdata);
344
    /* transform eigenvectors by measurement frame */
345
    ELL_3MV_MUL(tmpvec, msFr, evec + 0);
346
    ELL_3V_COPY_TT(evec + 0, float, tmpvec);
347
    ELL_3MV_MUL(tmpvec, msFr, evec + 3);
348
    ELL_3V_COPY_TT(evec + 3, float, tmpvec);
349
    ELL_3MV_MUL(tmpvec, msFr, evec + 6);
350
    ELL_3V_COPY_TT(evec + 6, float, tmpvec);
351
    ELL_3V_CROSS(tmpvec, evec + 0, evec + 3);
352
    if (0 > ELL_3V_DOT(tmpvec, evec + 6)) {
353
      ELL_3V_SCALE(evec + 6, -1, evec + 6);
354
    }
355
    ELL_3M_TRANSPOSE(rotEvec, evec);
356
    if (parm->doSlice
357
        && pI[parm->sliceAxis] == parm->slicePos) {
358
      /* set sliceGray */
359
      if (nslc) {
360
        /* we aren't masked by confidence, as anisotropy slice is */
361
        for (duh=0; duh<parm->sliceAxis; duh++) {
362
          slcCoord[duh] = (int)(pI[duh]);
363
        }
364
        for (duh=duh<parm->sliceAxis; duh<2; duh++) {
365
          slcCoord[duh] = (int)(pI[duh+1]);
366
        }
367
        /* HEY: GLK has no idea what's going here */
368
        slcCoord[0] = (int)(pI[0]);
369
        slcCoord[1] = (int)(pI[1]);
370
        slcCoord[2] = (int)(pI[2]);
371
        sliceGray =
372
          nrrdFLookup[nslc->type](nslc->data, slcCoord[0]
373
                                  + nslc->axis[0].size*slcCoord[1]);
374
      } else {
375
        if (!( tdata[0] >= parm->confThresh )) {
376
          if (parm->verbose >= 2) {
377
            fprintf(stderr, "%s: glyph %d/%d (slice): conf %g < thresh %g\n",
378
                    me, idx, numGlyphs, tdata[0], parm->confThresh);
379
          }
380
          continue;
381
        }
382
        sliceGray = tenAnisoEval_f(eval, parm->sliceAnisoType);
383
      }
384
      if (parm->sliceGamma > 0) {
385
        sliceGray = AIR_AFFINE(0, sliceGray, 1, parm->sliceBias, 1);
386
        sliceGray = pow(sliceGray, 1.0/parm->sliceGamma);
387
      } else {
388
        sliceGray = AIR_AFFINE(0, sliceGray, 1, 0, 1-parm->sliceBias);
389
        sliceGray = 1.0 - pow(sliceGray, -1.0/parm->sliceGamma);
390
      }
391
      /* make slice contribution */
392
      /* HEY: this is *NOT* aware of shape->fromOrientation */
393
      if (glyphsLimn) {
394
        lookIdx = limnObjectLookAdd(glyphsLimn);
395
        look = glyphsLimn->look + lookIdx;
396
        ELL_4V_SET_TT(look->rgba, float, sliceGray, sliceGray, sliceGray, 1);
397
        ELL_3V_SET(look->kads, 1, 0, 0);
398
        look->spow = 0;
399
        glyphIdx = limnObjectSquareAdd(glyphsLimn, lookIdx);
400
        ELL_4M_IDENTITY_SET(mA);
401
        ell_4m_post_mul_d(mA, sRot);
402
        if (!npos) {
403
          ELL_4M_SCALE_SET(mB,
404
                           shape->spacing[0],
405
                           shape->spacing[1],
406
                           shape->spacing[2]);
407
        }
408
        ell_4m_post_mul_d(mA, mB);
409
        ELL_4M_TRANSLATE_SET(mB, pW[0], pW[1], pW[2]);
410
        ell_4m_post_mul_d(mA, mB);
411
        ELL_4M_TRANSLATE_SET(mB,
412
                             originOffset[0],
413
                             originOffset[1],
414
                             originOffset[2]);
415
        ell_4m_post_mul_d(mA, mB);
416
        ELL_4M_COPY_TT(mA_f, float, mA);
417
        limnObjectPartTransform(glyphsLimn, glyphIdx, mA_f);
418
      }
419
      if (glyphsEcho) {
420
        esquare = echoObjectNew(glyphsEcho,echoTypeRectangle);
421
        ELL_3V_ADD2(((echoRectangle*)esquare)->origin, pW, originOffset);
422
        ELL_3V_COPY(((echoRectangle*)esquare)->edge0, edge0);
423
        ELL_3V_COPY(((echoRectangle*)esquare)->edge1, edge1);
424
        echoColorSet(esquare,
425
                     AIR_CAST(echoCol_t, sliceGray),
426
                     AIR_CAST(echoCol_t, sliceGray),
427
                     AIR_CAST(echoCol_t, sliceGray), 1);
428
        /* this is pretty arbitrary- but I want shadows to have some effect.
429
           Previously, the material was all ambient: (A,D,S) = (1,0,0),
430
           which avoided all shadow effects. */
431
        echoMatterPhongSet(glyphsEcho, esquare, 0.4f, 0.6f, 0, 40);
432
        echoListAdd(list, esquare);
433
      }
434
    }
435
    if (parm->onlyPositive) {
436
      if (eval[2] < 0) {
437
        /* didn't have all positive eigenvalues, its outta here */
438
        if (parm->verbose >= 2) {
439
          fprintf(stderr, "%s: glyph %d/%d: not all evals %g %g %g > 0\n",
440
                  me, idx, numGlyphs, eval[0], eval[1], eval[2]);
441
        }
442
        continue;
443
      }
444
    }
445
    if (!( tdata[0] >= parm->confThresh )) {
446
      if (parm->verbose >= 2) {
447
        fprintf(stderr, "%s: glyph %d/%d: conf %g < thresh %g\n",
448
                me, idx, numGlyphs, tdata[0], parm->confThresh);
449
      }
450
      continue;
451
    }
452
    if (!( tenAnisoEval_f(eval, parm->anisoType) >= parm->anisoThresh )) {
453
      if (parm->verbose >= 2) {
454
        fprintf(stderr, "%s: glyph %d/%d: aniso[%d] %g < thresh %g\n",
455
                me, idx, numGlyphs, parm->anisoType,
456
                tenAnisoEval_f(eval, parm->anisoType), parm->anisoThresh);
457
      }
458
      continue;
459
    }
460
    glyphAniso = tenAnisoEval_f(eval, parm->colAnisoType);
461
    /*
462
      fprintf(stderr, "%s: eret = %d; evals = %g %g %g\n", me,
463
      eret, eval[0], eval[1], eval[2]);
464
      ELL_3V_CROSS(tmp1, evec+0, evec+3); tmp2[0] = ELL_3V_LEN(tmp1);
465
      ELL_3V_CROSS(tmp1, evec+0, evec+6); tmp2[1] = ELL_3V_LEN(tmp1);
466
      ELL_3V_CROSS(tmp1, evec+3, evec+6); tmp2[2] = ELL_3V_LEN(tmp1);
467
      fprintf(stderr, "%s: crosses = %g %g %g\n", me,
468
      tmp2[0], tmp2[1], tmp2[2]);
469
    */
470
471
    /* set transform (in mA) */
472
    ELL_3V_ABS(absEval, eval);
473
    ELL_4M_IDENTITY_SET(mA);                        /* reset */
474
    ELL_3V_SCALE(glyphScl, parm->glyphScale, absEval); /* scale by evals */
475
    ELL_4M_SCALE_SET(mB, glyphScl[0], glyphScl[1], glyphScl[2]);
476
477
    ell_4m_post_mul_d(mA, mB);
478
    ELL_43M_INSET(mB, rotEvec);                     /* rotate by evecs */
479
    ell_4m_post_mul_d(mA, mB);
480
    ELL_4M_TRANSLATE_SET(mB, pW[0], pW[1], pW[2]);  /* translate */
481
    ell_4m_post_mul_d(mA, mB);
482
483
    /* set color (in R,G,B) */
484
    cvec = evec + 3*(AIR_CLAMP(0, parm->colEvec, 2));
485
    R = AIR_ABS(cvec[0]);                           /* standard mapping */
486
    G = AIR_ABS(cvec[1]);
487
    B = AIR_ABS(cvec[2]);
488
    /* desaturate by colMaxSat */
489
    R = AIR_AFFINE(0.0, parm->colMaxSat, 1.0, parm->colIsoGray, R);
490
    G = AIR_AFFINE(0.0, parm->colMaxSat, 1.0, parm->colIsoGray, G);
491
    B = AIR_AFFINE(0.0, parm->colMaxSat, 1.0, parm->colIsoGray, B);
492
    /* desaturate some by anisotropy */
493
    R = AIR_AFFINE(0.0, parm->colAnisoModulate, 1.0,
494
                   R, AIR_AFFINE(0.0, glyphAniso, 1.0, parm->colIsoGray, R));
495
    G = AIR_AFFINE(0.0, parm->colAnisoModulate, 1.0,
496
                   G, AIR_AFFINE(0.0, glyphAniso, 1.0, parm->colIsoGray, G));
497
    B = AIR_AFFINE(0.0, parm->colAnisoModulate, 1.0,
498
                   B, AIR_AFFINE(0.0, glyphAniso, 1.0, parm->colIsoGray, B));
499
    /* clamp and do gamma */
500
    R = AIR_CLAMP(0.0, R, 1.0);
501
    G = AIR_CLAMP(0.0, G, 1.0);
502
    B = AIR_CLAMP(0.0, B, 1.0);
503
    R = pow(R, parm->colGamma);
504
    G = pow(G, parm->colGamma);
505
    B = pow(B, parm->colGamma);
506
507
    /* find axis, and superquad exponents qA and qB */
508
    if (eval[2] > 0) {
509
      /* all evals positive */
510
      cl = AIR_MIN(0.99, tenAnisoEval_f(eval, tenAniso_Cl1));
511
      cp = AIR_MIN(0.99, tenAnisoEval_f(eval, tenAniso_Cp1));
512
      if (cl > cp) {
513
        axis = 0;
514
        qA = pow(1-cp, parm->sqdSharp);
515
        qB = pow(1-cl, parm->sqdSharp);
516
      } else {
517
        axis = 2;
518
        qA = pow(1-cl, parm->sqdSharp);
519
        qB = pow(1-cp, parm->sqdSharp);
520
      }
521
      qC = qB;
522
    } else if (eval[0] < 0) {
523
      /* all evals negative */
524
      float aef[3];
525
      aef[0] = absEval[2];
526
      aef[1] = absEval[1];
527
      aef[2] = absEval[0];
528
      cl = AIR_MIN(0.99, tenAnisoEval_f(aef, tenAniso_Cl1));
529
      cp = AIR_MIN(0.99, tenAnisoEval_f(aef, tenAniso_Cp1));
530
      if (cl > cp) {
531
        axis = 2;
532
        qA = pow(1-cp, parm->sqdSharp);
533
        qB = pow(1-cl, parm->sqdSharp);
534
      } else {
535
        axis = 0;
536
        qA = pow(1-cl, parm->sqdSharp);
537
        qB = pow(1-cp, parm->sqdSharp);
538
      }
539
      qC = qB;
540
    } else {
541
#define OOSQRT2 0.70710678118654752440
542
#define OOSQRT3 0.57735026918962576451
543
      /* double poleA[3]={OOSQRT3, OOSQRT3, OOSQRT3}; */
544
      double poleB[3]={1, 0, 0};
545
      double poleC[3]={OOSQRT2, OOSQRT2, 0};
546
      double poleD[3]={OOSQRT3, -OOSQRT3, -OOSQRT3};
547
      double poleE[3]={OOSQRT2, 0, -OOSQRT2};
548
      double poleF[3]={OOSQRT3, OOSQRT3, -OOSQRT3};
549
      double poleG[3]={0, -OOSQRT2, -OOSQRT2};
550
      double poleH[3]={0, 0, -1};
551
      /* double poleI[3]={-OOSQRT3, -OOSQRT3, -OOSQRT3}; */
552
      double funk[3]={0,4,2}, thrn[3]={1,4,4};
553
      double octa[3]={0,2,2}, cone[3]={1,2,2};
554
      double evalN[3], tmp, bary[3];
555
      double qq[3];
556
557
      ELL_3V_NORM(evalN, eval, tmp);
558
      if (eval[1] >= -eval[2]) {
559
        /* inside B-F-C */
560
        ell_3v_barycentric_spherical_d(bary, poleB, poleF, poleC, evalN);
561
        ELL_3V_SCALE_ADD3(qq, bary[0], octa, bary[1], thrn, bary[2], cone);
562
        axis = 2;
563
      } else if (eval[0] >= -eval[2]) {
564
        /* inside B-D-F */
565
        if (eval[1] >= 0) {
566
          /* inside B-E-F */
567
          ell_3v_barycentric_spherical_d(bary, poleB, poleE, poleF, evalN);
568
          ELL_3V_SCALE_ADD3(qq, bary[0], octa, bary[1], funk, bary[2], thrn);
569
          axis = 2;
570
        } else {
571
          /* inside B-D-E */
572
          ell_3v_barycentric_spherical_d(bary, poleB, poleD, poleE, evalN);
573
          ELL_3V_SCALE_ADD3(qq, bary[0], cone, bary[1], thrn, bary[2], funk);
574
          axis = 0;
575
        }
576
      } else if (eval[0] < -eval[1]) {
577
        /* inside D-G-H */
578
        ell_3v_barycentric_spherical_d(bary, poleD, poleG, poleH, evalN);
579
        ELL_3V_SCALE_ADD3(qq, bary[0], thrn, bary[1], cone, bary[2], octa);
580
        axis = 0;
581
      } else if (eval[1] < 0) {
582
        /* inside E-D-H */
583
        ell_3v_barycentric_spherical_d(bary, poleE, poleD, poleH, evalN);
584
        ELL_3V_SCALE_ADD3(qq, bary[0], funk, bary[1], thrn, bary[2], octa);
585
        axis = 0;
586
      } else {
587
        /* inside F-E-H */
588
        ell_3v_barycentric_spherical_d(bary, poleF, poleE, poleH, evalN);
589
        ELL_3V_SCALE_ADD3(qq, bary[0], thrn, bary[1], funk, bary[2], cone);
590
        axis = 2;
591
      }
592
      qA = qq[0];
593
      qB = qq[1];
594
      qC = qq[2];
595
#undef OOSQRT2
596
#undef OOSQRT3
597
    }
598
599
    /* add the glyph */
600
    if (parm->verbose >= 2) {
601
      fprintf(stderr, "%s: glyph %d/%d: the glyph stays!\n",
602
              me, idx, numGlyphs);
603
    }
604
    if (glyphsLimn) {
605
      lookIdx = limnObjectLookAdd(glyphsLimn);
606
      look = glyphsLimn->look + lookIdx;
607
      ELL_4V_SET_TT(look->rgba, float, R, G, B, 1);
608
      ELL_3V_SET(look->kads, parm->ADSP[0], parm->ADSP[1], parm->ADSP[2]);
609
      look->spow = 0;
610
      switch(parm->glyphType) {
611
      case tenGlyphTypeBox:
612
        glyphIdx = limnObjectCubeAdd(glyphsLimn, lookIdx);
613
        break;
614
      case tenGlyphTypeSphere:
615
        glyphIdx = limnObjectPolarSphereAdd(glyphsLimn, lookIdx, axis,
616
                                            2*parm->facetRes, parm->facetRes);
617
        break;
618
      case tenGlyphTypeCylinder:
619
        glyphIdx = limnObjectCylinderAdd(glyphsLimn, lookIdx, axis,
620
                                         parm->facetRes);
621
        break;
622
      case tenGlyphTypeSuperquad:
623
      default:
624
        glyphIdx =
625
          limnObjectPolarSuperquadFancyAdd(glyphsLimn, lookIdx, axis,
626
                                           AIR_CAST(float, qA),
627
                                           AIR_CAST(float, qB),
628
                                           AIR_CAST(float, qC), 0,
629
                                           2*parm->facetRes,
630
                                           parm->facetRes);
631
        break;
632
      }
633
      ELL_4M_COPY_TT(mA_f, float, mA);
634
      limnObjectPartTransform(glyphsLimn, glyphIdx, mA_f);
635
    }
636
    if (glyphsEcho) {
637
      switch(parm->glyphType) {
638
      case tenGlyphTypeBox:
639
        eglyph = echoObjectNew(glyphsEcho, echoTypeCube);
640
        /* nothing else to set */
641
        break;
642
      case tenGlyphTypeSphere:
643
        eglyph = echoObjectNew(glyphsEcho, echoTypeSphere);
644
        echoSphereSet(eglyph, 0, 0, 0, 1);
645
        break;
646
      case tenGlyphTypeCylinder:
647
        eglyph = echoObjectNew(glyphsEcho, echoTypeCylinder);
648
        echoCylinderSet(eglyph, axis);
649
        break;
650
      case tenGlyphTypeSuperquad:
651
      default:
652
        eglyph = echoObjectNew(glyphsEcho, echoTypeSuperquad);
653
        echoSuperquadSet(eglyph, axis, qA, qB);
654
        break;
655
      }
656
      echoColorSet(eglyph,
657
                   AIR_CAST(echoCol_t, R),
658
                   AIR_CAST(echoCol_t, G),
659
                   AIR_CAST(echoCol_t, B), 1);
660
      echoMatterPhongSet(glyphsEcho, eglyph,
661
                         parm->ADSP[0], parm->ADSP[1],
662
                         parm->ADSP[2], parm->ADSP[3]);
663
      inst = echoObjectNew(glyphsEcho, echoTypeInstance);
664
      ELL_4M_COPY(eM, mA);
665
      echoInstanceSet(inst, eM, eglyph);
666
      echoListAdd(list, inst);
667
    }
668
  }
669
  if (glyphsLimn) {
670
    glyphsLimn->setVertexRGBAFromLook = svRGBAfl;
671
  }
672
  if (glyphsEcho) {
673
    split = echoListSplit3(glyphsEcho, list, 10);
674
    echoObjectAdd(glyphsEcho, split);
675
  }
676
677
  airMopOkay(mop);
678
  return 0;
679
}
680
681
/*
682
** Zone from Eval
683
*/
684
unsigned int
685
tenGlyphBqdZoneEval(const double eval[3]) {
686
  double x, y, z;
687
  unsigned int zone;
688
689
  x = eval[0];
690
  y = eval[1];
691
  z = eval[2];
692
  if (y > 0) {   /* 0 1 2 3 4 */
693
    if (z > 0) { /* 0 1 */
694
      if (x - y > y - z) {
695
        zone = 0;
696
      } else {
697
        zone = 1;
698
      }
699
    } else {     /* 2 3 4 */
700
      if (y > -z) {
701
        zone = 2;
702
      } else if (x > -z) {
703
        zone = 3;
704
      } else {
705
        zone = 4;
706
      }
707
    }
708
  } else {       /* 5 6 7 8 9 */
709
    if (x > 0) { /* 5 6 7 */
710
      if (x > -z) {
711
        zone = 5;
712
      } else if (x > -y) {
713
        zone = 6;
714
      } else {
715
        zone = 7;
716
      }
717
    } else {     /* 8 9 */
718
      if (x - y > y - z) {
719
        zone = 8;
720
      } else {
721
        zone = 9;
722
      }
723
    }
724
  }
725
  return zone;
726
}
727
728
/*
729
** UV from Eval
730
*/
731
void
732
tenGlyphBqdUvEval(double uv[2], const double eval[3]) {
733
  double xx, yy, zz, ax, ay, az, mm;
734
735
  ax = AIR_ABS(eval[0]);
736
  ay = AIR_ABS(eval[1]);
737
  az = AIR_ABS(eval[2]);
738
  mm = AIR_MAX(ax, AIR_MAX(ay, az));
739
  if (mm==0) { /* do not divide */
740
    uv[0]=uv[1]=0;
741
    return;
742
  }
743
  xx = eval[0]/mm;
744
  yy = eval[1]/mm;
745
  zz = eval[2]/mm;
746
  uv[0] = AIR_AFFINE(-1, yy, 1, 0, 1);
747
  if (xx > -zz) {
748
    uv[1] = AIR_AFFINE(-1, zz, 1, 0, 1) - uv[0] + 1;
749
  } else {
750
    uv[1] = AIR_AFFINE(-1, xx, 1, -1, 0) - uv[0] + 1;
751
  }
752
  return;
753
}
754
755
/*
756
** Eval from UV
757
*/
758
void
759
tenGlyphBqdEvalUv(double eval[3], const double uv[2]) {
760
  double xx, yy, zz, ll;
761
762
  yy = AIR_AFFINE(0, uv[0], 1, -1, 1);
763
  if (uv[0] + uv[1] > 1) {
764
    zz = AIR_AFFINE(0, uv[1], 1, -1, 1) - 1 + yy;
765
    xx = 1;
766
  } else {
767
    xx = AIR_AFFINE(0, uv[1], 1, -1, 1) + yy + 1;
768
    zz = -1;
769
  }
770
  ELL_3V_SET(eval, xx, yy, zz);
771
  ELL_3V_NORM(eval, eval, ll);
772
  return;
773
}
774
775
/*
776
** Zone from UV
777
*/
778
unsigned int
779
tenGlyphBqdZoneUv(const double uv[2]) {
780
  /* the use of "volatile" here, as well as additional variables for
781
     expressions involving u and v, is based on browsing this summary of the
782
     subtleties of IEEE 754: http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323
783
     In this function, "if (u + v > 0.5)" returned one thing for cygwin, and
784
     something else for other platforms.  Adding volatile and more variables
785
     for expressions brings cygwin back into line with the other platforms */
786
40000
  volatile double u, v, upv, tupv;
787
  unsigned int zone;
788
789
20000
  u = uv[0];
790
20000
  v = uv[1];
791
20000
  upv = u + v;
792
20000
  tupv = 2*u + v;
793

40000
  if (u > 0.5) {       /* 0 1 2 3 4 */
794
30000
    if (upv > 1.5) { /* 0 1 */
795
2450
      if (u < v) {
796
        zone = 0;
797
1200
      } else {
798
        zone = 1;
799
      }
800
    } else {           /* 2 3 4 */
801
7550
      if (tupv > 2) {
802
        zone = 2;
803
7550
      } else if (upv > 1) {
804
        zone = 3;
805
2450
      } else {
806
        zone = 4;
807
      }
808
    }
809
  } else {             /* 5 6 7 8 9 */
810
10000
    if (upv > 0.5) { /* 5 6 7 */
811
7450
      if (upv > 1) {
812
        zone = 5;
813
7450
      } else if (tupv > 1) {
814
        zone = 6;
815
2550
      } else {
816
        zone = 7;
817
      }
818
    } else {           /* 8 9 */
819
2550
      if (u < v) {
820
        zone = 8;
821
1250
      } else {
822
        zone = 9;
823
      }
824
    }
825
  }
826
20000
  return zone;
827
20000
}
828
829
static void
830
baryFind(double bcoord[3], const double uvp[2],
831
         const double uv0[2],
832
         const double uv1[2],
833
         const double uv2[2]) {
834
  double mat[9], a, a01, a02, a12;
835
836
40000
  ELL_3M_SET(mat,
837
             uv0[0], uv0[1], 1,
838
             uv1[0], uv1[1], 1,
839
             uvp[0], uvp[1], 1);
840
20000
  a01 = ELL_3M_DET(mat); a01 = AIR_ABS(a01);
841
842
20000
  ELL_3M_SET(mat,
843
             uv0[0], uv0[1], 1,
844
             uv2[0], uv2[1], 1,
845
             uvp[0], uvp[1], 1);
846
20000
  a02 = ELL_3M_DET(mat); a02 = AIR_ABS(a02);
847
848
  ELL_3M_SET(mat,
849
             uv1[0], uv1[1], 1,
850
             uv2[0], uv2[1], 1,
851
             uvp[0], uvp[1], 1);
852
20000
  a12 = ELL_3M_DET(mat); a12 = AIR_ABS(a12);
853
854
20000
  a = a01 + a02 + a12;
855
20000
  ELL_3V_SET(bcoord, a12/a, a02/a, a01/a);
856
  return;
857
20000
}
858
859
static void
860
baryBlend(double abc[3], const double co[3],
861
          const double abc0[3],
862
          const double abc1[3],
863
          const double abc2[3]) {
864
  unsigned int ii;
865
866
180000
  for (ii=0; ii<3; ii++) {
867
60000
    abc[ii] = co[0]*abc0[ii] + co[1]*abc1[ii] + co[2]*abc2[ii];
868
  }
869
  return;
870
20000
}
871
872
void
873
tenGlyphBqdAbcUv(double abc[3], const double uv[2], double betaMax) {
874
  static const unsigned int vertsZone[10][3] = {{0, 1, 2},   /* 0 */
875
                                                {0, 2, 3},   /* 1 */
876
                                                {1, 3, 4},   /* 2 */
877
                                                {1, 4, 5},   /* 3 */
878
                                                {4, 5, 9},   /* 4 */
879
                                                {1, 5, 6},   /* 5 */
880
                                                {5, 6, 9},   /* 6 */
881
                                                {6, 7, 9},   /* 7 */
882
                                                {7, 8, 10},  /* 8 */
883
                                                {8, 9, 10}}; /* 9 */
884
  static const double uvVert[11][2] = {{1.00, 1.00},   /* 0 */
885
                                       {0.50, 1.00},   /* 1 */
886
                                       {0.75, 0.75},   /* 2 */
887
                                       {1.00, 0.50},   /* 3 */
888
                                       {1.00, 0.00},   /* 4 */
889
                                       {0.50, 0.50},   /* 5 */
890
                                       {0.00, 1.00},   /* 6 */
891
                                       {0.00, 0.50},   /* 7 */
892
                                       {0.25, 0.25},   /* 8 */
893
                                       {0.50, 0.00},   /* 9 */
894
                                       {0.00, 0.00}};  /* 10 */
895
40000
  double abcBall[3], abcCyli[3], abcFunk[3], abcThrn[3],
896
    abcOcta[3], abcCone[3], abcHalf[3];
897
  /* old compile-time setting
898
  const double *abcAll[10][11] = {
899
     zone \ vert 0      1        2        3        4        5        6        7        8        9       10
900
      0    {abcBall, abcCyli, abcHalf,  NULL,    NULL,    NULL,    NULL,    NULL,    NULL,    NULL,    NULL   },
901
      1    {abcBall,  NULL,   abcHalf, abcCyli,  NULL,    NULL,    NULL,    NULL,    NULL,    NULL,    NULL   },
902
      2    { NULL,   abcOcta,  NULL,   abcCone, abcThrn,  NULL,    NULL,    NULL,    NULL,    NULL,    NULL   },
903
      3    { NULL,   abcOcta,  NULL,    NULL,   abcThrn, abcFunk,  NULL,    NULL,    NULL,    NULL,    NULL   },
904
      4    { NULL,    NULL,    NULL,    NULL,   abcThrn, abcFunk,  NULL,    NULL,    NULL,   abcCone,  NULL   },
905
      5    { NULL,   abcCone,  NULL,    NULL,    NULL,   abcFunk, abcThrn,  NULL,    NULL,    NULL,    NULL   },
906
      6    { NULL,    NULL,    NULL,    NULL,    NULL,   abcFunk, abcThrn,  NULL,    NULL,   abcOcta,  NULL   },
907
      7    { NULL,    NULL,    NULL,    NULL,    NULL,    NULL,   abcThrn, abcCone,  NULL,   abcOcta,  NULL   },
908
      8    { NULL,    NULL,    NULL,    NULL,    NULL,    NULL,    NULL,   abcCyli, abcHalf,  NULL,   abcBall },
909
      9    { NULL,    NULL,    NULL,    NULL,    NULL,    NULL,    NULL,    NULL,   abcHalf, abcCyli, abcBall }};
910
  */
911
20000
  const double *abcAll[10][11];
912
  unsigned int pvi[3], zone, vert;
913
20000
  double bcoord[3];
914
915
20000
  ELL_3V_SET(abcBall, 1, 1, 1);
916
20000
  ELL_3V_SET(abcCyli, 1, 0, 0);
917
20000
  ELL_3V_SET(abcFunk, 0, betaMax, 2);  /* only one with c != b  */
918
20000
  ELL_3V_SET(abcThrn, 1, betaMax, 3);
919
20000
  ELL_3V_SET(abcOcta, 0, 2, 2);
920
20000
  ELL_3V_SET(abcCone, 1, 2, 2);
921
20000
  ELL_3V_SET(abcHalf, 0.5, 0.5, 0.5); /* alpha is half-way between alpha of
922
                                         octa and cone and beta has to be
923
                                         the same as alpha at for the
924
                                         seam to be shape-continuous */
925
  /* run-time setting of abcAll[][]; compile-time setting (comments above)
926
     gives "initializer element is not computable at load time" warnings */
927
440000
  for (zone=0; zone<10; zone++) {
928
4800000
    for (vert=0; vert<11; vert++) {
929
2200000
      abcAll[zone][vert]=NULL;
930
    }
931
  }
932
#define SET(zi, vi0, vi1, vi2, sh0, sh1, sh2) \
933
  abcAll[zi][vi0] = abc##sh0; \
934
  abcAll[zi][vi1] = abc##sh1; \
935
  abcAll[zi][vi2] = abc##sh2
936
937
20000
  SET(0, 0, 1, 2, Ball, Cyli, Half);
938
20000
  SET(1, 0, 2, 3, Ball, Half, Cyli);
939
20000
  SET(2, 1, 3, 4, Octa, Cone, Thrn);
940
20000
  SET(3, 1, 4, 5, Octa, Thrn, Funk);
941
20000
  SET(4, 4, 5, 9, Thrn, Funk, Cone);
942
20000
  SET(5, 1, 5, 6, Cone, Funk, Thrn);
943
20000
  SET(6, 5, 6, 9, Funk, Thrn, Octa);
944
20000
  SET(7, 6, 7, 9, Thrn, Cone, Octa);
945
20000
  SET(8, 7, 8,10, Cyli, Half, Ball);
946
20000
  SET(9, 8, 9,10, Half, Cyli, Ball);
947
948
#undef SET
949
950
20000
  zone = tenGlyphBqdZoneUv(uv);
951
20000
  ELL_3V_COPY(pvi, vertsZone[zone]);
952
20000
  baryFind(bcoord, uv, uvVert[pvi[0]], uvVert[pvi[1]], uvVert[pvi[2]]);
953
20000
  baryBlend(abc, bcoord,
954
20000
            abcAll[zone][pvi[0]],
955
20000
            abcAll[zone][pvi[1]],
956
20000
            abcAll[zone][pvi[2]]);
957
  return;
958
20000
}
959