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