Patch #3675 by Ed Halley
[blender.git] / source / blender / render / intern / source / ray.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
29 #include <math.h>
30 #include <string.h>
31 #include <stdlib.h>
32 #include <float.h>
33
34 #include "MEM_guardedalloc.h"
35
36 #include "DNA_material_types.h"
37 #include "DNA_lamp_types.h"
38
39 #include "BKE_utildefines.h"
40 #include "BKE_global.h"
41
42 #include "BLI_arithb.h"
43 #include "BLI_rand.h"
44 #include "BLI_jitter.h"
45
46 #include "PIL_time.h"
47
48 #include "render_types.h"
49 #include "renderpipeline.h"
50 #include "rendercore.h"
51 #include "pixelblending.h"
52 #include "pixelshading.h"
53 #include "texture.h"
54
55 #define DDA_SHADOW 0
56 #define DDA_MIRROR 1
57 #define DDA_SHADOW_TRA 2
58
59 #define RAY_TRA         1
60 #define RAY_TRAFLIP     2
61
62 #define DEPTH_SHADOW_TRA  10
63
64 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
65 /* defined in pipeline.c, is hardcopy of active dynamic allocated Render */
66 /* only to be used here in this file, it's for speed */
67 extern struct Render R;
68 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
69
70 /* ********** structs *************** */
71
72 #define BRANCH_ARRAY 1024
73 #define NODE_ARRAY 4096
74
75 typedef struct Isect {
76         float start[3], vec[3], end[3];         /* start+vec = end, in d3dda */
77         float labda, u, v;
78         struct VlakRen *vlr, *vlrcontr, *vlrorig;
79         short isect, mode;      /* isect: which half of quad, mode: DDA_SHADOW, DDA_MIRROR, DDA_SHADOW_TRA */
80         float ddalabda;
81         float col[4];           /* RGBA for shadow_tra */
82         int lay;                        /* -1 default, set for layer lamps */
83         short vlrisect;         /* flag whether vlrcontr was done or not */
84         /* for optimize, last intersected face */
85         VlakRen *vlr_last;
86 } Isect;
87
88 typedef struct Branch
89 {
90         struct Branch *b[8];
91 } Branch;
92
93 typedef struct OcVal 
94 {
95         short ocx, ocy, ocz;
96 } OcVal;
97
98 typedef struct Node
99 {
100         struct VlakRen *v[8];
101         struct OcVal ov[8];
102         struct Node *next;
103 } Node;
104
105
106 /* ******** globals ***************** */
107
108 /* just for statistics */
109 static int raycount;
110 static int accepted, rejected, coherent_ray;
111
112
113 /* **************** ocval method ******************* */
114 /* within one octree node, a set of 3x15 bits defines a 'boundbox' to OR with */
115
116 #define OCVALRES        15
117 #define BROW16(min, max)      (((max)>=OCVALRES? 0xFFFF: (1<<(max+1))-1) - ((min>0)? ((1<<(min))-1):0) )
118
119 static void calc_ocval_face(float *v1, float *v2, float *v3, float *v4, short x, short y, short z, OcVal *ov)
120 {
121         float min[3], max[3];
122         int ocmin, ocmax;
123         
124         VECCOPY(min, v1);
125         VECCOPY(max, v1);
126         DO_MINMAX(v2, min, max);
127         DO_MINMAX(v3, min, max);
128         if(v4) {
129                 DO_MINMAX(v4, min, max);
130         }
131         
132         ocmin= OCVALRES*(min[0]-x); 
133         ocmax= OCVALRES*(max[0]-x);
134         ov->ocx= BROW16(ocmin, ocmax);
135         
136         ocmin= OCVALRES*(min[1]-y); 
137         ocmax= OCVALRES*(max[1]-y);
138         ov->ocy= BROW16(ocmin, ocmax);
139         
140         ocmin= OCVALRES*(min[2]-z); 
141         ocmax= OCVALRES*(max[2]-z);
142         ov->ocz= BROW16(ocmin, ocmax);
143
144 }
145
146 static void calc_ocval_ray(OcVal *ov, float xo, float yo, float zo, float *vec1, float *vec2)
147 {
148         int ocmin, ocmax;
149         
150         if(vec1[0]<vec2[0]) {
151                 ocmin= OCVALRES*(vec1[0] - xo);
152                 ocmax= OCVALRES*(vec2[0] - xo);
153         } else {
154                 ocmin= OCVALRES*(vec2[0] - xo);
155                 ocmax= OCVALRES*(vec1[0] - xo);
156         }
157         ov->ocx= BROW16(ocmin, ocmax);
158
159         if(vec1[1]<vec2[1]) {
160                 ocmin= OCVALRES*(vec1[1] - yo);
161                 ocmax= OCVALRES*(vec2[1] - yo);
162         } else {
163                 ocmin= OCVALRES*(vec2[1] - yo);
164                 ocmax= OCVALRES*(vec1[1] - yo);
165         }
166         ov->ocy= BROW16(ocmin, ocmax);
167
168         if(vec1[2]<vec2[2]) {
169                 ocmin= OCVALRES*(vec1[2] - zo);
170                 ocmax= OCVALRES*(vec2[2] - zo);
171         } else {
172                 ocmin= OCVALRES*(vec2[2] - zo);
173                 ocmax= OCVALRES*(vec1[2] - zo);
174         }
175         ov->ocz= BROW16(ocmin, ocmax);
176 }
177
178 /* ************* octree ************** */
179
180 static Branch *addbranch(Octree *oc, Branch *br, short ocb)
181 {
182         int index;
183         
184         if(br->b[ocb]) return br->b[ocb];
185         
186         oc->branchcount++;
187         index= oc->branchcount>>12;
188         
189         if(oc->adrbranch[index]==NULL)
190                 oc->adrbranch[index]= MEM_callocN(4096*sizeof(Branch), "new oc branch");
191
192         if(oc->branchcount>= BRANCH_ARRAY*4096) {
193                 printf("error; octree branches full\n");
194                 oc->branchcount=0;
195         }
196         
197         return br->b[ocb]= oc->adrbranch[index]+(oc->branchcount & 4095);
198 }
199
200 static Node *addnode(Octree *oc)
201 {
202         int index;
203         
204         oc->nodecount++;
205         index= oc->nodecount>>12;
206         
207         if(oc->adrnode[index]==NULL)
208                 oc->adrnode[index]= MEM_callocN(4096*sizeof(Node),"addnode");
209
210         if(oc->nodecount> NODE_ARRAY*NODE_ARRAY) {
211                 printf("error; octree nodes full\n");
212                 oc->nodecount=0;
213         }
214         
215         return oc->adrnode[index]+(oc->nodecount & 4095);
216 }
217
218 static int face_in_node(VlakRen *vlr, short x, short y, short z, float rtf[][3])
219 {
220         static float nor[3], d;
221         float fx, fy, fz;
222         
223         // init static vars 
224         if(vlr) {
225                 CalcNormFloat(rtf[0], rtf[1], rtf[2], nor);
226                 d= -nor[0]*rtf[0][0] - nor[1]*rtf[0][1] - nor[2]*rtf[0][2];
227                 return 0;
228         }
229         
230         fx= x;
231         fy= y;
232         fz= z;
233         
234         if((x+0)*nor[0] + (y+0)*nor[1] + (z+0)*nor[2] + d > 0.0) {
235                 if((x+1)*nor[0] + (y+0)*nor[1] + (z+0)*nor[2] + d < 0.0) return 1;
236                 if((x+0)*nor[0] + (y+1)*nor[1] + (z+0)*nor[2] + d < 0.0) return 1;
237                 if((x+1)*nor[0] + (y+1)*nor[1] + (z+0)*nor[2] + d < 0.0) return 1;
238         
239                 if((x+0)*nor[0] + (y+0)*nor[1] + (z+1)*nor[2] + d < 0.0) return 1;
240                 if((x+1)*nor[0] + (y+0)*nor[1] + (z+1)*nor[2] + d < 0.0) return 1;
241                 if((x+0)*nor[0] + (y+1)*nor[1] + (z+1)*nor[2] + d < 0.0) return 1;
242                 if((x+1)*nor[0] + (y+1)*nor[1] + (z+1)*nor[2] + d < 0.0) return 1;
243         }
244         else {
245                 if((x+1)*nor[0] + (y+0)*nor[1] + (z+0)*nor[2] + d > 0.0) return 1;
246                 if((x+0)*nor[0] + (y+1)*nor[1] + (z+0)*nor[2] + d > 0.0) return 1;
247                 if((x+1)*nor[0] + (y+1)*nor[1] + (z+0)*nor[2] + d > 0.0) return 1;
248         
249                 if((x+0)*nor[0] + (y+0)*nor[1] + (z+1)*nor[2] + d > 0.0) return 1;
250                 if((x+1)*nor[0] + (y+0)*nor[1] + (z+1)*nor[2] + d > 0.0) return 1;
251                 if((x+0)*nor[0] + (y+1)*nor[1] + (z+1)*nor[2] + d > 0.0) return 1;
252                 if((x+1)*nor[0] + (y+1)*nor[1] + (z+1)*nor[2] + d > 0.0) return 1;
253         }
254
255         return 0;
256 }
257
258 static void ocwrite(Octree *oc, VlakRen *vlr, short x, short y, short z, float rtf[][3])
259 {
260         Branch *br;
261         Node *no;
262         short a, oc0, oc1, oc2, oc3, oc4, oc5;
263
264         x<<=2;
265         y<<=1;
266
267         br= oc->adrbranch[0];
268
269         if(oc->ocres==512) {
270                 oc0= ((x & 1024)+(y & 512)+(z & 256))>>8;
271                 br= addbranch(oc, br, oc0);
272         }
273         if(oc->ocres>=256) {
274                 oc0= ((x & 512)+(y & 256)+(z & 128))>>7;
275                 br= addbranch(oc, br, oc0);
276         }
277         if(oc->ocres>=128) {
278                 oc0= ((x & 256)+(y & 128)+(z & 64))>>6;
279                 br= addbranch(oc, br, oc0);
280         }
281
282         oc0= ((x & 128)+(y & 64)+(z & 32))>>5;
283         oc1= ((x & 64)+(y & 32)+(z & 16))>>4;
284         oc2= ((x & 32)+(y & 16)+(z & 8))>>3;
285         oc3= ((x & 16)+(y & 8)+(z & 4))>>2;
286         oc4= ((x & 8)+(y & 4)+(z & 2))>>1;
287         oc5= ((x & 4)+(y & 2)+(z & 1));
288
289         br= addbranch(oc, br,oc0);
290         br= addbranch(oc, br,oc1);
291         br= addbranch(oc, br,oc2);
292         br= addbranch(oc, br,oc3);
293         br= addbranch(oc, br,oc4);
294         no= (Node *)br->b[oc5];
295         if(no==NULL) br->b[oc5]= (Branch *)(no= addnode(oc));
296
297         while(no->next) no= no->next;
298
299         a= 0;
300         if(no->v[7]) {          /* node full */
301                 no->next= addnode(oc);
302                 no= no->next;
303         }
304         else {
305                 while(no->v[a]!=NULL) a++;
306         }
307         
308         no->v[a]= vlr;
309         
310         calc_ocval_face(rtf[0], rtf[1], rtf[2], rtf[3], x>>2, y>>1, z, &no->ov[a]);
311
312 }
313
314 static void d2dda(Octree *oc, short b1, short b2, short c1, short c2, char *ocface, short rts[][3], float rtf[][3])
315 {
316         int ocx1,ocx2,ocy1,ocy2;
317         int x,y,dx=0,dy=0;
318         float ox1,ox2,oy1,oy2;
319         float labda,labdao,labdax,labday,ldx,ldy;
320
321         ocx1= rts[b1][c1];
322         ocy1= rts[b1][c2];
323         ocx2= rts[b2][c1];
324         ocy2= rts[b2][c2];
325
326         if(ocx1==ocx2 && ocy1==ocy2) {
327                 ocface[oc->ocres*ocx1+ocy1]= 1;
328                 return;
329         }
330
331         ox1= rtf[b1][c1];
332         oy1= rtf[b1][c2];
333         ox2= rtf[b2][c1];
334         oy2= rtf[b2][c2];
335
336         if(ox1!=ox2) {
337                 if(ox2-ox1>0.0) {
338                         labdax= (ox1-ocx1-1.0)/(ox1-ox2);
339                         ldx= -1.0/(ox1-ox2);
340                         dx= 1;
341                 } else {
342                         labdax= (ox1-ocx1)/(ox1-ox2);
343                         ldx= 1.0/(ox1-ox2);
344                         dx= -1;
345                 }
346         } else {
347                 labdax=1.0;
348                 ldx=0;
349         }
350
351         if(oy1!=oy2) {
352                 if(oy2-oy1>0.0) {
353                         labday= (oy1-ocy1-1.0)/(oy1-oy2);
354                         ldy= -1.0/(oy1-oy2);
355                         dy= 1;
356                 } else {
357                         labday= (oy1-ocy1)/(oy1-oy2);
358                         ldy= 1.0/(oy1-oy2);
359                         dy= -1;
360                 }
361         } else {
362                 labday=1.0;
363                 ldy=0;
364         }
365         
366         x=ocx1; y=ocy1;
367         labda= MIN2(labdax, labday);
368         
369         while(TRUE) {
370                 
371                 if(x<0 || y<0 || x>=oc->ocres || y>=oc->ocres);
372                 else ocface[oc->ocres*x+y]= 1;
373                 
374                 labdao=labda;
375                 if(labdax==labday) {
376                         labdax+=ldx;
377                         x+=dx;
378                         labday+=ldy;
379                         y+=dy;
380                 } else {
381                         if(labdax<labday) {
382                                 labdax+=ldx;
383                                 x+=dx;
384                         } else {
385                                 labday+=ldy;
386                                 y+=dy;
387                         }
388                 }
389                 labda=MIN2(labdax,labday);
390                 if(labda==labdao) break;
391                 if(labda>=1.0) break;
392         }
393         ocface[oc->ocres*ocx2+ocy2]=1;
394 }
395
396 static void filltriangle(Octree *oc, short c1, short c2, char *ocface, short *ocmin)
397 {
398         short *ocmax;
399         int a, x, y, y1, y2;
400
401         ocmax=ocmin+3;
402
403         for(x=ocmin[c1];x<=ocmax[c1];x++) {
404                 a= oc->ocres*x;
405                 for(y=ocmin[c2];y<=ocmax[c2];y++) {
406                         if(ocface[a+y]) {
407                                 y++;
408                                 while(ocface[a+y] && y!=ocmax[c2]) y++;
409                                 for(y1=ocmax[c2];y1>y;y1--) {
410                                         if(ocface[a+y1]) {
411                                                 for(y2=y;y2<=y1;y2++) ocface[a+y2]=1;
412                                                 y1=0;
413                                         }
414                                 }
415                                 y=ocmax[c2];
416                         }
417                 }
418         }
419 }
420
421 void freeoctree(Render *re)
422 {
423         Octree *oc= &re->oc;
424         
425         if(G.f & G_DEBUG) {
426                 printf("branches %d nodes %d\n", oc->branchcount, oc->nodecount);
427                 printf("raycount %d \n", raycount);     
428                 printf("ray coherent %d \n", coherent_ray);
429                 printf("accepted %d rejected %d\n", accepted, rejected);
430         }
431         
432         if(oc->adrbranch) {
433                 int a= 0;
434                 while(oc->adrbranch[a]) {
435                         MEM_freeN(oc->adrbranch[a]);
436                         oc->adrbranch[a]= NULL;
437                         a++;
438                 }
439                 MEM_freeN(oc->adrbranch);
440                 oc->adrbranch= NULL;
441         }
442         oc->branchcount= 0;
443         
444         if(oc->adrnode) {
445                 int a= 0;
446                 while(oc->adrnode[a]) {
447                         MEM_freeN(oc->adrnode[a]);
448                         oc->adrnode[a]= NULL;
449                         a++;
450                 }
451                 MEM_freeN(oc->adrnode);
452                 oc->adrnode= NULL;
453         }
454         oc->nodecount= 0;
455 }
456
457 void makeoctree(Render *re)
458 {
459         Octree *oc;
460         VlakRen *vlr=NULL;
461         VertRen *v1, *v2, *v3, *v4;
462         float ocfac[3], t00, t01, t02;
463         float rtf[4][3];
464         int v;
465         int a, b, c, oc1, oc2, oc3, oc4, x, y, z, ocres2;
466         short rts[4][3], ocmin[6], *ocmax;
467         char *ocface;   // front, top, size view of face, to fill in
468         double lasttime= PIL_check_seconds_timer();
469         
470         oc= &re->oc;
471         oc->adrbranch= MEM_callocN(sizeof(void *)*BRANCH_ARRAY, "octree branches");
472         oc->adrnode= MEM_callocN(sizeof(void *)*NODE_ARRAY, "octree nodes");
473
474         ocmax= ocmin+3;
475         
476         /* only for debug info */
477         raycount=0;
478         accepted= 0;
479         rejected= 0;
480         coherent_ray= 0;
481         
482         /* fill main octree struct */
483         oc->ocres= re->r.ocres;
484         ocres2= oc->ocres*oc->ocres;
485         INIT_MINMAX(oc->min, oc->max);
486         
487         /* first min max octree space */
488         for(v=0;v<re->totvlak;v++) {
489                 if((v & 255)==0) vlr= re->blovl[v>>8];  
490                 else vlr++;
491                 if(vlr->mat->mode & MA_TRACEBLE) {      
492                         if((vlr->mat->mode & MA_WIRE)==0) {     
493                                 
494                                 DO_MINMAX(vlr->v1->co, oc->min, oc->max);
495                                 DO_MINMAX(vlr->v2->co, oc->min, oc->max);
496                                 DO_MINMAX(vlr->v3->co, oc->min, oc->max);
497                                 if(vlr->v4) {
498                                         DO_MINMAX(vlr->v4->co, oc->min, oc->max);
499                                 }
500                         }
501                 }
502         }
503
504         if(oc->min[0] > oc->max[0]) return;     /* empty octree */
505
506         oc->adrbranch[0]=(Branch *)MEM_callocN(4096*sizeof(Branch), "makeoctree");
507         
508         /* the lookup table, per face, for which nodes to fill in */
509         ocface= MEM_callocN( 3*ocres2 + 8, "ocface");
510         memset(ocface, 0, 3*ocres2);
511
512         for(c=0;c<3;c++) {      /* octree enlarge, still needed? */
513                 oc->min[c]-= 0.01;
514                 oc->max[c]+= 0.01;
515         }
516
517         t00= oc->max[0]-oc->min[0];
518         t01= oc->max[1]-oc->min[1];
519         t02= oc->max[2]-oc->min[2];
520         
521         /* this minus 0.1 is old safety... seems to be needed? */
522         oc->ocfacx=ocfac[0]= (oc->ocres-0.1)/t00;
523         oc->ocfacy=ocfac[1]= (oc->ocres-0.1)/t01;
524         oc->ocfacz=ocfac[2]= (oc->ocres-0.1)/t02;
525         
526         oc->ocsize= sqrt(t00*t00+t01*t01+t02*t02);      /* global, max size octree */
527
528         for(v=0; v<re->totvlak; v++) {
529                 if((v & 255)==0) {
530                         double time= PIL_check_seconds_timer();
531
532                         vlr= re->blovl[v>>8];
533                         if(re->test_break())
534                                 break;
535                         if(time-lasttime>1.0) {
536                                 char str[32];
537                                 sprintf(str, "Filling Octree: %d", v);
538                                 re->i.infostr= str;
539                                 re->stats_draw(&re->i);
540                                 lasttime= time;
541                         }
542                 }
543                 else vlr++;
544                 
545                 if(vlr->mat->mode & MA_TRACEBLE) {
546                         if((vlr->mat->mode & MA_WIRE)==0) {     
547                                 
548                                 v1= vlr->v1;
549                                 v2= vlr->v2;
550                                 v3= vlr->v3;
551                                 v4= vlr->v4;
552                                 
553                                 for(c=0;c<3;c++) {
554                                         rtf[0][c]= (v1->co[c]-oc->min[c])*ocfac[c] ;
555                                         rts[0][c]= (short)rtf[0][c];
556                                         rtf[1][c]= (v2->co[c]-oc->min[c])*ocfac[c] ;
557                                         rts[1][c]= (short)rtf[1][c];
558                                         rtf[2][c]= (v3->co[c]-oc->min[c])*ocfac[c] ;
559                                         rts[2][c]= (short)rtf[2][c];
560                                         if(v4) {
561                                                 rtf[3][c]= (v4->co[c]-oc->min[c])*ocfac[c] ;
562                                                 rts[3][c]= (short)rtf[3][c];
563                                         }
564                                 }
565                                 
566                                 for(c=0;c<3;c++) {
567                                         oc1= rts[0][c];
568                                         oc2= rts[1][c];
569                                         oc3= rts[2][c];
570                                         if(v4==NULL) {
571                                                 ocmin[c]= MIN3(oc1,oc2,oc3);
572                                                 ocmax[c]= MAX3(oc1,oc2,oc3);
573                                         }
574                                         else {
575                                                 oc4= rts[3][c];
576                                                 ocmin[c]= MIN4(oc1,oc2,oc3,oc4);
577                                                 ocmax[c]= MAX4(oc1,oc2,oc3,oc4);
578                                         }
579                                         if(ocmax[c]>oc->ocres-1) ocmax[c]=oc->ocres-1;
580                                         if(ocmin[c]<0) ocmin[c]=0;
581                                 }
582                                 
583                                 if(ocmin[0]==ocmax[0] && ocmin[1]==ocmax[1] && ocmin[2]==ocmax[2]) {
584                                         ocwrite(oc, vlr, ocmin[0], ocmin[1], ocmin[2], rtf);
585                                 }
586                                 else {
587                                         
588                                         d2dda(oc, 0,1,0,1,ocface+ocres2,rts,rtf);
589                                         d2dda(oc, 0,1,0,2,ocface,rts,rtf);
590                                         d2dda(oc, 0,1,1,2,ocface+2*ocres2,rts,rtf);
591                                         d2dda(oc, 1,2,0,1,ocface+ocres2,rts,rtf);
592                                         d2dda(oc, 1,2,0,2,ocface,rts,rtf);
593                                         d2dda(oc, 1,2,1,2,ocface+2*ocres2,rts,rtf);
594                                         if(v4==NULL) {
595                                                 d2dda(oc, 2,0,0,1,ocface+ocres2,rts,rtf);
596                                                 d2dda(oc, 2,0,0,2,ocface,rts,rtf);
597                                                 d2dda(oc, 2,0,1,2,ocface+2*ocres2,rts,rtf);
598                                         }
599                                         else {
600                                                 d2dda(oc, 2,3,0,1,ocface+ocres2,rts,rtf);
601                                                 d2dda(oc, 2,3,0,2,ocface,rts,rtf);
602                                                 d2dda(oc, 2,3,1,2,ocface+2*ocres2,rts,rtf);
603                                                 d2dda(oc, 3,0,0,1,ocface+ocres2,rts,rtf);
604                                                 d2dda(oc, 3,0,0,2,ocface,rts,rtf);
605                                                 d2dda(oc, 3,0,1,2,ocface+2*ocres2,rts,rtf);
606                                         }
607                                         /* nothing todo with triangle..., just fills :) */
608                                         filltriangle(oc, 0,1,ocface+ocres2,ocmin);
609                                         filltriangle(oc, 0,2,ocface,ocmin);
610                                         filltriangle(oc, 1,2,ocface+2*ocres2,ocmin);
611                                         
612                                         /* init static vars here */
613                                         face_in_node(vlr, 0,0,0, rtf);
614                                         
615                                         for(x=ocmin[0];x<=ocmax[0];x++) {
616                                                 a= oc->ocres*x;
617                                                 for(y=ocmin[1];y<=ocmax[1];y++) {
618                                                         if(ocface[a+y+ocres2]) {
619                                                                 b= oc->ocres*y+2*ocres2;
620                                                                 for(z=ocmin[2];z<=ocmax[2];z++) {
621                                                                         if(ocface[b+z] && ocface[a+z]) {
622                                                                                 if(face_in_node(NULL, x, y, z, rtf))
623                                                                                         ocwrite(oc, vlr, x,y,z, rtf);
624                                                                         }
625                                                                 }
626                                                         }
627                                                 }
628                                         }
629                                         
630                                         /* same loops to clear octree, doubt it can be done smarter */
631                                         for(x=ocmin[0];x<=ocmax[0];x++) {
632                                                 a= oc->ocres*x;
633                                                 for(y=ocmin[1];y<=ocmax[1];y++) {
634                                                         /* x-y */
635                                                         ocface[a+y+ocres2]= 0;
636
637                                                         b= oc->ocres*y + 2*ocres2;
638                                                         for(z=ocmin[2];z<=ocmax[2];z++) {
639                                                                 /* y-z */
640                                                                 ocface[b+z]= 0;
641                                                                 /* x-z */
642                                                                 ocface[a+z]= 0;
643                                                         }
644                                                 }
645                                         }
646                                 }
647                         }
648                 }
649         }
650         
651         MEM_freeN(ocface);
652         re->i.infostr= NULL;
653         re->stats_draw(&re->i);
654 }
655
656 /* ************ raytracer **************** */
657
658 /* only for self-intersecting test with current render face (where ray left) */
659 static int intersection2(VlakRen *vlr, float r0, float r1, float r2, float rx1, float ry1, float rz1)
660 {
661         VertRen *v1,*v2,*v3,*v4=NULL;
662         float x0,x1,x2,t00,t01,t02,t10,t11,t12,t20,t21,t22;
663         float m0, m1, m2, divdet, det, det1;
664         float u1, v, u2;
665
666         v1= vlr->v1; 
667         v2= vlr->v2; 
668         if(vlr->v4) {
669                 v3= vlr->v4;
670                 v4= vlr->v3;
671         }
672         else v3= vlr->v3;       
673
674         t00= v3->co[0]-v1->co[0];
675         t01= v3->co[1]-v1->co[1];
676         t02= v3->co[2]-v1->co[2];
677         t10= v3->co[0]-v2->co[0];
678         t11= v3->co[1]-v2->co[1];
679         t12= v3->co[2]-v2->co[2];
680         
681         x0= t11*r2-t12*r1;
682         x1= t12*r0-t10*r2;
683         x2= t10*r1-t11*r0;
684
685         divdet= t00*x0+t01*x1+t02*x2;
686
687         m0= rx1-v3->co[0];
688         m1= ry1-v3->co[1];
689         m2= rz1-v3->co[2];
690         det1= m0*x0+m1*x1+m2*x2;
691         
692         if(divdet!=0.0) {
693                 u1= det1/divdet;
694
695                 if(u1<=0.0) {
696                         det= t00*(m1*r2-m2*r1);
697                         det+= t01*(m2*r0-m0*r2);
698                         det+= t02*(m0*r1-m1*r0);
699                         v= det/divdet;
700
701                         if(v<=0.0 && (u1 + v) >= -1.0) {
702                                 return 1;
703                         }
704                 }
705         }
706
707         if(v4) {
708
709                 t20= v3->co[0]-v4->co[0];
710                 t21= v3->co[1]-v4->co[1];
711                 t22= v3->co[2]-v4->co[2];
712
713                 divdet= t20*x0+t21*x1+t22*x2;
714                 if(divdet!=0.0) {
715                         u2= det1/divdet;
716                 
717                         if(u2<=0.0) {
718                                 det= t20*(m1*r2-m2*r1);
719                                 det+= t21*(m2*r0-m0*r2);
720                                 det+= t22*(m0*r1-m1*r0);
721                                 v= det/divdet;
722         
723                                 if(v<=0.0 && (u2 + v) >= -1.0) {
724                                         return 2;
725                                 }
726                         }
727                 }
728         }
729         return 0;
730 }
731
732 #if 0
733 /* ray - line intersection */
734 /* disabled until i got real & fast cylinder checking, this code doesnt work proper
735 for faster strands */
736
737 static int intersection_strand(Isect *is)
738 {
739         float v1[3], v2[3];             /* length of strand */
740         float axis[3], rc[3], nor[3], radline, dist, len;
741         
742         /* radius strand */
743         radline= 0.5f*VecLenf(is->vlr->v1->co, is->vlr->v2->co);
744         
745         VecMidf(v1, is->vlr->v1->co, is->vlr->v2->co);
746         VecMidf(v2, is->vlr->v3->co, is->vlr->v4->co);
747         
748         VECSUB(rc, v1, is->start);      /* vector from base ray to base cylinder */
749         VECSUB(axis, v2, v1);           /* cylinder axis */
750         
751         CROSS(nor, is->vec, axis);
752         len= VecLength(nor);
753         
754         if(len<FLT_EPSILON)
755                 return 0;
756
757         dist= INPR(rc, nor)/len;        /* distance between ray and axis cylinder */
758         
759         if(dist<radline && dist>-radline) {
760                 float dot1, dot2, dot3, rlen, alen, div;
761                 float labda;
762                 
763                 /* calculating the intersection point of shortest distance */
764                 dot1 = INPR(rc, is->vec);
765                 dot2 = INPR(is->vec, axis);
766                 dot3 = INPR(rc, axis);
767                 rlen = INPR(is->vec, is->vec);
768                 alen = INPR(axis, axis);
769                 
770                 div = alen * rlen - dot2 * dot2;
771                 if (ABS(div) < FLT_EPSILON)
772                         return 0;
773                 
774                 labda = (dot1*dot2 - dot3*rlen)/div;
775                 
776                 radline/= sqrt(alen);
777                 
778                 /* labda: where on axis do we have closest intersection? */
779                 if(labda >= -radline && labda <= 1.0f+radline) {
780                         VlakRen *vlr= is->vlrorig;
781                         VertRen *v1= is->vlr->v1, *v2= is->vlr->v2, *v3= is->vlr->v3, *v4= is->vlr->v4;
782                                 /* but we dont do shadows from faces sharing edge */
783                         
784                         if(v1==vlr->v1 || v2==vlr->v1 || v3==vlr->v1 || v4==vlr->v1) return 0;
785                         if(v1==vlr->v2 || v2==vlr->v2 || v3==vlr->v2 || v4==vlr->v2) return 0;
786                         if(v1==vlr->v3 || v2==vlr->v3 || v3==vlr->v3 || v4==vlr->v3) return 0;
787                         if(vlr->v4) {
788                                 if(v1==vlr->v4 || v2==vlr->v4 || v3==vlr->v4 || v4==vlr->v4) return 0;
789                         }
790                         return 1;
791                 }
792         }
793         return 0;
794 }
795 #endif
796
797 /* ray - triangle or quad intersection */
798 static int intersection(Isect *is)
799 {
800         VertRen *v1,*v2,*v3,*v4=NULL;
801         float x0,x1,x2,t00,t01,t02,t10,t11,t12,t20,t21,t22,r0,r1,r2;
802         float m0, m1, m2, divdet, det1;
803         short ok=0;
804
805         /* disabled until i got real & fast cylinder checking, this code doesnt work proper
806            for faster strands */
807 //      if(is->mode==DDA_SHADOW && is->vlr->flag & R_STRAND) 
808 //              return intersection_strand(is);
809         
810         v1= is->vlr->v1; 
811         v2= is->vlr->v2; 
812         if(is->vlr->v4) {
813                 v3= is->vlr->v4;
814                 v4= is->vlr->v3;
815         }
816         else v3= is->vlr->v3;   
817
818         t00= v3->co[0]-v1->co[0];
819         t01= v3->co[1]-v1->co[1];
820         t02= v3->co[2]-v1->co[2];
821         t10= v3->co[0]-v2->co[0];
822         t11= v3->co[1]-v2->co[1];
823         t12= v3->co[2]-v2->co[2];
824         
825         r0= is->vec[0];
826         r1= is->vec[1];
827         r2= is->vec[2];
828         
829         x0= t12*r1-t11*r2;
830         x1= t10*r2-t12*r0;
831         x2= t11*r0-t10*r1;
832
833         divdet= t00*x0+t01*x1+t02*x2;
834
835         m0= is->start[0]-v3->co[0];
836         m1= is->start[1]-v3->co[1];
837         m2= is->start[2]-v3->co[2];
838         det1= m0*x0+m1*x1+m2*x2;
839         
840         if(divdet!=0.0) {
841                 float u;
842
843                 divdet= 1.0/divdet;
844                 u= det1*divdet;
845                 if(u<0.0 && u>-1.0) {
846                         float v, cros0, cros1, cros2;
847                         
848                         cros0= m1*t02-m2*t01;
849                         cros1= m2*t00-m0*t02;
850                         cros2= m0*t01-m1*t00;
851                         v= divdet*(cros0*r0 + cros1*r1 + cros2*r2);
852
853                         if(v<0.0 && (u + v) > -1.0) {
854                                 float labda;
855                                 labda= divdet*(cros0*t10 + cros1*t11 + cros2*t12);
856
857                                 if(labda>0.0 && labda<1.0) {
858                                         is->labda= labda;
859                                         is->u= u; is->v= v;
860                                         ok= 1;
861                                 }
862                         }
863                 }
864         }
865
866         if(ok==0 && v4) {
867
868                 t20= v3->co[0]-v4->co[0];
869                 t21= v3->co[1]-v4->co[1];
870                 t22= v3->co[2]-v4->co[2];
871
872                 divdet= t20*x0+t21*x1+t22*x2;
873                 if(divdet!=0.0) {
874                         float u;
875                         divdet= 1.0/divdet;
876                         u = det1*divdet;
877                         
878                         if(u<0.0 && u>-1.0) {
879                                 float v, cros0, cros1, cros2;
880                                 cros0= m1*t22-m2*t21;
881                                 cros1= m2*t20-m0*t22;
882                                 cros2= m0*t21-m1*t20;
883                                 v= divdet*(cros0*r0 + cros1*r1 + cros2*r2);
884         
885                                 if(v<0.0 && (u + v) > -1.0) {
886                                         float labda;
887                                         labda= divdet*(cros0*t10 + cros1*t11 + cros2*t12);
888                                         
889                                         if(labda>0.0 && labda<1.0) {
890                                                 ok= 2;
891                                                 is->labda= labda;
892                                                 is->u= u; is->v= v;
893                                         }
894                                 }
895                         }
896                 }
897         }
898
899         if(ok) {
900                 is->isect= ok;  // wich half of the quad
901                 
902                 if(is->mode!=DDA_SHADOW) {
903                         /* for mirror & tra-shadow: large faces can be filled in too often, this prevents
904                            a face being detected too soon... */
905                         if(is->labda > is->ddalabda) {
906                                 return 0;
907                         }
908                 }
909                 
910                 /* when a shadow ray leaves a face, it can be little outside the edges of it, causing
911                 intersection to be detected in its neighbour face */
912                 
913                 if(is->vlrcontr && is->vlrisect);       // optimizing, the tests below are not needed
914                 else if(is->labda< .1) {
915                         VlakRen *vlr= is->vlrorig;
916                         short de= 0;
917                         
918                         if(v1==vlr->v1 || v2==vlr->v1 || v3==vlr->v1 || v4==vlr->v1) de++;
919                         if(v1==vlr->v2 || v2==vlr->v2 || v3==vlr->v2 || v4==vlr->v2) de++;
920                         if(v1==vlr->v3 || v2==vlr->v3 || v3==vlr->v3 || v4==vlr->v3) de++;
921                         if(vlr->v4) {
922                                 if(v1==vlr->v4 || v2==vlr->v4 || v3==vlr->v4 || v4==vlr->v4) de++;
923                         }
924                         if(de) {
925                                 
926                                 /* so there's a shared edge or vertex, let's intersect ray with vlr
927                                 itself, if that's true we can safely return 1, otherwise we assume
928                                 the intersection is invalid, 0 */
929                                 
930                                 if(is->vlrcontr==NULL) {
931                                         is->vlrcontr= vlr;
932                                         is->vlrisect= intersection2(vlr, -r0, -r1, -r2, is->start[0], is->start[1], is->start[2]);
933                                 }
934
935                                 if(is->vlrisect) return 1;
936                                 return 0;
937                         }
938                 }
939                 
940                 return 1;
941         }
942
943         return 0;
944 }
945
946 /* check all faces in this node */
947 static int testnode(Isect *is, Node *no, OcVal ocval)
948 {
949         VlakRen *vlr;
950         short nr=0;
951         OcVal *ov;
952         
953         if(is->mode==DDA_SHADOW) {
954                 
955                 vlr= no->v[0];
956                 while(vlr) {
957                 
958                         if(is->vlrorig != vlr) {
959
960                                 if(is->lay & vlr->lay) {
961                                         
962                                         ov= no->ov+nr;
963                                         if( (ov->ocx & ocval.ocx) && (ov->ocy & ocval.ocy) && (ov->ocz & ocval.ocz) ) { 
964                                                 //accepted++;
965                                                 is->vlr= vlr;
966         
967                                                 if(intersection(is)) {
968                                                         is->vlr_last= vlr;
969                                                         return 1;
970                                                 }
971                                         }
972                                         //else rejected++;
973                                 }
974                         }
975                         
976                         nr++;
977                         if(nr==8) {
978                                 no= no->next;
979                                 if(no==0) return 0;
980                                 nr=0;
981                         }
982                         vlr= no->v[nr];
983                 }
984         }
985         else {                  /* else mirror and glass  */
986                 Isect isect;
987                 int found= 0;
988                 
989                 is->labda= 1.0; /* needed? */
990                 isect= *is;             /* copy for sorting */
991                 
992                 vlr= no->v[0];
993                 while(vlr) {
994                         
995                         if(is->vlrorig != vlr) {
996                                 /* I now... cpu cycle waste, might do smarter once */
997                                 if(is->mode==DDA_MIRROR && (vlr->mat->mode & MA_ONLYCAST));
998                                 else {
999                                         ov= no->ov+nr;
1000                                         if( (ov->ocx & ocval.ocx) && (ov->ocy & ocval.ocy) && (ov->ocz & ocval.ocz) ) { 
1001                                                 //accepted++;
1002
1003                                                 isect.vlr= vlr;
1004                                                 if(intersection(&isect)) {
1005                                                         if(isect.labda<is->labda) *is= isect;
1006                                                         found= 1;
1007                                                 }
1008                                         }
1009                                         //else rejected++;
1010                                 }
1011                         }
1012                         
1013                         nr++;
1014                         if(nr==8) {
1015                                 no= no->next;
1016                                 if(no==NULL) break;
1017                                 nr=0;
1018                         }
1019                         vlr= no->v[nr];
1020                 }
1021                 
1022                 return found;
1023         }
1024
1025         return 0;
1026 }
1027
1028 /* find the Node for the octree coord x y z */
1029 static Node *ocread(int x, int y, int z)
1030 {
1031         Branch *br;
1032         int oc1;
1033         
1034         x<<=2;
1035         y<<=1;
1036         
1037         br= R.oc.adrbranch[0];
1038         
1039         if(R.oc.ocres==512) {
1040                 oc1= ((x & 1024)+(y & 512)+(z & 256))>>8;
1041                 br= br->b[oc1];
1042                 if(br==NULL) {
1043                         return NULL;
1044                 }
1045         }
1046         if(R.oc.ocres>=256) {
1047                 oc1= ((x & 512)+(y & 256)+(z & 128))>>7;
1048                 br= br->b[oc1];
1049                 if(br==NULL) {
1050                         return NULL;
1051                 }
1052         }
1053         if(R.oc.ocres>=128) {
1054                 oc1= ((x & 256)+(y & 128)+(z & 64))>>6;
1055                 br= br->b[oc1];
1056                 if(br==NULL) {
1057                         return NULL;
1058                 }
1059         }
1060         
1061         oc1= ((x & 128)+(y & 64)+(z & 32))>>5;
1062         br= br->b[oc1];
1063         if(br) {
1064                 oc1= ((x & 64)+(y & 32)+(z & 16))>>4;
1065                 br= br->b[oc1];
1066                 if(br) {
1067                         oc1= ((x & 32)+(y & 16)+(z & 8))>>3;
1068                         br= br->b[oc1];
1069                         if(br) {
1070                                 oc1= ((x & 16)+(y & 8)+(z & 4))>>2;
1071                                 br= br->b[oc1];
1072                                 if(br) {
1073                                         oc1= ((x & 8)+(y & 4)+(z & 2))>>1;
1074                                         br= br->b[oc1];
1075                                         if(br) {
1076                                                 oc1= ((x & 4)+(y & 2)+(z & 1));
1077                                                 return (Node *)br->b[oc1];
1078                                         }
1079                                 }
1080                         }
1081                 }
1082         }
1083         
1084         return NULL;
1085 }
1086
1087 static int cliptest(float p, float q, float *u1, float *u2)
1088 {
1089         float r;
1090
1091         if(p<0.0) {
1092                 if(q<p) return 0;
1093                 else if(q<0.0) {
1094                         r= q/p;
1095                         if(r>*u2) return 0;
1096                         else if(r>*u1) *u1=r;
1097                 }
1098         }
1099         else {
1100                 if(p>0.0) {
1101                         if(q<0.0) return 0;
1102                         else if(q<p) {
1103                                 r= q/p;
1104                                 if(r<*u1) return 0;
1105                                 else if(r<*u2) *u2=r;
1106                         }
1107                 }
1108                 else if(q<0.0) return 0;
1109         }
1110         return 1;
1111 }
1112
1113 /* extensive coherence checks/storage cancels out the benefit of it, and gives errors... we
1114    need better methods, sample code commented out below (ton) */
1115  
1116 /*
1117
1118 in top: static int coh_nodes[16*16*16][6];
1119 in makeoctree: memset(coh_nodes, 0, sizeof(coh_nodes));
1120  
1121 static void add_coherence_test(int ocx1, int ocx2, int ocy1, int ocy2, int ocz1, int ocz2)
1122 {
1123         short *sp;
1124         
1125         sp= coh_nodes[ (ocx2 & 15) + 16*(ocy2 & 15) + 256*(ocz2 & 15) ];
1126         sp[0]= ocx1; sp[1]= ocy1; sp[2]= ocz1;
1127         sp[3]= ocx2; sp[4]= ocy2; sp[5]= ocz2;
1128         
1129 }
1130
1131 static int do_coherence_test(int ocx1, int ocx2, int ocy1, int ocy2, int ocz1, int ocz2)
1132 {
1133         short *sp;
1134         
1135         sp= coh_nodes[ (ocx2 & 15) + 16*(ocy2 & 15) + 256*(ocz2 & 15) ];
1136         if(sp[0]==ocx1 && sp[1]==ocy1 && sp[2]==ocz1 &&
1137            sp[3]==ocx2 && sp[4]==ocy2 && sp[5]==ocz2) return 1;
1138         return 0;
1139 }
1140
1141 */
1142
1143 /* return 1: found valid intersection */
1144 /* starts with is->vlrorig */
1145 static int d3dda(Isect *is)     
1146 {
1147         Node *no;
1148         OcVal ocval;
1149         float vec1[3], vec2[3];
1150         float u1,u2,ox1,ox2,oy1,oy2,oz1,oz2;
1151         float labdao,labdax,ldx,labday,ldy,labdaz,ldz, ddalabda;
1152         int dx,dy,dz;   
1153         int xo,yo,zo,c1=0;
1154         int ocx1,ocx2,ocy1, ocy2,ocz1,ocz2;
1155         
1156         /* clip with octree */
1157         if(R.oc.branchcount==0) return 0;
1158         
1159         /* do this before intersect calls */
1160         is->vlrcontr= NULL;     /*  to check shared edge */
1161
1162         /* only for shadow! */
1163         if(is->mode==DDA_SHADOW) {
1164         
1165                 /* check with last intersected shadow face */
1166                 if(is->vlr_last!=NULL && is->vlr_last!=is->vlrorig) {
1167                         if(is->lay & is->vlr_last->lay) {
1168                                 is->vlr= is->vlr_last;
1169                                 VECSUB(is->vec, is->end, is->start);
1170                                 if(intersection(is)) return 1;
1171                         }
1172                 }
1173         }
1174         
1175         ldx= is->end[0] - is->start[0];
1176         u1= 0.0;
1177         u2= 1.0;
1178
1179         /* clip with octree cube */
1180         if(cliptest(-ldx, is->start[0]-R.oc.min[0], &u1,&u2)) {
1181                 if(cliptest(ldx, R.oc.max[0]-is->start[0], &u1,&u2)) {
1182                         ldy= is->end[1] - is->start[1];
1183                         if(cliptest(-ldy, is->start[1]-R.oc.min[1], &u1,&u2)) {
1184                                 if(cliptest(ldy, R.oc.max[1]-is->start[1], &u1,&u2)) {
1185                                         ldz= is->end[2] - is->start[2];
1186                                         if(cliptest(-ldz, is->start[2]-R.oc.min[2], &u1,&u2)) {
1187                                                 if(cliptest(ldz, R.oc.max[2]-is->start[2], &u1,&u2)) {
1188                                                         c1=1;
1189                                                         if(u2<1.0) {
1190                                                                 is->end[0]= is->start[0]+u2*ldx;
1191                                                                 is->end[1]= is->start[1]+u2*ldy;
1192                                                                 is->end[2]= is->start[2]+u2*ldz;
1193                                                         }
1194                                                         if(u1>0.0) {
1195                                                                 is->start[0]+=u1*ldx;
1196                                                                 is->start[1]+=u1*ldy;
1197                                                                 is->start[2]+=u1*ldz;
1198                                                         }
1199                                                 }
1200                                         }
1201                                 }
1202                         }
1203                 }
1204         }
1205
1206         if(c1==0) return 0;
1207
1208         /* reset static variables in ocread */
1209         //ocread(R.oc.ocres, 0, 0);
1210
1211         /* setup 3dda to traverse octree */
1212         ox1= (is->start[0]-R.oc.min[0])*R.oc.ocfacx;
1213         oy1= (is->start[1]-R.oc.min[1])*R.oc.ocfacy;
1214         oz1= (is->start[2]-R.oc.min[2])*R.oc.ocfacz;
1215         ox2= (is->end[0]-R.oc.min[0])*R.oc.ocfacx;
1216         oy2= (is->end[1]-R.oc.min[1])*R.oc.ocfacy;
1217         oz2= (is->end[2]-R.oc.min[2])*R.oc.ocfacz;
1218
1219         ocx1= (int)ox1;
1220         ocy1= (int)oy1;
1221         ocz1= (int)oz1;
1222         ocx2= (int)ox2;
1223         ocy2= (int)oy2;
1224         ocz2= (int)oz2;
1225
1226         /* for intersection */
1227         VECSUB(is->vec, is->end, is->start);
1228
1229         if(ocx1==ocx2 && ocy1==ocy2 && ocz1==ocz2) {
1230                 no= ocread(ocx1, ocy1, ocz1);
1231                 if(no) {
1232                         /* exact intersection with node */
1233                         vec1[0]= ox1; vec1[1]= oy1; vec1[2]= oz1;
1234                         vec2[0]= ox2; vec2[1]= oy2; vec2[2]= oz2;
1235                         calc_ocval_ray(&ocval, (float)ocx1, (float)ocy1, (float)ocz1, vec1, vec2);
1236                         is->ddalabda= 1.0;
1237                         if( testnode(is, no, ocval) ) return 1;
1238                 }
1239         }
1240         else {
1241                 //static int coh_ocx1,coh_ocx2,coh_ocy1, coh_ocy2,coh_ocz1,coh_ocz2;
1242                 float dox, doy, doz;
1243                 int eqval;
1244                 
1245                 /* calc labda en ld */
1246                 dox= ox1-ox2;
1247                 doy= oy1-oy2;
1248                 doz= oz1-oz2;
1249
1250                 if(dox<-FLT_EPSILON) {
1251                         ldx= -1.0/dox;
1252                         labdax= (ocx1-ox1+1.0)*ldx;
1253                         dx= 1;
1254                 } else if(dox>FLT_EPSILON) {
1255                         ldx= 1.0/dox;
1256                         labdax= (ox1-ocx1)*ldx;
1257                         dx= -1;
1258                 } else {
1259                         labdax=1.0;
1260                         ldx=0;
1261                         dx= 0;
1262                 }
1263
1264                 if(doy<-FLT_EPSILON) {
1265                         ldy= -1.0/doy;
1266                         labday= (ocy1-oy1+1.0)*ldy;
1267                         dy= 1;
1268                 } else if(doy>FLT_EPSILON) {
1269                         ldy= 1.0/doy;
1270                         labday= (oy1-ocy1)*ldy;
1271                         dy= -1;
1272                 } else {
1273                         labday=1.0;
1274                         ldy=0;
1275                         dy= 0;
1276                 }
1277
1278                 if(doz<-FLT_EPSILON) {
1279                         ldz= -1.0/doz;
1280                         labdaz= (ocz1-oz1+1.0)*ldz;
1281                         dz= 1;
1282                 } else if(doz>FLT_EPSILON) {
1283                         ldz= 1.0/doz;
1284                         labdaz= (oz1-ocz1)*ldz;
1285                         dz= -1;
1286                 } else {
1287                         labdaz=1.0;
1288                         ldz=0;
1289                         dz= 0;
1290                 }
1291                 
1292                 xo=ocx1; yo=ocy1; zo=ocz1;
1293                 labdao= ddalabda= MIN3(labdax,labday,labdaz);
1294                 
1295                 vec2[0]= ox1;
1296                 vec2[1]= oy1;
1297                 vec2[2]= oz1;
1298                 
1299                 /* this loop has been constructed to make sure the first and last node of ray
1300                    are always included, even when ddalabda==1.0 or larger */
1301
1302                 while(TRUE) {
1303
1304                         no= ocread(xo, yo, zo);
1305                         if(no) {
1306                                 
1307                                 /* calculate ray intersection with octree node */
1308                                 VECCOPY(vec1, vec2);
1309                                 // dox,y,z is negative
1310                                 vec2[0]= ox1-ddalabda*dox;
1311                                 vec2[1]= oy1-ddalabda*doy;
1312                                 vec2[2]= oz1-ddalabda*doz;
1313                                 calc_ocval_ray(&ocval, (float)xo, (float)yo, (float)zo, vec1, vec2);
1314                                                            
1315                                 is->ddalabda= ddalabda;
1316                                 if( testnode(is, no, ocval) ) return 1;
1317                         }
1318
1319                         labdao= ddalabda;
1320                         
1321                         /* traversing ocree nodes need careful detection of smallest values, with proper
1322                            exceptions for equal labdas */
1323                         eqval= (labdax==labday);
1324                         if(labday==labdaz) eqval += 2;
1325                         if(labdax==labdaz) eqval += 4;
1326                         
1327                         if(eqval) {     // only 4 cases exist!
1328                                 if(eqval==7) {  // x=y=z
1329                                         xo+=dx; labdax+=ldx;
1330                                         yo+=dy; labday+=ldy;
1331                                         zo+=dz; labdaz+=ldz;
1332                                 }
1333                                 else if(eqval==1) { // x=y 
1334                                         if(labday < labdaz) {
1335                                                 xo+=dx; labdax+=ldx;
1336                                                 yo+=dy; labday+=ldy;
1337                                         }
1338                                         else {
1339                                                 zo+=dz; labdaz+=ldz;
1340                                         }
1341                                 }
1342                                 else if(eqval==2) { // y=z
1343                                         if(labdax < labday) {
1344                                                 xo+=dx; labdax+=ldx;
1345                                         }
1346                                         else {
1347                                                 yo+=dy; labday+=ldy;
1348                                                 zo+=dz; labdaz+=ldz;
1349                                         }
1350                                 }
1351                                 else { // x=z
1352                                         if(labday < labdax) {
1353                                                 yo+=dy; labday+=ldy;
1354                                         }
1355                                         else {
1356                                                 xo+=dx; labdax+=ldx;
1357                                                 zo+=dz; labdaz+=ldz;
1358                                         }
1359                                 }
1360                         }
1361                         else {  // all three different, just three cases exist
1362                                 eqval= (labdax<labday);
1363                                 if(labday<labdaz) eqval += 2;
1364                                 if(labdax<labdaz) eqval += 4;
1365                                 
1366                                 if(eqval==7 || eqval==5) { // x smallest
1367                                         xo+=dx; labdax+=ldx;
1368                                 }
1369                                 else if(eqval==2 || eqval==6) { // y smallest
1370                                         yo+=dy; labday+=ldy;
1371                                 }
1372                                 else { // z smallest
1373                                         zo+=dz; labdaz+=ldz;
1374                                 }
1375                                 
1376                         }
1377
1378                         ddalabda=MIN3(labdax,labday,labdaz);
1379                         if(ddalabda==labdao) break;
1380                         /* to make sure the last node is always checked */
1381                         if(labdao>=1.0) break;
1382                 }
1383         }
1384         
1385         /* reached end, no intersections found */
1386         is->vlr_last= NULL;
1387         return 0;
1388 }               
1389
1390
1391 static void shade_ray(Isect *is, ShadeInput *shi, ShadeResult *shr)
1392 {
1393         VlakRen *vlr= is->vlr;
1394         float l;
1395         int osatex= 0;
1396         
1397         /* set up view vector */
1398         VECCOPY(shi->view, is->vec);
1399
1400         /* render co */
1401         shi->co[0]= is->start[0]+is->labda*(shi->view[0]);
1402         shi->co[1]= is->start[1]+is->labda*(shi->view[1]);
1403         shi->co[2]= is->start[2]+is->labda*(shi->view[2]);
1404         
1405         Normalise(shi->view);
1406
1407         shi->vlr= vlr;
1408         shi->mat= vlr->mat;
1409         memcpy(&shi->r, &shi->mat->r, 23*sizeof(float));        // note, keep this synced with render_types.h
1410         shi->har= shi->mat->har;
1411         
1412         /* face normal, check for flip */
1413         l= vlr->n[0]*shi->view[0]+vlr->n[1]*shi->view[1]+vlr->n[2]*shi->view[2];
1414         if(l<0.0) {     
1415                 shi->facenor[0]= -vlr->n[0];
1416                 shi->facenor[1]= -vlr->n[1];
1417                 shi->facenor[2]= -vlr->n[2];
1418                 // only flip lower 4 bits
1419                 shi->puno= vlr->puno ^ 15;
1420         }
1421         else {
1422                 VECCOPY(shi->facenor, vlr->n);
1423                 shi->puno= vlr->puno;
1424         }
1425         
1426         // Osa structs we leave unchanged now
1427         SWAP(int, osatex, shi->osatex);
1428         
1429         shi->dxco[0]= shi->dxco[1]= shi->dxco[2]= 0.0;
1430         shi->dyco[0]= shi->dyco[1]= shi->dyco[2]= 0.0;
1431         
1432         // but, set Osa stuff to zero where it can confuse texture code
1433         if(shi->mat->texco & (TEXCO_NORM|TEXCO_REFL) ) {
1434                 shi->dxno[0]= shi->dxno[1]= shi->dxno[2]= 0.0;
1435                 shi->dyno[0]= shi->dyno[1]= shi->dyno[2]= 0.0;
1436         }
1437
1438         if(vlr->v4) {
1439                 if(is->isect==2) 
1440                         shade_input_set_coords(shi, is->u, is->v, 2, 1, 3);
1441                 else
1442                         shade_input_set_coords(shi, is->u, is->v, 0, 1, 3);
1443         }
1444         else {
1445                 shade_input_set_coords(shi, is->u, is->v, 0, 1, 2);
1446         }
1447         
1448         // SWAP(int, osatex, shi->osatex);  XXXXX!!!!
1449
1450         if(is->mode==DDA_SHADOW_TRA) shade_color(shi, shr);
1451         else {
1452
1453                 shade_lamp_loop(shi, shr);      
1454
1455                 if(shi->translucency!=0.0) {
1456                         ShadeResult shr_t;
1457                         VecMulf(shi->vn, -1.0);
1458                         VecMulf(shi->facenor, -1.0);
1459                         shade_lamp_loop(shi, &shr_t);
1460                         shr->diff[0]+= shi->translucency*shr_t.diff[0];
1461                         shr->diff[1]+= shi->translucency*shr_t.diff[1];
1462                         shr->diff[2]+= shi->translucency*shr_t.diff[2];
1463                         VecMulf(shi->vn, -1.0);
1464                         VecMulf(shi->facenor, -1.0);
1465                 }
1466         }
1467         
1468         SWAP(int, osatex, shi->osatex);  // XXXXX!!!!
1469
1470 }
1471
1472 static int refraction(float *refract, float *n, float *view, float index)
1473 {
1474         float dot, fac;
1475
1476         VECCOPY(refract, view);
1477         
1478         dot= view[0]*n[0] + view[1]*n[1] + view[2]*n[2];
1479
1480         if(dot>0.0) {
1481                 index = 1.0/index;
1482                 fac= 1.0 - (1.0 - dot*dot)*index*index;
1483                 if(fac<= 0.0) return 0;
1484                 fac= -dot*index + sqrt(fac);
1485         }
1486         else {
1487                 fac= 1.0 - (1.0 - dot*dot)*index*index;
1488                 if(fac<= 0.0) return 0;
1489                 fac= -dot*index - sqrt(fac);
1490         }
1491
1492         refract[0]= index*view[0] + fac*n[0];
1493         refract[1]= index*view[1] + fac*n[1];
1494         refract[2]= index*view[2] + fac*n[2];
1495
1496         return 1;
1497 }
1498
1499 /* orn = original face normal */
1500 static void reflection(float *ref, float *n, float *view, float *orn)
1501 {
1502         float f1;
1503         
1504         f1= -2.0*(n[0]*view[0]+ n[1]*view[1]+ n[2]*view[2]);
1505         
1506         ref[0]= (view[0]+f1*n[0]);
1507         ref[1]= (view[1]+f1*n[1]);
1508         ref[2]= (view[2]+f1*n[2]);
1509
1510         if(orn) {
1511                 /* test phong normals, then we should prevent vector going to the back */
1512                 f1= ref[0]*orn[0]+ ref[1]*orn[1]+ ref[2]*orn[2];
1513                 if(f1>0.0) {
1514                         f1+= .01;
1515                         ref[0]-= f1*orn[0];
1516                         ref[1]-= f1*orn[1];
1517                         ref[2]-= f1*orn[2];
1518                 }
1519         }
1520 }
1521
1522 #if 0
1523 static void color_combine(float *result, float fac1, float fac2, float *col1, float *col2)
1524 {
1525         float col1t[3], col2t[3];
1526         
1527         col1t[0]= sqrt(col1[0]);
1528         col1t[1]= sqrt(col1[1]);
1529         col1t[2]= sqrt(col1[2]);
1530         col2t[0]= sqrt(col2[0]);
1531         col2t[1]= sqrt(col2[1]);
1532         col2t[2]= sqrt(col2[2]);
1533
1534         result[0]= (fac1*col1t[0] + fac2*col2t[0]);
1535         result[0]*= result[0];
1536         result[1]= (fac1*col1t[1] + fac2*col2t[1]);
1537         result[1]*= result[1];
1538         result[2]= (fac1*col1t[2] + fac2*col2t[2]);
1539         result[2]*= result[2];
1540 }
1541 #endif
1542
1543 static float shade_by_transmission(Isect *is, ShadeInput *shi, ShadeResult *shr)
1544 {
1545         float dx, dy, dz, d;
1546
1547         if (0 == (shi->mat->mode & (MA_RAYTRANSP|MA_ZTRA)))
1548                 return -1;
1549            
1550         /* shi.co[] calculated by shade_ray() */
1551         dx= shi->co[0] - is->start[0];
1552         dy= shi->co[1] - is->start[1];
1553         dz= shi->co[2] - is->start[2];
1554         d= sqrt(dx*dx+dy*dy+dz*dz);
1555
1556         shr->alpha *= d;
1557         if (shr->alpha > 1.0) shr->alpha= 1.0;
1558
1559         return d;
1560 }
1561
1562 /* the main recursive tracer itself */
1563 static void traceray(ShadeInput *origshi, short depth, float *start, float *vec, float *col, VlakRen *vlr, int traflag)
1564 {
1565         ShadeInput shi;
1566         ShadeResult shr;
1567         Isect isec;
1568         float f, f1, fr, fg, fb;
1569         float ref[3];
1570         
1571         VECCOPY(isec.start, start);
1572         isec.end[0]= start[0]+R.oc.ocsize*vec[0];
1573         isec.end[1]= start[1]+R.oc.ocsize*vec[1];
1574         isec.end[2]= start[2]+R.oc.ocsize*vec[2];
1575         isec.mode= DDA_MIRROR;
1576         isec.vlrorig= vlr;
1577
1578         if( d3dda(&isec) ) {
1579                 float d= 1.0;
1580                 
1581                 shi.mask= origshi->mask;
1582                 shi.osatex= origshi->osatex;
1583                 shi.depth= 1;   // only now to indicate tracing
1584                 shi.thread= origshi->thread;
1585                 shi.xs= origshi->xs;
1586                 shi.ys= origshi->ys;
1587                 shi.lay= origshi->lay;
1588                 shi.do_preview= 0;
1589                 
1590                 shade_ray(&isec, &shi, &shr);
1591                 if (traflag & RAY_TRA)
1592                         d= shade_by_transmission(&isec, &shi, &shr);
1593                 
1594                 if(depth>0) {
1595
1596                         if(shi.mat->mode & (MA_RAYTRANSP|MA_ZTRA) && shr.alpha < 1.0) {
1597                                 float nf, f, f1, refract[3], tracol[4];
1598                                 
1599                                 tracol[0]= shi.r;
1600                                 tracol[1]= shi.g;
1601                                 tracol[2]= shi.b;
1602                                 tracol[3]= col[3];      // we pass on and accumulate alpha
1603                                 
1604                                 if(shi.mat->mode & MA_RAYTRANSP) {
1605                                         /* odd depths: use normal facing viewer, otherwise flip */
1606                                         if(traflag & RAY_TRAFLIP) {
1607                                                 float norm[3];
1608                                                 norm[0]= - shi.vn[0];
1609                                                 norm[1]= - shi.vn[1];
1610                                                 norm[2]= - shi.vn[2];
1611                                                 if (!refraction(refract, norm, shi.view, shi.ang))
1612                                                         reflection(refract, norm, shi.view, shi.vn);
1613                                         }
1614                                         else {
1615                                                 if (!refraction(refract, shi.vn, shi.view, shi.ang))
1616                                                         reflection(refract, shi.vn, shi.view, shi.vn);
1617                                         }
1618                                         traflag |= RAY_TRA;
1619                                         traceray(origshi, depth-1, shi.co, refract, tracol, shi.vlr, traflag ^ RAY_TRAFLIP);
1620                                 }
1621                                 else
1622                                         traceray(origshi, depth-1, shi.co, shi.view, tracol, shi.vlr, 0);
1623                                 
1624                                 f= shr.alpha; f1= 1.0-f;
1625                                 nf= d * shi.mat->filter;
1626                                 fr= 1.0+ nf*(shi.r-1.0);
1627                                 fg= 1.0+ nf*(shi.g-1.0);
1628                                 fb= 1.0+ nf*(shi.b-1.0);
1629                                 shr.diff[0]= f*shr.diff[0] + f1*fr*tracol[0];
1630                                 shr.diff[1]= f*shr.diff[1] + f1*fg*tracol[1];
1631                                 shr.diff[2]= f*shr.diff[2] + f1*fb*tracol[2];
1632                                 
1633                                 shr.spec[0] *=f;
1634                                 shr.spec[1] *=f;
1635                                 shr.spec[2] *=f;
1636
1637                                 col[3]= f1*tracol[3] + f;
1638                         }
1639                         else 
1640                                 col[3]= 1.0;
1641
1642                         if(shi.mat->mode & MA_RAYMIRROR) {
1643                                 f= shi.ray_mirror;
1644                                 if(f!=0.0) f*= fresnel_fac(shi.view, shi.vn, shi.mat->fresnel_mir_i, shi.mat->fresnel_mir);
1645                         }
1646                         else f= 0.0;
1647                         
1648                         if(f!=0.0) {
1649                                 float mircol[4];
1650                                 
1651                                 reflection(ref, shi.vn, shi.view, NULL);                        
1652                                 traceray(origshi, depth-1, shi.co, ref, mircol, shi.vlr, 0);
1653                         
1654                                 f1= 1.0-f;
1655
1656                                 /* combine */
1657                                 //color_combine(col, f*fr*(1.0-shr.spec[0]), f1, col, shr.diff);
1658                                 //col[0]+= shr.spec[0];
1659                                 //col[1]+= shr.spec[1];
1660                                 //col[2]+= shr.spec[2];
1661                                 
1662                                 fr= shi.mirr;
1663                                 fg= shi.mirg;
1664                                 fb= shi.mirb;
1665                 
1666                                 col[0]= f*fr*(1.0-shr.spec[0])*mircol[0] + f1*shr.diff[0] + shr.spec[0];
1667                                 col[1]= f*fg*(1.0-shr.spec[1])*mircol[1] + f1*shr.diff[1] + shr.spec[1];
1668                                 col[2]= f*fb*(1.0-shr.spec[2])*mircol[2] + f1*shr.diff[2] + shr.spec[2];
1669                         }
1670                         else {
1671                                 col[0]= shr.diff[0] + shr.spec[0];
1672                                 col[1]= shr.diff[1] + shr.spec[1];
1673                                 col[2]= shr.diff[2] + shr.spec[2];
1674                         }
1675                 }
1676                 else {
1677                         col[0]= shr.diff[0] + shr.spec[0];
1678                         col[1]= shr.diff[1] + shr.spec[1];
1679                         col[2]= shr.diff[2] + shr.spec[2];
1680                 }
1681                 
1682         }
1683         else {  /* sky */
1684                 VECCOPY(shi.view, vec);
1685                 Normalise(shi.view);
1686                 
1687                 shadeSkyPixelFloat(col, isec.start, shi.view, NULL);
1688         }
1689 }
1690
1691 /* **************** jitter blocks ********** */
1692
1693 /* calc distributed planar energy */
1694
1695 static void DP_energy(float *table, float *vec, int tot, float xsize, float ysize)
1696 {
1697         int x, y, a;
1698         float *fp, force[3], result[3];
1699         float dx, dy, dist, min;
1700         
1701         min= MIN2(xsize, ysize);
1702         min*= min;
1703         result[0]= result[1]= 0.0;
1704         
1705         for(y= -1; y<2; y++) {
1706                 dy= ysize*y;
1707                 for(x= -1; x<2; x++) {
1708                         dx= xsize*x;
1709                         fp= table;
1710                         for(a=0; a<tot; a++, fp+= 2) {
1711                                 force[0]= vec[0] - fp[0]-dx;
1712                                 force[1]= vec[1] - fp[1]-dy;
1713                                 dist= force[0]*force[0] + force[1]*force[1];
1714                                 if(dist < min && dist>0.0) {
1715                                         result[0]+= force[0]/dist;
1716                                         result[1]+= force[1]/dist;
1717                                 }
1718                         }
1719                 }
1720         }
1721         vec[0] += 0.1*min*result[0]/(float)tot;
1722         vec[1] += 0.1*min*result[1]/(float)tot;
1723         // cyclic clamping
1724         vec[0]= vec[0] - xsize*floor(vec[0]/xsize + 0.5);
1725         vec[1]= vec[1] - ysize*floor(vec[1]/ysize + 0.5);
1726 }
1727
1728 // random offset of 1 in 2
1729 static void jitter_plane_offset(float *jitter1, float *jitter2, int tot, float sizex, float sizey, float ofsx, float ofsy)
1730 {
1731         float dsizex= sizex*ofsx;
1732         float dsizey= sizey*ofsy;
1733         float hsizex= 0.5*sizex, hsizey= 0.5*sizey;
1734         int x;
1735         
1736         for(x=tot; x>0; x--, jitter1+=2, jitter2+=2) {
1737                 jitter2[0]= jitter1[0] + dsizex;
1738                 jitter2[1]= jitter1[1] + dsizey;
1739                 if(jitter2[0] > hsizex) jitter2[0]-= sizex;
1740                 if(jitter2[1] > hsizey) jitter2[1]-= sizey;
1741         }
1742 }
1743
1744 /* called from convertBlenderScene.c */
1745 /* we do this in advance to get consistant random, not alter the render seed, and be threadsafe */
1746 void init_jitter_plane(LampRen *lar)
1747 {
1748         float *fp;
1749         int x, iter=12, tot= lar->ray_totsamp;
1750         
1751         fp=lar->jitter= MEM_mallocN(4*tot*2*sizeof(float), "lamp jitter tab");
1752         
1753         /* set per-lamp fixed seed */
1754         BLI_srandom(tot);
1755         
1756         /* fill table with random locations, area_size large */
1757         for(x=0; x<tot; x++, fp+=2) {
1758                 fp[0]= (BLI_frand()-0.5)*lar->area_size;
1759                 fp[1]= (BLI_frand()-0.5)*lar->area_sizey;
1760         }
1761         
1762         while(iter--) {
1763                 fp= lar->jitter;
1764                 for(x=tot; x>0; x--, fp+=2) {
1765                         DP_energy(lar->jitter, fp, tot, lar->area_size, lar->area_sizey);
1766                 }
1767         }
1768         
1769         /* create the dithered tables */
1770         jitter_plane_offset(lar->jitter, lar->jitter+2*tot, tot, lar->area_size, lar->area_sizey, 0.5, 0.0);
1771         jitter_plane_offset(lar->jitter, lar->jitter+4*tot, tot, lar->area_size, lar->area_sizey, 0.5, 0.5);
1772         jitter_plane_offset(lar->jitter, lar->jitter+6*tot, tot, lar->area_size, lar->area_sizey, 0.0, 0.5);
1773 }
1774
1775 /* table around origin, -0.5*size to 0.5*size */
1776 static float *give_jitter_plane(LampRen *lar, int thread, int xs, int ys)
1777 {
1778         int tot;
1779         
1780         tot= lar->ray_totsamp;
1781                         
1782         if(lar->ray_samp_type & LA_SAMP_JITTER) {
1783                 /* made it threadsafe */
1784                 if(thread & 1) {
1785                         if(lar->xold1!=xs || lar->yold1!=ys) {
1786                                 jitter_plane_offset(lar->jitter, lar->jitter+2*tot, tot, lar->area_size, lar->area_sizey, BLI_thread_frand(1), BLI_thread_frand(1));
1787                                 lar->xold1= xs; lar->yold1= ys;
1788                         }
1789                         return lar->jitter+2*tot;
1790                 }
1791                 else {
1792                         if(lar->xold2!=xs || lar->yold2!=ys) {
1793                                 jitter_plane_offset(lar->jitter, lar->jitter+4*tot, tot, lar->area_size, lar->area_sizey, BLI_thread_frand(0), BLI_thread_frand(0));
1794                                 lar->xold2= xs; lar->yold2= ys;
1795                         }
1796                         return lar->jitter+4*tot;
1797                 }
1798         }
1799         if(lar->ray_samp_type & LA_SAMP_DITHER) {
1800                 return lar->jitter + 2*tot*((xs & 1)+2*(ys & 1));
1801         }
1802         
1803         return lar->jitter;
1804 }
1805
1806
1807 /* ***************** main calls ************** */
1808
1809
1810 /* extern call from render loop */
1811 void ray_trace(ShadeInput *shi, ShadeResult *shr)
1812 {
1813         VlakRen *vlr;
1814         float i, f, f1, fr, fg, fb, vec[3], mircol[4], tracol[4];
1815         int do_tra, do_mir;
1816         
1817         do_tra= ((shi->mat->mode & (MA_RAYTRANSP|MA_ZTRA)) && shr->alpha!=1.0);
1818         do_mir= ((shi->mat->mode & MA_RAYMIRROR) && shi->ray_mirror!=0.0);
1819         vlr= shi->vlr;
1820         
1821         if(do_tra) {
1822                 float refract[3];
1823                 
1824                 tracol[3]= shr->alpha;
1825                 
1826                 if(shi->mat->mode & MA_RAYTRANSP) {
1827                         refraction(refract, shi->vn, shi->view, shi->ang);
1828                         traceray(shi, shi->mat->ray_depth_tra, shi->co, refract, tracol, shi->vlr, RAY_TRA|RAY_TRAFLIP);
1829                 }
1830                 else
1831                         traceray(shi, shi->mat->ray_depth_tra, shi->co, shi->view, tracol, shi->vlr, 0);
1832                 
1833                 f= shr->alpha; f1= 1.0-f;
1834                 fr= 1.0+ shi->mat->filter*(shi->r-1.0);
1835                 fg= 1.0+ shi->mat->filter*(shi->g-1.0);
1836                 fb= 1.0+ shi->mat->filter*(shi->b-1.0);
1837                 shr->diff[0]= f*shr->diff[0] + f1*fr*tracol[0];
1838                 shr->diff[1]= f*shr->diff[1] + f1*fg*tracol[1];
1839                 shr->diff[2]= f*shr->diff[2] + f1*fb*tracol[2];
1840
1841                 shr->alpha= tracol[3];
1842         }
1843         
1844         if(do_mir) {
1845         
1846                 i= shi->ray_mirror*fresnel_fac(shi->view, shi->vn, shi->mat->fresnel_mir_i, shi->mat->fresnel_mir);
1847                 if(i!=0.0) {
1848                         fr= shi->mirr;
1849                         fg= shi->mirg;
1850                         fb= shi->mirb;
1851
1852                         if(vlr->flag & R_SMOOTH) 
1853                                 reflection(vec, shi->vn, shi->view, shi->facenor);
1854                         else
1855                                 reflection(vec, shi->vn, shi->view, NULL);
1856         
1857                         traceray(shi, shi->mat->ray_depth, shi->co, vec, mircol, shi->vlr, 0);
1858                         
1859                         f= i*fr*(1.0-shr->spec[0]);     f1= 1.0-i;
1860                         shr->diff[0]= f*mircol[0] + f1*shr->diff[0];
1861                         
1862                         f= i*fg*(1.0-shr->spec[1]);     f1= 1.0-i;
1863                         shr->diff[1]= f*mircol[1] + f1*shr->diff[1];
1864                         
1865                         f= i*fb*(1.0-shr->spec[2]);     f1= 1.0-i;
1866                         shr->diff[2]= f*mircol[2] + f1*shr->diff[2];
1867                 }
1868         }
1869 }
1870
1871 /* color 'shadfac' passes through 'col' with alpha and filter */
1872 /* filter is only applied on alpha defined transparent part */
1873 static void addAlphaLight(float *shadfac, float *col, float alpha, float filter)
1874 {
1875         float fr, fg, fb;
1876         
1877         fr= 1.0+ filter*(col[0]-1.0);
1878         fg= 1.0+ filter*(col[1]-1.0);
1879         fb= 1.0+ filter*(col[2]-1.0);
1880         
1881         shadfac[0]= alpha*col[0] + fr*(1.0-alpha)*shadfac[0];
1882         shadfac[1]= alpha*col[1] + fg*(1.0-alpha)*shadfac[1];
1883         shadfac[2]= alpha*col[2] + fb*(1.0-alpha)*shadfac[2];
1884         
1885         shadfac[3]= (1.0-alpha)*shadfac[3];
1886 }
1887
1888 static void ray_trace_shadow_tra(Isect *is, int depth, int traflag)
1889 {
1890         /* ray to lamp, find first face that intersects, check alpha properties,
1891            if it has col[3]>0.0  continue. so exit when alpha is full */
1892         ShadeInput shi;
1893         ShadeResult shr;
1894
1895         if( d3dda(is)) {
1896                 float d= 1.0;
1897                 /* we got a face */
1898                 
1899                 shi.mask= 1;
1900                 shi.osatex= 0;
1901                 shi.depth= 1;   // only now to indicate tracing
1902                 
1903                 shade_ray(is, &shi, &shr);
1904                 if (traflag & RAY_TRA)
1905                         d= shade_by_transmission(is, &shi, &shr);
1906                 
1907                 /* mix colors based on shadfac (rgb + amount of light factor) */
1908                 addAlphaLight(is->col, shr.diff, shr.alpha, d*shi.mat->filter);
1909                 
1910                 if(depth>0 && is->col[3]>0.0) {
1911                         
1912                         /* adapt isect struct */
1913                         VECCOPY(is->start, shi.co);
1914                         is->vlrorig= shi.vlr;
1915
1916                         ray_trace_shadow_tra(is, depth-1, traflag | RAY_TRA);
1917                 }
1918         }
1919 }
1920
1921 /* not used, test function for ambient occlusion (yaf: pathlight) */
1922 /* main problem; has to be called within shading loop, giving unwanted recursion */
1923 int ray_trace_shadow_rad(ShadeInput *ship, ShadeResult *shr)
1924 {
1925         static int counter=0, only_one= 0;
1926         extern float hashvectf[];
1927         Isect isec;
1928         ShadeInput shi;
1929         ShadeResult shr_t;
1930         float vec[3], accum[3], div= 0.0;
1931         int a;
1932         
1933         if(only_one) {
1934                 return 0;
1935         }
1936         only_one= 1;
1937         
1938         accum[0]= accum[1]= accum[2]= 0.0;
1939         isec.mode= DDA_MIRROR;
1940         isec.vlrorig= ship->vlr;
1941         
1942         for(a=0; a<8*8; a++) {
1943                 
1944                 counter+=3;
1945                 counter %= 768;
1946                 VECCOPY(vec, hashvectf+counter);
1947                 if(ship->vn[0]*vec[0]+ship->vn[1]*vec[1]+ship->vn[2]*vec[2]>0.0) {
1948                         vec[0]-= vec[0];
1949                         vec[1]-= vec[1];
1950                         vec[2]-= vec[2];
1951                 }
1952                 VECCOPY(isec.start, ship->co);
1953                 isec.end[0]= isec.start[0] + R.oc.ocsize*vec[0];
1954                 isec.end[1]= isec.start[1] + R.oc.ocsize*vec[1];
1955                 isec.end[2]= isec.start[2] + R.oc.ocsize*vec[2];
1956                 
1957                 if( d3dda(&isec)) {
1958                         float fac;
1959                         shade_ray(&isec, &shi, &shr_t);
1960                         fac= isec.labda*isec.labda;
1961                         fac= 1.0;
1962                         accum[0]+= fac*(shr_t.diff[0]+shr_t.spec[0]);
1963                         accum[1]+= fac*(shr_t.diff[1]+shr_t.spec[1]);
1964                         accum[2]+= fac*(shr_t.diff[2]+shr_t.spec[2]);
1965                         div+= fac;
1966                 }
1967                 else div+= 1.0;
1968         }
1969         
1970         if(div!=0.0) {
1971                 shr->diff[0]+= accum[0]/div;
1972                 shr->diff[1]+= accum[1]/div;
1973                 shr->diff[2]+= accum[2]/div;
1974         }
1975         shr->alpha= 1.0;
1976         
1977         only_one= 0;
1978         return 1;
1979 }
1980
1981 /* aolight: function to create random unit sphere vectors for total random sampling */
1982 static void RandomSpherical(float *v)
1983 {
1984         float r;
1985         v[2] = 2.f*BLI_frand()-1.f;
1986         if ((r = 1.f - v[2]*v[2])>0.f) {
1987                 float a = 6.283185307f*BLI_frand();
1988                 r = sqrt(r);
1989                 v[0] = r * cos(a);
1990                 v[1] = r * sin(a);
1991         }
1992         else v[2] = 1.f;
1993 }
1994
1995 /* calc distributed spherical energy */
1996 static void DS_energy(float *sphere, int tot, float *vec)
1997 {
1998         float *fp, fac, force[3], res[3];
1999         int a;
2000         
2001         res[0]= res[1]= res[2]= 0.0;
2002         
2003         for(a=0, fp=sphere; a<tot; a++, fp+=3) {
2004                 VecSubf(force, vec, fp);
2005                 fac= force[0]*force[0] + force[1]*force[1] + force[2]*force[2];
2006                 if(fac!=0.0) {
2007                         fac= 1.0/fac;
2008                         res[0]+= fac*force[0];
2009                         res[1]+= fac*force[1];
2010                         res[2]+= fac*force[2];
2011                 }
2012         }
2013
2014         VecMulf(res, 0.5);
2015         VecAddf(vec, vec, res);
2016         Normalise(vec);
2017         
2018 }
2019
2020 /* called from convertBlenderScene.c */
2021 /* creates an equally distributed spherical sample pattern */
2022 void init_ao_sphere(float *sphere, int tot, int iter)
2023 {
2024         float *fp;
2025         int a;
2026
2027         BLI_srandom(tot);
2028         
2029         /* init */
2030         fp= sphere;
2031         for(a=0; a<tot; a++, fp+= 3) {
2032                 RandomSpherical(fp);
2033         }
2034         
2035         while(iter--) {
2036                 for(a=0, fp= sphere; a<tot; a++, fp+= 3) {
2037                         DS_energy(sphere, tot, fp);
2038                 }
2039         }
2040 }
2041
2042
2043 static float *threadsafe_table_sphere(int test, int thread, int xs, int ys)
2044 {
2045         static float sphere1[2*3*256];
2046         static float sphere2[2*3*256];
2047         static int xs1=-1, xs2=-1, ys1=-1, ys2=-1;
2048         
2049         if(thread & 1) {
2050                 if(xs==xs1 && ys==ys1) return sphere1;
2051                 if(test) return NULL;
2052                 xs1= xs; ys1= ys;
2053                 return sphere1;
2054         }
2055         else  {
2056                 if(xs==xs2 && ys==ys2) return sphere2;
2057                 if(test) return NULL;
2058                 xs2= xs; ys2= ys;
2059                 return sphere2;
2060         }
2061 }
2062
2063 static float *sphere_sampler(int type, int resol, int thread, int xs, int ys)
2064 {
2065         int tot;
2066         float *vec;
2067         
2068         if(resol>16) resol= 16;
2069         
2070         tot= 2*resol*resol;
2071
2072         if (type & WO_AORNDSMP) {
2073                 static float sphere[2*3*256];
2074                 int a;
2075                 
2076                 /* total random sampling. NOT THREADSAFE! (should be removed, is not useful) */
2077                 vec= sphere;
2078                 for (a=0; a<tot; a++, vec+=3) {
2079                         RandomSpherical(vec);
2080                 }
2081                 
2082                 return sphere;
2083         } 
2084         else {
2085                 float *sphere;
2086                 float cosfi, sinfi, cost, sint;
2087                 float ang, *vec1;
2088                 int a;
2089                 
2090                 sphere= threadsafe_table_sphere(1, thread, xs, ys);     // returns table if xs and ys were equal to last call
2091                 if(sphere==NULL) {
2092                         sphere= threadsafe_table_sphere(0, thread, xs, ys);
2093                         
2094                         // random rotation
2095                         ang= BLI_thread_frand(thread);
2096                         sinfi= sin(ang); cosfi= cos(ang);
2097                         ang= BLI_thread_frand(thread);
2098                         sint= sin(ang); cost= cos(ang);
2099                         
2100                         vec= R.wrld.aosphere;
2101                         vec1= sphere;
2102                         for (a=0; a<tot; a++, vec+=3, vec1+=3) {
2103                                 vec1[0]= cost*cosfi*vec[0] - sinfi*vec[1] + sint*cosfi*vec[2];
2104                                 vec1[1]= cost*sinfi*vec[0] + cosfi*vec[1] + sint*sinfi*vec[2];
2105                                 vec1[2]= -sint*vec[0] + cost*vec[2];                    
2106                         }
2107                 }
2108                 return sphere;
2109         }
2110 }
2111
2112
2113 /* extern call from shade_lamp_loop, ambient occlusion calculus */
2114 void ray_ao(ShadeInput *shi, float *shadfac)
2115 {
2116         Isect isec;
2117         float *vec, *nrm, div, bias, sh=0;
2118         float maxdist = R.wrld.aodist;
2119         int j= -1, tot, actual=0, skyadded=0;
2120
2121         isec.vlrorig= shi->vlr;
2122         isec.vlr_last= NULL;
2123         isec.mode= DDA_SHADOW;
2124         isec.lay= -1;
2125
2126         shadfac[0]= shadfac[1]= shadfac[2]= 0.0;
2127
2128         // bias prevents smoothed faces to appear flat
2129         if(shi->vlr->flag & R_SMOOTH) {
2130                 bias= G.scene->world->aobias;
2131                 nrm= shi->vn;
2132         }
2133         else {
2134                 bias= 0.0;
2135                 nrm= shi->facenor;
2136         }
2137
2138         vec= sphere_sampler(R.wrld.aomode, R.wrld.aosamp, shi->thread, shi->xs, shi->ys);
2139         
2140         // warning: since we use full sphere now, and dotproduct is below, we do twice as much
2141         tot= 2*R.wrld.aosamp*R.wrld.aosamp;
2142
2143         while(tot--) {
2144                 
2145                 if ((vec[0]*nrm[0] + vec[1]*nrm[1] + vec[2]*nrm[2]) > bias) {
2146                         // only ao samples for mask
2147                         if(R.r.mode & R_OSA) {
2148                                 j++;
2149                                 if(j==R.osa) j= 0;
2150                                 if(!(shi->mask & (1<<j))) {
2151                                         vec+=3;
2152                                         continue;
2153                                 }
2154                         }
2155                         
2156                         actual++;
2157                         
2158                         /* always set start/end, 3dda clips it */
2159                         VECCOPY(isec.start, shi->co);
2160                         isec.end[0] = shi->co[0] - maxdist*vec[0];
2161                         isec.end[1] = shi->co[1] - maxdist*vec[1];
2162                         isec.end[2] = shi->co[2] - maxdist*vec[2];
2163                         
2164                         /* do the trace */
2165                         if (d3dda(&isec)) {
2166                                 if (R.wrld.aomode & WO_AODIST) sh+= exp(-isec.labda*R.wrld.aodistfac); 
2167                                 else sh+= 1.0;
2168                         }
2169                         else if(R.wrld.aocolor!=WO_AOPLAIN) {
2170                                 float skycol[4];
2171                                 float fac, view[3];
2172                                 
2173                                 view[0]= -vec[0];
2174                                 view[1]= -vec[1];
2175                                 view[2]= -vec[2];
2176                                 Normalise(view);
2177                                 
2178                                 if(R.wrld.aocolor==WO_AOSKYCOL) {
2179                                         fac= 0.5*(1.0+view[0]*R.grvec[0]+ view[1]*R.grvec[1]+ view[2]*R.grvec[2]);
2180                                         shadfac[0]+= (1.0-fac)*R.wrld.horr + fac*R.wrld.zenr;
2181                                         shadfac[1]+= (1.0-fac)*R.wrld.horg + fac*R.wrld.zeng;
2182                                         shadfac[2]+= (1.0-fac)*R.wrld.horb + fac*R.wrld.zenb;
2183                                 }
2184                                 else {
2185                                         shadeSkyPixelFloat(skycol, isec.start, view, NULL);
2186                                         shadfac[0]+= skycol[0];
2187                                         shadfac[1]+= skycol[1];
2188                                         shadfac[2]+= skycol[2];
2189                                 }
2190                                 skyadded++;
2191                         }
2192                 }
2193                 // samples
2194                 vec+= 3;
2195         }
2196         
2197         if(actual==0) shadfac[3]= 1.0;
2198         else shadfac[3] = 1.0 - sh/((float)actual);
2199         
2200         if(R.wrld.aocolor!=WO_AOPLAIN && skyadded) {
2201                 div= shadfac[3]/((float)skyadded);
2202                 
2203                 shadfac[0]*= div;       // average color times distances/hits formula
2204                 shadfac[1]*= div;       // average color times distances/hits formula
2205                 shadfac[2]*= div;       // average color times distances/hits formula
2206         }
2207 }
2208
2209
2210
2211 /* extern call from shade_lamp_loop */
2212 void ray_shadow(ShadeInput *shi, LampRen *lar, float *shadfac)
2213 {
2214         Isect isec;
2215         float lampco[3];
2216
2217         /* setup isec */
2218         if(shi->mat->mode & MA_SHADOW_TRA) isec.mode= DDA_SHADOW_TRA;
2219         else isec.mode= DDA_SHADOW;
2220         
2221         if(lar->mode & LA_LAYER) isec.lay= lar->lay; else isec.lay= -1;
2222
2223         /* only when not mir tracing, first hit optimm */
2224         if(shi->depth==0) isec.vlr_last= lar->vlr_last;
2225         else isec.vlr_last= NULL;
2226         
2227         
2228         if(lar->type==LA_SUN || lar->type==LA_HEMI) {
2229                 lampco[0]= shi->co[0] - R.oc.ocsize*lar->vec[0];
2230                 lampco[1]= shi->co[1] - R.oc.ocsize*lar->vec[1];
2231                 lampco[2]= shi->co[2] - R.oc.ocsize*lar->vec[2];
2232         }
2233         else {
2234                 VECCOPY(lampco, lar->co);
2235         }
2236         
2237         if(lar->ray_totsamp<2) {
2238                 
2239                 isec.vlrorig= shi->vlr;
2240                 shadfac[3]= 1.0; // 1.0=full light
2241                 
2242                 /* set up isec vec */
2243                 VECCOPY(isec.start, shi->co);
2244                 VECCOPY(isec.end, lampco);
2245
2246                 if(isec.mode==DDA_SHADOW_TRA) {
2247                         /* isec.col is like shadfac, so defines amount of light (0.0 is full shadow) */
2248                         isec.col[0]= isec.col[1]= isec.col[2]=  1.0;
2249                         isec.col[3]= 1.0;
2250
2251                         ray_trace_shadow_tra(&isec, DEPTH_SHADOW_TRA, 0);
2252                         QUATCOPY(shadfac, isec.col);
2253                         //printf("shadfac %f %f %f %f\n", shadfac[0], shadfac[1], shadfac[2], shadfac[3]);
2254                 }
2255                 else if( d3dda(&isec)) shadfac[3]= 0.0;
2256         }
2257         else {
2258                 /* area soft shadow */
2259                 float *jitlamp;
2260                 float fac=0.0, div=0.0, vec[3];
2261                 int a, j= -1, mask;
2262                 
2263                 if(isec.mode==DDA_SHADOW_TRA) {
2264                         shadfac[0]= shadfac[1]= shadfac[2]= shadfac[3]= 0.0;
2265                 }
2266                 else shadfac[3]= 1.0;                                                   // 1.0=full light
2267                 
2268                 fac= 0.0;
2269                 jitlamp= give_jitter_plane(lar, shi->thread, shi->xs, shi->ys);
2270
2271                 a= lar->ray_totsamp;
2272                 
2273                 /* this correction to make sure we always take at least 1 sample */
2274                 mask= shi->mask;
2275                 if(a==4) mask |= (mask>>4)|(mask>>8);
2276                 else if(a==9) mask |= (mask>>9);
2277                 
2278                 while(a--) {
2279                         
2280                         if(R.r.mode & R_OSA) {
2281                                 j++;
2282                                 if(j>=R.osa) j= 0;
2283                                 if(!(mask & (1<<j))) {
2284                                         jitlamp+= 2;
2285                                         continue;
2286                                 }
2287                         }
2288                         
2289                         isec.vlrorig= shi->vlr; // ray_trace_shadow_tra changes it
2290                         
2291                         vec[0]= jitlamp[0];
2292                         vec[1]= jitlamp[1];
2293                         vec[2]= 0.0;
2294                         Mat3MulVecfl(lar->mat, vec);
2295                         
2296                         /* set start and end, d3dda clips it */
2297                         VECCOPY(isec.start, shi->co);
2298                         isec.end[0]= lampco[0]+vec[0];
2299                         isec.end[1]= lampco[1]+vec[1];
2300                         isec.end[2]= lampco[2]+vec[2];
2301                         
2302                         if(isec.mode==DDA_SHADOW_TRA) {
2303                                 /* isec.col is like shadfac, so defines amount of light (0.0 is full shadow) */
2304                                 isec.col[0]= isec.col[1]= isec.col[2]=  1.0;
2305                                 isec.col[3]= 1.0;
2306                                 
2307                                 ray_trace_shadow_tra(&isec, DEPTH_SHADOW_TRA, 0);
2308                                 shadfac[0] += isec.col[0];
2309                                 shadfac[1] += isec.col[1];
2310                                 shadfac[2] += isec.col[2];
2311                                 shadfac[3] += isec.col[3];
2312                         }
2313                         else if( d3dda(&isec) ) fac+= 1.0;
2314                         
2315                         div+= 1.0;
2316                         jitlamp+= 2;
2317                 }
2318                 
2319                 if(isec.mode==DDA_SHADOW_TRA) {
2320                         shadfac[0] /= div;
2321                         shadfac[1] /= div;
2322                         shadfac[2] /= div;
2323                         shadfac[3] /= div;
2324                 }
2325                 else {
2326                         // sqrt makes nice umbra effect
2327                         if(lar->ray_samp_type & LA_SAMP_UMBRA)
2328                                 shadfac[3]= sqrt(1.0-fac/div);
2329                         else
2330                                 shadfac[3]= 1.0-fac/div;
2331                 }
2332         }
2333
2334         /* for first hit optim, set last interesected shadow face */
2335         if(shi->depth==0) lar->vlr_last= isec.vlr_last;
2336
2337 }
2338
2339 /* only when face points away from lamp, in direction of lamp, trace ray and find first exit point */
2340 void ray_translucent(ShadeInput *shi, LampRen *lar, float *distfac, float *co)
2341 {
2342         Isect isec;
2343         float lampco[3];
2344         
2345         /* setup isec */
2346         isec.mode= DDA_SHADOW_TRA;
2347         
2348         if(lar->mode & LA_LAYER) isec.lay= lar->lay; else isec.lay= -1;
2349         
2350         if(lar->type==LA_SUN || lar->type==LA_HEMI) {
2351                 lampco[0]= shi->co[0] - R.oc.ocsize*lar->vec[0];
2352                 lampco[1]= shi->co[1] - R.oc.ocsize*lar->vec[1];
2353                 lampco[2]= shi->co[2] - R.oc.ocsize*lar->vec[2];
2354         }
2355         else {
2356                 VECCOPY(lampco, lar->co);
2357         }
2358         
2359         isec.vlrorig= shi->vlr;
2360         
2361         /* set up isec vec */
2362         VECCOPY(isec.start, shi->co);
2363         VECCOPY(isec.end, lampco);
2364         
2365         if( d3dda(&isec)) {
2366                 /* we got a face */
2367                 
2368                 /* render co */
2369                 co[0]= isec.start[0]+isec.labda*(isec.vec[0]);
2370                 co[1]= isec.start[1]+isec.labda*(isec.vec[1]);
2371                 co[2]= isec.start[2]+isec.labda*(isec.vec[2]);
2372                 
2373                 *distfac= VecLength(isec.vec);
2374         }
2375         else
2376                 *distfac= 0.0f;
2377 }
2378
2379