text: /*
   * ***** 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
   *
   * The Original Code is Copyright (C) 2008 Blender Foundation.
   * All rights reserved.
   *
   * The Original Code is: all of this file.
   *
   * Contributor(s): Brecht Van Lommel.
   *
   * ***** END GPL LICENSE BLOCK *****
   */
  
  /** \file blender/render/intern/source/occlusion.c
   *  \ingroup render
   */
  
  #include <math.h>
  #include <stdio.h>
  #include <stdlib.h>
  #include <string.h>
  
  #include "MEM_guardedalloc.h"
  
  #include "DNA_material_types.h"
  
  #include "BLI_math.h"
  #include "BLI_memarena.h"
  #include "BLI_threads.h"
  #include "BLI_utildefines.h"
  
  #include "BLT_translation.h"
  
  #include "BKE_node.h"
  #include "BKE_scene.h"
  
  
  #include "RE_shader_ext.h"
  
  /* local includes */
  #include "occlusion.h"
  #include "render_types.h"
  #include "rendercore.h"
  #include "renderdatabase.h"
  #include "shading.h"
  
  /* ------------------------- Declarations --------------------------- */
  
  #define INVPI ((float)M_1_PI)
  #define TOTCHILD 8
  #define CACHE_STEP 3
  
  typedef struct OcclusionCacheSample {
  	float co[3], n[3], ao[3], env[3], indirect[3], intensity, dist2;
  	int x, y, filled;
  } OcclusionCacheSample;
  
  typedef struct OcclusionCache {
  	OcclusionCacheSample *sample;
  	int x, y, w, h, step;
  } OcclusionCache;
  
  typedef struct OccFace {
  	int obi;
  	int facenr;
  } OccFace;
  
  typedef struct OccNode {
  	float co[3], area;
  	float sh[9], dco;
  	float occlusion, rad[3];
  	int childflag;
  	union {
  		//OccFace face;
  		int face;
  		struct OccNode *node;
  	} child[TOTCHILD];
  } OccNode;
  
  typedef struct OcclusionTree {
  	MemArena *arena;
  
  	float (*co)[3];     /* temporary during build */
  
  	OccFace *face;      /* instance and face indices */
  	float *occlusion;   /* occlusion for faces */
  	float (*rad)[3];    /* radiance for faces */
  	
  	OccNode *root;
  
  	OccNode **stack[BLENDER_MAX_THREADS];
  	int maxdepth;
  
  	int totface;
  
  	float error;
  	float distfac;
  
  	int dothreadedbuild;
  	int totbuildthread;
  	int doindirect;
  
  	OcclusionCache *cache;
  
  	int num_threads;
  } OcclusionTree;
  
  typedef struct OcclusionThread {
  	Render *re;
  	StrandSurface *mesh;
  	float (*faceao)[3];
  	float (*faceenv)[3];
  	float (*faceindirect)[3];
  	int begin, end;
  	int thread;
  } OcclusionThread;
  
  typedef struct OcclusionBuildThread {
  	OcclusionTree *tree;
  	int begin, end, depth;
  	OccNode *node;
  } OcclusionBuildThread;
  
  /* ------------------------- Shading --------------------------- */
  
  extern Render R; /* meh */
  
  static void occ_shade(ShadeSample *ssamp, ObjectInstanceRen *obi, VlakRen *vlr, float *rad)
  {
  	ShadeInput *shi = ssamp->shi;
  	ShadeResult *shr = ssamp->shr;
  	float l, u, v, *v1, *v2, *v3;
  	
  	/* init */
  	if (vlr->v4) {
  		shi->u = u = 0.5f;
  		shi->v = v = 0.5f;
  	}
  	else {
  		shi->u = u = 1.0f / 3.0f;
  		shi->v = v = 1.0f / 3.0f;
  	}
  
  	/* setup render coordinates */
  	v1 = vlr->v1->co;
  	v2 = vlr->v2->co;
  	v3 = vlr->v3->co;
  	
  	/* renderco */
  	l = 1.0f - u - v;
  	
  	shi->co[0] = l * v3[0] + u * v1[0] + v * v2[0];
  	shi->co[1] = l * v3[1] + u * v1[1] + v * v2[1];
  	shi->co[2] = l * v3[2] + u * v1[2] + v * v2[2];
  
  	shade_input_set_triangle_i(shi, obi, vlr, 0, 1, 2);
  
  	/* set up view vector */
  	copy_v3_v3(shi->view, shi->co);
  	normalize_v3(shi->view);
  	
  	/* cache for shadow */
  	shi->samplenr++;
  
  	shi->xs = 0; /* TODO */
  	shi->ys = 0;
  	
  	shade_input_set_normals(shi);
  
  	/* no normal flip */
  	if (shi->flippednor)
  		shade_input_flip_normals(shi);
  
  	madd_v3_v3fl(shi->co, shi->facenor, -0.0001f); /* ugly.. */
  
  	/* not a pretty solution, but fixes common cases */
  	if (shi->obr->ob && shi->obr->ob->transflag & OB_NEG_SCALE) {
  		negate_v3(shi->vn);
  		negate_v3(shi->vno);
  		negate_v3(shi->nmapnorm);
  	}
  
  	/* init material vars */
  	shade_input_init_material(shi);
  
  	/* render */
  	shade_input_set_shade_texco(shi);
  
  	if (shi->mat->nodetree && shi->mat->use_nodes) {
  		ntreeShaderExecTree(shi->mat->nodetree, shi, shr);
  		shi->mat = vlr->mat;  /* shi->mat is being set in nodetree */
  	}
  	else {
  		shade_material_loop(shi, shr);
  	}
  
  	copy_v3_v3(rad, shr->combined);
  }
  
  static void occ_build_shade(Render *re, OcclusionTree *tree)
  {
  	ShadeSample ssamp;
  	ObjectInstanceRen *obi;
  	VlakRen *vlr;
  	int a;
  
  	R = *re;
  
  	/* setup shade sample with correct passes */
  	memset(&ssamp, 0, sizeof(ShadeSample));
  	ssamp.shi[0].lay = re->lay;
  	ssamp.shi[0].passflag = SCE_PASS_DIFFUSE | SCE_PASS_RGBA;
  	ssamp.shi[0].combinedflag = ~(SCE_PASS_SPEC);
  	ssamp.tot = 1;
  
  	for (a = 0; a < tree->totface; a++) {
  		obi = &R.objectinstance[tree->face[a].obi];
  		vlr = RE_findOrAddVlak(obi->obr, tree->face[a].facenr);
  
  		occ_shade(&ssamp, obi, vlr, tree->rad[a]);
  
  		if (re->test_break(re->tbh))
  			break;
  	}
  }
  
  /* ------------------------- Spherical Harmonics --------------------------- */
  
  /* Use 2nd order SH => 9 coefficients, stored in this order:
   * 0 = (0,0),
   * 1 = (1,-1), 2 = (1,0), 3 = (1,1),
   * 4 = (2,-2), 5 = (2,-1), 6 = (2,0), 7 = (2,1), 8 = (2,2) */
  
  static void sh_copy(float *shresult, float *sh)
  {
  	memcpy(shresult, sh, sizeof(float) * 9);
  }
  
  static void sh_mul(float *sh, float f)
  {
  	int i;
  
  	for (i = 0; i < 9; i++)
  		sh[i] *= f;
  }
  
  static void sh_add(float *shresult, float *sh1, float *sh2)
  {
  	int i;
  
  	for (i = 0; i < 9; i++)
  		shresult[i] = sh1[i] + sh2[i];
  }
  
  static void sh_from_disc(float *n, float area, float *shresult)
  {
  	/* See formula (3) in:
  	 * "An Efficient Representation for Irradiance Environment Maps" */
  	float sh[9], x, y, z;
  
  	x = n[0];
  	y = n[1];
  	z = n[2];
  
  	sh[0] = 0.282095f;
  
  	sh[1] = 0.488603f * y;
  	sh[2] = 0.488603f * z;
  	sh[3] = 0.488603f * x;
  
  	sh[4] = 1.092548f * x * y;
  	sh[5] = 1.092548f * y * z;
  	sh[6] = 0.315392f * (3.0f * z * z - 1.0f);
  	sh[7] = 1.092548f * x * z;
  	sh[8] = 0.546274f * (x * x - y * y);
  
  	sh_mul(sh, area);
  	sh_copy(shresult, sh);
  }
  
  static float sh_eval(float *sh, float *v)
  {
  	/* See formula (13) in:
  	 * "An Efficient Representation for Irradiance Environment Maps" */
  	static const float c1 = 0.429043f, c2 = 0.511664f, c3 = 0.743125f;
  	static const float c4 = 0.886227f, c5 = 0.247708f;
  	float x, y, z, sum;
  
  	x = v[0];
  	y = v[1];
  	z = v[2];
  
  	sum = c1 * sh[8] * (x * x - y * y);
  	sum += c3 * sh[6] * z * z;
  	sum += c4 * sh[0];
  	sum += -c5 * sh[6];
  	sum += 2.0f * c1 * (sh[4] * x * y + sh[7] * x * z + sh[5] * y * z);
  	sum += 2.0f * c2 * (sh[3] * x + sh[1] * y + sh[2] * z);
  
  	return sum;
  }
  
  /* ------------------------------ Building --------------------------------- */
  
  static void occ_face(const OccFace *face, float co[3], float normal[3], float *area)
  {
  	ObjectInstanceRen *obi;
  	VlakRen *vlr;
  	float v1[3], v2[3], v3[3], v4[3];
  
  	obi = &R.objectinstance[face->obi];
  	vlr = RE_findOrAddVlak(obi->obr, face->facenr);
  	
  	if (co) {
  		if (vlr->v4)
  			mid_v3_v3v3(co, vlr->v1->co, vlr->v3->co);
  		else
  			mid_v3_v3v3v3(co, vlr->v1->co, vlr->v2->co, vlr->v3->co);
  
  		if (obi->flag & R_TRANSFORMED)
  			mul_m4_v3(obi->mat, co);
  	}
  	
  	if (normal) {
  		normal[0] = -vlr->n[0];
  		normal[1] = -vlr->n[1];
  		normal[2] = -vlr->n[2];
  
  		if (obi->flag & R_TRANSFORMED)
  			mul_m3_v3(obi->nmat, normal);
  	}
  
  	if (area) {
  		copy_v3_v3(v1, vlr->v1->co);
  		copy_v3_v3(v2, vlr->v2->co);
  		copy_v3_v3(v3, vlr->v3->co);
  		if (vlr->v4) copy_v3_v3(v4, vlr->v4->co);
  
  		if (obi->flag & R_TRANSFORMED) {
  			mul_m4_v3(obi->mat, v1);
  			mul_m4_v3(obi->mat, v2);
  			mul_m4_v3(obi->mat, v3);
  			if (vlr->v4) mul_m4_v3(obi->mat, v4);
  		}
  
  		/* todo: correct area for instances */
  		if (vlr->v4)
  			*area = area_quad_v3(v1, v2, v3, v4);
  		else
  			*area = area_tri_v3(v1, v2, v3);
  	}
  }
  
  static void occ_sum_occlusion(OcclusionTree *tree, OccNode *node)
  {
  	OccNode *child;
  	float occ, area, totarea, rad[3];
  	int a, b, indirect = tree->doindirect;
  
  	occ = 0.0f;
  	totarea = 0.0f;
  	if (indirect) zero_v3(rad);
  
  	for (b = 0; b < TOTCHILD; b++) {
  		if (node->childflag & (1 << b)) {
  			a = node->child[b].face;
  			occ_face(&tree->face[a], NULL, NULL, &area);
  			occ += area * tree->occlusion[a];
  			if (indirect) madd_v3_v3fl(rad, tree->rad[a], area);
  			totarea += area;
  		}
  		else if (node->child[b].node) {
  			child = node->child[b].node;
  			occ_sum_occlusion(tree, child);
  
  			occ += child->area * child->occlusion;
  			if (indirect) madd_v3_v3fl(rad, child->rad, child->area);
  			totarea += child->area;
  		}
  	}
  
  	if (totarea != 0.0f) {
  		occ /= totarea;
  		if (indirect) mul_v3_fl(rad, 1.0f / totarea);
  	}
  	
  	node->occlusion = occ;
  	if (indirect) copy_v3_v3(node->rad, rad);
  }
  
  static int occ_find_bbox_axis(OcclusionTree *tree, int begin, int end, float *min, float *max)
  {
  	float len, maxlen = -1.0f;
  	int a, axis = 0;
  
  	INIT_MINMAX(min, max);
  
  	for (a = begin; a < end; a++) {
  		minmax_v3v3_v3(min, max, tree->co[a]);
  	}
  
  	for (a = 0; a < 3; a++) {
  		len = max[a] - min[a];
  
  		if (len > maxlen) {
  			maxlen = len;
  			axis = a;
  		}
  	}
  
  	return axis;
  }
  
  static void occ_node_from_face(OccFace *face, OccNode *node)
  {
  	float n[3];
  
  	occ_face(face, node->co, n, &node->area);
  	node->dco = 0.0f;
  	sh_from_disc(n, node->area, node->sh);
  }
  
  static void occ_build_dco(OcclusionTree *tree, OccNode *node, const float co[3], float *dco)
  {
  	int b;
  	for (b = 0; b < TOTCHILD; b++) {
  		float dist, d[3], nco[3];
  
  		if (node->childflag & (1 << b)) {
  			occ_face(tree->face + node->child[b].face, nco, NULL, NULL);
  		}
  		else if (node->child[b].node) {
  			OccNode *child = node->child[b].node;
  			occ_build_dco(tree, child, co, dco);
  			copy_v3_v3(nco, child->co);
  		}
  		else {
  			continue;
  		}
  
  		sub_v3_v3v3(d, nco, co);
  		dist = dot_v3v3(d, d);
  		if (dist > *dco)
  			*dco = dist;
  	}
  }
  
  static void occ_build_split(OcclusionTree *tree, int begin, int end, int *split)
  {
  	float min[3], max[3], mid;
  	int axis, a, enda;
  
  	/* split in middle of boundbox. this seems faster than median split
  	 * on complex scenes, possibly since it avoids two distant faces to
  	 * be in the same node better? */
  	axis = occ_find_bbox_axis(tree, begin, end, min, max);
  	mid = 0.5f * (min[axis] + max[axis]);
  
  	a = begin;
  	enda = end;
  	while (a < enda) {
  		if (tree->co[a][axis] > mid) {
  			enda--;
  			SWAP(OccFace, tree->face[a], tree->face[enda]);
  			swap_v3_v3(tree->co[a], tree->co[enda]);
  		}
  		else
  			a++;
  	}
  
  	*split = enda;
  }
  
  static void occ_build_8_split(OcclusionTree *tree, int begin, int end, int *offset, int *count)
  {
  	/* split faces into eight groups */
  	int b, splitx, splity[2], splitz[4];
  
  	occ_build_split(tree, begin, end, &splitx);
  
  	/* force split if none found, to deal with degenerate geometry */
  	if (splitx == begin || splitx == end)
  		splitx = (begin + end) / 2;
  
  	occ_build_split(tree, begin, splitx, &splity[0]);
  	occ_build_split(tree, splitx, end, &splity[1]);
  
  	occ_build_split(tree, begin, splity[0], &splitz[0]);
  	occ_build_split(tree, splity[0], splitx, &splitz[1]);
  	occ_build_split(tree, splitx, splity[1], &splitz[2]);
  	occ_build_split(tree, splity[1], end, &splitz[3]);
  
  	offset[0] = begin;
  	offset[1] = splitz[0];
  	offset[2] = splity[0];
  	offset[3] = splitz[1];
  	offset[4] = splitx;
  	offset[5] = splitz[2];
  	offset[6] = splity[1];
  	offset[7] = splitz[3];
  
  	for (b = 0; b < 7; b++)
  		count[b] = offset[b + 1] - offset[b];
  	count[7] = end - offset[7];
  }
  
  static void occ_build_recursive(OcclusionTree *tree, OccNode *node, int begin, int end, int depth);
  
  static void *exec_occ_build(void *data)
  {
  	OcclusionBuildThread *othread = (OcclusionBuildThread *)data;
  
  	occ_build_recursive(othread->tree, othread->node, othread->begin, othread->end, othread->depth);
  
  	return NULL;
  }
  
  static void occ_build_recursive(OcclusionTree *tree, OccNode *node, int begin, int end, int depth)
  {
  	ListBase threads;
  	OcclusionBuildThread othreads[BLENDER_MAX_THREADS];
  	OccNode *child, tmpnode;
  	/* OccFace *face; */
  	int a, b, totthread = 0, offset[TOTCHILD], count[TOTCHILD];
  
  	/* add a new node */
  	node->occlusion = 1.0f;
  
  	/* leaf node with only children */
  	if (end - begin <= TOTCHILD) {
  		for (a = begin, b = 0; a < end; a++, b++) {
  			/* face= &tree->face[a]; */
  			node->child[b].face = a;
  			node->childflag |= (1 << b);
  		}
  	}
  	else {
  		/* order faces */
  		occ_build_8_split(tree, begin, end, offset, count);
  
  		if (depth == 1 && tree->dothreadedbuild)
  			BLI_init_threads(&threads, exec_occ_build, tree->totbuildthread);
  
  		for (b = 0; b < TOTCHILD; b++) {
  			if (count[b] == 0) {
  				node->child[b].node = NULL;
  			}
  			else if (count[b] == 1) {
  				/* face= &tree->face[offset[b]]; */
  				node->child[b].face = offset[b];
  				node->childflag |= (1 << b);
  			}
  			else {
  				if (tree->dothreadedbuild)
  					BLI_lock_thread(LOCK_CUSTOM1);
  
  				child = BLI_memarena_alloc(tree->arena, sizeof(OccNode));
  				node->child[b].node = child;
  
  				/* keep track of maximum depth for stack */
  				if (depth >= tree->maxdepth)
  					tree->maxdepth = depth + 1;
  
  				if (tree->dothreadedbuild)
  					BLI_unlock_thread(LOCK_CUSTOM1);
  
  				if (depth == 1 && tree->dothreadedbuild) {
  					othreads[totthread].tree = tree;
  					othreads[totthread].node = child;
  					othreads[totthread].begin = offset[b];
  					othreads[totthread].end = offset[b] + count[b];
  					othreads[totthread].depth = depth + 1;
  					BLI_insert_thread(&threads, &othreads[totthread]);
  					totthread++;
  				}
  				else
  					occ_build_recursive(tree, child, offset[b], offset[b] + count[b], depth + 1);
  			}
  		}
  
  		if (depth == 1 && tree->dothreadedbuild)
  			BLI_end_threads(&threads);
  	}
  
  	/* combine area, position and sh */
  	for (b = 0; b < TOTCHILD; b++) {
  		if (node->childflag & (1 << b)) {
  			child = &tmpnode;
  			occ_node_from_face(tree->face + node->child[b].face, &tmpnode);
  		}
  		else {
  			child = node->child[b].node;
  		}
  
  		if (child) {
  			node->area += child->area;
  			sh_add(node->sh, node->sh, child->sh);
  			madd_v3_v3fl(node->co, child->co, child->area);
  		}
  	}
  
  	if (node->area != 0.0f)
  		mul_v3_fl(node->co, 1.0f / node->area);
  
  	/* compute maximum distance from center */
  	node->dco = 0.0f;
  	if (node->area > 0.0f)
  		occ_build_dco(tree, node, node->co, &node->dco);
  }
  
  static void occ_build_sh_normalize(OccNode *node)
  {
  	/* normalize spherical harmonics to not include area, so
  	 * we can clamp the dot product and then multiply by area */
  	int b;
  
  	if (node->area != 0.0f)
  		sh_mul(node->sh, 1.0f / node->area);
  
  	for (b = 0; b < TOTCHILD; b++) {
  		if (node->childflag & (1 << b)) {
  			/* pass */
  		}
  		else if (node->child[b].node) {
  			occ_build_sh_normalize(node->child[b].node);
  		}
  	}
  }
  
  static OcclusionTree *occ_tree_build(Render *re)
  {
  	const int num_threads = re->r.threads;
  	OcclusionTree *tree;
  	ObjectInstanceRen *obi;
  	ObjectRen *obr;
  	Material *ma;
  	VlakRen *vlr = NULL;
  	int a, b, c, totface;
  
  	/* count */
  	totface = 0;
  	for (obi = re->instancetable.first; obi; obi = obi->next) {
  		obr = obi->obr;
  		for (a = 0; a < obr->totvlak; a++) {
  			if ((a & 255) == 0) vlr = obr->vlaknodes[a >> 8].vlak;
  			else vlr++;
  
  			ma = vlr->mat;
  
  			if ((ma->shade_flag & MA_APPROX_OCCLUSION) && (ma->material_type == MA_TYPE_SURFACE))
  				totface++;
  		}
  	}
  
  	if (totface == 0)
  		return NULL;
  	
  	tree = MEM_callocN(sizeof(OcclusionTree), "OcclusionTree");
  	tree->totface = totface;
  
  	/* parameters */
  	tree->error = get_render_aosss_error(&re->r, re->wrld.ao_approx_error);
  	tree->distfac = (re->wrld.aomode & WO_AODIST) ? re->wrld.aodistfac : 0.0f;
  	tree->doindirect = (re->wrld.ao_indirect_energy > 0.0f && re->wrld.ao_indirect_bounces > 0);
  
  	/* allocation */
  	tree->arena = BLI_memarena_new(0x8000 * sizeof(OccNode), "occ tree arena");
  	BLI_memarena_use_calloc(tree->arena);
  
  	if (re->wrld.aomode & WO_AOCACHE)
  		tree->cache = MEM_callocN(sizeof(OcclusionCache) * num_threads, "OcclusionCache");
  
  	tree->face = MEM_callocN(sizeof(OccFace) * totface, "OcclusionFace");
  	tree->co = MEM_callocN(sizeof(float) * 3 * totface, "OcclusionCo");
  	tree->occlusion = MEM_callocN(sizeof(float) * totface, "OcclusionOcclusion");
  
  	if (tree->doindirect)
  		tree->rad = MEM_callocN(sizeof(float) * 3 * totface, "OcclusionRad");
  
  	/* make array of face pointers */
  	for (b = 0, c = 0, obi = re->instancetable.first; obi; obi = obi->next, c++) {
  		obr = obi->obr;
  		for (a = 0; a < obr->totvlak; a++) {
  			if ((a & 255) == 0) vlr = obr->vlaknodes[a >> 8].vlak;
  			else vlr++;
  
  			ma = vlr->mat;
  
  			if ((ma->shade_flag & MA_APPROX_OCCLUSION) && (ma->material_type == MA_TYPE_SURFACE)) {
  				tree->face[b].obi = c;
  				tree->face[b].facenr = a;
  				tree->occlusion[b] = 1.0f;
  				occ_face(&tree->face[b], tree->co[b], NULL, NULL); 
  				b++;
  			}
  		}
  	}
  
  	/* threads */
  	tree->totbuildthread = (re->r.threads > 1 && totface > 10000) ? 8 : 1;
  	tree->dothreadedbuild = (tree->totbuildthread > 1);
  
  	/* recurse */
  	tree->root = BLI_memarena_alloc(tree->arena, sizeof(OccNode));
  	tree->maxdepth = 1;
  	occ_build_recursive(tree, tree->root, 0, totface, 1);
  
  	if (tree->doindirect) {
  		if (!(re->test_break(re->tbh)))
  			occ_build_shade(re, tree);
  
  		if (!(re->test_break(re->tbh)))
  			occ_sum_occlusion(tree, tree->root);
  	}
  	
  	MEM_freeN(tree->co);
  	tree->co = NULL;
  
  	if (!(re->test_break(re->tbh)))
  		occ_build_sh_normalize(tree->root);
  
  	for (a = 0; a < num_threads; a++)
  		tree->stack[a] = MEM_callocN(sizeof(OccNode) * TOTCHILD * (tree->maxdepth + 1), "OccStack");
  
  	tree->num_threads = num_threads;
  
  	return tree;
  }
  
  static void occ_free_tree(OcclusionTree *tree)
  {
  	int a;
  
  	if (tree) {
  		if (tree->arena) BLI_memarena_free(tree->arena);
  		for (a = 0; a < tree->num_threads; a++)
  			if (tree->stack[a])
  				MEM_freeN(tree->stack[a]);
  		if (tree->occlusion) MEM_freeN(tree->occlusion);
  		if (tree->cache) MEM_freeN(tree->cache);
  		if (tree->face) MEM_freeN(tree->face);
  		if (tree->rad) MEM_freeN(tree->rad);
  		MEM_freeN(tree);
  	}
  }
  
  /* ------------------------- Traversal --------------------------- */
  
  static float occ_solid_angle(OccNode *node, const float v[3], float d2, float invd2, const float receivenormal[3])
  {
  	float dotreceive, dotemit;
  	float ev[3];
  
  	ev[0] = -v[0] * invd2;
  	ev[1] = -v[1] * invd2;
  	ev[2] = -v[2] * invd2;
  	dotemit = sh_eval(node->sh, ev);
  	dotreceive = dot_v3v3(receivenormal, v) * invd2;
  
  	CLAMP(dotemit, 0.0f, 1.0f);
  	CLAMP(dotreceive, 0.0f, 1.0f);
  	
  	return ((node->area * dotemit * dotreceive) / (d2 + node->area * INVPI)) * INVPI;
  }
  
  static float occ_form_factor(OccFace *face, float *p, float *n)
  {
  	ObjectInstanceRen *obi;
  	VlakRen *vlr;
  	float v1[3], v2[3], v3[3], v4[3], q0[3], q1[3], q2[3], q3[3], contrib = 0.0f;
  
  	obi = &R.objectinstance[face->obi];
  	vlr = RE_findOrAddVlak(obi->obr, face->facenr);
  
  	copy_v3_v3(v1, vlr->v1->co);
  	copy_v3_v3(v2, vlr->v2->co);
  	copy_v3_v3(v3, vlr->v3->co);
  
  	if (obi->flag & R_TRANSFORMED) {
  		mul_m4_v3(obi->mat, v1);
  		mul_m4_v3(obi->mat, v2);
  		mul_m4_v3(obi->mat, v3);
  	}
  
  	if (form_factor_visible_quad(p, n, v1, v2, v3, q0, q1, q2, q3))
  		contrib += form_factor_quad(p, n, q0, q1, q2, q3);
  
  	if (vlr->v4) {
  		copy_v3_v3(v4, vlr->v4->co);
  		if (obi->flag & R_TRANSFORMED)
  			mul_m4_v3(obi->mat, v4);
  
  		if (form_factor_visible_quad(p, n, v1, v3, v4, q0, q1, q2, q3))
  			contrib += form_factor_quad(p, n, q0, q1, q2, q3);
  	}
  
  	return contrib;
  }
  
  static void occ_lookup(OcclusionTree *tree, int thread, OccFace *exclude,
                         const float pp[3], const float pn[3], float *occ, float rad[3], float bentn[3])
  {
  	OccNode *node, **stack;
  	OccFace *face;
  	float resultocc, resultrad[3], v[3], p[3], n[3], co[3], invd2;
  	float distfac, fac, error, d2, weight, emitarea;
  	int b, f, totstack;
  
  	/* init variables */
  	copy_v3_v3(p, pp);
  	copy_v3_v3(n, pn);
  	madd_v3_v3fl(p, n, 1e-4f);
  
  	if (bentn)
  		copy_v3_v3(bentn, n);
  	
  	error = tree->error;
  	distfac = tree->distfac;
  
  	resultocc = 0.0f;
  	zero_v3(resultrad);
  
  	/* init stack */
  	stack = tree->stack[thread];
  	stack[0] = tree->root;
  	totstack = 1;
  
  	while (totstack) {
  		/* pop point off the stack */
  		node = stack[--totstack];
  
  		sub_v3_v3v3(v, node->co, p);
  		d2 = dot_v3v3(v, v) + 1e-16f;
  		emitarea = MAX2(node->area, node->dco);
  
  		if (d2 * error > emitarea) {
  			if (distfac != 0.0f) {
  				fac = 1.0f / (1.0f + distfac * d2);
  				if (fac < 0.01f)
  					continue;
  			}
  			else
  				fac = 1.0f;
  
  			/* accumulate occlusion from spherical harmonics */
  			invd2 = 1.0f / sqrtf(d2);
  			weight = occ_solid_angle(node, v, d2, invd2, n);
  
  			if (rad)
  				madd_v3_v3fl(resultrad, node->rad, weight * fac);
  
  			weight *= node->occlusion;
  
  			if (bentn) {
  				bentn[0] -= weight * invd2 * v[0];
  				bentn[1] -= weight * invd2 * v[1];
  				bentn[2] -= weight * invd2 * v[2];
  			}
  
  			resultocc += weight * fac;
  		}
  		else {
  			/* traverse into children */
  			for (b = 0; b < TOTCHILD; b++) {
  				if (node->childflag & (1 << b)) {
  					f = node->child[b].face;
  					face = &tree->face[f];
  
  					/* accumulate occlusion with face form factor */
  					if (!exclude || !(face->obi == exclude->obi && face->facenr == exclude->facenr)) {
  						if (bentn || distfac != 0.0f) {
  							occ_face(face, co, NULL, NULL); 
  							sub_v3_v3v3(v, co, p);
  							d2 = dot_v3v3(v, v) + 1e-16f;
  
  							fac = (distfac == 0.0f) ? 1.0f : 1.0f / (1.0f + distfac * d2);
  							if (fac < 0.01f)
  								continue;
  						}
  						else
  							fac = 1.0f;
  
  						weight = occ_form_factor(face, p, n);
  
  						if (rad)
  							madd_v3_v3fl(resultrad, tree->rad[f], weight * fac);
  
  						weight *= tree->occlusion[f];
  
  						if (bentn) {
  							invd2 = 1.0f / sqrtf(d2);
  							bentn[0] -= weight * invd2 * v[0];
  							bentn[1] -= weight * invd2 * v[1];
  							bentn[2] -= weight * invd2 * v[2];
  						}
  
  						resultocc += weight * fac;
  					}
  				}
  				else if (node->child[b].node) {
  					/* push child on the stack */
  					stack[totstack++] = node->child[b].node;
  				}
  			}
  		}
  	}
  
  	if (occ) *occ = resultocc;
  	if (rad) copy_v3_v3(rad, resultrad);
  #if 0
  	if (rad && exclude) {
  		int a;
  		for (a = 0; a < tree->totface; a++)
  			if ((tree->face[a].obi == exclude->obi && tree->face[a].facenr == exclude->facenr))
  				copy_v3_v3(rad, tree->rad[a]);
  	}
  #endif
  	if (bentn) normalize_v3(bentn);
  }
  
  static void occ_compute_bounces(Render *re, OcclusionTree *tree, int totbounce)
  {
  	float (*rad)[3], (*sum)[3], (*tmp)[3], co[3], n[3], occ;
  	int bounce, i;
  
  	rad = MEM_callocN(sizeof(float) * 3 * tree->totface, "OcclusionBounceRad");
  	sum = MEM_dupallocN(tree->rad);
  
  	for (bounce = 1; bounce < totbounce; bounce++) {
  		for (i = 0; i < tree->totface; i++) {
  			occ_face(&tree->face[i], co, n, NULL);
  			madd_v3_v3fl(co, n, 1e-8f);
  
  			occ_lookup(tree, 0, &tree->face[i], co, n, &occ, rad[i], NULL);
  			rad[i][0] = MAX2(rad[i][0], 0.0f);
  			rad[i][1] = MAX2(rad[i][1], 0.0f);
  			rad[i][2] = MAX2(rad[i][2], 0.0f);
  			add_v3_v3(sum[i], rad[i]);
  
  			if (re->test_break(re->tbh))
  				break;
  		}
  
  		if (re->test_break(re->tbh))
  			break;
  
  		tmp = tree->rad;
  		tree->rad = rad;
  		rad = tmp;
  
  		occ_sum_occlusion(tree, tree->root);
  	}
  
  	MEM_freeN(rad);
  	MEM_freeN(tree->rad);
  	tree->rad = sum;
  
  	if (!re->test_break(re->tbh))
  		occ_sum_occlusion(tree, tree->root);
  }
  
  static void occ_compute_passes(Render *re, OcclusionTree *tree, int totpass)
  {
  	float *occ, co[3], n[3];
  	int pass, i;
  	
  	occ = MEM_callocN(sizeof(float) * tree->totface, "OcclusionPassOcc");
  
  	for (pass = 0; pass < totpass; pass++) {
  		for (i = 0; i < tree->totface; i++) {
  			occ_face(&tree->face[i], co, n, NULL);
  			negate_v3(n);
  			madd_v3_v3fl(co, n, 1e-8f);
  
  			occ_lookup(tree, 0, &tree->face[i], co, n, &occ[i], NULL, NULL);
  			if (re->test_break(re->tbh))
  				break;
  		}
  
  		if (re->test_break(re->tbh))
  			break;
  
  		for (i = 0; i < tree->totface; i++) {
  			tree->occlusion[i] -= occ[i]; //MAX2(1.0f-occ[i], 0.0f);
  			if (tree->occlusion[i] < 0.0f)
  				tree->occlusion[i] = 0.0f;
  		}
  
  		occ_sum_occlusion(tree, tree->root);
  	}
  
  	MEM_freeN(occ);
  }
  
  static void sample_occ_tree(Render *re, OcclusionTree *tree, OccFace *exclude,
                              const float co[3], const float n[3], int thread, int onlyshadow,
                              float *ao, float *env, float *indirect)
  {
  	float nn[3], bn[3], fac, occ, occlusion, correction, rad[3];
  	int envcolor;
  
  	envcolor = re->wrld.aocolor;
  	if (onlyshadow)
  		envcolor = WO_AOPLAIN;
  
  	negate_v3_v3(nn, n);
  
  	occ_lookup(tree, thread, exclude, co, nn, &occ, (tree->doindirect) ? rad : NULL, (env && envcolor) ? bn : NULL);
  
  	correction = re->wrld.ao_approx_correction;
  
  	occlusion = (1.0f - correction) * (1.0f - occ);
  	CLAMP(occlusion, 0.0f, 1.0f);
  	if (correction != 0.0f)
  		occlusion += correction * expf(-occ);
  
  	if (env) {
  		/* sky shading using bent normal */
  		if (ELEM(envcolor, WO_AOSKYCOL, WO_AOSKYTEX)) {
  			fac = 0.5f * (1.0f + dot_v3v3(bn, re->grvec));
  			env[0] = (1.0f - fac) * re->wrld.horr + fac * re->wrld.zenr;
  			env[1] = (1.0f - fac) * re->wrld.horg + fac * re->wrld.zeng;
  			env[2] = (1.0f - fac) * re->wrld.horb + fac * re->wrld.zenb;
  
  			mul_v3_fl(env, occlusion);
  		}
  		else {
  			env[0] = occlusion;
  			env[1] = occlusion;
  			env[2] = occlusion;
  		}
  #if 0
  		else {  /* WO_AOSKYTEX */
  			float dxyview[3];
  			bn[0] = -bn[0];
  			bn[1] = -bn[1];
  			bn[2] = -bn[2];
  			dxyview[0] = 1.0f;
  			dxyview[1] = 1.0f;
  			dxyview[2] = 0.0f;
  			shadeSkyView(ao, co, bn, dxyview);
  		}
  #endif
  	}
  
  	if (ao) {
  		ao[0] = occlusion;
  		ao[1] = occlusion;
  		ao[2] = occlusion;
  	}
  
  	if (tree->doindirect) copy_v3_v3(indirect, rad);
  	else zero_v3(indirect);
  }
  
  /* ---------------------------- Caching ------------------------------- */
  
  static OcclusionCacheSample *find_occ_sample(OcclusionCache *cache, int x, int y)
  {
  	x -= cache->x;
  	y -= cache->y;
  
  	x /= cache->step;
  	y /= cache->step;
  	x *= cache->step;
  	y *= cache->step;
  
  	if (x < 0 || x >= cache->w || y < 0 || y >= cache->h)
  		return NULL;
  	else
  		return &cache->sample[y * cache->w + x];
  }
  
  static int sample_occ_cache(OcclusionTree *tree, float *co, float *n, int x, int y, int thread, float *ao, float *env, float *indirect)
  {
  	OcclusionCache *cache;
  	OcclusionCacheSample *samples[4], *sample;
  	float wn[4], wz[4], wb[4], tx, ty, w, totw, mino, maxo;
  	float d[3], dist2;
  	int i, x1, y1, x2, y2;
  
  	if (!tree->cache)
  		return 0;
  	
  	/* first try to find a sample in the same pixel */
  	cache = &tree->cache[thread];
  
  	if (cache->sample && cache->step) {
  		sample = &cache->sample[(y - cache->y) * cache->w + (x - cache->x)];
  		if (sample->filled) {
  			sub_v3_v3v3(d, sample->co, co);
  			dist2 = dot_v3v3(d, d);
  			if (dist2 < 0.5f * sample->dist2 && dot_v3v3(sample->n, n) > 0.98f) {
  				copy_v3_v3(ao, sample->ao);
  				copy_v3_v3(env, sample->env);
  				copy_v3_v3(indirect, sample->indirect);
  				return 1;
  			}
  		}
  	}
  	else
  		return 0;
  
  	/* try to interpolate between 4 neighboring pixels */
  	samples[0] = find_occ_sample(cache, x, y);
  	samples[1] = find_occ_sample(cache, x + cache->step, y);
  	samples[2] = find_occ_sample(cache, x, y + cache->step);
  	samples[3] = find_occ_sample(cache, x + cache->step, y + cache->step);
  
  	for (i = 0; i < 4; i++)
  		if (!samples[i] || !samples[i]->filled)
  			return 0;
  
  	/* require intensities not being too different */
  	mino = min_ffff(samples[0]->intensity, samples[1]->intensity, samples[2]->intensity, samples[3]->intensity);
  	maxo = max_ffff(samples[0]->intensity, samples[1]->intensity, samples[2]->intensity, samples[3]->intensity);
  
  	if (maxo - mino > 0.05f)
  		return 0;
  
  	/* compute weighted interpolation between samples */
  	zero_v3(ao);
  	zero_v3(env);
  	zero_v3(indirect);
  	totw = 0.0f;
  
  	x1 = samples[0]->x;
  	y1 = samples[0]->y;
  	x2 = samples[3]->x;
  	y2 = samples[3]->y;
  
  	tx = (float)(x2 - x) / (float)(x2 - x1);
  	ty = (float)(y2 - y) / (float)(y2 - y1);
  
  	wb[3] = (1.0f - tx) * (1.0f - ty);
  	wb[2] = (tx) * (1.0f - ty);
  	wb[1] = (1.0f - tx) * (ty);
  	wb[0] = tx * ty;
  
  	for (i = 0; i < 4; i++) {
  		sub_v3_v3v3(d, samples[i]->co, co);
  		//dist2 = dot_v3v3(d, d);
  
  		wz[i] = 1.0f; //(samples[i]->dist2/(1e-4f + dist2));
  		wn[i] = pow(dot_v3v3(samples[i]->n, n), 32.0f);
  
  		w = wb[i] * wn[i] * wz[i];
  
  		totw += w;
  		madd_v3_v3fl(ao, samples[i]->ao, w);
  		madd_v3_v3fl(env, samples[i]->env, w);
  		madd_v3_v3fl(indirect, samples[i]->indirect, w);
  	}
  
  	if (totw >= 0.9f) {
  		totw = 1.0f / totw;
  		mul_v3_fl(ao, totw);
  		mul_v3_fl(env, totw);
  		mul_v3_fl(indirect, totw);
  		return 1;
  	}
  
  	return 0;
  }
  
  static void sample_occ_surface(ShadeInput *shi)
  {
  	StrandRen *strand = shi->strand;
  	StrandSurface *mesh = strand->buffer->surface;
  	const int *face, *index = RE_strandren_get_face(shi->obr, strand, 0);
  	float w[4], *co1, *co2, *co3, *co4;
  
  	if (mesh && mesh->face && mesh->co && mesh->ao && index) {
  		face = mesh->face[*index];
  
  		co1 = mesh->co[face[0]];
  		co2 = mesh->co[face[1]];
  		co3 = mesh->co[face[2]];
  
  		if (face[3]) {
  			co4 = mesh->co[face[3]];
  			interp_weights_quad_v3(w, co1, co2, co3, co4, strand->vert->co);
  		}
  		else {
  			interp_weights_tri_v3(w, co1, co2, co3, strand->vert->co);
  		}
  
  		zero_v3(shi->ao);
  		zero_v3(shi->env);
  		zero_v3(shi->indirect);
  
  		madd_v3_v3fl(shi->ao, mesh->ao[face[0]], w[0]);
  		madd_v3_v3fl(shi->env, mesh->env[face[0]], w[0]);
  		madd_v3_v3fl(shi->indirect, mesh->indirect[face[0]], w[0]);
  		madd_v3_v3fl(shi->ao, mesh->ao[face[1]], w[1]);
  		madd_v3_v3fl(shi->env, mesh->env[face[1]], w[1]);
  		madd_v3_v3fl(shi->indirect, mesh->indirect[face[1]], w[1]);
  		madd_v3_v3fl(shi->ao, mesh->ao[face[2]], w[2]);
  		madd_v3_v3fl(shi->env, mesh->env[face[2]], w[2]);
  		madd_v3_v3fl(shi->indirect, mesh->indirect[face[2]], w[2]);
  		if (face[3]) {
  			madd_v3_v3fl(shi->ao, mesh->ao[face[3]], w[3]);
  			madd_v3_v3fl(shi->env, mesh->env[face[3]], w[3]);
  			madd_v3_v3fl(shi->indirect, mesh->indirect[face[3]], w[3]);
  		}
  	}
  	else {
  		shi->ao[0] = 1.0f;
  		shi->ao[1] = 1.0f;
  		shi->ao[2] = 1.0f;
  		zero_v3(shi->env);
  		zero_v3(shi->indirect);
  	}
  }
  
  /* ------------------------- External Functions --------------------------- */
  
  static void *exec_strandsurface_sample(void *data)
  {
  	OcclusionThread *othread = (OcclusionThread *)data;
  	Render *re = othread->re;
  	StrandSurface *mesh = othread->mesh;
  	float ao[3], env[3], indirect[3], co[3], n[3], *co1, *co2, *co3, *co4;
  	int a, *face;
  
  	for (a = othread->begin; a < othread->end; a++) {
  		face = mesh->face[a];
  		co1 = mesh->co[face[0]];
  		co2 = mesh->co[face[1]];
  		co3 = mesh->co[face[2]];
  
  		if (face[3]) {
  			co4 = mesh->co[face[3]];
  
  			mid_v3_v3v3(co, co1, co3);
  			normal_quad_v3(n, co1, co2, co3, co4);
  		}
  		else {
  			mid_v3_v3v3v3(co, co1, co2, co3);
  			normal_tri_v3(n, co1, co2, co3);
  		}
  		negate_v3(n);
  
  		sample_occ_tree(re, re->occlusiontree, NULL, co, n, othread->thread, 0, ao, env, indirect);
  		copy_v3_v3(othread->faceao[a], ao);
  		copy_v3_v3(othread->faceenv[a], env);
  		copy_v3_v3(othread->faceindirect[a], indirect);
  	}
  
  	return NULL;
  }
  
  void make_occ_tree(Render *re)
  {
  	OcclusionThread othreads[BLENDER_MAX_THREADS];
  	OcclusionTree *tree;
  	StrandSurface *mesh;
  	ListBase threads;
  	float ao[3], env[3], indirect[3], (*faceao)[3], (*faceenv)[3], (*faceindirect)[3];
  	int a, totface, totthread, *face, *count;
  
  	/* ugly, needed for occ_face */
  	R = *re;
  
  	re->i.infostr = IFACE_("Occlusion preprocessing");
  	re->stats_draw(re->sdh, &re->i);
  	
  	re->occlusiontree = tree = occ_tree_build(re);
  	
  	if (tree && !re->test_break(re->tbh)) {
  		if (re->wrld.ao_approx_passes > 0)
  			occ_compute_passes(re, tree, re->wrld.ao_approx_passes);
  		if (tree->doindirect && (re->wrld.mode & WO_INDIRECT_LIGHT))
  			occ_compute_bounces(re, tree, re->wrld.ao_indirect_bounces);
  
  		for (mesh = re->strandsurface.first; mesh; mesh = mesh->next) {
  			if (!mesh->face || !mesh->co || !mesh->ao)
  				continue;
  
  			count = MEM_callocN(sizeof(int) * mesh->totvert, "OcclusionCount");
  			faceao = MEM_callocN(sizeof(float) * 3 * mesh->totface, "StrandSurfFaceAO");
  			faceenv = MEM_callocN(sizeof(float) * 3 * mesh->totface, "StrandSurfFaceEnv");
  			faceindirect = MEM_callocN(sizeof(float) * 3 * mesh->totface, "StrandSurfFaceIndirect");
  
  			totthread = (mesh->totface > 10000) ? re->r.threads : 1;
  			totface = mesh->totface / totthread;
  			for (a = 0; a < totthread; a++) {
  				othreads[a].re = re;
  				othreads[a].faceao = faceao;
  				othreads[a].faceenv = faceenv;
  				othreads[a].faceindirect = faceindirect;
  				othreads[a].thread = a;
  				othreads[a].mesh = mesh;
  				othreads[a].begin = a * totface;
  				othreads[a].end = (a == totthread - 1) ? mesh->totface : (a + 1) * totface;
  			}
  
  			if (totthread == 1) {
  				exec_strandsurface_sample(&othreads[0]);
  			}
  			else {
  				BLI_init_threads(&threads, exec_strandsurface_sample, totthread);
  
  				for (a = 0; a < totthread; a++)
  					BLI_insert_thread(&threads, &othreads[a]);
  
  				BLI_end_threads(&threads);
  			}
  
  			for (a = 0; a < mesh->totface; a++) {
  				face = mesh->face[a];
  
  				copy_v3_v3(ao, faceao[a]);
  				copy_v3_v3(env, faceenv[a]);
  				copy_v3_v3(indirect, faceindirect[a]);
  
  				add_v3_v3(mesh->ao[face[0]], ao);
  				add_v3_v3(mesh->env[face[0]], env);
  				add_v3_v3(mesh->indirect[face[0]], indirect);
  				count[face[0]]++;
  				add_v3_v3(mesh->ao[face[1]], ao);
  				add_v3_v3(mesh->env[face[1]], env);
  				add_v3_v3(mesh->indirect[face[1]], indirect);
  				count[face[1]]++;
  				add_v3_v3(mesh->ao[face[2]], ao);
  				add_v3_v3(mesh->env[face[2]], env);
  				add_v3_v3(mesh->indirect[face[2]], indirect);
  				count[face[2]]++;
  
  				if (face[3]) {
  					add_v3_v3(mesh->ao[face[3]], ao);
  					add_v3_v3(mesh->env[face[3]], env);
  					add_v3_v3(mesh->indirect[face[3]], indirect);
  					count[face[3]]++;
  				}
  			}
  
  			for (a = 0; a < mesh->totvert; a++) {
  				if (count[a]) {
  					mul_v3_fl(mesh->ao[a], 1.0f / count[a]);
  					mul_v3_fl(mesh->env[a], 1.0f / count[a]);
  					mul_v3_fl(mesh->indirect[a], 1.0f / count[a]);
  				}
  			}
  
  			MEM_freeN(count);
  			MEM_freeN(faceao);
  			MEM_freeN(faceenv);
  			MEM_freeN(faceindirect);
  		}
  	}
  }
  
  void free_occ(Render *re)
  {
  	if (re->occlusiontree) {
  		occ_free_tree(re->occlusiontree);
  		re->occlusiontree = NULL;
  	}
  }
  
  void sample_occ(Render *re, ShadeInput *shi)
  {
  	OcclusionTree *tree = re->occlusiontree;
  	OcclusionCache *cache;
  	OcclusionCacheSample *sample;
  	OccFace exclude;
  	int onlyshadow;
  
  	if (tree) {
  		if (shi->strand) {
  			sample_occ_surface(shi);
  		}
  		/* try to get result from the cache if possible */
  		else if (shi->depth != 0 || !sample_occ_cache(tree, shi->co, shi->vno, shi->xs, shi->ys, shi->thread, shi->ao, shi->env, shi->indirect)) {
  			/* no luck, let's sample the occlusion */
  			exclude.obi = shi->obi - re->objectinstance;
  			exclude.facenr = shi->vlr->index;
  			onlyshadow = (shi->mat->mode & MA_ONLYSHADOW);
  			sample_occ_tree(re, tree, &exclude, shi->co, shi->vno, shi->thread, onlyshadow, shi->ao, shi->env, shi->indirect);
  
  			/* fill result into sample, each time */
  			if (tree->cache) {
  				cache = &tree->cache[shi->thread];
  
  				if (cache->sample && cache->step) {
  					sample = &cache->sample[(shi->ys - cache->y) * cache->w + (shi->xs - cache->x)];
  					copy_v3_v3(sample->co, shi->co);
  					copy_v3_v3(sample->n, shi->vno);
  					copy_v3_v3(sample->ao, shi->ao);
  					copy_v3_v3(sample->env, shi->env);
  					copy_v3_v3(sample->indirect, shi->indirect);
  					sample->intensity = max_fff(sample->ao[0], sample->ao[1], sample->ao[2]);
  					sample->intensity = max_ff(sample->intensity, max_fff(sample->env[0], sample->env[1], sample->env[2]));
  					sample->intensity = max_ff(sample->intensity, max_fff(sample->indirect[0], sample->indirect[1], sample->indirect[2]));
  					sample->dist2 = dot_v3v3(shi->dxco, shi->dxco) + dot_v3v3(shi->dyco, shi->dyco);
  					sample->filled = 1;
  				}
  			}
  		}
  	}
  	else {
  		shi->ao[0] = 1.0f;
  		shi->ao[1] = 1.0f;
  		shi->ao[2] = 1.0f;
  
  		shi->env[0] = 0.0f;
  		shi->env[1] = 0.0f;
  		shi->env[2] = 0.0f;
  
  		shi->indirect[0] = 0.0f;
  		shi->indirect[1] = 0.0f;
  		shi->indirect[2] = 0.0f;
  	}
  }
  
  void cache_occ_samples(Render *re, RenderPart *pa, ShadeSample *ssamp)
  {
  	OcclusionTree *tree = re->occlusiontree;
  	PixStr ps;
  	OcclusionCache *cache;
  	OcclusionCacheSample *sample;
  	OccFace exclude;
  	ShadeInput *shi;
  	intptr_t *rd = NULL;
  	int *ro = NULL, *rp = NULL, *rz = NULL, onlyshadow;
  	int x, y, step = CACHE_STEP;
  
  	if (!tree->cache)
  		return;
  
  	cache = &tree->cache[pa->thread];
  	cache->w = pa->rectx;
  	cache->h = pa->recty;
  	cache->x = pa->disprect.xmin;
  	cache->y = pa->disprect.ymin;
  	cache->step = step;
  	cache->sample = MEM_callocN(sizeof(OcclusionCacheSample) * cache->w * cache->h, "OcclusionCacheSample");
  	sample = cache->sample;
  
  	if (re->osa) {
  		rd = pa->rectdaps;
  	}
  	else {
  		/* fake pixel struct for non-osa */
  		ps.next = NULL;
  		ps.mask = 0xFFFF;
  
  		ro = pa->recto;
  		rp = pa->rectp;
  		rz = pa->rectz;
  	}
  
  	/* compute a sample at every step pixels */
  	for (y = pa->disprect.ymin; y < pa->disprect.ymax; y++) {
  		for (x = pa->disprect.xmin; x < pa->disprect.xmax; x++, sample++, rd++, ro++, rp++, rz++) {
  			if (!(((x - pa->disprect.xmin + step) % step) == 0 || x == pa->disprect.xmax - 1))
  				continue;
  			if (!(((y - pa->disprect.ymin + step) % step) == 0 || y == pa->disprect.ymax - 1))
  				continue;
  
  			if (re->osa) {
  				if (!*rd) continue;
  
  				shade_samples_fill_with_ps(ssamp, (PixStr *)(*rd), x, y);
  			}
  			else {
  				if (!*rp) continue;
  
  				ps.obi = *ro;
  				ps.facenr = *rp;
  				ps.z = *rz;
  				shade_samples_fill_with_ps(ssamp, &ps, x, y);
  			}
  
  			shi = ssamp->shi;
  			if (shi->vlr) {
  				onlyshadow = (shi->mat->mode & MA_ONLYSHADOW);
  				exclude.obi = shi->obi - re->objectinstance;
  				exclude.facenr = shi->vlr->index;
  				sample_occ_tree(re, tree, &exclude, shi->co, shi->vno, shi->thread, onlyshadow, shi->ao, shi->env, shi->indirect);
  
  				copy_v3_v3(sample->co, shi->co);
  				copy_v3_v3(sample->n, shi->vno);
  				copy_v3_v3(sample->ao, shi->ao);
  				copy_v3_v3(sample->env, shi->env);
  				copy_v3_v3(sample->indirect, shi->indirect);
  				sample->intensity = max_fff(sample->ao[0], sample->ao[1], sample->ao[2]);
  				sample->intensity = max_ff(sample->intensity, max_fff(sample->env[0], sample->env[1], sample->env[2]));
  				sample->intensity = max_ff(sample->intensity, max_fff(sample->indirect[0], sample->indirect[1], sample->indirect[2]));
  				sample->dist2 = dot_v3v3(shi->dxco, shi->dxco) + dot_v3v3(shi->dyco, shi->dyco);
  				sample->x = shi->xs;
  				sample->y = shi->ys;
  				sample->filled = 1;
  			}
  
  			if (re->test_break(re->tbh))
  				break;
  		}
  	}
  }
  
  void free_occ_samples(Render *re, RenderPart *pa)
  {
  	OcclusionTree *tree = re->occlusiontree;
  	OcclusionCache *cache;
  
  	if (tree->cache) {
  		cache = &tree->cache[pa->thread];
  
  		if (cache->sample)
  			MEM_freeN(cache->sample);
  
  		cache->w = 0;
  		cache->h = 0;
  		cache->step = 0;
  	}
  }
  
  