f86032b5d686380f6bd7dad83ef34dc8d80c4b68
[blender.git] / source / blender / render / intern / source / raytrace.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
36 #include "MEM_guardedalloc.h"
37
38 #include "DNA_material_types.h"
39
40 #include "BKE_utildefines.h"
41
42 #include "BLI_arithb.h"
43
44 #include "RE_raytrace.h"
45
46 /* ********** structs *************** */
47
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         struct Branch **adrbranch;
70         struct Node **adrnode;
71         float ocsize;   /* ocsize: mult factor,  max size octree */
72         float ocfacx,ocfacy,ocfacz;
73         float min[3], max[3];
74         int ocres;
75         int branchcount, nodecount;
76         char *ocface; /* during building only */
77         RayCoordsFunc coordsfunc;
78         RayCheckFunc checkfunc;
79 } Octree;
80
81 /* ******** globals ***************** */
82
83 /* just for statistics */
84 static int raycount;
85 static int accepted, rejected, coherent_ray;
86
87
88 /* **************** ocval method ******************* */
89 /* within one octree node, a set of 3x15 bits defines a 'boundbox' to OR with */
90
91 #define OCVALRES        15
92 #define BROW16(min, max)      (((max)>=OCVALRES? 0xFFFF: (1<<(max+1))-1) - ((min>0)? ((1<<(min))-1):0) )
93
94 static void calc_ocval_face(float *v1, float *v2, float *v3, float *v4, short x, short y, short z, OcVal *ov)
95 {
96         float min[3], max[3];
97         int ocmin, ocmax;
98         
99         VECCOPY(min, v1);
100         VECCOPY(max, v1);
101         DO_MINMAX(v2, min, max);
102         DO_MINMAX(v3, min, max);
103         if(v4) {
104                 DO_MINMAX(v4, min, max);
105         }
106         
107         ocmin= OCVALRES*(min[0]-x); 
108         ocmax= OCVALRES*(max[0]-x);
109         ov->ocx= BROW16(ocmin, ocmax);
110         
111         ocmin= OCVALRES*(min[1]-y); 
112         ocmax= OCVALRES*(max[1]-y);
113         ov->ocy= BROW16(ocmin, ocmax);
114         
115         ocmin= OCVALRES*(min[2]-z); 
116         ocmax= OCVALRES*(max[2]-z);
117         ov->ocz= BROW16(ocmin, ocmax);
118
119 }
120
121 static void calc_ocval_ray(OcVal *ov, float xo, float yo, float zo, float *vec1, float *vec2)
122 {
123         int ocmin, ocmax;
124         
125         if(vec1[0]<vec2[0]) {
126                 ocmin= OCVALRES*(vec1[0] - xo);
127                 ocmax= OCVALRES*(vec2[0] - xo);
128         } else {
129                 ocmin= OCVALRES*(vec2[0] - xo);
130                 ocmax= OCVALRES*(vec1[0] - xo);
131         }
132         ov->ocx= BROW16(ocmin, ocmax);
133
134         if(vec1[1]<vec2[1]) {
135                 ocmin= OCVALRES*(vec1[1] - yo);
136                 ocmax= OCVALRES*(vec2[1] - yo);
137         } else {
138                 ocmin= OCVALRES*(vec2[1] - yo);
139                 ocmax= OCVALRES*(vec1[1] - yo);
140         }
141         ov->ocy= BROW16(ocmin, ocmax);
142
143         if(vec1[2]<vec2[2]) {
144                 ocmin= OCVALRES*(vec1[2] - zo);
145                 ocmax= OCVALRES*(vec2[2] - zo);
146         } else {
147                 ocmin= OCVALRES*(vec2[2] - zo);
148                 ocmax= OCVALRES*(vec1[2] - zo);
149         }
150         ov->ocz= BROW16(ocmin, ocmax);
151 }
152
153 /* ************* octree ************** */
154
155 static Branch *addbranch(Octree *oc, Branch *br, short ocb)
156 {
157         int index;
158         
159         if(br->b[ocb]) return br->b[ocb];
160         
161         oc->branchcount++;
162         index= oc->branchcount>>12;
163         
164         if(oc->adrbranch[index]==NULL)
165                 oc->adrbranch[index]= MEM_callocN(4096*sizeof(Branch), "new oc branch");
166
167         if(oc->branchcount>= BRANCH_ARRAY*4096) {
168                 printf("error; octree branches full\n");
169                 oc->branchcount=0;
170         }
171         
172         return br->b[ocb]= oc->adrbranch[index]+(oc->branchcount & 4095);
173 }
174
175 static Node *addnode(Octree *oc)
176 {
177         int index;
178         
179         oc->nodecount++;
180         index= oc->nodecount>>12;
181         
182         if(oc->adrnode[index]==NULL)
183                 oc->adrnode[index]= MEM_callocN(4096*sizeof(Node),"addnode");
184
185         if(oc->nodecount> NODE_ARRAY*NODE_ARRAY) {
186                 printf("error; octree nodes full\n");
187                 oc->nodecount=0;
188         }
189         
190         return oc->adrnode[index]+(oc->nodecount & 4095);
191 }
192
193 static int face_in_node(RayFace *face, short x, short y, short z, float rtf[][3])
194 {
195         static float nor[3], d;
196         float fx, fy, fz;
197         
198         // init static vars 
199         if(face) {
200                 CalcNormFloat(rtf[0], rtf[1], rtf[2], nor);
201                 d= -nor[0]*rtf[0][0] - nor[1]*rtf[0][1] - nor[2]*rtf[0][2];
202                 return 0;
203         }
204         
205         fx= x;
206         fy= y;
207         fz= z;
208         
209         if((fx)*nor[0] + (fy)*nor[1] + (fz)*nor[2] + d > 0.0f) {
210                 if((fx+1)*nor[0] + (fy  )*nor[1] + (fz  )*nor[2] + d < 0.0f) return 1;
211                 if((fx  )*nor[0] + (fy+1)*nor[1] + (fz  )*nor[2] + d < 0.0f) return 1;
212                 if((fx+1)*nor[0] + (fy+1)*nor[1] + (fz  )*nor[2] + d < 0.0f) return 1;
213         
214                 if((fx  )*nor[0] + (fy  )*nor[1] + (fz+1)*nor[2] + d < 0.0f) return 1;
215                 if((fx+1)*nor[0] + (fy  )*nor[1] + (fz+1)*nor[2] + d < 0.0f) return 1;
216                 if((fx  )*nor[0] + (fy+1)*nor[1] + (fz+1)*nor[2] + d < 0.0f) return 1;
217                 if((fx+1)*nor[0] + (fy+1)*nor[1] + (fz+1)*nor[2] + d < 0.0f) return 1;
218         }
219         else {
220                 if((fx+1)*nor[0] + (fy  )*nor[1] + (fz  )*nor[2] + d > 0.0f) return 1;
221                 if((fx  )*nor[0] + (fy+1)*nor[1] + (fz  )*nor[2] + d > 0.0f) return 1;
222                 if((fx+1)*nor[0] + (fy+1)*nor[1] + (fz  )*nor[2] + d > 0.0f) return 1;
223         
224                 if((fx  )*nor[0] + (fy  )*nor[1] + (fz+1)*nor[2] + d > 0.0f) return 1;
225                 if((fx+1)*nor[0] + (fy  )*nor[1] + (fz+1)*nor[2] + d > 0.0f) return 1;
226                 if((fx  )*nor[0] + (fy+1)*nor[1] + (fz+1)*nor[2] + d > 0.0f) return 1;
227                 if((fx+1)*nor[0] + (fy+1)*nor[1] + (fz+1)*nor[2] + d > 0.0f) return 1;
228         }
229         
230         return 0;
231 }
232
233 static void ocwrite(Octree *oc, RayFace *face, int quad, short x, short y, short z, float rtf[][3])
234 {
235         Branch *br;
236         Node *no;
237         short a, oc0, oc1, oc2, oc3, oc4, oc5;
238
239         x<<=2;
240         y<<=1;
241
242         br= oc->adrbranch[0];
243
244         if(oc->ocres==512) {
245                 oc0= ((x & 1024)+(y & 512)+(z & 256))>>8;
246                 br= addbranch(oc, br, oc0);
247         }
248         if(oc->ocres>=256) {
249                 oc0= ((x & 512)+(y & 256)+(z & 128))>>7;
250                 br= addbranch(oc, br, oc0);
251         }
252         if(oc->ocres>=128) {
253                 oc0= ((x & 256)+(y & 128)+(z & 64))>>6;
254                 br= addbranch(oc, br, oc0);
255         }
256
257         oc0= ((x & 128)+(y & 64)+(z & 32))>>5;
258         oc1= ((x & 64)+(y & 32)+(z & 16))>>4;
259         oc2= ((x & 32)+(y & 16)+(z & 8))>>3;
260         oc3= ((x & 16)+(y & 8)+(z & 4))>>2;
261         oc4= ((x & 8)+(y & 4)+(z & 2))>>1;
262         oc5= ((x & 4)+(y & 2)+(z & 1));
263
264         br= addbranch(oc, br,oc0);
265         br= addbranch(oc, br,oc1);
266         br= addbranch(oc, br,oc2);
267         br= addbranch(oc, br,oc3);
268         br= addbranch(oc, br,oc4);
269         no= (Node *)br->b[oc5];
270         if(no==NULL) br->b[oc5]= (Branch *)(no= addnode(oc));
271
272         while(no->next) no= no->next;
273
274         a= 0;
275         if(no->v[7]) {          /* node full */
276                 no->next= addnode(oc);
277                 no= no->next;
278         }
279         else {
280                 while(no->v[a]!=NULL) a++;
281         }
282         
283         no->v[a]= face;
284         
285         if(quad)
286                 calc_ocval_face(rtf[0], rtf[1], rtf[2], rtf[3], x>>2, y>>1, z, &no->ov[a]);
287         else
288                 calc_ocval_face(rtf[0], rtf[1], rtf[2], NULL, x>>2, y>>1, z, &no->ov[a]);
289 }
290
291 static void d2dda(Octree *oc, short b1, short b2, short c1, short c2, char *ocface, short rts[][3], float rtf[][3])
292 {
293         int ocx1,ocx2,ocy1,ocy2;
294         int x,y,dx=0,dy=0;
295         float ox1,ox2,oy1,oy2;
296         float labda,labdao,labdax,labday,ldx,ldy;
297
298         ocx1= rts[b1][c1];
299         ocy1= rts[b1][c2];
300         ocx2= rts[b2][c1];
301         ocy2= rts[b2][c2];
302
303         if(ocx1==ocx2 && ocy1==ocy2) {
304                 ocface[oc->ocres*ocx1+ocy1]= 1;
305                 return;
306         }
307
308         ox1= rtf[b1][c1];
309         oy1= rtf[b1][c2];
310         ox2= rtf[b2][c1];
311         oy2= rtf[b2][c2];
312
313         if(ox1!=ox2) {
314                 if(ox2-ox1>0.0f) {
315                         labdax= (ox1-ocx1-1.0f)/(ox1-ox2);
316                         ldx= -1.0f/(ox1-ox2);
317                         dx= 1;
318                 } else {
319                         labdax= (ox1-ocx1)/(ox1-ox2);
320                         ldx= 1.0f/(ox1-ox2);
321                         dx= -1;
322                 }
323         } else {
324                 labdax=1.0f;
325                 ldx=0;
326         }
327
328         if(oy1!=oy2) {
329                 if(oy2-oy1>0.0f) {
330                         labday= (oy1-ocy1-1.0f)/(oy1-oy2);
331                         ldy= -1.0f/(oy1-oy2);
332                         dy= 1;
333                 } else {
334                         labday= (oy1-ocy1)/(oy1-oy2);
335                         ldy= 1.0f/(oy1-oy2);
336                         dy= -1;
337                 }
338         } else {
339                 labday=1.0f;
340                 ldy=0;
341         }
342         
343         x=ocx1; y=ocy1;
344         labda= MIN2(labdax, labday);
345         
346         while(TRUE) {
347                 
348                 if(x<0 || y<0 || x>=oc->ocres || y>=oc->ocres);
349                 else ocface[oc->ocres*x+y]= 1;
350                 
351                 labdao=labda;
352                 if(labdax==labday) {
353                         labdax+=ldx;
354                         x+=dx;
355                         labday+=ldy;
356                         y+=dy;
357                 } else {
358                         if(labdax<labday) {
359                                 labdax+=ldx;
360                                 x+=dx;
361                         } else {
362                                 labday+=ldy;
363                                 y+=dy;
364                         }
365                 }
366                 labda=MIN2(labdax,labday);
367                 if(labda==labdao) break;
368                 if(labda>=1.0f) break;
369         }
370         ocface[oc->ocres*ocx2+ocy2]=1;
371 }
372
373 static void filltriangle(Octree *oc, short c1, short c2, char *ocface, short *ocmin)
374 {
375         short *ocmax;
376         int a, x, y, y1, y2;
377
378         ocmax=ocmin+3;
379
380         for(x=ocmin[c1];x<=ocmax[c1];x++) {
381                 a= oc->ocres*x;
382                 for(y=ocmin[c2];y<=ocmax[c2];y++) {
383                         if(ocface[a+y]) {
384                                 y++;
385                                 while(ocface[a+y] && y!=ocmax[c2]) y++;
386                                 for(y1=ocmax[c2];y1>y;y1--) {
387                                         if(ocface[a+y1]) {
388                                                 for(y2=y;y2<=y1;y2++) ocface[a+y2]=1;
389                                                 y1=0;
390                                         }
391                                 }
392                                 y=ocmax[c2];
393                         }
394                 }
395         }
396 }
397
398 void RE_ray_tree_free(RayTree *tree)
399 {
400         Octree *oc= (Octree*)tree;
401
402 #if 0
403         printf("branches %d nodes %d\n", oc->branchcount, oc->nodecount);
404         printf("raycount %d \n", raycount);     
405         printf("ray coherent %d \n", coherent_ray);
406         printf("accepted %d rejected %d\n", accepted, rejected);
407 #endif
408         if(oc->ocface)
409                 MEM_freeN(oc->ocface);
410
411         if(oc->adrbranch) {
412                 int a= 0;
413                 while(oc->adrbranch[a]) {
414                         MEM_freeN(oc->adrbranch[a]);
415                         oc->adrbranch[a]= NULL;
416                         a++;
417                 }
418                 MEM_freeN(oc->adrbranch);
419                 oc->adrbranch= NULL;
420         }
421         oc->branchcount= 0;
422         
423         if(oc->adrnode) {
424                 int a= 0;
425                 while(oc->adrnode[a]) {
426                         MEM_freeN(oc->adrnode[a]);
427                         oc->adrnode[a]= NULL;
428                         a++;
429                 }
430                 MEM_freeN(oc->adrnode);
431                 oc->adrnode= NULL;
432         }
433         oc->nodecount= 0;
434
435         MEM_freeN(oc);
436 }
437
438 RayTree *RE_ray_tree_create(int ocres, int totface, float *min, float *max, RayCoordsFunc coordsfunc, RayCheckFunc checkfunc)
439 {
440         Octree *oc;
441         float t00, t01, t02;
442         int c, ocres2;
443         
444         oc= MEM_callocN(sizeof(Octree), "Octree");
445         oc->adrbranch= MEM_callocN(sizeof(void *)*BRANCH_ARRAY, "octree branches");
446         oc->adrnode= MEM_callocN(sizeof(void *)*NODE_ARRAY, "octree nodes");
447
448         oc->coordsfunc= coordsfunc;
449         oc->checkfunc= checkfunc;
450
451         /* only for debug info */
452         raycount=0;
453         accepted= 0;
454         rejected= 0;
455         coherent_ray= 0;
456         
457         /* fill main octree struct */
458         oc->ocres= ocres;
459         ocres2= oc->ocres*oc->ocres;
460         
461         VECCOPY(oc->min, min);
462         VECCOPY(oc->max, max);
463
464         oc->adrbranch[0]=(Branch *)MEM_callocN(4096*sizeof(Branch), "makeoctree");
465         
466         /* the lookup table, per face, for which nodes to fill in */
467         oc->ocface= MEM_callocN( 3*ocres2 + 8, "ocface");
468         memset(oc->ocface, 0, 3*ocres2);
469
470         for(c=0;c<3;c++) {      /* octree enlarge, still needed? */
471                 oc->min[c]-= 0.01f;
472                 oc->max[c]+= 0.01f;
473         }
474
475         t00= oc->max[0]-oc->min[0];
476         t01= oc->max[1]-oc->min[1];
477         t02= oc->max[2]-oc->min[2];
478         
479         /* this minus 0.1 is old safety... seems to be needed? */
480         oc->ocfacx= (oc->ocres-0.1)/t00;
481         oc->ocfacy= (oc->ocres-0.1)/t01;
482         oc->ocfacz= (oc->ocres-0.1)/t02;
483         
484         oc->ocsize= sqrt(t00*t00+t01*t01+t02*t02);      /* global, max size octree */
485
486         return (RayTree*)oc;
487 }
488
489 void RE_ray_tree_add_face(RayTree *tree, RayFace *face)
490 {
491         Octree *oc = (Octree*)tree;
492         float *v1, *v2, *v3, *v4, ocfac[3], rtf[4][3];
493         short rts[4][3], ocmin[6], *ocmax;
494         char *ocface= oc->ocface;       // front, top, size view of face, to fill in
495         int a, b, c, oc1, oc2, oc3, oc4, x, y, z, ocres2;
496
497         ocfac[0]= oc->ocfacx;
498         ocfac[1]= oc->ocfacy;
499         ocfac[2]= oc->ocfacz;
500
501         ocres2= oc->ocres*oc->ocres;
502
503         ocmax= ocmin+3;
504
505         oc->coordsfunc(face, &v1, &v2, &v3, &v4);
506         
507         for(c=0;c<3;c++) {
508                 rtf[0][c]= (v1[c]-oc->min[c])*ocfac[c] ;
509                 rts[0][c]= (short)rtf[0][c];
510                 rtf[1][c]= (v2[c]-oc->min[c])*ocfac[c] ;
511                 rts[1][c]= (short)rtf[1][c];
512                 rtf[2][c]= (v3[c]-oc->min[c])*ocfac[c] ;
513                 rts[2][c]= (short)rtf[2][c];
514                 if(v4) {
515                         rtf[3][c]= (v4[c]-oc->min[c])*ocfac[c] ;
516                         rts[3][c]= (short)rtf[3][c];
517                 }
518         }
519         
520         for(c=0;c<3;c++) {
521                 oc1= rts[0][c];
522                 oc2= rts[1][c];
523                 oc3= rts[2][c];
524                 if(v4==NULL) {
525                         ocmin[c]= MIN3(oc1,oc2,oc3);
526                         ocmax[c]= MAX3(oc1,oc2,oc3);
527                 }
528                 else {
529                         oc4= rts[3][c];
530                         ocmin[c]= MIN4(oc1,oc2,oc3,oc4);
531                         ocmax[c]= MAX4(oc1,oc2,oc3,oc4);
532                 }
533                 if(ocmax[c]>oc->ocres-1) ocmax[c]=oc->ocres-1;
534                 if(ocmin[c]<0) ocmin[c]=0;
535         }
536         
537         if(ocmin[0]==ocmax[0] && ocmin[1]==ocmax[1] && ocmin[2]==ocmax[2]) {
538                 ocwrite(oc, face, (v4 != NULL), ocmin[0], ocmin[1], ocmin[2], rtf);
539         }
540         else {
541                 
542                 d2dda(oc, 0,1,0,1,ocface+ocres2,rts,rtf);
543                 d2dda(oc, 0,1,0,2,ocface,rts,rtf);
544                 d2dda(oc, 0,1,1,2,ocface+2*ocres2,rts,rtf);
545                 d2dda(oc, 1,2,0,1,ocface+ocres2,rts,rtf);
546                 d2dda(oc, 1,2,0,2,ocface,rts,rtf);
547                 d2dda(oc, 1,2,1,2,ocface+2*ocres2,rts,rtf);
548                 if(v4==NULL) {
549                         d2dda(oc, 2,0,0,1,ocface+ocres2,rts,rtf);
550                         d2dda(oc, 2,0,0,2,ocface,rts,rtf);
551                         d2dda(oc, 2,0,1,2,ocface+2*ocres2,rts,rtf);
552                 }
553                 else {
554                         d2dda(oc, 2,3,0,1,ocface+ocres2,rts,rtf);
555                         d2dda(oc, 2,3,0,2,ocface,rts,rtf);
556                         d2dda(oc, 2,3,1,2,ocface+2*ocres2,rts,rtf);
557                         d2dda(oc, 3,0,0,1,ocface+ocres2,rts,rtf);
558                         d2dda(oc, 3,0,0,2,ocface,rts,rtf);
559                         d2dda(oc, 3,0,1,2,ocface+2*ocres2,rts,rtf);
560                 }
561                 /* nothing todo with triangle..., just fills :) */
562                 filltriangle(oc, 0,1,ocface+ocres2,ocmin);
563                 filltriangle(oc, 0,2,ocface,ocmin);
564                 filltriangle(oc, 1,2,ocface+2*ocres2,ocmin);
565                 
566                 /* init static vars here */
567                 face_in_node(face, 0,0,0, rtf);
568                 
569                 for(x=ocmin[0];x<=ocmax[0];x++) {
570                         a= oc->ocres*x;
571                         for(y=ocmin[1];y<=ocmax[1];y++) {
572                                 if(ocface[a+y+ocres2]) {
573                                         b= oc->ocres*y+2*ocres2;
574                                         for(z=ocmin[2];z<=ocmax[2];z++) {
575                                                 if(ocface[b+z] && ocface[a+z]) {
576                                                         if(face_in_node(NULL, x, y, z, rtf))
577                                                                 ocwrite(oc, face, (v4 != NULL), x,y,z, rtf);
578                                                 }
579                                         }
580                                 }
581                         }
582                 }
583                 
584                 /* same loops to clear octree, doubt it can be done smarter */
585                 for(x=ocmin[0];x<=ocmax[0];x++) {
586                         a= oc->ocres*x;
587                         for(y=ocmin[1];y<=ocmax[1];y++) {
588                                 /* x-y */
589                                 ocface[a+y+ocres2]= 0;
590
591                                 b= oc->ocres*y + 2*ocres2;
592                                 for(z=ocmin[2];z<=ocmax[2];z++) {
593                                         /* y-z */
594                                         ocface[b+z]= 0;
595                                         /* x-z */
596                                         ocface[a+z]= 0;
597                                 }
598                         }
599                 }
600         }
601 }
602
603 void RE_ray_tree_done(RayTree *tree)
604 {
605         Octree *oc= (Octree*)tree;
606
607         MEM_freeN(oc->ocface);
608         oc->ocface= NULL;
609 }
610
611 /* ************ raytracer **************** */
612
613 /* only for self-intersecting test with current render face (where ray left) */
614 static int intersection2(RayFace *face, RayCoordsFunc coordsfunc, float r0, float r1, float r2, float rx1, float ry1, float rz1)
615 {
616         float *v1, *v2, *v3, *v4;
617         float x0,x1,x2,t00,t01,t02,t10,t11,t12,t20,t21,t22;
618         float m0, m1, m2, divdet, det, det1;
619         float u1, v, u2;
620
621         coordsfunc(face, &v1, &v2, &v3, &v4);
622
623         /* happens for baking with non existing face */
624         if(v1==NULL)
625                 return 1;
626         
627         if (v4) {
628                 SWAP(float*, v3, v4);
629         }
630         
631         t00= v3[0]-v1[0];
632         t01= v3[1]-v1[1];
633         t02= v3[2]-v1[2];
634         t10= v3[0]-v2[0];
635         t11= v3[1]-v2[1];
636         t12= v3[2]-v2[2];
637         
638         x0= t11*r2-t12*r1;
639         x1= t12*r0-t10*r2;
640         x2= t10*r1-t11*r0;
641
642         divdet= t00*x0+t01*x1+t02*x2;
643
644         m0= rx1-v3[0];
645         m1= ry1-v3[1];
646         m2= rz1-v3[2];
647         det1= m0*x0+m1*x1+m2*x2;
648         
649         if(divdet!=0.0f) {
650                 u1= det1/divdet;
651
652                 if(u1<=0.0f) {
653                         det= t00*(m1*r2-m2*r1);
654                         det+= t01*(m2*r0-m0*r2);
655                         det+= t02*(m0*r1-m1*r0);
656                         v= det/divdet;
657
658                         if(v<=0.0f && (u1 + v) >= -1.0f) {
659                                 return 1;
660                         }
661                 }
662         }
663
664         if(v4) {
665
666                 t20= v3[0]-v4[0];
667                 t21= v3[1]-v4[1];
668                 t22= v3[2]-v4[2];
669
670                 divdet= t20*x0+t21*x1+t22*x2;
671                 if(divdet!=0.0f) {
672                         u2= det1/divdet;
673                 
674                         if(u2<=0.0f) {
675                                 det= t20*(m1*r2-m2*r1);
676                                 det+= t21*(m2*r0-m0*r2);
677                                 det+= t22*(m0*r1-m1*r0);
678                                 v= det/divdet;
679         
680                                 if(v<=0.0f && (u2 + v) >= -1.0f) {
681                                         return 2;
682                                 }
683                         }
684                 }
685         }
686         return 0;
687 }
688
689 #if 0
690 /* ray - line intersection */
691 /* disabled until i got real & fast cylinder checking, this code doesnt work proper
692 for faster strands */
693
694 static int intersection_strand(Isect *is)
695 {
696         float v1[3], v2[3];             /* length of strand */
697         float axis[3], rc[3], nor[3], radline, dist, len;
698         
699         /* radius strand */
700         radline= 0.5f*VecLenf(is->vlr->v1->co, is->vlr->v2->co);
701         
702         VecMidf(v1, is->vlr->v1->co, is->vlr->v2->co);
703         VecMidf(v2, is->vlr->v3->co, is->vlr->v4->co);
704         
705         VECSUB(rc, v1, is->start);      /* vector from base ray to base cylinder */
706         VECSUB(axis, v2, v1);           /* cylinder axis */
707         
708         CROSS(nor, is->vec, axis);
709         len= VecLength(nor);
710         
711         if(len<FLT_EPSILON)
712                 return 0;
713
714         dist= INPR(rc, nor)/len;        /* distance between ray and axis cylinder */
715         
716         if(dist<radline && dist>-radline) {
717                 float dot1, dot2, dot3, rlen, alen, div;
718                 float labda;
719                 
720                 /* calculating the intersection point of shortest distance */
721                 dot1 = INPR(rc, is->vec);
722                 dot2 = INPR(is->vec, axis);
723                 dot3 = INPR(rc, axis);
724                 rlen = INPR(is->vec, is->vec);
725                 alen = INPR(axis, axis);
726                 
727                 div = alen * rlen - dot2 * dot2;
728                 if (ABS(div) < FLT_EPSILON)
729                         return 0;
730                 
731                 labda = (dot1*dot2 - dot3*rlen)/div;
732                 
733                 radline/= sqrt(alen);
734                 
735                 /* labda: where on axis do we have closest intersection? */
736                 if(labda >= -radline && labda <= 1.0f+radline) {
737                         VlakRen *vlr= is->faceorig;
738                         VertRen *v1= is->vlr->v1, *v2= is->vlr->v2, *v3= is->vlr->v3, *v4= is->vlr->v4;
739                                 /* but we dont do shadows from faces sharing edge */
740                         
741                         if(v1==vlr->v1 || v2==vlr->v1 || v3==vlr->v1 || v4==vlr->v1) return 0;
742                         if(v1==vlr->v2 || v2==vlr->v2 || v3==vlr->v2 || v4==vlr->v2) return 0;
743                         if(v1==vlr->v3 || v2==vlr->v3 || v3==vlr->v3 || v4==vlr->v3) return 0;
744                         if(vlr->v4) {
745                                 if(v1==vlr->v4 || v2==vlr->v4 || v3==vlr->v4 || v4==vlr->v4) return 0;
746                         }
747                         return 1;
748                 }
749         }
750         return 0;
751 }
752 #endif
753
754 /* ray - triangle or quad intersection */
755 int RE_ray_face_intersection(Isect *is, RayCoordsFunc coordsfunc)
756 {
757         RayFace *face= is->face;
758         float *v1,*v2,*v3,*v4;
759         float x0,x1,x2,t00,t01,t02,t10,t11,t12,t20,t21,t22,r0,r1,r2;
760         float m0, m1, m2, divdet, det1;
761         short ok=0;
762
763         /* disabled until i got real & fast cylinder checking, this code doesnt work proper
764            for faster strands */
765 //      if(is->mode==RE_RAY_SHADOW && is->vlr->flag & R_STRAND) 
766 //              return intersection_strand(is);
767         
768         coordsfunc(face, &v1, &v2, &v3, &v4);
769
770         if (v4) {
771                 SWAP(float*, v3, v4);
772         }
773
774         t00= v3[0]-v1[0];
775         t01= v3[1]-v1[1];
776         t02= v3[2]-v1[2];
777         t10= v3[0]-v2[0];
778         t11= v3[1]-v2[1];
779         t12= v3[2]-v2[2];
780         
781         r0= is->vec[0];
782         r1= is->vec[1];
783         r2= is->vec[2];
784         
785         x0= t12*r1-t11*r2;
786         x1= t10*r2-t12*r0;
787         x2= t11*r0-t10*r1;
788
789         divdet= t00*x0+t01*x1+t02*x2;
790
791         m0= is->start[0]-v3[0];
792         m1= is->start[1]-v3[1];
793         m2= is->start[2]-v3[2];
794         det1= m0*x0+m1*x1+m2*x2;
795         
796         if(divdet!=0.0f) {
797                 float u;
798
799                 divdet= 1.0f/divdet;
800                 u= det1*divdet;
801                 if(u<0.0f && u>-1.0f) {
802                         float v, cros0, cros1, cros2;
803                         
804                         cros0= m1*t02-m2*t01;
805                         cros1= m2*t00-m0*t02;
806                         cros2= m0*t01-m1*t00;
807                         v= divdet*(cros0*r0 + cros1*r1 + cros2*r2);
808
809                         if(v<0.0f && (u + v) > -1.0f) {
810                                 float labda;
811                                 labda= divdet*(cros0*t10 + cros1*t11 + cros2*t12);
812
813                                 if(labda>0.0f && labda<1.0f) {
814                                         is->labda= labda;
815                                         is->u= u; is->v= v;
816                                         ok= 1;
817                                 }
818                         }
819                 }
820         }
821
822         if(ok==0 && v4) {
823
824                 t20= v3[0]-v4[0];
825                 t21= v3[1]-v4[1];
826                 t22= v3[2]-v4[2];
827
828                 divdet= t20*x0+t21*x1+t22*x2;
829                 if(divdet!=0.0f) {
830                         float u;
831                         divdet= 1.0f/divdet;
832                         u = det1*divdet;
833                         
834                         if(u<0.0f && u>-1.0f) {
835                                 float v, cros0, cros1, cros2;
836                                 cros0= m1*t22-m2*t21;
837                                 cros1= m2*t20-m0*t22;
838                                 cros2= m0*t21-m1*t20;
839                                 v= divdet*(cros0*r0 + cros1*r1 + cros2*r2);
840         
841                                 if(v<0.0f && (u + v) > -1.0f) {
842                                         float labda;
843                                         labda= divdet*(cros0*t10 + cros1*t11 + cros2*t12);
844                                         
845                                         if(labda>0.0f && labda<1.0f) {
846                                                 ok= 2;
847                                                 is->labda= labda;
848                                                 is->u= u; is->v= v;
849                                         }
850                                 }
851                         }
852                 }
853         }
854
855         if(ok) {
856                 is->isect= ok;  // wich half of the quad
857                 
858                 if(is->mode!=RE_RAY_SHADOW) {
859                         /* for mirror & tra-shadow: large faces can be filled in too often, this prevents
860                            a face being detected too soon... */
861                         if(is->labda > is->ddalabda) {
862                                 return 0;
863                         }
864                 }
865                 
866                 /* when a shadow ray leaves a face, it can be little outside the edges of it, causing
867                 intersection to be detected in its neighbour face */
868                 
869                 if(is->facecontr && is->faceisect);     // optimizing, the tests below are not needed
870                 else if(is->labda< .1) {
871                         RayFace *face= is->faceorig;
872                         float *origv1, *origv2, *origv3, *origv4;
873                         short de= 0;
874
875                         coordsfunc(face, &origv1, &origv2, &origv3, &origv4);
876                         
877                         if(v1==origv1 || v2==origv1 || v3==origv1 || v4==origv1) de++;
878                         if(v1==origv2 || v2==origv2 || v3==origv2 || v4==origv2) de++;
879                         if(v1==origv3 || v2==origv3 || v3==origv3 || v4==origv3) de++;
880                         if(origv4) {
881                                 if(v1==origv4 || v2==origv4 || v3==origv4 || v4==origv4) de++;
882                         }
883                         if(de) {
884                                 /* so there's a shared edge or vertex, let's intersect ray with face
885                                 itself, if that's true we can safely return 1, otherwise we assume
886                                 the intersection is invalid, 0 */
887                                 
888                                 if(is->facecontr==NULL) {
889                                         is->facecontr= face;
890                                         is->faceisect= intersection2(face, coordsfunc, -r0, -r1, -r2, is->start[0], is->start[1], is->start[2]);
891                                 }
892
893                                 if(is->faceisect) return 1;
894                                 return 0;
895                         }
896                 }
897                 
898                 return 1;
899         }
900
901         return 0;
902 }
903
904 /* check all faces in this node */
905 static int testnode(Octree *oc, Isect *is, Node *no, OcVal ocval)
906 {
907         RayFace *face;
908         short nr=0;
909         OcVal *ov;
910
911         /* return on any first hit */
912         if(is->mode==RE_RAY_SHADOW) {
913                 
914                 face= no->v[0];
915                 while(face) {
916                 
917                         if(is->faceorig != face) {
918
919                                 if(oc->checkfunc(is, face)) {
920                                         
921                                         ov= no->ov+nr;
922                                         if( (ov->ocx & ocval.ocx) && (ov->ocy & ocval.ocy) && (ov->ocz & ocval.ocz) ) { 
923                                                 //accepted++;
924                                                 is->face= face;
925         
926                                                 if(RE_ray_face_intersection(is, oc->coordsfunc)) {
927                                                         is->face_last= face;
928                                                         return 1;
929                                                 }
930                                         }
931                                         //else rejected++;
932                                 }
933                         }
934                         
935                         nr++;
936                         if(nr==8) {
937                                 no= no->next;
938                                 if(no==0) return 0;
939                                 nr=0;
940                         }
941                         face= no->v[nr];
942                 }
943         }
944         else {                  /* else mirror or glass or shadowtra, return closest face  */
945                 Isect isect;
946                 int found= 0;
947                 
948                 is->labda= 1.0f;        /* needed? */
949                 isect= *is;             /* copy for sorting */
950                 
951                 face= no->v[0];
952                 while(face) {
953
954                         if(is->faceorig != face) {
955                                 if(oc->checkfunc(is, face)) {
956                                         ov= no->ov+nr;
957                                         if( (ov->ocx & ocval.ocx) && (ov->ocy & ocval.ocy) && (ov->ocz & ocval.ocz) ) { 
958                                                 //accepted++;
959
960                                                 isect.face= face;
961                                                 if(RE_ray_face_intersection(&isect, oc->coordsfunc)) {
962                                                         if(isect.labda<is->labda) *is= isect;
963                                                         found= 1;
964                                                 }
965                                         }
966                                         //else rejected++;
967                                 }
968                         }
969                         
970                         nr++;
971                         if(nr==8) {
972                                 no= no->next;
973                                 if(no==NULL) break;
974                                 nr=0;
975                         }
976                         face= no->v[nr];
977                 }
978                 
979                 return found;
980         }
981
982         return 0;
983 }
984
985 /* find the Node for the octree coord x y z */
986 static Node *ocread(Octree *oc, int x, int y, int z)
987 {
988         Branch *br;
989         int oc1;
990         
991         x<<=2;
992         y<<=1;
993         
994         br= oc->adrbranch[0];
995         
996         if(oc->ocres==512) {
997                 oc1= ((x & 1024)+(y & 512)+(z & 256))>>8;
998                 br= br->b[oc1];
999                 if(br==NULL) {
1000                         return NULL;
1001                 }
1002         }
1003         if(oc->ocres>=256) {
1004                 oc1= ((x & 512)+(y & 256)+(z & 128))>>7;
1005                 br= br->b[oc1];
1006                 if(br==NULL) {
1007                         return NULL;
1008                 }
1009         }
1010         if(oc->ocres>=128) {
1011                 oc1= ((x & 256)+(y & 128)+(z & 64))>>6;
1012                 br= br->b[oc1];
1013                 if(br==NULL) {
1014                         return NULL;
1015                 }
1016         }
1017         
1018         oc1= ((x & 128)+(y & 64)+(z & 32))>>5;
1019         br= br->b[oc1];
1020         if(br) {
1021                 oc1= ((x & 64)+(y & 32)+(z & 16))>>4;
1022                 br= br->b[oc1];
1023                 if(br) {
1024                         oc1= ((x & 32)+(y & 16)+(z & 8))>>3;
1025                         br= br->b[oc1];
1026                         if(br) {
1027                                 oc1= ((x & 16)+(y & 8)+(z & 4))>>2;
1028                                 br= br->b[oc1];
1029                                 if(br) {
1030                                         oc1= ((x & 8)+(y & 4)+(z & 2))>>1;
1031                                         br= br->b[oc1];
1032                                         if(br) {
1033                                                 oc1= ((x & 4)+(y & 2)+(z & 1));
1034                                                 return (Node *)br->b[oc1];
1035                                         }
1036                                 }
1037                         }
1038                 }
1039         }
1040         
1041         return NULL;
1042 }
1043
1044 static int cliptest(float p, float q, float *u1, float *u2)
1045 {
1046         float r;
1047
1048         if(p<0.0f) {
1049                 if(q<p) return 0;
1050                 else if(q<0.0f) {
1051                         r= q/p;
1052                         if(r>*u2) return 0;
1053                         else if(r>*u1) *u1=r;
1054                 }
1055         }
1056         else {
1057                 if(p>0.0f) {
1058                         if(q<0.0f) return 0;
1059                         else if(q<p) {
1060                                 r= q/p;
1061                                 if(r<*u1) return 0;
1062                                 else if(r<*u2) *u2=r;
1063                         }
1064                 }
1065                 else if(q<0.0f) return 0;
1066         }
1067         return 1;
1068 }
1069
1070 /* extensive coherence checks/storage cancels out the benefit of it, and gives errors... we
1071    need better methods, sample code commented out below (ton) */
1072  
1073 /*
1074
1075 in top: static int coh_nodes[16*16*16][6];
1076 in makeoctree: memset(coh_nodes, 0, sizeof(coh_nodes));
1077  
1078 static void add_coherence_test(int ocx1, int ocx2, int ocy1, int ocy2, int ocz1, int ocz2)
1079 {
1080         short *sp;
1081         
1082         sp= coh_nodes[ (ocx2 & 15) + 16*(ocy2 & 15) + 256*(ocz2 & 15) ];
1083         sp[0]= ocx1; sp[1]= ocy1; sp[2]= ocz1;
1084         sp[3]= ocx2; sp[4]= ocy2; sp[5]= ocz2;
1085         
1086 }
1087
1088 static int do_coherence_test(int ocx1, int ocx2, int ocy1, int ocy2, int ocz1, int ocz2)
1089 {
1090         short *sp;
1091         
1092         sp= coh_nodes[ (ocx2 & 15) + 16*(ocy2 & 15) + 256*(ocz2 & 15) ];
1093         if(sp[0]==ocx1 && sp[1]==ocy1 && sp[2]==ocz1 &&
1094            sp[3]==ocx2 && sp[4]==ocy2 && sp[5]==ocz2) return 1;
1095         return 0;
1096 }
1097
1098 */
1099
1100 /* return 1: found valid intersection */
1101 /* starts with is->faceorig */
1102 int RE_ray_tree_intersect(RayTree *tree, Isect *is)     
1103 {
1104         Octree *oc= (Octree*)tree;
1105         Node *no;
1106         OcVal ocval;
1107         float vec1[3], vec2[3];
1108         float u1,u2,ox1,ox2,oy1,oy2,oz1,oz2;
1109         float labdao,labdax,ldx,labday,ldy,labdaz,ldz, ddalabda;
1110         int dx,dy,dz;   
1111         int xo,yo,zo,c1=0;
1112         int ocx1,ocx2,ocy1, ocy2,ocz1,ocz2;
1113         
1114         /* clip with octree */
1115         if(oc->branchcount==0) return 0;
1116         
1117         /* do this before intersect calls */
1118         is->facecontr= NULL;                            /* to check shared edge */
1119         is->faceisect= is->isect= 0;            /* shared edge, quad half flag */
1120
1121         /* only for shadow! */
1122         if(is->mode==RE_RAY_SHADOW) {
1123         
1124                 /* check with last intersected shadow face */
1125                 if(is->face_last!=NULL && is->face_last!=is->faceorig) {
1126                         if(oc->checkfunc(is, is->face_last)) {
1127                                 is->face= is->face_last;
1128                                 VECSUB(is->vec, is->end, is->start);
1129                                 if(RE_ray_face_intersection(is, oc->coordsfunc)) return 1;
1130                         }
1131                 }
1132         }
1133         
1134         ldx= is->end[0] - is->start[0];
1135         u1= 0.0f;
1136         u2= 1.0f;
1137
1138         /* clip with octree cube */
1139         if(cliptest(-ldx, is->start[0]-oc->min[0], &u1,&u2)) {
1140                 if(cliptest(ldx, oc->max[0]-is->start[0], &u1,&u2)) {
1141                         ldy= is->end[1] - is->start[1];
1142                         if(cliptest(-ldy, is->start[1]-oc->min[1], &u1,&u2)) {
1143                                 if(cliptest(ldy, oc->max[1]-is->start[1], &u1,&u2)) {
1144                                         ldz= is->end[2] - is->start[2];
1145                                         if(cliptest(-ldz, is->start[2]-oc->min[2], &u1,&u2)) {
1146                                                 if(cliptest(ldz, oc->max[2]-is->start[2], &u1,&u2)) {
1147                                                         c1=1;
1148                                                         if(u2<1.0f) {
1149                                                                 is->end[0]= is->start[0]+u2*ldx;
1150                                                                 is->end[1]= is->start[1]+u2*ldy;
1151                                                                 is->end[2]= is->start[2]+u2*ldz;
1152                                                         }
1153                                                         if(u1>0.0f) {
1154                                                                 is->start[0]+=u1*ldx;
1155                                                                 is->start[1]+=u1*ldy;
1156                                                                 is->start[2]+=u1*ldz;
1157                                                         }
1158                                                 }
1159                                         }
1160                                 }
1161                         }
1162                 }
1163         }
1164
1165         if(c1==0) return 0;
1166
1167         /* reset static variables in ocread */
1168         //ocread(oc, oc->ocres, 0, 0);
1169
1170         /* setup 3dda to traverse octree */
1171         ox1= (is->start[0]-oc->min[0])*oc->ocfacx;
1172         oy1= (is->start[1]-oc->min[1])*oc->ocfacy;
1173         oz1= (is->start[2]-oc->min[2])*oc->ocfacz;
1174         ox2= (is->end[0]-oc->min[0])*oc->ocfacx;
1175         oy2= (is->end[1]-oc->min[1])*oc->ocfacy;
1176         oz2= (is->end[2]-oc->min[2])*oc->ocfacz;
1177
1178         ocx1= (int)ox1;
1179         ocy1= (int)oy1;
1180         ocz1= (int)oz1;
1181         ocx2= (int)ox2;
1182         ocy2= (int)oy2;
1183         ocz2= (int)oz2;
1184
1185         /* for intersection */
1186         VECSUB(is->vec, is->end, is->start);
1187
1188         if(ocx1==ocx2 && ocy1==ocy2 && ocz1==ocz2) {
1189                 no= ocread(oc, ocx1, ocy1, ocz1);
1190                 if(no) {
1191                         /* exact intersection with node */
1192                         vec1[0]= ox1; vec1[1]= oy1; vec1[2]= oz1;
1193                         vec2[0]= ox2; vec2[1]= oy2; vec2[2]= oz2;
1194                         calc_ocval_ray(&ocval, (float)ocx1, (float)ocy1, (float)ocz1, vec1, vec2);
1195                         is->ddalabda= 1.0f;
1196                         if( testnode(oc, is, no, ocval) ) return 1;
1197                 }
1198         }
1199         else {
1200                 //static int coh_ocx1,coh_ocx2,coh_ocy1, coh_ocy2,coh_ocz1,coh_ocz2;
1201                 float dox, doy, doz;
1202                 int eqval;
1203                 
1204                 /* calc labda en ld */
1205                 dox= ox1-ox2;
1206                 doy= oy1-oy2;
1207                 doz= oz1-oz2;
1208
1209                 if(dox<-FLT_EPSILON) {
1210                         ldx= -1.0f/dox;
1211                         labdax= (ocx1-ox1+1.0f)*ldx;
1212                         dx= 1;
1213                 } else if(dox>FLT_EPSILON) {
1214                         ldx= 1.0f/dox;
1215                         labdax= (ox1-ocx1)*ldx;
1216                         dx= -1;
1217                 } else {
1218                         labdax=1.0f;
1219                         ldx=0;
1220                         dx= 0;
1221                 }
1222
1223                 if(doy<-FLT_EPSILON) {
1224                         ldy= -1.0f/doy;
1225                         labday= (ocy1-oy1+1.0f)*ldy;
1226                         dy= 1;
1227                 } else if(doy>FLT_EPSILON) {
1228                         ldy= 1.0f/doy;
1229                         labday= (oy1-ocy1)*ldy;
1230                         dy= -1;
1231                 } else {
1232                         labday=1.0f;
1233                         ldy=0;
1234                         dy= 0;
1235                 }
1236
1237                 if(doz<-FLT_EPSILON) {
1238                         ldz= -1.0f/doz;
1239                         labdaz= (ocz1-oz1+1.0f)*ldz;
1240                         dz= 1;
1241                 } else if(doz>FLT_EPSILON) {
1242                         ldz= 1.0f/doz;
1243                         labdaz= (oz1-ocz1)*ldz;
1244                         dz= -1;
1245                 } else {
1246                         labdaz=1.0f;
1247                         ldz=0;
1248                         dz= 0;
1249                 }
1250                 
1251                 xo=ocx1; yo=ocy1; zo=ocz1;
1252                 labdao= ddalabda= MIN3(labdax,labday,labdaz);
1253                 
1254                 vec2[0]= ox1;
1255                 vec2[1]= oy1;
1256                 vec2[2]= oz1;
1257                 
1258                 /* this loop has been constructed to make sure the first and last node of ray
1259                    are always included, even when ddalabda==1.0f or larger */
1260
1261                 while(TRUE) {
1262
1263                         no= ocread(oc, xo, yo, zo);
1264                         if(no) {
1265                                 
1266                                 /* calculate ray intersection with octree node */
1267                                 VECCOPY(vec1, vec2);
1268                                 // dox,y,z is negative
1269                                 vec2[0]= ox1-ddalabda*dox;
1270                                 vec2[1]= oy1-ddalabda*doy;
1271                                 vec2[2]= oz1-ddalabda*doz;
1272                                 calc_ocval_ray(&ocval, (float)xo, (float)yo, (float)zo, vec1, vec2);
1273                                                            
1274                                 is->ddalabda= ddalabda;
1275                                 if( testnode(oc, is, no, ocval) ) return 1;
1276                         }
1277
1278                         labdao= ddalabda;
1279                         
1280                         /* traversing ocree nodes need careful detection of smallest values, with proper
1281                            exceptions for equal labdas */
1282                         eqval= (labdax==labday);
1283                         if(labday==labdaz) eqval += 2;
1284                         if(labdax==labdaz) eqval += 4;
1285                         
1286                         if(eqval) {     // only 4 cases exist!
1287                                 if(eqval==7) {  // x=y=z
1288                                         xo+=dx; labdax+=ldx;
1289                                         yo+=dy; labday+=ldy;
1290                                         zo+=dz; labdaz+=ldz;
1291                                 }
1292                                 else if(eqval==1) { // x=y 
1293                                         if(labday < labdaz) {
1294                                                 xo+=dx; labdax+=ldx;
1295                                                 yo+=dy; labday+=ldy;
1296                                         }
1297                                         else {
1298                                                 zo+=dz; labdaz+=ldz;
1299                                         }
1300                                 }
1301                                 else if(eqval==2) { // y=z
1302                                         if(labdax < labday) {
1303                                                 xo+=dx; labdax+=ldx;
1304                                         }
1305                                         else {
1306                                                 yo+=dy; labday+=ldy;
1307                                                 zo+=dz; labdaz+=ldz;
1308                                         }
1309                                 }
1310                                 else { // x=z
1311                                         if(labday < labdax) {
1312                                                 yo+=dy; labday+=ldy;
1313                                         }
1314                                         else {
1315                                                 xo+=dx; labdax+=ldx;
1316                                                 zo+=dz; labdaz+=ldz;
1317                                         }
1318                                 }
1319                         }
1320                         else {  // all three different, just three cases exist
1321                                 eqval= (labdax<labday);
1322                                 if(labday<labdaz) eqval += 2;
1323                                 if(labdax<labdaz) eqval += 4;
1324                                 
1325                                 if(eqval==7 || eqval==5) { // x smallest
1326                                         xo+=dx; labdax+=ldx;
1327                                 }
1328                                 else if(eqval==2 || eqval==6) { // y smallest
1329                                         yo+=dy; labday+=ldy;
1330                                 }
1331                                 else { // z smallest
1332                                         zo+=dz; labdaz+=ldz;
1333                                 }
1334                                 
1335                         }
1336
1337                         ddalabda=MIN3(labdax,labday,labdaz);
1338                         if(ddalabda==labdao) break;
1339                         /* to make sure the last node is always checked */
1340                         if(labdao>=1.0f) break;
1341                 }
1342         }
1343         
1344         /* reached end, no intersections found */
1345         is->face_last= NULL;
1346         return 0;
1347 }       
1348
1349 float RE_ray_tree_max_size(RayTree *tree)
1350 {
1351         return ((Octree*)tree)->ocsize;
1352 }
1353