Added neighbour test on detected ray hit
[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 RayObject *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         RayObject **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]= (RayObject*)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 *RayObject_octree_create(int ocres, int size)
451 {
452         Octree *oc= MEM_callocN(sizeof(Octree), "Octree");
453         
454         oc->rayobj.api = &octree_api;
455         
456         oc->ocres = ocres;
457         
458         oc->ro_nodes = MEM_callocN(sizeof(RayObject*)*size, "octree rayobject nodes");
459         oc->ro_nodes_size = size;
460         oc->ro_nodes_used = 0;
461
462         
463         return RayObject_unalign((RayObject*) oc);
464 }
465
466
467 static void RayObject_octree_add(RayObject *tree, RayObject *node)
468 {
469         Octree *oc = (Octree*)tree;
470
471         assert( oc->ro_nodes_used < oc->ro_nodes_size );
472         oc->ro_nodes[ oc->ro_nodes_used++ ] = node;     
473 }
474
475 static void octree_fill_rayface(Octree *oc, RayFace *face)
476 {
477         float ocfac[3], rtf[4][3];
478         float co1[3], co2[3], co3[3], co4[3];
479         short rts[4][3];
480         short ocmin[3], ocmax[3];
481         char *ocface= oc->ocface;       // front, top, size view of face, to fill in
482         int a, b, c, oc1, oc2, oc3, oc4, x, y, z, ocres2;
483
484         ocfac[0]= oc->ocfacx;
485         ocfac[1]= oc->ocfacy;
486         ocfac[2]= oc->ocfacz;
487
488         ocres2= oc->ocres*oc->ocres;
489
490         VECCOPY(co1, face->v1);
491         VECCOPY(co2, face->v2);
492         VECCOPY(co3, face->v3);
493         if(face->v4)
494                 VECCOPY(co4, face->v4);
495
496         for(c=0;c<3;c++) {
497                 rtf[0][c]= (co1[c]-oc->min[c])*ocfac[c] ;
498                 rts[0][c]= (short)rtf[0][c];
499                 rtf[1][c]= (co2[c]-oc->min[c])*ocfac[c] ;
500                 rts[1][c]= (short)rtf[1][c];
501                 rtf[2][c]= (co3[c]-oc->min[c])*ocfac[c] ;
502                 rts[2][c]= (short)rtf[2][c];
503                 if(face->v4) {
504                         rtf[3][c]= (co4[c]-oc->min[c])*ocfac[c] ;
505                         rts[3][c]= (short)rtf[3][c];
506                 }
507         }
508         
509         for(c=0;c<3;c++) {
510                 oc1= rts[0][c];
511                 oc2= rts[1][c];
512                 oc3= rts[2][c];
513                 if(face->v4==NULL) {
514                         ocmin[c]= MIN3(oc1,oc2,oc3);
515                         ocmax[c]= MAX3(oc1,oc2,oc3);
516                 }
517                 else {
518                         oc4= rts[3][c];
519                         ocmin[c]= MIN4(oc1,oc2,oc3,oc4);
520                         ocmax[c]= MAX4(oc1,oc2,oc3,oc4);
521                 }
522                 if(ocmax[c]>oc->ocres-1) ocmax[c]=oc->ocres-1;
523                 if(ocmin[c]<0) ocmin[c]=0;
524         }
525         
526         if(ocmin[0]==ocmax[0] && ocmin[1]==ocmax[1] && ocmin[2]==ocmax[2]) {
527                 ocwrite(oc, face, (face->v4 != NULL), ocmin[0], ocmin[1], ocmin[2], rtf);
528         }
529         else {
530                 
531                 d2dda(oc, 0,1,0,1,ocface+ocres2,rts,rtf);
532                 d2dda(oc, 0,1,0,2,ocface,rts,rtf);
533                 d2dda(oc, 0,1,1,2,ocface+2*ocres2,rts,rtf);
534                 d2dda(oc, 1,2,0,1,ocface+ocres2,rts,rtf);
535                 d2dda(oc, 1,2,0,2,ocface,rts,rtf);
536                 d2dda(oc, 1,2,1,2,ocface+2*ocres2,rts,rtf);
537                 if(face->v4==NULL) {
538                         d2dda(oc, 2,0,0,1,ocface+ocres2,rts,rtf);
539                         d2dda(oc, 2,0,0,2,ocface,rts,rtf);
540                         d2dda(oc, 2,0,1,2,ocface+2*ocres2,rts,rtf);
541                 }
542                 else {
543                         d2dda(oc, 2,3,0,1,ocface+ocres2,rts,rtf);
544                         d2dda(oc, 2,3,0,2,ocface,rts,rtf);
545                         d2dda(oc, 2,3,1,2,ocface+2*ocres2,rts,rtf);
546                         d2dda(oc, 3,0,0,1,ocface+ocres2,rts,rtf);
547                         d2dda(oc, 3,0,0,2,ocface,rts,rtf);
548                         d2dda(oc, 3,0,1,2,ocface+2*ocres2,rts,rtf);
549                 }
550                 /* nothing todo with triangle..., just fills :) */
551                 filltriangle(oc, 0,1,ocface+ocres2,ocmin,ocmax);
552                 filltriangle(oc, 0,2,ocface,ocmin,ocmax);
553                 filltriangle(oc, 1,2,ocface+2*ocres2,ocmin,ocmax);
554                 
555                 /* init static vars here */
556                 face_in_node(face, 0,0,0, rtf);
557                 
558                 for(x=ocmin[0];x<=ocmax[0];x++) {
559                         a= oc->ocres*x;
560                         for(y=ocmin[1];y<=ocmax[1];y++) {
561                                 if(ocface[a+y+ocres2]) {
562                                         b= oc->ocres*y+2*ocres2;
563                                         for(z=ocmin[2];z<=ocmax[2];z++) {
564                                                 if(ocface[b+z] && ocface[a+z]) {
565                                                         if(face_in_node(NULL, x, y, z, rtf))
566                                                                 ocwrite(oc, face, (face->v4 != NULL), x,y,z, rtf);
567                                                 }
568                                         }
569                                 }
570                         }
571                 }
572                 
573                 /* same loops to clear octree, doubt it can be done smarter */
574                 for(x=ocmin[0];x<=ocmax[0];x++) {
575                         a= oc->ocres*x;
576                         for(y=ocmin[1];y<=ocmax[1];y++) {
577                                 /* x-y */
578                                 ocface[a+y+ocres2]= 0;
579
580                                 b= oc->ocres*y + 2*ocres2;
581                                 for(z=ocmin[2];z<=ocmax[2];z++) {
582                                         /* y-z */
583                                         ocface[b+z]= 0;
584                                         /* x-z */
585                                         ocface[a+z]= 0;
586                                 }
587                         }
588                 }
589         }
590 }
591
592 static void RayObject_octree_done(RayObject *tree)
593 {
594         Octree *oc = (Octree*)tree;
595         int c;
596         float t00, t01, t02;
597         int ocres2 = oc->ocres*oc->ocres;
598         
599         INIT_MINMAX(oc->min, oc->max);
600         
601         /* Calculate Bounding Box */
602         for(c=0; c<oc->ro_nodes_used; c++)
603                 RayObject_merge_bb(oc->ro_nodes[c], oc->min, oc->max);
604                 
605         /* Alloc memory */
606         oc->adrbranch= MEM_callocN(sizeof(void *)*BRANCH_ARRAY, "octree branches");
607         oc->adrnode= MEM_callocN(sizeof(void *)*NODE_ARRAY, "octree nodes");
608         
609         oc->adrbranch[0]=(Branch *)MEM_callocN(4096*sizeof(Branch), "makeoctree");
610         
611         /* the lookup table, per face, for which nodes to fill in */
612         oc->ocface= MEM_callocN( 3*ocres2 + 8, "ocface");
613         memset(oc->ocface, 0, 3*ocres2);
614
615         for(c=0;c<3;c++) {      /* octree enlarge, still needed? */
616                 oc->min[c]-= 0.01f;
617                 oc->max[c]+= 0.01f;
618         }
619
620         t00= oc->max[0]-oc->min[0];
621         t01= oc->max[1]-oc->min[1];
622         t02= oc->max[2]-oc->min[2];
623         
624         /* this minus 0.1 is old safety... seems to be needed? */
625         oc->ocfacx= (oc->ocres-0.1)/t00;
626         oc->ocfacy= (oc->ocres-0.1)/t01;
627         oc->ocfacz= (oc->ocres-0.1)/t02;
628         
629         oc->ocsize= sqrt(t00*t00+t01*t01+t02*t02);      /* global, max size octree */
630
631         for(c=0; c<oc->ro_nodes_used; c++)
632         {
633                 assert( RayObject_isFace(oc->ro_nodes[c]) );
634                 octree_fill_rayface(oc, (RayFace*)oc->ro_nodes[c]);
635         }
636
637         MEM_freeN(oc->ocface);
638         MEM_freeN(oc->ro_nodes);
639         
640         printf("%f %f - %f\n", oc->min[0], oc->max[0], oc->ocfacx );
641         printf("%f %f - %f\n", oc->min[1], oc->max[1], oc->ocfacy );
642         printf("%f %f - %f\n", oc->min[2], oc->max[2], oc->ocfacz );
643 }
644
645 static void RayObject_octree_bb(RayObject *tree, float *min, float *max)
646 {
647         Octree *oc = (Octree*)tree;
648         DO_MINMAX(oc->min, min, max);
649         DO_MINMAX(oc->max, min, max);
650 }
651
652 /* check all faces in this node */
653 static int testnode(Octree *oc, Isect *is, Node *no, OcVal ocval)
654 {
655         short nr=0;
656
657         /* return on any first hit */
658         if(is->mode==RE_RAY_SHADOW) {
659         
660                 for(; no; no = no->next)
661                 for(nr=0; nr<8; nr++)
662                 {
663                         RayObject *face = no->v[nr];
664                         OcVal           *ov = no->ov+nr;
665                         
666                         if(!face) break; //TODO? return 0;
667                         
668                         if( (ov->ocx & ocval.ocx) && (ov->ocy & ocval.ocy) && (ov->ocz & ocval.ocz) )
669                         {
670                                 if( RayObject_intersect(face,is) )
671                                         return 1;
672                         }
673                 }
674         }
675         else
676         {                       /* else mirror or glass or shadowtra, return closest face  */
677                 Isect isect;
678                 int found= 0;
679                 
680                 assert(0);
681                 
682                 is->labda= 1.0f;        /* needed? */
683                 isect= *is;                     /* copy for sorting */
684                 
685                 for(; no; no = no->next)
686                 for(nr=0; nr<8; nr++)
687                 {
688                         RayObject *face = no->v[nr];
689                         OcVal           *ov = no->ov+nr;
690                         
691                         if(!face) return 0;
692                         
693                         if( (ov->ocx & ocval.ocx) && (ov->ocy & ocval.ocy) && (ov->ocz & ocval.ocz) )
694                         { 
695
696                                 if( RayObject_raycast(face,is) )
697                                         if(isect.labda<is->labda) {
698                                                 *is= isect;
699                                                 found= 1;
700                                         }
701                         }
702                 }
703                 
704                 return found;
705         }
706
707         return 0;
708 }
709
710 /* find the Node for the octree coord x y z */
711 static Node *ocread(Octree *oc, int x, int y, int z)
712 {
713         Branch *br;
714         int oc1;
715         
716         x<<=2;
717         y<<=1;
718         
719         br= oc->adrbranch[0];
720         
721         if(oc->ocres==512) {
722                 oc1= ((x & 1024)+(y & 512)+(z & 256))>>8;
723                 br= br->b[oc1];
724                 if(br==NULL) {
725                         return NULL;
726                 }
727         }
728         if(oc->ocres>=256) {
729                 oc1= ((x & 512)+(y & 256)+(z & 128))>>7;
730                 br= br->b[oc1];
731                 if(br==NULL) {
732                         return NULL;
733                 }
734         }
735         if(oc->ocres>=128) {
736                 oc1= ((x & 256)+(y & 128)+(z & 64))>>6;
737                 br= br->b[oc1];
738                 if(br==NULL) {
739                         return NULL;
740                 }
741         }
742         
743         oc1= ((x & 128)+(y & 64)+(z & 32))>>5;
744         br= br->b[oc1];
745         if(br) {
746                 oc1= ((x & 64)+(y & 32)+(z & 16))>>4;
747                 br= br->b[oc1];
748                 if(br) {
749                         oc1= ((x & 32)+(y & 16)+(z & 8))>>3;
750                         br= br->b[oc1];
751                         if(br) {
752                                 oc1= ((x & 16)+(y & 8)+(z & 4))>>2;
753                                 br= br->b[oc1];
754                                 if(br) {
755                                         oc1= ((x & 8)+(y & 4)+(z & 2))>>1;
756                                         br= br->b[oc1];
757                                         if(br) {
758                                                 oc1= ((x & 4)+(y & 2)+(z & 1));
759                                                 return (Node *)br->b[oc1];
760                                         }
761                                 }
762                         }
763                 }
764         }
765         
766         return NULL;
767 }
768
769 static int cliptest(float p, float q, float *u1, float *u2)
770 {
771         float r;
772
773         if(p<0.0f) {
774                 if(q<p) return 0;
775                 else if(q<0.0f) {
776                         r= q/p;
777                         if(r>*u2) return 0;
778                         else if(r>*u1) *u1=r;
779                 }
780         }
781         else {
782                 if(p>0.0f) {
783                         if(q<0.0f) return 0;
784                         else if(q<p) {
785                                 r= q/p;
786                                 if(r<*u1) return 0;
787                                 else if(r<*u2) *u2=r;
788                         }
789                 }
790                 else if(q<0.0f) return 0;
791         }
792         return 1;
793 }
794
795 /* extensive coherence checks/storage cancels out the benefit of it, and gives errors... we
796    need better methods, sample code commented out below (ton) */
797  
798 /*
799
800 in top: static int coh_nodes[16*16*16][6];
801 in makeoctree: memset(coh_nodes, 0, sizeof(coh_nodes));
802  
803 static void add_coherence_test(int ocx1, int ocx2, int ocy1, int ocy2, int ocz1, int ocz2)
804 {
805         short *sp;
806         
807         sp= coh_nodes[ (ocx2 & 15) + 16*(ocy2 & 15) + 256*(ocz2 & 15) ];
808         sp[0]= ocx1; sp[1]= ocy1; sp[2]= ocz1;
809         sp[3]= ocx2; sp[4]= ocy2; sp[5]= ocz2;
810         
811 }
812
813 static int do_coherence_test(int ocx1, int ocx2, int ocy1, int ocy2, int ocz1, int ocz2)
814 {
815         short *sp;
816         
817         sp= coh_nodes[ (ocx2 & 15) + 16*(ocy2 & 15) + 256*(ocz2 & 15) ];
818         if(sp[0]==ocx1 && sp[1]==ocy1 && sp[2]==ocz1 &&
819            sp[3]==ocx2 && sp[4]==ocy2 && sp[5]==ocz2) return 1;
820         return 0;
821 }
822
823 */
824
825 /* return 1: found valid intersection */
826 /* starts with is->orig.face */
827 static int RayObject_octree_intersect(RayObject *tree, Isect *is)
828 {
829         Octree *oc= (Octree*)tree;
830         Node *no;
831         OcVal ocval;
832         float vec1[3], vec2[3], end[3];
833         float u1,u2,ox1,ox2,oy1,oy2,oz1,oz2;
834         float labdao,labdax,ldx,labday,ldy,labdaz,ldz, ddalabda;
835         int dx,dy,dz;   
836         int xo,yo,zo,c1=0;
837         int ocx1,ocx2,ocy1, ocy2,ocz1,ocz2;
838         
839         /* clip with octree */
840         if(oc->branchcount==0) return 0;
841         
842         /* do this before intersect calls */
843 #if 0
844         is->facecontr= NULL;                            /* to check shared edge */
845         is->obcontr= 0;
846         is->faceisect= is->isect= 0;            /* shared edge, quad half flag */
847         is->userdata= oc->userdata;
848 #endif
849
850 #if 0
851         /* only for shadow! */
852         if(is->mode==RE_RAY_SHADOW) {
853         
854                 /* TODO check with last intersected shadow face */
855                 if(is->last.face!=NULL && !(is->last.face==is->orig.face && is->last.ob==is->orig.ob)) {
856                         if(checkfunc(is, is->last.ob, is->last.face)) {
857                                 is->ob= is->last.ob;
858                                 is->face= is->last.face;
859                                 VECSUB(is->vec, is->end, is->start);
860                                 if(RE_ray_face_intersection(is, oc->transformfunc, oc->coordsfunc)) return 1;
861                         }
862                 }
863         }
864 #endif
865         
866         VECADDFAC( end, is->start, is->vec, is->labda );
867         ldx= end[0] - is->start[0];
868         u1= 0.0f;
869         u2= 1.0f;
870         
871         /* clip with octree cube */
872         if(cliptest(-ldx, is->start[0]-oc->min[0], &u1,&u2)) {
873                 if(cliptest(ldx, oc->max[0]-is->start[0], &u1,&u2)) {
874                         ldy= end[1] - is->start[1];
875                         if(cliptest(-ldy, is->start[1]-oc->min[1], &u1,&u2)) {
876                                 if(cliptest(ldy, oc->max[1]-is->start[1], &u1,&u2)) {
877                                         ldz = end[2] - is->start[2];
878                                         if(cliptest(-ldz, is->start[2]-oc->min[2], &u1,&u2)) {
879                                                 if(cliptest(ldz, oc->max[2]-is->start[2], &u1,&u2)) {
880                                                         c1=1;
881                                                         if(u2<1.0f) {
882                                                                 end[0]= is->start[0]+u2*ldx;
883                                                                 end[1]= is->start[1]+u2*ldy;
884                                                                 end[2]= is->start[2]+u2*ldz;
885                                                         }
886                                                         if(u1>0.0f) {
887                                                                 assert( 0 );
888                                                                 is->start[0]+=u1*ldx;
889                                                                 is->start[1]+=u1*ldy;
890                                                                 is->start[2]+=u1*ldz;
891                                                         }
892                                                 }
893                                         }
894                                 }
895                         }
896                 }
897         }
898
899         if(c1==0) return 0;
900
901         /* reset static variables in ocread */
902         //ocread(oc, oc->ocres, 0, 0);
903
904         /* setup 3dda to traverse octree */
905         ox1= (is->start[0]-oc->min[0])*oc->ocfacx;
906         oy1= (is->start[1]-oc->min[1])*oc->ocfacy;
907         oz1= (is->start[2]-oc->min[2])*oc->ocfacz;
908         ox2= (end[0]-oc->min[0])*oc->ocfacx;
909         oy2= (end[1]-oc->min[1])*oc->ocfacy;
910         oz2= (end[2]-oc->min[2])*oc->ocfacz;
911
912         ocx1= (int)ox1;
913         ocy1= (int)oy1;
914         ocz1= (int)oz1;
915         ocx2= (int)ox2;
916         ocy2= (int)oy2;
917         ocz2= (int)oz2;
918         
919 //      for(ocx1=0; ocx1<oc->ocres; ocx1++)
920 //      for(ocy1=0; ocy1<oc->ocres; ocy1++)
921 //      for(ocz1=0; ocz1<oc->ocres; ocz1++)
922 //      {
923 //              no= ocread(oc, ocx1, ocy1, ocz1);
924 //              if( testnode(oc, is, no, ocval) ) return 1;
925 //      }
926
927 //      return 0;
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 //                      is->ddalabda= 1.0f;
937                         if( testnode(oc, is, no, ocval) ) return 1;
938                 }
939         }
940         else {
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 labda 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                         labdax= (ocx1-ox1+1.0f)*ldx;
953                         dx= 1;
954                 } else if(dox>FLT_EPSILON) {
955                         ldx= 1.0f/dox;
956                         labdax= (ox1-ocx1)*ldx;
957                         dx= -1;
958                 } else {
959                         labdax=1.0f;
960                         ldx=0;
961                         dx= 0;
962                 }
963
964                 if(doy<-FLT_EPSILON) {
965                         ldy= -1.0f/doy;
966                         labday= (ocy1-oy1+1.0f)*ldy;
967                         dy= 1;
968                 } else if(doy>FLT_EPSILON) {
969                         ldy= 1.0f/doy;
970                         labday= (oy1-ocy1)*ldy;
971                         dy= -1;
972                 } else {
973                         labday=1.0f;
974                         ldy=0;
975                         dy= 0;
976                 }
977
978                 if(doz<-FLT_EPSILON) {
979                         ldz= -1.0f/doz;
980                         labdaz= (ocz1-oz1+1.0f)*ldz;
981                         dz= 1;
982                 } else if(doz>FLT_EPSILON) {
983                         ldz= 1.0f/doz;
984                         labdaz= (oz1-ocz1)*ldz;
985                         dz= -1;
986                 } else {
987                         labdaz=1.0f;
988                         ldz=0;
989                         dz= 0;
990                 }
991                 
992                 xo=ocx1; yo=ocy1; zo=ocz1;
993                 labdao= ddalabda= MIN3(labdax,labday,labdaz);
994                 
995                 vec2[0]= ox1;
996                 vec2[1]= oy1;
997                 vec2[2]= oz1;
998                 
999                 /* this loop has been constructed to make sure the first and last node of ray
1000                    are always included, even when ddalabda==1.0f or larger */
1001
1002                 while(TRUE) {
1003
1004                         no= ocread(oc, xo, yo, zo);
1005                         if(no) {
1006                                 
1007                                 /* calculate ray intersection with octree node */
1008                                 VECCOPY(vec1, vec2);
1009                                 // dox,y,z is negative
1010                                 vec2[0]= ox1-ddalabda*dox;
1011                                 vec2[1]= oy1-ddalabda*doy;
1012                                 vec2[2]= oz1-ddalabda*doz;
1013                                 calc_ocval_ray(&ocval, (float)xo, (float)yo, (float)zo, vec1, vec2);
1014                                                            
1015 //                              is->ddalabda= ddalabda;
1016                                 if( testnode(oc, is, no, ocval) ) return 1;
1017                         }
1018
1019                         labdao= ddalabda;
1020                         
1021                         /* traversing ocree nodes need careful detection of smallest values, with proper
1022                            exceptions for equal labdas */
1023                         eqval= (labdax==labday);
1024                         if(labday==labdaz) eqval += 2;
1025                         if(labdax==labdaz) eqval += 4;
1026                         
1027                         if(eqval) {     // only 4 cases exist!
1028                                 if(eqval==7) {  // x=y=z
1029                                         xo+=dx; labdax+=ldx;
1030                                         yo+=dy; labday+=ldy;
1031                                         zo+=dz; labdaz+=ldz;
1032                                 }
1033                                 else if(eqval==1) { // x=y 
1034                                         if(labday < labdaz) {
1035                                                 xo+=dx; labdax+=ldx;
1036                                                 yo+=dy; labday+=ldy;
1037                                         }
1038                                         else {
1039                                                 zo+=dz; labdaz+=ldz;
1040                                         }
1041                                 }
1042                                 else if(eqval==2) { // y=z
1043                                         if(labdax < labday) {
1044                                                 xo+=dx; labdax+=ldx;
1045                                         }
1046                                         else {
1047                                                 yo+=dy; labday+=ldy;
1048                                                 zo+=dz; labdaz+=ldz;
1049                                         }
1050                                 }
1051                                 else { // x=z
1052                                         if(labday < labdax) {
1053                                                 yo+=dy; labday+=ldy;
1054                                         }
1055                                         else {
1056                                                 xo+=dx; labdax+=ldx;
1057                                                 zo+=dz; labdaz+=ldz;
1058                                         }
1059                                 }
1060                         }
1061                         else {  // all three different, just three cases exist
1062                                 eqval= (labdax<labday);
1063                                 if(labday<labdaz) eqval += 2;
1064                                 if(labdax<labdaz) eqval += 4;
1065                                 
1066                                 if(eqval==7 || eqval==5) { // x smallest
1067                                         xo+=dx; labdax+=ldx;
1068                                 }
1069                                 else if(eqval==2 || eqval==6) { // y smallest
1070                                         yo+=dy; labday+=ldy;
1071                                 }
1072                                 else { // z smallest
1073                                         zo+=dz; labdaz+=ldz;
1074                                 }
1075                                 
1076                         }
1077
1078                         ddalabda=MIN3(labdax,labday,labdaz);
1079                         if(ddalabda==labdao) break;
1080                         /* to make sure the last node is always checked */
1081                         if(labdao>=1.0f) break;
1082                 }
1083         }
1084         
1085         /* reached end, no intersections found */
1086         is->hit.ob   = 0;
1087         is->hit.face = NULL;
1088         return 0;
1089 }       
1090
1091
1092