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