path: root/source/blender/render/intern/source/volume_precache.c

/**
 *
 * ***** 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 "PIL_time.h"

#include "RE_shader_ext.h"
#include "RE_raytrace.h"

#include "DNA_material_types.h"

#include "render_types.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;
/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */


/* 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 */
int intersect_outside_volume(RayTree *tree, Isect *isect, float *offset, int limit, int depth)
{
	if (limit == 0) return depth;
	
	if (RE_ray_tree_intersect(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 */
int point_inside_obi(RayTree *tree, ObjectInstanceRen *obi, float *co)
{
	float maxsize = RE_ray_tree_max_size(tree);
	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);
	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;
	isect.face_last= 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;
}

static int inside_check_func(Isect *is, int ob, RayFace *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;
}

RayTree *create_raytree_obi(ObjectInstanceRen *obi, float *bbmin, float *bbmax)
{
	int v;
	VlakRen *vlr= NULL;
	
	/* create empty raytree */
	RayTree *tree = RE_ray_tree_create(64, obi->obr->totvlak, bbmin, bbmax,
		vlr_face_coords, inside_check_func, NULL, NULL);
	
	/* 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;
}

static float get_avg_surrounds(float *cache, int res, int res_2, int res_3, int rgb, int xx, int yy, int zz)
{
	int x, y, z, x_, y_, z_;
	int added=0;
	float tot=0.0f;
	int i;
	
	for (x=-1; x <= 1; x++) {
		x_ = xx+x;
		if (x_ >= 0 && x_ <= res-1) {
		
			for (y=-1; y <= 1; y++) {
				y_ = yy+y;
				if (y_ >= 0 && y_ <= res-1) {
				
					for (z=-1; z <= 1; z++) {
						z_ = zz+z;
						if (z_ >= 0 && z_ <= res-1) {
						
							i = rgb*res_3 + x_*res_2 + y_*res + z_;
							if (cache[i] > 0.0f) {
								tot += cache[i];
								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(float *cache, int res)
{
	int x, y, z, rgb;
	int res_2, res_3;
	int i;
	
	res_2 = res*res;
	res_3 = res*res*res;

	for (x=0; x < res; x++) {
		for (y=0; y < res; y++) {
			for (z=0; z < res; z++) {
				for (rgb=0; rgb < 3; rgb++) {
					i = rgb*res_3 + x*res_2 + y*res + z;

					/* trigger for outside mesh */
					if (cache[i] < -0.5f) cache[i] = get_avg_surrounds(cache, res, res_2, res_3, rgb, x, y, z);
				}
			}
		}
	}
}

static inline int I(int x,int y,int z,int n) //has a pad of 1 voxel surrounding the core for boundary simulation
{ 
	return (x*(n+2)+y)*(n+2)+z;
}


/* get the total amount of light energy in the light cache. used to normalise after multiple scattering */
static float total_ss_energy(float *cache, const int res)
{
	int x, y, z, rgb;
	int res_2, res_3;
	int i;
	float energy=0.f;
	
	res_2 = res*res;
	res_3 = res*res*res;

	for (x=0; x < res; x++) {
		for (y=0; y < res; y++) {
			for (z=0; z < res; z++) {
				for (rgb=0; rgb < 3; rgb++) {
					i = rgb*res_3 + x*res_2 + y*res + z;
					
					if (cache[i] > 0.f) energy += cache[i];
				}
			}
		}
	}
	
	return energy;
}

static float total_ms_energy(float *sr, float *sg, float *sb, const int res)
{
	int x, y, z, i;
	float energy=0.f;
	
	for (z=1;z<=res;z++) {
		for (y=1;y<=res;y++) {
			for (x=1;x<=res;x++) {
			
				i = 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*n*n;
	
	for (l=0; l<20; l++)
	{
		for (k=1; k<=n; k++)
		{
			for (j=1; j<=n; j++)
			{
				for (i=1; i<=n; i++)
				{
					x[I(i,j,k,n)] = (x0[I(i,j,k,n)] + a*(
														 x[I(i-1,j,k,n)]+x[I(i+1,j,k,n)]+
														 x[I(i,j-1,k,n)]+x[I(i,j+1,k,n)]+
														 x[I(i,j,k-1,n)]+x[I(i,j,k+1,n)]))/(1+6*a);
				}
			}
		}
	}
}

void multiple_scattering_diffusion(Render *re, float *cache, int res, 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 i, j, k, m;
	int n = res;
	const int size = (n+2)*(n+2)*(n+2);
	double time, lasttime= PIL_check_seconds_timer();
	float total;
	float c=1.0f;
	int index;
	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*n*n*simframes);
	
	energy_ss = total_ss_energy(cache, res);
	
	/* Scattering as diffusion pass */
	for (m=0; m<simframes; m++)
	{
		/* add sources */
		for (k=1; k<=n; k++)
		{
			for (j=1; j<=n; j++)
			{
				for (i=1; i<=n; i++)
				{
					time= PIL_check_seconds_timer();
					c++;
					
					index=(i-1)*n*n + (j-1)*n + k-1;
					
					if (cache[index] > 0.0f)
						sr[I(i,j,k,n)] += cache[index];
					if (cache[1*n*n*n + index] > 0.0f)
						sg[I(i,j,k,n)] += cache[1*n*n*n + index];
					if (cache[2*n*n*n + index] > 0.0f)
						sb[I(i,j,k,n)] += cache[2*n*n*n + index];


					/* 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->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()) break;
	}
	
	/* normalisation factor to conserve energy */
	energy_ms = total_ms_energy(sr, sg, sb, res);
	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 (k=1;k<=n;k++)
	{
		for (j=1;j<=n;j++)
		{
			for (i=1;i<=n;i++)
			{
				index=(i-1)*n*n + (j-1)*n + k-1;
				cache[index]			= origf * cache[index]  + fac * sr[I(i,j,k,n)];
				cache[1*n*n*n + index]	= origf * cache[1*n*n*n + index] + fac * sg[I(i,j,k,n)];
				cache[2*n*n*n + index]	= origf * cache[2*n*n*n + index] + fac * sb[I(i,j,k,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;
	RayTree *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=pa->res, res_2=pa->res*pa->res, res_3=pa->res*pa->res*pa->res;
	const float stepsize = vol_get_stepsize(shi, STEPSIZE_VIEW);

	for (x= pa->minx; x < pa->maxx; x++) {
		co[0] = pa->bbmin[0] + (pa->voxel[0] * x);
		
		for (y= pa->miny; y < pa->maxy; y++) {
			co[1] = pa->bbmin[1] + (pa->voxel[1] * y);
			
			for (z=pa->minz; z < pa->maxz; z++) {
				co[2] = pa->bbmin[2] + (pa->voxel[2] * z);
				
				// don't bother if the point is not inside the volume mesh
				if (!point_inside_obi(tree, obi, co)) {
					obi->volume_precache[0*res_3 + x*res_2 + y*res + z] = -1.0f;
					obi->volume_precache[1*res_3 + x*res_2 + y*res + z] = -1.0f;
					obi->volume_precache[2*res_3 + x*res_2 + y*res + z] = -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[0*res_3 + x*res_2 + y*res + z] = scatter_col[0];
				obi->volume_precache[1*res_3 + x*res_2 + y*res + z] = scatter_col[1];
				obi->volume_precache[2*res_3 + x*res_2 + y*res + z] = 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;
}

static void precache_init_parts(Render *re, RayTree *tree, ShadeInput *shi, ObjectInstanceRen *obi, float *bbmin, float *bbmax, int res, int totthread, int *parts)
{
	int i=0, x, y, z;
	float voxel[3];
	int sizex, sizey, sizez;
	int minx, maxx;
	int miny, maxy;
	int minz, maxz;
	
	BLI_freelistN(&re->volume_precache_parts);
	
	/* currently we just subdivide the box, number of threads per side */
	parts[0] = parts[1] = parts[2] = totthread;
	
	VecSubf(voxel, bbmax, bbmin);
	if ((voxel[0] < FLT_EPSILON) || (voxel[1] < FLT_EPSILON) || (voxel[2] < FLT_EPSILON))
		return;
	VecMulf(voxel, 1.0f/res);

	for (x=0; x < parts[0]; x++) {
		sizex = ceil(res / (float)parts[0]);
		minx = x * sizex;
		maxx = minx + sizex;
		maxx = (maxx>res)?res:maxx;
		
		for (y=0; y < parts[1]; y++) {
			sizey = ceil(res / (float)parts[1]);
			miny = y * sizey;
			maxy = miny + sizey;
			maxy = (maxy>res)?res:maxy;
			
			for (z=0; z < parts[2]; z++) {
				VolPrecachePart *pa= MEM_callocN(sizeof(VolPrecachePart), "new precache part");
				
				sizez = ceil(res / (float)parts[2]);
				minz = z * sizez;
				maxz = minz + sizez;
				maxz = (maxz>res)?res: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);
				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;
}

/* 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, float *bbmin, float *bbmax)
{
	VolPrecachePart *nextpa, *pa;
	RayTree *tree;
	ShadeInput shi;
	ListBase threads;
	const int res = ma->vol_precache_resolution;
	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. */
	tree = create_raytree_obi(obi, bbmin, bbmax);
	if (!tree) return;
	
	obi->volume_precache = MEM_callocN(sizeof(float)*res*res*res*3, "volume light cache");

	/* Need a shadeinput to calculate scattering */
	precache_setup_shadeinput(re, obi, ma, &shi);
	
	precache_init_parts(re, tree, &shi, obi, bbmin, bbmax, res, 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())) {
			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() && 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->i);
			re->i.infostr= NULL;
			lasttime= time;
		}
	}
	
	BLI_end_threads(&threads);
	BLI_freelistN(&re->volume_precache_parts);
	
	if(tree) {
		RE_ray_tree_free(tree);
		tree= NULL;
	}
	
	lightcache_filter(obi->volume_precache, res);
	
	if (ELEM(ma->vol_shade_type, MA_VOL_SHADE_MULTIPLE, MA_VOL_SHADE_SINGLEPLUSMULTIPLE))
	{
		multiple_scattering_diffusion(re, obi->volume_precache, res, ma);
	}
}

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;

	/* ignore light cache and multiple scattering for preview render .. for now? */
	if (re->r.scemode & R_PREVIEWBUTS) return;

	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, obi->obr->boundbox[0], obi->obr->boundbox[1]);
				}
			}
		}
	}
	
	re->i.infostr= NULL;
	re->stats_draw(&re->i);
}

void free_volume_precache(Render *re)
{
	ObjectInstanceRen *obi;
	
	for(obi= re->instancetable.first; obi; obi= obi->next) {
		if (obi->volume_precache)
			MEM_freeN(obi->volume_precache);
	}
	
	BLI_freelistN(&re->volumes);
}

int point_inside_volume_objectinstance(ObjectInstanceRen *obi, float *co)
{
	RayTree *tree;
	int inside=0;
	
	tree = create_raytree_obi(obi, obi->obr->boundbox[0], obi->obr->boundbox[1]);
	if (!tree) return 0;
	
	inside = point_inside_obi(tree, obi, co);
	
	RE_ray_tree_free(tree);
	tree= NULL;
	
	return inside;
}