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