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