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