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