Blender  V2.59
rayshade.c
Go to the documentation of this file.
00001 /*
00002  * $Id: rayshade.c 37667 2011-06-20 15:17:02Z campbellbarton $
00003  *
00004  * ***** BEGIN GPL LICENSE BLOCK *****
00005  *
00006  * This program is free software; you can redistribute it and/or
00007  * modify it under the terms of the GNU General Public License
00008  * as published by the Free Software Foundation; either version 2
00009  * of the License, or (at your option) any later version. 
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software Foundation,
00018  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
00019  *
00020  * The Original Code is Copyright (C) 1990-1998 NeoGeo BV.
00021  * All rights reserved.
00022  *
00023  * Contributors: 2004/2005 Blender Foundation, full recode
00024  *
00025  * ***** END GPL LICENSE BLOCK *****
00026  */
00027 
00033 #include <stdio.h>
00034 #include <math.h>
00035 #include <string.h>
00036 #include <stdlib.h>
00037 #include <float.h>
00038 #include <assert.h>
00039 
00040 #include "MEM_guardedalloc.h"
00041 
00042 #include "DNA_material_types.h"
00043 #include "DNA_lamp_types.h"
00044 
00045 #include "BLI_blenlib.h"
00046 #include "BLI_cpu.h"
00047 #include "BLI_jitter.h"
00048 #include "BLI_math.h"
00049 #include "BLI_rand.h"
00050 #include "BLI_utildefines.h"
00051 
00052 #include "BKE_global.h"
00053 #include "BKE_node.h"
00054 
00055 
00056 #include "PIL_time.h"
00057 
00058 #include "render_types.h"
00059 #include "renderpipeline.h"
00060 #include "rendercore.h"
00061 #include "renderdatabase.h"
00062 #include "pixelblending.h"
00063 #include "pixelshading.h"
00064 #include "shading.h"
00065 #include "texture.h"
00066 #include "volumetric.h"
00067 
00068 #include "rayintersection.h"
00069 #include "rayobject.h"
00070 #include "raycounter.h"
00071 
00072 
00073 #define RAY_TRA         1
00074 #define RAY_INSIDE      2
00075 
00076 #define DEPTH_SHADOW_TRA  10
00077 
00078 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
00079 /* defined in pipeline.c, is hardcopy of active dynamic allocated Render */
00080 /* only to be used here in this file, it's for speed */
00081 extern struct Render R;
00082 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
00083 static int test_break(void *data)
00084 {
00085         Render *re = (Render*)data;
00086         return re->test_break(re->tbh);
00087 }
00088 
00089 static void RE_rayobject_config_control(RayObject *r, Render *re)
00090 {
00091         if(RE_rayobject_isRayAPI(r))
00092         {
00093                 r = RE_rayobject_align( r );
00094                 r->control.data = re;
00095                 r->control.test_break = test_break;
00096         }
00097 }
00098 
00099 static RayObject*  RE_rayobject_create(Render *re, int type, int size)
00100 {
00101         RayObject * res = NULL;
00102 
00103         if(type == R_RAYSTRUCTURE_AUTO)
00104         {
00105                 //TODO
00106                 //if(detect_simd())
00107 #ifdef __SSE__
00108                 type = BLI_cpu_support_sse2()? R_RAYSTRUCTURE_SIMD_SVBVH: R_RAYSTRUCTURE_VBVH;
00109 #else
00110                 type = R_RAYSTRUCTURE_VBVH;
00111 #endif
00112         }
00113         
00114 #ifndef __SSE__
00115         if(type == R_RAYSTRUCTURE_SIMD_SVBVH || type == R_RAYSTRUCTURE_SIMD_QBVH)
00116         {
00117                 puts("Warning: Using VBVH (SSE was disabled at compile time)");
00118                 type = R_RAYSTRUCTURE_VBVH;
00119         }
00120 #endif
00121         
00122                 
00123         if(type == R_RAYSTRUCTURE_OCTREE) //TODO dynamic ocres
00124                 res = RE_rayobject_octree_create(re->r.ocres, size);
00125         else if(type == R_RAYSTRUCTURE_BLIBVH)
00126                 res = RE_rayobject_blibvh_create(size);
00127         else if(type == R_RAYSTRUCTURE_VBVH)
00128                 res = RE_rayobject_vbvh_create(size);
00129         else if(type == R_RAYSTRUCTURE_SIMD_SVBVH)
00130                 res = RE_rayobject_svbvh_create(size);
00131         else if(type == R_RAYSTRUCTURE_SIMD_QBVH)
00132                 res = RE_rayobject_qbvh_create(size);
00133         else
00134                 res = RE_rayobject_vbvh_create(size);   //Fallback
00135         
00136         
00137         if(res)
00138                 RE_rayobject_config_control( res, re );
00139         
00140         return res;
00141 }
00142 
00143 #ifdef RE_RAYCOUNTER
00144 RayCounter re_rc_counter[BLENDER_MAX_THREADS];
00145 #endif
00146 
00147 
00148 void freeraytree(Render *re)
00149 {
00150         ObjectInstanceRen *obi;
00151         
00152         if(re->raytree)
00153         {
00154                 RE_rayobject_free(re->raytree);
00155                 re->raytree = NULL;
00156         }
00157         if(re->rayfaces)
00158         {
00159                 MEM_freeN(re->rayfaces);
00160                 re->rayfaces = NULL;
00161         }
00162         if(re->rayprimitives)
00163         {
00164                 MEM_freeN(re->rayprimitives);
00165                 re->rayprimitives = NULL;
00166         }
00167 
00168         for(obi=re->instancetable.first; obi; obi=obi->next)
00169         {
00170                 ObjectRen *obr = obi->obr;
00171                 if(obr->raytree)
00172                 {
00173                         RE_rayobject_free(obr->raytree);
00174                         obr->raytree = NULL;
00175                 }
00176                 if(obr->rayfaces)
00177                 {
00178                         MEM_freeN(obr->rayfaces);
00179                         obr->rayfaces = NULL;
00180                 }
00181                 if(obi->raytree)
00182                 {
00183                         RE_rayobject_free(obi->raytree);
00184                         obi->raytree = NULL;
00185                 }
00186         }
00187         
00188 #ifdef RE_RAYCOUNTER
00189         {
00190                 RayCounter sum;
00191                 memset( &sum, 0, sizeof(sum) );
00192                 int i;
00193                 for(i=0; i<BLENDER_MAX_THREADS; i++)
00194                         RE_RC_MERGE(&sum, re_rc_counter+i);
00195                 RE_RC_INFO(&sum);
00196         }
00197 #endif
00198 }
00199 
00200 static int is_raytraceable_vlr(Render *re, VlakRen *vlr)
00201 {
00202         /* note: volumetric must be tracable, wire must not */
00203         if((re->flag & R_BAKE_TRACE) || (vlr->flag & R_TRACEBLE) || (vlr->mat->material_type == MA_TYPE_VOLUME))
00204                 if(vlr->mat->material_type != MA_TYPE_WIRE)
00205                         return 1;
00206         return 0;
00207 }
00208 
00209 static int is_raytraceable(Render *re, ObjectInstanceRen *obi)
00210 {
00211         int v;
00212         ObjectRen *obr = obi->obr;
00213 
00214         if(re->excludeob && obr->ob == re->excludeob)
00215                 return 0;
00216 
00217         for(v=0;v<obr->totvlak;v++) {
00218                 VlakRen *vlr = obr->vlaknodes[v>>8].vlak + (v&255);
00219 
00220                 if(is_raytraceable_vlr(re, vlr))
00221                         return 1;
00222         }
00223 
00224         return 0;
00225 }
00226 
00227 
00228 RayObject* makeraytree_object(Render *re, ObjectInstanceRen *obi)
00229 {
00230         //TODO
00231         // out-of-memory safeproof
00232         // break render
00233         // update render stats
00234         ObjectRen *obr = obi->obr;
00235         
00236         if(obr->raytree == NULL)
00237         {
00238                 RayObject *raytree;
00239                 RayFace *face = NULL;
00240                 VlakPrimitive *vlakprimitive = NULL;
00241                 int v;
00242                 
00243                 //Count faces
00244                 int faces = 0;
00245                 for(v=0;v<obr->totvlak;v++)
00246                 {
00247                         VlakRen *vlr = obr->vlaknodes[v>>8].vlak + (v&255);
00248                         if(is_raytraceable_vlr(re, vlr))
00249                                 faces++;
00250                 }
00251                 
00252                 if (faces == 0)
00253                         return NULL;
00254 
00255                 //Create Ray cast accelaration structure                
00256                 raytree = RE_rayobject_create( re,  re->r.raytrace_structure, faces );
00257                 if(  (re->r.raytrace_options & R_RAYTRACE_USE_LOCAL_COORDS) )
00258                         vlakprimitive = obr->rayprimitives = (VlakPrimitive*)MEM_callocN(faces*sizeof(VlakPrimitive), "ObjectRen primitives");
00259                 else
00260                         face = obr->rayfaces = (RayFace*)MEM_callocN(faces*sizeof(RayFace), "ObjectRen faces");
00261 
00262                 obr->rayobi = obi;
00263                 
00264                 for(v=0;v<obr->totvlak;v++)
00265                 {
00266                         VlakRen *vlr = obr->vlaknodes[v>>8].vlak + (v&255);
00267                         if(is_raytraceable_vlr(re, vlr))
00268                         {
00269                                 if(  (re->r.raytrace_options & R_RAYTRACE_USE_LOCAL_COORDS) )
00270                                 {
00271                                         RE_rayobject_add( raytree, RE_vlakprimitive_from_vlak( vlakprimitive, obi, vlr ) );
00272                                         vlakprimitive++;
00273                                 }
00274                                 else
00275                                 {
00276                                         RE_rayface_from_vlak( face, obi, vlr );                         
00277                                         RE_rayobject_add( raytree, RE_rayobject_unalignRayFace(face) );
00278                                         face++;
00279                                 }
00280                         }
00281                 }
00282                 RE_rayobject_done( raytree );
00283 
00284                 /* in case of cancel during build, raytree is not usable */
00285                 if(test_break(re))
00286                         RE_rayobject_free(raytree);
00287                 else
00288                         obr->raytree= raytree;
00289         }
00290 
00291         if(obr->raytree) {
00292                 if((obi->flag & R_TRANSFORMED) && obi->raytree == NULL)
00293                 {
00294                         obi->transform_primitives = 0;
00295                         obi->raytree = RE_rayobject_instance_create( obr->raytree, obi->mat, obi, obi->obr->rayobi );
00296                 }
00297         }
00298         
00299         if(obi->raytree) return obi->raytree;
00300         return obi->obr->raytree;
00301 }
00302 
00303 static int has_special_rayobject(Render *re, ObjectInstanceRen *obi)
00304 {
00305         if( (obi->flag & R_TRANSFORMED) && (re->r.raytrace_options & R_RAYTRACE_USE_INSTANCES) )
00306         {
00307                 ObjectRen *obr = obi->obr;
00308                 int v, faces = 0;
00309                 
00310                 for(v=0;v<obr->totvlak;v++)
00311                 {
00312                         VlakRen *vlr = obr->vlaknodes[v>>8].vlak + (v&255);
00313                         if(is_raytraceable_vlr(re, vlr))
00314                         {
00315                                 faces++;
00316                                 if(faces > 4)
00317                                         return 1;
00318                         }
00319                 }
00320         }
00321         return 0;
00322 }
00323 /*
00324  * create a single raytrace structure with all faces
00325  */
00326 static void makeraytree_single(Render *re)
00327 {
00328         ObjectInstanceRen *obi;
00329         RayObject *raytree;
00330         RayFace *face = NULL;
00331         VlakPrimitive *vlakprimitive = NULL;
00332         int faces = 0, obs = 0, special = 0;
00333 
00334         for(obi=re->instancetable.first; obi; obi=obi->next)
00335         if(is_raytraceable(re, obi))
00336         {
00337                 int v;
00338                 ObjectRen *obr = obi->obr;
00339                 obs++;
00340                 
00341                 if(has_special_rayobject(re, obi))
00342                 {
00343                         special++;
00344                 }
00345                 else
00346                 {
00347                         for(v=0;v<obr->totvlak;v++)
00348                         {
00349                                 VlakRen *vlr = obr->vlaknodes[v>>8].vlak + (v&255);
00350                                 if(is_raytraceable_vlr(re, vlr))
00351                                         faces++;
00352                         }
00353                 }
00354         }
00355         
00356         if(faces + special == 0)
00357         {
00358                 re->raytree = RE_rayobject_empty_create();
00359                 return;
00360         }
00361         
00362         //Create raytree
00363         raytree = re->raytree = RE_rayobject_create( re, re->r.raytrace_structure, faces+special );
00364 
00365         if( (re->r.raytrace_options & R_RAYTRACE_USE_LOCAL_COORDS) )
00366         {
00367                 vlakprimitive = re->rayprimitives = (VlakPrimitive*)MEM_callocN(faces*sizeof(VlakPrimitive), "Raytrace vlak-primitives");
00368         }
00369         else
00370         {
00371                 face = re->rayfaces     = (RayFace*)MEM_callocN(faces*sizeof(RayFace), "Render ray faces");
00372         }
00373         
00374         for(obi=re->instancetable.first; obi; obi=obi->next)
00375         if(is_raytraceable(re, obi))
00376         {
00377                 if(test_break(re))
00378                         break;
00379 
00380                 if(has_special_rayobject(re, obi))
00381                 {
00382                         RayObject *obj = makeraytree_object(re, obi);
00383 
00384                         if(test_break(re))
00385                                 break;
00386 
00387                         if (obj)
00388                                 RE_rayobject_add( re->raytree, obj );
00389                 }
00390                 else
00391                 {
00392                         int v;
00393                         ObjectRen *obr = obi->obr;
00394                         
00395                         if(obi->flag & R_TRANSFORMED)
00396                         {
00397                                 obi->transform_primitives = 1;
00398                         }
00399 
00400                         for(v=0;v<obr->totvlak;v++)
00401                         {
00402                                 VlakRen *vlr = obr->vlaknodes[v>>8].vlak + (v&255);
00403                                 if(is_raytraceable_vlr(re, vlr))
00404                                 {
00405                                         if( (re->r.raytrace_options & R_RAYTRACE_USE_LOCAL_COORDS) )
00406                                         {
00407                                                 RayObject *obj = RE_vlakprimitive_from_vlak( vlakprimitive, obi, vlr );
00408                                                 RE_rayobject_add( raytree, obj );
00409                                                 vlakprimitive++;
00410                                         }
00411                                         else
00412                                         {
00413                                                 RE_rayface_from_vlak(face, obi, vlr);
00414                                                 if((obi->flag & R_TRANSFORMED))
00415                                                 {
00416                                                         mul_m4_v3(obi->mat, face->v1);
00417                                                         mul_m4_v3(obi->mat, face->v2);
00418                                                         mul_m4_v3(obi->mat, face->v3);
00419                                                         if(RE_rayface_isQuad(face))
00420                                                                 mul_m4_v3(obi->mat, face->v4);
00421                                                 }
00422 
00423                                                 RE_rayobject_add( raytree, RE_rayobject_unalignRayFace(face) );
00424                                                 face++;
00425                                         }
00426                                 }
00427                         }
00428                 }
00429         }
00430         
00431         if(!test_break(re))
00432         {       
00433                 re->i.infostr= "Raytree.. building";
00434                 re->stats_draw(re->sdh, &re->i);
00435 
00436                 RE_rayobject_done( raytree );   
00437         }
00438 }
00439 
00440 void makeraytree(Render *re)
00441 {
00442         float min[3], max[3], sub[3];
00443         int i;
00444         
00445         re->i.infostr= "Raytree.. preparing";
00446         re->stats_draw(re->sdh, &re->i);
00447 
00448         /* disable options not yet supported by octree,
00449            they might actually never be supported (unless people really need it) */
00450         if(re->r.raytrace_structure == R_RAYSTRUCTURE_OCTREE)
00451                 re->r.raytrace_options &= ~( R_RAYTRACE_USE_INSTANCES | R_RAYTRACE_USE_LOCAL_COORDS);
00452 
00453         makeraytree_single(re);
00454 
00455         if(test_break(re))
00456         {
00457                 freeraytree(re);
00458 
00459                 re->i.infostr= "Raytree building canceled";
00460                 re->stats_draw(re->sdh, &re->i);
00461         }
00462         else
00463         {
00464                 //Calculate raytree max_size
00465                 //This is ONLY needed to kept a bogus behaviour of SUN and HEMI lights
00466                 INIT_MINMAX(min, max);
00467                 RE_rayobject_merge_bb( re->raytree, min, max );
00468                 for(i=0; i<3; i++)
00469                 {
00470                         min[i] += 0.01f;
00471                         max[i] += 0.01f;
00472                         sub[i] = max[i]-min[i];
00473                 }
00474 
00475                 re->maxdist= sub[0]*sub[0] + sub[1]*sub[1] + sub[2]*sub[2];
00476                 if(re->maxdist > 0.0f) re->maxdist= sqrt(re->maxdist);
00477 
00478                 re->i.infostr= "Raytree finished";
00479                 re->stats_draw(re->sdh, &re->i);
00480         }
00481 
00482 #ifdef RE_RAYCOUNTER
00483         memset( re_rc_counter, 0, sizeof(re_rc_counter) );
00484 #endif
00485 }
00486 
00487 /*      if(shi->osatex)  */
00488 static void shade_ray_set_derivative(ShadeInput *shi)
00489 {
00490         float detsh, t00, t10, t01, t11, xn, yn, zn;
00491         int axis1, axis2;
00492 
00493         /* find most stable axis to project */
00494         xn= fabs(shi->facenor[0]);
00495         yn= fabs(shi->facenor[1]);
00496         zn= fabs(shi->facenor[2]);
00497         
00498         if(zn>=xn && zn>=yn) { axis1= 0; axis2= 1; }
00499         else if(yn>=xn && yn>=zn) { axis1= 0; axis2= 2; }
00500         else { axis1= 1; axis2= 2; }
00501         
00502         /* compute u,v and derivatives */
00503         if(shi->obi->flag & R_TRANSFORMED) {
00504                 float v1[3], v2[3], v3[3];
00505 
00506                 mul_v3_m3v3(v1, shi->obi->nmat, shi->v1->co);
00507                 mul_v3_m3v3(v2, shi->obi->nmat, shi->v2->co);
00508                 mul_v3_m3v3(v3, shi->obi->nmat, shi->v3->co);
00509 
00510                 /* same as below */
00511                 t00= v3[axis1]-v1[axis1]; t01= v3[axis2]-v1[axis2];
00512                 t10= v3[axis1]-v2[axis1]; t11= v3[axis2]-v2[axis2];
00513         }
00514         else {
00515                 float *v1= shi->v1->co;
00516                 float *v2= shi->v2->co;
00517                 float *v3= shi->v3->co;
00518 
00519                 /* same as above */
00520                 t00= v3[axis1]-v1[axis1]; t01= v3[axis2]-v1[axis2];
00521                 t10= v3[axis1]-v2[axis1]; t11= v3[axis2]-v2[axis2];
00522         }
00523 
00524         detsh= 1.0f/(t00*t11-t10*t01);
00525         t00*= detsh; t01*=detsh; 
00526         t10*=detsh; t11*=detsh;
00527         
00528         shi->dx_u=  shi->dxco[axis1]*t11- shi->dxco[axis2]*t10;
00529         shi->dx_v=  shi->dxco[axis2]*t00- shi->dxco[axis1]*t01;
00530         shi->dy_u=  shi->dyco[axis1]*t11- shi->dyco[axis2]*t10;
00531         shi->dy_v=  shi->dyco[axis2]*t00- shi->dyco[axis1]*t01;
00532         
00533 }
00534 
00535 
00536 void shade_ray(Isect *is, ShadeInput *shi, ShadeResult *shr)
00537 {
00538         ObjectInstanceRen *obi= (ObjectInstanceRen*)is->hit.ob;
00539         VlakRen *vlr= (VlakRen*)is->hit.face;
00540         
00541         /* set up view vector */
00542         VECCOPY(shi->view, is->dir);
00543 
00544         /* render co */
00545         shi->co[0]= is->start[0]+is->dist*(shi->view[0]);
00546         shi->co[1]= is->start[1]+is->dist*(shi->view[1]);
00547         shi->co[2]= is->start[2]+is->dist*(shi->view[2]);
00548         
00549         normalize_v3(shi->view);
00550 
00551         shi->obi= obi;
00552         shi->obr= obi->obr;
00553         shi->vlr= vlr;
00554         shi->mat= vlr->mat;
00555         shade_input_init_material(shi);
00556         
00557         if(is->isect==2) 
00558                 shade_input_set_triangle_i(shi, obi, vlr, 0, 2, 3);
00559         else
00560                 shade_input_set_triangle_i(shi, obi, vlr, 0, 1, 2);
00561 
00562         shi->u= is->u;
00563         shi->v= is->v;
00564         shi->dx_u= shi->dx_v= shi->dy_u= shi->dy_v=  0.0f;
00565 
00566         if(shi->osatex)
00567                 shade_ray_set_derivative(shi);
00568         shade_input_set_normals(shi);
00569 
00570         shade_input_set_shade_texco(shi);
00571         if (shi->mat->material_type == MA_TYPE_VOLUME) {
00572                 if(ELEM(is->mode, RE_RAY_SHADOW, RE_RAY_SHADOW_TRA)) {
00573                         shade_volume_shadow(shi, shr, is);
00574                 } else {
00575                         shade_volume_outside(shi, shr);
00576                 }
00577         }
00578         else if(is->mode==RE_RAY_SHADOW_TRA) {
00579                 /* temp hack to prevent recursion */
00580                 if(shi->nodes==0 && shi->mat->nodetree && shi->mat->use_nodes) {
00581                         ntreeShaderExecTree(shi->mat->nodetree, shi, shr);
00582                         shi->mat= vlr->mat;             /* shi->mat is being set in nodetree */
00583                 }
00584                 else
00585                         shade_color(shi, shr);
00586         }
00587         else {
00588                 if(shi->mat->nodetree && shi->mat->use_nodes) {
00589                         ntreeShaderExecTree(shi->mat->nodetree, shi, shr);
00590                         shi->mat= vlr->mat;             /* shi->mat is being set in nodetree */
00591                 }
00592                 else {
00593                         shade_material_loop(shi, shr);
00594                 }
00595                 
00596                 /* raytrace likes to separate the spec color */
00597                 VECSUB(shr->diff, shr->combined, shr->spec);
00598         }       
00599 
00600 }
00601 
00602 static int refraction(float *refract, float *n, float *view, float index)
00603 {
00604         float dot, fac;
00605 
00606         VECCOPY(refract, view);
00607         
00608         dot= view[0]*n[0] + view[1]*n[1] + view[2]*n[2];
00609 
00610         if(dot>0.0f) {
00611                 index = 1.0f/index;
00612                 fac= 1.0f - (1.0f - dot*dot)*index*index;
00613                 if(fac<= 0.0f) return 0;
00614                 fac= -dot*index + sqrt(fac);
00615         }
00616         else {
00617                 fac= 1.0f - (1.0f - dot*dot)*index*index;
00618                 if(fac<= 0.0f) return 0;
00619                 fac= -dot*index - sqrt(fac);
00620         }
00621 
00622         refract[0]= index*view[0] + fac*n[0];
00623         refract[1]= index*view[1] + fac*n[1];
00624         refract[2]= index*view[2] + fac*n[2];
00625 
00626         return 1;
00627 }
00628 
00629 /* orn = original face normal */
00630 static void reflection(float *ref, float *n, float *view, float *orn)
00631 {
00632         float f1;
00633         
00634         f1= -2.0f*(n[0]*view[0]+ n[1]*view[1]+ n[2]*view[2]);
00635         
00636         ref[0]= (view[0]+f1*n[0]);
00637         ref[1]= (view[1]+f1*n[1]);
00638         ref[2]= (view[2]+f1*n[2]);
00639 
00640         if(orn) {
00641                 /* test phong normals, then we should prevent vector going to the back */
00642                 f1= ref[0]*orn[0]+ ref[1]*orn[1]+ ref[2]*orn[2];
00643                 if(f1>0.0f) {
00644                         f1+= .01f;
00645                         ref[0]-= f1*orn[0];
00646                         ref[1]-= f1*orn[1];
00647                         ref[2]-= f1*orn[2];
00648                 }
00649         }
00650 }
00651 
00652 #if 0
00653 static void color_combine(float *result, float fac1, float fac2, float *col1, float *col2)
00654 {
00655         float col1t[3], col2t[3];
00656         
00657         col1t[0]= sqrt(col1[0]);
00658         col1t[1]= sqrt(col1[1]);
00659         col1t[2]= sqrt(col1[2]);
00660         col2t[0]= sqrt(col2[0]);
00661         col2t[1]= sqrt(col2[1]);
00662         col2t[2]= sqrt(col2[2]);
00663 
00664         result[0]= (fac1*col1t[0] + fac2*col2t[0]);
00665         result[0]*= result[0];
00666         result[1]= (fac1*col1t[1] + fac2*col2t[1]);
00667         result[1]*= result[1];
00668         result[2]= (fac1*col1t[2] + fac2*col2t[2]);
00669         result[2]*= result[2];
00670 }
00671 #endif
00672 
00673 static float shade_by_transmission(Isect *is, ShadeInput *shi, ShadeResult *shr)
00674 {
00675         float dx, dy, dz, d, p;
00676 
00677         if (0 == (shi->mat->mode & MA_TRANSP))
00678                 return -1;
00679            
00680         if (shi->mat->tx_limit <= 0.0f) {
00681                 d= 1.0f;
00682         } 
00683         else {
00684                 /* shi.co[] calculated by shade_ray() */
00685                 dx= shi->co[0] - is->start[0];
00686                 dy= shi->co[1] - is->start[1];
00687                 dz= shi->co[2] - is->start[2];
00688                 d= sqrt(dx*dx+dy*dy+dz*dz);
00689                 if (d > shi->mat->tx_limit)
00690                         d= shi->mat->tx_limit;
00691 
00692                 p = shi->mat->tx_falloff;
00693                 if(p < 0.0f) p= 0.0f;
00694                 else if (p > 10.0f) p= 10.0f;
00695 
00696                 shr->alpha *= pow(d, p);
00697                 if (shr->alpha > 1.0f)
00698                         shr->alpha= 1.0f;
00699         }
00700 
00701         return d;
00702 }
00703 
00704 static void ray_fadeout_endcolor(float *col, ShadeInput *origshi, ShadeInput *shi, ShadeResult *shr, Isect *isec, float *vec)
00705 {
00706         /* un-intersected rays get either rendered material color or sky color */
00707         if (origshi->mat->fadeto_mir == MA_RAYMIR_FADETOMAT) {
00708                 VECCOPY(col, shr->combined);
00709         } else if (origshi->mat->fadeto_mir == MA_RAYMIR_FADETOSKY) {
00710                 VECCOPY(shi->view, vec);
00711                 normalize_v3(shi->view);
00712                 
00713                 shadeSkyView(col, isec->start, shi->view, NULL, shi->thread);
00714                 shadeSunView(col, shi->view);
00715         }
00716 }
00717 
00718 static void ray_fadeout(Isect *is, ShadeInput *shi, float *col, float *blendcol, float dist_mir)
00719 {
00720         /* if fading out, linear blend against fade color */
00721         float blendfac;
00722 
00723         blendfac = 1.0 - len_v3v3(shi->co, is->start)/dist_mir;
00724         
00725         col[0] = col[0]*blendfac + (1.0 - blendfac)*blendcol[0];
00726         col[1] = col[1]*blendfac + (1.0 - blendfac)*blendcol[1];
00727         col[2] = col[2]*blendfac + (1.0 - blendfac)*blendcol[2];
00728 }
00729 
00730 /* the main recursive tracer itself
00731  * note: 'col' must be initialized */
00732 static void traceray(ShadeInput *origshi, ShadeResult *origshr, short depth, float *start, float *dir, float *col, ObjectInstanceRen *obi, VlakRen *vlr, int traflag)
00733 {
00734         ShadeInput shi= {0};
00735         Isect isec;
00736         float dist_mir = origshi->mat->dist_mir;
00737         
00738         VECCOPY(isec.start, start);
00739         VECCOPY(isec.dir, dir );
00740         isec.dist = dist_mir > 0 ? dist_mir : RE_RAYTRACE_MAXDIST;
00741         isec.mode= RE_RAY_MIRROR;
00742         isec.check = RE_CHECK_VLR_RENDER;
00743         isec.skip = RE_SKIP_VLR_NEIGHBOUR;
00744         isec.hint = 0;
00745 
00746         isec.orig.ob   = obi;
00747         isec.orig.face = vlr;
00748         RE_RC_INIT(isec, shi);
00749 
00750         if(RE_rayobject_raycast(R.raytree, &isec)) {
00751                 ShadeResult shr= {{0}};
00752                 float d= 1.0f;
00753 
00754                 /* for as long we don't have proper dx/dy transform for rays we copy over original */
00755                 VECCOPY(shi.dxco, origshi->dxco);
00756                 VECCOPY(shi.dyco, origshi->dyco);
00757                 
00758                 shi.mask= origshi->mask;
00759                 shi.osatex= origshi->osatex;
00760                 shi.depth= origshi->depth + 1;                                  /* only used to indicate tracing */
00761                 shi.thread= origshi->thread;
00762                 //shi.sample= 0; // memset above, so dont need this
00763                 shi.xs= origshi->xs;
00764                 shi.ys= origshi->ys;
00765                 shi.lay= origshi->lay;
00766                 shi.passflag= SCE_PASS_COMBINED; /* result of tracing needs no pass info */
00767                 shi.combinedflag= 0xFFFFFF;              /* ray trace does all options */
00768                 //shi.do_preview= 0; // memset above, so dont need this
00769                 shi.light_override= origshi->light_override;
00770                 shi.mat_override= origshi->mat_override;
00771                 
00772                 shade_ray(&isec, &shi, &shr);
00773                 /* ray has traveled inside the material, so shade by transmission */
00774                 if (traflag & RAY_INSIDE)
00775                         d= shade_by_transmission(&isec, &shi, &shr);
00776                 
00777                 if(depth>0) {
00778                         float fr, fg, fb, f, f1;
00779 
00780                         if((shi.mat->mode_l & MA_TRANSP) && shr.alpha < 1.0f && (shi.mat->mode_l & (MA_ZTRANSP | MA_RAYTRANSP))) { 
00781                                 float nf, f, refract[3], tracol[4];
00782                                 
00783                                 tracol[0]= shi.r;
00784                                 tracol[1]= shi.g;
00785                                 tracol[2]= shi.b;
00786                                 tracol[3]= col[3];      // we pass on and accumulate alpha
00787                                 
00788                                 if((shi.mat->mode & MA_TRANSP) && (shi.mat->mode & MA_RAYTRANSP)) {
00789                                         if(traflag & RAY_INSIDE) {
00790                                                 /* inside the material, so use inverse normal */
00791                                                 float norm[3];
00792                                                 norm[0]= - shi.vn[0];
00793                                                 norm[1]= - shi.vn[1];
00794                                                 norm[2]= - shi.vn[2];
00795 
00796                                                 if (refraction(refract, norm, shi.view, shi.ang)) {
00797                                                         /* ray comes out from the material into air */
00798                                                         traflag &= ~RAY_INSIDE;
00799                                                 }
00800                                                 else {
00801                                                         /* total internal reflection (ray stays inside the material) */
00802                                                         reflection(refract, norm, shi.view, shi.vn);
00803                                                 }
00804                                         }
00805                                         else {
00806                                                 if (refraction(refract, shi.vn, shi.view, shi.ang)) {
00807                                                         /* ray goes in to the material from air */
00808                                                         traflag |= RAY_INSIDE;
00809                                                 }
00810                                                 else {
00811                                                         /* total external reflection (ray doesn't enter the material) */
00812                                                         reflection(refract, shi.vn, shi.view, shi.vn);
00813                                                 }
00814                                         }
00815                                         traceray(origshi, origshr, depth-1, shi.co, refract, tracol, shi.obi, shi.vlr, traflag);
00816                                 }
00817                                 else
00818                                         traceray(origshi, origshr, depth-1, shi.co, shi.view, tracol, shi.obi, shi.vlr, 0);
00819                                 
00820                                 f= shr.alpha; f1= 1.0f-f;
00821                                 nf= d * shi.mat->filter;
00822                                 fr= 1.0f+ nf*(shi.r-1.0f);
00823                                 fg= 1.0f+ nf*(shi.g-1.0f);
00824                                 fb= 1.0f+ nf*(shi.b-1.0f);
00825                                 shr.diff[0]= f*shr.diff[0] + f1*fr*tracol[0];
00826                                 shr.diff[1]= f*shr.diff[1] + f1*fg*tracol[1];
00827                                 shr.diff[2]= f*shr.diff[2] + f1*fb*tracol[2];
00828                                 
00829                                 shr.spec[0] *=f;
00830                                 shr.spec[1] *=f;
00831                                 shr.spec[2] *=f;
00832 
00833                                 col[3]= f1*tracol[3] + f;
00834                         }
00835                         else 
00836                                 col[3]= 1.0f;
00837 
00838                         if(shi.mat->mode_l & MA_RAYMIRROR) {
00839                                 f= shi.ray_mirror;
00840                                 if(f!=0.0f) f*= fresnel_fac(shi.view, shi.vn, shi.mat->fresnel_mir_i, shi.mat->fresnel_mir);
00841                         }
00842                         else f= 0.0f;
00843                         
00844                         if(f!=0.0f) {
00845                                 float mircol[4];
00846                                 float ref[3];
00847                                 
00848                                 reflection(ref, shi.vn, shi.view, NULL);                        
00849                                 traceray(origshi, origshr, depth-1, shi.co, ref, mircol, shi.obi, shi.vlr, 0);
00850                         
00851                                 f1= 1.0f-f;
00852 
00853                                 /* combine */
00854                                 //color_combine(col, f*fr*(1.0f-shr.spec[0]), f1, col, shr.diff);
00855                                 //col[0]+= shr.spec[0];
00856                                 //col[1]+= shr.spec[1];
00857                                 //col[2]+= shr.spec[2];
00858                                 
00859                                 fr= shi.mirr;
00860                                 fg= shi.mirg;
00861                                 fb= shi.mirb;
00862                 
00863                                 col[0]= f*fr*(1.0f-shr.spec[0])*mircol[0] + f1*shr.diff[0] + shr.spec[0];
00864                                 col[1]= f*fg*(1.0f-shr.spec[1])*mircol[1] + f1*shr.diff[1] + shr.spec[1];
00865                                 col[2]= f*fb*(1.0f-shr.spec[2])*mircol[2] + f1*shr.diff[2] + shr.spec[2];
00866                         }
00867                         else {
00868                                 col[0]= shr.diff[0] + shr.spec[0];
00869                                 col[1]= shr.diff[1] + shr.spec[1];
00870                                 col[2]= shr.diff[2] + shr.spec[2];
00871                         }
00872                         
00873                         if (dist_mir > 0.0) {
00874                                 float blendcol[3];
00875                                 
00876                                 /* max ray distance set, but found an intersection, so fade this color
00877                                  * out towards the sky/material color for a smooth transition */
00878                                 ray_fadeout_endcolor(blendcol, origshi, &shi, origshr, &isec, dir);
00879                                 ray_fadeout(&isec, &shi, col, blendcol, dist_mir);
00880                         }
00881                 }
00882                 else {
00883                         col[0]= shr.diff[0] + shr.spec[0];
00884                         col[1]= shr.diff[1] + shr.spec[1];
00885                         col[2]= shr.diff[2] + shr.spec[2];
00886                 }
00887                 
00888         }
00889         else {
00890                 ray_fadeout_endcolor(col, origshi, &shi, origshr, &isec, dir);
00891         }
00892         RE_RC_MERGE(&origshi->raycounter, &shi.raycounter);
00893 }
00894 
00895 /* **************** jitter blocks ********** */
00896 
00897 /* calc distributed planar energy */
00898 
00899 static void DP_energy(float *table, float *vec, int tot, float xsize, float ysize)
00900 {
00901         int x, y, a;
00902         float *fp, force[3], result[3];
00903         float dx, dy, dist, min;
00904         
00905         min= MIN2(xsize, ysize);
00906         min*= min;
00907         result[0]= result[1]= 0.0f;
00908         
00909         for(y= -1; y<2; y++) {
00910                 dy= ysize*y;
00911                 for(x= -1; x<2; x++) {
00912                         dx= xsize*x;
00913                         fp= table;
00914                         for(a=0; a<tot; a++, fp+= 2) {
00915                                 force[0]= vec[0] - fp[0]-dx;
00916                                 force[1]= vec[1] - fp[1]-dy;
00917                                 dist= force[0]*force[0] + force[1]*force[1];
00918                                 if(dist < min && dist>0.0f) {
00919                                         result[0]+= force[0]/dist;
00920                                         result[1]+= force[1]/dist;
00921                                 }
00922                         }
00923                 }
00924         }
00925         vec[0] += 0.1*min*result[0]/(float)tot;
00926         vec[1] += 0.1*min*result[1]/(float)tot;
00927         // cyclic clamping
00928         vec[0]= vec[0] - xsize*floor(vec[0]/xsize + 0.5);
00929         vec[1]= vec[1] - ysize*floor(vec[1]/ysize + 0.5);
00930 }
00931 
00932 // random offset of 1 in 2
00933 static void jitter_plane_offset(float *jitter1, float *jitter2, int tot, float sizex, float sizey, float ofsx, float ofsy)
00934 {
00935         float dsizex= sizex*ofsx;
00936         float dsizey= sizey*ofsy;
00937         float hsizex= 0.5*sizex, hsizey= 0.5*sizey;
00938         int x;
00939         
00940         for(x=tot; x>0; x--, jitter1+=2, jitter2+=2) {
00941                 jitter2[0]= jitter1[0] + dsizex;
00942                 jitter2[1]= jitter1[1] + dsizey;
00943                 if(jitter2[0] > hsizex) jitter2[0]-= sizex;
00944                 if(jitter2[1] > hsizey) jitter2[1]-= sizey;
00945         }
00946 }
00947 
00948 /* called from convertBlenderScene.c */
00949 /* we do this in advance to get consistent random, not alter the render seed, and be threadsafe */
00950 void init_jitter_plane(LampRen *lar)
00951 {
00952         float *fp;
00953         int x, iter=12, tot= lar->ray_totsamp;
00954         
00955         /* test if already initialized */
00956         if(lar->jitter) return;
00957         
00958         /* at least 4, or max threads+1 tables */
00959         if(BLENDER_MAX_THREADS < 4) x= 4;
00960         else x= BLENDER_MAX_THREADS+1;
00961         fp= lar->jitter= MEM_callocN(x*tot*2*sizeof(float), "lamp jitter tab");
00962         
00963         /* if 1 sample, we leave table to be zero's */
00964         if(tot>1) {
00965                 
00966                 /* set per-lamp fixed seed */
00967                 BLI_srandom(tot);
00968                 
00969                 /* fill table with random locations, area_size large */
00970                 for(x=0; x<tot; x++, fp+=2) {
00971                         fp[0]= (BLI_frand()-0.5)*lar->area_size;
00972                         fp[1]= (BLI_frand()-0.5)*lar->area_sizey;
00973                 }
00974                 
00975                 while(iter--) {
00976                         fp= lar->jitter;
00977                         for(x=tot; x>0; x--, fp+=2) {
00978                                 DP_energy(lar->jitter, fp, tot, lar->area_size, lar->area_sizey);
00979                         }
00980                 }
00981         }       
00982         /* create the dithered tables (could just check lamp type!) */
00983         jitter_plane_offset(lar->jitter, lar->jitter+2*tot, tot, lar->area_size, lar->area_sizey, 0.5f, 0.0f);
00984         jitter_plane_offset(lar->jitter, lar->jitter+4*tot, tot, lar->area_size, lar->area_sizey, 0.5f, 0.5f);
00985         jitter_plane_offset(lar->jitter, lar->jitter+6*tot, tot, lar->area_size, lar->area_sizey, 0.0f, 0.5f);
00986 }
00987 
00988 /* table around origin, -0.5*size to 0.5*size */
00989 static float *give_jitter_plane(LampRen *lar, int thread, int xs, int ys)
00990 {
00991         int tot;
00992         
00993         tot= lar->ray_totsamp;
00994                         
00995         if(lar->ray_samp_type & LA_SAMP_JITTER) {
00996                 /* made it threadsafe */
00997                 
00998                 if(lar->xold[thread]!=xs || lar->yold[thread]!=ys) {
00999                         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));
01000                         lar->xold[thread]= xs; 
01001                         lar->yold[thread]= ys;
01002                 }
01003                 return lar->jitter+2*(thread+1)*tot;
01004         }
01005         if(lar->ray_samp_type & LA_SAMP_DITHER) {
01006                 return lar->jitter + 2*tot*((xs & 1)+2*(ys & 1));
01007         }
01008         
01009         return lar->jitter;
01010 }
01011 
01012 
01013 /* **************** QMC sampling *************** */
01014 
01015 static void halton_sample(double *ht_invprimes, double *ht_nums, double *v)
01016 {
01017         // incremental halton sequence generator, from:
01018         // "Instant Radiosity", Keller A.
01019         unsigned int i;
01020         
01021         for (i = 0; i < 2; i++)
01022         {
01023                 double r = fabs((1.0 - ht_nums[i]) - 1e-10);
01024                 
01025                 if (ht_invprimes[i] >= r)
01026                 {
01027                         double lasth;
01028                         double h = ht_invprimes[i];
01029                         
01030                         do {
01031                                 lasth = h;
01032                                 h *= ht_invprimes[i];
01033                         } while (h >= r);
01034                         
01035                         ht_nums[i] += ((lasth + h) - 1.0);
01036                 }
01037                 else
01038                         ht_nums[i] += ht_invprimes[i];
01039                 
01040                 v[i] = (float)ht_nums[i];
01041         }
01042 }
01043 
01044 /* Generate Hammersley points in [0,1)^2
01045  * From Lucille renderer */
01046 static void hammersley_create(double *out, int n)
01047 {
01048         double p, t;
01049         int k, kk;
01050 
01051         for (k = 0; k < n; k++) {
01052                 t = 0;
01053                 for (p = 0.5, kk = k; kk; p *= 0.5, kk >>= 1) {
01054                         if (kk & 1) {           /* kk mod 2 = 1         */
01055                                 t += p;
01056                         }
01057                 }
01058         
01059                 out[2 * k + 0] = (double)k / (double)n;
01060                 out[2 * k + 1] = t;
01061         }
01062 }
01063 
01064 static struct QMCSampler *QMC_initSampler(int type, int tot)
01065 {       
01066         QMCSampler *qsa = MEM_callocN(sizeof(QMCSampler), "qmc sampler");
01067         qsa->samp2d = MEM_callocN(2*sizeof(double)*tot, "qmc sample table");
01068 
01069         qsa->tot = tot;
01070         qsa->type = type;
01071         
01072         if (qsa->type==SAMP_TYPE_HAMMERSLEY) 
01073                 hammersley_create(qsa->samp2d, qsa->tot);
01074                 
01075         return qsa;
01076 }
01077 
01078 static void QMC_initPixel(QMCSampler *qsa, int thread)
01079 {
01080         if (qsa->type==SAMP_TYPE_HAMMERSLEY)
01081         {
01082                 /* hammersley sequence is fixed, already created in QMCSampler init.
01083                  * per pixel, gets a random offset. We create separate offsets per thread, for write-safety */
01084                 qsa->offs[thread][0] = 0.5 * BLI_thread_frand(thread);
01085                 qsa->offs[thread][1] = 0.5 * BLI_thread_frand(thread);
01086         }
01087         else {  /* SAMP_TYPE_HALTON */
01088                 
01089                 /* generate a new randomised halton sequence per pixel
01090                  * to alleviate qmc artifacts and make it reproducable 
01091                  * between threads/frames */
01092                 double ht_invprimes[2], ht_nums[2];
01093                 double r[2];
01094                 int i;
01095         
01096                 ht_nums[0] = BLI_thread_frand(thread);
01097                 ht_nums[1] = BLI_thread_frand(thread);
01098                 ht_invprimes[0] = 0.5;
01099                 ht_invprimes[1] = 1.0/3.0;
01100                 
01101                 for (i=0; i< qsa->tot; i++) {
01102                         halton_sample(ht_invprimes, ht_nums, r);
01103                         qsa->samp2d[2*i+0] = r[0];
01104                         qsa->samp2d[2*i+1] = r[1];
01105                 }
01106         }
01107 }
01108 
01109 static void QMC_freeSampler(QMCSampler *qsa)
01110 {
01111         MEM_freeN(qsa->samp2d);
01112         MEM_freeN(qsa);
01113 }
01114 
01115 static void QMC_getSample(double *s, QMCSampler *qsa, int thread, int num)
01116 {
01117         if (qsa->type == SAMP_TYPE_HAMMERSLEY) {
01118                 s[0] = fmod(qsa->samp2d[2*num+0] + qsa->offs[thread][0], 1.0f);
01119                 s[1] = fmod(qsa->samp2d[2*num+1] + qsa->offs[thread][1], 1.0f);
01120         }
01121         else { /* SAMP_TYPE_HALTON */
01122                 s[0] = qsa->samp2d[2*num+0];
01123                 s[1] = qsa->samp2d[2*num+1];
01124         }
01125 }
01126 
01127 /* phong weighted disc using 'blur' for exponent, centred on 0,0 */
01128 static void QMC_samplePhong(float *vec, QMCSampler *qsa, int thread, int num, float blur)
01129 {
01130         double s[2];
01131         float phi, pz, sqr;
01132         
01133         QMC_getSample(s, qsa, thread, num);
01134 
01135         phi = s[0]*2*M_PI;
01136         pz = pow(s[1], blur);
01137         sqr = sqrt(1.0f-pz*pz);
01138 
01139         vec[0] = cos(phi)*sqr;
01140         vec[1] = sin(phi)*sqr;
01141         vec[2] = 0.0f;
01142 }
01143 
01144 /* rect of edge lengths sizex, sizey, centred on 0.0,0.0 i.e. ranging from -sizex/2 to +sizey/2 */
01145 static void QMC_sampleRect(float *vec, QMCSampler *qsa, int thread, int num, float sizex, float sizey)
01146 {
01147         double s[2];
01148 
01149         QMC_getSample(s, qsa, thread, num);
01150                 
01151         vec[0] = (s[0] - 0.5) * sizex;
01152         vec[1] = (s[1] - 0.5) * sizey;
01153         vec[2] = 0.0f;
01154 }
01155 
01156 /* disc of radius 'radius', centred on 0,0 */
01157 static void QMC_sampleDisc(float *vec, QMCSampler *qsa, int thread, int num, float radius)
01158 {
01159         double s[2];
01160         float phi, sqr;
01161         
01162         QMC_getSample(s, qsa, thread, num);
01163         
01164         phi = s[0]*2*M_PI;
01165         sqr = sqrt(s[1]);
01166 
01167         vec[0] = cos(phi)*sqr* radius/2.0;
01168         vec[1] = sin(phi)*sqr* radius/2.0;
01169         vec[2] = 0.0f;
01170 }
01171 
01172 /* uniform hemisphere sampling */
01173 static void QMC_sampleHemi(float *vec, QMCSampler *qsa, int thread, int num)
01174 {
01175         double s[2];
01176         float phi, sqr;
01177         
01178         QMC_getSample(s, qsa, thread, num);
01179         
01180         phi = s[0]*2.f*M_PI;    
01181         sqr = sqrt(s[1]);
01182 
01183         vec[0] = cos(phi)*sqr;
01184         vec[1] = sin(phi)*sqr;
01185         vec[2] = 1.f - s[1]*s[1];
01186 }
01187 
01188 #if 0 /* currently not used */
01189 /* cosine weighted hemisphere sampling */
01190 static void QMC_sampleHemiCosine(float *vec, QMCSampler *qsa, int thread, int num)
01191 {
01192         double s[2];
01193         float phi, sqr;
01194         
01195         QMC_getSample(s, qsa, thread, num);
01196         
01197         phi = s[0]*2.f*M_PI;    
01198         sqr = s[1]*sqrt(2-s[1]*s[1]);
01199 
01200         vec[0] = cos(phi)*sqr;
01201         vec[1] = sin(phi)*sqr;
01202         vec[2] = 1.f - s[1]*s[1];
01203 
01204 }
01205 #endif
01206 
01207 /* called from convertBlenderScene.c */
01208 void init_render_qmcsampler(Render *re)
01209 {
01210         re->qmcsamplers= MEM_callocN(sizeof(ListBase)*BLENDER_MAX_THREADS, "QMCListBase");
01211 }
01212 
01213 static QMCSampler *get_thread_qmcsampler(Render *re, int thread, int type, int tot)
01214 {
01215         QMCSampler *qsa;
01216 
01217         /* create qmc samplers as needed, since recursion makes it hard to
01218          * predict how many are needed */
01219 
01220         for(qsa=re->qmcsamplers[thread].first; qsa; qsa=qsa->next) {
01221                 if(qsa->type == type && qsa->tot == tot && !qsa->used) {
01222                         qsa->used= 1;
01223                         return qsa;
01224                 }
01225         }
01226 
01227         qsa= QMC_initSampler(type, tot);
01228         qsa->used= 1;
01229         BLI_addtail(&re->qmcsamplers[thread], qsa);
01230 
01231         return qsa;
01232 }
01233 
01234 static void release_thread_qmcsampler(Render *UNUSED(re), int UNUSED(thread), QMCSampler *qsa)
01235 {
01236         qsa->used= 0;
01237 }
01238 
01239 void free_render_qmcsampler(Render *re)
01240 {
01241         QMCSampler *qsa, *next;
01242         int a;
01243 
01244         if(re->qmcsamplers) {
01245                 for(a=0; a<BLENDER_MAX_THREADS; a++) {
01246                         for(qsa=re->qmcsamplers[a].first; qsa; qsa=next) {
01247                                 next= qsa->next;
01248                                 QMC_freeSampler(qsa);
01249                         }
01250 
01251                         re->qmcsamplers[a].first= re->qmcsamplers[a].last= NULL;
01252                 }
01253 
01254                 MEM_freeN(re->qmcsamplers);
01255                 re->qmcsamplers= NULL;
01256         }
01257 }
01258 
01259 static int adaptive_sample_variance(int samples, float *col, float *colsq, float thresh)
01260 {
01261         float var[3], mean[3];
01262 
01263         /* scale threshold just to give a bit more precision in input rather than dealing with 
01264          * tiny tiny numbers in the UI */
01265         thresh /= 2;
01266         
01267         mean[0] = col[0] / (float)samples;
01268         mean[1] = col[1] / (float)samples;
01269         mean[2] = col[2] / (float)samples;
01270 
01271         var[0] = (colsq[0] / (float)samples) - (mean[0]*mean[0]);
01272         var[1] = (colsq[1] / (float)samples) - (mean[1]*mean[1]);
01273         var[2] = (colsq[2] / (float)samples) - (mean[2]*mean[2]);
01274         
01275         if ((var[0] * 0.4 < thresh) && (var[1] * 0.3 < thresh) && (var[2] * 0.6 < thresh))
01276                 return 1;
01277         else
01278                 return 0;
01279 }
01280 
01281 static int adaptive_sample_contrast_val(int samples, float prev, float val, float thresh)
01282 {
01283         /* if the last sample's contribution to the total value was below a small threshold
01284          * (i.e. the samples taken are very similar), then taking more samples that are probably 
01285          * going to be the same is wasting effort */
01286         if (fabs( prev/(float)(samples-1) - val/(float)samples ) < thresh) {
01287                 return 1;
01288         } else
01289                 return 0;
01290 }
01291 
01292 static float get_avg_speed(ShadeInput *shi)
01293 {
01294         float pre_x, pre_y, post_x, post_y, speedavg;
01295         
01296         pre_x = (shi->winspeed[0] == PASS_VECTOR_MAX)?0.0:shi->winspeed[0];
01297         pre_y = (shi->winspeed[1] == PASS_VECTOR_MAX)?0.0:shi->winspeed[1];
01298         post_x = (shi->winspeed[2] == PASS_VECTOR_MAX)?0.0:shi->winspeed[2];
01299         post_y = (shi->winspeed[3] == PASS_VECTOR_MAX)?0.0:shi->winspeed[3];
01300         
01301         speedavg = (sqrt(pre_x*pre_x + pre_y*pre_y) + sqrt(post_x*post_x + post_y*post_y)) / 2.0;
01302         
01303         return speedavg;
01304 }
01305 
01306 /* ***************** main calls ************** */
01307 
01308 
01309 static void trace_refract(float *col, ShadeInput *shi, ShadeResult *shr)
01310 {
01311         QMCSampler *qsa=NULL;
01312         int samp_type;
01313         int traflag=0;
01314         
01315         float samp3d[3], orthx[3], orthy[3];
01316         float v_refract[3], v_refract_new[3];
01317         float sampcol[4], colsq[4];
01318         
01319         float blur = pow(1.0 - shi->mat->gloss_tra, 3);
01320         short max_samples = shi->mat->samp_gloss_tra;
01321         float adapt_thresh = shi->mat->adapt_thresh_tra;
01322         
01323         int samples=0;
01324         
01325         colsq[0] = colsq[1] = colsq[2] = 0.0;
01326         col[0] = col[1] = col[2] = 0.0;
01327         col[3]= shr->alpha;
01328         
01329         if (blur > 0.0) {
01330                 if (adapt_thresh != 0.0) samp_type = SAMP_TYPE_HALTON;
01331                 else samp_type = SAMP_TYPE_HAMMERSLEY;
01332                         
01333                 /* all samples are generated per pixel */
01334                 qsa = get_thread_qmcsampler(&R, shi->thread, samp_type, max_samples);
01335                 QMC_initPixel(qsa, shi->thread);
01336         } else 
01337                 max_samples = 1;
01338         
01339 
01340         while (samples < max_samples) {         
01341                 if(refraction(v_refract, shi->vn, shi->view, shi->ang)) {
01342                         traflag |= RAY_INSIDE;
01343                 } else {
01344                         /* total external reflection can happen for materials with IOR < 1.0 */
01345                         if((shi->vlr->flag & R_SMOOTH)) 
01346                                 reflection(v_refract, shi->vn, shi->view, shi->facenor);
01347                         else
01348                                 reflection(v_refract, shi->vn, shi->view, NULL);
01349 
01350                         /* can't blur total external reflection */
01351                         max_samples = 1;
01352                 }
01353                 
01354                 if (max_samples > 1) {
01355                         /* get a quasi-random vector from a phong-weighted disc */
01356                         QMC_samplePhong(samp3d, qsa, shi->thread, samples, blur);
01357                                                 
01358                         ortho_basis_v3v3_v3( orthx, orthy,v_refract);
01359                         mul_v3_fl(orthx, samp3d[0]);
01360                         mul_v3_fl(orthy, samp3d[1]);
01361                                 
01362                         /* and perturb the refraction vector in it */
01363                         add_v3_v3v3(v_refract_new, v_refract, orthx);
01364                         add_v3_v3(v_refract_new, orthy);
01365                         
01366                         normalize_v3(v_refract_new);
01367                 } else {
01368                         /* no blurriness, use the original normal */
01369                         VECCOPY(v_refract_new, v_refract);
01370                 }
01371                 
01372                 sampcol[0]= sampcol[1]= sampcol[2]= sampcol[3]= 0.0f;
01373 
01374                 traceray(shi, shr, shi->mat->ray_depth_tra, shi->co, v_refract_new, sampcol, shi->obi, shi->vlr, traflag);
01375         
01376                 col[0] += sampcol[0];
01377                 col[1] += sampcol[1];
01378                 col[2] += sampcol[2];
01379                 col[3] += sampcol[3];
01380                 
01381                 /* for variance calc */
01382                 colsq[0] += sampcol[0]*sampcol[0];
01383                 colsq[1] += sampcol[1]*sampcol[1];
01384                 colsq[2] += sampcol[2]*sampcol[2];
01385                 
01386                 samples++;
01387                 
01388                 /* adaptive sampling */
01389                 if (adapt_thresh < 1.0 && samples > max_samples/2) 
01390                 {
01391                         if (adaptive_sample_variance(samples, col, colsq, adapt_thresh))
01392                                 break;
01393                         
01394                         /* if the pixel so far is very dark, we can get away with less samples */
01395                         if ( (col[0] + col[1] + col[2])/3.0/(float)samples < 0.01 )
01396                                 max_samples--;
01397                 }
01398         }
01399         
01400         col[0] /= (float)samples;
01401         col[1] /= (float)samples;
01402         col[2] /= (float)samples;
01403         col[3] /= (float)samples;
01404         
01405         if (qsa)
01406                 release_thread_qmcsampler(&R, shi->thread, qsa);
01407 }
01408 
01409 static void trace_reflect(float *col, ShadeInput *shi, ShadeResult *shr, float fresnelfac)
01410 {
01411         QMCSampler *qsa=NULL;
01412         int samp_type;
01413         
01414         float samp3d[3], orthx[3], orthy[3];
01415         float v_nor_new[3], v_reflect[3];
01416         float sampcol[4], colsq[4];
01417                 
01418         float blur = pow(1.0 - shi->mat->gloss_mir, 3);
01419         short max_samples = shi->mat->samp_gloss_mir;
01420         float adapt_thresh = shi->mat->adapt_thresh_mir;
01421         float aniso = 1.0 - shi->mat->aniso_gloss_mir;
01422         
01423         int samples=0;
01424         
01425         col[0] = col[1] = col[2] = 0.0;
01426         colsq[0] = colsq[1] = colsq[2] = 0.0;
01427         
01428         if (blur > 0.0) {
01429                 if (adapt_thresh != 0.0) samp_type = SAMP_TYPE_HALTON;
01430                 else samp_type = SAMP_TYPE_HAMMERSLEY;
01431                         
01432                 /* all samples are generated per pixel */
01433                 qsa = get_thread_qmcsampler(&R, shi->thread, samp_type, max_samples);
01434                 QMC_initPixel(qsa, shi->thread);
01435         } else 
01436                 max_samples = 1;
01437         
01438         while (samples < max_samples) {
01439                                 
01440                 if (max_samples > 1) {
01441                         /* get a quasi-random vector from a phong-weighted disc */
01442                         QMC_samplePhong(samp3d, qsa, shi->thread, samples, blur);
01443 
01444                         /* find the normal's perpendicular plane, blurring along tangents
01445                          * if tangent shading enabled */
01446                         if (shi->mat->mode & (MA_TANGENT_V)) {
01447                                 cross_v3_v3v3(orthx, shi->vn, shi->tang);      // bitangent
01448                                 VECCOPY(orthy, shi->tang);
01449                                 mul_v3_fl(orthx, samp3d[0]);
01450                                 mul_v3_fl(orthy, samp3d[1]*aniso);
01451                         } else {
01452                                 ortho_basis_v3v3_v3( orthx, orthy,shi->vn);
01453                                 mul_v3_fl(orthx, samp3d[0]);
01454                                 mul_v3_fl(orthy, samp3d[1]);
01455                         }
01456 
01457                         /* and perturb the normal in it */
01458                         add_v3_v3v3(v_nor_new, shi->vn, orthx);
01459                         add_v3_v3(v_nor_new, orthy);
01460                         normalize_v3(v_nor_new);
01461                 } else {
01462                         /* no blurriness, use the original normal */
01463                         VECCOPY(v_nor_new, shi->vn);
01464                 }
01465                 
01466                 if((shi->vlr->flag & R_SMOOTH)) 
01467                         reflection(v_reflect, v_nor_new, shi->view, shi->facenor);
01468                 else
01469                         reflection(v_reflect, v_nor_new, shi->view, NULL);
01470                 
01471                 sampcol[0]= sampcol[1]= sampcol[2]= sampcol[3]= 0.0f;
01472 
01473                 traceray(shi, shr, shi->mat->ray_depth, shi->co, v_reflect, sampcol, shi->obi, shi->vlr, 0);
01474 
01475                 
01476                 col[0] += sampcol[0];
01477                 col[1] += sampcol[1];
01478                 col[2] += sampcol[2];
01479         
01480                 /* for variance calc */
01481                 colsq[0] += sampcol[0]*sampcol[0];
01482                 colsq[1] += sampcol[1]*sampcol[1];
01483                 colsq[2] += sampcol[2]*sampcol[2];
01484                 
01485                 samples++;
01486 
01487                 /* adaptive sampling */
01488                 if (adapt_thresh > 0.0 && samples > max_samples/3) 
01489                 {
01490                         if (adaptive_sample_variance(samples, col, colsq, adapt_thresh))
01491                                 break;
01492                         
01493                         /* if the pixel so far is very dark, we can get away with less samples */
01494                         if ( (col[0] + col[1] + col[2])/3.0/(float)samples < 0.01 )
01495                                 max_samples--;
01496                 
01497                         /* reduce samples when reflection is dim due to low ray mirror blend value or fresnel factor
01498                          * and when reflection is blurry */
01499                         if (fresnelfac < 0.1 * (blur+1)) {
01500                                 max_samples--;
01501                                 
01502                                 /* even more for very dim */
01503                                 if (fresnelfac < 0.05 * (blur+1)) 
01504                                         max_samples--;
01505                         }
01506                 }
01507         }
01508         
01509         col[0] /= (float)samples;
01510         col[1] /= (float)samples;
01511         col[2] /= (float)samples;
01512         
01513         if (qsa)
01514                 release_thread_qmcsampler(&R, shi->thread, qsa);
01515 }
01516 
01517 /* extern call from render loop */
01518 void ray_trace(ShadeInput *shi, ShadeResult *shr)
01519 {
01520         float i, f, f1, fr, fg, fb;
01521         float mircol[4], tracol[4];
01522         float diff[3];
01523         int do_tra, do_mir;
01524         
01525         do_tra= ((shi->mode & MA_TRANSP) && (shi->mode & MA_RAYTRANSP) && shr->alpha!=1.0f && (shi->depth <= shi->mat->ray_depth_tra));
01526         do_mir= ((shi->mat->mode & MA_RAYMIRROR) && shi->ray_mirror!=0.0f && (shi->depth <= shi->mat->ray_depth));
01527         
01528         /* raytrace mirror amd refract like to separate the spec color */
01529         if(shi->combinedflag & SCE_PASS_SPEC)
01530                 VECSUB(diff, shr->combined, shr->spec) /* no ; */
01531         else
01532                 VECCOPY(diff, shr->combined);
01533         
01534         if(do_tra) {
01535                 float olddiff[3];
01536                 
01537                 trace_refract(tracol, shi, shr);
01538                 
01539                 f= shr->alpha; f1= 1.0f-f;
01540                 fr= 1.0f+ shi->mat->filter*(shi->r-1.0f);
01541                 fg= 1.0f+ shi->mat->filter*(shi->g-1.0f);
01542                 fb= 1.0f+ shi->mat->filter*(shi->b-1.0f);
01543                 
01544                 /* for refract pass */
01545                 VECCOPY(olddiff, diff);
01546                 
01547                 diff[0]= f*diff[0] + f1*fr*tracol[0];
01548                 diff[1]= f*diff[1] + f1*fg*tracol[1];
01549                 diff[2]= f*diff[2] + f1*fb*tracol[2];
01550                 
01551                 if(shi->passflag & SCE_PASS_REFRACT)
01552                         VECSUB(shr->refr, diff, olddiff);
01553                 
01554                 if(!(shi->combinedflag & SCE_PASS_REFRACT))
01555                         VECSUB(diff, diff, shr->refr);
01556                 
01557                 shr->alpha= MIN2(1.0f, tracol[3]);
01558         }
01559         
01560         if(do_mir) {
01561         
01562                 i= shi->ray_mirror*fresnel_fac(shi->view, shi->vn, shi->mat->fresnel_mir_i, shi->mat->fresnel_mir);
01563                 if(i!=0.0f) {
01564                 
01565                         trace_reflect(mircol, shi, shr, i);
01566                         
01567                         fr= i*shi->mirr;
01568                         fg= i*shi->mirg;
01569                         fb= i*shi->mirb;
01570 
01571                         if(shi->passflag & SCE_PASS_REFLECT) {
01572                                 /* mirror pass is not blocked out with spec */
01573                                 shr->refl[0]= fr*mircol[0] - fr*diff[0];
01574                                 shr->refl[1]= fg*mircol[1] - fg*diff[1];
01575                                 shr->refl[2]= fb*mircol[2] - fb*diff[2];
01576                         }
01577                         
01578                         if(shi->combinedflag & SCE_PASS_REFLECT) {
01579                                 /* values in shr->spec can be greater then 1.0.
01580                                  * In this case the mircol uses a zero blending factor, so ignoring it is ok.
01581                                  * Fixes bug #18837 - when the spec is higher then 1.0,
01582                                  * diff can become a negative color - Campbell  */
01583                                 
01584                                 f1= 1.0f-i;
01585                                 
01586                                 diff[0] *= f1;
01587                                 diff[1] *= f1;
01588                                 diff[2] *= f1;
01589                                 
01590                                 if(shr->spec[0]<1.0f)   diff[0] += mircol[0] * (fr*(1.0f-shr->spec[0]));
01591                                 if(shr->spec[1]<1.0f)   diff[1] += mircol[1] * (fg*(1.0f-shr->spec[1]));
01592                                 if(shr->spec[2]<1.0f)   diff[2] += mircol[2] * (fb*(1.0f-shr->spec[2]));
01593                         }
01594                 }
01595         }
01596         /* put back together */
01597         if(shi->combinedflag & SCE_PASS_SPEC)
01598                 VECADD(shr->combined, diff, shr->spec) /* no ; */
01599         else
01600                 VECCOPY(shr->combined, diff);
01601 }
01602 
01603 /* color 'shadfac' passes through 'col' with alpha and filter */
01604 /* filter is only applied on alpha defined transparent part */
01605 static void addAlphaLight(float *shadfac, float *col, float alpha, float filter)
01606 {
01607         float fr, fg, fb;
01608         
01609         fr= 1.0f+ filter*(col[0]-1.0f);
01610         fg= 1.0f+ filter*(col[1]-1.0f);
01611         fb= 1.0f+ filter*(col[2]-1.0f);
01612         
01613         shadfac[0]= alpha*col[0] + fr*(1.0f-alpha)*shadfac[0];
01614         shadfac[1]= alpha*col[1] + fg*(1.0f-alpha)*shadfac[1];
01615         shadfac[2]= alpha*col[2] + fb*(1.0f-alpha)*shadfac[2];
01616         
01617         shadfac[3]= (1.0f-alpha)*shadfac[3];
01618 }
01619 
01620 static void ray_trace_shadow_tra(Isect *is, ShadeInput *origshi, int depth, int traflag, float col[4])
01621 {
01622         /* ray to lamp, find first face that intersects, check alpha properties,
01623            if it has col[3]>0.0f  continue. so exit when alpha is full */
01624         ShadeInput shi;
01625         ShadeResult shr;
01626         float initial_dist = is->dist;
01627         
01628         if(RE_rayobject_raycast(R.raytree, is)) {
01629                 float d= 1.0f;
01630                 /* we got a face */
01631                 
01632                 /* Warning, This is not that nice, and possibly a bit slow for every ray,
01633                 however some variables were not initialized properly in, unless using shade_input_initialize(...), we need to do a memset */
01634                 memset(&shi, 0, sizeof(ShadeInput)); 
01635                 /* end warning! - Campbell */
01636                 
01637                 shi.depth= origshi->depth + 1;                                  /* only used to indicate tracing */
01638                 shi.mask= origshi->mask;
01639                 shi.thread= origshi->thread;
01640                 shi.passflag= SCE_PASS_COMBINED;
01641                 shi.combinedflag= 0xFFFFFF;              /* ray trace does all options */
01642         
01643                 shi.xs= origshi->xs;
01644                 shi.ys= origshi->ys;
01645                 shi.lay= origshi->lay;
01646                 shi.nodes= origshi->nodes;
01647                 
01648                 shade_ray(is, &shi, &shr);
01649                 if (shi.mat->material_type == MA_TYPE_SURFACE) {
01650                         if (traflag & RAY_TRA)
01651                                 d= shade_by_transmission(is, &shi, &shr);
01652                         
01653                         /* mix colors based on shadfac (rgb + amount of light factor) */
01654                         addAlphaLight(col, shr.diff, shr.alpha, d*shi.mat->filter);
01655                 } else if (shi.mat->material_type == MA_TYPE_VOLUME) {
01656                         const float a = col[3];
01657                         
01658                         col[0] = a*col[0] + shr.alpha*shr.combined[0];
01659                         col[1] = a*col[1] + shr.alpha*shr.combined[1];
01660                         col[2] = a*col[2] + shr.alpha*shr.combined[2];
01661                         
01662                         col[3] = (1.0 - shr.alpha)*a;
01663                 }
01664                 
01665                 if(depth>0 && col[3]>0.0f) {
01666                         
01667                         /* adapt isect struct */
01668                         VECCOPY(is->start, shi.co);
01669                         is->dist = initial_dist-is->dist;
01670                         is->orig.ob   = shi.obi;
01671                         is->orig.face = shi.vlr;
01672 
01673                         ray_trace_shadow_tra(is, origshi, depth-1, traflag | RAY_TRA, col);
01674                 }
01675                 
01676                 RE_RC_MERGE(&origshi->raycounter, &shi.raycounter);
01677         }
01678 }
01679 
01680 /* not used, test function for ambient occlusion (yaf: pathlight) */
01681 /* main problem; has to be called within shading loop, giving unwanted recursion */
01682 static int ray_trace_shadow_rad(ShadeInput *ship, ShadeResult *shr)
01683 {
01684         static int counter=0, only_one= 0;
01685         extern float hashvectf[];
01686         Isect isec;
01687         ShadeInput shi;
01688         ShadeResult shr_t;
01689         float vec[3], accum[3], div= 0.0f;
01690         int a;
01691         
01692         assert(0);
01693         
01694         if(only_one) {
01695                 return 0;
01696         }
01697         only_one= 1;
01698         
01699         accum[0]= accum[1]= accum[2]= 0.0f;
01700         isec.mode= RE_RAY_MIRROR;
01701         isec.orig.ob   = ship->obi;
01702         isec.orig.face = ship->vlr;
01703         isec.hint = 0;
01704 
01705         VECCOPY(isec.start, ship->co);
01706         
01707         RE_RC_INIT(isec, shi);
01708         
01709         for(a=0; a<8*8; a++) {
01710                 
01711                 counter+=3;
01712                 counter %= 768;
01713                 VECCOPY(vec, hashvectf+counter);
01714                 if(ship->vn[0]*vec[0]+ship->vn[1]*vec[1]+ship->vn[2]*vec[2]>0.0f) {
01715                         vec[0]-= vec[0];
01716                         vec[1]-= vec[1];
01717                         vec[2]-= vec[2];
01718                 }
01719 
01720                 VECCOPY(isec.dir, vec );
01721                 isec.dist = RE_RAYTRACE_MAXDIST;
01722 
01723                 if(RE_rayobject_raycast(R.raytree, &isec)) {
01724                         float fac;
01725                         
01726                         /* Warning, This is not that nice, and possibly a bit slow for every ray,
01727                         however some variables were not initialized properly in, unless using shade_input_initialize(...), we need to do a memset */
01728                         memset(&shi, 0, sizeof(ShadeInput)); 
01729                         /* end warning! - Campbell */
01730                         
01731                         shade_ray(&isec, &shi, &shr_t);
01732                         /* fac= isec.dist*isec.dist; */
01733                         fac= 1.0f;
01734                         accum[0]+= fac*(shr_t.diff[0]+shr_t.spec[0]);
01735                         accum[1]+= fac*(shr_t.diff[1]+shr_t.spec[1]);
01736                         accum[2]+= fac*(shr_t.diff[2]+shr_t.spec[2]);
01737                         div+= fac;
01738                 }
01739                 else div+= 1.0f;
01740         }
01741         
01742         if(div!=0.0f) {
01743                 shr->diff[0]+= accum[0]/div;
01744                 shr->diff[1]+= accum[1]/div;
01745                 shr->diff[2]+= accum[2]/div;
01746         }
01747         shr->alpha= 1.0f;
01748         
01749         only_one= 0;
01750         return 1;
01751 }
01752 
01753 /* aolight: function to create random unit sphere vectors for total random sampling */
01754 static void RandomSpherical(float *v)
01755 {
01756         float r;
01757         v[2] = 2.f*BLI_frand()-1.f;
01758         if ((r = 1.f - v[2]*v[2])>0.f) {
01759                 float a = 6.283185307f*BLI_frand();
01760                 r = sqrt(r);
01761                 v[0] = r * cos(a);
01762                 v[1] = r * sin(a);
01763         }
01764         else v[2] = 1.f;
01765 }
01766 
01767 /* calc distributed spherical energy */
01768 static void DS_energy(float *sphere, int tot, float *vec)
01769 {
01770         float *fp, fac, force[3], res[3];
01771         int a;
01772         
01773         res[0]= res[1]= res[2]= 0.0f;
01774         
01775         for(a=0, fp=sphere; a<tot; a++, fp+=3) {
01776                 sub_v3_v3v3(force, vec, fp);
01777                 fac= force[0]*force[0] + force[1]*force[1] + force[2]*force[2];
01778                 if(fac!=0.0f) {
01779                         fac= 1.0f/fac;
01780                         res[0]+= fac*force[0];
01781                         res[1]+= fac*force[1];
01782                         res[2]+= fac*force[2];
01783                 }
01784         }
01785 
01786         mul_v3_fl(res, 0.5);
01787         add_v3_v3(vec, res);
01788         normalize_v3(vec);
01789         
01790 }
01791 
01792 /* called from convertBlenderScene.c */
01793 /* creates an equally distributed spherical sample pattern */
01794 /* and allocates threadsafe memory */
01795 void init_ao_sphere(World *wrld)
01796 {
01797         float *fp;
01798         int a, tot, iter= 16;
01799 
01800         /* we make twice the amount of samples, because only a hemisphere is used */
01801         tot= 2*wrld->aosamp*wrld->aosamp;
01802         
01803         wrld->aosphere= MEM_mallocN(3*tot*sizeof(float), "AO sphere");
01804         
01805         /* fixed random */
01806         BLI_srandom(tot);
01807         
01808         /* init */
01809         fp= wrld->aosphere;
01810         for(a=0; a<tot; a++, fp+= 3) {
01811                 RandomSpherical(fp);
01812         }
01813         
01814         while(iter--) {
01815                 for(a=0, fp= wrld->aosphere; a<tot; a++, fp+= 3) {
01816                         DS_energy(wrld->aosphere, tot, fp);
01817                 }
01818         }
01819         
01820         /* tables */
01821         wrld->aotables= MEM_mallocN(BLENDER_MAX_THREADS*3*tot*sizeof(float), "AO tables");
01822 }
01823 
01824 /* give per thread a table, we have to compare xs ys because of way OSA works... */
01825 static float *threadsafe_table_sphere(int test, int thread, int xs, int ys, int tot)
01826 {
01827         static int xso[BLENDER_MAX_THREADS], yso[BLENDER_MAX_THREADS];
01828         static int firsttime= 1;
01829         
01830         if(firsttime) {
01831                 memset(xso, 255, sizeof(xso));
01832                 memset(yso, 255, sizeof(yso));
01833                 firsttime= 0;
01834         }
01835         
01836         if(xs==xso[thread] && ys==yso[thread]) return R.wrld.aotables+ thread*tot*3;
01837         if(test) return NULL;
01838         xso[thread]= xs; yso[thread]= ys;
01839         return R.wrld.aotables+ thread*tot*3;
01840 }
01841 
01842 static float *sphere_sampler(int type, int resol, int thread, int xs, int ys, int reset)
01843 {
01844         int tot;
01845         float *vec;
01846         
01847         tot= 2*resol*resol;
01848 
01849         if (type & WO_AORNDSMP) {
01850                 float *sphere;
01851                 int a;
01852                 
01853                 // always returns table
01854                 sphere= threadsafe_table_sphere(0, thread, xs, ys, tot);
01855 
01856                 /* total random sampling. NOT THREADSAFE! (should be removed, is not useful) */
01857                 vec= sphere;
01858                 for (a=0; a<tot; a++, vec+=3) {
01859                         RandomSpherical(vec);
01860                 }
01861                 
01862                 return sphere;
01863         } 
01864         else {
01865                 float *sphere;
01866                 float cosfi, sinfi, cost, sint;
01867                 float ang, *vec1;
01868                 int a;
01869                 
01870                 // returns table if xs and ys were equal to last call, and not resetting
01871                 sphere= (reset)? NULL: threadsafe_table_sphere(1, thread, xs, ys, tot);
01872                 if(sphere==NULL) {
01873                         sphere= threadsafe_table_sphere(0, thread, xs, ys, tot);
01874                         
01875                         // random rotation
01876                         ang= BLI_thread_frand(thread);
01877                         sinfi= sin(ang); cosfi= cos(ang);
01878                         ang= BLI_thread_frand(thread);
01879                         sint= sin(ang); cost= cos(ang);
01880                         
01881                         vec= R.wrld.aosphere;
01882                         vec1= sphere;
01883                         for (a=0; a<tot; a++, vec+=3, vec1+=3) {
01884                                 vec1[0]= cost*cosfi*vec[0] - sinfi*vec[1] + sint*cosfi*vec[2];
01885                                 vec1[1]= cost*sinfi*vec[0] + cosfi*vec[1] + sint*sinfi*vec[2];
01886                                 vec1[2]= -sint*vec[0] + cost*vec[2];                    
01887                         }
01888                 }
01889                 return sphere;
01890         }
01891 }
01892 
01893 static void ray_ao_qmc(ShadeInput *shi, float *ao, float *env)
01894 {
01895         Isect isec;
01896         RayHint point_hint;
01897         QMCSampler *qsa=NULL;
01898         float samp3d[3];
01899         float up[3], side[3], dir[3], nrm[3];
01900         
01901         float maxdist = R.wrld.aodist;
01902         float fac=0.0f, prev=0.0f;
01903         float adapt_thresh = R.wrld.ao_adapt_thresh;
01904         float adapt_speed_fac = R.wrld.ao_adapt_speed_fac;
01905         
01906         int samples=0;
01907         int max_samples = R.wrld.aosamp*R.wrld.aosamp;
01908         
01909         float dxyview[3], skyadded=0;
01910         int envcolor;
01911         
01912         RE_RC_INIT(isec, *shi);
01913         isec.orig.ob   = shi->obi;
01914         isec.orig.face = shi->vlr;
01915         isec.check = RE_CHECK_VLR_NON_SOLID_MATERIAL;
01916         isec.skip = RE_SKIP_VLR_NEIGHBOUR;
01917         isec.hint = 0;
01918 
01919         isec.hit.ob   = 0;
01920         isec.hit.face = 0;
01921 
01922         isec.last_hit = NULL;
01923         
01924         isec.mode= (R.wrld.aomode & WO_AODIST)?RE_RAY_SHADOW_TRA:RE_RAY_SHADOW;
01925         isec.lay= -1;
01926         
01927         VECCOPY(isec.start, shi->co);           
01928         RE_rayobject_hint_bb( R.raytree, &point_hint, isec.start, isec.start );
01929         isec.hint = &point_hint;
01930 
01931         zero_v3(ao);
01932         zero_v3(env);
01933         
01934         /* prevent sky colors to be added for only shadow (shadow becomes alpha) */
01935         envcolor= R.wrld.aocolor;
01936         if(shi->mat->mode & MA_ONLYSHADOW)
01937                 envcolor= WO_AOPLAIN;
01938         
01939         if(envcolor == WO_AOSKYTEX) {
01940                 dxyview[0]= 1.0f/(float)R.wrld.aosamp;
01941                 dxyview[1]= 1.0f/(float)R.wrld.aosamp;
01942                 dxyview[2]= 0.0f;
01943         }
01944         
01945         if(shi->vlr->flag & R_SMOOTH) {
01946                 VECCOPY(nrm, shi->vn);
01947         }
01948         else {
01949                 VECCOPY(nrm, shi->facenor);
01950         }
01951 
01952         ortho_basis_v3v3_v3( up, side,nrm);
01953         
01954         /* sampling init */
01955         if (R.wrld.ao_samp_method==WO_AOSAMP_HALTON) {
01956                 float speedfac;
01957                 
01958                 speedfac = get_avg_speed(shi) * adapt_speed_fac;
01959                 CLAMP(speedfac, 1.0, 1000.0);
01960                 max_samples /= speedfac;
01961                 if (max_samples < 5) max_samples = 5;
01962                 
01963                 qsa = get_thread_qmcsampler(&R, shi->thread, SAMP_TYPE_HALTON, max_samples);
01964         } else if (R.wrld.ao_samp_method==WO_AOSAMP_HAMMERSLEY)
01965                 qsa = get_thread_qmcsampler(&R, shi->thread, SAMP_TYPE_HAMMERSLEY, max_samples);
01966 
01967         QMC_initPixel(qsa, shi->thread);
01968         
01969         while (samples < max_samples) {
01970 
01971                 /* sampling, returns quasi-random vector in unit hemisphere */
01972                 QMC_sampleHemi(samp3d, qsa, shi->thread, samples);
01973 
01974                 dir[0] = (samp3d[0]*up[0] + samp3d[1]*side[0] + samp3d[2]*nrm[0]);
01975                 dir[1] = (samp3d[0]*up[1] + samp3d[1]*side[1] + samp3d[2]*nrm[1]);
01976                 dir[2] = (samp3d[0]*up[2] + samp3d[1]*side[2] + samp3d[2]*nrm[2]);
01977                 
01978                 normalize_v3(dir);
01979                         
01980                 isec.dir[0] = -dir[0];
01981                 isec.dir[1] = -dir[1];
01982                 isec.dir[2] = -dir[2];
01983                 isec.dist = maxdist;
01984                 
01985                 prev = fac;
01986                 
01987                 if(RE_rayobject_raycast(R.raytree, &isec)) {
01988                         if (R.wrld.aomode & WO_AODIST) fac+= exp(-isec.dist*R.wrld.aodistfac); 
01989                         else fac+= 1.0f;
01990                 }
01991                 else if(envcolor!=WO_AOPLAIN) {
01992                         float skycol[4];
01993                         float skyfac, view[3];
01994                         
01995                         view[0]= -dir[0];
01996                         view[1]= -dir[1];
01997                         view[2]= -dir[2];
01998                         normalize_v3(view);
01999                         
02000                         if(envcolor==WO_AOSKYCOL) {
02001                                 skyfac= 0.5*(1.0f+view[0]*R.grvec[0]+ view[1]*R.grvec[1]+ view[2]*R.grvec[2]);
02002                                 env[0]+= (1.0f-skyfac)*R.wrld.horr + skyfac*R.wrld.zenr;
02003                                 env[1]+= (1.0f-skyfac)*R.wrld.horg + skyfac*R.wrld.zeng;
02004                                 env[2]+= (1.0f-skyfac)*R.wrld.horb + skyfac*R.wrld.zenb;
02005                         }
02006                         else {  /* WO_AOSKYTEX */
02007                                 shadeSkyView(skycol, isec.start, view, dxyview, shi->thread);
02008                                 shadeSunView(skycol, shi->view);
02009                                 env[0]+= skycol[0];
02010                                 env[1]+= skycol[1];
02011                                 env[2]+= skycol[2];
02012                         }
02013                         skyadded++;
02014                 }
02015                 
02016                 samples++;
02017                 
02018                 if (qsa->type == SAMP_TYPE_HALTON) {
02019                         /* adaptive sampling - consider samples below threshold as in shadow (or vice versa) and exit early */          
02020                         if (adapt_thresh > 0.0 && (samples > max_samples/2) ) {
02021                                 
02022                                 if (adaptive_sample_contrast_val(samples, prev, fac, adapt_thresh)) {
02023                                         break;
02024                                 }
02025                         }
02026                 }
02027         }
02028         
02029         /* average color times distances/hits formula */
02030         ao[0]= ao[1]= ao[2]= 1.0f - fac/(float)samples;
02031 
02032         if(envcolor!=WO_AOPLAIN && skyadded)
02033                 mul_v3_fl(env, (1.0f - fac/(float)samples)/((float)skyadded));
02034         else
02035                 copy_v3_v3(env, ao);
02036         
02037         if (qsa)
02038                 release_thread_qmcsampler(&R, shi->thread, qsa);
02039 }
02040 
02041 /* extern call from shade_lamp_loop, ambient occlusion calculus */
02042 static void ray_ao_spheresamp(ShadeInput *shi, float *ao, float *env)
02043 {
02044         Isect isec;
02045         RayHint point_hint;
02046         float *vec, *nrm, bias, sh=0.0f;
02047         float maxdist = R.wrld.aodist;
02048         float dxyview[3];
02049         int j= -1, tot, actual=0, skyadded=0, envcolor, resol= R.wrld.aosamp;
02050         
02051         RE_RC_INIT(isec, *shi);
02052         isec.orig.ob   = shi->obi;
02053         isec.orig.face = shi->vlr;
02054         isec.check = RE_CHECK_VLR_RENDER;
02055         isec.skip = RE_SKIP_VLR_NEIGHBOUR;
02056         isec.hint = 0;
02057 
02058         isec.hit.ob   = 0;
02059         isec.hit.face = 0;
02060         
02061         isec.last_hit = NULL;
02062         
02063         isec.mode= (R.wrld.aomode & WO_AODIST)?RE_RAY_SHADOW_TRA:RE_RAY_SHADOW;
02064         isec.lay= -1;
02065 
02066         VECCOPY(isec.start, shi->co);           
02067         RE_rayobject_hint_bb( R.raytree, &point_hint, isec.start, isec.start );
02068         isec.hint = &point_hint;
02069 
02070         zero_v3(ao);
02071         zero_v3(env);
02072 
02073         /* bias prevents smoothed faces to appear flat */
02074         if(shi->vlr->flag & R_SMOOTH) {
02075                 bias= R.wrld.aobias;
02076                 nrm= shi->vn;
02077         }
02078         else {
02079                 bias= 0.0f;
02080                 nrm= shi->facenor;
02081         }
02082 
02083         /* prevent sky colors to be added for only shadow (shadow becomes alpha) */
02084         envcolor= R.wrld.aocolor;
02085         if(shi->mat->mode & MA_ONLYSHADOW)
02086                 envcolor= WO_AOPLAIN;
02087         
02088         if(resol>32) resol= 32;
02089 
02090         /* get sphere samples. for faces we get the same samples for sample x/y values,
02091            for strand render we always require a new sampler because x/y are not set */
02092         vec= sphere_sampler(R.wrld.aomode, resol, shi->thread, shi->xs, shi->ys, shi->strand != NULL);
02093         
02094         // warning: since we use full sphere now, and dotproduct is below, we do twice as much
02095         tot= 2*resol*resol;
02096 
02097         if(envcolor == WO_AOSKYTEX) {
02098                 dxyview[0]= 1.0f/(float)resol;
02099                 dxyview[1]= 1.0f/(float)resol;
02100                 dxyview[2]= 0.0f;
02101         }
02102         
02103         while(tot--) {
02104                 
02105                 if ((vec[0]*nrm[0] + vec[1]*nrm[1] + vec[2]*nrm[2]) > bias) {
02106                         /* only ao samples for mask */
02107                         if(R.r.mode & R_OSA) {
02108                                 j++;
02109                                 if(j==R.osa) j= 0;
02110                                 if(!(shi->mask & (1<<j))) {
02111                                         vec+=3;
02112                                         continue;
02113                                 }
02114                         }
02115                         
02116                         actual++;
02117                         
02118                         /* always set start/vec/dist */
02119                         isec.dir[0] = -vec[0];
02120                         isec.dir[1] = -vec[1];
02121                         isec.dir[2] = -vec[2];
02122                         isec.dist = maxdist;
02123                         
02124                         /* do the trace */
02125                         if(RE_rayobject_raycast(R.raytree, &isec)) {
02126                                 if (R.wrld.aomode & WO_AODIST) sh+= exp(-isec.dist*R.wrld.aodistfac); 
02127                                 else sh+= 1.0f;
02128                         }
02129                         else if(envcolor!=WO_AOPLAIN) {
02130                                 float skycol[4];
02131                                 float fac, view[3];
02132                                 
02133                                 view[0]= -vec[0];
02134                                 view[1]= -vec[1];
02135                                 view[2]= -vec[2];
02136                                 normalize_v3(view);
02137                                 
02138                                 if(envcolor==WO_AOSKYCOL) {
02139                                         fac= 0.5*(1.0f+view[0]*R.grvec[0]+ view[1]*R.grvec[1]+ view[2]*R.grvec[2]);
02140                                         env[0]+= (1.0f-fac)*R.wrld.horr + fac*R.wrld.zenr;
02141                                         env[1]+= (1.0f-fac)*R.wrld.horg + fac*R.wrld.zeng;
02142                                         env[2]+= (1.0f-fac)*R.wrld.horb + fac*R.wrld.zenb;
02143                                 }
02144                                 else {  /* WO_AOSKYTEX */
02145                                         shadeSkyView(skycol, isec.start, view, dxyview, shi->thread);
02146                                         shadeSunView(skycol, shi->view);
02147                                         env[0]+= skycol[0];
02148                                         env[1]+= skycol[1];
02149                                         env[2]+= skycol[2];
02150                                 }
02151                                 skyadded++;
02152                         }
02153                 }
02154                 // samples
02155                 vec+= 3;
02156         }
02157         
02158         if(actual==0) sh= 1.0f;
02159         else sh = 1.0f - sh/((float)actual);
02160         
02161         /* average color times distances/hits formula */
02162         ao[0]= ao[1]= ao[2]= sh;
02163 
02164         if(envcolor!=WO_AOPLAIN && skyadded)
02165                 mul_v3_fl(env, sh/((float)skyadded));
02166         else
02167                 copy_v3_v3(env, ao);
02168 }
02169 
02170 void ray_ao(ShadeInput *shi, float *ao, float *env)
02171 {
02172         /* Unfortunately, the unusual way that the sphere sampler calculates roughly twice as many
02173          * samples as are actually traced, and skips them based on bias and OSA settings makes it very difficult
02174          * to reuse code between these two functions. This is the easiest way I can think of to do it
02175          * --broken */
02176         if (ELEM(R.wrld.ao_samp_method, WO_AOSAMP_HAMMERSLEY, WO_AOSAMP_HALTON))
02177                 ray_ao_qmc(shi, ao, env);
02178         else if (R.wrld.ao_samp_method == WO_AOSAMP_CONSTANT)
02179                 ray_ao_spheresamp(shi, ao, env);
02180 }
02181 
02182 static void ray_shadow_jittered_coords(ShadeInput *shi, int max, float jitco[RE_MAX_OSA][3], int *totjitco)
02183 {
02184         /* magic numbers for reordering sample positions to give better
02185          * results with adaptive sample, when it usually only takes 4 samples */
02186         int order8[8] = {0, 1, 5, 6, 2, 3, 4, 7};
02187         int order11[11] = {1, 3, 8, 10, 0, 2, 4, 5, 6, 7, 9};
02188         int order16[16] = {1, 3, 9, 12, 0, 6, 7, 8, 13, 2, 4, 5, 10, 11, 14, 15};
02189         int count = count_mask(shi->mask);
02190 
02191         /* for better antialising shadow samples are distributed over the subpixel
02192          * sample coordinates, this only works for raytracing depth 0 though */
02193         if(!shi->strand && shi->depth == 0 && count > 1 && count <= max) {
02194                 float xs, ys, zs, view[3];
02195                 int samp, ordsamp, tot= 0;
02196 
02197                 for(samp=0; samp<R.osa; samp++) {
02198                         if(R.osa == 8) ordsamp = order8[samp];
02199                         else if(R.osa == 11) ordsamp = order11[samp];
02200                         else if(R.osa == 16) ordsamp = order16[samp];
02201                         else ordsamp = samp;
02202 
02203                         if(shi->mask & (1<<ordsamp)) {
02204                                 /* zbuffer has this inverse corrected, ensures xs,ys are inside pixel */
02205                                 xs= (float)shi->scanco[0] + R.jit[ordsamp][0] + 0.5f;
02206                                 ys= (float)shi->scanco[1] + R.jit[ordsamp][1] + 0.5f;
02207                                 zs= shi->scanco[2];
02208 
02209                                 shade_input_calc_viewco(shi, xs, ys, zs, view, NULL, jitco[tot], NULL, NULL);
02210                                 tot++;
02211                         }
02212                 }
02213 
02214                 *totjitco= tot;
02215         }
02216         else {
02217                 VECCOPY(jitco[0], shi->co);
02218                 *totjitco= 1;
02219         }
02220 }
02221 
02222 static void ray_shadow_qmc(ShadeInput *shi, LampRen *lar, float *lampco, float *shadfac, Isect *isec)
02223 {
02224         QMCSampler *qsa=NULL;
02225         int samples=0;
02226         float samp3d[3];
02227 
02228         float fac=0.0f, vec[3], end[3];
02229         float colsq[4];
02230         float adapt_thresh = lar->adapt_thresh;
02231         int min_adapt_samples=4, max_samples = lar->ray_totsamp;
02232         float *co;
02233         int do_soft=1, full_osa=0, i;
02234 
02235         float min[3], max[3];
02236         RayHint bb_hint;
02237 
02238         float jitco[RE_MAX_OSA][3];
02239         int totjitco;
02240 
02241         colsq[0] = colsq[1] = colsq[2] = 0.0;
02242         if(isec->mode==RE_RAY_SHADOW_TRA) {
02243                 shadfac[0]= shadfac[1]= shadfac[2]= shadfac[3]= 0.0f;
02244         } else
02245                 shadfac[3]= 1.0f;
02246         
02247         if (lar->ray_totsamp < 2) do_soft = 0;
02248         if ((R.r.mode & R_OSA) && (R.osa > 0) && (shi->vlr->flag & R_FULL_OSA)) full_osa = 1;
02249         
02250         if (full_osa) {
02251                 if (do_soft) max_samples  = max_samples/R.osa + 1;
02252                 else max_samples = 1;
02253         } else {
02254                 if (do_soft) max_samples = lar->ray_totsamp;
02255                 else if (shi->depth == 0) max_samples = (R.osa > 4)?R.osa:5;
02256                 else max_samples = 1;
02257         }
02258         
02259         ray_shadow_jittered_coords(shi, max_samples, jitco, &totjitco);
02260 
02261         /* sampling init */
02262         if (lar->ray_samp_method==LA_SAMP_HALTON)
02263                 qsa = get_thread_qmcsampler(&R, shi->thread, SAMP_TYPE_HALTON, max_samples);
02264         else if (lar->ray_samp_method==LA_SAMP_HAMMERSLEY)
02265                 qsa = get_thread_qmcsampler(&R, shi->thread, SAMP_TYPE_HAMMERSLEY, max_samples);
02266         
02267         QMC_initPixel(qsa, shi->thread);
02268 
02269         INIT_MINMAX(min, max);
02270         for(i=0; i<totjitco; i++)
02271         {
02272                 DO_MINMAX(jitco[i], min, max);
02273         }
02274         RE_rayobject_hint_bb( R.raytree, &bb_hint, min, max);
02275         
02276         isec->hint = &bb_hint;
02277         isec->check = RE_CHECK_VLR_RENDER;
02278         isec->skip = RE_SKIP_VLR_NEIGHBOUR;
02279         VECCOPY(vec, lampco);
02280         
02281         while (samples < max_samples) {
02282 
02283                 isec->orig.ob   = shi->obi;
02284                 isec->orig.face = shi->vlr;
02285 
02286                 /* manually jitter the start shading co-ord per sample
02287                  * based on the pre-generated OSA texture sampling offsets, 
02288                  * for anti-aliasing sharp shadow edges. */
02289                 co = jitco[samples % totjitco];
02290 
02291                 if (do_soft) {
02292                         /* sphere shadow source */
02293                         if (lar->type == LA_LOCAL) {
02294                                 float ru[3], rv[3], v[3], s[3];
02295                                 
02296                                 /* calc tangent plane vectors */
02297                                 v[0] = co[0] - lampco[0];
02298                                 v[1] = co[1] - lampco[1];
02299                                 v[2] = co[2] - lampco[2];
02300                                 normalize_v3(v);
02301                                 ortho_basis_v3v3_v3( ru, rv,v);
02302                                 
02303                                 /* sampling, returns quasi-random vector in area_size disc */
02304                                 QMC_sampleDisc(samp3d, qsa, shi->thread, samples,lar->area_size);
02305 
02306                                 /* distribute disc samples across the tangent plane */
02307                                 s[0] = samp3d[0]*ru[0] + samp3d[1]*rv[0];
02308                                 s[1] = samp3d[0]*ru[1] + samp3d[1]*rv[1];
02309                                 s[2] = samp3d[0]*ru[2] + samp3d[1]*rv[2];
02310                                 
02311                                 VECCOPY(samp3d, s);
02312                         }
02313                         else {
02314                                 /* sampling, returns quasi-random vector in [sizex,sizey]^2 plane */
02315                                 QMC_sampleRect(samp3d, qsa, shi->thread, samples, lar->area_size, lar->area_sizey);
02316                                                                 
02317                                 /* align samples to lamp vector */
02318                                 mul_m3_v3(lar->mat, samp3d);
02319                         }
02320                         end[0] = vec[0]+samp3d[0];
02321                         end[1] = vec[1]+samp3d[1];
02322                         end[2] = vec[2]+samp3d[2];
02323                 } else {
02324                         VECCOPY(end, vec);
02325                 }
02326 
02327                 if(shi->strand) {
02328                         /* bias away somewhat to avoid self intersection */
02329                         float jitbias= 0.5f*(len_v3(shi->dxco) + len_v3(shi->dyco));
02330                         float v[3];
02331 
02332                         VECSUB(v, co, end);
02333                         normalize_v3(v);
02334 
02335                         co[0] -= jitbias*v[0];
02336                         co[1] -= jitbias*v[1];
02337                         co[2] -= jitbias*v[2];
02338                 }
02339 
02340                 VECCOPY(isec->start, co);
02341                 isec->dir[0] = end[0]-isec->start[0];
02342                 isec->dir[1] = end[1]-isec->start[1];
02343                 isec->dir[2] = end[2]-isec->start[2];
02344                 isec->dist = normalize_v3(isec->dir);
02345                 
02346                 /* trace the ray */
02347                 if(isec->mode==RE_RAY_SHADOW_TRA) {
02348                         float col[4] = {1.0f, 1.0f, 1.0f, 1.0f};
02349                         
02350                         ray_trace_shadow_tra(isec, shi, DEPTH_SHADOW_TRA, 0, col);
02351                         shadfac[0] += col[0];
02352                         shadfac[1] += col[1];
02353                         shadfac[2] += col[2];
02354                         shadfac[3] += col[3];
02355                         
02356                         /* for variance calc */
02357                         colsq[0] += col[0]*col[0];
02358                         colsq[1] += col[1]*col[1];
02359                         colsq[2] += col[2]*col[2];
02360                 }
02361                 else {
02362                         if( RE_rayobject_raycast(R.raytree, isec) ) fac+= 1.0f;
02363                 }
02364                 
02365                 samples++;
02366                 
02367                 if (lar->ray_samp_method == LA_SAMP_HALTON) {
02368                 
02369                         /* adaptive sampling - consider samples below threshold as in shadow (or vice versa) and exit early */
02370                         if ((max_samples > min_adapt_samples) && (adapt_thresh > 0.0) && (samples > max_samples / 3)) {
02371                                 if (isec->mode==RE_RAY_SHADOW_TRA) {
02372                                         if ((shadfac[3] / samples > (1.0-adapt_thresh)) || (shadfac[3] / samples < adapt_thresh))
02373                                                 break;
02374                                         else if (adaptive_sample_variance(samples, shadfac, colsq, adapt_thresh))
02375                                                 break;
02376                                 } else {
02377                                         if ((fac / samples > (1.0-adapt_thresh)) || (fac / samples < adapt_thresh))
02378                                                 break;
02379                                 }
02380                         }
02381                 }
02382         }
02383         
02384         if(isec->mode==RE_RAY_SHADOW_TRA) {
02385                 shadfac[0] /= samples;
02386                 shadfac[1] /= samples;
02387                 shadfac[2] /= samples;
02388                 shadfac[3] /= samples;
02389         } else
02390                 shadfac[3]= 1.0f-fac/samples;
02391 
02392         if (qsa)
02393                 release_thread_qmcsampler(&R, shi->thread, qsa);
02394 }
02395 
02396 static void ray_shadow_jitter(ShadeInput *shi, LampRen *lar, float *lampco, float *shadfac, Isect *isec)
02397 {
02398         /* area soft shadow */
02399         float *jitlamp;
02400         float fac=0.0f, div=0.0f, vec[3];
02401         int a, j= -1, mask;
02402         RayHint point_hint;
02403         
02404         if(isec->mode==RE_RAY_SHADOW_TRA) {
02405                 shadfac[0]= shadfac[1]= shadfac[2]= shadfac[3]= 0.0f;
02406         }
02407         else shadfac[3]= 1.0f;
02408         
02409         fac= 0.0f;
02410         jitlamp= give_jitter_plane(lar, shi->thread, shi->xs, shi->ys);
02411 
02412         a= lar->ray_totsamp;
02413         
02414         /* this correction to make sure we always take at least 1 sample */
02415         mask= shi->mask;
02416         if(a==4) mask |= (mask>>4)|(mask>>8);
02417         else if(a==9) mask |= (mask>>9);
02418         
02419         VECCOPY(isec->start, shi->co);          
02420         isec->orig.ob   = shi->obi;
02421         isec->orig.face = shi->vlr;
02422         RE_rayobject_hint_bb( R.raytree, &point_hint, isec->start, isec->start );
02423         isec->hint = &point_hint;
02424         
02425         while(a--) {
02426                 
02427                 if(R.r.mode & R_OSA) {
02428                         j++;
02429                         if(j>=R.osa) j= 0;
02430                         if(!(mask & (1<<j))) {
02431                                 jitlamp+= 2;
02432                                 continue;
02433                         }
02434                 }
02435                 
02436                 vec[0]= jitlamp[0];
02437                 vec[1]= jitlamp[1];
02438                 vec[2]= 0.0f;
02439                 mul_m3_v3(lar->mat, vec);
02440                 
02441                 /* set start and vec */
02442                 isec->dir[0] = vec[0]+lampco[0]-isec->start[0];
02443                 isec->dir[1] = vec[1]+lampco[1]-isec->start[1];
02444                 isec->dir[2] = vec[2]+lampco[2]-isec->start[2];
02445                 isec->dist = 1.0f;
02446                 isec->check = RE_CHECK_VLR_RENDER;
02447                 isec->skip = RE_SKIP_VLR_NEIGHBOUR;
02448                 
02449                 if(isec->mode==RE_RAY_SHADOW_TRA) {
02450                         /* isec.col is like shadfac, so defines amount of light (0.0 is full shadow) */
02451                         float col[4] = {1.0f, 1.0f, 1.0f, 1.0f};
02452                         
02453                         ray_trace_shadow_tra(isec, shi, DEPTH_SHADOW_TRA, 0, col);
02454                         shadfac[0] += col[0];
02455                         shadfac[1] += col[1];
02456                         shadfac[2] += col[2];
02457                         shadfac[3] += col[3];
02458                 }
02459                 else if( RE_rayobject_raycast(R.raytree, isec) ) fac+= 1.0f;
02460                 
02461                 div+= 1.0f;
02462                 jitlamp+= 2;
02463         }
02464         
02465         if(isec->mode==RE_RAY_SHADOW_TRA) {
02466                 shadfac[0] /= div;
02467                 shadfac[1] /= div;
02468                 shadfac[2] /= div;
02469                 shadfac[3] /= div;
02470         }
02471         else {
02472                 // sqrt makes nice umbra effect
02473                 if(lar->ray_samp_type & LA_SAMP_UMBRA)
02474                         shadfac[3]= sqrt(1.0f-fac/div);
02475                 else
02476                         shadfac[3]= 1.0f-fac/div;
02477         }
02478 }
02479 /* extern call from shade_lamp_loop */
02480 void ray_shadow(ShadeInput *shi, LampRen *lar, float *shadfac)
02481 {
02482         Isect isec;
02483         float lampco[3];
02484 
02485         /* setup isec */
02486         RE_RC_INIT(isec, *shi);
02487         if(shi->mat->mode & MA_SHADOW_TRA) isec.mode= RE_RAY_SHADOW_TRA;
02488         else isec.mode= RE_RAY_SHADOW;
02489         isec.hint = 0;
02490         
02491         if(lar->mode & (LA_LAYER|LA_LAYER_SHADOW))
02492                 isec.lay= lar->lay;
02493         else
02494                 isec.lay= -1;
02495 
02496         /* only when not mir tracing, first hit optimm */
02497         if(shi->depth==0) {
02498                 isec.last_hit = lar->last_hit[shi->thread];
02499         }
02500         else {
02501                 isec.last_hit = NULL;
02502         }
02503         
02504         if(lar->type==LA_SUN || lar->type==LA_HEMI) {
02505                 /* jitter and QMC sampling add a displace vector to the lamp position
02506                  * that's incorrect because a SUN lamp does not has an exact position
02507                  * and the displace should be done at the ray vector instead of the
02508                  * lamp position.
02509                  * This is easily verified by noticing that shadows of SUN lights change
02510                  * with the scene BB.
02511                  * 
02512                  * This was detected during SoC 2009 - Raytrace Optimization, but to keep
02513                  * consistency with older render code it wasn't removed.
02514                  * 
02515                  * If the render code goes through some recode/serious bug-fix then this
02516                  * is something to consider!
02517                  */
02518                 lampco[0]= shi->co[0] - R.maxdist*lar->vec[0];
02519                 lampco[1]= shi->co[1] - R.maxdist*lar->vec[1];
02520                 lampco[2]= shi->co[2] - R.maxdist*lar->vec[2];
02521         }
02522         else {
02523                 VECCOPY(lampco, lar->co);
02524         }
02525         
02526         if (ELEM(lar->ray_samp_method, LA_SAMP_HALTON, LA_SAMP_HAMMERSLEY)) {
02527                 
02528                 ray_shadow_qmc(shi, lar, lampco, shadfac, &isec);
02529                 
02530         } else {
02531                 if(lar->ray_totsamp<2) {
02532                         
02533                         isec.orig.ob   = shi->obi;
02534                         isec.orig.face = shi->vlr;
02535                         
02536                         shadfac[3]= 1.0f; // 1.0=full light
02537                         
02538                         /* set up isec.dir */
02539                         VECCOPY(isec.start, shi->co);
02540                         VECSUB(isec.dir, lampco, isec.start);
02541                         isec.dist = normalize_v3(isec.dir);
02542 
02543                         if(isec.mode==RE_RAY_SHADOW_TRA) {
02544                                 /* isec.col is like shadfac, so defines amount of light (0.0 is full shadow) */
02545                                 float col[4] = {1.0f, 1.0f, 1.0f, 1.0f};
02546 
02547                                 ray_trace_shadow_tra(&isec, shi, DEPTH_SHADOW_TRA, 0, col);
02548                                 QUATCOPY(shadfac, col);
02549                         }
02550                         else if(RE_rayobject_raycast(R.raytree, &isec))
02551                                 shadfac[3]= 0.0f;
02552                 }
02553                 else {
02554                         ray_shadow_jitter(shi, lar, lampco, shadfac, &isec);
02555                 }
02556         }
02557                 
02558         /* for first hit optim, set last interesected shadow face */
02559         if(shi->depth==0) {
02560                 lar->last_hit[shi->thread] = isec.last_hit;
02561         }
02562 
02563 }
02564 
02565 #if 0
02566 /* only when face points away from lamp, in direction of lamp, trace ray and find first exit point */
02567 static void ray_translucent(ShadeInput *shi, LampRen *lar, float *distfac, float *co)
02568 {
02569         Isect isec;
02570         float lampco[3];
02571         
02572         assert(0);
02573         
02574         /* setup isec */
02575         RE_RC_INIT(isec, *shi);
02576         isec.mode= RE_RAY_SHADOW_TRA;
02577         isec.hint = 0;
02578         
02579         if(lar->mode & LA_LAYER) isec.lay= lar->lay; else isec.lay= -1;
02580         
02581         if(lar->type==LA_SUN || lar->type==LA_HEMI) {
02582                 lampco[0]= shi->co[0] - RE_RAYTRACE_MAXDIST*lar->vec[0];
02583                 lampco[1]= shi->co[1] - RE_RAYTRACE_MAXDIST*lar->vec[1];
02584                 lampco[2]= shi->co[2] - RE_RAYTRACE_MAXDIST*lar->vec[2];
02585         }
02586         else {
02587                 VECCOPY(lampco, lar->co);
02588         }
02589         
02590         isec.orig.ob   = shi->obi;
02591         isec.orig.face = shi->vlr;
02592         
02593         /* set up isec.dir */
02594         VECCOPY(isec.start, shi->co);
02595         VECCOPY(isec.end, lampco);
02596         
02597         if(RE_rayobject_raycast(R.raytree, &isec)) {
02598                 /* we got a face */
02599                 
02600                 /* render co */
02601                 co[0]= isec.start[0]+isec.dist*(isec.dir[0]);
02602                 co[1]= isec.start[1]+isec.dist*(isec.dir[1]);
02603                 co[2]= isec.start[2]+isec.dist*(isec.dir[2]);
02604                 
02605                 *distfac= len_v3(isec.dir);
02606         }
02607         else
02608                 *distfac= 0.0f;
02609 }
02610 
02611 #endif
02612