source/blender/render/intern/source/rayshade.c

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