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 "echo.h"
25 #include "privateEcho.h"
26
27 char _echoBuff[128] = "";
28
29 char *
_echoDot(int depth)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
echoTextureLookup(echoCol_t rgba[4],Nrrd * ntext,echoPos_t u,echoPos_t v,echoRTParm * parm)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
echoIntxMaterialColor(echoCol_t rgba[4],echoIntx * intx,echoRTParm * parm)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
echoIntxLightColor(echoCol_t ambi[3],echoCol_t diff[3],echoCol_t spec[3],echoCol_t sp,echoIntx * intx,echoScene * scene,echoRTParm * parm,echoThreadState * tstate)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
_echoIntxColorPhong(INTXCOLOR_ARGS)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
_echoIntxColorMetal(INTXCOLOR_ARGS)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
_echoRefract(echoPos_t T[3],echoPos_t V[3],echoPos_t N[3],echoCol_t indexr,echoThreadState * tstate)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
_echoIntxColorGlass(INTXCOLOR_ARGS)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
_echoIntxColorLight(INTXCOLOR_ARGS)404 _echoIntxColorLight(INTXCOLOR_ARGS) {
405
406 AIR_UNUSED(scene);
407 AIR_UNUSED(tstate);
408 echoIntxMaterialColor(rgba, intx, parm);
409 }
410
411 void
_echoIntxColorUnknown(INTXCOLOR_ARGS)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
echoIntxFuzzify(echoIntx * intx,echoCol_t fuzz,echoThreadState * tstate)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
echoIntxColor(echoCol_t rgba[4],echoIntx * intx,echoScene * scene,echoRTParm * parm,echoThreadState * tstate)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
495