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