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