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