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