1 /* Copyright (C) 1992-1998 The Geometry Center
2 * Copyright (C) 1998-2000 Stuart Levy, Tamara Munzner, Mark Phillips
3 *
4 * This file is part of Geomview.
5 *
6 * Geomview is free software; you can redistribute it and/or modify it
7 * under the terms of the GNU Lesser General Public License as published
8 * by the Free Software Foundation; either version 2, or (at your option)
9 * any later version.
10 *
11 * Geomview is distributed in the hope that it will be useful, but
12 * WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with Geomview; see the file COPYING. If not, write
18 * to the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139,
19 * USA, or visit http://www.gnu.org.
20 */
21
22 #if HAVE_CONFIG_H
23 # include "config.h"
24 #endif
25
26 #if 0
27 static char copyright[] = "Copyright (C) 1992-1998 The Geometry Center\n\
28 Copyright (C) 1998-2000 Stuart Levy, Tamara Munzner, Mark Phillips";
29 #endif
30
31 /*
32 * Software shaders.
33 * These operate in the mg environment.
34 */
35 #include "drawer.h"
36 #include "mgP.h"
37 #include "3d.h"
38
39 #define DOT(a,b) ((a).x*(b).x + (a).y*(b).y + (a).z*(b).z)
40
41 #define CCMUL(A,B,C) ((C).r=(A).r*(B).r, (C).g=(A).g*(B).g, (C).b=(A).b*(B).b)
42 #define CCADD(A,B,C) ((C).r=(A).r+(B).r, (C).g=(A).g+(B).g, (C).b=(A).b+(B).b)
43
44 #define CSCALE(s, A, B) ( (B).r = (A).r * s, (B).g = (A).g * s, (B).b = (A).b * s, (B).a = (A).b * s)
45
46
47 #define CCSADD(A,B,s,C) \
48 ((C).r += (s)*(A).r*(B).r, \
49 (C).g += (s)*(A).g*(B).g, \
50 (C).b += (s)*(A).b*(B).b)
51
52 /* Eye position as an HPoint3 */
53 #define EYE ((HPoint3 *) (_mgc->C2W[3]))
54
55
56
57 /* Shading parameters, left here for easy dbx viewing. */
58 static float fog = 0.0; /* fog control */
59 static float hypk1 = .5; /* controls exponential falloff */
60 static float hypk2 = .0; /* controls exponential falloff */
61 static float sphk1 = .5; /* controls exponential falloff */
62 static float sphk2 = 0.0; /* safety zone around light */
63
64 /*
65 * Note: code for the Euclidean shader now lives in src/lib/mg/common/mgshade.c.
66 * This file contains only the non-Euclidean ones.
67 */
68
69 #define HPt3R31GramSchmidt HPt3HypGramSchmidt
70 #define HPt3R31Distance HPt3HypDistance
71
72 #define HPt3R40GramSchmidt HPt3SphGramSchmidt
73 #define HPt3R40Distance HPt3SphDistance
74
75 /*
76 * A routine to determine the bisector of two tangent vectors
77 */
HPt3R31Bisector(HPoint3 * Lt,HPoint3 * I,HPoint3 * H)78 void HPt3R31Bisector(HPoint3 *Lt, HPoint3 *I, HPoint3 *H)
79 {
80 HPoint3 tmp1, tmp2;
81 float a,b;
82
83 HPt3Sub(I, Lt, &tmp1);
84 a = HPt3R31Dot(I, &tmp1);
85 b = HPt3R31Dot(Lt, &tmp1);
86 if ( b == 0) b = 1;
87 HPt3Scale(-(a/b), Lt, &tmp2);
88 HPt3Add(I, &tmp2, H);
89 if (H->w < 0) HPt3Scale(-1, H, H);
90 return;
91 }
92
HPt3R40Bisector(HPoint3 * Lt,HPoint3 * I,HPoint3 * H)93 void HPt3R40Bisector(HPoint3 *Lt, HPoint3 *I, HPoint3 *H)
94 {
95 HPoint3 tmp1, tmp2;
96 float a,b;
97
98 HPt3Sub(I, Lt, &tmp1);
99 a = HPt3R40Dot(I, &tmp1);
100 b = HPt3R40Dot(Lt, &tmp1);
101 if ( b == 0) b = 1;
102 HPt3Scale(-(a/b), Lt, &tmp2);
103 HPt3Add(I, &tmp2, H);
104 return;
105 }
106
107 /*
108 * Hyperbolic shader
109 */
110
hypshade(int nv,HPoint3 * v,Point3 * n,ColorA * c,ColorA * cs)111 int hypshade(int nv, HPoint3 *v, Point3 *n, ColorA *c, ColorA *cs)
112 {
113 struct mgxstk *mx = _mgc->xstk;
114 struct mgastk *ma = _mgc->astk;
115 struct LtLight *l;
116 int i, lno;
117 HPoint3 V, n4, N, I;
118 float s, ni;
119 Color Ci;
120 static LtLight *cachedlights[AP_MAXLIGHTS];
121 static int ncached = 0;
122 static short cached_light_seq = -1;
123
124 hypk1 = ma->lighting.attenmult;
125 hypk2 = ma->lighting.attenmult2;
126 fog = ma->lighting.attenconst>1 ? ma->lighting.attenconst-1 : 0;
127
128 Ci.r = Ci.g = Ci.b = 0;
129 CCSADD(ma->mat.ambient, ma->lighting.ambient, ma->mat.ka, Ci);
130
131 /* for now, move lights to be inside hyperbolic space */
132
133 if(ma->light_seq != cached_light_seq) {
134 LtLight **lp;
135 float d;
136 HPoint3 *L;
137
138 for(i = 0; i < ncached; i++)
139 LtDelete(cachedlights[i]);
140 ncached = 0;
141 LM_FOR_ALL_LIGHTS(&ma->lighting, i,lp) {
142 if(ncached >= AP_MAXLIGHTS)
143 break;
144
145 cachedlights[ncached++] = l = LtCopy(*lp, NULL);
146
147 /* use euclidean distance */
148
149 L = &l->globalposition;
150 d = HPt3R30Dot(L, L);
151 if (d >= 1) { /* move light to sphere of radius .8 */
152 d = .8/sqrt(d);
153 L->x *= d;
154 L->y *= d;
155 L->z *= d;
156 L->w = 1.0;
157 }
158 }
159 cached_light_seq = ma->light_seq;
160 }
161
162 for(i = 0; i < nv; i++, v++, n++, c++, cs++) {
163
164 /* Transform point by T */
165 HPt3Transform(mx->T, v, &V); /* V = v in world coords */
166 if (HPt3R31Dot(&V, &V) >= 0) /* point at infinity, or beyond */
167 continue; /* outside hyperbolic space */
168 HPt3R31Normalize(&V); /* lie on surface < , > = -1 */
169
170 HPt3Sub( (HPoint3 *) EYE, &V, &I ); /* I = EYE - V (point-to-eye vector) */
171 HPt3R31GramSchmidt(&V, &I);
172 HPt3R31Normalize(&I);
173
174 n4.x = n->x; n4.y = n->y; n4.z = n->z; n4.w = 1.0;
175 HPt3R31GramSchmidt( v, &n4);
176 HPt3R31Normalize(&n4);
177
178 HPt3Transform(mx->T, &n4, &N);/* transform by inverse adjoint */
179
180 /* tangent space at V has a Euclidean (positive definite) structure */
181 ni = HPt3R31Dot(&N,&I);
182 if (ni < 0)
183 {
184 float s = -1.0;
185 HPt3Scale(s, &N, &N);
186 }
187
188
189 for(lno = ncached; --lno >= 0; ) {
190 HPoint3 *L;
191 HPoint3 Lt;
192 float bright, ll, ln, light_intensity, d;
193
194 l = cachedlights[lno];
195 L = &l->globalposition;
196 /* Diffuse term */
197 /* first compute the color of arriving light */
198 bright = l->intensity;
199 ll = HPt3R31Dot(L,L);
200 if (ll >= 0) /* ignore lights outside of hspace */
201 continue;
202 ln = HPt3R31Dot(L,&N);
203 if(ln <= 0) /* Ignore lights shining from behind */
204 continue;
205 /* compute the distance for computation of falloff of light */
206 d = HPt3R31Distance(L, &V);
207 /* the following models the exponential falloff of intensity */
208 light_intensity = exp(-hypk1*d);
209
210 /* now compute the cosine for diffuse shading */
211 Lt = l->globalposition; /* we'll change the values here */
212 HPt3Sub(&Lt, &V, &Lt); /* make it a difference vector */
213 HPt3R31GramSchmidt(&V, &Lt);
214 HPt3R31Normalize(&Lt);
215 d = HPt3R31Dot(&Lt, &N); /* cos of angle between L and N */
216
217 s = ma->mat.kd * bright * d * light_intensity;
218
219 CCSADD(l->color, *c, s, Ci);
220
221 if(ma->mat.ks > 0) { /* Specular term */
222 HPoint3 H; /* H: halfway between L and I */
223
224 HPt3R31Bisector(&Lt, &I, &H);
225 HPt3R31Normalize(&H);
226
227 s = HPt3R31Dot(&H, &N); /* cos of angle between H and N */
228 if (s < 0) continue;
229
230 s = ma->mat.ks * bright * pow(s, ma->mat.shininess);
231
232 CCSADD(l->color, ma->mat.specular, s, Ci);
233 }
234 }
235 /* insert fog code */
236 if (fog) {
237 float k1, k2;
238 float d = HPt3HypDistance(&V, EYE);
239 ColorA surfcolor, fogcolor;
240 d = d - hypk2; /* fog-free sphere of radius euck2 */
241 if (d < 0) d = 0;
242 k1 = exp( -fog * d);
243 k2 = 1.0 - k1;
244 CSCALE(k1, Ci, surfcolor);
245 CSCALE(k2, _mgc->background, fogcolor);
246 CCADD(surfcolor, fogcolor, Ci);
247 }
248 if(Ci.r < 0) Ci.r = 0; else if(Ci.r > 1) Ci.r = 1;
249 if(Ci.g < 0) Ci.g = 0; else if(Ci.g > 1) Ci.g = 1;
250 if(Ci.b < 0) Ci.b = 0; else if(Ci.b > 1) Ci.b = 1;
251 *(Color *)cs = Ci;
252 cs->a = c->a;
253 }
254 return 0;
255 }
256
sphshade(int nv,HPoint3 * v,Point3 * n,ColorA * c,ColorA * cs)257 int sphshade(int nv, HPoint3 *v, Point3 *n, ColorA *c, ColorA *cs)
258 {
259 struct mgxstk *mx = _mgc->xstk;
260 struct mgastk *ma = _mgc->astk;
261 struct LtLight *l;
262 int i, lno;
263 HPoint3 V, n4, N, I;
264 float s, ni;
265 Color Ci;
266 static LtLight *cachedlights[AP_MAXLIGHTS];
267 static int ncached = 0;
268 static short cached_light_seq = -1;
269
270 sphk1 = ma->lighting.attenmult;
271 sphk2 = ma->lighting.attenmult2;
272 fog = ma->lighting.attenconst>1 ? ma->lighting.attenconst-1 : 0;
273
274 /* "ambient" color */
275
276 Ci.r = Ci.g = Ci.b = 0;
277 CCSADD(ma->mat.ambient, ma->lighting.ambient, ma->mat.ka, Ci);
278
279 /* for now, move lights to be inside spherical space */
280 if(ma->light_seq != cached_light_seq) {
281 float d;
282 HPoint3 *L;
283 LtLight **lp;
284
285 for(i = 0; i < ncached; i++)
286 LtDelete(cachedlights[i]);
287 ncached = 0;
288 LM_FOR_ALL_LIGHTS(&ma->lighting, i,lp) {
289 if(ncached >= AP_MAXLIGHTS)
290 break;
291
292 cachedlights[ncached++] = l = LtCopy(*lp, NULL);
293
294 /* use euclidean distance */
295
296 L = &l->globalposition;
297 /* use euclidean distance */
298 d = HPt3R40Dot(L, L);
299 if (d >= 1.0) { /* move light to be in S3 */
300 d = 1.0/sqrt(d);
301 L->x *= d;
302 L->y *= d;
303 L->z *= d;
304 L->w *= d;
305 }
306 }
307 cached_light_seq = ma->light_seq;
308 }
309
310 for(i = 0; i < nv; i++, v++, n++, c++, cs++) {
311
312 /* Transform point by T */
313 HPt3Transform(mx->T, v, &V); /* V = v in world coords */
314 HPt3R40Normalize(&V); /* lie on surface < , > = -1 */
315 /* I = EYE - V (point-to-eye vector) */
316 HPt3Sub( (HPoint3 *) EYE, &V, &I );
317 HPt3R40GramSchmidt(&V, &I);
318 HPt3R40Normalize(&I);
319
320 n4.x = n->x; n4.y = n->y; n4.z = n->z; n4.w = 1.0;
321 HPt3R40GramSchmidt( v, &n4);
322 HPt3R40Normalize(&n4);
323
324 HPt3Transform(mx->T, &n4, &N);/* transform by inverse adjoint */
325
326 /* tangent space at V has a Euclidean (positive definite) structure */
327 ni = HPt3R40Dot(&N,&I);
328 if (ni < 0)
329 {
330 float s = -1.0;
331 HPt3Scale(s, &N, &N);
332 }
333
334 for(lno = ncached; --lno >= 0; ) {
335 HPoint3 *L;
336 HPoint3 Lt;
337 float bright, ln, d, radius, light_intensity;
338
339 l = cachedlights[lno];
340 L = &l->globalposition;
341 /* Diffuse term */
342 /* first compute the color of arriving light */
343 bright = l->intensity;
344 /*ll = HPt3R40Dot(L,L);*/ /* this should be 1.0 by above */
345 ln = HPt3R40Dot(L,&N);
346 if(ln <= 0) /* Ignore lights shining from behind */
347 continue;
348 /* compute the distance for computation of falloff of light */
349 d = HPt3R40Distance(L, &V);
350 /* the amount of light from the source which arrives at
351 the surface is inversely proportional to the area
352 of the sphere centered at the light which contains
353 the surface point. The radius of this sphere is the sin
354 of the distance between L and V */
355 radius = sin((double) d);
356 if (radius > sphk2)
357 {
358 float kk = (sphk2/radius); /* 0 <= sphk2 <= 1.0 */
359 light_intensity = pow(kk, -sphk1);
360 }
361 else light_intensity = 1.0;
362 /* we model the atmosphere as a murky fluid that has exponential
363 dissipation; sphk1 controls exponential drop off; value of 0
364 avoids any drop off */
365 d = HPt3R40Distance((HPoint3 *) EYE, &V);
366 light_intensity *= exp(-sphk1 * d);
367
368 /* now compute the cosine for diffuse shading */
369 Lt = l->globalposition; /* we'll change the values here */
370 HPt3Sub(&Lt, &V, &Lt); /* make it a difference vector */
371 HPt3R40GramSchmidt(&V, &Lt);
372 HPt3R40Normalize(&Lt);
373 d = HPt3R40Dot(&Lt, &N); /* cos of angle between L and N */
374
375 s = ma->mat.kd * bright * d * light_intensity;
376
377 CCSADD(l->color, *c, s, Ci);
378
379 if(ma->mat.ks > 0) { /* Specular term */
380 HPoint3 H; /* H: halfway between L and I */
381
382 HPt3R40Bisector(&Lt, &I, &H);
383 HPt3R40Normalize(&H);
384
385 s = HPt3R40Dot(&H, &N); /* cos of angle between H and N */
386 if (s < 0) continue;
387
388 s = ma->mat.ks * bright * pow(s, ma->mat.shininess);
389
390 CCSADD(l->color, ma->mat.specular, s, Ci);
391 }
392 }
393 /* insert fog code */
394 if (fog) {
395 float k1, k2;
396 float d = HPt3SphDistance(&V, EYE);
397 ColorA surfcolor, fogcolor;
398
399 d = d - sphk2; /* fog-free sphere of radius euck2 */
400 if (d < 0) d = 0;
401 k1 = exp( -fog * d);
402 k2 = 1.0 - k1;
403 CSCALE(k1, Ci, surfcolor);
404 CSCALE(k2, _mgc->background, fogcolor);
405 CCADD(surfcolor, fogcolor, Ci);
406 }
407 if(Ci.r < 0) Ci.r = 0; else if(Ci.r > 1) Ci.r = 1;
408 if(Ci.g < 0) Ci.g = 0; else if(Ci.g > 1) Ci.g = 1;
409 if(Ci.b < 0) Ci.b = 0; else if(Ci.b > 1) Ci.b = 1;
410 *(Color *)cs = Ci;
411 cs->a = c->a;
412 }
413 return 0;
414 }
415
416
417 mgshadefunc
softshader(int camid)418 softshader(int camid)
419 {
420 switch(spaceof(camid)) {
421 case TM_HYPERBOLIC: return hypshade;
422 case TM_EUCLIDEAN: return mg_eucshade;
423 case TM_SPHERICAL: return sphshade;
424 default: return mg_eucshade;
425 }
426 }
427