1 /*
2     Copyright (c) 2005-2021 Intel Corporation
3 
4     Licensed under the Apache License, Version 2.0 (the "License");
5     you may not use this file except in compliance with the License.
6     You may obtain a copy of the License at
7 
8         http://www.apache.org/licenses/LICENSE-2.0
9 
10     Unless required by applicable law or agreed to in writing, software
11     distributed under the License is distributed on an "AS IS" BASIS,
12     WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13     See the License for the specific language governing permissions and
14     limitations under the License.
15 */
16 
17 /*
18     The original source for this example is
19     Copyright (c) 1994-2008 John E. Stone
20     All rights reserved.
21 
22     Redistribution and use in source and binary forms, with or without
23     modification, are permitted provided that the following conditions
24     are met:
25     1. Redistributions of source code must retain the above copyright
26        notice, this list of conditions and the following disclaimer.
27     2. Redistributions in binary form must reproduce the above copyright
28        notice, this list of conditions and the following disclaimer in the
29        documentation and/or other materials provided with the distribution.
30     3. The name of the author may not be used to endorse or promote products
31        derived from this software without specific prior written permission.
32 
33     THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS
34     OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
35     WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
36     ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
37     DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
38     DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
39     OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
40     HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
41     LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
42     OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
43     SUCH DAMAGE.
44 */
45 
46 /*
47  * extvol.cpp - Volume rendering helper routines etc.
48  */
49 
50 #include <cstdio>
51 
52 #include "machine.hpp"
53 #include "types.hpp"
54 #include "macros.hpp"
55 #include "vector.hpp"
56 #include "util.hpp"
57 #include "box.hpp"
58 #include "extvol.hpp"
59 #include "trace.hpp"
60 #include "sphere.hpp"
61 #include "light.hpp"
62 #include "shade.hpp"
63 #include "global.hpp"
64 
65 int extvol_bbox(void *obj, vector *min, vector *max) {
66     box *b = (box *)obj;
67 
68     *min = b->min;
69     *max = b->max;
70 
71     return 1;
72 }
73 
74 static object_methods extvol_methods = { (void (*)(void *, void *))(box_intersect),
75                                          (void (*)(void *, void *, void *, void *))(box_normal),
76                                          extvol_bbox,
77                                          free };
78 
79 extvol *newextvol(void *voidtex,
80                   vector min,
81                   vector max,
82                   int samples,
83                   flt (*evaluator)(flt, flt, flt)) {
84     extvol *xvol;
85     texture *tex;
86 
87     tex = (texture *)voidtex;
88 
89     xvol = (extvol *)rt_getmem(sizeof(extvol));
90     memset(xvol, 0, sizeof(extvol));
91 
92     xvol->methods = &extvol_methods;
93 
94     xvol->min = min;
95     xvol->max = max;
96     xvol->evaluator = evaluator;
97     xvol->ambient = tex->ambient;
98     xvol->diffuse = tex->diffuse;
99     xvol->opacity = tex->opacity;
100     xvol->samples = samples;
101 
102     xvol->tex = (texture *)rt_getmem(sizeof(texture));
103     memset(xvol->tex, 0, sizeof(texture));
104 
105     xvol->tex->ctr.x = 0.0;
106     xvol->tex->ctr.y = 0.0;
107     xvol->tex->ctr.z = 0.0;
108     xvol->tex->rot = xvol->tex->ctr;
109     xvol->tex->scale = xvol->tex->ctr;
110     xvol->tex->uaxs = xvol->tex->ctr;
111     xvol->tex->vaxs = xvol->tex->ctr;
112     xvol->tex->islight = 0;
113     xvol->tex->shadowcast = 0;
114 
115     xvol->tex->col = tex->col;
116     xvol->tex->ambient = 1.0;
117     xvol->tex->diffuse = 0.0;
118     xvol->tex->specular = 0.0;
119     xvol->tex->opacity = 1.0;
120     xvol->tex->img = nullptr;
121     xvol->tex->texfunc = (color(*)(void *, void *, void *))(ext_volume_texture);
122     xvol->tex->obj = (void *)xvol; /* XXX hack! */
123 
124     return xvol;
125 }
126 
127 color ExtVoxelColor(flt scalar) {
128     color col;
129 
130     if (scalar > 1.0)
131         scalar = 1.0;
132 
133     if (scalar < 0.0)
134         scalar = 0.0;
135 
136     if (scalar < 0.5) {
137         col.g = 0.0;
138     }
139     else {
140         col.g = (scalar - 0.5) * 2.0;
141     }
142 
143     col.r = scalar;
144     col.b = 1.0 - (scalar / 2.0);
145 
146     return col;
147 }
148 
149 color ext_volume_texture(vector *hit, texture *tex, ray *ry) {
150     color col, col2;
151     box *bx;
152     extvol *xvol;
153     flt a, tx1, tx2, ty1, ty2, tz1, tz2;
154     flt tnear, tfar;
155     flt t, tdist, dt, ddt, sum, tt;
156     vector pnt, bln;
157     flt scalar, transval;
158     int i;
159     point_light *li;
160     color diffint;
161     vector N, L;
162     flt inten;
163 
164     col.r = 0.0;
165     col.g = 0.0;
166     col.b = 0.0;
167 
168     bx = (box *)tex->obj;
169     xvol = (extvol *)tex->obj;
170 
171     tnear = -FHUGE;
172     tfar = FHUGE;
173 
174     if (ry->d.x == 0.0) {
175         if ((ry->o.x < bx->min.x) || (ry->o.x > bx->max.x))
176             return col;
177     }
178     else {
179         tx1 = (bx->min.x - ry->o.x) / ry->d.x;
180         tx2 = (bx->max.x - ry->o.x) / ry->d.x;
181         if (tx1 > tx2) {
182             a = tx1;
183             tx1 = tx2;
184             tx2 = a;
185         }
186         if (tx1 > tnear)
187             tnear = tx1;
188         if (tx2 < tfar)
189             tfar = tx2;
190     }
191     if (tnear > tfar)
192         return col;
193     if (tfar < 0.0)
194         return col;
195 
196     if (ry->d.y == 0.0) {
197         if ((ry->o.y < bx->min.y) || (ry->o.y > bx->max.y))
198             return col;
199     }
200     else {
201         ty1 = (bx->min.y - ry->o.y) / ry->d.y;
202         ty2 = (bx->max.y - ry->o.y) / ry->d.y;
203         if (ty1 > ty2) {
204             a = ty1;
205             ty1 = ty2;
206             ty2 = a;
207         }
208         if (ty1 > tnear)
209             tnear = ty1;
210         if (ty2 < tfar)
211             tfar = ty2;
212     }
213     if (tnear > tfar)
214         return col;
215     if (tfar < 0.0)
216         return col;
217 
218     if (ry->d.z == 0.0) {
219         if ((ry->o.z < bx->min.z) || (ry->o.z > bx->max.z))
220             return col;
221     }
222     else {
223         tz1 = (bx->min.z - ry->o.z) / ry->d.z;
224         tz2 = (bx->max.z - ry->o.z) / ry->d.z;
225         if (tz1 > tz2) {
226             a = tz1;
227             tz1 = tz2;
228             tz2 = a;
229         }
230         if (tz1 > tnear)
231             tnear = tz1;
232         if (tz2 < tfar)
233             tfar = tz2;
234     }
235     if (tnear > tfar)
236         return col;
237     if (tfar < 0.0)
238         return col;
239 
240     if (tnear < 0.0)
241         tnear = 0.0;
242 
243     tdist = xvol->samples;
244 
245     tt = (xvol->opacity / tdist);
246 
247     bln.x = fabs(bx->min.x - bx->max.x);
248     bln.y = fabs(bx->min.y - bx->max.y);
249     bln.z = fabs(bx->min.z - bx->max.z);
250 
251     dt = 1.0 / tdist;
252     sum = 0.0;
253 
254     /* Accumulate color as the ray passes through the voxels */
255     for (t = tnear; t <= tfar; t += dt) {
256         if (sum < 1.0) {
257             pnt.x = ((ry->o.x + (ry->d.x * t)) - bx->min.x) / bln.x;
258             pnt.y = ((ry->o.y + (ry->d.y * t)) - bx->min.y) / bln.y;
259             pnt.z = ((ry->o.z + (ry->d.z * t)) - bx->min.z) / bln.z;
260 
261             /* call external evaluator assume 0.0 -> 1.0 range.. */
262             scalar = xvol->evaluator(pnt.x, pnt.y, pnt.z);
263 
264             transval = tt * scalar;
265             sum += transval;
266 
267             col2 = ExtVoxelColor(scalar);
268 
269             col.r += transval * col2.r * xvol->ambient;
270             col.g += transval * col2.g * xvol->ambient;
271             col.b += transval * col2.b * xvol->ambient;
272 
273             ddt = dt;
274 
275             /* Add in diffuse shaded light sources (no shadows) */
276             if (xvol->diffuse > 0.0) {
277                 /* Calculate the Volume gradient at the voxel */
278                 N.x = (xvol->evaluator(pnt.x - ddt, pnt.y, pnt.z) -
279                        xvol->evaluator(pnt.x + ddt, pnt.y, pnt.z)) *
280                       8.0 * tt;
281 
282                 N.y = (xvol->evaluator(pnt.x, pnt.y - ddt, pnt.z) -
283                        xvol->evaluator(pnt.x, pnt.y + ddt, pnt.z)) *
284                       8.0 * tt;
285 
286                 N.z = (xvol->evaluator(pnt.x, pnt.y, pnt.z - ddt) -
287                        xvol->evaluator(pnt.x, pnt.y, pnt.z + ddt)) *
288                       8.0 * tt;
289 
290                 /* only light surfaces with enough of a normal.. */
291                 if ((N.x * N.x + N.y * N.y + N.z * N.z) > 0.0) {
292                     diffint.r = 0.0;
293                     diffint.g = 0.0;
294                     diffint.b = 0.0;
295 
296                     /* add the contribution of each of the lights.. */
297                     for (i = 0; i < numlights; i++) {
298                         li = lightlist[i];
299                         VSUB(li->ctr, (*hit), L)
300                         VNorm(&L);
301                         VDOT(inten, N, L)
302 
303                         /* only add light if its from the front of the surface */
304                         /* could add back-lighting if we wanted to later.. */
305                         if (inten > 0.0) {
306                             diffint.r += inten * li->tex->col.r;
307                             diffint.g += inten * li->tex->col.g;
308                             diffint.b += inten * li->tex->col.b;
309                         }
310                     }
311                     col.r += col2.r * diffint.r * xvol->diffuse;
312                     col.g += col2.g * diffint.g * xvol->diffuse;
313                     col.b += col2.b * diffint.b * xvol->diffuse;
314                 }
315             }
316         }
317         else {
318             sum = 1.0;
319         }
320     }
321 
322     /* Add in transmitted ray from outside environment */
323     if (sum < 1.0) { /* spawn transmission rays / refraction */
324         color transcol;
325 
326         transcol = shade_transmission(ry, hit, 1.0 - sum);
327 
328         col.r += transcol.r; /* add the transmitted ray  */
329         col.g += transcol.g; /* to the diffuse and       */
330         col.b += transcol.b; /* transmission total..     */
331     }
332 
333     return col;
334 }
335