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