1 /*
2 Copyright (c) 2005-2020 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 * grid.cpp - spatial subdivision efficiency structures
48 */
49
50 #include "machine.h"
51 #include "types.h"
52 #include "macros.h"
53 #include "vector.h"
54 #include "intersect.h"
55 #include "util.h"
56
57 #define GRID_PRIVATE
58 #include "grid.h"
59
60 #ifndef cbrt
61 #define cbrt(x) ((x) > 0.0 ? pow((double)(x), 1.0/3.0) : \
62 ((x) < 0.0 ? -pow((double)-(x), 1.0/3.0) : 0.0))
63
64 #define qbrt(x) ((x) > 0.0 ? pow((double)(x), 1.0/4.0) : \
65 ((x) < 0.0 ? -pow((double)-(x), 1.0/4.0) : 0.0))
66
67 #endif
68
69 static object_methods grid_methods = {
70 (void (*)(void *, void *))(grid_intersect),
71 (void (*)(void *, void *, void *, void *))(NULL),
72 grid_bbox,
73 grid_free
74 };
75
76 extern bool silent_mode;
77
newgrid(int xsize,int ysize,int zsize,vector min,vector max)78 object * newgrid(int xsize, int ysize, int zsize, vector min, vector max) {
79 grid * g;
80
81 g = (grid *) rt_getmem(sizeof(grid));
82 memset(g, 0, sizeof(grid));
83
84 g->methods = &grid_methods;
85 g->id = new_objectid();
86
87 g->xsize = xsize;
88 g->ysize = ysize;
89 g->zsize = zsize;
90
91 g->min = min;
92 g->max = max;
93
94 VSub(&g->max, &g->min, &g->voxsize);
95 g->voxsize.x /= (flt) g->xsize;
96 g->voxsize.y /= (flt) g->ysize;
97 g->voxsize.z /= (flt) g->zsize;
98
99 g->cells = (objectlist **) rt_getmem(xsize*ysize*zsize*sizeof(objectlist *));
100 memset(g->cells, 0, xsize*ysize*zsize * sizeof(objectlist *));
101
102 /* fprintf(stderr, "New grid, size: %8d %8d %8d\n", g->xsize, g->ysize, g->zsize); */
103
104 return (object *) g;
105 }
106
grid_bbox(void * obj,vector * min,vector * max)107 static int grid_bbox(void * obj, vector * min, vector * max) {
108 grid * g = (grid *) obj;
109
110 *min = g->min;
111 *max = g->max;
112
113 return 1;
114 }
115
grid_free(void * v)116 static void grid_free(void * v) {
117 int i, numvoxels;
118 grid * g = (grid *) v;
119
120 /* loop through all voxels and free the object lists */
121 numvoxels = g->xsize * g->ysize * g->zsize;
122 for (i=0; i<numvoxels; i++) {
123 objectlist * lcur, * lnext;
124
125 lcur = g->cells[i];
126 while (lcur != NULL) {
127 lnext = lcur->next;
128 free(lcur);
129 }
130 }
131
132 /* free the grid cells */
133 free(g->cells);
134
135 /* free all objects on the grid object list */
136 free_objects(g->objects);
137
138 free(g);
139 }
140
globalbound(object ** rootlist,vector * gmin,vector * gmax)141 static void globalbound(object ** rootlist, vector * gmin, vector * gmax) {
142 vector min, max;
143 object * cur;
144
145 if (*rootlist == NULL) /* don't bound non-existent objects */
146 return;
147
148 gmin->x = FHUGE; gmin->y = FHUGE; gmin->z = FHUGE;
149 gmax->x = -FHUGE; gmax->y = -FHUGE; gmax->z = -FHUGE;
150
151 cur=*rootlist;
152 while (cur != NULL) { /* Go! */
153 min.x = -FHUGE; min.y = -FHUGE; min.z = -FHUGE;
154 max.x = FHUGE; max.y = FHUGE; max.z = FHUGE;
155
156 if (cur->methods->bbox((void *) cur, &min, &max)) {
157 gmin->x = MYMIN( gmin->x , min.x);
158 gmin->y = MYMIN( gmin->y , min.y);
159 gmin->z = MYMIN( gmin->z , min.z);
160
161 gmax->x = MYMAX( gmax->x , max.x);
162 gmax->y = MYMAX( gmax->y , max.y);
163 gmax->z = MYMAX( gmax->z , max.z);
164 }
165
166 cur=(object *)cur->nextobj;
167 }
168 }
169
170
cellbound(grid * g,gridindex * index,vector * cmin,vector * cmax)171 static int cellbound(grid *g, gridindex *index, vector * cmin, vector * cmax) {
172 vector min, max, cellmin, cellmax;
173 objectlist * cur;
174 int numinbounds = 0;
175
176 cur = g->cells[index->z*g->xsize*g->ysize + index->y*g->xsize + index->x];
177
178 if (cur == NULL) /* don't bound non-existent objects */
179 return 0;
180
181 cellmin.x = voxel2x(g, index->x);
182 cellmin.y = voxel2y(g, index->y);
183 cellmin.z = voxel2z(g, index->z);
184
185 cellmax.x = cellmin.x + g->voxsize.x;
186 cellmax.y = cellmin.y + g->voxsize.y;
187 cellmax.z = cellmin.z + g->voxsize.z;
188
189 cmin->x = FHUGE; cmin->y = FHUGE; cmin->z = FHUGE;
190 cmax->x = -FHUGE; cmax->y = -FHUGE; cmax->z = -FHUGE;
191
192 while (cur != NULL) { /* Go! */
193 min.x = -FHUGE; min.y = -FHUGE; min.z = -FHUGE;
194 max.x = FHUGE; max.y = FHUGE; max.z = FHUGE;
195
196 if (cur->obj->methods->bbox((void *) cur->obj, &min, &max)) {
197 if ((min.x >= cellmin.x) && (max.x <= cellmax.x) &&
198 (min.y >= cellmin.y) && (max.y <= cellmax.y) &&
199 (min.z >= cellmin.z) && (max.z <= cellmax.z)) {
200
201 cmin->x = MYMIN( cmin->x , min.x);
202 cmin->y = MYMIN( cmin->y , min.y);
203 cmin->z = MYMIN( cmin->z , min.z);
204
205 cmax->x = MYMAX( cmax->x , max.x);
206 cmax->y = MYMAX( cmax->y , max.y);
207 cmax->z = MYMAX( cmax->z , max.z);
208
209 numinbounds++;
210 }
211 }
212
213 cur=cur->next;
214 }
215
216 /* in case we get a 0.0 sized axis on the cell bounds, we'll */
217 /* use the original cell bounds */
218 if ((cmax->x - cmin->x) < EPSILON) {
219 cmax->x += EPSILON;
220 cmin->x -= EPSILON;
221 }
222 if ((cmax->y - cmin->y) < EPSILON) {
223 cmax->y += EPSILON;
224 cmin->y -= EPSILON;
225 }
226 if ((cmax->z - cmin->z) < EPSILON) {
227 cmax->z += EPSILON;
228 cmin->z -= EPSILON;
229 }
230
231 return numinbounds;
232 }
233
countobj(object * root)234 static int countobj(object * root) {
235 object * cur; /* counts the number of objects on a list */
236 int numobj;
237
238 numobj=0;
239 cur=root;
240
241 while (cur != NULL) {
242 cur=(object *)cur->nextobj;
243 numobj++;
244 }
245 return numobj;
246 }
247
countobjlist(objectlist * root)248 static int countobjlist(objectlist * root) {
249 objectlist * cur;
250 int numobj;
251
252 numobj=0;
253 cur = root;
254
255 while (cur != NULL) {
256 cur = cur->next;
257 numobj++;
258 }
259 return numobj;
260 }
261
engrid_scene(object ** list)262 int engrid_scene(object ** list) {
263 grid * g;
264 int numobj, numcbrt;
265 vector gmin, gmax;
266 gridindex index;
267
268 if (*list == NULL)
269 return 0;
270
271 numobj = countobj(*list);
272
273 if ( !silent_mode )
274 fprintf(stderr, "Scene contains %d bounded objects.\n", numobj);
275
276 if (numobj > 16) {
277 numcbrt = (int) cbrt(4*numobj);
278 globalbound(list, &gmin, &gmax);
279
280 g = (grid *) newgrid(numcbrt, numcbrt, numcbrt, gmin, gmax);
281 engrid_objlist(g, list);
282
283 numobj = countobj(*list);
284 g->nextobj = *list;
285 *list = (object *) g;
286
287 /* now create subgrids.. */
288 for (index.z=0; index.z<g->zsize; index.z++) {
289 for (index.y=0; index.y<g->ysize; index.y++) {
290 for (index.x=0; index.x<g->xsize; index.x++) {
291 engrid_cell(g, &index);
292 }
293 }
294 }
295 }
296
297 return 1;
298 }
299
300
engrid_objlist(grid * g,object ** list)301 void engrid_objlist(grid * g, object ** list) {
302 object * cur, * next, **prev;
303
304 if (*list == NULL)
305 return;
306
307 prev = list;
308 cur = *list;
309
310 while (cur != NULL) {
311 next = (object *)cur->nextobj;
312
313 if (engrid_object(g, cur))
314 *prev = next;
315 else
316 prev = (object **) &cur->nextobj;
317
318 cur = next;
319 }
320 }
321
engrid_cell(grid * gold,gridindex * index)322 static int engrid_cell(grid * gold, gridindex *index) {
323 vector gmin, gmax, gsize;
324 flt len;
325 int numobj, numcbrt, xs, ys, zs;
326 grid * g;
327 objectlist **list;
328 objectlist * newobj;
329
330 list = &gold->cells[index->z*gold->xsize*gold->ysize +
331 index->y*gold->xsize + index->x];
332
333 if (*list == NULL)
334 return 0;
335
336 numobj = cellbound(gold, index, &gmin, &gmax);
337
338 VSub(&gmax, &gmin, &gsize);
339 len = 1.0 / (MYMAX( MYMAX(gsize.x, gsize.y), gsize.z ));
340 gsize.x *= len;
341 gsize.y *= len;
342 gsize.z *= len;
343
344 if (numobj > 16) {
345 numcbrt = (int) cbrt(2*numobj);
346
347 xs = (int) ((flt) numcbrt * gsize.x);
348 if (xs < 1) xs = 1;
349 ys = (int) ((flt) numcbrt * gsize.y);
350 if (ys < 1) ys = 1;
351 zs = (int) ((flt) numcbrt * gsize.z);
352 if (zs < 1) zs = 1;
353
354 g = (grid *) newgrid(xs, ys, zs, gmin, gmax);
355 engrid_objectlist(g, list);
356
357 newobj = (objectlist *) rt_getmem(sizeof(objectlist));
358 newobj->obj = (object *) g;
359 newobj->next = *list;
360 *list = newobj;
361
362 g->nextobj = gold->objects;
363 gold->objects = (object *) g;
364 }
365
366 return 1;
367 }
368
engrid_objectlist(grid * g,objectlist ** list)369 static int engrid_objectlist(grid * g, objectlist ** list) {
370 objectlist * cur, * next, **prev;
371 int numsucceeded = 0;
372
373 if (*list == NULL)
374 return 0;
375
376 prev = list;
377 cur = *list;
378
379 while (cur != NULL) {
380 next = cur->next;
381
382 if (engrid_object(g, cur->obj)) {
383 *prev = next;
384 free(cur);
385 numsucceeded++;
386 }
387 else {
388 prev = &cur->next;
389 }
390
391 cur = next;
392 }
393
394 return numsucceeded;
395 }
396
397
398
engrid_object(grid * g,object * obj)399 static int engrid_object(grid * g, object * obj) {
400 vector omin, omax;
401 gridindex low, high;
402 int x, y, z, zindex, yindex, voxindex;
403 objectlist * tmp;
404
405 if (obj->methods->bbox(obj, &omin, &omax)) {
406 if (!pos2grid(g, &omin, &low) || !pos2grid(g, &omax, &high)) {
407 return 0; /* object is not wholly contained in the grid */
408 }
409 }
410 else {
411 return 0; /* object is unbounded */
412 }
413
414 /* add the object to the complete list of objects in the grid */
415 obj->nextobj = g->objects;
416 g->objects = obj;
417
418 /* add this object to all voxels it inhabits */
419 for (z=low.z; z<=high.z; z++) {
420 zindex = z * g->xsize * g->ysize;
421 for (y=low.y; y<=high.y; y++) {
422 yindex = y * g->xsize;
423 for (x=low.x; x<=high.x; x++) {
424 voxindex = x + yindex + zindex;
425 tmp = (objectlist *) rt_getmem(sizeof(objectlist));
426 tmp->next = g->cells[voxindex];
427 tmp->obj = obj;
428 g->cells[voxindex] = tmp;
429 }
430 }
431 }
432
433 return 1;
434 }
435
pos2grid(grid * g,vector * pos,gridindex * index)436 static int pos2grid(grid * g, vector * pos, gridindex * index) {
437 index->x = (int) ((pos->x - g->min.x) / g->voxsize.x);
438 index->y = (int) ((pos->y - g->min.y) / g->voxsize.y);
439 index->z = (int) ((pos->z - g->min.z) / g->voxsize.z);
440
441 if (index->x == g->xsize)
442 index->x--;
443 if (index->y == g->ysize)
444 index->y--;
445 if (index->z == g->zsize)
446 index->z--;
447
448 if (index->x < 0 || index->x > g->xsize ||
449 index->y < 0 || index->y > g->ysize ||
450 index->z < 0 || index->z > g->zsize)
451 return 0;
452
453 if (pos->x < g->min.x || pos->x > g->max.x ||
454 pos->y < g->min.y || pos->y > g->max.y ||
455 pos->z < g->min.z || pos->z > g->max.z)
456 return 0;
457
458 return 1;
459 }
460
461
462 /* the real thing */
grid_intersect(grid * g,ray * ry)463 static void grid_intersect(grid * g, ray * ry) {
464 flt tnear, tfar, offset;
465 vector curpos, tmax, tdelta, pdeltaX, pdeltaY, pdeltaZ, nXp, nYp, nZp;
466 gridindex curvox, step, out;
467 int voxindex;
468 objectlist * cur;
469
470 if (ry->flags & RT_RAY_FINISHED)
471 return;
472
473 if (!grid_bounds_intersect(g, ry, &tnear, &tfar))
474 return;
475
476 if (ry->maxdist < tnear)
477 return;
478
479 curpos = Raypnt(ry, tnear);
480 pos2grid(g, &curpos, &curvox);
481 offset = tnear;
482
483 /* Setup X iterator stuff */
484 if (fabs(ry->d.x) < EPSILON) {
485 tmax.x = FHUGE;
486 tdelta.x = 0.0;
487 step.x = 0;
488 out.x = 0; /* never goes out of bounds on this axis */
489 }
490 else if (ry->d.x < 0.0) {
491 tmax.x = offset + ((voxel2x(g, curvox.x) - curpos.x) / ry->d.x);
492 tdelta.x = g->voxsize.x / - ry->d.x;
493 step.x = out.x = -1;
494 }
495 else {
496 tmax.x = offset + ((voxel2x(g, curvox.x + 1) - curpos.x) / ry->d.x);
497 tdelta.x = g->voxsize.x / ry->d.x;
498 step.x = 1;
499 out.x = g->xsize;
500 }
501
502 /* Setup Y iterator stuff */
503 if (fabs(ry->d.y) < EPSILON) {
504 tmax.y = FHUGE;
505 tdelta.y = 0.0;
506 step.y = 0;
507 out.y = 0; /* never goes out of bounds on this axis */
508 }
509 else if (ry->d.y < 0.0) {
510 tmax.y = offset + ((voxel2y(g, curvox.y) - curpos.y) / ry->d.y);
511 tdelta.y = g->voxsize.y / - ry->d.y;
512 step.y = out.y = -1;
513 }
514 else {
515 tmax.y = offset + ((voxel2y(g, curvox.y + 1) - curpos.y) / ry->d.y);
516 tdelta.y = g->voxsize.y / ry->d.y;
517 step.y = 1;
518 out.y = g->ysize;
519 }
520
521 /* Setup Z iterator stuff */
522 if (fabs(ry->d.z) < EPSILON) {
523 tmax.z = FHUGE;
524 tdelta.z = 0.0;
525 step.z = 0;
526 out.z = 0; /* never goes out of bounds on this axis */
527 }
528 else if (ry->d.z < 0.0) {
529 tmax.z = offset + ((voxel2z(g, curvox.z) - curpos.z) / ry->d.z);
530 tdelta.z = g->voxsize.z / - ry->d.z;
531 step.z = out.z = -1;
532 }
533 else {
534 tmax.z = offset + ((voxel2z(g, curvox.z + 1) - curpos.z) / ry->d.z);
535 tdelta.z = g->voxsize.z / ry->d.z;
536 step.z = 1;
537 out.z = g->zsize;
538 }
539
540 pdeltaX = ry->d;
541 VScale(&pdeltaX, tdelta.x);
542 pdeltaY = ry->d;
543 VScale(&pdeltaY, tdelta.y);
544 pdeltaZ = ry->d;
545 VScale(&pdeltaZ, tdelta.z);
546
547 nXp = Raypnt(ry, tmax.x);
548 nYp = Raypnt(ry, tmax.y);
549 nZp = Raypnt(ry, tmax.z);
550
551 voxindex = curvox.z*g->xsize*g->ysize + curvox.y*g->xsize + curvox.x;
552 while (1) {
553 if (tmax.x < tmax.y && tmax.x < tmax.z) {
554 cur = g->cells[voxindex];
555 while (cur != NULL) {
556 if (ry->mbox[cur->obj->id] != ry->serial) {
557 ry->mbox[cur->obj->id] = ry->serial;
558 cur->obj->methods->intersect(cur->obj, ry);
559 }
560 cur = cur->next;
561 }
562 curvox.x += step.x;
563 if (ry->maxdist < tmax.x || curvox.x == out.x)
564 break;
565 voxindex += step.x;
566 tmax.x += tdelta.x;
567 curpos = nXp;
568 nXp.x += pdeltaX.x;
569 nXp.y += pdeltaX.y;
570 nXp.z += pdeltaX.z;
571 }
572 else if (tmax.z < tmax.y) {
573 cur = g->cells[voxindex];
574 while (cur != NULL) {
575 if (ry->mbox[cur->obj->id] != ry->serial) {
576 ry->mbox[cur->obj->id] = ry->serial;
577 cur->obj->methods->intersect(cur->obj, ry);
578 }
579 cur = cur->next;
580 }
581 curvox.z += step.z;
582 if (ry->maxdist < tmax.z || curvox.z == out.z)
583 break;
584 voxindex += step.z*g->xsize*g->ysize;
585 tmax.z += tdelta.z;
586 curpos = nZp;
587 nZp.x += pdeltaZ.x;
588 nZp.y += pdeltaZ.y;
589 nZp.z += pdeltaZ.z;
590 }
591 else {
592 cur = g->cells[voxindex];
593 while (cur != NULL) {
594 if (ry->mbox[cur->obj->id] != ry->serial) {
595 ry->mbox[cur->obj->id] = ry->serial;
596 cur->obj->methods->intersect(cur->obj, ry);
597 }
598 cur = cur->next;
599 }
600 curvox.y += step.y;
601 if (ry->maxdist < tmax.y || curvox.y == out.y)
602 break;
603 voxindex += step.y*g->xsize;
604 tmax.y += tdelta.y;
605 curpos = nYp;
606 nYp.x += pdeltaY.x;
607 nYp.y += pdeltaY.y;
608 nYp.z += pdeltaY.z;
609 }
610
611 if (ry->flags & RT_RAY_FINISHED)
612 break;
613 }
614 }
615
voxel_intersect(grid * g,ray * ry,int voxindex)616 static void voxel_intersect(grid * g, ray * ry, int voxindex) {
617 objectlist * cur;
618
619 cur = g->cells[voxindex];
620 while (cur != NULL) {
621 cur->obj->methods->intersect(cur->obj, ry);
622 cur = cur->next;
623 }
624 }
625
grid_bounds_intersect(grid * g,ray * ry,flt * nr,flt * fr)626 static int grid_bounds_intersect(grid * g, ray * ry, flt *nr, flt *fr) {
627 flt a, tx1, tx2, ty1, ty2, tz1, tz2;
628 flt tnear, tfar;
629
630 tnear= -FHUGE;
631 tfar= FHUGE;
632
633 if (ry->d.x == 0.0) {
634 if ((ry->o.x < g->min.x) || (ry->o.x > g->max.x)) return 0;
635 }
636 else {
637 tx1 = (g->min.x - ry->o.x) / ry->d.x;
638 tx2 = (g->max.x - ry->o.x) / ry->d.x;
639 if (tx1 > tx2) { a=tx1; tx1=tx2; tx2=a; }
640 if (tx1 > tnear) tnear=tx1;
641 if (tx2 < tfar) tfar=tx2;
642 }
643 if (tnear > tfar) return 0;
644 if (tfar < 0.0) return 0;
645
646 if (ry->d.y == 0.0) {
647 if ((ry->o.y < g->min.y) || (ry->o.y > g->max.y)) return 0;
648 }
649 else {
650 ty1 = (g->min.y - ry->o.y) / ry->d.y;
651 ty2 = (g->max.y - ry->o.y) / ry->d.y;
652 if (ty1 > ty2) { a=ty1; ty1=ty2; ty2=a; }
653 if (ty1 > tnear) tnear=ty1;
654 if (ty2 < tfar) tfar=ty2;
655 }
656 if (tnear > tfar) return 0;
657 if (tfar < 0.0) return 0;
658
659 if (ry->d.z == 0.0) {
660 if ((ry->o.z < g->min.z) || (ry->o.z > g->max.z)) return 0;
661 }
662 else {
663 tz1 = (g->min.z - ry->o.z) / ry->d.z;
664 tz2 = (g->max.z - ry->o.z) / ry->d.z;
665 if (tz1 > tz2) { a=tz1; tz1=tz2; tz2=a; }
666 if (tz1 > tnear) tnear=tz1;
667 if (tz2 < tfar) tfar=tz2;
668 }
669 if (tnear > tfar) return 0;
670 if (tfar < 0.0) return 0;
671
672 *nr = tnear;
673 *fr = tfar;
674 return 1;
675 }
676