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