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