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