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