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

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