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