svn merge -r 22571:22800 https://svn.blender.org/svnroot/bf-blender/trunk/blender
[blender-staging.git] / source / blender / render / intern / source / rayshade.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 #include <math.h>
29 #include <string.h>
30 #include <stdlib.h>
31 #include <float.h>
32 #include <assert.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_blenlib.h"
45 #include "BLI_jitter.h"
46 #include "BLI_rand.h"
47
48 #include "PIL_time.h"
49
50 #include "render_types.h"
51 #include "renderpipeline.h"
52 #include "rendercore.h"
53 #include "renderdatabase.h"
54 #include "pixelblending.h"
55 #include "pixelshading.h"
56 #include "shading.h"
57 #include "texture.h"
58 #include "volumetric.h"
59
60 #include "RE_raytrace.h"
61 #include "rayobject.h"
62 #include "raycounter.h"
63
64 #define RAY_TRA         1
65 #define RAY_TRAFLIP     2
66
67 #define DEPTH_SHADOW_TRA  10
68
69 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
70 /* defined in pipeline.c, is hardcopy of active dynamic allocated Render */
71 /* only to be used here in this file, it's for speed */
72 extern struct Render R;
73 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
74 RayObject *  RE_rayobject_create(int type, int size)
75 {
76         if(type == R_RAYSTRUCTURE_AUTO)
77         {
78                 //TODO
79 //              if(detect_simd())
80 //                      type = R_RAYSTRUCTURE_SIMD_SVBVH;
81 //              else
82 //                      type = R_RAYSTRUCTURE_VBVH;
83
84                         type = R_RAYSTRUCTURE_SIMD_QBVH;
85         }
86         
87         
88         if(type == R_RAYSTRUCTURE_OCTREE)
89         {
90                 //TODO dynamic ocres
91                 return RE_rayobject_octree_create(R.r.ocres, size);
92         }
93         if(type == R_RAYSTRUCTURE_BLIBVH)
94         {
95                 return RE_rayobject_blibvh_create(size);
96         }
97         if(type == R_RAYSTRUCTURE_VBVH)
98         {
99                 return RE_rayobject_vbvh_create(size);
100         }
101         if(type == R_RAYSTRUCTURE_SIMD_SVBVH)
102         {
103                 return RE_rayobject_svbvh_create(size);
104         }
105         if(type == R_RAYSTRUCTURE_SIMD_QBVH)
106         {
107                 return RE_rayobject_qbvh_create(size);
108         }
109         if(type == R_RAYSTRUCTURE_BIH)
110         {
111 //              return RE_rayobject_bih_create(size);
112         }
113         
114         return NULL;
115 }
116
117 #ifdef RE_RAYCOUNTER
118 RayCounter re_rc_counter[BLENDER_MAX_THREADS] = {};
119 #endif
120
121
122 void freeraytree(Render *re)
123 {
124         ObjectInstanceRen *obi;
125         
126         if(re->raytree)
127         {
128                 RE_rayobject_free(re->raytree);
129                 re->raytree = NULL;
130         }
131         if(re->rayfaces)
132         {
133                 MEM_freeN(re->rayfaces);
134                 re->rayfaces = NULL;
135         }
136
137         for(obi=re->instancetable.first; obi; obi=obi->next)
138         {
139                 ObjectRen *obr = obi->obr;
140                 if(obr->raytree)
141                 {
142                         RE_rayobject_free(obr->raytree);
143                         obr->raytree = NULL;
144                 }
145                 if(obr->rayfaces)
146                 {
147                         MEM_freeN(obr->rayfaces);
148                         obr->rayfaces = NULL;
149                 }
150                 if(obi->raytree)
151                 {
152                         RE_rayobject_free(obi->raytree);
153                         obi->raytree = NULL;
154                 }
155         }
156         
157 #ifdef RE_RAYCOUNTER
158         {
159                 RayCounter sum = {};
160                 int i;
161                 for(i=0; i<BLENDER_MAX_THREADS; i++)
162                         RE_RC_MERGE(&sum, re_rc_counter+i);
163                 RE_RC_INFO(&sum);
164         }
165 #endif
166 }
167
168 static int is_raytraceable_vlr(Render *re, VlakRen *vlr)
169 {
170         if((re->flag & R_BAKE_TRACE) || (vlr->mat->mode & MA_TRACEBLE))
171         if(vlr->mat->material_type != MA_TYPE_WIRE)
172                 return 1;
173         return 0;
174 }
175
176 static int is_raytraceable(Render *re, ObjectInstanceRen *obi)
177 {
178         int v;
179         ObjectRen *obr = obi->obr;
180
181         if(re->excludeob && obr->ob == re->excludeob)
182                 return 0;
183
184         for(v=0;v<obr->totvlak;v++)
185         {
186                 VlakRen *vlr = obr->vlaknodes[v>>8].vlak + (v&255);
187                 if(is_raytraceable_vlr(re, vlr))
188                         return 1;
189         }
190         return 0;
191 }
192
193
194 RayObject* makeraytree_object(Render *re, ObjectInstanceRen *obi)
195 {
196         //TODO
197         // out-of-memory safeproof
198         // break render
199         // update render stats
200         ObjectRen *obr = obi->obr;
201         
202         if(obr->raytree == NULL)
203         {
204                 RayObject *raytree;
205                 RayFace *face;
206                 int v;
207                 
208                 //Count faces
209                 int faces = 0;
210                 for(v=0;v<obr->totvlak;v++)
211                 {
212                         VlakRen *vlr = obr->vlaknodes[v>>8].vlak + (v&255);
213                         if(is_raytraceable_vlr(re, vlr))
214                                 faces++;
215                 }
216                 assert( faces > 0 );
217
218                 //Create Ray cast accelaration structure                
219                 raytree = obr->raytree = RE_rayobject_create( re->r.raytrace_tree_type, faces );
220                 face = obr->rayfaces = (RayFace*)MEM_callocN(faces*sizeof(RayFace), "ObjectRen faces");
221                 obr->rayobi = obi;
222                 
223                 for(v=0;v<obr->totvlak;v++)
224                 {
225                         VlakRen *vlr = obr->vlaknodes[v>>8].vlak + (v&255);
226                         if(is_raytraceable_vlr(re, vlr))
227                         {
228                                 RE_rayface_from_vlak( face, obi, vlr );                         
229                                 RE_rayobject_add( raytree, RE_rayobject_unalignRayFace(face) );
230                                 face++;
231                         }
232                 }
233                 RE_rayobject_done( raytree );
234         }
235
236
237         if(obi->flag & R_TRANSFORMED)
238         {
239                 obi->raytree = RE_rayobject_instance_create( obr->raytree, obi->mat, obi, obi->obr->rayobi );
240         }
241         
242         if(obi->raytree) return obi->raytree;
243         return obi->obr->raytree;
244 }
245
246 static int has_special_rayobject(Render *re, ObjectInstanceRen *obi)
247 {
248         if( (obi->flag & R_TRANSFORMED) )
249         {
250                 ObjectRen *obr = obi->obr;
251                 int v, faces = 0;
252                 
253                 for(v=0;v<obr->totvlak;v++)
254                 {
255                         VlakRen *vlr = obr->vlaknodes[v>>8].vlak + (v&255);
256                         if(is_raytraceable_vlr(re, vlr))
257                         {
258                                 faces++;
259                                 if(faces > 4)
260                                         return 1;
261                         }
262                 }
263         }
264         return 0;
265 }
266 /*
267  * create a single raytrace structure with all faces
268  */
269 static void makeraytree_single(Render *re)
270 {
271         ObjectInstanceRen *obi;
272         RayObject *raytree;
273         RayFace *face;
274         int faces = 0, obs = 0;
275
276         for(obi=re->instancetable.first; obi; obi=obi->next)
277         if(is_raytraceable(re, obi))
278         {
279                 int v;
280                 ObjectRen *obr = obi->obr;
281                 obs++;
282                 
283                 if(has_special_rayobject(re, obi))
284                 {
285                         faces++;
286                 }
287                 else
288                 {
289                         for(v=0;v<obr->totvlak;v++)
290                         {
291                                 VlakRen *vlr = obr->vlaknodes[v>>8].vlak + (v&255);
292                                 if(is_raytraceable_vlr(re, vlr))
293                                         faces++;
294                         }
295                 }
296         }
297         
298         //Create raytree
299         raytree = re->raytree = RE_rayobject_create( re->r.raytrace_tree_type, faces );
300
301         face    = re->rayfaces  = (RayFace*)MEM_callocN(faces*sizeof(RayFace), "Render ray faces");
302         
303         for(obi=re->instancetable.first; obi; obi=obi->next)
304         if(is_raytraceable(re, obi))
305         {
306                 if(has_special_rayobject(re, obi))
307                 {
308                         RayObject *obj = makeraytree_object(re, obi);
309                         RE_rayobject_add( re->raytree, obj );
310                 }
311                 else
312                 {
313                         int v;
314                         ObjectRen *obr = obi->obr;
315
316                         for(v=0;v<obr->totvlak;v++)
317                         {
318                                 VlakRen *vlr = obr->vlaknodes[v>>8].vlak + (v&255);
319                                 if(is_raytraceable_vlr(re, vlr))
320                                 {
321                                         RE_rayface_from_vlak(face, obi, vlr);
322                                         if((obi->flag & R_TRANSFORMED))
323                                         {
324                                                 Mat4MulVecfl(obi->mat, face->v1);
325                                                 Mat4MulVecfl(obi->mat, face->v2);
326                                                 Mat4MulVecfl(obi->mat, face->v3);
327                                                 if(RE_rayface_isQuad(face))
328                                                         Mat4MulVecfl(obi->mat, face->v4);
329                                         }
330
331                                         RE_rayobject_add( raytree, RE_rayobject_unalignRayFace(face) );
332                                         face++;
333                                 }
334                         }
335                 }
336         }
337         RE_rayobject_done( raytree );   
338 }
339
340 void makeraytree(Render *re)
341 {
342         float min[3], max[3], sub[3];
343         int i;
344         
345         re->i.infostr= "Make raytree";
346         re->stats_draw(re->sdh, &re->i);
347
348         BENCH(makeraytree_single(re), tree_build);
349                 
350         //Calculate raytree max_size
351         //This is ONLY needed to kept a bogus behaviour of SUN and HEMI lights
352         RE_rayobject_merge_bb( re->raytree, min, max );
353         for(i=0; i<3; i++)
354         {
355                 min[i] += 0.01f;
356                 max[i] += 0.01f;
357                 sub[i] = max[i]-min[i];
358         }
359         re->maxdist = sqrt( sub[0]*sub[0] + sub[1]*sub[1] + sub[2]*sub[2] );
360
361         re->i.infostr= "Raytree finished";
362 }
363
364 void shade_ray(Isect *is, ShadeInput *shi, ShadeResult *shr)
365 {
366         ObjectInstanceRen *obi= (ObjectInstanceRen*)is->hit.ob;
367         VlakRen *vlr= (VlakRen*)is->hit.face;
368         int osatex= 0;
369         
370         /* set up view vector */
371         VECCOPY(shi->view, is->vec);
372
373         /* render co */
374         shi->co[0]= is->start[0]+is->labda*(shi->view[0]);
375         shi->co[1]= is->start[1]+is->labda*(shi->view[1]);
376         shi->co[2]= is->start[2]+is->labda*(shi->view[2]);
377         
378         Normalize(shi->view);
379
380         shi->obi= obi;
381         shi->obr= obi->obr;
382         shi->vlr= vlr;
383         shi->mat= vlr->mat;
384         shade_input_init_material(shi);
385         
386         // Osa structs we leave unchanged now
387         SWAP(int, osatex, shi->osatex);
388         
389         shi->dxco[0]= shi->dxco[1]= shi->dxco[2]= 0.0f;
390         shi->dyco[0]= shi->dyco[1]= shi->dyco[2]= 0.0f;
391         
392         // but, set Osa stuff to zero where it can confuse texture code
393         if(shi->mat->texco & (TEXCO_NORM|TEXCO_REFL) ) {
394                 shi->dxno[0]= shi->dxno[1]= shi->dxno[2]= 0.0f;
395                 shi->dyno[0]= shi->dyno[1]= shi->dyno[2]= 0.0f;
396         }
397
398         if(vlr->v4) {
399                 if(is->isect==2) 
400                         shade_input_set_triangle_i(shi, obi, vlr, 2, 1, 3);
401                 else
402                         shade_input_set_triangle_i(shi, obi, vlr, 0, 1, 3);
403         }
404         else {
405                 shade_input_set_triangle_i(shi, obi, vlr, 0, 1, 2);
406         }
407
408         shi->u= is->u;
409         shi->v= is->v;
410         shi->dx_u= shi->dx_v= shi->dy_u= shi->dy_v=  0.0f;
411
412         shade_input_set_normals(shi);
413
414         /* point normals to viewing direction */
415         if(INPR(shi->facenor, shi->view) < 0.0f)
416                 shade_input_flip_normals(shi);
417
418         shade_input_set_shade_texco(shi);
419         if (shi->mat->material_type == MA_TYPE_VOLUME) {
420                 if(ELEM(is->mode, RE_RAY_SHADOW, RE_RAY_SHADOW_TRA)) {
421                         shade_volume_shadow(shi, shr, is);
422                 } else {
423                         shade_volume_outside(shi, shr);
424                 }
425         }
426         else if(is->mode==RE_RAY_SHADOW_TRA) {
427                 /* temp hack to prevent recursion */
428                 if(shi->nodes==0 && shi->mat->nodetree && shi->mat->use_nodes) {
429                         ntreeShaderExecTree(shi->mat->nodetree, shi, shr);
430                         shi->mat= vlr->mat;             /* shi->mat is being set in nodetree */
431                 }
432                 else
433                         shade_color(shi, shr);
434         }
435         else {
436                 if(shi->mat->nodetree && shi->mat->use_nodes) {
437                         ntreeShaderExecTree(shi->mat->nodetree, shi, shr);
438                         shi->mat= vlr->mat;             /* shi->mat is being set in nodetree */
439                 }
440                 else {
441                         int tempdepth;
442                         /* XXX dodgy business here, set ray depth to -1
443                          * to ignore raytrace in shade_material_loop()
444                          * this could really use a refactor --Matt */
445                         if (shi->volume_depth == 0) {
446                                 tempdepth = shi->depth;
447                                 shi->depth = -1;
448                                 shade_material_loop(shi, shr);
449                                 shi->depth = tempdepth;
450                         } else {
451                                 shade_material_loop(shi, shr);
452                         }
453                 }
454                 /* raytrace likes to separate the spec color */
455                 VECSUB(shr->diff, shr->combined, shr->spec);
456         }       
457         
458         SWAP(int, osatex, shi->osatex);  // XXXXX!!!!
459
460 }
461
462 static int refraction(float *refract, float *n, float *view, float index)
463 {
464         float dot, fac;
465
466         VECCOPY(refract, view);
467         
468         dot= view[0]*n[0] + view[1]*n[1] + view[2]*n[2];
469
470         if(dot>0.0f) {
471                 index = 1.0f/index;
472                 fac= 1.0f - (1.0f - dot*dot)*index*index;
473                 if(fac<= 0.0f) return 0;
474                 fac= -dot*index + sqrt(fac);
475         }
476         else {
477                 fac= 1.0f - (1.0f - dot*dot)*index*index;
478                 if(fac<= 0.0f) return 0;
479                 fac= -dot*index - sqrt(fac);
480         }
481
482         refract[0]= index*view[0] + fac*n[0];
483         refract[1]= index*view[1] + fac*n[1];
484         refract[2]= index*view[2] + fac*n[2];
485
486         return 1;
487 }
488
489 /* orn = original face normal */
490 static void reflection(float *ref, float *n, float *view, float *orn)
491 {
492         float f1;
493         
494         f1= -2.0f*(n[0]*view[0]+ n[1]*view[1]+ n[2]*view[2]);
495         
496         ref[0]= (view[0]+f1*n[0]);
497         ref[1]= (view[1]+f1*n[1]);
498         ref[2]= (view[2]+f1*n[2]);
499
500         if(orn) {
501                 /* test phong normals, then we should prevent vector going to the back */
502                 f1= ref[0]*orn[0]+ ref[1]*orn[1]+ ref[2]*orn[2];
503                 if(f1>0.0f) {
504                         f1+= .01f;
505                         ref[0]-= f1*orn[0];
506                         ref[1]-= f1*orn[1];
507                         ref[2]-= f1*orn[2];
508                 }
509         }
510 }
511
512 #if 0
513 static void color_combine(float *result, float fac1, float fac2, float *col1, float *col2)
514 {
515         float col1t[3], col2t[3];
516         
517         col1t[0]= sqrt(col1[0]);
518         col1t[1]= sqrt(col1[1]);
519         col1t[2]= sqrt(col1[2]);
520         col2t[0]= sqrt(col2[0]);
521         col2t[1]= sqrt(col2[1]);
522         col2t[2]= sqrt(col2[2]);
523
524         result[0]= (fac1*col1t[0] + fac2*col2t[0]);
525         result[0]*= result[0];
526         result[1]= (fac1*col1t[1] + fac2*col2t[1]);
527         result[1]*= result[1];
528         result[2]= (fac1*col1t[2] + fac2*col2t[2]);
529         result[2]*= result[2];
530 }
531 #endif
532
533 static float shade_by_transmission(Isect *is, ShadeInput *shi, ShadeResult *shr)
534 {
535         float dx, dy, dz, d, p;
536
537         if (0 == (shi->mat->mode & MA_TRANSP))
538                 return -1;
539            
540         if (shi->mat->tx_limit <= 0.0f) {
541                 d= 1.0f;
542         } 
543         else {
544                 /* shi.co[] calculated by shade_ray() */
545                 dx= shi->co[0] - is->start[0];
546                 dy= shi->co[1] - is->start[1];
547                 dz= shi->co[2] - is->start[2];
548                 d= sqrt(dx*dx+dy*dy+dz*dz);
549                 if (d > shi->mat->tx_limit)
550                         d= shi->mat->tx_limit;
551
552                 p = shi->mat->tx_falloff;
553                 if(p < 0.0f) p= 0.0f;
554                 else if (p > 10.0f) p= 10.0f;
555
556                 shr->alpha *= pow(d, p);
557                 if (shr->alpha > 1.0f)
558                         shr->alpha= 1.0f;
559         }
560
561         return d;
562 }
563
564 static void ray_fadeout_endcolor(float *col, ShadeInput *origshi, ShadeInput *shi, ShadeResult *shr, Isect *isec, float *vec)
565 {
566         /* un-intersected rays get either rendered material color or sky color */
567         if (origshi->mat->fadeto_mir == MA_RAYMIR_FADETOMAT) {
568                 VECCOPY(col, shr->combined);
569         } else if (origshi->mat->fadeto_mir == MA_RAYMIR_FADETOSKY) {
570                 VECCOPY(shi->view, vec);
571                 Normalize(shi->view);
572                 
573                 shadeSkyView(col, isec->start, shi->view, NULL, shi->thread);
574                 shadeSunView(col, shi->view);
575         }
576 }
577
578 static void ray_fadeout(Isect *is, ShadeInput *shi, float *col, float *blendcol, float dist_mir)
579 {
580         /* if fading out, linear blend against fade color */
581         float blendfac;
582
583         blendfac = 1.0 - VecLenf(shi->co, is->start)/dist_mir;
584         
585         col[0] = col[0]*blendfac + (1.0 - blendfac)*blendcol[0];
586         col[1] = col[1]*blendfac + (1.0 - blendfac)*blendcol[1];
587         col[2] = col[2]*blendfac + (1.0 - blendfac)*blendcol[2];
588 }
589
590 /* the main recursive tracer itself */
591 static void traceray(ShadeInput *origshi, ShadeResult *origshr, short depth, float *start, float *vec, float *col, ObjectInstanceRen *obi, VlakRen *vlr, int traflag)
592 {
593         ShadeInput shi;
594         ShadeResult shr;
595         Isect isec;
596         float f, f1, fr, fg, fb;
597         float ref[3];
598         float dist_mir = origshi->mat->dist_mir;
599
600         /* Warning, This is not that nice, and possibly a bit slow for every ray,
601         however some variables were not initialized properly in, unless using shade_input_initialize(...), we need to do a memset */
602         memset(&shi, 0, sizeof(ShadeInput)); 
603         /* end warning! - Campbell */
604         
605         VECCOPY(isec.start, start);
606         VECCOPY(isec.vec, vec );
607         isec.labda = dist_mir > 0 ? dist_mir : RE_RAYTRACE_MAXDIST;
608         isec.mode= RE_RAY_MIRROR;
609         isec.skip = RE_SKIP_VLR_NEIGHBOUR | RE_SKIP_VLR_RENDER_CHECK;
610         isec.hint = 0;
611
612         isec.orig.ob   = obi;
613         isec.orig.face = vlr;
614         RE_RC_INIT(isec, shi);
615
616         if(RE_rayobject_raycast(R.raytree, &isec)) {
617                 float d= 1.0f;
618                 
619                 shi.mask= origshi->mask;
620                 shi.osatex= origshi->osatex;
621                 shi.depth= 1;                                   /* only used to indicate tracing */
622                 shi.thread= origshi->thread;
623                 //shi.sample= 0; // memset above, so dont need this
624                 shi.xs= origshi->xs;
625                 shi.ys= origshi->ys;
626                 shi.lay= origshi->lay;
627                 shi.passflag= SCE_PASS_COMBINED; /* result of tracing needs no pass info */
628                 shi.combinedflag= 0xFFFFFF;              /* ray trace does all options */
629                 //shi.do_preview= 0; // memset above, so dont need this
630                 shi.light_override= origshi->light_override;
631                 shi.mat_override= origshi->mat_override;
632                 
633                 memset(&shr, 0, sizeof(ShadeResult));
634                 
635                 shade_ray(&isec, &shi, &shr);
636                 if (traflag & RAY_TRA)
637                         d= shade_by_transmission(&isec, &shi, &shr);
638                 
639                 if(depth>0) {
640
641                         if((shi.mat->mode_l & MA_TRANSP) && shr.alpha < 1.0f) {
642                                 float nf, f, f1, refract[3], tracol[4];
643                                 
644                                 tracol[0]= shi.r;
645                                 tracol[1]= shi.g;
646                                 tracol[2]= shi.b;
647                                 tracol[3]= col[3];      // we pass on and accumulate alpha
648                                 
649                                 if((shi.mat->mode & MA_TRANSP) && (shi.mat->mode & MA_RAYTRANSP)) {
650                                         /* odd depths: use normal facing viewer, otherwise flip */
651                                         if(traflag & RAY_TRAFLIP) {
652                                                 float norm[3];
653                                                 norm[0]= - shi.vn[0];
654                                                 norm[1]= - shi.vn[1];
655                                                 norm[2]= - shi.vn[2];
656                                                 if (!refraction(refract, norm, shi.view, shi.ang))
657                                                         reflection(refract, norm, shi.view, shi.vn);
658                                         }
659                                         else {
660                                                 if (!refraction(refract, shi.vn, shi.view, shi.ang))
661                                                         reflection(refract, shi.vn, shi.view, shi.vn);
662                                         }
663                                         traflag |= RAY_TRA;
664                                         traceray(origshi, origshr, depth-1, shi.co, refract, tracol, shi.obi, shi.vlr, traflag ^ RAY_TRAFLIP);
665                                 }
666                                 else
667                                         traceray(origshi, origshr, depth-1, shi.co, shi.view, tracol, shi.obi, shi.vlr, 0);
668                                 
669                                 f= shr.alpha; f1= 1.0f-f;
670                                 nf= d * shi.mat->filter;
671                                 fr= 1.0f+ nf*(shi.r-1.0f);
672                                 fg= 1.0f+ nf*(shi.g-1.0f);
673                                 fb= 1.0f+ nf*(shi.b-1.0f);
674                                 shr.diff[0]= f*shr.diff[0] + f1*fr*tracol[0];
675                                 shr.diff[1]= f*shr.diff[1] + f1*fg*tracol[1];
676                                 shr.diff[2]= f*shr.diff[2] + f1*fb*tracol[2];
677                                 
678                                 shr.spec[0] *=f;
679                                 shr.spec[1] *=f;
680                                 shr.spec[2] *=f;
681
682                                 col[3]= f1*tracol[3] + f;
683                         }
684                         else 
685                                 col[3]= 1.0f;
686
687                         if(shi.mat->mode_l & MA_RAYMIRROR) {
688                                 f= shi.ray_mirror;
689                                 if(f!=0.0f) f*= fresnel_fac(shi.view, shi.vn, shi.mat->fresnel_mir_i, shi.mat->fresnel_mir);
690                         }
691                         else f= 0.0f;
692                         
693                         if(f!=0.0f) {
694                                 float mircol[4];
695                                 
696                                 reflection(ref, shi.vn, shi.view, NULL);                        
697                                 traceray(origshi, origshr, depth-1, shi.co, ref, mircol, shi.obi, shi.vlr, 0);
698                         
699                                 f1= 1.0f-f;
700
701                                 /* combine */
702                                 //color_combine(col, f*fr*(1.0f-shr.spec[0]), f1, col, shr.diff);
703                                 //col[0]+= shr.spec[0];
704                                 //col[1]+= shr.spec[1];
705                                 //col[2]+= shr.spec[2];
706                                 
707                                 fr= shi.mirr;
708                                 fg= shi.mirg;
709                                 fb= shi.mirb;
710                 
711                                 col[0]= f*fr*(1.0f-shr.spec[0])*mircol[0] + f1*shr.diff[0] + shr.spec[0];
712                                 col[1]= f*fg*(1.0f-shr.spec[1])*mircol[1] + f1*shr.diff[1] + shr.spec[1];
713                                 col[2]= f*fb*(1.0f-shr.spec[2])*mircol[2] + f1*shr.diff[2] + shr.spec[2];
714                         }
715                         else {
716                                 col[0]= shr.diff[0] + shr.spec[0];
717                                 col[1]= shr.diff[1] + shr.spec[1];
718                                 col[2]= shr.diff[2] + shr.spec[2];
719                         }
720                         
721                         if (dist_mir > 0.0) {
722                                 float blendcol[3];
723                                 
724                                 /* max ray distance set, but found an intersection, so fade this color
725                                  * out towards the sky/material color for a smooth transition */
726                                 ray_fadeout_endcolor(blendcol, origshi, &shi, origshr, &isec, vec);
727                                 ray_fadeout(&isec, &shi, col, blendcol, dist_mir);
728                         }
729                 }
730                 else {
731                         col[0]= shr.diff[0] + shr.spec[0];
732                         col[1]= shr.diff[1] + shr.spec[1];
733                         col[2]= shr.diff[2] + shr.spec[2];
734                 }
735                 
736         }
737         else {
738                 ray_fadeout_endcolor(col, origshi, &shi, origshr, &isec, vec);
739         }
740         RE_RC_MERGE(&origshi->raycounter, &shi.raycounter);
741 }
742
743 /* **************** jitter blocks ********** */
744
745 /* calc distributed planar energy */
746
747 static void DP_energy(float *table, float *vec, int tot, float xsize, float ysize)
748 {
749         int x, y, a;
750         float *fp, force[3], result[3];
751         float dx, dy, dist, min;
752         
753         min= MIN2(xsize, ysize);
754         min*= min;
755         result[0]= result[1]= 0.0f;
756         
757         for(y= -1; y<2; y++) {
758                 dy= ysize*y;
759                 for(x= -1; x<2; x++) {
760                         dx= xsize*x;
761                         fp= table;
762                         for(a=0; a<tot; a++, fp+= 2) {
763                                 force[0]= vec[0] - fp[0]-dx;
764                                 force[1]= vec[1] - fp[1]-dy;
765                                 dist= force[0]*force[0] + force[1]*force[1];
766                                 if(dist < min && dist>0.0f) {
767                                         result[0]+= force[0]/dist;
768                                         result[1]+= force[1]/dist;
769                                 }
770                         }
771                 }
772         }
773         vec[0] += 0.1*min*result[0]/(float)tot;
774         vec[1] += 0.1*min*result[1]/(float)tot;
775         // cyclic clamping
776         vec[0]= vec[0] - xsize*floor(vec[0]/xsize + 0.5);
777         vec[1]= vec[1] - ysize*floor(vec[1]/ysize + 0.5);
778 }
779
780 // random offset of 1 in 2
781 static void jitter_plane_offset(float *jitter1, float *jitter2, int tot, float sizex, float sizey, float ofsx, float ofsy)
782 {
783         float dsizex= sizex*ofsx;
784         float dsizey= sizey*ofsy;
785         float hsizex= 0.5*sizex, hsizey= 0.5*sizey;
786         int x;
787         
788         for(x=tot; x>0; x--, jitter1+=2, jitter2+=2) {
789                 jitter2[0]= jitter1[0] + dsizex;
790                 jitter2[1]= jitter1[1] + dsizey;
791                 if(jitter2[0] > hsizex) jitter2[0]-= sizex;
792                 if(jitter2[1] > hsizey) jitter2[1]-= sizey;
793         }
794 }
795
796 /* called from convertBlenderScene.c */
797 /* we do this in advance to get consistant random, not alter the render seed, and be threadsafe */
798 void init_jitter_plane(LampRen *lar)
799 {
800         float *fp;
801         int x, iter=12, tot= lar->ray_totsamp;
802         
803         /* test if already initialized */
804         if(lar->jitter) return;
805         
806         /* at least 4, or max threads+1 tables */
807         if(BLENDER_MAX_THREADS < 4) x= 4;
808         else x= BLENDER_MAX_THREADS+1;
809         fp= lar->jitter= MEM_callocN(x*tot*2*sizeof(float), "lamp jitter tab");
810         
811         /* if 1 sample, we leave table to be zero's */
812         if(tot>1) {
813                 
814                 /* set per-lamp fixed seed */
815                 BLI_srandom(tot);
816                 
817                 /* fill table with random locations, area_size large */
818                 for(x=0; x<tot; x++, fp+=2) {
819                         fp[0]= (BLI_frand()-0.5)*lar->area_size;
820                         fp[1]= (BLI_frand()-0.5)*lar->area_sizey;
821                 }
822                 
823                 while(iter--) {
824                         fp= lar->jitter;
825                         for(x=tot; x>0; x--, fp+=2) {
826                                 DP_energy(lar->jitter, fp, tot, lar->area_size, lar->area_sizey);
827                         }
828                 }
829         }       
830         /* create the dithered tables (could just check lamp type!) */
831         jitter_plane_offset(lar->jitter, lar->jitter+2*tot, tot, lar->area_size, lar->area_sizey, 0.5f, 0.0f);
832         jitter_plane_offset(lar->jitter, lar->jitter+4*tot, tot, lar->area_size, lar->area_sizey, 0.5f, 0.5f);
833         jitter_plane_offset(lar->jitter, lar->jitter+6*tot, tot, lar->area_size, lar->area_sizey, 0.0f, 0.5f);
834 }
835
836 /* table around origin, -0.5*size to 0.5*size */
837 static float *give_jitter_plane(LampRen *lar, int thread, int xs, int ys)
838 {
839         int tot;
840         
841         tot= lar->ray_totsamp;
842                         
843         if(lar->ray_samp_type & LA_SAMP_JITTER) {
844                 /* made it threadsafe */
845                 
846                 if(lar->xold[thread]!=xs || lar->yold[thread]!=ys) {
847                         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));
848                         lar->xold[thread]= xs; 
849                         lar->yold[thread]= ys;
850                 }
851                 return lar->jitter+2*(thread+1)*tot;
852         }
853         if(lar->ray_samp_type & LA_SAMP_DITHER) {
854                 return lar->jitter + 2*tot*((xs & 1)+2*(ys & 1));
855         }
856         
857         return lar->jitter;
858 }
859
860
861 /* **************** QMC sampling *************** */
862
863 static void halton_sample(double *ht_invprimes, double *ht_nums, double *v)
864 {
865         // incremental halton sequence generator, from:
866         // "Instant Radiosity", Keller A.
867         unsigned int i;
868         
869         for (i = 0; i < 2; i++)
870         {
871                 double r = fabs((1.0 - ht_nums[i]) - 1e-10);
872                 
873                 if (ht_invprimes[i] >= r)
874                 {
875                         double lasth;
876                         double h = ht_invprimes[i];
877                         
878                         do {
879                                 lasth = h;
880                                 h *= ht_invprimes[i];
881                         } while (h >= r);
882                         
883                         ht_nums[i] += ((lasth + h) - 1.0);
884                 }
885                 else
886                         ht_nums[i] += ht_invprimes[i];
887                 
888                 v[i] = (float)ht_nums[i];
889         }
890 }
891
892 /* Generate Hammersley points in [0,1)^2
893  * From Lucille renderer */
894 static void hammersley_create(double *out, int n)
895 {
896         double p, t;
897         int k, kk;
898
899         for (k = 0; k < n; k++) {
900                 t = 0;
901                 for (p = 0.5, kk = k; kk; p *= 0.5, kk >>= 1) {
902                         if (kk & 1) {           /* kk mod 2 = 1         */
903                                 t += p;
904                         }
905                 }
906         
907                 out[2 * k + 0] = (double)k / (double)n;
908                 out[2 * k + 1] = t;
909         }
910 }
911
912 static struct QMCSampler *QMC_initSampler(int type, int tot)
913 {       
914         QMCSampler *qsa = MEM_callocN(sizeof(QMCSampler), "qmc sampler");
915         qsa->samp2d = MEM_callocN(2*sizeof(double)*tot, "qmc sample table");
916
917         qsa->tot = tot;
918         qsa->type = type;
919         
920         if (qsa->type==SAMP_TYPE_HAMMERSLEY) 
921                 hammersley_create(qsa->samp2d, qsa->tot);
922                 
923         return qsa;
924 }
925
926 static void QMC_initPixel(QMCSampler *qsa, int thread)
927 {
928         if (qsa->type==SAMP_TYPE_HAMMERSLEY)
929         {
930                 /* hammersley sequence is fixed, already created in QMCSampler init.
931                  * per pixel, gets a random offset. We create separate offsets per thread, for write-safety */
932                 qsa->offs[thread][0] = 0.5 * BLI_thread_frand(thread);
933                 qsa->offs[thread][1] = 0.5 * BLI_thread_frand(thread);
934         }
935         else {  /* SAMP_TYPE_HALTON */
936                 
937                 /* generate a new randomised halton sequence per pixel
938                  * to alleviate qmc artifacts and make it reproducable 
939                  * between threads/frames */
940                 double ht_invprimes[2], ht_nums[2];
941                 double r[2];
942                 int i;
943         
944                 ht_nums[0] = BLI_thread_frand(thread);
945                 ht_nums[1] = BLI_thread_frand(thread);
946                 ht_invprimes[0] = 0.5;
947                 ht_invprimes[1] = 1.0/3.0;
948                 
949                 for (i=0; i< qsa->tot; i++) {
950                         halton_sample(ht_invprimes, ht_nums, r);
951                         qsa->samp2d[2*i+0] = r[0];
952                         qsa->samp2d[2*i+1] = r[1];
953                 }
954         }
955 }
956
957 static void QMC_freeSampler(QMCSampler *qsa)
958 {
959         MEM_freeN(qsa->samp2d);
960         MEM_freeN(qsa);
961 }
962
963 static void QMC_getSample(double *s, QMCSampler *qsa, int thread, int num)
964 {
965         if (qsa->type == SAMP_TYPE_HAMMERSLEY) {
966                 s[0] = fmod(qsa->samp2d[2*num+0] + qsa->offs[thread][0], 1.0f);
967                 s[1] = fmod(qsa->samp2d[2*num+1] + qsa->offs[thread][1], 1.0f);
968         }
969         else { /* SAMP_TYPE_HALTON */
970                 s[0] = qsa->samp2d[2*num+0];
971                 s[1] = qsa->samp2d[2*num+1];
972         }
973 }
974
975 /* phong weighted disc using 'blur' for exponent, centred on 0,0 */
976 static void QMC_samplePhong(float *vec, QMCSampler *qsa, int thread, int num, float blur)
977 {
978         double s[2];
979         float phi, pz, sqr;
980         
981         QMC_getSample(s, qsa, thread, num);
982
983         phi = s[0]*2*M_PI;
984         pz = pow(s[1], blur);
985         sqr = sqrt(1.0f-pz*pz);
986
987         vec[0] = cos(phi)*sqr;
988         vec[1] = sin(phi)*sqr;
989         vec[2] = 0.0f;
990 }
991
992 /* rect of edge lengths sizex, sizey, centred on 0.0,0.0 i.e. ranging from -sizex/2 to +sizey/2 */
993 static void QMC_sampleRect(float *vec, QMCSampler *qsa, int thread, int num, float sizex, float sizey)
994 {
995         double s[2];
996
997         QMC_getSample(s, qsa, thread, num);
998                 
999         vec[0] = (s[0] - 0.5) * sizex;
1000         vec[1] = (s[1] - 0.5) * sizey;
1001         vec[2] = 0.0f;
1002 }
1003
1004 /* disc of radius 'radius', centred on 0,0 */
1005 static void QMC_sampleDisc(float *vec, QMCSampler *qsa, int thread, int num, float radius)
1006 {
1007         double s[2];
1008         float phi, sqr;
1009         
1010         QMC_getSample(s, qsa, thread, num);
1011         
1012         phi = s[0]*2*M_PI;
1013         sqr = sqrt(s[1]);
1014
1015         vec[0] = cos(phi)*sqr* radius/2.0;
1016         vec[1] = sin(phi)*sqr* radius/2.0;
1017         vec[2] = 0.0f;
1018 }
1019
1020 /* uniform hemisphere sampling */
1021 static void QMC_sampleHemi(float *vec, QMCSampler *qsa, int thread, int num)
1022 {
1023         double s[2];
1024         float phi, sqr;
1025         
1026         QMC_getSample(s, qsa, thread, num);
1027         
1028         phi = s[0]*2.f*M_PI;    
1029         sqr = sqrt(s[1]);
1030
1031         vec[0] = cos(phi)*sqr;
1032         vec[1] = sin(phi)*sqr;
1033         vec[2] = 1.f - s[1]*s[1];
1034 }
1035
1036 #if 0 /* currently not used */
1037 /* cosine weighted hemisphere sampling */
1038 static void QMC_sampleHemiCosine(float *vec, QMCSampler *qsa, int thread, int num)
1039 {
1040         double s[2];
1041         float phi, sqr;
1042         
1043         QMC_getSample(s, qsa, thread, num);
1044         
1045         phi = s[0]*2.f*M_PI;    
1046         sqr = s[1]*sqrt(2-s[1]*s[1]);
1047
1048         vec[0] = cos(phi)*sqr;
1049         vec[1] = sin(phi)*sqr;
1050         vec[2] = 1.f - s[1]*s[1];
1051
1052 }
1053 #endif
1054
1055 /* called from convertBlenderScene.c */
1056 void init_render_qmcsampler(Render *re)
1057 {
1058         re->qmcsamplers= MEM_callocN(sizeof(ListBase)*BLENDER_MAX_THREADS, "QMCListBase");
1059 }
1060
1061 static QMCSampler *get_thread_qmcsampler(Render *re, int thread, int type, int tot)
1062 {
1063         QMCSampler *qsa;
1064
1065         /* create qmc samplers as needed, since recursion makes it hard to
1066          * predict how many are needed */
1067
1068         for(qsa=re->qmcsamplers[thread].first; qsa; qsa=qsa->next) {
1069                 if(qsa->type == type && qsa->tot == tot && !qsa->used) {
1070                         qsa->used= 1;
1071                         return qsa;
1072                 }
1073         }
1074
1075         qsa= QMC_initSampler(type, tot);
1076         qsa->used= 1;
1077         BLI_addtail(&re->qmcsamplers[thread], qsa);
1078
1079         return qsa;
1080 }
1081
1082 static void release_thread_qmcsampler(Render *re, int thread, QMCSampler *qsa)
1083 {
1084         qsa->used= 0;
1085 }
1086
1087 void free_render_qmcsampler(Render *re)
1088 {
1089         QMCSampler *qsa, *next;
1090         int a;
1091
1092         if(re->qmcsamplers) {
1093                 for(a=0; a<BLENDER_MAX_THREADS; a++) {
1094                         for(qsa=re->qmcsamplers[a].first; qsa; qsa=next) {
1095                                 next= qsa->next;
1096                                 QMC_freeSampler(qsa);
1097                         }
1098
1099                         re->qmcsamplers[a].first= re->qmcsamplers[a].last= NULL;
1100                 }
1101
1102                 MEM_freeN(re->qmcsamplers);
1103                 re->qmcsamplers= NULL;
1104         }
1105 }
1106
1107 static int adaptive_sample_variance(int samples, float *col, float *colsq, float thresh)
1108 {
1109         float var[3], mean[3];
1110
1111         /* scale threshold just to give a bit more precision in input rather than dealing with 
1112          * tiny tiny numbers in the UI */
1113         thresh /= 2;
1114         
1115         mean[0] = col[0] / (float)samples;
1116         mean[1] = col[1] / (float)samples;
1117         mean[2] = col[2] / (float)samples;
1118
1119         var[0] = (colsq[0] / (float)samples) - (mean[0]*mean[0]);
1120         var[1] = (colsq[1] / (float)samples) - (mean[1]*mean[1]);
1121         var[2] = (colsq[2] / (float)samples) - (mean[2]*mean[2]);
1122         
1123         if ((var[0] * 0.4 < thresh) && (var[1] * 0.3 < thresh) && (var[2] * 0.6 < thresh))
1124                 return 1;
1125         else
1126                 return 0;
1127 }
1128
1129 static int adaptive_sample_contrast_val(int samples, float prev, float val, float thresh)
1130 {
1131         /* if the last sample's contribution to the total value was below a small threshold
1132          * (i.e. the samples taken are very similar), then taking more samples that are probably 
1133          * going to be the same is wasting effort */
1134         if (fabs( prev/(float)(samples-1) - val/(float)samples ) < thresh) {
1135                 return 1;
1136         } else
1137                 return 0;
1138 }
1139
1140 static float get_avg_speed(ShadeInput *shi)
1141 {
1142         float pre_x, pre_y, post_x, post_y, speedavg;
1143         
1144         pre_x = (shi->winspeed[0] == PASS_VECTOR_MAX)?0.0:shi->winspeed[0];
1145         pre_y = (shi->winspeed[1] == PASS_VECTOR_MAX)?0.0:shi->winspeed[1];
1146         post_x = (shi->winspeed[2] == PASS_VECTOR_MAX)?0.0:shi->winspeed[2];
1147         post_y = (shi->winspeed[3] == PASS_VECTOR_MAX)?0.0:shi->winspeed[3];
1148         
1149         speedavg = (sqrt(pre_x*pre_x + pre_y*pre_y) + sqrt(post_x*post_x + post_y*post_y)) / 2.0;
1150         
1151         return speedavg;
1152 }
1153
1154 /* ***************** main calls ************** */
1155
1156
1157 static void trace_refract(float *col, ShadeInput *shi, ShadeResult *shr)
1158 {
1159         QMCSampler *qsa=NULL;
1160         int samp_type;
1161         
1162         float samp3d[3], orthx[3], orthy[3];
1163         float v_refract[3], v_refract_new[3];
1164         float sampcol[4], colsq[4];
1165         
1166         float blur = pow(1.0 - shi->mat->gloss_tra, 3);
1167         short max_samples = shi->mat->samp_gloss_tra;
1168         float adapt_thresh = shi->mat->adapt_thresh_tra;
1169         
1170         int samples=0;
1171         
1172         colsq[0] = colsq[1] = colsq[2] = 0.0;
1173         col[0] = col[1] = col[2] = 0.0;
1174         col[3]= shr->alpha;
1175         
1176         if (blur > 0.0) {
1177                 if (adapt_thresh != 0.0) samp_type = SAMP_TYPE_HALTON;
1178                 else samp_type = SAMP_TYPE_HAMMERSLEY;
1179                         
1180                 /* all samples are generated per pixel */
1181                 qsa = get_thread_qmcsampler(&R, shi->thread, samp_type, max_samples);
1182                 QMC_initPixel(qsa, shi->thread);
1183         } else 
1184                 max_samples = 1;
1185         
1186
1187         while (samples < max_samples) {         
1188                 refraction(v_refract, shi->vn, shi->view, shi->ang);
1189                 
1190                 if (max_samples > 1) {
1191                         /* get a quasi-random vector from a phong-weighted disc */
1192                         QMC_samplePhong(samp3d, qsa, shi->thread, samples, blur);
1193                                                 
1194                         VecOrthoBasisf(v_refract, orthx, orthy);
1195                         VecMulf(orthx, samp3d[0]);
1196                         VecMulf(orthy, samp3d[1]);
1197                                 
1198                         /* and perturb the refraction vector in it */
1199                         VecAddf(v_refract_new, v_refract, orthx);
1200                         VecAddf(v_refract_new, v_refract_new, orthy);
1201                         
1202                         Normalize(v_refract_new);
1203                 } else {
1204                         /* no blurriness, use the original normal */
1205                         VECCOPY(v_refract_new, v_refract);
1206                 }
1207         
1208                 traceray(shi, shr, shi->mat->ray_depth_tra, shi->co, v_refract_new, sampcol, shi->obi, shi->vlr, RAY_TRA|RAY_TRAFLIP);
1209         
1210                 col[0] += sampcol[0];
1211                 col[1] += sampcol[1];
1212                 col[2] += sampcol[2];
1213                 col[3] += sampcol[3];
1214                 
1215                 /* for variance calc */
1216                 colsq[0] += sampcol[0]*sampcol[0];
1217                 colsq[1] += sampcol[1]*sampcol[1];
1218                 colsq[2] += sampcol[2]*sampcol[2];
1219                 
1220                 samples++;
1221                 
1222                 /* adaptive sampling */
1223                 if (adapt_thresh < 1.0 && samples > max_samples/2) 
1224                 {
1225                         if (adaptive_sample_variance(samples, col, colsq, adapt_thresh))
1226                                 break;
1227                         
1228                         /* if the pixel so far is very dark, we can get away with less samples */
1229                         if ( (col[0] + col[1] + col[2])/3.0/(float)samples < 0.01 )
1230                                 max_samples--;
1231                 }
1232         }
1233         
1234         col[0] /= (float)samples;
1235         col[1] /= (float)samples;
1236         col[2] /= (float)samples;
1237         col[3] /= (float)samples;
1238         
1239         if (qsa)
1240                 release_thread_qmcsampler(&R, shi->thread, qsa);
1241 }
1242
1243 static void trace_reflect(float *col, ShadeInput *shi, ShadeResult *shr, float fresnelfac)
1244 {
1245         QMCSampler *qsa=NULL;
1246         int samp_type;
1247         
1248         float samp3d[3], orthx[3], orthy[3];
1249         float v_nor_new[3], v_reflect[3];
1250         float sampcol[4], colsq[4];
1251                 
1252         float blur = pow(1.0 - shi->mat->gloss_mir, 3);
1253         short max_samples = shi->mat->samp_gloss_mir;
1254         float adapt_thresh = shi->mat->adapt_thresh_mir;
1255         float aniso = 1.0 - shi->mat->aniso_gloss_mir;
1256         
1257         int samples=0;
1258         
1259         col[0] = col[1] = col[2] = 0.0;
1260         colsq[0] = colsq[1] = colsq[2] = 0.0;
1261         
1262         if (blur > 0.0) {
1263                 if (adapt_thresh != 0.0) samp_type = SAMP_TYPE_HALTON;
1264                 else samp_type = SAMP_TYPE_HAMMERSLEY;
1265                         
1266                 /* all samples are generated per pixel */
1267                 qsa = get_thread_qmcsampler(&R, shi->thread, samp_type, max_samples);
1268                 QMC_initPixel(qsa, shi->thread);
1269         } else 
1270                 max_samples = 1;
1271         
1272         while (samples < max_samples) {
1273                                 
1274                 if (max_samples > 1) {
1275                         /* get a quasi-random vector from a phong-weighted disc */
1276                         QMC_samplePhong(samp3d, qsa, shi->thread, samples, blur);
1277
1278                         /* find the normal's perpendicular plane, blurring along tangents
1279                          * if tangent shading enabled */
1280                         if (shi->mat->mode & (MA_TANGENT_V)) {
1281                                 Crossf(orthx, shi->vn, shi->tang);      // bitangent
1282                                 VECCOPY(orthy, shi->tang);
1283                                 VecMulf(orthx, samp3d[0]);
1284                                 VecMulf(orthy, samp3d[1]*aniso);
1285                         } else {
1286                                 VecOrthoBasisf(shi->vn, orthx, orthy);
1287                                 VecMulf(orthx, samp3d[0]);
1288                                 VecMulf(orthy, samp3d[1]);
1289                         }
1290
1291                         /* and perturb the normal in it */
1292                         VecAddf(v_nor_new, shi->vn, orthx);
1293                         VecAddf(v_nor_new, v_nor_new, orthy);
1294                         Normalize(v_nor_new);
1295                 } else {
1296                         /* no blurriness, use the original normal */
1297                         VECCOPY(v_nor_new, shi->vn);
1298                 }
1299                 
1300                 if((shi->vlr->flag & R_SMOOTH)) 
1301                         reflection(v_reflect, v_nor_new, shi->view, shi->facenor);
1302                 else
1303                         reflection(v_reflect, v_nor_new, shi->view, NULL);
1304                 
1305                 traceray(shi, shr, shi->mat->ray_depth, shi->co, v_reflect, sampcol, shi->obi, shi->vlr, 0);
1306
1307                 
1308                 col[0] += sampcol[0];
1309                 col[1] += sampcol[1];
1310                 col[2] += sampcol[2];
1311         
1312                 /* for variance calc */
1313                 colsq[0] += sampcol[0]*sampcol[0];
1314                 colsq[1] += sampcol[1]*sampcol[1];
1315                 colsq[2] += sampcol[2]*sampcol[2];
1316                 
1317                 samples++;
1318
1319                 /* adaptive sampling */
1320                 if (adapt_thresh > 0.0 && samples > max_samples/3) 
1321                 {
1322                         if (adaptive_sample_variance(samples, col, colsq, adapt_thresh))
1323                                 break;
1324                         
1325                         /* if the pixel so far is very dark, we can get away with less samples */
1326                         if ( (col[0] + col[1] + col[2])/3.0/(float)samples < 0.01 )
1327                                 max_samples--;
1328                 
1329                         /* reduce samples when reflection is dim due to low ray mirror blend value or fresnel factor
1330                          * and when reflection is blurry */
1331                         if (fresnelfac < 0.1 * (blur+1)) {
1332                                 max_samples--;
1333                                 
1334                                 /* even more for very dim */
1335                                 if (fresnelfac < 0.05 * (blur+1)) 
1336                                         max_samples--;
1337                         }
1338                 }
1339         }
1340         
1341         col[0] /= (float)samples;
1342         col[1] /= (float)samples;
1343         col[2] /= (float)samples;
1344         
1345         if (qsa)
1346                 release_thread_qmcsampler(&R, shi->thread, qsa);
1347 }
1348
1349 /* extern call from render loop */
1350 void ray_trace(ShadeInput *shi, ShadeResult *shr)
1351 {
1352         float i, f, f1, fr, fg, fb;
1353         float mircol[4], tracol[4];
1354         float diff[3];
1355         int do_tra, do_mir;
1356         
1357         do_tra= ((shi->mat->mode & MA_TRANSP) && (shi->mat->mode & MA_RAYTRANSP) && shr->alpha!=1.0f);
1358         do_mir= ((shi->mat->mode & MA_RAYMIRROR) && shi->ray_mirror!=0.0f);
1359         
1360         /* raytrace mirror amd refract like to separate the spec color */
1361         if(shi->combinedflag & SCE_PASS_SPEC)
1362                 VECSUB(diff, shr->combined, shr->spec) /* no ; */
1363         else
1364                 VECCOPY(diff, shr->combined);
1365         
1366         if(do_tra) {
1367                 float olddiff[3];
1368                 
1369                 trace_refract(tracol, shi, shr);
1370                 
1371                 f= shr->alpha; f1= 1.0f-f;
1372                 fr= 1.0f+ shi->mat->filter*(shi->r-1.0f);
1373                 fg= 1.0f+ shi->mat->filter*(shi->g-1.0f);
1374                 fb= 1.0f+ shi->mat->filter*(shi->b-1.0f);
1375                 
1376                 /* for refract pass */
1377                 VECCOPY(olddiff, diff);
1378                 
1379                 diff[0]= f*diff[0] + f1*fr*tracol[0];
1380                 diff[1]= f*diff[1] + f1*fg*tracol[1];
1381                 diff[2]= f*diff[2] + f1*fb*tracol[2];
1382                 
1383                 if(shi->passflag & SCE_PASS_REFRACT)
1384                         VECSUB(shr->refr, diff, olddiff);
1385                 
1386                 if(!(shi->combinedflag & SCE_PASS_REFRACT))
1387                         VECSUB(diff, diff, shr->refr);
1388                 
1389                 shr->alpha= tracol[3];
1390         }
1391         
1392         if(do_mir) {
1393         
1394                 i= shi->ray_mirror*fresnel_fac(shi->view, shi->vn, shi->mat->fresnel_mir_i, shi->mat->fresnel_mir);
1395                 if(i!=0.0f) {
1396                 
1397                         trace_reflect(mircol, shi, shr, i);
1398                         
1399                         fr= i*shi->mirr;
1400                         fg= i*shi->mirg;
1401                         fb= i*shi->mirb;
1402
1403                         if(shi->passflag & SCE_PASS_REFLECT) {
1404                                 /* mirror pass is not blocked out with spec */
1405                                 shr->refl[0]= fr*mircol[0] - fr*diff[0];
1406                                 shr->refl[1]= fg*mircol[1] - fg*diff[1];
1407                                 shr->refl[2]= fb*mircol[2] - fb*diff[2];
1408                         }
1409                         
1410                         if(shi->combinedflag & SCE_PASS_REFLECT) {
1411                                 /* values in shr->spec can be greater then 1.0.
1412                                  * In this case the mircol uses a zero blending factor, so ignoring it is ok.
1413                                  * Fixes bug #18837 - when the spec is higher then 1.0,
1414                                  * diff can become a negative color - Campbell  */
1415                                 
1416                                 f1= 1.0f-i;
1417                                 
1418                                 diff[0] *= f1;
1419                                 diff[1] *= f1;
1420                                 diff[2] *= f1;
1421                                 
1422                                 if(shr->spec[0]<1.0f)   diff[0] += mircol[0] * (fr*(1.0f-shr->spec[0]));
1423                                 if(shr->spec[1]<1.0f)   diff[1] += mircol[1] * (fg*(1.0f-shr->spec[1]));
1424                                 if(shr->spec[2]<1.0f)   diff[2] += mircol[2] * (fb*(1.0f-shr->spec[2]));
1425                         }
1426                 }
1427         }
1428         /* put back together */
1429         if(shi->combinedflag & SCE_PASS_SPEC)
1430                 VECADD(shr->combined, diff, shr->spec) /* no ; */
1431         else
1432                 VECCOPY(shr->combined, diff);
1433 }
1434
1435 /* color 'shadfac' passes through 'col' with alpha and filter */
1436 /* filter is only applied on alpha defined transparent part */
1437 static void addAlphaLight(float *shadfac, float *col, float alpha, float filter)
1438 {
1439         float fr, fg, fb;
1440         
1441         fr= 1.0f+ filter*(col[0]-1.0f);
1442         fg= 1.0f+ filter*(col[1]-1.0f);
1443         fb= 1.0f+ filter*(col[2]-1.0f);
1444         
1445         shadfac[0]= alpha*col[0] + fr*(1.0f-alpha)*shadfac[0];
1446         shadfac[1]= alpha*col[1] + fg*(1.0f-alpha)*shadfac[1];
1447         shadfac[2]= alpha*col[2] + fb*(1.0f-alpha)*shadfac[2];
1448         
1449         shadfac[3]= (1.0f-alpha)*shadfac[3];
1450 }
1451
1452 static void ray_trace_shadow_tra(Isect *is, ShadeInput *origshi, int depth, int traflag)
1453 {
1454         /* ray to lamp, find first face that intersects, check alpha properties,
1455            if it has col[3]>0.0f  continue. so exit when alpha is full */
1456         ShadeInput shi;
1457         ShadeResult shr;
1458         float initial_labda = is->labda;
1459         
1460         if(RE_rayobject_raycast(R.raytree, is)) {
1461                 float d= 1.0f;
1462                 /* we got a face */
1463                 
1464                 /* Warning, This is not that nice, and possibly a bit slow for every ray,
1465                 however some variables were not initialized properly in, unless using shade_input_initialize(...), we need to do a memset */
1466                 memset(&shi, 0, sizeof(ShadeInput)); 
1467                 /* end warning! - Campbell */
1468                 
1469                 shi.depth= 1;                                   /* only used to indicate tracing */
1470                 shi.mask= origshi->mask;
1471                 shi.thread= origshi->thread;
1472                 shi.passflag= SCE_PASS_COMBINED;
1473                 shi.combinedflag= 0xFFFFFF;              /* ray trace does all options */
1474         
1475                 shi.xs= origshi->xs;
1476                 shi.ys= origshi->ys;
1477                 shi.lay= origshi->lay;
1478                 shi.nodes= origshi->nodes;
1479                 
1480                 shade_ray(is, &shi, &shr);
1481                 if (shi.mat->material_type == MA_TYPE_SURFACE) {
1482                         if (traflag & RAY_TRA)
1483                                 d= shade_by_transmission(is, &shi, &shr);
1484                         
1485                         /* mix colors based on shadfac (rgb + amount of light factor) */
1486                         addAlphaLight(is->col, shr.diff, shr.alpha, d*shi.mat->filter);
1487                 } else if (shi.mat->material_type == MA_TYPE_VOLUME) {
1488                         addAlphaLight(is->col, shr.combined, shr.alpha, 1.0f);
1489                 }
1490                 
1491                 if(depth>0 && is->col[3]>0.0f) {
1492                         
1493                         /* adapt isect struct */
1494                         VECCOPY(is->start, shi.co);
1495                         is->labda = initial_labda-is->labda;
1496                         is->orig.ob   = shi.obi;
1497                         is->orig.face = shi.vlr;
1498
1499                         ray_trace_shadow_tra(is, origshi, depth-1, traflag | RAY_TRA);
1500                 }
1501                 
1502                 RE_RC_MERGE(&origshi->raycounter, &shi.raycounter);
1503         }
1504 }
1505
1506 /* not used, test function for ambient occlusion (yaf: pathlight) */
1507 /* main problem; has to be called within shading loop, giving unwanted recursion */
1508 int ray_trace_shadow_rad(ShadeInput *ship, ShadeResult *shr)
1509 {
1510         static int counter=0, only_one= 0;
1511         extern float hashvectf[];
1512         Isect isec;
1513         ShadeInput shi;
1514         ShadeResult shr_t;
1515         float vec[3], accum[3], div= 0.0f;
1516         int a;
1517         
1518         assert(0);
1519         
1520         if(only_one) {
1521                 return 0;
1522         }
1523         only_one= 1;
1524         
1525         accum[0]= accum[1]= accum[2]= 0.0f;
1526         isec.mode= RE_RAY_MIRROR;
1527         isec.orig.ob   = ship->obi;
1528         isec.orig.face = ship->vlr;
1529         isec.hint = 0;
1530
1531         VECCOPY(isec.start, ship->co);
1532         
1533         RE_RC_INIT(isec, shi);
1534         
1535         for(a=0; a<8*8; a++) {
1536                 
1537                 counter+=3;
1538                 counter %= 768;
1539                 VECCOPY(vec, hashvectf+counter);
1540                 if(ship->vn[0]*vec[0]+ship->vn[1]*vec[1]+ship->vn[2]*vec[2]>0.0f) {
1541                         vec[0]-= vec[0];
1542                         vec[1]-= vec[1];
1543                         vec[2]-= vec[2];
1544                 }
1545
1546                 VECCOPY(isec.vec, vec );
1547                 isec.labda = RE_RAYTRACE_MAXDIST;
1548
1549                 if(RE_rayobject_raycast(R.raytree, &isec)) {
1550                         float fac;
1551                         
1552                         /* Warning, This is not that nice, and possibly a bit slow for every ray,
1553                         however some variables were not initialized properly in, unless using shade_input_initialize(...), we need to do a memset */
1554                         memset(&shi, 0, sizeof(ShadeInput)); 
1555                         /* end warning! - Campbell */
1556                         
1557                         shade_ray(&isec, &shi, &shr_t);
1558                         fac= isec.labda*isec.labda;
1559                         fac= 1.0f;
1560                         accum[0]+= fac*(shr_t.diff[0]+shr_t.spec[0]);
1561                         accum[1]+= fac*(shr_t.diff[1]+shr_t.spec[1]);
1562                         accum[2]+= fac*(shr_t.diff[2]+shr_t.spec[2]);
1563                         div+= fac;
1564                 }
1565                 else div+= 1.0f;
1566         }
1567         
1568         if(div!=0.0f) {
1569                 shr->diff[0]+= accum[0]/div;
1570                 shr->diff[1]+= accum[1]/div;
1571                 shr->diff[2]+= accum[2]/div;
1572         }
1573         shr->alpha= 1.0f;
1574         
1575         only_one= 0;
1576         return 1;
1577 }
1578
1579 /* aolight: function to create random unit sphere vectors for total random sampling */
1580 static void RandomSpherical(float *v)
1581 {
1582         float r;
1583         v[2] = 2.f*BLI_frand()-1.f;
1584         if ((r = 1.f - v[2]*v[2])>0.f) {
1585                 float a = 6.283185307f*BLI_frand();
1586                 r = sqrt(r);
1587                 v[0] = r * cos(a);
1588                 v[1] = r * sin(a);
1589         }
1590         else v[2] = 1.f;
1591 }
1592
1593 /* calc distributed spherical energy */
1594 static void DS_energy(float *sphere, int tot, float *vec)
1595 {
1596         float *fp, fac, force[3], res[3];
1597         int a;
1598         
1599         res[0]= res[1]= res[2]= 0.0f;
1600         
1601         for(a=0, fp=sphere; a<tot; a++, fp+=3) {
1602                 VecSubf(force, vec, fp);
1603                 fac= force[0]*force[0] + force[1]*force[1] + force[2]*force[2];
1604                 if(fac!=0.0f) {
1605                         fac= 1.0f/fac;
1606                         res[0]+= fac*force[0];
1607                         res[1]+= fac*force[1];
1608                         res[2]+= fac*force[2];
1609                 }
1610         }
1611
1612         VecMulf(res, 0.5);
1613         VecAddf(vec, vec, res);
1614         Normalize(vec);
1615         
1616 }
1617
1618 /* called from convertBlenderScene.c */
1619 /* creates an equally distributed spherical sample pattern */
1620 /* and allocates threadsafe memory */
1621 void init_ao_sphere(World *wrld)
1622 {
1623         float *fp;
1624         int a, tot, iter= 16;
1625
1626         /* we make twice the amount of samples, because only a hemisphere is used */
1627         tot= 2*wrld->aosamp*wrld->aosamp;
1628         
1629         wrld->aosphere= MEM_mallocN(3*tot*sizeof(float), "AO sphere");
1630         
1631         /* fixed random */
1632         BLI_srandom(tot);
1633         
1634         /* init */
1635         fp= wrld->aosphere;
1636         for(a=0; a<tot; a++, fp+= 3) {
1637                 RandomSpherical(fp);
1638         }
1639         
1640         while(iter--) {
1641                 for(a=0, fp= wrld->aosphere; a<tot; a++, fp+= 3) {
1642                         DS_energy(wrld->aosphere, tot, fp);
1643                 }
1644         }
1645         
1646         /* tables */
1647         wrld->aotables= MEM_mallocN(BLENDER_MAX_THREADS*3*tot*sizeof(float), "AO tables");
1648 }
1649
1650 /* give per thread a table, we have to compare xs ys because of way OSA works... */
1651 static float *threadsafe_table_sphere(int test, int thread, int xs, int ys, int tot)
1652 {
1653         static int xso[BLENDER_MAX_THREADS], yso[BLENDER_MAX_THREADS];
1654         static int firsttime= 1;
1655         
1656         if(firsttime) {
1657                 memset(xso, 255, sizeof(xso));
1658                 memset(yso, 255, sizeof(yso));
1659                 firsttime= 0;
1660         }
1661         
1662         if(xs==xso[thread] && ys==yso[thread]) return R.wrld.aotables+ thread*tot*3;
1663         if(test) return NULL;
1664         xso[thread]= xs; yso[thread]= ys;
1665         return R.wrld.aotables+ thread*tot*3;
1666 }
1667
1668 static float *sphere_sampler(int type, int resol, int thread, int xs, int ys)
1669 {
1670         int tot;
1671         float *vec;
1672         
1673         tot= 2*resol*resol;
1674
1675         if (type & WO_AORNDSMP) {
1676                 float *sphere;
1677                 int a;
1678                 
1679                 // always returns table
1680                 sphere= threadsafe_table_sphere(0, thread, xs, ys, tot);
1681
1682                 /* total random sampling. NOT THREADSAFE! (should be removed, is not useful) */
1683                 vec= sphere;
1684                 for (a=0; a<tot; a++, vec+=3) {
1685                         RandomSpherical(vec);
1686                 }
1687                 
1688                 return sphere;
1689         } 
1690         else {
1691                 float *sphere;
1692                 float cosfi, sinfi, cost, sint;
1693                 float ang, *vec1;
1694                 int a;
1695                 
1696                 // returns table if xs and ys were equal to last call
1697                 sphere= threadsafe_table_sphere(1, thread, xs, ys, tot);
1698                 if(sphere==NULL) {
1699                         sphere= threadsafe_table_sphere(0, thread, xs, ys, tot);
1700                         
1701                         // random rotation
1702                         ang= BLI_thread_frand(thread);
1703                         sinfi= sin(ang); cosfi= cos(ang);
1704                         ang= BLI_thread_frand(thread);
1705                         sint= sin(ang); cost= cos(ang);
1706                         
1707                         vec= R.wrld.aosphere;
1708                         vec1= sphere;
1709                         for (a=0; a<tot; a++, vec+=3, vec1+=3) {
1710                                 vec1[0]= cost*cosfi*vec[0] - sinfi*vec[1] + sint*cosfi*vec[2];
1711                                 vec1[1]= cost*sinfi*vec[0] + cosfi*vec[1] + sint*sinfi*vec[2];
1712                                 vec1[2]= -sint*vec[0] + cost*vec[2];                    
1713                         }
1714                 }
1715                 return sphere;
1716         }
1717 }
1718
1719 static void ray_ao_qmc(ShadeInput *shi, float *shadfac)
1720 {
1721         Isect isec;
1722         RayHint point_hint;
1723         QMCSampler *qsa=NULL;
1724         float samp3d[3];
1725         float up[3], side[3], dir[3], nrm[3];
1726         
1727         float maxdist = R.wrld.aodist;
1728         float fac=0.0f, prev=0.0f;
1729         float adapt_thresh = R.wrld.ao_adapt_thresh;
1730         float adapt_speed_fac = R.wrld.ao_adapt_speed_fac;
1731         
1732         int samples=0;
1733         int max_samples = R.wrld.aosamp*R.wrld.aosamp;
1734         
1735         float dxyview[3], skyadded=0, div;
1736         int aocolor;
1737         
1738         RE_RC_INIT(isec, *shi);
1739         isec.orig.ob   = shi->obi;
1740         isec.orig.face = shi->vlr;
1741         isec.skip = RE_SKIP_VLR_NEIGHBOUR | RE_SKIP_VLR_RENDER_CHECK | RE_SKIP_VLR_NON_SOLID_MATERIAL;
1742         isec.hint = 0;
1743
1744         isec.hit.ob   = 0;
1745         isec.hit.face = 0;
1746
1747         isec.last_hit = NULL;
1748         
1749         isec.mode= (R.wrld.aomode & WO_AODIST)?RE_RAY_SHADOW_TRA:RE_RAY_SHADOW;
1750         isec.lay= -1;
1751         
1752         VECCOPY(isec.start, shi->co);           
1753         RE_rayobject_hint_bb( R.raytree, &point_hint, isec.start, isec.start );
1754         isec.hint = &point_hint;
1755
1756         
1757         shadfac[0]= shadfac[1]= shadfac[2]= 0.0f;
1758         
1759         /* prevent sky colors to be added for only shadow (shadow becomes alpha) */
1760         aocolor= R.wrld.aocolor;
1761         if(shi->mat->mode & MA_ONLYSHADOW)
1762                 aocolor= WO_AOPLAIN;
1763         
1764         if(aocolor == WO_AOSKYTEX) {
1765                 dxyview[0]= 1.0f/(float)R.wrld.aosamp;
1766                 dxyview[1]= 1.0f/(float)R.wrld.aosamp;
1767                 dxyview[2]= 0.0f;
1768         }
1769         
1770         if(shi->vlr->flag & R_SMOOTH) {
1771                 VECCOPY(nrm, shi->vn);
1772         }
1773         else {
1774                 VECCOPY(nrm, shi->facenor);
1775         }
1776         
1777         VecOrthoBasisf(nrm, up, side);
1778         
1779         /* sampling init */
1780         if (R.wrld.ao_samp_method==WO_AOSAMP_HALTON) {
1781                 float speedfac;
1782                 
1783                 speedfac = get_avg_speed(shi) * adapt_speed_fac;
1784                 CLAMP(speedfac, 1.0, 1000.0);
1785                 max_samples /= speedfac;
1786                 if (max_samples < 5) max_samples = 5;
1787                 
1788                 qsa = get_thread_qmcsampler(&R, shi->thread, SAMP_TYPE_HALTON, max_samples);
1789         } else if (R.wrld.ao_samp_method==WO_AOSAMP_HAMMERSLEY)
1790                 qsa = get_thread_qmcsampler(&R, shi->thread, SAMP_TYPE_HAMMERSLEY, max_samples);
1791
1792         QMC_initPixel(qsa, shi->thread);
1793         
1794         
1795         while (samples < max_samples) {
1796
1797                 /* sampling, returns quasi-random vector in unit hemisphere */
1798                 QMC_sampleHemi(samp3d, qsa, shi->thread, samples);
1799
1800                 dir[0] = (samp3d[0]*up[0] + samp3d[1]*side[0] + samp3d[2]*nrm[0]);
1801                 dir[1] = (samp3d[0]*up[1] + samp3d[1]*side[1] + samp3d[2]*nrm[1]);
1802                 dir[2] = (samp3d[0]*up[2] + samp3d[1]*side[2] + samp3d[2]*nrm[2]);
1803                 
1804                 Normalize(dir);
1805                         
1806                 isec.vec[0] = -dir[0];
1807                 isec.vec[1] = -dir[1];
1808                 isec.vec[2] = -dir[2];
1809                 isec.labda = maxdist;
1810                 
1811                 prev = fac;
1812                 
1813                 if(RE_rayobject_raycast(R.raytree, &isec)) {
1814                         if (R.wrld.aomode & WO_AODIST) fac+= exp(-isec.labda*R.wrld.aodistfac); 
1815                         else fac+= 1.0f;
1816                 }
1817                 else if(aocolor!=WO_AOPLAIN) {
1818                         float skycol[4];
1819                         float skyfac, view[3];
1820                         
1821                         view[0]= -dir[0];
1822                         view[1]= -dir[1];
1823                         view[2]= -dir[2];
1824                         Normalize(view);
1825                         
1826                         if(aocolor==WO_AOSKYCOL) {
1827                                 skyfac= 0.5*(1.0f+view[0]*R.grvec[0]+ view[1]*R.grvec[1]+ view[2]*R.grvec[2]);
1828                                 shadfac[0]+= (1.0f-skyfac)*R.wrld.horr + skyfac*R.wrld.zenr;
1829                                 shadfac[1]+= (1.0f-skyfac)*R.wrld.horg + skyfac*R.wrld.zeng;
1830                                 shadfac[2]+= (1.0f-skyfac)*R.wrld.horb + skyfac*R.wrld.zenb;
1831                         }
1832                         else {  /* WO_AOSKYTEX */
1833                                 shadeSkyView(skycol, isec.start, view, dxyview, shi->thread);
1834                                 shadeSunView(skycol, shi->view);
1835                                 shadfac[0]+= skycol[0];
1836                                 shadfac[1]+= skycol[1];
1837                                 shadfac[2]+= skycol[2];
1838                         }
1839                         skyadded++;
1840                 }
1841                 
1842                 samples++;
1843                 
1844                 if (qsa->type == SAMP_TYPE_HALTON) {
1845                         /* adaptive sampling - consider samples below threshold as in shadow (or vice versa) and exit early */          
1846                         if (adapt_thresh > 0.0 && (samples > max_samples/2) ) {
1847                                 
1848                                 if (adaptive_sample_contrast_val(samples, prev, fac, adapt_thresh)) {
1849                                         break;
1850                                 }
1851                         }
1852                 }
1853         }
1854         
1855         if(aocolor!=WO_AOPLAIN && skyadded) {
1856                 div= (1.0f - fac/(float)samples)/((float)skyadded);
1857                 
1858                 shadfac[0]*= div;       // average color times distances/hits formula
1859                 shadfac[1]*= div;       // average color times distances/hits formula
1860                 shadfac[2]*= div;       // average color times distances/hits formula
1861         } else {
1862                 shadfac[0]= shadfac[1]= shadfac[2]= 1.0f - fac/(float)samples;
1863         }
1864         
1865         if (qsa)
1866                 release_thread_qmcsampler(&R, shi->thread, qsa);
1867 }
1868
1869 /* extern call from shade_lamp_loop, ambient occlusion calculus */
1870 static void ray_ao_spheresamp(ShadeInput *shi, float *shadfac)
1871 {
1872         Isect isec;
1873         RayHint point_hint;
1874         float *vec, *nrm, div, bias, sh=0.0f;
1875         float maxdist = R.wrld.aodist;
1876         float dxyview[3];
1877         int j= -1, tot, actual=0, skyadded=0, aocolor, resol= R.wrld.aosamp;
1878         
1879         RE_RC_INIT(isec, *shi);
1880         isec.orig.ob   = shi->obi;
1881         isec.orig.face = shi->vlr;
1882         isec.skip = RE_SKIP_VLR_NEIGHBOUR | RE_SKIP_VLR_RENDER_CHECK;
1883         isec.hint = 0;
1884
1885         isec.hit.ob   = 0;
1886         isec.hit.face = 0;
1887         
1888         isec.last_hit = NULL;
1889         
1890         isec.mode= (R.wrld.aomode & WO_AODIST)?RE_RAY_SHADOW_TRA:RE_RAY_SHADOW;
1891         isec.lay= -1;
1892
1893         VECCOPY(isec.start, shi->co);           
1894         RE_rayobject_hint_bb( R.raytree, &point_hint, isec.start, isec.start );
1895         isec.hint = &point_hint;
1896
1897         shadfac[0]= shadfac[1]= shadfac[2]= 0.0f;
1898
1899         /* bias prevents smoothed faces to appear flat */
1900         if(shi->vlr->flag & R_SMOOTH) {
1901                 bias= R.wrld.aobias;
1902                 nrm= shi->vn;
1903         }
1904         else {
1905                 bias= 0.0f;
1906                 nrm= shi->facenor;
1907         }
1908
1909         /* prevent sky colors to be added for only shadow (shadow becomes alpha) */
1910         aocolor= R.wrld.aocolor;
1911         if(shi->mat->mode & MA_ONLYSHADOW)
1912                 aocolor= WO_AOPLAIN;
1913         
1914         if(resol>32) resol= 32;
1915         
1916         vec= sphere_sampler(R.wrld.aomode, resol, shi->thread, shi->xs, shi->ys);
1917         
1918         // warning: since we use full sphere now, and dotproduct is below, we do twice as much
1919         tot= 2*resol*resol;
1920
1921         if(aocolor == WO_AOSKYTEX) {
1922                 dxyview[0]= 1.0f/(float)resol;
1923                 dxyview[1]= 1.0f/(float)resol;
1924                 dxyview[2]= 0.0f;
1925         }
1926         
1927         while(tot--) {
1928                 
1929                 if ((vec[0]*nrm[0] + vec[1]*nrm[1] + vec[2]*nrm[2]) > bias) {
1930                         /* only ao samples for mask */
1931                         if(R.r.mode & R_OSA) {
1932                                 j++;
1933                                 if(j==R.osa) j= 0;
1934                                 if(!(shi->mask & (1<<j))) {
1935                                         vec+=3;
1936                                         continue;
1937                                 }
1938                         }
1939                         
1940                         actual++;
1941                         
1942                         /* always set start/vec/labda */
1943                         isec.vec[0] = -vec[0];
1944                         isec.vec[1] = -vec[1];
1945                         isec.vec[2] = -vec[2];
1946                         isec.labda = maxdist;
1947                         
1948                         /* do the trace */
1949                         if(RE_rayobject_raycast(R.raytree, &isec)) {
1950                                 if (R.wrld.aomode & WO_AODIST) sh+= exp(-isec.labda*R.wrld.aodistfac); 
1951                                 else sh+= 1.0f;
1952                         }
1953                         else if(aocolor!=WO_AOPLAIN) {
1954                                 float skycol[4];
1955                                 float fac, view[3];
1956                                 
1957                                 view[0]= -vec[0];
1958                                 view[1]= -vec[1];
1959                                 view[2]= -vec[2];
1960                                 Normalize(view);
1961                                 
1962                                 if(aocolor==WO_AOSKYCOL) {
1963                                         fac= 0.5*(1.0f+view[0]*R.grvec[0]+ view[1]*R.grvec[1]+ view[2]*R.grvec[2]);
1964                                         shadfac[0]+= (1.0f-fac)*R.wrld.horr + fac*R.wrld.zenr;
1965                                         shadfac[1]+= (1.0f-fac)*R.wrld.horg + fac*R.wrld.zeng;
1966                                         shadfac[2]+= (1.0f-fac)*R.wrld.horb + fac*R.wrld.zenb;
1967                                 }
1968                                 else {  /* WO_AOSKYTEX */
1969                                         shadeSkyView(skycol, isec.start, view, dxyview, shi->thread);
1970                                         shadeSunView(skycol, shi->view);
1971                                         shadfac[0]+= skycol[0];
1972                                         shadfac[1]+= skycol[1];
1973                                         shadfac[2]+= skycol[2];
1974                                 }
1975                                 skyadded++;
1976                         }
1977                 }
1978                 // samples
1979                 vec+= 3;
1980         }
1981         
1982         if(actual==0) sh= 1.0f;
1983         else sh = 1.0f - sh/((float)actual);
1984         
1985         if(aocolor!=WO_AOPLAIN && skyadded) {
1986                 div= sh/((float)skyadded);
1987                 
1988                 shadfac[0]*= div;       // average color times distances/hits formula
1989                 shadfac[1]*= div;       // average color times distances/hits formula
1990                 shadfac[2]*= div;       // average color times distances/hits formula
1991         }
1992         else {
1993                 shadfac[0]= shadfac[1]= shadfac[2]= sh;
1994         }
1995 }
1996
1997 void ray_ao(ShadeInput *shi, float *shadfac)
1998 {
1999         /* Unfortunately, the unusual way that the sphere sampler calculates roughly twice as many
2000          * samples as are actually traced, and skips them based on bias and OSA settings makes it very difficult
2001          * to reuse code between these two functions. This is the easiest way I can think of to do it
2002          * --broken */
2003         if (ELEM(R.wrld.ao_samp_method, WO_AOSAMP_HAMMERSLEY, WO_AOSAMP_HALTON))
2004                 ray_ao_qmc(shi, shadfac);
2005         else if (R.wrld.ao_samp_method == WO_AOSAMP_CONSTANT)
2006                 ray_ao_spheresamp(shi, shadfac);
2007 }
2008
2009 static void ray_shadow_jittered_coords(ShadeInput *shi, int max, float jitco[RE_MAX_OSA][3], int *totjitco)
2010 {
2011         /* magic numbers for reordering sample positions to give better
2012          * results with adaptive sample, when it usually only takes 4 samples */
2013         int order8[8] = {0, 1, 5, 6, 2, 3, 4, 7};
2014         int order11[11] = {1, 3, 8, 10, 0, 2, 4, 5, 6, 7, 9};
2015         int order16[16] = {1, 3, 9, 12, 0, 6, 7, 8, 13, 2, 4, 5, 10, 11, 14, 15};
2016         int count = count_mask(shi->mask);
2017
2018         /* for better antialising shadow samples are distributed over the subpixel
2019          * sample coordinates, this only works for raytracing depth 0 though */
2020         if(!shi->strand && shi->depth == 0 && count > 1 && count <= max) {
2021                 float xs, ys, zs, view[3];
2022                 int samp, ordsamp, tot= 0;
2023
2024                 for(samp=0; samp<R.osa; samp++) {
2025                         if(R.osa == 8) ordsamp = order8[samp];
2026                         else if(R.osa == 11) ordsamp = order11[samp];
2027                         else if(R.osa == 16) ordsamp = order16[samp];
2028                         else ordsamp = samp;
2029
2030                         if(shi->mask & (1<<ordsamp)) {
2031                                 /* zbuffer has this inverse corrected, ensures xs,ys are inside pixel */
2032                                 xs= (float)shi->scanco[0] + R.jit[ordsamp][0] + 0.5f;
2033                                 ys= (float)shi->scanco[1] + R.jit[ordsamp][1] + 0.5f;
2034                                 zs= shi->scanco[2];
2035
2036                                 shade_input_calc_viewco(shi, xs, ys, zs, view, NULL, jitco[tot], NULL, NULL);
2037                                 tot++;
2038                         }
2039                 }
2040
2041                 *totjitco= tot;
2042         }
2043         else {
2044                 VECCOPY(jitco[0], shi->co);
2045                 *totjitco= 1;
2046         }
2047 }
2048
2049 static void ray_shadow_qmc(ShadeInput *shi, LampRen *lar, float *lampco, float *shadfac, Isect *isec)
2050 {
2051         QMCSampler *qsa=NULL;
2052         int samples=0;
2053         float samp3d[3];
2054
2055         float fac=0.0f, vec[3], end[3];
2056         float colsq[4];
2057         float adapt_thresh = lar->adapt_thresh;
2058         int min_adapt_samples=4, max_samples = lar->ray_totsamp;
2059         float *co;
2060         int do_soft=1, full_osa=0, i;
2061
2062         float min[3], max[3];
2063         RayHint bb_hint;
2064
2065         float jitco[RE_MAX_OSA][3];
2066         int totjitco;
2067
2068         colsq[0] = colsq[1] = colsq[2] = 0.0;
2069         if(isec->mode==RE_RAY_SHADOW_TRA) {
2070                 shadfac[0]= shadfac[1]= shadfac[2]= shadfac[3]= 0.0f;
2071         } else
2072                 shadfac[3]= 1.0f;
2073         
2074         if (lar->ray_totsamp < 2) do_soft = 0;
2075         if ((R.r.mode & R_OSA) && (R.osa > 0) && (shi->vlr->flag & R_FULL_OSA)) full_osa = 1;
2076         
2077         if (full_osa) {
2078                 if (do_soft) max_samples  = max_samples/R.osa + 1;
2079                 else max_samples = 1;
2080         } else {
2081                 if (do_soft) max_samples = lar->ray_totsamp;
2082                 else max_samples = (R.osa > 4)?R.osa:5;
2083         }
2084         
2085         ray_shadow_jittered_coords(shi, max_samples, jitco, &totjitco);
2086
2087         /* sampling init */
2088         if (lar->ray_samp_method==LA_SAMP_HALTON)
2089                 qsa = get_thread_qmcsampler(&R, shi->thread, SAMP_TYPE_HALTON, max_samples);
2090         else if (lar->ray_samp_method==LA_SAMP_HAMMERSLEY)
2091                 qsa = get_thread_qmcsampler(&R, shi->thread, SAMP_TYPE_HAMMERSLEY, max_samples);
2092         
2093         QMC_initPixel(qsa, shi->thread);
2094
2095         INIT_MINMAX(min, max);
2096         for(i=0; i<totjitco; i++)
2097         {
2098                 DO_MINMAX(jitco[i], min, max);
2099         }
2100         RE_rayobject_hint_bb( R.raytree, &bb_hint, min, max);
2101         
2102         isec->hint = &bb_hint;
2103         isec->skip = RE_SKIP_VLR_NEIGHBOUR | RE_SKIP_VLR_RENDER_CHECK;
2104         VECCOPY(vec, lampco);
2105         
2106         while (samples < max_samples) {
2107
2108                 isec->orig.ob   = shi->obi;
2109                 isec->orig.face = shi->vlr;
2110
2111                 /* manually jitter the start shading co-ord per sample
2112                  * based on the pre-generated OSA texture sampling offsets, 
2113                  * for anti-aliasing sharp shadow edges. */
2114                 co = jitco[samples % totjitco];
2115
2116                 if (do_soft) {
2117                         /* sphere shadow source */
2118                         if (lar->type == LA_LOCAL) {
2119                                 float ru[3], rv[3], v[3], s[3];
2120                                 
2121                                 /* calc tangent plane vectors */
2122                                 v[0] = co[0] - lampco[0];
2123                                 v[1] = co[1] - lampco[1];
2124                                 v[2] = co[2] - lampco[2];
2125                                 Normalize(v);
2126                                 VecOrthoBasisf(v, ru, rv);
2127                                 
2128                                 /* sampling, returns quasi-random vector in area_size disc */
2129                                 QMC_sampleDisc(samp3d, qsa, shi->thread, samples,lar->area_size);
2130
2131                                 /* distribute disc samples across the tangent plane */
2132                                 s[0] = samp3d[0]*ru[0] + samp3d[1]*rv[0];
2133                                 s[1] = samp3d[0]*ru[1] + samp3d[1]*rv[1];
2134                                 s[2] = samp3d[0]*ru[2] + samp3d[1]*rv[2];
2135                                 
2136                                 VECCOPY(samp3d, s);
2137                         }
2138                         else {
2139                                 /* sampling, returns quasi-random vector in [sizex,sizey]^2 plane */
2140                                 QMC_sampleRect(samp3d, qsa, shi->thread, samples, lar->area_size, lar->area_sizey);
2141                                                                 
2142                                 /* align samples to lamp vector */
2143                                 Mat3MulVecfl(lar->mat, samp3d);
2144                         }
2145                         end[0] = vec[0]+samp3d[0];
2146                         end[1] = vec[1]+samp3d[1];
2147                         end[2] = vec[2]+samp3d[2];
2148                 } else {
2149                         VECCOPY(end, vec);
2150                 }
2151
2152                 if(shi->strand) {
2153                         /* bias away somewhat to avoid self intersection */
2154                         float jitbias= 0.5f*(VecLength(shi->dxco) + VecLength(shi->dyco));
2155                         float v[3];
2156
2157                         VECSUB(v, co, end);
2158                         Normalize(v);
2159
2160                         co[0] -= jitbias*v[0];
2161                         co[1] -= jitbias*v[1];
2162                         co[2] -= jitbias*v[2];
2163                 }
2164
2165                 VECCOPY(isec->start, co);
2166                 isec->vec[0] = end[0]-isec->start[0];
2167                 isec->vec[1] = end[1]-isec->start[1];
2168                 isec->vec[2] = end[2]-isec->start[2];
2169                 isec->labda = 1.0f; // * Normalize(isec->vec);
2170                 
2171                 /* trace the ray */
2172                 if(isec->mode==RE_RAY_SHADOW_TRA) {
2173                         isec->col[0]= isec->col[1]= isec->col[2]=  1.0f;
2174                         isec->col[3]= 1.0f;
2175                         
2176                         ray_trace_shadow_tra(isec, shi, DEPTH_SHADOW_TRA, 0);
2177                         shadfac[0] += isec->col[0];
2178                         shadfac[1] += isec->col[1];
2179                         shadfac[2] += isec->col[2];
2180                         shadfac[3] += isec->col[3];
2181                         
2182                         /* for variance calc */
2183                         colsq[0] += isec->col[0]*isec->col[0];
2184                         colsq[1] += isec->col[1]*isec->col[1];
2185                         colsq[2] += isec->col[2]*isec->col[2];
2186                 }
2187                 else {
2188                         if( RE_rayobject_raycast(R.raytree, isec) ) fac+= 1.0f;
2189                 }
2190                 
2191                 samples++;
2192                 
2193                 if ((lar->ray_samp_method == LA_SAMP_HALTON)) {
2194                 
2195                         /* adaptive sampling - consider samples below threshold as in shadow (or vice versa) and exit early */
2196                         if ((max_samples > min_adapt_samples) && (adapt_thresh > 0.0) && (samples > max_samples / 3)) {
2197                                 if (isec->mode==RE_RAY_SHADOW_TRA) {
2198                                         if ((shadfac[3] / samples > (1.0-adapt_thresh)) || (shadfac[3] / samples < adapt_thresh))
2199                                                 break;
2200                                         else if (adaptive_sample_variance(samples, shadfac, colsq, adapt_thresh))
2201                                                 break;
2202                                 } else {
2203                                         if ((fac / samples > (1.0-adapt_thresh)) || (fac / samples < adapt_thresh))
2204                                                 break;
2205                                 }
2206                         }
2207                 }
2208         }
2209         
2210         if(isec->mode==RE_RAY_SHADOW_TRA) {
2211                 shadfac[0] /= samples;
2212                 shadfac[1] /= samples;
2213                 shadfac[2] /= samples;
2214                 shadfac[3] /= samples;
2215         } else
2216                 shadfac[3]= 1.0f-fac/samples;
2217
2218         if (qsa)
2219                 release_thread_qmcsampler(&R, shi->thread, qsa);
2220 }
2221
2222 static void ray_shadow_jitter(ShadeInput *shi, LampRen *lar, float *lampco, float *shadfac, Isect *isec)
2223 {
2224         /* area soft shadow */
2225         float *jitlamp;
2226         float fac=0.0f, div=0.0f, vec[3];
2227         int a, j= -1, mask;
2228         RayHint point_hint;
2229         
2230         if(isec->mode==RE_RAY_SHADOW_TRA) {
2231                 shadfac[0]= shadfac[1]= shadfac[2]= shadfac[3]= 0.0f;
2232         }
2233         else shadfac[3]= 1.0f;
2234         
2235         fac= 0.0f;
2236         jitlamp= give_jitter_plane(lar, shi->thread, shi->xs, shi->ys);
2237
2238         a= lar->ray_totsamp;
2239         
2240         /* this correction to make sure we always take at least 1 sample */
2241         mask= shi->mask;
2242         if(a==4) mask |= (mask>>4)|(mask>>8);
2243         else if(a==9) mask |= (mask>>9);
2244         
2245         VECCOPY(isec->start, shi->co);          
2246         isec->orig.ob   = shi->obi;
2247         isec->orig.face = shi->vlr;
2248         RE_rayobject_hint_bb( R.raytree, &point_hint, isec->start, isec->start );
2249         isec->hint = &point_hint;
2250         
2251         while(a--) {
2252                 
2253                 if(R.r.mode & R_OSA) {
2254                         j++;
2255                         if(j>=R.osa) j= 0;
2256                         if(!(mask & (1<<j))) {
2257                                 jitlamp+= 2;
2258                                 continue;
2259                         }
2260                 }
2261                 
2262                 vec[0]= jitlamp[0];
2263                 vec[1]= jitlamp[1];
2264                 vec[2]= 0.0f;
2265                 Mat3MulVecfl(lar->mat, vec);
2266                 
2267                 /* set start and vec */
2268                 isec->vec[0] = vec[0]+lampco[0]-isec->start[0];
2269                 isec->vec[1] = vec[1]+lampco[1]-isec->start[1];
2270                 isec->vec[2] = vec[2]+lampco[2]-isec->start[2];
2271                 isec->labda = 1.0f;
2272                 isec->skip = RE_SKIP_VLR_NEIGHBOUR | RE_SKIP_VLR_RENDER_CHECK;
2273                 
2274                 if(isec->mode==RE_RAY_SHADOW_TRA) {
2275                         /* isec.col is like shadfac, so defines amount of light (0.0 is full shadow) */
2276                         isec->col[0]= isec->col[1]= isec->col[2]=  1.0f;
2277                         isec->col[3]= 1.0f;
2278                         
2279                         ray_trace_shadow_tra(isec, shi, DEPTH_SHADOW_TRA, 0);
2280                         shadfac[0] += isec->col[0];
2281                         shadfac[1] += isec->col[1];
2282                         shadfac[2] += isec->col[2];
2283                         shadfac[3] += isec->col[3];
2284                 }
2285                 else if( RE_rayobject_raycast(R.raytree, isec) ) fac+= 1.0f;
2286                 
2287                 div+= 1.0f;
2288                 jitlamp+= 2;
2289         }
2290         
2291         if(isec->mode==RE_RAY_SHADOW_TRA) {
2292                 shadfac[0] /= div;
2293                 shadfac[1] /= div;
2294                 shadfac[2] /= div;
2295                 shadfac[3] /= div;
2296         }
2297         else {
2298                 // sqrt makes nice umbra effect
2299                 if(lar->ray_samp_type & LA_SAMP_UMBRA)
2300                         shadfac[3]= sqrt(1.0f-fac/div);
2301                 else
2302                         shadfac[3]= 1.0f-fac/div;
2303         }
2304 }
2305 /* extern call from shade_lamp_loop */
2306 void ray_shadow(ShadeInput *shi, LampRen *lar, float *shadfac)
2307 {
2308         Isect isec;
2309         float lampco[3];
2310
2311         /* setup isec */
2312         RE_RC_INIT(isec, *shi);
2313         if(shi->mat->mode & MA_SHADOW_TRA) isec.mode= RE_RAY_SHADOW_TRA;
2314         else isec.mode= RE_RAY_SHADOW;
2315         isec.hint = 0;
2316         
2317         if(lar->mode & (LA_LAYER|LA_LAYER_SHADOW))
2318                 isec.lay= lar->lay;
2319         else
2320                 isec.lay= -1;
2321
2322         /* only when not mir tracing, first hit optimm */
2323         if(shi->depth==0) {
2324                 isec.last_hit = lar->last_hit[shi->thread];
2325         }
2326         else {
2327                 isec.last_hit = NULL;
2328         }
2329         
2330         if(lar->type==LA_SUN || lar->type==LA_HEMI) {
2331                 /* jitter and QMC sampling add a displace vector to the lamp position
2332                  * that's incorrect because a SUN lamp does not has an exact position
2333                  * and the displace should be done at the ray vector instead of the
2334                  * lamp position.
2335                  * This is easily verified by noticing that shadows of SUN lights change
2336                  * with the scene BB.
2337                  * 
2338                  * This was detected during SoC 2009 - Raytrace Optimization, but to keep
2339                  * consistency with older render code it wasn't removed.
2340                  * 
2341                  * If the render code goes through some recode/serious bug-fix then this
2342                  * is something to consider!
2343                  */
2344                 lampco[0]= shi->co[0] - R.maxdist*lar->vec[0];
2345                 lampco[1]= shi->co[1] - R.maxdist*lar->vec[1];
2346                 lampco[2]= shi->co[2] - R.maxdist*lar->vec[2];
2347         }
2348         else {
2349                 VECCOPY(lampco, lar->co);
2350         }
2351         
2352         if (ELEM(lar->ray_samp_method, LA_SAMP_HALTON, LA_SAMP_HAMMERSLEY)) {
2353                 
2354                 ray_shadow_qmc(shi, lar, lampco, shadfac, &isec);
2355                 
2356         } else {
2357                 if(lar->ray_totsamp<2) {
2358                         
2359                         isec.orig.ob   = shi->obi;
2360                         isec.orig.face = shi->vlr;
2361                         
2362                         shadfac[3]= 1.0f; // 1.0=full light
2363                         
2364                         /* set up isec vec */
2365                         VECCOPY(isec.start, shi->co);
2366                         VECSUB(isec.vec, lampco, isec.start);
2367                         isec.labda = 1.0f;
2368
2369                         if(isec.mode==RE_RAY_SHADOW_TRA) {
2370                                 /* isec.col is like shadfac, so defines amount of light (0.0 is full shadow) */
2371                                 isec.col[0]= isec.col[1]= isec.col[2]=  1.0f;
2372                                 isec.col[3]= 1.0f;
2373
2374                                 ray_trace_shadow_tra(&isec, shi, DEPTH_SHADOW_TRA, 0);
2375                                 QUATCOPY(shadfac, isec.col);
2376                         }
2377                         else if(RE_rayobject_raycast(R.raytree, &isec))
2378                                 shadfac[3]= 0.0f;
2379                 }
2380                 else {
2381                         ray_shadow_jitter(shi, lar, lampco, shadfac, &isec);
2382                 }
2383         }
2384                 
2385         /* for first hit optim, set last interesected shadow face */
2386         if(shi->depth==0) {
2387                 lar->last_hit[shi->thread] = isec.last_hit;
2388         }
2389
2390 }
2391
2392 #if 0
2393 /* only when face points away from lamp, in direction of lamp, trace ray and find first exit point */
2394 static void ray_translucent(ShadeInput *shi, LampRen *lar, float *distfac, float *co)
2395 {
2396         Isect isec;
2397         float lampco[3];
2398         
2399         assert(0);
2400         
2401         /* setup isec */
2402         RE_RC_INIT(isec, *shi);
2403         isec.mode= RE_RAY_SHADOW_TRA;
2404         isec.hint = 0;
2405         
2406         if(lar->mode & LA_LAYER) isec.lay= lar->lay; else isec.lay= -1;
2407         
2408         if(lar->type==LA_SUN || lar->type==LA_HEMI) {
2409                 lampco[0]= shi->co[0] - RE_RAYTRACE_MAXDIST*lar->vec[0];
2410                 lampco[1]= shi->co[1] - RE_RAYTRACE_MAXDIST*lar->vec[1];
2411                 lampco[2]= shi->co[2] - RE_RAYTRACE_MAXDIST*lar->vec[2];
2412         }
2413         else {
2414                 VECCOPY(lampco, lar->co);
2415         }
2416         
2417         isec.orig.ob   = shi->obi;
2418         isec.orig.face = shi->vlr;
2419         
2420         /* set up isec vec */
2421         VECCOPY(isec.start, shi->co);
2422         VECCOPY(isec.end, lampco);
2423         
2424         if(RE_rayobject_raycast(R.raytree, &isec)) {
2425                 /* we got a face */
2426                 
2427                 /* render co */
2428                 co[0]= isec.start[0]+isec.labda*(isec.vec[0]);
2429                 co[1]= isec.start[1]+isec.labda*(isec.vec[1]);
2430                 co[2]= isec.start[2]+isec.labda*(isec.vec[2]);
2431                 
2432                 *distfac= VecLength(isec.vec);
2433         }
2434         else
2435                 *distfac= 0.0f;
2436 }
2437
2438 #endif
2439