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