svn merge -r 22571:22800 https://svn.blender.org/svnroot/bf-blender/trunk/blender
authorAndre Susano Pinto <andresusanopinto@gmail.com>
Tue, 15 Sep 2009 15:15:43 +0000 (15:15 +0000)
committerAndre Susano Pinto <andresusanopinto@gmail.com>
Tue, 15 Sep 2009 15:15:43 +0000 (15:15 +0000)
svn merge -r 22800:23207 https://svn.blender.org/svnroot/bf-blender/trunk/blender

Merged volumetric with new raytrace code (it compiles and rendered volume-cube.blend withouth problems)

Part1:
source/blender

18 files changed:
1  2 
source/blender/blenkernel/BKE_utildefines.h
source/blender/blenlib/intern/BLI_kdopbvh.c
source/blender/editors/armature/meshlaplacian.c
source/blender/makesdna/DNA_scene_types.h
source/blender/makesrna/intern/rna_scene.c
source/blender/render/SConscript
source/blender/render/extern/include/RE_raytrace.h
source/blender/render/extern/include/RE_shader_ext.h
source/blender/render/intern/include/rayobject.h
source/blender/render/intern/include/render_types.h
source/blender/render/intern/include/rendercore.h
source/blender/render/intern/source/rayobject.c
source/blender/render/intern/source/rayshade.c
source/blender/render/intern/source/rendercore.c
source/blender/render/intern/source/shadeinput.c
source/blender/render/intern/source/volume_precache.c
source/blender/render/intern/source/volumetric.c
source/blender/windowmanager/intern/wm_init_exit.c

@@@ -1,10 -1,8 +1,10 @@@
  #!/usr/bin/python
  Import ('env')
  
- cflags = ['-O3','-msse2','-mfpmath=sse']
- cxxflags = ['-O3','-msse2','-mfpmath=sse']
 -cflags=''
++cflags = ['-msse2','-mfpmath=sse']
++cxxflags = ['-msse2','-mfpmath=sse']
  sources = env.Glob('intern/source/*.c')
 +raysources = env.Glob('intern/raytrace/*.cpp')
  
  incs = 'intern/include #/intern/guardedalloc ../blenlib ../makesdna ../makesrna'
  incs += ' extern/include ../blenkernel ../radiosity/extern/include ../imbuf'
  #ifndef RE_RAYTRACE_H
  #define RE_RAYTRACE_H
  
 -/* ray types */
 -#define RE_RAY_SHADOW 0
 -#define RE_RAY_MIRROR 1
 -#define RE_RAY_SHADOW_TRA 2
 +#ifdef __cplusplus
 +extern "C" {
 +#endif
  
 -/* spatial tree for raytracing acceleration */
 -typedef void RayTree;
 -/* abstraction of face type */
 -typedef void RayFace;
 +//#define RE_RAYCOUNTER                       /* enable counters per ray, usefull for measuring raytrace structures performance */
  
 -/* object numbers above this are transformed */
 -#define RE_RAY_TRANSFORM_OFFS 0x8000000
 +#define RE_RAY_LCTS_MAX_SIZE  256
 +#define RT_USE_LAST_HIT                       /* last shadow hit is reused before raycasting on whole tree */
 +//#define RT_USE_HINT                 /* last hit object is reused before raycasting on whole tree */
  
 -/* convert from pointer to index in array and back, with offset if the
 - * instance is transformed */
 -#define RAY_OBJECT_SET(re, obi) \
 -      ((obi == NULL)? 0: \
 -      ((obi - (re)->objectinstance) + ((obi->flag & R_TRANSFORMED)? RE_RAY_TRANSFORM_OFFS: 0)))
  
 -#define RAY_OBJECT_GET(re, i) \
 -      ((re)->objectinstance + ((i >= RE_RAY_TRANSFORM_OFFS)? i-RE_RAY_TRANSFORM_OFFS: i))
  
 +/* Internals about raycasting structures can be found on intern/raytree.h */
 +typedef struct RayObject RayObject;
 +typedef struct Isect Isect;
 +typedef struct RayHint RayHint;
 +typedef struct RayTraceHint RayTraceHint;
 +typedef struct RayCounter RayCounter;
  
 -/* struct for intersection data */
 -typedef struct Isect {
 -      float start[3];                 /* start+vec = end, in ray_tree_intersect */
 -      float vec[3];
 -      float end[3];                   
 +struct DerivedMesh;
 +struct Mesh;
 +
 +int  RE_rayobject_raycast(RayObject *r, Isect *i);
 +void RE_rayobject_add    (RayObject *r, RayObject *);
 +void RE_rayobject_done(RayObject *r);
 +void RE_rayobject_free(RayObject *r);
 +
++/* Extend min/max coords so that the rayobject is inside them */
++void RE_rayobject_merge_bb(RayObject *ob, float *min, float *max);
++
 +/* initializes an hint for optiming raycast where it is know that a ray will pass by the given BB often the origin point */
 +void RE_rayobject_hint_bb(RayObject *r, RayHint *hint, float *min, float *max);
 +
 +/* initializes an hint for optiming raycast where it is know that a ray will be contained inside the given cone*/
 +/* void RE_rayobject_hint_cone(RayObject *r, RayHint *hint, float *); */
 +
 +/* RayObject constructors */
 +RayObject* RE_rayobject_octree_create(int ocres, int size);
 +RayObject* RE_rayobject_instance_create(RayObject *target, float transform[][4], void *ob, void *target_ob);
 +
 +RayObject* RE_rayobject_blibvh_create(int size);      /* BLI_kdopbvh.c   */
 +RayObject* RE_rayobject_bvh_create(int size);         /* raytrace/rayobject_bvh.c */
 +RayObject* RE_rayobject_vbvh_create(int size);                /* raytrace/rayobject_vbvh.c */
 +RayObject* RE_rayobject_qbvh_create(int size);                /* raytrace/rayobject_qbvh.c */
 +RayObject* RE_rayobject_svbvh_create(int size);               /* raytrace/rayobject_svbvh.c */
 +RayObject* RE_rayobject_bih_create(int size);         /* rayobject_bih.c */
 +
 +
 +typedef struct LCTSHint LCTSHint;
 +struct LCTSHint
 +{
 +      int size;
 +      RayObject *stack[RE_RAY_LCTS_MAX_SIZE];
 +};
  
 -      float labda, u, v;              /* distance to hitpoint, uv weights */
 +struct RayHint
 +{
 +      union
 +      {
 +              LCTSHint lcts;
 +      } data;
 +};
  
 -      RayFace *face;                  /* face is where to intersect with */
 -      int ob;
 -      RayFace *faceorig;              /* start face */
 -      int oborig;
 -      RayFace *face_last;             /* for shadow optimize, last intersected face */
 -      int ob_last;
  
 +/* Ray Intersection */
 +struct Isect
 +{
 +      float start[3];
 +      float vec[3];
 +      float labda;
 +
 +      /* length of vec, configured by RE_rayobject_raycast */
 +      int   bv_index[6];
 +      float idot_axis[3];
 +      float dist;
 +
 +/*    float end[3];                    - not used */
 +
 +      float u, v;
 +      
 +      struct
 +      {
 +              void *ob;
 +              void *face;
 +      }
 +      hit, orig;
 +      
 +      RayObject *last_hit;    /* last hit optimization */
 +
 +#ifdef RT_USE_HINT
 +      RayTraceHint *hint, *hit_hint;
 +#endif
 +      
        short isect;                    /* which half of quad */
        short mode;                             /* RE_RAY_SHADOW, RE_RAY_MIRROR, RE_RAY_SHADOW_TRA */
        int lay;                                /* -1 default, set for layer lamps */
 +      
 +      int skip;                               /* RE_SKIP_CULLFACE */
  
 -      /* only used externally */
        float col[4];                   /* RGBA for shadow_tra */
  
 -      /* octree only */
 -      RayFace *facecontr;
 -      int obcontr;
 -      float ddalabda;
 -      short faceisect;                /* flag if facecontr was done or not */
 -
 -      /* custom pointer to be used in the RayCheckFunc */
        void *userdata;
 -} Isect;
 -
 -/* function callbacks for face type abstraction */
 -typedef void (*RayCoordsFunc)(RayFace *face,
 -      float **v1, float **v2, float **v3, float **v4);
 -typedef int (*RayCheckFunc)(Isect *is, int ob, RayFace *face);
 -typedef float *(*RayObjectTransformFunc)(void *userdata, int ob);
 -
 -/* tree building and freeing */
 -RayTree *RE_ray_tree_create(int ocres, int totface, float *min, float *max,
 -      RayCoordsFunc coordfunc, RayCheckFunc checkfunc,
 -      RayObjectTransformFunc transformfunc, void *userdata);
 -void RE_ray_tree_add_face(RayTree *tree, int ob, RayFace *face);
 -void RE_ray_tree_done(RayTree *tree);
 -void RE_ray_tree_free(RayTree *tree);
 -
 -/* intersection with full tree and single face */
 -int RE_ray_tree_intersect(RayTree *tree, Isect *is);
 -int RE_ray_tree_intersect_check(RayTree *tree, Isect *is, RayCheckFunc check);
 -int RE_ray_face_intersection(Isect *is, RayObjectTransformFunc transformfunc,
 -      RayCoordsFunc coordsfunc);
 -
 -/* retrieve the diameter of the tree structure, for setting intersection
 -   end distance */
 -float RE_ray_tree_max_size(RayTree *tree);
 +      
 +      RayHint *hint;
 +      
 +#ifdef RE_RAYCOUNTER
 +      RayCounter *raycounter;
 +#endif
 +};
  
 -#endif /*__RE_RAYTRACE_H__*/
 +/* ray types */
 +#define RE_RAY_SHADOW 0
 +#define RE_RAY_MIRROR 1
 +#define RE_RAY_SHADOW_TRA 2
  
- #define RE_SKIP_VLR_NEIGHBOUR         (1 << 1)
- #define RE_SKIP_VLR_RENDER_CHECK      (1 << 2)
 +/* skip options */
 +#define RE_SKIP_CULLFACE              (1 << 0)
 +
 +/* if using this flag then *face should be a pointer to a VlakRen */
++#define RE_SKIP_VLR_NEIGHBOUR                 (1 << 1)
++#define RE_SKIP_VLR_RENDER_CHECK              (1 << 2)
++#define RE_SKIP_VLR_NON_SOLID_MATERIAL        (1 << 3)
 +
 +/* TODO use: FLT_MAX? */
 +#define RE_RAYTRACE_MAXDIST   1e33
 +
 +#ifdef __cplusplus
 +}
 +#endif
 +
 +
 +#endif /*__RE_RAYTRACE_H__*/
index 2a9b42d,0000000..19608fb
mode 100644,000000..100644
--- /dev/null
@@@ -1,219 -1,0 +1,212 @@@
- /*
-  * Extend min/max coords so that the rayobject is inside them
-  */
- void RE_rayobject_merge_bb(RayObject *ob, float *min, float *max);
 +/**
 + * $Id$
 + *
 + * ***** BEGIN GPL LICENSE BLOCK *****
 + *
 + * This program is free software; you can redistribute it and/or
 + * modify it under the terms of the GNU General Public License
 + * as published by the Free Software Foundation; either version 2
 + * of the License, or (at your option) any later version. 
 + *
 + * This program is distributed in the hope that it will be useful,
 + * but WITHOUT ANY WARRANTY; without even the implied warranty of
 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 + * GNU General Public License for more details.
 + *
 + * You should have received a copy of the GNU General Public License
 + * along with this program; if not, write to the Free Software Foundation,
 + * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 + *
 + * The Original Code is Copyright (C) 2009 Blender Foundation.
 + * All rights reserved.
 + *
 + * The Original Code is: all of this file.
 + *
 + * Contributor(s): AndrĂ© Pinto.
 + *
 + * ***** END GPL LICENSE BLOCK *****
 + */
 +#ifndef RE_RAYOBJECT_H
 +#define RE_RAYOBJECT_H
 +
 +#ifdef __cplusplus
 +extern "C" {
 +#endif
 +
 +#include "RE_raytrace.h"
 +#include "render_types.h"
 +#include <float.h>
 +
 +
 +/* RayObject
 +      
 +      A ray object is everything where we can cast rays like:
 +              * a face/triangle
 +              * an octree
 +              * a bvh tree
 +              * an octree of bvh's
 +              * a bvh of bvh's
 +      
 +              
 +      All types of RayObjects can be created by implementing the
 +      callbacks of the RayObject.
 +
 +      Due to high computing time evolved with casting on faces
 +      there is a special type of RayObject (named RayFace)
 +      which won't use callbacks like other generic nodes.
 +      
 +      In order to allow a mixture of RayFace+RayObjects,
 +      all RayObjects must be 4byte aligned, allowing us to use the
 +      2 least significant bits (with the mask 0x02) to define the
 +      type of RayObject.
 +      
 +      This leads to 4 possible types of RayObject, but at the moment
 +      only 2 are used:
 +
 +       addr&2  - type of object
 +              0               Self (reserved for each structure)
 +              1       RayFace
 +              2               RayObject (generic with API callbacks)
 +              3               RayObject_Vlak
 +
 +      0 means it's reserved and has it own meaning inside each ray acceleration structure
 +      (this way each structure can use the allign offset to determine if a node represents a
 +       RayObject primitive, which can be used to save memory)
 +
 +      You actually don't need to care about this if you are only using the API
 +      described on RE_raytrace.h
 + */
 +
 +/* used to align a given ray object */
 +#define RE_rayobject_align(o)                         ((RayObject*)(((intptr_t)o)&(~3)))
 +
 +/* used to unalign a given ray object */
 +#define RE_rayobject_unalignRayFace(o)                ((RayObject*)(((intptr_t)o)|1))
 +#define RE_rayobject_unalignRayAPI(o)         ((RayObject*)(((intptr_t)o)|2))
 +#define RE_rayobject_unalignRayVlak(o)                ((RayObject*)(((intptr_t)o)|3))
 +
 +/* used to test the type of ray object */
 +#define RE_rayobject_isAligned(o)     ((((intptr_t)o)&3) == 0)
 +#define RE_rayobject_isRayFace(o)     ((((intptr_t)o)&3) == 1)
 +#define RE_rayobject_isRayAPI(o)      ((((intptr_t)o)&3) == 2)
 +#define RE_rayobject_isRayVlak(o)     ((((intptr_t)o)&3) == 3)
 +
 +
 +/*
 + * This ray object represents faces directly from a given VlakRen structure.
 + * Thus allowing to save memory, but making code dependant on render structures
 +typedef struct RayVlak
 +{
 +      struct ObjectInstanceRen *ob;
 +      struct VlakRen *face;
 +} RayVlak;
 + */
 +
 +/*
 + * This ray object represents a triangle or a quad face.
 + * All data needed to realize intersection is "localy" available.
 + */
 +typedef struct RayFace
 +{
 +      float v1[4], v2[4], v3[4], v4[3];
 +      int quad;
 +      void *ob;
 +      void *face;
 +      
 +} RayFace;
 +
 +#define RE_rayface_isQuad(a) ((a)->quad)
 +/* Loads a VlakRen on a RayFace */
 +void RE_rayface_from_vlak(RayFace *face, ObjectInstanceRen *obi, VlakRen *vlr);
 +
 +
 +
 +/*
 + * This rayobject represents a generic object. With it's own callbacks for raytrace operations.
 + * It's suitable to implement things like LOD.
 + */
 +struct RayObject
 +{
 +      struct RayObjectAPI *api;
 +};
 +
 +
 +typedef int  (*RE_rayobject_raycast_callback)(RayObject *, Isect *);
 +typedef void (*RE_rayobject_add_callback)(RayObject *raytree, RayObject *rayobject);
 +typedef void (*RE_rayobject_done_callback)(RayObject *);
 +typedef void (*RE_rayobject_free_callback)(RayObject *);
 +typedef void (*RE_rayobject_merge_bb_callback)(RayObject *, float *min, float *max);
 +typedef float (*RE_rayobject_cost_callback)(RayObject *);
 +typedef void (*RE_rayobject_hint_bb_callback)(RayObject *, RayHint *, float *, float *);
 +
 +typedef struct RayObjectAPI
 +{
 +      RE_rayobject_raycast_callback   raycast;
 +      RE_rayobject_add_callback               add;
 +      RE_rayobject_done_callback              done;
 +      RE_rayobject_free_callback              free;
 +      RE_rayobject_merge_bb_callback  bb;
 +      RE_rayobject_cost_callback              cost;
 +      RE_rayobject_hint_bb_callback   hint_bb;
 +      
 +} RayObjectAPI;
 +
 +
 +/*
 + * This function differs from RE_rayobject_raycast
 + * RE_rayobject_intersect does NOT perform last-hit optimization
 + * So this is probably a function to call inside raytrace structures
 + */
 +int RE_rayobject_intersect(RayObject *r, Isect *i);
 +
 +/*
 + * Returns distance ray must travel to hit the given bounding box
 + * BB should be in format [2][3]
 + */
 +/* float RE_rayobject_bb_intersect(const Isect *i, const float *bb); */
 +int RE_rayobject_bb_intersect_test(const Isect *i, const float *bb); /* same as bb_intersect but doens't calculates distance */
 +
 +/*
 + * Returns the expected cost of raycast on this node, primitives have a cost of 1
 + */
 +float RE_rayobject_cost(RayObject *r);
 +
 +
 +
 +
 +#define ISECT_EPSILON ((float)FLT_EPSILON)
 +
 +
 +
 +#if !defined(_WIN32)
 +
 +#include <sys/time.h>
 +#include <time.h>
 +#include <stdio.h>
 +
 +#define BENCH(a,name) \
 +      do {                    \
 +              double _t1, _t2;                                \
 +              struct timeval _tstart, _tend;  \
 +              clock_t _clock_init = clock();  \
 +              gettimeofday ( &_tstart, NULL); \
 +              (a);                                                    \
 +              gettimeofday ( &_tend, NULL);   \
 +              _t1 = ( double ) _tstart.tv_sec + ( double ) _tstart.tv_usec/ ( 1000*1000 );    \
 +              _t2 = ( double )   _tend.tv_sec + ( double )   _tend.tv_usec/ ( 1000*1000 );    \
 +              printf("BENCH:%s: %fs (real) %fs (cpu)\n", #name, _t2-_t1, (float)(clock()-_clock_init)/CLOCKS_PER_SEC);\
 +      } while(0)
 +#else
 +
 +#define BENCH(a)      (a)
 +
 +#endif
 +
 +
 +
 +#ifdef __cplusplus
 +}
 +#endif
 +
 +
 +#endif
@@@ -302,13 -295,11 +305,15 @@@ typedef struct ObjectInstanceRen 
  
        float dupliorco[3], dupliuv[2];
        float (*duplitexmat)[4];
+       
+       struct VolumePrecache *volume_precache;
+       
        float *vectors;
        int totvector;
 +      
 +      /* used on makeraytree */
 +      struct RayObject *raytree;
 +
  } ObjectInstanceRen;
  
  /* ------------------------------------------------------------------------- */
@@@ -416,6 -407,46 +421,46 @@@ typedef struct StrandRen 
        float orco[3];
  } StrandRen;
  
 -      struct RayTree *tree;
+ /* ------------------------------------------------------------------------- */
+ typedef struct VolumeOb
+ {
+       struct VolumeOb *next, *prev;
+       struct Material *ma;
+       struct ObjectRen *obr;
+ } VolumeOb;
+ typedef struct MatInside {
+       struct MatInside *next, *prev;
+       struct Material *ma;
+       struct ObjectInstanceRen *obi;
+ } MatInside;
+ typedef struct VolPrecachePart
+ {
+       struct VolPrecachePart *next, *prev;
++      struct RayObject *tree;
+       struct ShadeInput *shi;
+       struct ObjectInstanceRen *obi;
+       int num;
+       int minx, maxx;
+       int miny, maxy;
+       int minz, maxz;
+       int res[3];
+       float bbmin[3];
+       float voxel[3];
+       int working, done;
+ } VolPrecachePart;
+ typedef struct VolumePrecache
+ {
+       int res[3];
+       float *data_r;
+       float *data_g;
+       float *data_b;
+ } VolumePrecache;
+ /* ------------------------------------------------------------------------- */
  
  struct LampRen;
  struct MTex;
@@@ -95,6 -95,6 +95,7 @@@ int get_sample_layers(struct RenderPar
  
  extern void freeraytree(Render *re);
  extern void makeraytree(Render *re);
++RayObject* makeraytree_object(Render *re, ObjectInstanceRen *obi);
  
  extern void ray_shadow(ShadeInput *, LampRen *, float *);
  extern void ray_trace(ShadeInput *, ShadeResult *);
index 0530836,0000000..4bd8a12
mode 100644,000000..100644
--- /dev/null
@@@ -1,457 -1,0 +1,471 @@@
 +/**
 + * $Id$
 + *
 + * ***** BEGIN GPL LICENSE BLOCK *****
 + *
 + * This program is free software; you can redistribute it and/or
 + * modify it under the terms of the GNU General Public License
 + * as published by the Free Software Foundation; either version 2
 + * of the License, or (at your option) any later version. 
 + *
 + * This program is distributed in the hope that it will be useful,
 + * but WITHOUT ANY WARRANTY; without even the implied warranty of
 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 + * GNU General Public License for more details.
 + *
 + * You should have received a copy of the GNU General Public License
 + * along with this program; if not, write to the Free Software Foundation,
 + * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 + *
 + * The Original Code is Copyright (C) 2009 Blender Foundation.
 + * All rights reserved.
 + *
 + * The Original Code is: all of this file.
 + *
 + * Contributor(s): AndrĂ© Pinto.
 + *
 + * ***** END GPL LICENSE BLOCK *****
 + */
 +#include <assert.h>
 +
 +#include "BKE_utildefines.h"
 +#include "BLI_arithb.h"
 +
 +#include "RE_raytrace.h"
 +#include "render_types.h"
 +#include "rayobject.h"
 +#include "raycounter.h"
 +
 +/*
 + * Determines the distance that the ray must travel to hit the bounding volume of the given node
 + * Based on Tactical Optimization of Ray/Box Intersection, by Graham Fyffe
 + *  [http://tog.acm.org/resources/RTNews/html/rtnv21n1.html#art9]
 + */
 +int RE_rayobject_bb_intersect_test(const Isect *isec, const float *_bb)
 +{
 +      const float *bb = _bb;
 +      
 +      float t1x = (bb[isec->bv_index[0]] - isec->start[0]) * isec->idot_axis[0];
 +      float t2x = (bb[isec->bv_index[1]] - isec->start[0]) * isec->idot_axis[0];
 +      float t1y = (bb[isec->bv_index[2]] - isec->start[1]) * isec->idot_axis[1];
 +      float t2y = (bb[isec->bv_index[3]] - isec->start[1]) * isec->idot_axis[1];
 +      float t1z = (bb[isec->bv_index[4]] - isec->start[2]) * isec->idot_axis[2];
 +      float t2z = (bb[isec->bv_index[5]] - isec->start[2]) * isec->idot_axis[2];
 +
 +      RE_RC_COUNT(isec->raycounter->bb.test);
 +      
 +      if(t1x > t2y || t2x < t1y || t1x > t2z || t2x < t1z || t1y > t2z || t2y < t1z) return 0;
 +      if(t2x < 0.0 || t2y < 0.0 || t2z < 0.0) return 0;
 +      if(t1x > isec->labda || t1y > isec->labda || t1z > isec->labda) return 0;
 +      RE_RC_COUNT(isec->raycounter->bb.hit);  
 +
 +      return 1;
 +}
 +
 +
 +/* only for self-intersecting test with current render face (where ray left) */
 +static int intersection2(VlakRen *face, float r0, float r1, float r2, float rx1, float ry1, float rz1)
 +{
 +      float co1[3], co2[3], co3[3], co4[3];
 +      float x0,x1,x2,t00,t01,t02,t10,t11,t12,t20,t21,t22;
 +      float m0, m1, m2, divdet, det, det1;
 +      float u1, v, u2;
 +
 +      VECCOPY(co1, face->v1->co);
 +      VECCOPY(co2, face->v2->co);
 +      if(face->v4)
 +      {
 +              VECCOPY(co3, face->v4->co);
 +              VECCOPY(co4, face->v3->co);
 +      }
 +      else
 +      {
 +              VECCOPY(co3, face->v3->co);
 +      }
 +
 +      t00= co3[0]-co1[0];
 +      t01= co3[1]-co1[1];
 +      t02= co3[2]-co1[2];
 +      t10= co3[0]-co2[0];
 +      t11= co3[1]-co2[1];
 +      t12= co3[2]-co2[2];
 +      
 +      x0= t11*r2-t12*r1;
 +      x1= t12*r0-t10*r2;
 +      x2= t10*r1-t11*r0;
 +
 +      divdet= t00*x0+t01*x1+t02*x2;
 +
 +      m0= rx1-co3[0];
 +      m1= ry1-co3[1];
 +      m2= rz1-co3[2];
 +      det1= m0*x0+m1*x1+m2*x2;
 +      
 +      if(divdet!=0.0f) {
 +              u1= det1/divdet;
 +
 +              if(u1<ISECT_EPSILON) {
 +                      det= t00*(m1*r2-m2*r1);
 +                      det+= t01*(m2*r0-m0*r2);
 +                      det+= t02*(m0*r1-m1*r0);
 +                      v= det/divdet;
 +
 +                      if(v<ISECT_EPSILON && (u1 + v) > -(1.0f+ISECT_EPSILON)) {
 +                              return 1;
 +                      }
 +              }
 +      }
 +
 +      if(face->v4) {
 +
 +              t20= co3[0]-co4[0];
 +              t21= co3[1]-co4[1];
 +              t22= co3[2]-co4[2];
 +
 +              divdet= t20*x0+t21*x1+t22*x2;
 +              if(divdet!=0.0f) {
 +                      u2= det1/divdet;
 +              
 +                      if(u2<ISECT_EPSILON) {
 +                              det= t20*(m1*r2-m2*r1);
 +                              det+= t21*(m2*r0-m0*r2);
 +                              det+= t22*(m0*r1-m1*r0);
 +                              v= det/divdet;
 +      
 +                              if(v<ISECT_EPSILON && (u2 + v) >= -(1.0f+ISECT_EPSILON)) {
 +                                      return 2;
 +                              }
 +                      }
 +              }
 +      }
 +      return 0;
 +}
 +
 +#include "DNA_material_types.h"
 +static int vlr_check_intersect(Isect *is, ObjectInstanceRen *obi, VlakRen *vlr)
 +{
 +      /* for baking selected to active non-traceable materials might still
 +       * be in the raytree */
 +      if(!(vlr->mat->mode & MA_TRACEBLE))
 +              return 0;
 +
 +      /* I know... cpu cycle waste, might do smarter once */
 +      if(is->mode==RE_RAY_MIRROR)
 +              return !(vlr->mat->mode & MA_ONLYCAST);
 +      else
 +              return (is->lay & obi->lay);
 +}
 +
++static int vlr_check_intersect_solid(Isect *is, ObjectInstanceRen* obi, VlakRen *vlr)
++{
++      /* solid material types only */
++      if (vlr->mat->material_type == MA_TYPE_SURFACE)
++              return 1;
++      else
++              return 0;
++}
++
 +/* ray - triangle or quad intersection */
 +/* this function shall only modify Isect if it detects an hit */
 +static int intersect_rayface(RayFace *face, Isect *is)
 +{
 +      float co1[3],co2[3],co3[3],co4[3];
 +      float x0,x1,x2,t00,t01,t02,t10,t11,t12,t20,t21,t22,r0,r1,r2;
 +      float m0, m1, m2, divdet, det1;
 +      float labda, u, v;
 +      short ok=0;
 +      
 +      if(is->orig.ob == face->ob && is->orig.face == face->face)
 +              return 0;
 +              
 +      if(is->skip & RE_SKIP_VLR_RENDER_CHECK)
 +      {
 +              if(vlr_check_intersect(is, (ObjectInstanceRen*)face->ob, (VlakRen*)face->face ) == 0)
 +                      return 0;
 +      }
++      if(is->skip & RE_SKIP_VLR_NON_SOLID_MATERIAL)
++      {
++              if(vlr_check_intersect_solid(is, (ObjectInstanceRen*)face->ob, (VlakRen*)face->face) == 0)
++                      return 0;
++      }
 +
 +      RE_RC_COUNT(is->raycounter->faces.test);
 +
 +      //Load coords
 +      VECCOPY(co1, face->v1);
 +      VECCOPY(co2, face->v2);
 +      if(RE_rayface_isQuad(face))
 +      {
 +              VECCOPY(co3, face->v4);
 +              VECCOPY(co4, face->v3);
 +      }
 +      else
 +      {
 +              VECCOPY(co3, face->v3);
 +      }
 +
 +      t00= co3[0]-co1[0];
 +      t01= co3[1]-co1[1];
 +      t02= co3[2]-co1[2];
 +      t10= co3[0]-co2[0];
 +      t11= co3[1]-co2[1];
 +      t12= co3[2]-co2[2];
 +      
 +      r0= is->vec[0];
 +      r1= is->vec[1];
 +      r2= is->vec[2];
 +      
 +      x0= t12*r1-t11*r2;
 +      x1= t10*r2-t12*r0;
 +      x2= t11*r0-t10*r1;
 +
 +      divdet= t00*x0+t01*x1+t02*x2;
 +
 +      m0= is->start[0]-co3[0];
 +      m1= is->start[1]-co3[1];
 +      m2= is->start[2]-co3[2];
 +      det1= m0*x0+m1*x1+m2*x2;
 +      
 +      if(divdet!=0.0f) {
 +
 +              divdet= 1.0f/divdet;
 +              u= det1*divdet;
 +              if(u<ISECT_EPSILON && u>-(1.0f+ISECT_EPSILON)) {
 +                      float cros0, cros1, cros2;
 +                      
 +                      cros0= m1*t02-m2*t01;
 +                      cros1= m2*t00-m0*t02;
 +                      cros2= m0*t01-m1*t00;
 +                      v= divdet*(cros0*r0 + cros1*r1 + cros2*r2);
 +
 +                      if(v<ISECT_EPSILON && (u + v) > -(1.0f+ISECT_EPSILON)) {
 +                              labda= divdet*(cros0*t10 + cros1*t11 + cros2*t12);
 +
 +                              if(labda>-ISECT_EPSILON && labda<is->labda) {
 +                                      ok= 1;
 +                              }
 +                      }
 +              }
 +      }
 +
 +      if(ok==0 && RE_rayface_isQuad(face)) {
 +
 +              t20= co3[0]-co4[0];
 +              t21= co3[1]-co4[1];
 +              t22= co3[2]-co4[2];
 +
 +              divdet= t20*x0+t21*x1+t22*x2;
 +              if(divdet!=0.0f) {
 +                      divdet= 1.0f/divdet;
 +                      u = det1*divdet;
 +                      
 +                      if(u<ISECT_EPSILON && u>-(1.0f+ISECT_EPSILON)) {
 +                              float cros0, cros1, cros2;
 +                              cros0= m1*t22-m2*t21;
 +                              cros1= m2*t20-m0*t22;
 +                              cros2= m0*t21-m1*t20;
 +                              v= divdet*(cros0*r0 + cros1*r1 + cros2*r2);
 +      
 +                              if(v<ISECT_EPSILON && (u + v) >-(1.0f+ISECT_EPSILON)) {
 +                                      labda= divdet*(cros0*t10 + cros1*t11 + cros2*t12);
 +                                      
 +                                      if(labda>-ISECT_EPSILON && labda<is->labda) {
 +                                              ok= 2;
 +                                      }
 +                              }
 +                      }
 +              }
 +      }
 +
 +      if(ok) {
 +      
 +              /* when a shadow ray leaves a face, it can be little outside the edges of it, causing
 +              intersection to be detected in its neighbour face */
 +              if(is->skip & RE_SKIP_VLR_NEIGHBOUR)
 +              {
 +                      if(labda < 0.1f && is->orig.ob == face->ob)
 +                      {
 +                              VlakRen * a = is->orig.face;
 +                              VlakRen * b = face->face;
 +
 +                              /* so there's a shared edge or vertex, let's intersect ray with face
 +                              itself, if that's true we can safely return 1, otherwise we assume
 +                              the intersection is invalid, 0 */
 +                              if(a->v1==b->v1 || a->v2==b->v1 || a->v3==b->v1 || a->v4==b->v1
 +                              || a->v1==b->v2 || a->v2==b->v2 || a->v3==b->v2 || a->v4==b->v2
 +                              || a->v1==b->v3 || a->v2==b->v3 || a->v3==b->v3 || a->v4==b->v3
 +                              || (b->v4 && (a->v1==b->v4 || a->v2==b->v4 || a->v3==b->v4 || a->v4==b->v4)))
 +                              if(!intersection2((VlakRen*)a, -r0, -r1, -r2, is->start[0], is->start[1], is->start[2]))
 +                              {
 +                                      return 0;
 +                              }
 +                      }
 +              }
 +
 +              RE_RC_COUNT(is->raycounter->faces.hit);
 +
 +              is->isect= ok;  // wich half of the quad
 +              is->labda= labda;
 +              is->u= u; is->v= v;
 +
 +              is->hit.ob   = face->ob;
 +              is->hit.face = face->face;
 +#ifdef RT_USE_LAST_HIT
 +              is->last_hit = (RayObject*) RE_rayobject_unalignRayFace(face);
 +#endif
 +              return 1;
 +      }
 +
 +      return 0;
 +}
 +
 +void RE_rayface_from_vlak(RayFace *face, ObjectInstanceRen *obi, VlakRen *vlr)
 +{
 +      VECCOPY(face->v1, vlr->v1->co);
 +      VECCOPY(face->v2, vlr->v2->co);
 +      VECCOPY(face->v3, vlr->v3->co);
 +      if(vlr->v4)
 +      {
 +              VECCOPY(face->v4, vlr->v4->co);
 +              face->quad = 1;
 +      }
 +      else
 +      {
 +              face->quad = 0;
 +      }
 +
 +      face->ob   = obi;
 +      face->face = vlr;
 +}
 +
 +
 +int RE_rayobject_raycast(RayObject *r, Isect *isec)
 +{
 +      int i;
 +      RE_RC_COUNT(isec->raycounter->raycast.test);
 +
 +      /* Setup vars used on raycast */
 +      isec->labda *= Normalize(isec->vec);
 +      isec->dist = VecLength(isec->vec);
 +      
 +      for(i=0; i<3; i++)
 +      {
 +              isec->idot_axis[i]              = 1.0f / isec->vec[i];
 +              
 +              isec->bv_index[2*i]             = isec->idot_axis[i] < 0.0 ? 1 : 0;
 +              isec->bv_index[2*i+1]   = 1 - isec->bv_index[2*i];
 +              
 +              isec->bv_index[2*i]             = i+3*isec->bv_index[2*i];
 +              isec->bv_index[2*i+1]   = i+3*isec->bv_index[2*i+1];
 +      }
 +
 +#ifdef RT_USE_LAST_HIT        
 +      /* Last hit heuristic */
 +      if(isec->mode==RE_RAY_SHADOW && isec->last_hit)
 +      {
 +              RE_RC_COUNT(isec->raycounter->rayshadow_last_hit.test);
 +              
 +              if(RE_rayobject_intersect(isec->last_hit, isec))
 +              {
 +                      RE_RC_COUNT(isec->raycounter->raycast.hit);
 +                      RE_RC_COUNT(isec->raycounter->rayshadow_last_hit.hit);
 +                      return 1;
 +              }
 +      }
 +#endif
 +
 +#ifdef RT_USE_HINT
 +      isec->hit_hint = 0;
 +#endif
 +
 +      if(RE_rayobject_intersect(r, isec))
 +      {
 +              RE_RC_COUNT(isec->raycounter->raycast.hit);
 +
 +#ifdef RT_USE_HINT
 +              isec->hint = isec->hit_hint;
 +#endif
 +              return 1;
 +      }
 +      return 0;
 +}
 +
 +int RE_rayobject_intersect(RayObject *r, Isect *i)
 +{
 +      if(RE_rayobject_isRayFace(r))
 +      {
 +              return intersect_rayface( (RayFace*) RE_rayobject_align(r), i);
 +      }
 +      else if(RE_rayobject_isRayAPI(r))
 +      {
 +              r = RE_rayobject_align( r );
 +              return r->api->raycast( r, i );
 +      }
 +      else assert(0);
 +}
 +
 +void RE_rayobject_add(RayObject *r, RayObject *o)
 +{
 +      r = RE_rayobject_align( r );
 +      return r->api->add( r, o );
 +}
 +
 +void RE_rayobject_done(RayObject *r)
 +{
 +      r = RE_rayobject_align( r );
 +      r->api->done( r );
 +}
 +
 +void RE_rayobject_free(RayObject *r)
 +{
 +      r = RE_rayobject_align( r );
 +      r->api->free( r );
 +}
 +
 +void RE_rayobject_merge_bb(RayObject *r, float *min, float *max)
 +{
 +      if(RE_rayobject_isRayFace(r))
 +      {
 +              RayFace *face = (RayFace*) RE_rayobject_align(r);
 +              
 +              DO_MINMAX( face->v1, min, max );
 +              DO_MINMAX( face->v2, min, max );
 +              DO_MINMAX( face->v3, min, max );
 +              if(RE_rayface_isQuad(face)) DO_MINMAX( face->v4, min, max );
 +      }
 +      else if(RE_rayobject_isRayAPI(r))
 +      {
 +              r = RE_rayobject_align( r );
 +              r->api->bb( r, min, max );
 +      }
 +      else assert(0);
 +}
 +
 +float RE_rayobject_cost(RayObject *r)
 +{
 +      if(RE_rayobject_isRayFace(r))
 +      {
 +              return 1.0;
 +      }
 +      else if(RE_rayobject_isRayAPI(r))
 +      {
 +              r = RE_rayobject_align( r );
 +              return r->api->cost( r );
 +      }
 +      else assert(0);
 +}
 +
 +void RE_rayobject_hint_bb(RayObject *r, RayHint *hint, float *min, float *max)
 +{
 +      if(RE_rayobject_isRayFace(r))
 +      {
 +              return;
 +      }
 +      else if(RE_rayobject_isRayAPI(r))
 +      {
 +              r = RE_rayobject_align( r );
 +              return r->api->hint_bb( r, hint, min, max );
 +      }
 +      else assert(0);
 +}
 +
  #include "pixelshading.h"
  #include "shading.h"
  #include "texture.h"
+ #include "volumetric.h"
  
  #include "RE_raytrace.h"
 +#include "rayobject.h"
 +#include "raycounter.h"
  
  #define RAY_TRA               1
  #define RAY_TRAFLIP   2
  /* only to be used here in this file, it's for speed */
  extern struct Render R;
  /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
 -
 -static void vlr_face_coords(RayFace *face, float **v1, float **v2, float **v3, float **v4)
 +RayObject *  RE_rayobject_create(int type, int size)
  {
 -      VlakRen *vlr= (VlakRen*)face;
 +      if(type == R_RAYSTRUCTURE_AUTO)
 +      {
 +              //TODO
 +//            if(detect_simd())
 +//                    type = R_RAYSTRUCTURE_SIMD_SVBVH;
 +//            else
 +//                    type = R_RAYSTRUCTURE_VBVH;
  
 -      *v1 = (vlr->v1)? vlr->v1->co: NULL;
 -      *v2 = (vlr->v2)? vlr->v2->co: NULL;
 -      *v3 = (vlr->v3)? vlr->v3->co: NULL;
 -      *v4 = (vlr->v4)? vlr->v4->co: NULL;
 +                      type = R_RAYSTRUCTURE_SIMD_QBVH;
 +      }
 +      
 +      
 +      if(type == R_RAYSTRUCTURE_OCTREE)
 +      {
 +              //TODO dynamic ocres
 +              return RE_rayobject_octree_create(R.r.ocres, size);
 +      }
 +      if(type == R_RAYSTRUCTURE_BLIBVH)
 +      {
 +              return RE_rayobject_blibvh_create(size);
 +      }
 +      if(type == R_RAYSTRUCTURE_VBVH)
 +      {
 +              return RE_rayobject_vbvh_create(size);
 +      }
 +      if(type == R_RAYSTRUCTURE_SIMD_SVBVH)
 +      {
 +              return RE_rayobject_svbvh_create(size);
 +      }
 +      if(type == R_RAYSTRUCTURE_SIMD_QBVH)
 +      {
 +              return RE_rayobject_qbvh_create(size);
 +      }
 +      if(type == R_RAYSTRUCTURE_BIH)
 +      {
 +//            return RE_rayobject_bih_create(size);
 +      }
 +      
 +      return NULL;
  }
  
 -static int vlr_check_intersect(Isect *is, int ob, RayFace *face)
 -{
 -      ObjectInstanceRen *obi= RAY_OBJECT_GET((Render*)is->userdata, ob);
 -      VlakRen *vlr = (VlakRen*)face;
 +#ifdef RE_RAYCOUNTER
 +RayCounter re_rc_counter[BLENDER_MAX_THREADS] = {};
 +#endif
  
 -      /* for baking selected to active non-traceable materials might still
 -       * be in the raytree */
 -      if(!(vlr->mat->mode & MA_TRACEBLE))
 -              return 0;
 -      /* I know... cpu cycle waste, might do smarter once */
 -      if(is->mode==RE_RAY_MIRROR)
 -              return !(vlr->mat->mode & MA_ONLYCAST);
 -      else
 -              return (is->lay & obi->lay);
 +void freeraytree(Render *re)
 +{
 +      ObjectInstanceRen *obi;
 +      
 +      if(re->raytree)
 +      {
 +              RE_rayobject_free(re->raytree);
 +              re->raytree = NULL;
 +      }
 +      if(re->rayfaces)
 +      {
 +              MEM_freeN(re->rayfaces);
 +              re->rayfaces = NULL;
 +      }
 +
 +      for(obi=re->instancetable.first; obi; obi=obi->next)
 +      {
 +              ObjectRen *obr = obi->obr;
 +              if(obr->raytree)
 +              {
 +                      RE_rayobject_free(obr->raytree);
 +                      obr->raytree = NULL;
 +              }
 +              if(obr->rayfaces)
 +              {
 +                      MEM_freeN(obr->rayfaces);
 +                      obr->rayfaces = NULL;
 +              }
 +              if(obi->raytree)
 +              {
 +                      RE_rayobject_free(obi->raytree);
 +                      obi->raytree = NULL;
 +              }
 +      }
 +      
 +#ifdef RE_RAYCOUNTER
 +      {
 +              RayCounter sum = {};
 +              int i;
 +              for(i=0; i<BLENDER_MAX_THREADS; i++)
 +                      RE_RC_MERGE(&sum, re_rc_counter+i);
 +              RE_RC_INFO(&sum);
 +      }
 +#endif
  }
  
 -static int vlr_check_intersect_solid(Isect *is, int ob, RayFace *face)
 +static int is_raytraceable_vlr(Render *re, VlakRen *vlr)
  {
 -      VlakRen *vlr = (VlakRen*)face;
 -
 -      /* solid material types only */
 -      if (vlr->mat->material_type == MA_TYPE_SURFACE)
 +      if((re->flag & R_BAKE_TRACE) || (vlr->mat->mode & MA_TRACEBLE))
 +      if(vlr->mat->material_type != MA_TYPE_WIRE)
                return 1;
 -      else
 -              return 0;
 +      return 0;
  }
  
 -static float *vlr_get_transform(void *userdata, int i)
 +static int is_raytraceable(Render *re, ObjectInstanceRen *obi)
  {
 -      ObjectInstanceRen *obi= RAY_OBJECT_GET((Render*)userdata, i);
 +      int v;
 +      ObjectRen *obr = obi->obr;
  
 -      return (obi->flag & R_TRANSFORMED)? (float*)obi->mat: NULL;
 -}
 +      if(re->excludeob && obr->ob == re->excludeob)
 +              return 0;
  
 -void freeraytree(Render *re)
 -{
 -      if(re->raytree) {
 -              RE_ray_tree_free(re->raytree);
 -              re->raytree= NULL;
 +      for(v=0;v<obr->totvlak;v++)
 +      {
 +              VlakRen *vlr = obr->vlaknodes[v>>8].vlak + (v&255);
 +              if(is_raytraceable_vlr(re, vlr))
 +                      return 1;
        }
 +      return 0;
  }
  
 -void makeraytree(Render *re)
 -{
 -      ObjectInstanceRen *obi;
 -      ObjectRen *obr;
 -      VlakRen *vlr= NULL;
 -      float min[3], max[3], co1[3], co2[3], co3[3], co4[3];
 -      double lasttime= PIL_check_seconds_timer();
 -      int v, totv = 0, totface = 0;
 -
 -      INIT_MINMAX(min, max);
 -
 -      /* first min max raytree space */
 -      for(obi=re->instancetable.first; obi; obi=obi->next) {
 -              obr= obi->obr;
 -
 -              if(re->excludeob && obr->ob == re->excludeob)
 -                      continue;
 -
 -              for(v=0;v<obr->totvlak;v++) {
 -                      if((v & 255)==0) vlr= obr->vlaknodes[v>>8].vlak;
 -                      else vlr++;
 -                      /* baking selected to active needs non-traceable too */
 -                      if((re->flag & R_BAKE_TRACE) || (vlr->mat->mode & MA_TRACEBLE)) {       
 -                              if(vlr->mat->material_type != MA_TYPE_WIRE) {
 -                                      VECCOPY(co1, vlr->v1->co);
 -                                      VECCOPY(co2, vlr->v2->co);
 -                                      VECCOPY(co3, vlr->v3->co);
 -
 -                                      if(obi->flag & R_TRANSFORMED) {
 -                                              Mat4MulVecfl(obi->mat, co1);
 -                                              Mat4MulVecfl(obi->mat, co2);
 -                                              Mat4MulVecfl(obi->mat, co3);
 -                                      }
 -
 -                                      DO_MINMAX(co1, min, max);
 -                                      DO_MINMAX(co2, min, max);
 -                                      DO_MINMAX(co3, min, max);
  
 -                                      if(vlr->v4) {
 -                                              VECCOPY(co4, vlr->v4->co);
 -                                              if(obi->flag & R_TRANSFORMED)
 -                                                      Mat4MulVecfl(obi->mat, co4);
 -                                              DO_MINMAX(co4, min, max);
 -                                      }
 +RayObject* makeraytree_object(Render *re, ObjectInstanceRen *obi)
 +{
 +      //TODO
 +      // out-of-memory safeproof
 +      // break render
 +      // update render stats
 +      ObjectRen *obr = obi->obr;
 +      
 +      if(obr->raytree == NULL)
 +      {
 +              RayObject *raytree;
 +              RayFace *face;
 +              int v;
 +              
 +              //Count faces
 +              int faces = 0;
 +              for(v=0;v<obr->totvlak;v++)
 +              {
 +                      VlakRen *vlr = obr->vlaknodes[v>>8].vlak + (v&255);
 +                      if(is_raytraceable_vlr(re, vlr))
 +                              faces++;
 +              }
 +              assert( faces > 0 );
  
 -                                      totface++;
 -                              }
 +              //Create Ray cast accelaration structure                
 +              raytree = obr->raytree = RE_rayobject_create( re->r.raytrace_tree_type, faces );
 +              face = obr->rayfaces = (RayFace*)MEM_callocN(faces*sizeof(RayFace), "ObjectRen faces");
 +              obr->rayobi = obi;
 +              
 +              for(v=0;v<obr->totvlak;v++)
 +              {
 +                      VlakRen *vlr = obr->vlaknodes[v>>8].vlak + (v&255);
 +                      if(is_raytraceable_vlr(re, vlr))
 +                      {
 +                              RE_rayface_from_vlak( face, obi, vlr );                         
 +                              RE_rayobject_add( raytree, RE_rayobject_unalignRayFace(face) );
 +                              face++;
                        }
                }
 +              RE_rayobject_done( raytree );
        }
  
 -      re->raytree= RE_ray_tree_create(re->r.ocres, totface, min, max,
 -              vlr_face_coords, vlr_check_intersect, vlr_get_transform, re);
  
 -      if(min[0] > max[0]) { /* empty raytree */
 -              RE_ray_tree_done(re->raytree);
 -              return; 
 +      if(obi->flag & R_TRANSFORMED)
 +      {
 +              obi->raytree = RE_rayobject_instance_create( obr->raytree, obi->mat, obi, obi->obr->rayobi );
        }
 +      
 +      if(obi->raytree) return obi->raytree;
 +      return obi->obr->raytree;
 +}
  
 -      for(obi=re->instancetable.first; obi; obi=obi->next) {
 -              obr= obi->obr;
 +static int has_special_rayobject(Render *re, ObjectInstanceRen *obi)
 +{
 +      if( (obi->flag & R_TRANSFORMED) )
 +      {
 +              ObjectRen *obr = obi->obr;
 +              int v, faces = 0;
 +              
 +              for(v=0;v<obr->totvlak;v++)
 +              {
 +                      VlakRen *vlr = obr->vlaknodes[v>>8].vlak + (v&255);
 +                      if(is_raytraceable_vlr(re, vlr))
 +                      {
 +                              faces++;
 +                              if(faces > 4)
 +                                      return 1;
 +                      }
 +              }
 +      }
 +      return 0;
 +}
 +/*
 + * create a single raytrace structure with all faces
 + */
 +static void makeraytree_single(Render *re)
 +{
 +      ObjectInstanceRen *obi;
 +      RayObject *raytree;
 +      RayFace *face;
 +      int faces = 0, obs = 0;
  
 -              if(re->excludeob && obr->ob == re->excludeob)
 -                      continue;
 +      for(obi=re->instancetable.first; obi; obi=obi->next)
 +      if(is_raytraceable(re, obi))
 +      {
 +              int v;
 +              ObjectRen *obr = obi->obr;
 +              obs++;
 +              
 +              if(has_special_rayobject(re, obi))
 +              {
 +                      faces++;
 +              }
 +              else
 +              {
 +                      for(v=0;v<obr->totvlak;v++)
 +                      {
 +                              VlakRen *vlr = obr->vlaknodes[v>>8].vlak + (v&255);
 +                              if(is_raytraceable_vlr(re, vlr))
 +                                      faces++;
 +                      }
 +              }
 +      }
 +      
 +      //Create raytree
 +      raytree = re->raytree = RE_rayobject_create( re->r.raytrace_tree_type, faces );
  
 -              for(v=0; v<obr->totvlak; v++, totv++) {
 -                      if((v & 255)==0) {
 -                              double time= PIL_check_seconds_timer();
 +      face    = re->rayfaces  = (RayFace*)MEM_callocN(faces*sizeof(RayFace), "Render ray faces");
 +      
 +      for(obi=re->instancetable.first; obi; obi=obi->next)
 +      if(is_raytraceable(re, obi))
 +      {
 +              if(has_special_rayobject(re, obi))
 +              {
 +                      RayObject *obj = makeraytree_object(re, obi);
 +                      RE_rayobject_add( re->raytree, obj );
 +              }
 +              else
 +              {
 +                      int v;
 +                      ObjectRen *obr = obi->obr;
 +
 +                      for(v=0;v<obr->totvlak;v++)
 +                      {
 +                              VlakRen *vlr = obr->vlaknodes[v>>8].vlak + (v&255);
 +                              if(is_raytraceable_vlr(re, vlr))
 +                              {
 +                                      RE_rayface_from_vlak(face, obi, vlr);
 +                                      if((obi->flag & R_TRANSFORMED))
 +                                      {
 +                                              Mat4MulVecfl(obi->mat, face->v1);
 +                                              Mat4MulVecfl(obi->mat, face->v2);
 +                                              Mat4MulVecfl(obi->mat, face->v3);
 +                                              if(RE_rayface_isQuad(face))
 +                                                      Mat4MulVecfl(obi->mat, face->v4);
 +                                      }
  
 -                              vlr= obr->vlaknodes[v>>8].vlak;
 -                              if(re->test_break(re->tbh))
 -                                      break;
 -                              if(time-lasttime>1.0f) {
 -                                      char str[32];
 -                                      sprintf(str, "Filling Octree: %d", totv);
 -                                      re->i.infostr= str;
 -                                      re->stats_draw(re->sdh, &re->i);
 -                                      re->i.infostr= NULL;
 -                                      lasttime= time;
 +                                      RE_rayobject_add( raytree, RE_rayobject_unalignRayFace(face) );
 +                                      face++;
                                }
                        }
 -                      else vlr++;
 -                      
 -                      if((re->flag & R_BAKE_TRACE) || (vlr->mat->mode & MA_TRACEBLE))
 -                              if(vlr->mat->material_type != MA_TYPE_WIRE)
 -                                      RE_ray_tree_add_face(re->raytree, RAY_OBJECT_SET(re, obi), vlr);
                }
        }
 +      RE_rayobject_done( raytree );   
 +}
  
 -      RE_ray_tree_done(re->raytree);
 +void makeraytree(Render *re)
 +{
 +      float min[3], max[3], sub[3];
 +      int i;
        
 -      re->i.infostr= NULL;
 +      re->i.infostr= "Make raytree";
        re->stats_draw(re->sdh, &re->i);
 +
 +      BENCH(makeraytree_single(re), tree_build);
 +              
 +      //Calculate raytree max_size
 +      //This is ONLY needed to kept a bogus behaviour of SUN and HEMI lights
 +      RE_rayobject_merge_bb( re->raytree, min, max );
 +      for(i=0; i<3; i++)
 +      {
 +              min[i] += 0.01f;
 +              max[i] += 0.01f;
 +              sub[i] = max[i]-min[i];
 +      }
 +      re->maxdist = sqrt( sub[0]*sub[0] + sub[1]*sub[1] + sub[2]*sub[2] );
 +
 +      re->i.infostr= "Raytree finished";
  }
  
- static void shade_ray(Isect *is, ShadeInput *shi, ShadeResult *shr)
+ void shade_ray(Isect *is, ShadeInput *shi, ShadeResult *shr)
  {
 -      VlakRen *vlr= (VlakRen*)is->face;
 -      ObjectInstanceRen *obi= RAY_OBJECT_GET(&R, is->ob);
 +      ObjectInstanceRen *obi= (ObjectInstanceRen*)is->hit.ob;
 +      VlakRen *vlr= (VlakRen*)is->hit.face;
        int osatex= 0;
        
        /* set up view vector */
@@@ -1709,17 -1583,10 +1735,17 @@@ static void ray_ao_qmc(ShadeInput *shi
        float dxyview[3], skyadded=0, div;
        int aocolor;
        
 -      isec.faceorig= (RayFace*)shi->vlr;
 -      isec.oborig= RAY_OBJECT_SET(&R, shi->obi);
 -      isec.face_last= NULL;
 -      isec.ob_last= 0;
 +      RE_RC_INIT(isec, *shi);
 +      isec.orig.ob   = shi->obi;
 +      isec.orig.face = shi->vlr;
-       isec.skip = RE_SKIP_VLR_NEIGHBOUR | RE_SKIP_VLR_RENDER_CHECK;
++      isec.skip = RE_SKIP_VLR_NEIGHBOUR | RE_SKIP_VLR_RENDER_CHECK | RE_SKIP_VLR_NON_SOLID_MATERIAL;
 +      isec.hint = 0;
 +
 +      isec.hit.ob   = 0;
 +      isec.hit.face = 0;
 +
 +      isec.last_hit = NULL;
 +      
        isec.mode= (R.wrld.aomode & WO_AODIST)?RE_RAY_SHADOW_TRA:RE_RAY_SHADOW;
        isec.lay= -1;
        
@@@ -2246,15 -2234,11 +2246,15 @@@ static int bake_check_intersect(Isect *
        /* no direction checking for now, doesn't always improve the result
         * (INPR(shi->facenor, bs->dir) > 0.0f); */
  
-       return (R.objectinstance[ob].obr->ob != bs->actob);
+       return (R.objectinstance[ob & ~RE_RAY_TRANSFORM_OFFS].obr->ob != bs->actob);
  }
 +#endif
  
 -static int bake_intersect_tree(RayTree* raytree, Isect* isect, float *start, float *dir, float sign, float *hitco, float *dist)
 +static int bake_intersect_tree(RayObject* raytree, Isect* isect, float *start, float *dir, float sign, float *hitco, float *dist)
  {
 +      //TODO
 +      assert( 0 );
 +#if 0
        float maxdist;
        int hit;
  
index 0000000,15d8643..ca19f25
mode 000000,100644..100644
--- /dev/null
@@@ -1,0 -1,746 +1,759 @@@
 -int intersect_outside_volume(RayTree *tree, Isect *isect, float *offset, int limit, int depth)
+ /**
+  *
+  * ***** BEGIN GPL LICENSE BLOCK *****
+  *
+  * This program is free software; you can redistribute it and/or
+  * modify it under the terms of the GNU General Public License
+  * as published by the Free Software Foundation; either version 2
+  * of the License, or (at your option) any later version.
+  *
+  * This program is distributed in the hope that it will be useful,
+  * but WITHOUT ANY WARRANTY; without even the implied warranty of
+  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  * GNU General Public License for more details.
+  *
+  * You should have received a copy of the GNU General Public License
+  * along with this program; if not, write to the Free Software Foundation,
+  * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+  *
+  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
+  * All rights reserved.
+  *
+  * The Original Code is: all of this file.
+  *
+  * Contributor(s): Matt Ebb.
+  *
+  * ***** END GPL LICENSE BLOCK *****
+  */
+ #include <math.h>
+ #include <stdlib.h>
+ #include <string.h>
+ #include <float.h>
+ #include "MEM_guardedalloc.h"
+ #include "BLI_blenlib.h"
+ #include "BLI_arithb.h"
+ #include "BLI_threads.h"
+ #include "BLI_voxel.h"
+ #include "PIL_time.h"
+ #include "RE_shader_ext.h"
+ #include "RE_raytrace.h"
+ #include "DNA_material_types.h"
+ #include "render_types.h"
++#include "rendercore.h"
+ #include "renderdatabase.h"
+ #include "volumetric.h"
+ #include "volume_precache.h"
+ #if defined( _MSC_VER ) && !defined( __cplusplus )
+ # define inline __inline
+ #endif // defined( _MSC_VER ) && !defined( __cplusplus )
+ #include "BKE_global.h"
+ /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+ /* defined in pipeline.c, is hardcopy of active dynamic allocated Render */
+ /* only to be used here in this file, it's for speed */
+ extern struct Render R;
+ /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+ /* *** utility code to set up an individual raytree for objectinstance, for checking inside/outside *** */
+ /* Recursive test for intersections, from a point inside the mesh, to outside
+  * Number of intersections (depth) determine if a point is inside or outside the mesh */
 -      if (RE_ray_tree_intersect(tree, isect)) {
++int intersect_outside_volume(RayObject *tree, Isect *isect, float *offset, int limit, int depth)
+ {
+       if (limit == 0) return depth;
+       
 -int point_inside_obi(RayTree *tree, ObjectInstanceRen *obi, float *co)
++      if (RE_rayobject_raycast(tree, isect)) {
+               float hitco[3];
+               
+               hitco[0] = isect->start[0] + isect->labda*isect->vec[0];
+               hitco[1] = isect->start[1] + isect->labda*isect->vec[1];
+               hitco[2] = isect->start[2] + isect->labda*isect->vec[2];
+               VecAddf(isect->start, hitco, offset);
+               return intersect_outside_volume(tree, isect, offset, limit-1, depth+1);
+       } else {
+               return depth;
+       }
+ }
+ /* Uses ray tracing to check if a point is inside or outside an ObjectInstanceRen */
 -      float maxsize = RE_ray_tree_max_size(tree);
++int point_inside_obi(RayObject *tree, ObjectInstanceRen *obi, float *co)
+ {
 -      isect.face_last= NULL;
+       Isect isect;
+       float vec[3] = {0.0f,0.0f,1.0f};
+       int final_depth=0, depth=0, limit=20;
+       
+       /* set up the isect */
+       memset(&isect, 0, sizeof(isect));
+       VECCOPY(isect.start, co);
++      VECCOPY(isect.vec, vec);
++      isect.labda = FLT_MAX;
++      
++      /*
+       isect.end[0] = co[0] + vec[0] * maxsize;
+       isect.end[1] = co[1] + vec[1] * maxsize;
+       isect.end[2] = co[2] + vec[2] * maxsize;
++      */
+       
+       /* and give it a little offset to prevent self-intersections */
+       VecMulf(vec, 1e-5);
+       VecAddf(isect.start, isect.start, vec);
+       
+       isect.mode= RE_RAY_MIRROR;
 -static int inside_check_func(Isect *is, int ob, RayFace *face)
++      isect.last_hit= NULL;
+       isect.lay= -1;
+       
+       final_depth = intersect_outside_volume(tree, &isect, vec, limit, depth);
+       
+       /* even number of intersections: point is outside
+        * odd number: point is inside */
+       if (final_depth % 2 == 0) return 0;
+       else return 1;
+ }
 -RayTree *create_raytree_obi(ObjectInstanceRen *obi, float *bbmin, float *bbmax)
++/*
++static int inside_check_func(Isect *is, int ob, RayObject *face)
+ {
+       return 1;
+ }
++
+ static void vlr_face_coords(RayFace *face, float **v1, float **v2, float **v3, float **v4)
+ {
+       VlakRen *vlr= (VlakRen*)face;
+       *v1 = (vlr->v1)? vlr->v1->co: NULL;
+       *v2 = (vlr->v2)? vlr->v2->co: NULL;
+       *v3 = (vlr->v3)? vlr->v3->co: NULL;
+       *v4 = (vlr->v4)? vlr->v4->co: NULL;
+ }
 -      /* create empty raytree */
++RayObject *create_raytree_obi(ObjectInstanceRen *obi, float *bbmin, float *bbmax)
+ {
+       int v;
+       VlakRen *vlr= NULL;
+       
 -      /* fill it with faces */
++      / * create empty raytree * /
+       RayTree *tree = RE_ray_tree_create(64, obi->obr->totvlak, bbmin, bbmax,
+               vlr_face_coords, inside_check_func, NULL, NULL);
+       
 -      RayTree *tree = pa->tree;
++      / * fill it with faces * /
+       for(v=0; v<obi->obr->totvlak; v++) {
+               if((v & 255)==0)
+                       vlr= obi->obr->vlaknodes[v>>8].vlak;
+               else
+                       vlr++;
+       
+               RE_ray_tree_add_face(tree, 0, vlr);
+       }
+       
+       RE_ray_tree_done(tree);
+       
+       return tree;
+ }
++*/
+ /* *** light cache filtering *** */
+ static float get_avg_surrounds(float *cache, int *res, int xx, int yy, int zz)
+ {
+       int x, y, z, x_, y_, z_;
+       int added=0;
+       float tot=0.0f;
+       
+       for (z=-1; z <= 1; z++) {
+               z_ = zz+z;
+               if (z_ >= 0 && z_ <= res[2]-1) {
+               
+                       for (y=-1; y <= 1; y++) {
+                               y_ = yy+y;
+                               if (y_ >= 0 && y_ <= res[1]-1) {
+                               
+                                       for (x=-1; x <= 1; x++) {
+                                               x_ = xx+x;
+                                               if (x_ >= 0 && x_ <= res[0]-1) {
+                                               
+                                                       if (cache[ V_I(x_, y_, z_, res) ] > 0.0f) {
+                                                               tot += cache[ V_I(x_, y_, z_, res) ];
+                                                               added++;
+                                                       }
+                                                       
+                                               }
+                                       }
+                               }
+                       }
+               }
+       }
+       
+       tot /= added;
+       
+       return ((added>0)?tot:0.0f);
+ }
+ /* function to filter the edges of the light cache, where there was no volume originally.
+  * For each voxel which was originally external to the mesh, it finds the average values of
+  * the surrounding internal voxels and sets the original external voxel to that average amount.
+  * Works almost a bit like a 'dilate' filter */
+ static void lightcache_filter(VolumePrecache *vp)
+ {
+       int x, y, z;
+       for (z=0; z < vp->res[2]; z++) {
+               for (y=0; y < vp->res[1]; y++) {
+                       for (x=0; x < vp->res[0]; x++) {
+                               /* trigger for outside mesh */
+                               if (vp->data_r[ V_I(x, y, z, vp->res) ] < -0.5f)
+                                       vp->data_r[ V_I(x, y, z, vp->res) ] = get_avg_surrounds(vp->data_r, vp->res, x, y, z);
+                               if (vp->data_g[ V_I(x, y, z, vp->res) ] < -0.5f)
+                                       vp->data_g[ V_I(x, y, z, vp->res) ] = get_avg_surrounds(vp->data_g, vp->res, x, y, z);
+                               if (vp->data_b[ V_I(x, y, z, vp->res) ] < -0.5f)
+                                       vp->data_b[ V_I(x, y, z, vp->res) ] = get_avg_surrounds(vp->data_b, vp->res, x, y, z);
+                       }
+               }
+       }
+ }
+ static inline int ms_I(int x, int y, int z, int *n) //has a pad of 1 voxel surrounding the core for boundary simulation
+ { 
+       return z*(n[1]+2)*(n[0]+2) + y*(n[0]+2) + x;
+ }
+ /* *** multiple scattering approximation *** */
+ /* get the total amount of light energy in the light cache. used to normalise after multiple scattering */
+ static float total_ss_energy(VolumePrecache *vp)
+ {
+       int x, y, z;
+       int *res = vp->res;
+       float energy=0.f;
+       
+       for (z=0; z < res[2]; z++) {
+               for (y=0; y < res[1]; y++) {
+                       for (x=0; x < res[0]; x++) {
+                               if (vp->data_r[ V_I(x, y, z, res) ] > 0.f) energy += vp->data_r[ V_I(x, y, z, res) ];
+                               if (vp->data_g[ V_I(x, y, z, res) ] > 0.f) energy += vp->data_g[ V_I(x, y, z, res) ];
+                               if (vp->data_b[ V_I(x, y, z, res) ] > 0.f) energy += vp->data_b[ V_I(x, y, z, res) ];
+                       }
+               }
+       }
+       
+       return energy;
+ }
+ static float total_ms_energy(float *sr, float *sg, float *sb, int *res)
+ {
+       int x, y, z, i;
+       float energy=0.f;
+       
+       for (z=1;z<=res[2];z++) {
+               for (y=1;y<=res[1];y++) {
+                       for (x=1;x<=res[0];x++) {
+                       
+                               i = ms_I(x,y,z,res);
+                               if (sr[i] > 0.f) energy += sr[i];
+                               if (sg[i] > 0.f) energy += sg[i];
+                               if (sb[i] > 0.f) energy += sb[i];
+                       }
+               }
+       }
+       
+       return energy;
+ }
+ static void ms_diffuse(int b, float* x0, float* x, float diff, int *n)
+ {
+       int i, j, k, l;
+       const float dt = VOL_MS_TIMESTEP;
+       const float a = dt*diff*n[0]*n[1]*n[2];
+       
+       for (l=0; l<20; l++)
+       {
+               for (k=1; k<=n[2]; k++)
+               {
+                       for (j=1; j<=n[1]; j++)
+                       {
+                               for (i=1; i<=n[0]; i++)
+                               {
+                                       x[ms_I(i,j,k,n)] = (x0[ms_I(i,j,k,n)] + a*(
+                                                x[ms_I(i-1,j,k,n)]+x[ms_I(i+1,j,k,n)]+
+                                                x[ms_I(i,j-1,k,n)]+x[ms_I(i,j+1,k,n)]+
+                                                x[ms_I(i,j,k-1,n)]+x[ms_I(i,j,k+1,n)]))/(1+6*a);
+                               }
+                       }
+               }
+       }
+ }
+ void multiple_scattering_diffusion(Render *re, VolumePrecache *vp, Material *ma)
+ {
+       const float diff = ma->vol.ms_diff * 0.001f;    /* compensate for scaling for a nicer UI range */
+       const float simframes = ma->vol.ms_steps;
+       const int shade_type = ma->vol.shade_type;
+       float fac = ma->vol.ms_intensity;
+       
+       int x, y, z, m;
+       int *n = vp->res;
+       const int size = (n[0]+2)*(n[1]+2)*(n[2]+2);
+       double time, lasttime= PIL_check_seconds_timer();
+       float total;
+       float c=1.0f;
+       int i;
+       float origf;    /* factor for blending in original light cache */
+       float energy_ss, energy_ms;
+       float *sr0=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
+       float *sr=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
+       float *sg0=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
+       float *sg=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
+       float *sb0=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
+       float *sb=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
+       total = (float)(n[0]*n[1]*n[2]*simframes);
+       
+       energy_ss = total_ss_energy(vp);
+       
+       /* Scattering as diffusion pass */
+       for (m=0; m<simframes; m++)
+       {
+               /* add sources */
+               for (z=1; z<=n[2]; z++)
+               {
+                       for (y=1; y<=n[1]; y++)
+                       {
+                               for (x=1; x<=n[0]; x++)
+                               {
+                                       i = V_I((x-1), (y-1), (z-1), n);
+                                       time= PIL_check_seconds_timer();
+                                       c++;
+                                                                               
+                                       if (vp->data_r[i] > 0.f)
+                                               sr[ms_I(x,y,z,n)] += vp->data_r[i];
+                                       if (vp->data_g[i] > 0.f)
+                                               sg[ms_I(x,y,z,n)] += vp->data_g[i];
+                                       if (vp->data_b[i] > 0.f)
+                                               sb[ms_I(x,y,z,n)] += vp->data_b[i];
+                                       
+                                       /* Displays progress every second */
+                                       if(time-lasttime>1.0f) {
+                                               char str[64];
+                                               sprintf(str, "Simulating multiple scattering: %d%%", (int)
+                                                               (100.0f * (c / total)));
+                                               re->i.infostr= str;
+                                               re->stats_draw(re->sdh, &re->i);
+                                               re->i.infostr= NULL;
+                                               lasttime= time;
+                                       }
+                               }
+                       }
+               }
+               SWAP(float *, sr, sr0);
+               SWAP(float *, sg, sg0);
+               SWAP(float *, sb, sb0);
+               /* main diffusion simulation */
+               ms_diffuse(0, sr0, sr, diff, n);
+               ms_diffuse(0, sg0, sg, diff, n);
+               ms_diffuse(0, sb0, sb, diff, n);
+               
+               if (re->test_break(re->tbh)) break;
+       }
+       
+       /* normalisation factor to conserve energy */
+       energy_ms = total_ms_energy(sr, sg, sb, n);
+       fac *= (energy_ss / energy_ms);
+       
+       /* blend multiple scattering back in the light cache */
+       if (shade_type == MA_VOL_SHADE_SINGLEPLUSMULTIPLE) {
+               /* conserve energy - half single, half multiple */
+               origf = 0.5f;
+               fac *= 0.5f;
+       } else {
+               origf = 0.0f;
+       }
+       for (z=1;z<=n[2];z++)
+       {
+               for (y=1;y<=n[1];y++)
+               {
+                       for (x=1;x<=n[0];x++)
+                       {
+                               int index=(x-1)*n[1]*n[2] + (y-1)*n[2] + z-1;
+                               vp->data_r[index] = origf * vp->data_r[index] + fac * sr[ms_I(x,y,z,n)];
+                               vp->data_g[index] = origf * vp->data_g[index] + fac * sg[ms_I(x,y,z,n)];
+                               vp->data_b[index] = origf * vp->data_b[index] + fac * sb[ms_I(x,y,z,n)];
+                       }
+               }
+       }
+       MEM_freeN(sr0);
+       MEM_freeN(sr);
+       MEM_freeN(sg0);
+       MEM_freeN(sg);
+       MEM_freeN(sb0);
+       MEM_freeN(sb);
+ }
+ #if 0 // debug stuff
+ static void *vol_precache_part_test(void *data)
+ {
+       VolPrecachePart *pa = data;
+       printf("part number: %d \n", pa->num);
+       printf("done: %d \n", pa->done);
+       printf("x min: %d   x max: %d \n", pa->minx, pa->maxx);
+       printf("y min: %d   y max: %d \n", pa->miny, pa->maxy);
+       printf("z min: %d   z max: %d \n", pa->minz, pa->maxz);
+       return NULL;
+ }
+ #endif
+ /* Iterate over the 3d voxel grid, and fill the voxels with scattering information
+  *
+  * It's stored in memory as 3 big float grids next to each other, one for each RGB channel.
+  * I'm guessing the memory alignment may work out better this way for the purposes
+  * of doing linear interpolation, but I haven't actually tested this theory! :)
+  */
+ static void *vol_precache_part(void *data)
+ {
+       VolPrecachePart *pa =  (VolPrecachePart *)data;
+       ObjectInstanceRen *obi = pa->obi;
 -static void precache_init_parts(Render *re, RayTree *tree, ShadeInput *shi, ObjectInstanceRen *obi, int totthread, int *parts)
++      RayObject *tree = pa->tree;
+       ShadeInput *shi = pa->shi;
+       float density, scatter_col[3] = {0.f, 0.f, 0.f};
+       float co[3];
+       int x, y, z;
+       const int res[3]= {pa->res[0], pa->res[1], pa->res[2]};
+       const float stepsize = vol_get_stepsize(shi, STEPSIZE_VIEW);
+       for (z= pa->minz; z < pa->maxz; z++) {
+               co[2] = pa->bbmin[2] + (pa->voxel[2] * (z + 0.5f));
+               
+               for (y= pa->miny; y < pa->maxy; y++) {
+                       co[1] = pa->bbmin[1] + (pa->voxel[1] * (y + 0.5f));
+                       
+                       for (x=pa->minx; x < pa->maxx; x++) {
+                               co[0] = pa->bbmin[0] + (pa->voxel[0] * (x + 0.5f));
+                               
+                               // don't bother if the point is not inside the volume mesh
+                               if (!point_inside_obi(tree, obi, co)) {
+                                       obi->volume_precache->data_r[ V_I(x, y, z, res) ] = -1.0f;
+                                       obi->volume_precache->data_g[ V_I(x, y, z, res) ] = -1.0f;
+                                       obi->volume_precache->data_b[ V_I(x, y, z, res) ] = -1.0f;
+                                       continue;
+                               }
+                               
+                               VecCopyf(shi->view, co);
+                               Normalize(shi->view);
+                               density = vol_get_density(shi, co);
+                               vol_get_scattering(shi, scatter_col, co, stepsize, density);
+                       
+                               obi->volume_precache->data_r[ V_I(x, y, z, res) ] = scatter_col[0];
+                               obi->volume_precache->data_g[ V_I(x, y, z, res) ] = scatter_col[1];
+                               obi->volume_precache->data_b[ V_I(x, y, z, res) ] = scatter_col[2];
+                       }
+               }
+       }
+       
+       pa->done = 1;
+       
+       return 0;
+ }
+ static void precache_setup_shadeinput(Render *re, ObjectInstanceRen *obi, Material *ma, ShadeInput *shi)
+ {
+       memset(shi, 0, sizeof(ShadeInput)); 
+       shi->depth= 1;
+       shi->mask= 1;
+       shi->mat = ma;
+       shi->vlr = NULL;
+       memcpy(&shi->r, &shi->mat->r, 23*sizeof(float));        // note, keep this synced with render_types.h
+       shi->har= shi->mat->har;
+       shi->obi= obi;
+       shi->obr= obi->obr;
+       shi->lay = re->scene->lay;
+ }
 -      RayTree *tree;
++static void precache_init_parts(Render *re, RayObject *tree, ShadeInput *shi, ObjectInstanceRen *obi, int totthread, int *parts)
+ {
+       VolumePrecache *vp = obi->volume_precache;
+       int i=0, x, y, z;
+       float voxel[3];
+       int sizex, sizey, sizez;
+       float *bbmin=obi->obr->boundbox[0], *bbmax=obi->obr->boundbox[1];
+       int *res;
+       int minx, maxx;
+       int miny, maxy;
+       int minz, maxz;
+       
+       if (!vp) return;
+       BLI_freelistN(&re->volume_precache_parts);
+       
+       /* currently we just subdivide the box, number of threads per side */
+       parts[0] = parts[1] = parts[2] = totthread;
+       res = vp->res;
+       
+       VecSubf(voxel, bbmax, bbmin);
+       
+       voxel[0] /= res[0];
+       voxel[1] /= res[1];
+       voxel[2] /= res[2];
+       for (x=0; x < parts[0]; x++) {
+               sizex = ceil(res[0] / (float)parts[0]);
+               minx = x * sizex;
+               maxx = minx + sizex;
+               maxx = (maxx>res[0])?res[0]:maxx;
+               
+               for (y=0; y < parts[1]; y++) {
+                       sizey = ceil(res[1] / (float)parts[1]);
+                       miny = y * sizey;
+                       maxy = miny + sizey;
+                       maxy = (maxy>res[1])?res[1]:maxy;
+                       
+                       for (z=0; z < parts[2]; z++) {
+                               VolPrecachePart *pa= MEM_callocN(sizeof(VolPrecachePart), "new precache part");
+                               
+                               sizez = ceil(res[2] / (float)parts[2]);
+                               minz = z * sizez;
+                               maxz = minz + sizez;
+                               maxz = (maxz>res[2])?res[2]:maxz;
+                                               
+                               pa->done = 0;
+                               pa->working = 0;
+                               
+                               pa->num = i;
+                               pa->tree = tree;
+                               pa->shi = shi;
+                               pa->obi = obi;
+                               VECCOPY(pa->bbmin, bbmin);
+                               VECCOPY(pa->voxel, voxel);
+                               VECCOPY(pa->res, res);
+                               
+                               pa->minx = minx; pa->maxx = maxx;
+                               pa->miny = miny; pa->maxy = maxy;
+                               pa->minz = minz; pa->maxz = maxz;
+                               
+                               
+                               BLI_addtail(&re->volume_precache_parts, pa);
+                               
+                               i++;
+                       }
+               }
+       }
+ }
+ static VolPrecachePart *precache_get_new_part(Render *re)
+ {
+       VolPrecachePart *pa, *nextpa=NULL;
+       
+       for (pa = re->volume_precache_parts.first; pa; pa=pa->next)
+       {
+               if (pa->done==0 && pa->working==0) {
+                       nextpa = pa;
+                       break;
+               }
+       }
+       return nextpa;
+ }
+ static int precache_resolution(VolumePrecache *vp, float *bbmin, float *bbmax, int res)
+ {
+       float dim[3], div;
+       
+       VecSubf(dim, bbmax, bbmin);
+       
+       div = MAX3(dim[0], dim[1], dim[2]);
+       dim[0] /= div;
+       dim[1] /= div;
+       dim[2] /= div;
+       
+       vp->res[0] = dim[0] * (float)res;
+       vp->res[1] = dim[1] * (float)res;
+       vp->res[2] = dim[2] * (float)res;
+       
+       if ((vp->res[0] < 1) || (vp->res[1] < 1) || (vp->res[2] < 1))
+               return 0;
+       
+       return 1;
+ }
+ /* Precache a volume into a 3D voxel grid.
+  * The voxel grid is stored in the ObjectInstanceRen, 
+  * in camera space, aligned with the ObjectRen's bounding box.
+  * Resolution is defined by the user.
+  */
+ void vol_precache_objectinstance_threads(Render *re, ObjectInstanceRen *obi, Material *ma)
+ {
+       VolumePrecache *vp;
+       VolPrecachePart *nextpa, *pa;
 -      tree = create_raytree_obi(obi, bbmin, bbmax);
++      RayObject *tree;
+       ShadeInput shi;
+       ListBase threads;
+       float *bbmin=obi->obr->boundbox[0], *bbmax=obi->obr->boundbox[1];
+       int parts[3], totparts;
+       
+       int caching=1, counter=0;
+       int totthread = re->r.threads;
+       
+       double time, lasttime= PIL_check_seconds_timer();
+       
+       R = *re;
+       /* create a raytree with just the faces of the instanced ObjectRen, 
+        * used for checking if the cached point is inside or outside. */
 -              RE_ray_tree_free(tree);
 -              tree= NULL;
++      //tree = create_raytree_obi(obi, bbmin, bbmax);
++      tree = makeraytree_object(&R, obi);
+       if (!tree) return;
++      INIT_MINMAX(bbmin, bbmax);
++      RE_rayobject_merge_bb( tree, bbmin, bbmax);
+       vp = MEM_callocN(sizeof(VolumePrecache), "volume light cache");
+       
+       if (!precache_resolution(vp, bbmin, bbmax, ma->vol.precache_resolution)) {
+               MEM_freeN(vp);
+               vp = NULL;
+               return;
+       }
+       vp->data_r = MEM_callocN(sizeof(float)*vp->res[0]*vp->res[1]*vp->res[2], "volume light cache data red channel");
+       vp->data_g = MEM_callocN(sizeof(float)*vp->res[0]*vp->res[1]*vp->res[2], "volume light cache data green channel");
+       vp->data_b = MEM_callocN(sizeof(float)*vp->res[0]*vp->res[1]*vp->res[2], "volume light cache data blue channel");
+       obi->volume_precache = vp;
+       /* Need a shadeinput to calculate scattering */
+       precache_setup_shadeinput(re, obi, ma, &shi);
+       
+       precache_init_parts(re, tree, &shi, obi, totthread, parts);
+       totparts = parts[0] * parts[1] * parts[2];
+       
+       BLI_init_threads(&threads, vol_precache_part, totthread);
+       
+       while(caching) {
+               if(BLI_available_threads(&threads) && !(re->test_break(re->tbh))) {
+                       nextpa = precache_get_new_part(re);
+                       if (nextpa) {
+                               nextpa->working = 1;
+                               BLI_insert_thread(&threads, nextpa);
+                       }
+               }
+               else PIL_sleep_ms(50);
+               caching=0;
+               counter=0;
+               for(pa= re->volume_precache_parts.first; pa; pa= pa->next) {
+                       
+                       if(pa->done) {
+                               counter++;
+                               BLI_remove_thread(&threads, pa);
+                       } else
+                               caching = 1;
+               }
+               
+               if (re->test_break(re->tbh) && BLI_available_threads(&threads)==totthread)
+                       caching=0;
+               
+               time= PIL_check_seconds_timer();
+               if(time-lasttime>1.0f) {
+                       char str[64];
+                       sprintf(str, "Precaching volume: %d%%", (int)(100.0f * ((float)counter / (float)totparts)));
+                       re->i.infostr= str;
+                       re->stats_draw(re->sdh, &re->i);
+                       re->i.infostr= NULL;
+                       lasttime= time;
+               }
+       }
+       
+       BLI_end_threads(&threads);
+       BLI_freelistN(&re->volume_precache_parts);
+       
+       if(tree) {
 -      RayTree *tree;
++              //TODO: makeraytree_object creates a tree and saves it on OBI, if we free this tree we should also clear other pointers to it
++              //RE_rayobject_free(tree);
++              //tree= NULL;
+       }
+       
+       lightcache_filter(obi->volume_precache);
+       
+       if (ELEM(ma->vol.shade_type, MA_VOL_SHADE_MULTIPLE, MA_VOL_SHADE_SINGLEPLUSMULTIPLE))
+       {
+               multiple_scattering_diffusion(re, vp, ma);
+       }
+ }
+ static int using_lightcache(Material *ma)
+ {
+       return (((ma->vol.shadeflag & MA_VOL_PRECACHESHADING) && (ma->vol.shade_type == MA_VOL_SHADE_SINGLE))
+               || (ELEM(ma->vol.shade_type, MA_VOL_SHADE_MULTIPLE, MA_VOL_SHADE_SINGLEPLUSMULTIPLE)));
+ }
+ /* loop through all objects (and their associated materials)
+  * marked for pre-caching in convertblender.c, and pre-cache them */
+ void volume_precache(Render *re)
+ {
+       ObjectInstanceRen *obi;
+       VolumeOb *vo;
+       for(vo= re->volumes.first; vo; vo= vo->next) {
+               if (using_lightcache(vo->ma)) {
+                       for(obi= re->instancetable.first; obi; obi= obi->next) {
+                               if (obi->obr == vo->obr) {
+                                       vol_precache_objectinstance_threads(re, obi, vo->ma);
+                               }
+                       }
+               }
+       }
+       
+       re->i.infostr= NULL;
+       re->stats_draw(re->sdh, &re->i);
+ }
+ void free_volume_precache(Render *re)
+ {
+       ObjectInstanceRen *obi;
+       
+       for(obi= re->instancetable.first; obi; obi= obi->next) {
+               if (obi->volume_precache != NULL) {
+                       MEM_freeN(obi->volume_precache->data_r);
+                       MEM_freeN(obi->volume_precache->data_g);
+                       MEM_freeN(obi->volume_precache->data_b);
+                       MEM_freeN(obi->volume_precache);
+                       obi->volume_precache = NULL;
+               }
+       }
+       
+       BLI_freelistN(&re->volumes);
+ }
+ int point_inside_volume_objectinstance(ObjectInstanceRen *obi, float *co)
+ {
 -      tree = create_raytree_obi(obi, obi->obr->boundbox[0], obi->obr->boundbox[1]);
++      RayObject *tree;
+       int inside=0;
+       
 -      RE_ray_tree_free(tree);
 -      tree= NULL;
++      tree = makeraytree_object(&R, obi); //create_raytree_obi(obi, obi->obr->boundbox[0], obi->obr->boundbox[1]);
+       if (!tree) return 0;
+       
+       inside = point_inside_obi(tree, obi, co);
+       
++      //TODO: makeraytree_object creates a tree and saves it on OBI, if we free this tree we should also clear other pointers to it
++      //RE_rayobject_free(tree);
++      //tree= NULL;
+       
+       return inside;
+ }
index 0000000,bc425c8..2808d26
mode 000000,100644..100644
--- /dev/null
@@@ -1,0 -1,752 +1,749 @@@
 -      float maxsize = RE_ray_tree_max_size(R.raytree);
 -      
+ /**
+  *
+  * ***** BEGIN GPL LICENSE BLOCK *****
+  *
+  * This program is free software; you can redistribute it and/or
+  * modify it under the terms of the GNU General Public License
+  * as published by the Free Software Foundation; either version 2
+  * of the License, or (at your option) any later version.
+  *
+  * This program is distributed in the hope that it will be useful,
+  * but WITHOUT ANY WARRANTY; without even the implied warranty of
+  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  * GNU General Public License for more details.
+  *
+  * You should have received a copy of the GNU General Public License
+  * along with this program; if not, write to the Free Software Foundation,
+  * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+  *
+  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
+  * All rights reserved.
+  *
+  * The Original Code is: all of this file.
+  *
+  * Contributor(s): Matt Ebb, Raul Fernandez Hernandez (Farsthary)
+  *
+  * ***** END GPL LICENSE BLOCK *****
+  */
+ #include <math.h>
+ #include <stdlib.h>
+ #include <string.h>
+ #include <float.h>
+ #include "MEM_guardedalloc.h"
+ #include "BLI_blenlib.h"
+ #include "BLI_arithb.h"
+ #include "BLI_rand.h"
+ #include "BLI_voxel.h"
+ #include "RE_shader_ext.h"
+ #include "RE_raytrace.h"
+ #include "DNA_material_types.h"
+ #include "DNA_group_types.h"
+ #include "DNA_lamp_types.h"
+ #include "DNA_meta_types.h"
+ #include "BKE_global.h"
+ #include "render_types.h"
+ #include "pixelshading.h"
+ #include "shading.h"
+ #include "texture.h"
+ #include "volumetric.h"
+ #include "volume_precache.h"
+ #if defined( _MSC_VER ) && !defined( __cplusplus )
+ # define inline __inline
+ #endif // defined( _MSC_VER ) && !defined( __cplusplus )
+ /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+ /* defined in pipeline.c, is hardcopy of active dynamic allocated Render */
+ /* only to be used here in this file, it's for speed */
+ extern struct Render R;
+ /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+ /* luminance rec. 709 */
+ inline float luminance(float* col)
+ {
+       return (0.212671f*col[0] + 0.71516f*col[1] + 0.072169f*col[2]);
+ }
+ /* tracing */
+ static int vol_get_bounds(ShadeInput *shi, float *co, float *vec, float *hitco, Isect *isect, int intersect_type)
+ {
 -      isect->oborig= RAY_OBJECT_SET(&R, shi->obi);
 -      isect->face_last= NULL;
 -      isect->ob_last= 0;
+       /* XXX TODO - get raytrace max distance from object instance's bounding box */
+       /* need to account for scaling only, but keep coords in camera space...
+        * below code is WIP and doesn't work!
+       VecSubf(bb_dim, shi->obi->obr->boundbox[1], shi->obi->obr->boundbox[2]);
+       Mat3MulVecfl(shi->obi->nmat, bb_dim);
+       maxsize = VecLength(bb_dim);
+       */
+       
+       VECCOPY(isect->start, co);
++      VECCOPY(isect->vec, vec );
++      isect->labda = FLT_MAX;
++      /*
+       isect->end[0] = co[0] + vec[0] * maxsize;
+       isect->end[1] = co[1] + vec[1] * maxsize;
+       isect->end[2] = co[2] + vec[2] * maxsize;
++      */
+       
+       isect->mode= RE_RAY_MIRROR;
 -      if (intersect_type == VOL_BOUNDS_DEPTH) isect->faceorig= (RayFace*)shi->vlr;
 -      else if (intersect_type == VOL_BOUNDS_SS) isect->faceorig= NULL;
++      isect->orig.ob = (void*)shi->obi;
++      isect->last_hit = NULL;
+       isect->lay= -1;
+       
 -      if(RE_ray_tree_intersect(R.raytree, isect))
++      if (intersect_type == VOL_BOUNDS_DEPTH) isect->orig.face = (void*)shi->vlr;
++      else if (intersect_type == VOL_BOUNDS_SS) isect->orig.face= NULL;
+       
 -      float maxsize = RE_ray_tree_max_size(R.raytree);
++      if(RE_rayobject_raycast(R.raytree, isect))
+       {
+               hitco[0] = isect->start[0] + isect->labda*isect->vec[0];
+               hitco[1] = isect->start[1] + isect->labda*isect->vec[1];
+               hitco[2] = isect->start[2] + isect->labda*isect->vec[2];
+               return 1;
+       } else {
+               return 0;
+       }
+ }
+ static void shade_intersection(ShadeInput *shi, float *col, Isect *is)
+ {
+       ShadeInput shi_new;
+       ShadeResult shr_new;
+       
+       memset(&shi_new, 0, sizeof(ShadeInput)); 
+       
+       shi_new.mask= shi->mask;
+       shi_new.osatex= shi->osatex;
+       shi_new.thread= shi->thread;
+       shi_new.depth = shi->depth + 1;
+       shi_new.volume_depth= shi->volume_depth + 1;
+       shi_new.xs= shi->xs;
+       shi_new.ys= shi->ys;
+       shi_new.lay= shi->lay;
+       shi_new.passflag= SCE_PASS_COMBINED; /* result of tracing needs no pass info */
+       shi_new.combinedflag= 0xFFFFFF;          /* ray trace does all options */
+       shi_new.light_override= shi->light_override;
+       shi_new.mat_override= shi->mat_override;
+       
+       VECCOPY(shi_new.camera_co, is->start);
+       
+       memset(&shr_new, 0, sizeof(ShadeResult));
+       
+       /* hardcoded limit of 100 for now - prevents problems in weird geometry */
+       if (shi->volume_depth < 100) {
+               shade_ray(is, &shi_new, &shr_new);
+       }
+       
+       VecCopyf(col, shr_new.combined);
+       col[3] = shr_new.alpha;
+ }
+ static void vol_trace_behind(ShadeInput *shi, VlakRen *vlr, float *co, float *col)
+ {
+       Isect isect;
 -      isect.end[0] = isect.start[0] + shi->view[0] * maxsize;
 -      isect.end[1] = isect.start[1] + shi->view[1] * maxsize;
 -      isect.end[2] = isect.start[2] + shi->view[2] * maxsize;
 -      
 -      isect.faceorig= (RayFace *)vlr;
+       
+       VECCOPY(isect.start, co);
 -      isect.oborig= RAY_OBJECT_SET(&R, shi->obi);
 -      isect.face_last= NULL;
 -      isect.ob_last= 0;
++      VECCOPY(isect.vec, shi->view);
++      isect.labda = FLT_MAX;
+       
+       isect.mode= RE_RAY_MIRROR;
 -      if(RE_ray_tree_intersect(R.raytree, &isect)) {
++      isect.orig.ob = (void*) shi->obi;
++      isect.orig.face = (void*)vlr;
++      isect.last_hit = NULL;
+       isect.lay= -1;
+       
+       /* check to see if there's anything behind the volume, otherwise shade the sky */
 -                      VlakRen *vlr = (VlakRen *)is.face;
++      if(RE_rayobject_raycast(R.raytree, &isect)) {
+               shade_intersection(shi, col, &isect);
+       } else {
+               shadeSkyView(col, co, shi->view, NULL, shi->thread);
+               shadeSunView(col, shi->view);
+       }
+ }
+ /* input shader data */
+ float vol_get_stepsize(struct ShadeInput *shi, int context)
+ {
+       if (shi->mat->vol.stepsize_type == MA_VOL_STEP_RANDOMIZED) {
+               /* range between 0.75 and 1.25 */
+               const float rnd = 0.5f * BLI_thread_frand(shi->thread) + 0.75f;
+       
+               if (context == STEPSIZE_VIEW)
+                       return shi->mat->vol.stepsize * rnd;
+               else if (context == STEPSIZE_SHADE)
+                       return shi->mat->vol.shade_stepsize * rnd;
+       }
+       else {  // MA_VOL_STEP_CONSTANT
+               
+               if (context == STEPSIZE_VIEW)
+                       return shi->mat->vol.stepsize;
+               else if (context == STEPSIZE_SHADE)
+                       return shi->mat->vol.shade_stepsize;
+       }
+       
+       return shi->mat->vol.stepsize;
+ }
+ /* trilinear interpolation */
+ static void vol_get_precached_scattering(ShadeInput *shi, float *scatter_col, float *co)
+ {
+       VolumePrecache *vp = shi->obi->volume_precache;
+       float bbmin[3], bbmax[3], dim[3];
+       float sample_co[3];
+       
+       if (!vp) return;
+       
+       /* convert input coords to 0.0, 1.0 */
+       VECCOPY(bbmin, shi->obi->obr->boundbox[0]);
+       VECCOPY(bbmax, shi->obi->obr->boundbox[1]);
+       VecSubf(dim, bbmax, bbmin);
+       sample_co[0] = ((co[0] - bbmin[0]) / dim[0]);
+       sample_co[1] = ((co[1] - bbmin[1]) / dim[1]);
+       sample_co[2] = ((co[2] - bbmin[2]) / dim[2]);
+       scatter_col[0] = voxel_sample_trilinear(vp->data_r, vp->res, sample_co);
+       scatter_col[1] = voxel_sample_trilinear(vp->data_g, vp->res, sample_co);
+       scatter_col[2] = voxel_sample_trilinear(vp->data_b, vp->res, sample_co);
+ }
+ /* Meta object density, brute force for now 
+  * (might be good enough anyway, don't need huge number of metaobs to model volumetric objects */
+ static float metadensity(Object* ob, float* co)
+ {
+       float mat[4][4], imat[4][4], dens = 0.f;
+       MetaBall* mb = (MetaBall*)ob->data;
+       MetaElem* ml;
+       
+       /* transform co to meta-element */
+       float tco[3] = {co[0], co[1], co[2]};
+       Mat4MulMat4(mat, ob->obmat, R.viewmat);
+       Mat4Invert(imat, mat);
+       Mat4MulVecfl(imat, tco);
+       
+       for (ml = mb->elems.first; ml; ml=ml->next) {
+               float bmat[3][3], dist2;
+               
+               /* element rotation transform */
+               float tp[3] = {ml->x - tco[0], ml->y - tco[1], ml->z - tco[2]};
+               QuatToMat3(ml->quat, bmat);
+               Mat3Transp(bmat);       // rot.only, so inverse == transpose
+               Mat3MulVecfl(bmat, tp);
+               
+               /* MB_BALL default */
+               switch (ml->type) {
+                       case MB_ELIPSOID:
+                               tp[0] /= ml->expx, tp[1] /= ml->expy, tp[2] /= ml->expz;
+                               break;
+                       case MB_CUBE:
+                               tp[2] = (tp[2] > ml->expz) ? (tp[2] - ml->expz) : ((tp[2] < -ml->expz) ? (tp[2] + ml->expz) : 0.f);
+                               // no break, xy as plane
+                       case MB_PLANE:
+                               tp[1] = (tp[1] > ml->expy) ? (tp[1] - ml->expy) : ((tp[1] < -ml->expy) ? (tp[1] + ml->expy) : 0.f);
+                               // no break, x as tube
+                       case MB_TUBE:
+                               tp[0] = (tp[0] > ml->expx) ? (tp[0] - ml->expx) : ((tp[0] < -ml->expx) ? (tp[0] + ml->expx) : 0.f);
+               }
+               
+               /* ml->rad2 is not set */
+               dist2 = 1.f - ((tp[0]*tp[0] + tp[1]*tp[1] + tp[2]*tp[2]) / (ml->rad*ml->rad));
+               if (dist2 > 0.f)
+                       dens += (ml->flag & MB_NEGATIVE) ? -ml->s*dist2*dist2*dist2 : ml->s*dist2*dist2*dist2;
+       }
+       
+       dens -= mb->thresh;
+       return (dens < 0.f) ? 0.f : dens;
+ }
+ float vol_get_density(struct ShadeInput *shi, float *co)
+ {
+       float density = shi->mat->vol.density;
+       float density_scale = shi->mat->vol.density_scale;
+               
+       do_volume_tex(shi, co, MAP_DENSITY, NULL, &density);
+       
+       // if meta-object, modulate by metadensity without increasing it
+       if (shi->obi->obr->ob->type == OB_MBALL) {
+               const float md = metadensity(shi->obi->obr->ob, co);
+               if (md < 1.f) density *= md;
+       }
+       
+       return density * density_scale;
+ }
+ /* scattering multiplier, values above 1.0 are non-physical, 
+  * but can be useful to tweak lighting */
+ float vol_get_scattering_fac(ShadeInput *shi, float *co)
+ {
+       float scatter = shi->mat->vol.scattering;
+       float col[3] = {0.0, 0.0, 0.0};
+       
+       do_volume_tex(shi, co, MAP_SCATTERING, col, &scatter);
+       
+       return scatter;
+ }
+ /* compute emission component, amount of radiance to add per segment
+  * can be textured with 'emit' */
+ void vol_get_emission(ShadeInput *shi, float *emission_col, float *co, float density)
+ {
+       float emission = shi->mat->vol.emission;
+       VECCOPY(emission_col, shi->mat->vol.emission_col);
+       
+       do_volume_tex(shi, co, MAP_EMISSION+MAP_EMISSION_COL, emission_col, &emission);
+       
+       emission_col[0] = emission_col[0] * emission * density;
+       emission_col[1] = emission_col[1] * emission * density;
+       emission_col[2] = emission_col[2] * emission * density;
+ }
+ void vol_get_absorption(ShadeInput *shi, float *absorb_col, float *co)
+ {
+       float absorption = shi->mat->vol.absorption;
+       VECCOPY(absorb_col, shi->mat->vol.absorption_col);
+       
+       do_volume_tex(shi, co, MAP_ABSORPTION+MAP_ABSORPTION_COL, absorb_col, &absorption);
+       
+       absorb_col[0] = (1.0f - absorb_col[0]) * absorption;
+       absorb_col[1] = (1.0f - absorb_col[1]) * absorption;
+       absorb_col[2] = (1.0f - absorb_col[2]) * absorption;
+ }
+ /* phase function - determines in which directions the light 
+  * is scattered in the volume relative to incoming direction 
+  * and view direction */
+ float vol_get_phasefunc(ShadeInput *shi, short phasefunc_type, float g, float *w, float *wp)
+ {
+       const float costheta = Inpf(w, wp);
+       const float scale = M_PI;
+       
+       /*
+        * Scale constant is required, since Blender's shading system doesn't normalise for
+        * energy conservation - eg. scaling by 1/pi for a lambert shader.
+        * This makes volumes darker than other solid objects, for the same lighting intensity.
+        * To correct this, scale up the phase function values
+        * until Blender's shading system supports this better. --matt
+        */
+       
+       switch (phasefunc_type) {
+               case MA_VOL_PH_MIEHAZY:
+                       return scale * (0.5f + 4.5f * powf(0.5 * (1.f + costheta), 8.f)) / (4.f*M_PI);
+               case MA_VOL_PH_MIEMURKY:
+                       return scale * (0.5f + 16.5f * powf(0.5 * (1.f + costheta), 32.f)) / (4.f*M_PI);
+               case MA_VOL_PH_RAYLEIGH:
+                       return scale * 3.f/(16.f*M_PI) * (1 + costheta * costheta);
+               case MA_VOL_PH_HG:
+                       return scale * (1.f / (4.f * M_PI) * (1.f - g*g) / powf(1.f + g*g - 2.f * g * costheta, 1.5f));
+               case MA_VOL_PH_SCHLICK:
+               {
+                       const float k = 1.55f * g - .55f * g * g * g;
+                       const float kcostheta = k * costheta;
+                       return scale * (1.f / (4.f * M_PI) * (1.f - k*k) / ((1.f - kcostheta) * (1.f - kcostheta)));
+               }
+               case MA_VOL_PH_ISOTROPIC:
+               default:
+                       return scale * (1.f / (4.f * M_PI));
+       }
+ }
+ /* Compute transmittance = e^(-attenuation) */
+ void vol_get_transmittance_seg(ShadeInput *shi, float *tr, float stepsize, float *co, float density)
+ {
+       /* input density = density at co */
+       float tau[3] = {0.f, 0.f, 0.f};
+       float absorb[3];
+       const float scatter_dens = vol_get_scattering_fac(shi, co) * density * stepsize;
+       vol_get_absorption(shi, absorb, co);
+       
+       /* homogenous volume within the sampled distance */
+       tau[0] += scatter_dens * absorb[0];
+       tau[1] += scatter_dens * absorb[1];
+       tau[2] += scatter_dens * absorb[2];
+       
+       tr[0] *= exp(-tau[0]);
+       tr[1] *= exp(-tau[1]);
+       tr[2] *= exp(-tau[2]);
+ }
+ /* Compute transmittance = e^(-attenuation) */
+ static void vol_get_transmittance(ShadeInput *shi, float *tr, float *co, float *endco)
+ {
+       float p[3] = {co[0], co[1], co[2]};
+       float step_vec[3] = {endco[0] - co[0], endco[1] - co[1], endco[2] - co[2]};
+       //const float ambtau = -logf(shi->mat->vol.depth_cutoff);       // never zero
+       float tau[3] = {0.f, 0.f, 0.f};
+       float t0 = 0.f;
+       float t1 = Normalize(step_vec);
+       float pt0 = t0;
+       
+       t0 += shi->mat->vol.shade_stepsize * ((shi->mat->vol.stepsize_type == MA_VOL_STEP_CONSTANT) ? 0.5f : BLI_thread_frand(shi->thread));
+       p[0] += t0 * step_vec[0];
+       p[1] += t0 * step_vec[1];
+       p[2] += t0 * step_vec[2];
+       VecMulf(step_vec, shi->mat->vol.shade_stepsize);
+       for (; t0 < t1; pt0 = t0, t0 += shi->mat->vol.shade_stepsize) {
+               float absorb[3];
+               const float d = vol_get_density(shi, p);
+               const float stepd = (t0 - pt0) * d;
+               const float scatter_dens = vol_get_scattering_fac(shi, p) * stepd;
+               vol_get_absorption(shi, absorb, p);
+               
+               tau[0] += scatter_dens * absorb[0];
+               tau[1] += scatter_dens * absorb[1];
+               tau[2] += scatter_dens * absorb[2];
+               
+               //if (luminance(tau) >= ambtau) break;
+               VecAddf(p, p, step_vec);
+       }
+       
+       /* return transmittance */
+       tr[0] = expf(-tau[0]);
+       tr[1] = expf(-tau[1]);
+       tr[2] = expf(-tau[2]);
+ }
+ void vol_shade_one_lamp(struct ShadeInput *shi, float *co, LampRen *lar, float *lacol)
+ {
+       float visifac, lv[3], lampdist;
+       float tr[3]={1.0,1.0,1.0};
+       float hitco[3], *atten_co;
+       float p;
+       float scatter_fac;
+       float shade_stepsize = vol_get_stepsize(shi, STEPSIZE_SHADE);
+       
+       if (lar->mode & LA_LAYER) if((lar->lay & shi->obi->lay)==0) return;
+       if ((lar->lay & shi->lay)==0) return;
+       if (lar->energy == 0.0) return;
+       
+       if ((visifac= lamp_get_visibility(lar, co, lv, &lampdist)) == 0.f) return;
+       
+       VecCopyf(lacol, &lar->r);
+       
+       if(lar->mode & LA_TEXTURE) {
+               shi->osatex= 0;
+               do_lamp_tex(lar, lv, shi, lacol, LA_TEXTURE);
+       }
+       VecMulf(lacol, visifac);
+       if (ELEM(lar->type, LA_SUN, LA_HEMI))
+               VECCOPY(lv, lar->vec);
+       VecMulf(lv, -1.0f);
+       
+       if (shi->mat->vol.shade_type != MA_VOL_SHADE_NONE) {
+               Isect is;
+               
+               /* find minimum of volume bounds, or lamp coord */
+               if (vol_get_bounds(shi, co, lv, hitco, &is, VOL_BOUNDS_SS)) {
+                       float dist = VecLenf(co, hitco);
 -              VlakRen *vlr = (VlakRen *)is.face;
++                      VlakRen *vlr = (VlakRen *)is.hit.face;
+                       
+                       /* simple internal shadowing */
+                       if (vlr->mat->material_type == MA_TYPE_SURFACE) {
+                               lacol[0] = lacol[1] = lacol[2] = 0.0f;
+                               return;
+                       }
+                       if (ELEM(lar->type, LA_SUN, LA_HEMI))
+                               /* infinite lights, can never be inside volume */
+                               atten_co = hitco;
+                       else if ( lampdist < dist ) {
+                               atten_co = lar->co;
+                       } else
+                               atten_co = hitco;
+                       
+                       vol_get_transmittance(shi, tr, co, atten_co);
+                       
+                       VecMulVecf(lacol, lacol, tr);
+               }
+               else {
+                       /* Point is on the outside edge of the volume,
+                        * therefore no attenuation, full transmission.
+                        * Radiance from lamp remains unchanged */
+               }
+       }
+       
+       p = vol_get_phasefunc(shi, shi->mat->vol.phasefunc_type, shi->mat->vol.phasefunc_g, shi->view, lv);
+       VecMulf(lacol, p);
+       
+       scatter_fac = vol_get_scattering_fac(shi, co);
+       VecMulf(lacol, scatter_fac);
+ }
+ /* single scattering only for now */
+ void vol_get_scattering(ShadeInput *shi, float *scatter_col, float *co, float stepsize, float density)
+ {
+       ListBase *lights;
+       GroupObject *go;
+       LampRen *lar;
+       
+       scatter_col[0] = scatter_col[1] = scatter_col[2] = 0.f;
+       
+       lights= get_lights(shi);
+       for(go=lights->first; go; go= go->next)
+       {
+               float lacol[3] = {0.f, 0.f, 0.f};
+               lar= go->lampren;
+               
+               if (lar) {
+                       vol_shade_one_lamp(shi, co, lar, lacol);
+                       VecAddf(scatter_col, scatter_col, lacol);
+               }
+       }
+ }
+       
+ /*
+ The main volumetric integrator, using an emission/absorption/scattering model.
+ Incoming radiance = 
+ outgoing radiance from behind surface * beam transmittance/attenuation
+ + added radiance from all points along the ray due to participating media
+       --> radiance for each segment = 
+               (radiance added by scattering + radiance added by emission) * beam transmittance/attenuation
+ */
+ static void volumeintegrate(struct ShadeInput *shi, float *col, float *co, float *endco)
+ {
+       float tr[3] = {1.0f, 1.0f, 1.0f};
+       float radiance[3] = {0.f, 0.f, 0.f}, d_radiance[3] = {0.f, 0.f, 0.f};
+       float stepsize = vol_get_stepsize(shi, STEPSIZE_VIEW);
+       int nsteps, s;
+       float emit_col[3], scatter_col[3] = {0.0, 0.0, 0.0};
+       float stepvec[3], step_sta[3], step_end[3], step_mid[3];
+       float density;
+       const float depth_cutoff = shi->mat->vol.depth_cutoff;
+       /* ray marching */
+       nsteps = (int)((VecLenf(co, endco) / stepsize) + 0.5);
+       
+       VecSubf(stepvec, endco, co);
+       VecMulf(stepvec, 1.0f / nsteps);
+       VecCopyf(step_sta, co);
+       VecAddf(step_end, step_sta, stepvec);
+       
+       /* get radiance from all points along the ray due to participating media */
+       for (s = 0; s < nsteps; s++) {
+               density = vol_get_density(shi, step_sta);
+               
+               /* there's only any use in shading here if there's actually some density to shade! */
+               if (density > 0.01f) {
+               
+                       /* transmittance component (alpha) */
+                       vol_get_transmittance_seg(shi, tr, stepsize, co, density);
+                       step_mid[0] = step_sta[0] + (stepvec[0] * 0.5);
+                       step_mid[1] = step_sta[1] + (stepvec[1] * 0.5);
+                       step_mid[2] = step_sta[2] + (stepvec[2] * 0.5);
+               
+                       /* incoming light via emission or scattering (additive) */
+                       vol_get_emission(shi, emit_col, step_mid, density);
+                       
+                       if (shi->obi->volume_precache)
+                               vol_get_precached_scattering(shi, scatter_col, step_mid);
+                       else
+                               vol_get_scattering(shi, scatter_col, step_mid, stepsize, density);
+                       
+                       VecMulf(scatter_col, density);
+                       VecAddf(d_radiance, emit_col, scatter_col);
+                       
+                       /*   Lv += Tr * (Lve() + Ld) */
+                       VecMulVecf(d_radiance, tr, d_radiance);
+                       VecMulf(d_radiance, stepsize);
+                       
+                       VecAddf(radiance, radiance, d_radiance);        
+               }
+               VecCopyf(step_sta, step_end);
+               VecAddf(step_end, step_end, stepvec);
+               
+               /* luminance rec. 709 */
+               if ((0.2126*tr[0] + 0.7152*tr[1] + 0.0722*tr[2]) < depth_cutoff) break; 
+       }
+       
+       /* multiply original color (behind volume) with beam transmittance over entire distance */
+       VecMulVecf(col, tr, col);       
+       VecAddf(col, col, radiance);
+       
+       /* alpha <-- transmission luminance */
+       col[3] = 1.0f -(0.2126*tr[0] + 0.7152*tr[1] + 0.0722*tr[2]);
+ }
+ /* the main entry point for volume shading */
+ static void volume_trace(struct ShadeInput *shi, struct ShadeResult *shr, int inside_volume)
+ {
+       float hitco[3], col[4] = {0.f,0.f,0.f,0.f};
+       float *startco, *endco;
+       int trace_behind = 1;
+       const int ztransp= ((shi->depth==0) && (shi->mat->mode & MA_TRANSP) && (shi->mat->mode & MA_ZTRANSP));
+       Isect is;
+       /* check for shading an internal face a volume object directly */
+       if (inside_volume == VOL_SHADE_INSIDE)
+               trace_behind = 0;
+       else if (inside_volume == VOL_SHADE_OUTSIDE) {
+               if (shi->flippednor)
+                       inside_volume = VOL_SHADE_INSIDE;
+       }
+       
+       if (ztransp && inside_volume == VOL_SHADE_INSIDE) {
+               MatInside *mi;
+               int render_this=0;
+               
+               /* don't render the backfaces of ztransp volume materials.
+                
+                * volume shading renders the internal volume from between the
+                * near view intersection of the solid volume to the
+                * intersection on the other side, as part of the shading of
+                * the front face.
+                
+                * Because ztransp renders both front and back faces independently
+                * this will double up, so here we prevent rendering the backface as well, 
+                * which would otherwise render the volume in between the camera and the backface
+                * --matt */
+               
+               for (mi=R.render_volumes_inside.first; mi; mi=mi->next) {
+                       /* weak... */
+                       if (mi->ma == shi->mat) render_this=1;
+               }
+               if (!render_this) return;
+       }
+       
+       if (inside_volume == VOL_SHADE_INSIDE)
+       {
+               startco = shi->camera_co;
+               endco = shi->co;
+               
+               if (trace_behind) {
+                       if (!ztransp)
+                               /* trace behind the volume object */
+                               vol_trace_behind(shi, shi->vlr, endco, col);
+               } else {
+                       /* we're tracing through the volume between the camera 
+                        * and a solid surface, so use that pre-shaded radiance */
+                       QUATCOPY(col, shr->combined);
+               }
+               
+               /* shade volume from 'camera' to 1st hit point */
+               volumeintegrate(shi, col, startco, endco);
+       }
+       /* trace to find a backface, the other side bounds of the volume */
+       /* (ray intersect ignores front faces here) */
+       else if (vol_get_bounds(shi, shi->co, shi->view, hitco, &is, VOL_BOUNDS_DEPTH))
+       {
 -                              vol_trace_behind(shi, (VlakRen *)is.face, endco, col);
++              VlakRen *vlr = (VlakRen *)is.hit.face;
+               
+               startco = shi->co;
+               endco = hitco;
+               
+               if (!ztransp) {
+                       /* if it's another face in the same material */
+                       if (vlr->mat == shi->mat) {
+                               /* trace behind the 2nd (raytrace) hit point */
++                              vol_trace_behind(shi, (VlakRen *)is.hit.face, endco, col);
+                       } else {
+                               shade_intersection(shi, col, &is);
+                       }
+               }
+               
+               /* shade volume from 1st hit point to 2nd hit point */
+               volumeintegrate(shi, col, startco, endco);
+       }
+       
+       if (ztransp)
+               col[3] = col[3]>1.f?1.f:col[3];
+       else
+               col[3] = 1.f;
+       
+       VecCopyf(shr->combined, col);
+       shr->alpha = col[3];
+       
+       VECCOPY(shr->diff, shr->combined);
+ }
+ /* Traces a shadow through the object, 
+  * pretty much gets the transmission over a ray path */
+ void shade_volume_shadow(struct ShadeInput *shi, struct ShadeResult *shr, struct Isect *last_is)
+ {
+       float hitco[3];
+       float tr[3] = {1.0,1.0,1.0};
+       Isect is;
+       float shade_stepsize = vol_get_stepsize(shi, STEPSIZE_SHADE);
+       float *startco, *endco;
+       float density=0.f;
+       memset(shr, 0, sizeof(ShadeResult));
+       
+       /* if 1st hit normal is facing away from the camera, 
+        * then we're inside the volume already. */
+       if (shi->flippednor) {
+               startco = last_is->start;
+               endco = shi->co;
+       }
+       /* trace to find a backface, the other side bounds of the volume */
+       /* (ray intersect ignores front faces here) */
+       else if (vol_get_bounds(shi, shi->co, shi->view, hitco, &is, VOL_BOUNDS_DEPTH)) {
+               startco = shi->co;
+               endco = hitco;
+       }
+       else {
+               shr->combined[0] = shr->combined[1] = shr->combined[2] = 0.f;
+               shr->alpha = shr->combined[3] = 1.f;
+               return;
+       }
+       
+       density = vol_get_density(shi, startco);
+       vol_get_transmittance(shi, tr, startco, endco);
+       
+       VecCopyf(shr->combined, tr);
+       shr->combined[3] = 1.0f -(0.2126*tr[0] + 0.7152*tr[1] + 0.0722*tr[2]);
+       shr->alpha = shr->combined[3];
+ }
+ /* delivers a fully filled in ShadeResult, for all passes */
+ void shade_volume_outside(ShadeInput *shi, ShadeResult *shr)
+ {
+       memset(shr, 0, sizeof(ShadeResult));
+       volume_trace(shi, shr, VOL_SHADE_OUTSIDE);
+ }
+ void shade_volume_inside(ShadeInput *shi, ShadeResult *shr)
+ {
+       MatInside *m;
+       Material *mat_backup;
+       ObjectInstanceRen *obi_backup;
+       float prev_alpha = shr->alpha;
+       
+       //if (BLI_countlist(&R.render_volumes_inside) == 0) return;
+       
+       /* XXX: extend to multiple volumes perhaps later */
+       mat_backup = shi->mat;
+       obi_backup = shi->obi;
+       
+       m = R.render_volumes_inside.first;
+       shi->mat = m->ma;
+       shi->obi = m->obi;
+       shi->obr = m->obi->obr;
+       
+       volume_trace(shi, shr, VOL_SHADE_INSIDE);
+       shr->alpha += prev_alpha;
+       CLAMP(shr->alpha, 0.f, 1.f);
+       
+       shi->mat = mat_backup;
+       shi->obi = obi_backup;
+       shi->obr = obi_backup->obr;
+ }