GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: src/echo/color.c Lines: 0 234 0.0 %
Date: 2017-05-26 Branches: 0 81 0.0 %

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 "echo.h"
25
#include "privateEcho.h"
26
27
char _echoBuff[128] = "";
28
29
char *
30
_echoDot(int depth) {
31
  int i;
32
33
  _echoBuff[0] = '\0';
34
  for (i=1; i<=depth; i++) {
35
    strcat(_echoBuff, ".  ");
36
  }
37
  return _echoBuff;
38
}
39
40
void
41
echoTextureLookup(echoCol_t rgba[4], Nrrd *ntext,
42
                  echoPos_t u, echoPos_t v, echoRTParm *parm) {
43
  int su, sv, ui, vi;
44
  float uf, vf;
45
  unsigned char *tdata00, *tdata10, *tdata01, *tdata11;
46
47
  su = ntext->axis[1].size;
48
  sv = ntext->axis[2].size;
49
  if (parm->textureNN) {
50
    ui = airIndex(0.0, u, 1.0, su);
51
    vi = airIndex(0.0, v, 1.0, sv);
52
    tdata00 = (unsigned char*)(ntext->data) + 4*(ui + su*vi);
53
    ELL_4V_SET_TT(rgba, echoCol_t,
54
                  tdata00[0]/255.0, tdata00[1]/255.0,
55
                  tdata00[2]/255.0, tdata00[3]/255.0);
56
  } else {
57
    u = AIR_AFFINE(0.0, u, 1.0, 0.0, su-1);  u = AIR_CLAMP(0, u, su-1);
58
    v = AIR_AFFINE(0.0, v, 1.0, 0.0, sv-1);  v = AIR_CLAMP(0, v, sv-1);
59
    u -= (u == su-1);  ui = (int)u;  uf = AIR_CAST(float, u - ui);
60
    v -= (v == sv-1);  vi = (int)v;  vf = AIR_CAST(float, v - vi);
61
    tdata00 = (unsigned char*)(ntext->data) + 4*(ui + su*vi);
62
    tdata01 = tdata00 + 4;
63
    tdata10 = tdata00 + 4*su;
64
    tdata11 = tdata10 + 4;
65
    ELL_4V_SET_TT(rgba, echoCol_t,
66
                  ((1-vf)*(1-uf)*tdata00[0] + (1-vf)*uf*tdata01[0] +
67
                   vf*(1-uf)*tdata10[0] + vf*uf*tdata11[0])/255.0,
68
                  ((1-vf)*(1-uf)*tdata00[1] + (1-vf)*uf*tdata01[1] +
69
                   vf*(1-uf)*tdata10[1] + vf*uf*tdata11[1])/255.0,
70
                  ((1-vf)*(1-uf)*tdata00[2] + (1-vf)*uf*tdata01[2] +
71
                   vf*(1-uf)*tdata10[2] + vf*uf*tdata11[2])/255.0,
72
                  ((1-vf)*(1-uf)*tdata00[3] + (1-vf)*uf*tdata01[3] +
73
                   vf*(1-uf)*tdata10[3] + vf*uf*tdata11[3])/255.0);
74
  }
75
}
76
77
void
78
echoIntxMaterialColor(echoCol_t rgba[4], echoIntx *intx, echoRTParm *parm) {
79
80
  if (intx->obj->ntext) {
81
    _echoRayIntxUV[intx->obj->type](intx);
82
    echoTextureLookup(rgba, intx->obj->ntext, intx->u, intx->v, parm);
83
    rgba[0] *= intx->obj->rgba[0];
84
    rgba[1] *= intx->obj->rgba[1];
85
    rgba[2] *= intx->obj->rgba[2];
86
    rgba[3] *= intx->obj->rgba[3];
87
  } else {
88
    ELL_4V_COPY(rgba, intx->obj->rgba);
89
  }
90
}
91
92
/*
93
******** echoIntxLightColor
94
**
95
** determines ambient, diffuse, and (Phong) specular contributions to lighting
96
** a given intersection.  "ambi" and "diff" MUST be given as non-NULL,
97
** "spec" can be given as NULL in order to bypass specular computations
98
*/
99
void
100
echoIntxLightColor(echoCol_t ambi[3], echoCol_t diff[3], echoCol_t spec[3],
101
                   echoCol_t sp,
102
                   echoIntx *intx, echoScene *scene, echoRTParm *parm,
103
                   echoThreadState *tstate) {
104
  unsigned int Lidx;
105
  int blocked;
106
  echoRay shadRay;
107
  echoIntx shadIntx;
108
  echoPos_t Ldist, Ldir[3], Lpos[3], Ldot;
109
  echoCol_t Lcol[3], fracseen;
110
111
  if (parm->shadow) {
112
    /* from, neer, shadow */
113
    shadRay.shadow = AIR_TRUE;
114
    ELL_3V_COPY(shadRay.from, intx->pos);
115
    /* HEY: this has to be fixed: shadow rays were getting aborted
116
       because of this in a scene of instanced superquadrics, and
117
       so epsilon had to be increased.  Something about this epsilon
118
       has to be adjusted according to the instancing matrix */
119
    shadRay.neer = 30*ECHO_EPSILON;
120
  }
121
122
  /* set ambient (easy) */
123
  ELL_3V_COPY(ambi, scene->ambi);
124
125
  /* environment contributes only to diffuse */
126
  if (scene->envmap) {
127
    echoEnvmapLookup(diff, intx->norm, scene->envmap);
128
  } else {
129
    ELL_3V_SET(diff, 0, 0, 0);
130
  }
131
132
  /* lights contribute to diffuse and specular */
133
  if (spec) {
134
    ELL_3V_SET(spec, 0, 0, 0);
135
  }
136
  for (Lidx=0; Lidx<scene->lightArr->len; Lidx++) {
137
    echoLightPosition(Lpos, scene->light[Lidx], tstate);
138
    ELL_3V_SUB(Ldir, Lpos, intx->pos);
139
    ELL_3V_NORM(Ldir, Ldir, Ldist);
140
    Ldot = ELL_3V_DOT(Ldir, intx->norm);
141
    /* HEY: HACK: we have to general per-object-type flag that says,
142
       this kind of object has no notion of in-versus-out facing . . . */
143
    if (echoTypeRectangle == intx->obj->type) {
144
      Ldot = AIR_ABS(Ldot);
145
    }
146
    if (Ldot <= 0) {
147
      continue;
148
      /* to next light, we aren't facing this one.  NB: this means
149
         that there aren't diffuse or specular contributions on the
150
         backsides of even semi-transparent surfaces */
151
    }
152
    if (parm->shadow) {
153
      ELL_3V_COPY(shadRay.dir, Ldir);
154
      shadRay.faar = Ldist;
155
      if (echoRayIntx(&shadIntx, &shadRay, scene, parm, tstate)) {
156
        if (1.0 == parm->shadow) {
157
          /* skip to next light, this one is obscured by something,
158
             and we don't do any partial shadowing */
159
          continue;
160
        }
161
        blocked = AIR_TRUE;
162
      } else {
163
        blocked = AIR_FALSE;
164
      }
165
    } else {
166
      blocked = AIR_FALSE;
167
    }
168
    fracseen = AIR_CAST(echoCol_t, blocked ? 1.0 - parm->shadow : 1.0);
169
    echoLightColor(Lcol, Ldist, scene->light[Lidx], parm, tstate);
170
    ELL_3V_SCALE_INCR_TT(diff, echoCol_t, fracseen*Ldot, Lcol);
171
    if (spec) {
172
      Ldot = ELL_3V_DOT(Ldir, intx->refl);
173
      if (echoTypeRectangle == intx->obj->type) {
174
        Ldot = AIR_ABS(Ldot);
175
      }
176
      if (Ldot > 0) {
177
        Ldot = pow(Ldot, sp);
178
        ELL_3V_SCALE_INCR_TT(spec, echoCol_t, fracseen*Ldot, Lcol);
179
      }
180
    }
181
  }
182
  return;
183
}
184
185
void
186
_echoIntxColorPhong(INTXCOLOR_ARGS) {
187
  echoCol_t ambi[3], diff[3], spec[3], ka, kd, ks, sp;
188
189
  ka = intx->obj->mat[echoMatterPhongKa];
190
  kd = intx->obj->mat[echoMatterPhongKd];
191
  ks = intx->obj->mat[echoMatterPhongKs];
192
  sp = intx->obj->mat[echoMatterPhongSp];
193
194
  echoIntxMaterialColor(rgba, intx, parm);
195
  ELL_3V_SET(spec, 0, 0, 0);
196
  echoIntxLightColor(ambi, diff, ks ? spec : NULL, sp,
197
                     intx, scene, parm, tstate);
198
  rgba[0] = rgba[0]*(ka*ambi[0] + kd*diff[0]) + ks*spec[0];
199
  rgba[1] = rgba[1]*(ka*ambi[1] + kd*diff[1]) + ks*spec[1];
200
  rgba[2] = rgba[2]*(ka*ambi[2] + kd*diff[2]) + ks*spec[2];
201
  return;
202
}
203
204
void
205
_echoIntxColorMetal(INTXCOLOR_ARGS) {
206
  echoCol_t ka, kd, kp, RA, RD, RS, ambi[3], diff[3], spec[4];
207
  echoPos_t c;
208
  echoRay reflRay;
209
210
  if (0 && tstate->verbose) {
211
    fprintf(stderr, "%s%s: t = %g\n",
212
            _echoDot(tstate->depth), "_echoIntxColorMetal", intx->t);
213
  }
214
215
  ELL_3V_SET(spec, 0, 0, 0);
216
  echoIntxMaterialColor(rgba, intx, parm);
217
  c = ELL_3V_DOT(intx->view, intx->norm);
218
  if (c <= 0) {
219
    /* see only surface color on backside of metal */
220
    return;
221
  }
222
  c = 1 - c;
223
  c = c*c*c*c*c;
224
  RS = intx->obj->mat[echoMatterMetalR0];
225
  RS = AIR_CAST(echoCol_t, RS + (1 - RS)*c);
226
  ka = intx->obj->mat[echoMatterMetalKa];
227
  kd = intx->obj->mat[echoMatterMetalKd];
228
  kp = ka + kd;
229
  /* neer, faar, shadow, from, dir */
230
  ELL_3V_COPY(reflRay.from, intx->pos);
231
  ELL_3V_COPY(reflRay.dir, intx->refl);
232
  reflRay.neer = ECHO_EPSILON;
233
  reflRay.faar = ECHO_POS_MAX;
234
  reflRay.shadow = AIR_FALSE;
235
  echoRayColor(spec, &reflRay, scene, parm, tstate);
236
  if (kp) {
237
    RA = (1 - RS)*ka/kp;
238
    RD = (1 - RS)*kd/kp;
239
    echoIntxLightColor(ambi, diff, NULL, 0.0,
240
                       intx, scene, parm, tstate);
241
    /* NB: surface color does attenuate reflected color (unlike phong) */
242
    rgba[0] *= RA*ambi[0] + RD*diff[0] + RS*spec[0];
243
    rgba[1] *= RA*ambi[1] + RD*diff[1] + RS*spec[1];
244
    rgba[2] *= RA*ambi[2] + RD*diff[2] + RS*spec[2];
245
  } else {
246
    rgba[0] *= RS*spec[0];
247
    rgba[1] *= RS*spec[1];
248
    rgba[2] *= RS*spec[2];
249
  }
250
251
  return;
252
}
253
254
/*
255
** "th" = theta = angle of incidence
256
** "ph" = phi = angle of refraction
257
** "index" = (index of outgoing material)/(index of incoming material)
258
*/
259
int
260
_echoRefract(echoPos_t T[3], echoPos_t V[3],
261
             echoPos_t N[3], echoCol_t indexr, echoThreadState *tstate) {
262
  echoPos_t cosTh, cosPh, sinPhSq, cosPhSq, tmp1, tmp2;
263
264
  cosTh = ELL_3V_DOT(V, N);
265
  sinPhSq = (1 - cosTh*cosTh)/(indexr*indexr);
266
  cosPhSq = 1 - sinPhSq;
267
  if (cosPhSq < 0) {
268
    if (tstate->verbose) {
269
      fprintf(stderr, "%s%s: cosTh = %g --%g--> TIR!!\n",
270
              _echoDot(tstate->depth), "_echoRefract",
271
              cosTh, indexr);
272
    }
273
    return AIR_FALSE;
274
  }
275
  /* else we do not have total internal reflection */
276
  cosPh = sqrt(cosPhSq);
277
  if (tstate->verbose) {
278
    fprintf(stderr, "%s%s: cosTh = %g --%g--> cosPh = %g\n",
279
            _echoDot(tstate->depth), "_echoRefract",
280
            cosTh, indexr, cosPh);
281
  }
282
  tmp1 = -1.0/indexr; tmp2 = cosTh/indexr - cosPh;
283
  ELL_3V_SCALE_ADD2(T, tmp1, V, tmp2, N);
284
  ELL_3V_NORM(T, T, tmp1);
285
  return AIR_TRUE;
286
}
287
288
void
289
_echoIntxColorGlass(INTXCOLOR_ARGS) {
290
  char me[]="_echoIntxColorGlass";
291
  echoCol_t
292
    ambi[3], diff[3],
293
    ka, kd, RP, RS, RT, R0,
294
    indexr,        /* (index of material we're going into) /
295
                      (index of material we're leaving) */
296
    k[3],          /* attenuation of color due to travel through medium */
297
    matlCol[4],    /* inherent color */
298
    reflCol[4],    /* color from reflected ray */
299
    tranCol[4];    /* color from transmitted ray */
300
  echoPos_t tmp,
301
    negnorm[3];
302
  double c;
303
  echoRay tranRay, reflRay;
304
305
  echoIntxMaterialColor(matlCol, intx, parm);
306
307
  /* from, neer, faar, shadow */
308
  ELL_3V_COPY(tranRay.from, intx->pos);
309
  ELL_3V_COPY(reflRay.from, intx->pos);
310
  tranRay.neer = reflRay.neer = ECHO_EPSILON;
311
  tranRay.faar = reflRay.faar = ECHO_POS_MAX;
312
  tranRay.shadow = reflRay.shadow = AIR_FALSE;
313
  ELL_3V_COPY(reflRay.dir, intx->refl);
314
  /* tranRay.dir set below */
315
  indexr = intx->obj->mat[echoMatterGlassIndex];
316
317
  RS = 0.0;  /* this is a flag meaning: "AFAIK, there's no total int refl" */
318
  tmp = ELL_3V_DOT(intx->norm, intx->view);
319
  if (tmp > 0) {
320
    /* "d.n < 0": we're coming from outside the glass, and we
321
       assume this means that we're going into a HIGHER index material,
322
       which means there is NO total internal reflection */
323
    _echoRefract(tranRay.dir, intx->view, intx->norm, indexr, tstate);
324
    if (tstate->verbose) {
325
      fprintf(stderr, "%s%s: V=(%g,%g,%g),N=(%g,%g,%g),n=%g -> T=(%g,%g,%g)\n",
326
              _echoDot(tstate->depth), me,
327
              intx->view[0], intx->view[1], intx->view[2],
328
              intx->norm[0], intx->norm[1], intx->norm[2], indexr,
329
              tranRay.dir[0], tranRay.dir[1], tranRay.dir[2]);
330
    }
331
    c = tmp;
332
    ELL_3V_SET(k, 1, 1, 1);
333
  } else {
334
    /* we're coming from inside the glass */
335
    /* the reasoning for my Beer's law implementation is this: if a
336
       channel (r, g, or b) is full on (1.0), then there should be no
337
       attenuation in its color.  The more the color is below 1.0, the
338
       more it should be damped with distance. */
339
    k[0] = AIR_CAST(echoCol_t, exp(parm->glassC*(matlCol[0]-1)*intx->t));
340
    k[1] = AIR_CAST(echoCol_t, exp(parm->glassC*(matlCol[1]-1)*intx->t));
341
    k[2] = AIR_CAST(echoCol_t, exp(parm->glassC*(matlCol[2]-1)*intx->t));
342
    if (tstate->verbose) {
343
      fprintf(stderr, "%s%s: internal refl @ t = %g -> k = %g %g %g\n",
344
              _echoDot(tstate->depth), me, intx->t, k[0], k[1], k[2]);
345
    }
346
    ELL_3V_SCALE(negnorm, -1, intx->norm);
347
    if (_echoRefract(tranRay.dir, intx->view, negnorm, 1/indexr, tstate)) {
348
      c = -ELL_3V_DOT(tranRay.dir, negnorm);
349
    } else {
350
      /* its total internal reflection time! */
351
      c = 0.0;
352
      RS = 1.0;
353
    }
354
  }
355
356
  if (RS) {
357
    /* total internal reflection */
358
    RT = 0;
359
  } else {
360
    R0 = (indexr - 1)/(indexr + 1);
361
    R0 *= R0;
362
    c = 1 - c;
363
    c = c*c*c*c*c;
364
    RS = AIR_CAST(echoCol_t, R0 + (1-R0)*c);
365
    RT = 1 - RS;
366
  }
367
  ka = intx->obj->mat[echoMatterMetalKa];
368
  kd = intx->obj->mat[echoMatterMetalKd];
369
  RP = ka + kd;
370
  if (RP) {
371
    RS *= 1 - RP;
372
    RT *= 1 - RP;
373
    echoIntxLightColor(ambi, diff, NULL, 0.0,
374
                       intx, scene, parm, tstate);
375
  } else {
376
    ELL_3V_SET(ambi, 0, 0, 0);
377
    ELL_3V_SET(diff, 0, 0, 0);
378
  }
379
  if (tstate->verbose) {
380
    fprintf(stderr, "%s%s: --- reflRay (reflected)\n",
381
            _echoDot(tstate->depth), me);
382
  }
383
  echoRayColor(reflCol, &reflRay, scene, parm, tstate);
384
  if (RT) {
385
    if (tstate->verbose) {
386
      fprintf(stderr, "%s%s: --- tranRay (refracted)\n",
387
              _echoDot(tstate->depth), me);
388
    }
389
    echoRayColor(tranCol, &tranRay, scene, parm, tstate);
390
  } else {
391
    ELL_3V_SET(tranCol, 0, 0, 0);
392
  }
393
  rgba[0] = (matlCol[0]*(ka*ambi[0] + kd*diff[0])  +
394
             k[0]*(RS*reflCol[0] + RT*tranCol[0]));
395
  rgba[1] = (matlCol[1]*(ka*ambi[1] + kd*diff[1])  +
396
             k[1]*(RS*reflCol[1] + RT*tranCol[1]));
397
  rgba[2] = (matlCol[2]*(ka*ambi[2] + kd*diff[2])  +
398
             k[2]*(RS*reflCol[2] + RT*tranCol[2]));
399
  rgba[3] = 1.0;
400
  return;
401
}
402
403
void
404
_echoIntxColorLight(INTXCOLOR_ARGS) {
405
406
  AIR_UNUSED(scene);
407
  AIR_UNUSED(tstate);
408
  echoIntxMaterialColor(rgba, intx, parm);
409
}
410
411
void
412
_echoIntxColorUnknown(INTXCOLOR_ARGS) {
413
414
  AIR_UNUSED(rgba);
415
  AIR_UNUSED(intx);
416
  AIR_UNUSED(scene);
417
  AIR_UNUSED(parm);
418
  fprintf(stderr, "%s%s: can't color intx with object with unset material\n",
419
          _echoDot(tstate->depth), "_echoIntxColorNone");
420
}
421
422
_echoIntxColor_t
423
_echoIntxColor[ECHO_MATTER_MAX+1] = {
424
  _echoIntxColorUnknown,
425
  _echoIntxColorPhong,
426
  _echoIntxColorGlass,
427
  _echoIntxColorMetal,
428
  _echoIntxColorLight,
429
};
430
431
/*
432
******** echoIntxFuzzify()
433
**
434
** this modifies the refl and norm fields of a given intx
435
*/
436
void
437
echoIntxFuzzify(echoIntx *intx, echoCol_t fuzz, echoThreadState *tstate) {
438
  echoPos_t tmp, *jitt, oldNorm[3], perp0[3], perp1[3], jj0, jj1;
439
  int side;
440
441
  /* at some point I thought this was important to avoid bias when
442
     going through glass, but now I'm not so sure . . . It is likely
443
     totally moot if jitter vectors are NOT reused between pixels. */
444
  if (ELL_3V_DOT(intx->norm, intx->view) > 0) {
445
    jitt = tstate->jitt + 2*echoJittableNormalA;
446
  } else {
447
    jitt = tstate->jitt + 2*echoJittableNormalB;
448
  }
449
  jj0 = fuzz*jitt[0];
450
  jj1 = fuzz*jitt[1];
451
  ELL_3V_COPY(oldNorm, intx->norm);
452
  side = ELL_3V_DOT(intx->refl, oldNorm) > 0;
453
  ell_3v_PERP(perp0, oldNorm);
454
  ELL_3V_NORM(perp0, perp0, tmp);
455
  ELL_3V_CROSS(perp1, perp0, oldNorm);
456
  ELL_3V_SCALE_ADD3(intx->norm, 1, oldNorm, jj0, perp0, jj1, perp1);
457
  ELL_3V_NORM(intx->norm, intx->norm, tmp);
458
  _ECHO_REFLECT(intx->refl, intx->norm, intx->view, tmp);
459
  if (side != (ELL_3V_DOT(intx->refl, oldNorm) > 0)) {
460
    ELL_3V_SCALE_ADD3(intx->norm, 1, oldNorm, -jj0, perp0, -jj1, perp1);
461
    ELL_3V_NORM(intx->norm, intx->norm, tmp);
462
    _ECHO_REFLECT(intx->refl, intx->norm, intx->view, tmp);
463
  }
464
  if (tstate->verbose) {
465
    fprintf(stderr, "%s%s: fuzz[%g](%g,%g,%g) --> (%g,%g,%g)\n",
466
            _echoDot(tstate->depth), "echoIntxFuzzify", fuzz,
467
            oldNorm[0], oldNorm[1], oldNorm[2],
468
            intx->norm[0], intx->norm[1], intx->norm[2]);
469
  }
470
  return;
471
}
472
473
void
474
echoIntxColor(echoCol_t rgba[4], echoIntx *intx,
475
              echoScene *scene, echoRTParm *parm, echoThreadState *tstate) {
476
  echoCol_t fuzz;
477
478
  switch(intx->obj->matter) {
479
  case echoMatterGlass:
480
    fuzz = intx->obj->mat[echoMatterGlassFuzzy];
481
    break;
482
  case echoMatterMetal:
483
    fuzz = intx->obj->mat[echoMatterMetalFuzzy];
484
    break;
485
  default:
486
    fuzz = 0;
487
    break;
488
  }
489
  if (fuzz) {
490
    echoIntxFuzzify(intx, fuzz, tstate);
491
  }
492
  _echoIntxColor[intx->obj->matter](rgba, intx, scene, parm, tstate);
493
}
494