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