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  * vol.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 "vol.hpp"
58 #include "box.hpp"
59 #include "trace.hpp"
60 #include "ui.hpp"
61 #include "light.hpp"
62 #include "shade.hpp"
63 
64 int scalarvol_bbox(void *obj, vector *min, vector *max) {
65     box *b = (box *)obj;
66 
67     *min = b->min;
68     *max = b->max;
69 
70     return 1;
71 }
72 
73 void *newscalarvol(void *intex,
74                    vector min,
75                    vector max,
76                    int xs,
77                    int ys,
78                    int zs,
79                    char *fname,
80                    scalarvol *invol) {
81     box *bx;
82     texture *tx, *tex;
83     scalarvol *vol;
84 
85     tex = (texture *)intex;
86     tex->shadowcast = 0; /* doesn't cast a shadow */
87 
88     tx = (texture *)rt_getmem(sizeof(texture));
89 
90     /* is the volume data already loaded? */
91     if (invol == nullptr) {
92         vol = (scalarvol *)rt_getmem(sizeof(scalarvol));
93         vol->loaded = 0;
94         vol->data = nullptr;
95     }
96     else
97         vol = invol;
98 
99     vol->opacity = tex->opacity;
100     vol->xres = xs;
101     vol->yres = ys;
102     vol->zres = zs;
103     strcpy(vol->name, fname);
104 
105     tx->ctr.x = 0.0;
106     tx->ctr.y = 0.0;
107     tx->ctr.z = 0.0;
108     tx->rot = tx->ctr;
109     tx->scale = tx->ctr;
110     tx->uaxs = tx->ctr;
111     tx->vaxs = tx->ctr;
112 
113     tx->islight = 0;
114     tx->shadowcast = 0; /* doesn't cast a shadow */
115 
116     tx->col = tex->col;
117     tx->ambient = 1.0;
118     tx->diffuse = 0.0;
119     tx->specular = 0.0;
120     tx->opacity = 1.0;
121     tx->img = vol;
122     tx->texfunc = (color(*)(void *, void *, void *))(scalar_volume_texture);
123 
124     bx = newbox(tx, min, max);
125     tx->obj = (void *)bx; /* XXX hack! */
126 
127     return (void *)bx;
128 }
129 
130 color VoxelColor(flt scalar) {
131     color col;
132 
133     if (scalar > 1.0)
134         scalar = 1.0;
135 
136     if (scalar < 0.0)
137         scalar = 0.0;
138 
139     if (scalar < 0.25) {
140         col.r = scalar * 4.0;
141         col.g = 0.0;
142         col.b = 0.0;
143     }
144     else {
145         if (scalar < 0.75) {
146             col.r = 1.0;
147             col.g = (scalar - 0.25) * 2.0;
148             col.b = 0.0;
149         }
150         else {
151             col.r = 1.0;
152             col.g = 1.0;
153             col.b = (scalar - 0.75) * 4.0;
154         }
155     }
156 
157     return col;
158 }
159 
160 color scalar_volume_texture(vector *hit, texture *tex, ray *ry) {
161     color col, col2;
162     box *bx;
163     flt a, tx1, tx2, ty1, ty2, tz1, tz2;
164     flt tnear, tfar;
165     flt t, tdist, dt, sum, tt;
166     vector pnt, bln;
167     scalarvol *vol;
168     flt scalar, transval;
169     int x, y, z;
170     unsigned char *ptr;
171 
172     bx = (box *)tex->obj;
173     vol = (scalarvol *)bx->tex->img;
174 
175     col.r = 0.0;
176     col.g = 0.0;
177     col.b = 0.0;
178 
179     tnear = -FHUGE;
180     tfar = FHUGE;
181 
182     if (ry->d.x == 0.0) {
183         if ((ry->o.x < bx->min.x) || (ry->o.x > bx->max.x))
184             return col;
185     }
186     else {
187         tx1 = (bx->min.x - ry->o.x) / ry->d.x;
188         tx2 = (bx->max.x - ry->o.x) / ry->d.x;
189         if (tx1 > tx2) {
190             a = tx1;
191             tx1 = tx2;
192             tx2 = a;
193         }
194         if (tx1 > tnear)
195             tnear = tx1;
196         if (tx2 < tfar)
197             tfar = tx2;
198     }
199     if (tnear > tfar)
200         return col;
201     if (tfar < 0.0)
202         return col;
203 
204     if (ry->d.y == 0.0) {
205         if ((ry->o.y < bx->min.y) || (ry->o.y > bx->max.y))
206             return col;
207     }
208     else {
209         ty1 = (bx->min.y - ry->o.y) / ry->d.y;
210         ty2 = (bx->max.y - ry->o.y) / ry->d.y;
211         if (ty1 > ty2) {
212             a = ty1;
213             ty1 = ty2;
214             ty2 = a;
215         }
216         if (ty1 > tnear)
217             tnear = ty1;
218         if (ty2 < tfar)
219             tfar = ty2;
220     }
221     if (tnear > tfar)
222         return col;
223     if (tfar < 0.0)
224         return col;
225 
226     if (ry->d.z == 0.0) {
227         if ((ry->o.z < bx->min.z) || (ry->o.z > bx->max.z))
228             return col;
229     }
230     else {
231         tz1 = (bx->min.z - ry->o.z) / ry->d.z;
232         tz2 = (bx->max.z - ry->o.z) / ry->d.z;
233         if (tz1 > tz2) {
234             a = tz1;
235             tz1 = tz2;
236             tz2 = a;
237         }
238         if (tz1 > tnear)
239             tnear = tz1;
240         if (tz2 < tfar)
241             tfar = tz2;
242     }
243     if (tnear > tfar)
244         return col;
245     if (tfar < 0.0)
246         return col;
247 
248     if (tnear < 0.0)
249         tnear = 0.0;
250 
251     tdist = sqrt((flt)(vol->xres * vol->xres + vol->yres * vol->yres + vol->zres * vol->zres));
252     tt = (vol->opacity / tdist);
253 
254     bln.x = fabs(bx->min.x - bx->max.x);
255     bln.y = fabs(bx->min.y - bx->max.y);
256     bln.z = fabs(bx->min.z - bx->max.z);
257 
258     dt = sqrt(bln.x * bln.x + bln.y * bln.y + bln.z * bln.z) / tdist;
259     sum = 0.0;
260 
261     /* move the volume residency check out of loop.. */
262     if (!vol->loaded) {
263         LoadVol(vol);
264         vol->loaded = 1;
265     }
266 
267     for (t = tnear; t <= tfar; t += dt) {
268         pnt.x = ((ry->o.x + (ry->d.x * t)) - bx->min.x) / bln.x;
269         pnt.y = ((ry->o.y + (ry->d.y * t)) - bx->min.y) / bln.y;
270         pnt.z = ((ry->o.z + (ry->d.z * t)) - bx->min.z) / bln.z;
271 
272         x = (int)((vol->xres - 1.5) * pnt.x + 0.5);
273         y = (int)((vol->yres - 1.5) * pnt.y + 0.5);
274         z = (int)((vol->zres - 1.5) * pnt.z + 0.5);
275 
276         ptr = vol->data + ((vol->xres * vol->yres * z) + (vol->xres * y) + x);
277 
278         scalar = (flt)((flt)1.0 * ((int)ptr[0])) / 255.0;
279 
280         sum += tt * scalar;
281 
282         transval = tt * scalar;
283 
284         col2 = VoxelColor(scalar);
285 
286         if (sum < 1.0) {
287             col.r += transval * col2.r;
288             col.g += transval * col2.g;
289             col.b += transval * col2.b;
290             if (sum < 0.0)
291                 sum = 0.0;
292         }
293         else {
294             sum = 1.0;
295         }
296     }
297 
298     if (sum < 1.0) { /* spawn transmission rays / refraction */
299         color transcol;
300 
301         transcol = shade_transmission(ry, hit, 1.0 - sum);
302 
303         col.r += transcol.r; /* add the transmitted ray  */
304         col.g += transcol.g; /* to the diffuse and       */
305         col.b += transcol.b; /* transmission total..     */
306     }
307 
308     return col;
309 }
310 
311 void LoadVol(scalarvol *vol) {
312     FILE *dfile;
313     std::size_t status;
314     char msgtxt[2048];
315 
316     dfile = fopen(vol->name, "r");
317     if (dfile == nullptr) {
318         char msgtxt[2048];
319         sprintf(msgtxt, "Vol: can't open %s for input!!! Aborting\n", vol->name);
320         rt_ui_message(MSG_ERR, msgtxt);
321         rt_ui_message(MSG_ABORT, "Rendering Aborted.");
322         std::exit(-1);
323     }
324 
325     sprintf(
326         msgtxt, "loading %dx%dx%d volume set from %s", vol->xres, vol->yres, vol->zres, vol->name);
327     rt_ui_message(MSG_0, msgtxt);
328 
329     vol->data = (unsigned char *)rt_getmem(vol->xres * vol->yres * vol->zres);
330 
331     status = fread(vol->data, 1, (vol->xres * vol->yres * vol->zres), dfile);
332     fclose(dfile);
333 }
334