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