GHOST: Fix uninitialized var
[blender.git] / source / blender / render / intern / raytrace / rayobject_octree.cpp
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * The Original Code is Copyright (C) 1990-1998 NeoGeo BV.
19  * All rights reserved.
20  *
21  * Contributors: 2004/2005 Blender Foundation, full recode
22  *
23  * ***** END GPL LICENSE BLOCK *****
24  */
25
26 /** \file blender/render/intern/raytrace/rayobject_octree.cpp
27  *  \ingroup render
28  */
29
30
31 /* IMPORTANT NOTE: this code must be independent of any other render code
32  * to use it outside the renderer! */
33
34 #include <math.h>
35 #include <string.h>
36 #include <stdlib.h>
37 #include <float.h>
38 #include <assert.h>
39
40 #include "MEM_guardedalloc.h"
41
42 #include "DNA_material_types.h"
43
44 #include "BLI_math.h"
45 #include "BLI_utildefines.h"
46
47 #include "rayintersection.h"
48 #include "rayobject.h"
49
50 /* ********** structs *************** */
51 #define BRANCH_ARRAY 1024
52 #define NODE_ARRAY 4096
53
54 typedef struct Branch {
55         struct Branch *b[8];
56 } Branch;
57
58 typedef struct OcVal {
59         short ocx, ocy, ocz;
60 } OcVal;
61
62 typedef struct Node {
63         struct RayFace *v[8];
64         struct OcVal ov[8];
65         struct Node *next;
66 } Node;
67
68 typedef struct Octree {
69         RayObject rayobj;
70
71         struct Branch **adrbranch;
72         struct Node **adrnode;
73         float ocsize;   /* ocsize: mult factor,  max size octree */
74         float ocfacx, ocfacy, ocfacz;
75         float min[3], max[3];
76         int ocres;
77         int branchcount, nodecount;
78
79         /* during building only */
80         char *ocface;
81
82         RayFace **ro_nodes;
83         int ro_nodes_size, ro_nodes_used;
84
85 } Octree;
86
87 static int  RE_rayobject_octree_intersect(RayObject *o, Isect *isec);
88 static void RE_rayobject_octree_add(RayObject *o, RayObject *ob);
89 static void RE_rayobject_octree_done(RayObject *o);
90 static void RE_rayobject_octree_free(RayObject *o);
91 static void RE_rayobject_octree_bb(RayObject *o, float *min, float *max);
92
93 /*
94  * This function is not expected to be called by current code state.
95  */
96 static float RE_rayobject_octree_cost(RayObject *UNUSED(o))
97 {
98         return 1.0;
99 }
100
101 static void RE_rayobject_octree_hint_bb(RayObject *UNUSED(o), RayHint *UNUSED(hint),
102                                         float *UNUSED(min), float *UNUSED(max))
103 {
104         return;
105 }
106
107 static RayObjectAPI octree_api =
108 {
109         RE_rayobject_octree_intersect,
110         RE_rayobject_octree_add,
111         RE_rayobject_octree_done,
112         RE_rayobject_octree_free,
113         RE_rayobject_octree_bb,
114         RE_rayobject_octree_cost,
115         RE_rayobject_octree_hint_bb
116 };
117
118 /* **************** ocval method ******************* */
119 /* within one octree node, a set of 3x15 bits defines a 'boundbox' to OR with */
120
121 #define OCVALRES    15
122 #define BROW16(min, max) \
123         (((max) >= OCVALRES ? 0xFFFF : (1 << ((max) + 1)) - 1) - (((min) > 0) ? ((1 << (min)) - 1) : 0))
124
125 static void calc_ocval_face(float *v1, float *v2, float *v3, float *v4, short x, short y, short z, OcVal *ov)
126 {
127         float min[3], max[3];
128         int ocmin, ocmax;
129
130         copy_v3_v3(min, v1);
131         copy_v3_v3(max, v1);
132         DO_MINMAX(v2, min, max);
133         DO_MINMAX(v3, min, max);
134         if (v4) {
135                 DO_MINMAX(v4, min, max);
136         }
137
138         ocmin = OCVALRES * (min[0] - x);
139         ocmax = OCVALRES * (max[0] - x);
140         ov->ocx = BROW16(ocmin, ocmax);
141
142         ocmin = OCVALRES * (min[1] - y);
143         ocmax = OCVALRES * (max[1] - y);
144         ov->ocy = BROW16(ocmin, ocmax);
145
146         ocmin = OCVALRES * (min[2] - z);
147         ocmax = OCVALRES * (max[2] - z);
148         ov->ocz = BROW16(ocmin, ocmax);
149
150 }
151
152 static void calc_ocval_ray(OcVal *ov, float xo, float yo, float zo, float *vec1, float *vec2)
153 {
154         int ocmin, ocmax;
155
156         if (vec1[0] < vec2[0]) {
157                 ocmin = OCVALRES * (vec1[0] - xo);
158                 ocmax = OCVALRES * (vec2[0] - xo);
159         }
160         else {
161                 ocmin = OCVALRES * (vec2[0] - xo);
162                 ocmax = OCVALRES * (vec1[0] - xo);
163         }
164         ov->ocx = BROW16(ocmin, ocmax);
165
166         if (vec1[1] < vec2[1]) {
167                 ocmin = OCVALRES * (vec1[1] - yo);
168                 ocmax = OCVALRES * (vec2[1] - yo);
169         }
170         else {
171                 ocmin = OCVALRES * (vec2[1] - yo);
172                 ocmax = OCVALRES * (vec1[1] - yo);
173         }
174         ov->ocy = BROW16(ocmin, ocmax);
175
176         if (vec1[2] < vec2[2]) {
177                 ocmin = OCVALRES * (vec1[2] - zo);
178                 ocmax = OCVALRES * (vec2[2] - zo);
179         }
180         else {
181                 ocmin = OCVALRES * (vec2[2] - zo);
182                 ocmax = OCVALRES * (vec1[2] - zo);
183         }
184         ov->ocz = BROW16(ocmin, ocmax);
185 }
186
187 /* ************* octree ************** */
188
189 static Branch *addbranch(Octree *oc, Branch *br, short ocb)
190 {
191         int index;
192
193         if (br->b[ocb]) return br->b[ocb];
194
195         oc->branchcount++;
196         index = oc->branchcount >> 12;
197
198         if (oc->adrbranch[index] == NULL)
199                 oc->adrbranch[index] = (Branch *)MEM_callocN(4096 * sizeof(Branch), "new oc branch");
200
201         if (oc->branchcount >= BRANCH_ARRAY * 4096) {
202                 printf("error; octree branches full\n");
203                 oc->branchcount = 0;
204         }
205
206         return br->b[ocb] = oc->adrbranch[index] + (oc->branchcount & 4095);
207 }
208
209 static Node *addnode(Octree *oc)
210 {
211         int index;
212
213         oc->nodecount++;
214         index = oc->nodecount >> 12;
215
216         if (oc->adrnode[index] == NULL)
217                 oc->adrnode[index] = (Node *)MEM_callocN(4096 * sizeof(Node), "addnode");
218
219         if (oc->nodecount > NODE_ARRAY * NODE_ARRAY) {
220                 printf("error; octree nodes full\n");
221                 oc->nodecount = 0;
222         }
223
224         return oc->adrnode[index] + (oc->nodecount & 4095);
225 }
226
227 static bool face_in_node(RayFace *face, short x, short y, short z, float rtf[4][3])
228 {
229         static float nor[3], d;
230         float fx, fy, fz;
231
232         // init static vars
233         if (face) {
234                 normal_tri_v3(nor, rtf[0], rtf[1], rtf[2]);
235                 d = -nor[0] * rtf[0][0] - nor[1] * rtf[0][1] - nor[2] * rtf[0][2];
236                 return 0;
237         }
238
239         fx = x;
240         fy = y;
241         fz = z;
242
243         if ((fx) * nor[0] + (fy) * nor[1] + (fz) * nor[2] + d > 0.0f) {
244                 if ((fx + 1) * nor[0] + (fy    ) * nor[1] + (fz    ) * nor[2] + d < 0.0f) return 1;
245                 if ((fx    ) * nor[0] + (fy + 1) * nor[1] + (fz    ) * nor[2] + d < 0.0f) return 1;
246                 if ((fx + 1) * nor[0] + (fy + 1) * nor[1] + (fz    ) * nor[2] + d < 0.0f) return 1;
247
248                 if ((fx    ) * nor[0] + (fy    ) * nor[1] + (fz + 1) * nor[2] + d < 0.0f) return 1;
249                 if ((fx + 1) * nor[0] + (fy    ) * nor[1] + (fz + 1) * nor[2] + d < 0.0f) return 1;
250                 if ((fx    ) * nor[0] + (fy + 1) * nor[1] + (fz + 1) * nor[2] + d < 0.0f) return 1;
251                 if ((fx + 1) * nor[0] + (fy + 1) * nor[1] + (fz + 1) * nor[2] + d < 0.0f) return 1;
252         }
253         else {
254                 if ((fx + 1) * nor[0] + (fy    ) * nor[1] + (fz    ) * nor[2] + d > 0.0f) return 1;
255                 if ((fx    ) * nor[0] + (fy + 1) * nor[1] + (fz    ) * nor[2] + d > 0.0f) return 1;
256                 if ((fx + 1) * nor[0] + (fy + 1) * nor[1] + (fz    ) * nor[2] + d > 0.0f) return 1;
257
258                 if ((fx    ) * nor[0] + (fy    ) * nor[1] + (fz + 1) * nor[2] + d > 0.0f) return 1;
259                 if ((fx + 1) * nor[0] + (fy    ) * nor[1] + (fz + 1) * nor[2] + d > 0.0f) return 1;
260                 if ((fx    ) * nor[0] + (fy + 1) * nor[1] + (fz + 1) * nor[2] + d > 0.0f) return 1;
261                 if ((fx + 1) * nor[0] + (fy + 1) * nor[1] + (fz + 1) * nor[2] + d > 0.0f) return 1;
262         }
263
264         return 0;
265 }
266
267 static void ocwrite(Octree *oc, RayFace *face, int quad, short x, short y, short z, float rtf[4][3])
268 {
269         Branch *br;
270         Node *no;
271         short a, oc0, oc1, oc2, oc3, oc4, oc5;
272
273         x <<= 2;
274         y <<= 1;
275
276         br = oc->adrbranch[0];
277
278         if (oc->ocres == 512) {
279                 oc0 = ((x & 1024) + (y & 512) + (z & 256)) >> 8;
280                 br = addbranch(oc, br, oc0);
281         }
282         if (oc->ocres >= 256) {
283                 oc0 = ((x & 512) + (y & 256) + (z & 128)) >> 7;
284                 br = addbranch(oc, br, oc0);
285         }
286         if (oc->ocres >= 128) {
287                 oc0 = ((x & 256) + (y & 128) + (z & 64)) >> 6;
288                 br = addbranch(oc, br, oc0);
289         }
290
291         oc0 = ((x & 128) + (y & 64) + (z & 32)) >> 5;
292         oc1 = ((x & 64) + (y & 32) + (z & 16)) >> 4;
293         oc2 = ((x & 32) + (y & 16) + (z & 8)) >> 3;
294         oc3 = ((x & 16) + (y & 8) + (z & 4)) >> 2;
295         oc4 = ((x & 8) + (y & 4) + (z & 2)) >> 1;
296         oc5 = ((x & 4) + (y & 2) + (z & 1));
297
298         br = addbranch(oc, br, oc0);
299         br = addbranch(oc, br, oc1);
300         br = addbranch(oc, br, oc2);
301         br = addbranch(oc, br, oc3);
302         br = addbranch(oc, br, oc4);
303         no = (Node *)br->b[oc5];
304         if (no == NULL) br->b[oc5] = (Branch *)(no = addnode(oc));
305
306         while (no->next) no = no->next;
307
308         a = 0;
309         if (no->v[7]) {     /* node full */
310                 no->next = addnode(oc);
311                 no = no->next;
312         }
313         else {
314                 while (no->v[a] != NULL) a++;
315         }
316
317         no->v[a] = (RayFace *) RE_rayobject_align(face);
318
319         if (quad)
320                 calc_ocval_face(rtf[0], rtf[1], rtf[2], rtf[3], x >> 2, y >> 1, z, &no->ov[a]);
321         else
322                 calc_ocval_face(rtf[0], rtf[1], rtf[2], NULL, x >> 2, y >> 1, z, &no->ov[a]);
323 }
324
325 static void d2dda(Octree *oc, short b1, short b2, short c1, short c2, char *ocface, short rts[4][3], float rtf[4][3])
326 {
327         int ocx1, ocx2, ocy1, ocy2;
328         int x, y, dx = 0, dy = 0;
329         float ox1, ox2, oy1, oy2;
330         float lambda, lambda_o, lambda_x, lambda_y, ldx, ldy;
331
332         ocx1 = rts[b1][c1];
333         ocy1 = rts[b1][c2];
334         ocx2 = rts[b2][c1];
335         ocy2 = rts[b2][c2];
336
337         if (ocx1 == ocx2 && ocy1 == ocy2) {
338                 ocface[oc->ocres * ocx1 + ocy1] = 1;
339                 return;
340         }
341
342         ox1 = rtf[b1][c1];
343         oy1 = rtf[b1][c2];
344         ox2 = rtf[b2][c1];
345         oy2 = rtf[b2][c2];
346
347         if (ox1 != ox2) {
348                 if (ox2 - ox1 > 0.0f) {
349                         lambda_x = (ox1 - ocx1 - 1.0f) / (ox1 - ox2);
350                         ldx = -1.0f / (ox1 - ox2);
351                         dx = 1;
352                 }
353                 else {
354                         lambda_x = (ox1 - ocx1) / (ox1 - ox2);
355                         ldx = 1.0f / (ox1 - ox2);
356                         dx = -1;
357                 }
358         }
359         else {
360                 lambda_x = 1.0f;
361                 ldx = 0;
362         }
363
364         if (oy1 != oy2) {
365                 if (oy2 - oy1 > 0.0f) {
366                         lambda_y = (oy1 - ocy1 - 1.0f) / (oy1 - oy2);
367                         ldy = -1.0f / (oy1 - oy2);
368                         dy = 1;
369                 }
370                 else {
371                         lambda_y = (oy1 - ocy1) / (oy1 - oy2);
372                         ldy = 1.0f / (oy1 - oy2);
373                         dy = -1;
374                 }
375         }
376         else {
377                 lambda_y = 1.0f;
378                 ldy = 0;
379         }
380
381         x = ocx1; y = ocy1;
382         lambda = MIN2(lambda_x, lambda_y);
383
384         while (true) {
385
386                 if (x < 0 || y < 0 || x >= oc->ocres || y >= oc->ocres) {
387                         /* pass*/
388                 }
389                 else {
390                         ocface[oc->ocres * x + y] = 1;
391                 }
392
393                 lambda_o = lambda;
394                 if (lambda_x == lambda_y) {
395                         lambda_x += ldx;
396                         x += dx;
397                         lambda_y += ldy;
398                         y += dy;
399                 }
400                 else {
401                         if (lambda_x < lambda_y) {
402                                 lambda_x += ldx;
403                                 x += dx;
404                         }
405                         else {
406                                 lambda_y += ldy;
407                                 y += dy;
408                         }
409                 }
410                 lambda = MIN2(lambda_x, lambda_y);
411                 if (lambda == lambda_o) break;
412                 if (lambda >= 1.0f) break;
413         }
414         ocface[oc->ocres * ocx2 + ocy2] = 1;
415 }
416
417 static void filltriangle(Octree *oc, short c1, short c2, char *ocface, short *ocmin, short *ocmax)
418 {
419         int a, x, y, y1, y2;
420
421         for (x = ocmin[c1]; x <= ocmax[c1]; x++) {
422                 a = oc->ocres * x;
423                 for (y = ocmin[c2]; y <= ocmax[c2]; y++) {
424                         if (ocface[a + y]) {
425                                 y++;
426                                 while (ocface[a + y] && y != ocmax[c2]) y++;
427                                 for (y1 = ocmax[c2]; y1 > y; y1--) {
428                                         if (ocface[a + y1]) {
429                                                 for (y2 = y; y2 <= y1; y2++) ocface[a + y2] = 1;
430                                                 y1 = 0;
431                                         }
432                                 }
433                                 y = ocmax[c2];
434                         }
435                 }
436         }
437 }
438
439 static void RE_rayobject_octree_free(RayObject *tree)
440 {
441         Octree *oc = (Octree *)tree;
442
443 #if 0
444         printf("branches %d nodes %d\n", oc->branchcount, oc->nodecount);
445         printf("raycount %d\n", raycount);
446         printf("ray coherent %d\n", coherent_ray);
447         printf("accepted %d rejected %d\n", accepted, rejected);
448 #endif
449         if (oc->ocface)
450                 MEM_freeN(oc->ocface);
451
452         if (oc->adrbranch) {
453                 int a = 0;
454                 while (oc->adrbranch[a]) {
455                         MEM_freeN(oc->adrbranch[a]);
456                         oc->adrbranch[a] = NULL;
457                         a++;
458                 }
459                 MEM_freeN(oc->adrbranch);
460                 oc->adrbranch = NULL;
461         }
462         oc->branchcount = 0;
463
464         if (oc->adrnode) {
465                 int a = 0;
466                 while (oc->adrnode[a]) {
467                         MEM_freeN(oc->adrnode[a]);
468                         oc->adrnode[a] = NULL;
469                         a++;
470                 }
471                 MEM_freeN(oc->adrnode);
472                 oc->adrnode = NULL;
473         }
474         oc->nodecount = 0;
475
476         MEM_freeN(oc);
477 }
478
479
480 RayObject *RE_rayobject_octree_create(int ocres, int size)
481 {
482         Octree *oc = (Octree *)MEM_callocN(sizeof(Octree), "Octree");
483         assert(RE_rayobject_isAligned(oc) );  /* RayObject API assumes real data to be 4-byte aligned */
484
485         oc->rayobj.api = &octree_api;
486
487         oc->ocres = ocres;
488
489         oc->ro_nodes = (RayFace **)MEM_callocN(sizeof(RayFace *) * size, "octree rayobject nodes");
490         oc->ro_nodes_size = size;
491         oc->ro_nodes_used = 0;
492
493
494         return RE_rayobject_unalignRayAPI((RayObject *) oc);
495 }
496
497
498 static void RE_rayobject_octree_add(RayObject *tree, RayObject *node)
499 {
500         Octree *oc = (Octree *)tree;
501
502         assert(RE_rayobject_isRayFace(node) );
503         assert(oc->ro_nodes_used < oc->ro_nodes_size);
504         oc->ro_nodes[oc->ro_nodes_used++] = (RayFace *)RE_rayobject_align(node);
505 }
506
507 static void octree_fill_rayface(Octree *oc, RayFace *face)
508 {
509         float ocfac[3], rtf[4][3];
510         float co1[3], co2[3], co3[3], co4[3];
511         short rts[4][3];
512         short ocmin[3], ocmax[3];
513         char *ocface = oc->ocface;   // front, top, size view of face, to fill in
514         int a, b, c, oc1, oc2, oc3, oc4, x, y, z, ocres2;
515
516         ocfac[0] = oc->ocfacx;
517         ocfac[1] = oc->ocfacy;
518         ocfac[2] = oc->ocfacz;
519
520         ocres2 = oc->ocres * oc->ocres;
521
522         copy_v3_v3(co1, face->v1);
523         copy_v3_v3(co2, face->v2);
524         copy_v3_v3(co3, face->v3);
525         if (RE_rayface_isQuad(face))
526                 copy_v3_v3(co4, face->v4);
527
528         for (c = 0; c < 3; c++) {
529                 rtf[0][c] = (co1[c] - oc->min[c]) * ocfac[c];
530                 rts[0][c] = (short)rtf[0][c];
531                 rtf[1][c] = (co2[c] - oc->min[c]) * ocfac[c];
532                 rts[1][c] = (short)rtf[1][c];
533                 rtf[2][c] = (co3[c] - oc->min[c]) * ocfac[c];
534                 rts[2][c] = (short)rtf[2][c];
535                 if (RE_rayface_isQuad(face)) {
536                         rtf[3][c] = (co4[c] - oc->min[c]) * ocfac[c];
537                         rts[3][c] = (short)rtf[3][c];
538                 }
539         }
540
541         for (c = 0; c < 3; c++) {
542                 oc1 = rts[0][c];
543                 oc2 = rts[1][c];
544                 oc3 = rts[2][c];
545                 if (!RE_rayface_isQuad(face)) {
546                         ocmin[c] = min_iii(oc1, oc2, oc3);
547                         ocmax[c] = max_iii(oc1, oc2, oc3);
548                 }
549                 else {
550                         oc4 = rts[3][c];
551                         ocmin[c] = min_iiii(oc1, oc2, oc3, oc4);
552                         ocmax[c] = max_iiii(oc1, oc2, oc3, oc4);
553                 }
554                 if (ocmax[c] > oc->ocres - 1) ocmax[c] = oc->ocres - 1;
555                 if (ocmin[c] < 0) ocmin[c] = 0;
556         }
557
558         if (ocmin[0] == ocmax[0] && ocmin[1] == ocmax[1] && ocmin[2] == ocmax[2]) {
559                 ocwrite(oc, face, RE_rayface_isQuad(face), ocmin[0], ocmin[1], ocmin[2], rtf);
560         }
561         else {
562
563                 d2dda(oc, 0, 1, 0, 1, ocface + ocres2, rts, rtf);
564                 d2dda(oc, 0, 1, 0, 2, ocface, rts, rtf);
565                 d2dda(oc, 0, 1, 1, 2, ocface + 2 * ocres2, rts, rtf);
566                 d2dda(oc, 1, 2, 0, 1, ocface + ocres2, rts, rtf);
567                 d2dda(oc, 1, 2, 0, 2, ocface, rts, rtf);
568                 d2dda(oc, 1, 2, 1, 2, ocface + 2 * ocres2, rts, rtf);
569                 if (!RE_rayface_isQuad(face)) {
570                         d2dda(oc, 2, 0, 0, 1, ocface + ocres2, rts, rtf);
571                         d2dda(oc, 2, 0, 0, 2, ocface, rts, rtf);
572                         d2dda(oc, 2, 0, 1, 2, ocface + 2 * ocres2, rts, rtf);
573                 }
574                 else {
575                         d2dda(oc, 2, 3, 0, 1, ocface + ocres2, rts, rtf);
576                         d2dda(oc, 2, 3, 0, 2, ocface, rts, rtf);
577                         d2dda(oc, 2, 3, 1, 2, ocface + 2 * ocres2, rts, rtf);
578                         d2dda(oc, 3, 0, 0, 1, ocface + ocres2, rts, rtf);
579                         d2dda(oc, 3, 0, 0, 2, ocface, rts, rtf);
580                         d2dda(oc, 3, 0, 1, 2, ocface + 2 * ocres2, rts, rtf);
581                 }
582                 /* nothing todo with triangle..., just fills :) */
583                 filltriangle(oc, 0, 1, ocface + ocres2, ocmin, ocmax);
584                 filltriangle(oc, 0, 2, ocface, ocmin, ocmax);
585                 filltriangle(oc, 1, 2, ocface + 2 * ocres2, ocmin, ocmax);
586
587                 /* init static vars here */
588                 face_in_node(face, 0, 0, 0, rtf);
589
590                 for (x = ocmin[0]; x <= ocmax[0]; x++) {
591                         a = oc->ocres * x;
592                         for (y = ocmin[1]; y <= ocmax[1]; y++) {
593                                 if (ocface[a + y + ocres2]) {
594                                         b = oc->ocres * y + 2 * ocres2;
595                                         for (z = ocmin[2]; z <= ocmax[2]; z++) {
596                                                 if (ocface[b + z] && ocface[a + z]) {
597                                                         if (face_in_node(NULL, x, y, z, rtf))
598                                                                 ocwrite(oc, face, RE_rayface_isQuad(face), x, y, z, rtf);
599                                                 }
600                                         }
601                                 }
602                         }
603                 }
604
605                 /* same loops to clear octree, doubt it can be done smarter */
606                 for (x = ocmin[0]; x <= ocmax[0]; x++) {
607                         a = oc->ocres * x;
608                         for (y = ocmin[1]; y <= ocmax[1]; y++) {
609                                 /* x-y */
610                                 ocface[a + y + ocres2] = 0;
611
612                                 b = oc->ocres * y + 2 * ocres2;
613                                 for (z = ocmin[2]; z <= ocmax[2]; z++) {
614                                         /* y-z */
615                                         ocface[b + z] = 0;
616                                         /* x-z */
617                                         ocface[a + z] = 0;
618                                 }
619                         }
620                 }
621         }
622 }
623
624 static void RE_rayobject_octree_done(RayObject *tree)
625 {
626         Octree *oc = (Octree *)tree;
627         int c;
628         float t00, t01, t02;
629         int ocres2 = oc->ocres * oc->ocres;
630
631         INIT_MINMAX(oc->min, oc->max);
632
633         /* Calculate Bounding Box */
634         for (c = 0; c < oc->ro_nodes_used; c++)
635                 RE_rayobject_merge_bb(RE_rayobject_unalignRayFace(oc->ro_nodes[c]), oc->min, oc->max);
636
637         /* Alloc memory */
638         oc->adrbranch = (Branch **)MEM_callocN(sizeof(void *) * BRANCH_ARRAY, "octree branches");
639         oc->adrnode = (Node **)MEM_callocN(sizeof(void *) * NODE_ARRAY, "octree nodes");
640
641         oc->adrbranch[0] = (Branch *)MEM_callocN(4096 * sizeof(Branch), "makeoctree");
642
643         /* the lookup table, per face, for which nodes to fill in */
644         oc->ocface = (char *)MEM_callocN(3 * ocres2 + 8, "ocface");
645         memset(oc->ocface, 0, 3 * ocres2);
646
647         for (c = 0; c < 3; c++) { /* octree enlarge, still needed? */
648                 oc->min[c] -= 0.01f;
649                 oc->max[c] += 0.01f;
650         }
651
652         t00 = oc->max[0] - oc->min[0];
653         t01 = oc->max[1] - oc->min[1];
654         t02 = oc->max[2] - oc->min[2];
655
656         /* this minus 0.1 is old safety... seems to be needed? */
657         oc->ocfacx = (oc->ocres - 0.1f) / t00;
658         oc->ocfacy = (oc->ocres - 0.1f) / t01;
659         oc->ocfacz = (oc->ocres - 0.1f) / t02;
660
661         oc->ocsize = sqrtf(t00 * t00 + t01 * t01 + t02 * t02);  /* global, max size octree */
662
663         for (c = 0; c < oc->ro_nodes_used; c++) {
664                 octree_fill_rayface(oc, oc->ro_nodes[c]);
665         }
666
667         MEM_freeN(oc->ocface);
668         oc->ocface = NULL;
669         MEM_freeN(oc->ro_nodes);
670         oc->ro_nodes = NULL;
671
672 #if 0
673         printf("%f %f - %f\n", oc->min[0], oc->max[0], oc->ocfacx);
674         printf("%f %f - %f\n", oc->min[1], oc->max[1], oc->ocfacy);
675         printf("%f %f - %f\n", oc->min[2], oc->max[2], oc->ocfacz);
676 #endif
677 }
678
679 static void RE_rayobject_octree_bb(RayObject *tree, float *min, float *max)
680 {
681         Octree *oc = (Octree *)tree;
682         DO_MINMAX(oc->min, min, max);
683         DO_MINMAX(oc->max, min, max);
684 }
685
686 /* check all faces in this node */
687 static int testnode(Octree *UNUSED(oc), Isect *is, Node *no, OcVal ocval)
688 {
689         short nr = 0;
690
691         /* return on any first hit */
692         if (is->mode == RE_RAY_SHADOW) {
693
694                 for (; no; no = no->next) {
695                         for (nr = 0; nr < 8; nr++) {
696                                 RayFace *face = no->v[nr];
697                                 OcVal     *ov = no->ov + nr;
698
699                                 if (!face) break;
700
701                                 if ( (ov->ocx & ocval.ocx) && (ov->ocy & ocval.ocy) && (ov->ocz & ocval.ocz) ) {
702                                         if (RE_rayobject_intersect(RE_rayobject_unalignRayFace(face), is) )
703                                                 return 1;
704                                 }
705                         }
706                 }
707         }
708         else {
709                 /* else mirror or glass or shadowtra, return closest face  */
710                 int found = 0;
711
712                 for (; no; no = no->next) {
713                         for (nr = 0; nr < 8; nr++) {
714                                 RayFace *face = no->v[nr];
715                                 OcVal     *ov = no->ov + nr;
716
717                                 if (!face) break;
718
719                                 if ( (ov->ocx & ocval.ocx) && (ov->ocy & ocval.ocy) && (ov->ocz & ocval.ocz) ) {
720                                         if (RE_rayobject_intersect(RE_rayobject_unalignRayFace(face), is) ) {
721                                                 found = 1;
722                                         }
723                                 }
724                         }
725                 }
726
727                 return found;
728         }
729
730         return 0;
731 }
732
733 /* find the Node for the octree coord x y z */
734 static Node *ocread(Octree *oc, int x, int y, int z)
735 {
736         Branch *br;
737         int oc1;
738
739         x <<= 2;
740         y <<= 1;
741
742         br = oc->adrbranch[0];
743
744         if (oc->ocres == 512) {
745                 oc1 = ((x & 1024) + (y & 512) + (z & 256)) >> 8;
746                 br = br->b[oc1];
747                 if (br == NULL) {
748                         return NULL;
749                 }
750         }
751         if (oc->ocres >= 256) {
752                 oc1 = ((x & 512) + (y & 256) + (z & 128)) >> 7;
753                 br = br->b[oc1];
754                 if (br == NULL) {
755                         return NULL;
756                 }
757         }
758         if (oc->ocres >= 128) {
759                 oc1 = ((x & 256) + (y & 128) + (z & 64)) >> 6;
760                 br = br->b[oc1];
761                 if (br == NULL) {
762                         return NULL;
763                 }
764         }
765
766         oc1 = ((x & 128) + (y & 64) + (z & 32)) >> 5;
767         br = br->b[oc1];
768         if (br) {
769                 oc1 = ((x & 64) + (y & 32) + (z & 16)) >> 4;
770                 br = br->b[oc1];
771                 if (br) {
772                         oc1 = ((x & 32) + (y & 16) + (z & 8)) >> 3;
773                         br = br->b[oc1];
774                         if (br) {
775                                 oc1 = ((x & 16) + (y & 8) + (z & 4)) >> 2;
776                                 br = br->b[oc1];
777                                 if (br) {
778                                         oc1 = ((x & 8) + (y & 4) + (z & 2)) >> 1;
779                                         br = br->b[oc1];
780                                         if (br) {
781                                                 oc1 = ((x & 4) + (y & 2) + (z & 1));
782                                                 return (Node *)br->b[oc1];
783                                         }
784                                 }
785                         }
786                 }
787         }
788
789         return NULL;
790 }
791
792 static int cliptest(float p, float q, float *u1, float *u2)
793 {
794         float r;
795
796         if (p < 0.0f) {
797                 if (q < p) return 0;
798                 else if (q < 0.0f) {
799                         r = q / p;
800                         if (r > *u2) return 0;
801                         else if (r > *u1) *u1 = r;
802                 }
803         }
804         else {
805                 if (p > 0.0f) {
806                         if (q < 0.0f) return 0;
807                         else if (q < p) {
808                                 r = q / p;
809                                 if (r < *u1) return 0;
810                                 else if (r < *u2) *u2 = r;
811                         }
812                 }
813                 else if (q < 0.0f) return 0;
814         }
815         return 1;
816 }
817
818 /* extensive coherence checks/storage cancels out the benefit of it, and gives errors... we
819  * need better methods, sample code commented out below (ton) */
820
821 #if 0
822
823 in top : static int coh_nodes[16 * 16 * 16][6];
824 in makeoctree : memset(coh_nodes, 0, sizeof(coh_nodes));
825
826 static void add_coherence_test(int ocx1, int ocx2, int ocy1, int ocy2, int ocz1, int ocz2)
827 {
828         short *sp;
829
830         sp = coh_nodes[(ocx2 & 15) + 16 * (ocy2 & 15) + 256 * (ocz2 & 15)];
831         sp[0] = ocx1; sp[1] = ocy1; sp[2] = ocz1;
832         sp[3] = ocx2; sp[4] = ocy2; sp[5] = ocz2;
833
834 }
835
836 static int do_coherence_test(int ocx1, int ocx2, int ocy1, int ocy2, int ocz1, int ocz2)
837 {
838         short *sp;
839
840         sp = coh_nodes[(ocx2 & 15) + 16 * (ocy2 & 15) + 256 * (ocz2 & 15)];
841         if (sp[0] == ocx1 && sp[1] == ocy1 && sp[2] == ocz1 &&
842             sp[3] == ocx2 && sp[4] == ocy2 && sp[5] == ocz2) return 1;
843         return 0;
844 }
845
846 #endif
847
848 /* return 1: found valid intersection */
849 /* starts with is->orig.face */
850 static int RE_rayobject_octree_intersect(RayObject *tree, Isect *is)
851 {
852         Octree *oc = (Octree *)tree;
853         Node *no;
854         OcVal ocval;
855         float vec1[3], vec2[3], start[3], end[3];
856         float u1, u2, ox1, ox2, oy1, oy2, oz1, oz2;
857         float lambda_o, lambda_x, ldx, lambda_y, ldy, lambda_z, ldz, dda_lambda;
858         float o_lambda = 0;
859         int dx, dy, dz;
860         int xo, yo, zo, c1 = 0;
861         int ocx1, ocx2, ocy1, ocy2, ocz1, ocz2;
862
863         /* clip with octree */
864         if (oc->branchcount == 0) return 0;
865
866         /* do this before intersect calls */
867 #if 0
868         is->facecontr = NULL;                /* to check shared edge */
869         is->obcontr = 0;
870         is->faceisect = is->isect = 0;        /* shared edge, quad half flag */
871         is->userdata = oc->userdata;
872 #endif
873
874         copy_v3_v3(start, is->start);
875         madd_v3_v3v3fl(end, is->start, is->dir, is->dist);
876         ldx = is->dir[0] * is->dist;
877         o_lambda = is->dist;
878         u1 = 0.0f;
879         u2 = 1.0f;
880
881         /* clip with octree cube */
882         if (cliptest(-ldx, start[0] - oc->min[0], &u1, &u2)) {
883                 if (cliptest(ldx, oc->max[0] - start[0], &u1, &u2)) {
884                         ldy = is->dir[1] * is->dist;
885                         if (cliptest(-ldy, start[1] - oc->min[1], &u1, &u2)) {
886                                 if (cliptest(ldy, oc->max[1] - start[1], &u1, &u2)) {
887                                         ldz = is->dir[2] * is->dist;
888                                         if (cliptest(-ldz, start[2] - oc->min[2], &u1, &u2)) {
889                                                 if (cliptest(ldz, oc->max[2] - start[2], &u1, &u2)) {
890                                                         c1 = 1;
891                                                         if (u2 < 1.0f) {
892                                                                 end[0] = start[0] + u2 * ldx;
893                                                                 end[1] = start[1] + u2 * ldy;
894                                                                 end[2] = start[2] + u2 * ldz;
895                                                         }
896
897                                                         if (u1 > 0.0f) {
898                                                                 start[0] += u1 * ldx;
899                                                                 start[1] += u1 * ldy;
900                                                                 start[2] += u1 * ldz;
901                                                         }
902                                                 }
903                                         }
904                                 }
905                         }
906                 }
907         }
908
909         if (c1 == 0) return 0;
910
911         /* reset static variables in ocread */
912         //ocread(oc, oc->ocres, 0, 0);
913
914         /* setup 3dda to traverse octree */
915         ox1 = (start[0] - oc->min[0]) * oc->ocfacx;
916         oy1 = (start[1] - oc->min[1]) * oc->ocfacy;
917         oz1 = (start[2] - oc->min[2]) * oc->ocfacz;
918         ox2 = (end[0] - oc->min[0]) * oc->ocfacx;
919         oy2 = (end[1] - oc->min[1]) * oc->ocfacy;
920         oz2 = (end[2] - oc->min[2]) * oc->ocfacz;
921
922         ocx1 = (int)ox1;
923         ocy1 = (int)oy1;
924         ocz1 = (int)oz1;
925         ocx2 = (int)ox2;
926         ocy2 = (int)oy2;
927         ocz2 = (int)oz2;
928
929         if (ocx1 == ocx2 && ocy1 == ocy2 && ocz1 == ocz2) {
930                 no = ocread(oc, ocx1, ocy1, ocz1);
931                 if (no) {
932                         /* exact intersection with node */
933                         vec1[0] = ox1; vec1[1] = oy1; vec1[2] = oz1;
934                         vec2[0] = ox2; vec2[1] = oy2; vec2[2] = oz2;
935                         calc_ocval_ray(&ocval, (float)ocx1, (float)ocy1, (float)ocz1, vec1, vec2);
936                         if (testnode(oc, is, no, ocval) ) return 1;
937                 }
938         }
939         else {
940                 int found = 0;
941                 //static int coh_ocx1, coh_ocx2, coh_ocy1, coh_ocy2, coh_ocz1, coh_ocz2;
942                 float dox, doy, doz;
943                 int eqval;
944
945                 /* calc lambda en ld */
946                 dox = ox1 - ox2;
947                 doy = oy1 - oy2;
948                 doz = oz1 - oz2;
949
950                 if (dox < -FLT_EPSILON) {
951                         ldx = -1.0f / dox;
952                         lambda_x = (ocx1 - ox1 + 1.0f) * ldx;
953                         dx = 1;
954                 }
955                 else if (dox > FLT_EPSILON) {
956                         ldx = 1.0f / dox;
957                         lambda_x = (ox1 - ocx1) * ldx;
958                         dx = -1;
959                 }
960                 else {
961                         lambda_x = 1.0f;
962                         ldx = 0;
963                         dx = 0;
964                 }
965
966                 if (doy < -FLT_EPSILON) {
967                         ldy = -1.0f / doy;
968                         lambda_y = (ocy1 - oy1 + 1.0f) * ldy;
969                         dy = 1;
970                 }
971                 else if (doy > FLT_EPSILON) {
972                         ldy = 1.0f / doy;
973                         lambda_y = (oy1 - ocy1) * ldy;
974                         dy = -1;
975                 }
976                 else {
977                         lambda_y = 1.0f;
978                         ldy = 0;
979                         dy = 0;
980                 }
981
982                 if (doz < -FLT_EPSILON) {
983                         ldz = -1.0f / doz;
984                         lambda_z = (ocz1 - oz1 + 1.0f) * ldz;
985                         dz = 1;
986                 }
987                 else if (doz > FLT_EPSILON) {
988                         ldz = 1.0f / doz;
989                         lambda_z = (oz1 - ocz1) * ldz;
990                         dz = -1;
991                 }
992                 else {
993                         lambda_z = 1.0f;
994                         ldz = 0;
995                         dz = 0;
996                 }
997
998                 xo = ocx1; yo = ocy1; zo = ocz1;
999                 dda_lambda = min_fff(lambda_x, lambda_y, lambda_z);
1000
1001                 vec2[0] = ox1;
1002                 vec2[1] = oy1;
1003                 vec2[2] = oz1;
1004
1005                 /* this loop has been constructed to make sure the first and last node of ray
1006                  * are always included, even when dda_lambda==1.0f or larger */
1007
1008                 while (true) {
1009
1010                         no = ocread(oc, xo, yo, zo);
1011                         if (no) {
1012
1013                                 /* calculate ray intersection with octree node */
1014                                 copy_v3_v3(vec1, vec2);
1015                                 // dox, y, z is negative
1016                                 vec2[0] = ox1 - dda_lambda * dox;
1017                                 vec2[1] = oy1 - dda_lambda * doy;
1018                                 vec2[2] = oz1 - dda_lambda * doz;
1019                                 calc_ocval_ray(&ocval, (float)xo, (float)yo, (float)zo, vec1, vec2);
1020
1021                                 //is->dist = (u1 + dda_lambda * (u2 - u1)) * o_lambda;
1022                                 if (testnode(oc, is, no, ocval) )
1023                                         found = 1;
1024
1025                                 if (is->dist < (u1 + dda_lambda * (u2 - u1)) * o_lambda)
1026                                         return found;
1027                         }
1028
1029
1030                         lambda_o = dda_lambda;
1031
1032                         /* traversing octree nodes need careful detection of smallest values, with proper
1033                          * exceptions for equal lambdas */
1034                         eqval = (lambda_x == lambda_y);
1035                         if (lambda_y == lambda_z) eqval += 2;
1036                         if (lambda_x == lambda_z) eqval += 4;
1037
1038                         if (eqval) {    // only 4 cases exist!
1039                                 if (eqval == 7) { // x=y=z
1040                                         xo += dx; lambda_x += ldx;
1041                                         yo += dy; lambda_y += ldy;
1042                                         zo += dz; lambda_z += ldz;
1043                                 }
1044                                 else if (eqval == 1) { // x=y
1045                                         if (lambda_y < lambda_z) {
1046                                                 xo += dx; lambda_x += ldx;
1047                                                 yo += dy; lambda_y += ldy;
1048                                         }
1049                                         else {
1050                                                 zo += dz; lambda_z += ldz;
1051                                         }
1052                                 }
1053                                 else if (eqval == 2) { // y=z
1054                                         if (lambda_x < lambda_y) {
1055                                                 xo += dx; lambda_x += ldx;
1056                                         }
1057                                         else {
1058                                                 yo += dy; lambda_y += ldy;
1059                                                 zo += dz; lambda_z += ldz;
1060                                         }
1061                                 }
1062                                 else { // x=z
1063                                         if (lambda_y < lambda_x) {
1064                                                 yo += dy; lambda_y += ldy;
1065                                         }
1066                                         else {
1067                                                 xo += dx; lambda_x += ldx;
1068                                                 zo += dz; lambda_z += ldz;
1069                                         }
1070                                 }
1071                         }
1072                         else {  // all three different, just three cases exist
1073                                 eqval = (lambda_x < lambda_y);
1074                                 if (lambda_y < lambda_z) eqval += 2;
1075                                 if (lambda_x < lambda_z) eqval += 4;
1076
1077                                 if (eqval == 7 || eqval == 5) { // x smallest
1078                                         xo += dx; lambda_x += ldx;
1079                                 }
1080                                 else if (eqval == 2 || eqval == 6) { // y smallest
1081                                         yo += dy; lambda_y += ldy;
1082                                 }
1083                                 else { // z smallest
1084                                         zo += dz; lambda_z += ldz;
1085                                 }
1086
1087                         }
1088
1089                         dda_lambda = min_fff(lambda_x, lambda_y, lambda_z);
1090                         if (dda_lambda == lambda_o) break;
1091                         /* to make sure the last node is always checked */
1092                         if (lambda_o >= 1.0f) break;
1093                 }
1094         }
1095
1096         /* reached end, no intersections found */
1097         return 0;
1098 }
1099
1100
1101