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