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