Blender V2.61 - r43446

rayshade.c

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