Subsurface scattering:

Documentation on the settings, known limitations and implementation
info can be found here:
http://www.blender.org/development/current-projects/changes-since-243/subsurface-scattering/

20 files changed, 1794 insertions, 37 deletions

diff --git a/source/blender/blenkernel/BKE_material.h b/source/blender/blenkernel/BKE_material.h
index 93260bac73d..2add4b95080 100644
--- a/source/blender/blenkernel/BKE_material.h
+++ b/source/blender/blenkernel/BKE_material.h
@@ -62,6 +62,8 @@ void init_render_materials(int, float *);
 void end_render_material(struct Material *);
 void end_render_materials(void);
 
+int material_in_material(struct Material *parmat, struct Material *mat);
+
 void automatname(struct Material *);
 void delete_material_index(void);            
 
diff --git a/source/blender/blenkernel/intern/material.c b/source/blender/blenkernel/intern/material.c
index d33cc7bd1c4..56b8307020a 100644
--- a/source/blender/blenkernel/intern/material.c
+++ b/source/blender/blenkernel/intern/material.c
@@ -144,6 +144,20 @@ void init_material(Material *ma)
 	ma->pr_lamp= 3;			/* two lamps, is bits */
 	ma->pr_type= MA_SPHERE;
 
+	ma->sss_radius[0]= 1.0f;
+	ma->sss_radius[1]= 1.0f;
+	ma->sss_radius[2]= 1.0f;
+	ma->sss_col[0]= 0.8f;
+	ma->sss_col[1]= 0.8f;
+	ma->sss_col[2]= 0.8f;
+	ma->sss_error= 0.05f;
+	ma->sss_scale= 0.1f;
+	ma->sss_ior= 1.3f;
+	ma->sss_colfac= 1.0f;
+	ma->sss_texfac= 0.0f;
+	ma->sss_front= 1.0f;
+	ma->sss_back= 1.0f;
+
 	ma->mode= MA_TRACEBLE|MA_SHADBUF|MA_SHADOW|MA_RADIO|MA_RAYBIAS|MA_TANGENT_STR;
 }
 
@@ -702,6 +716,33 @@ void end_render_materials(void)
 			end_render_material(ma);
 }
 
+static int material_in_nodetree(bNodeTree *ntree, Material *mat)
+{
+	bNode *node;
+
+	for(node=ntree->nodes.first; node; node= node->next) {
+		if(node->id && GS(node->id->name)==ID_MA) {
+			if(node->id==(ID*)mat)
+				return 1;
+		}
+		else if(node->type==NODE_GROUP)
+			if(material_in_nodetree((bNodeTree*)node->id, mat))
+				return 1;
+	}
+
+	return 0;
+}
+
+int material_in_material(Material *parmat, Material *mat)
+{
+	if(parmat==mat)
+		return 1;
+	else if(parmat->nodetree && parmat->use_nodes)
+		return material_in_nodetree(parmat->nodetree, mat);
+	else
+		return 0;
+}
+	
 /* ****************** */
 
 char colname_array[125][20]= {
diff --git a/source/blender/blenloader/intern/readfile.c b/source/blender/blenloader/intern/readfile.c
index ace436cf04c..c283e5c9fcf 100644
--- a/source/blender/blenloader/intern/readfile.c
+++ b/source/blender/blenloader/intern/readfile.c
@@ -6387,10 +6387,36 @@ static void do_versions(FileData *fd, Library *lib, Main *main)
 	if(main->versionfile <= 243) {
 		Object *ob= main->object.first;
 		Camera *cam = main->camera.first;
+		Material *ma;
 		
 		for(; cam; cam= cam->id.next) {
 			cam->angle= 360.0f * atan(16.0f/cam->lens) / M_PI;
 		}
+
+		for(ma=main->mat.first; ma; ma= ma->id.next) {
+			if(ma->sss_scale==0.0f) {
+				ma->sss_radius[0]= 1.0f;
+				ma->sss_radius[1]= 1.0f;
+				ma->sss_radius[2]= 1.0f;
+				ma->sss_col[0]= 0.8f;
+				ma->sss_col[1]= 0.8f;
+				ma->sss_col[2]= 0.8f;
+				ma->sss_error= 0.05f;
+				ma->sss_scale= 0.1f;
+				ma->sss_ior= 1.3f;
+				ma->sss_colfac= 1.0f;
+				ma->sss_texfac= 0.0f;
+			}
+			if(ma->sss_front==0 && ma->sss_back==0) {
+				ma->sss_front= 1.0f;
+				ma->sss_back= 1.0f;
+			}
+			if(ma->sss_col[0]==0 && ma->sss_col[1]==0 && ma->sss_col[2]==0) {
+				ma->sss_col[0]= ma->r;
+				ma->sss_col[1]= ma->g;
+				ma->sss_col[2]= ma->b;
+			}
+		}
 		
 		for(; ob; ob= ob->id.next) {
 			bDeformGroup *curdef;
diff --git a/source/blender/makesdna/DNA_material_types.h b/source/blender/makesdna/DNA_material_types.h
index 7b948b93b10..0fad110d64b 100644
--- a/source/blender/makesdna/DNA_material_types.h
+++ b/source/blender/makesdna/DNA_material_types.h
@@ -111,6 +111,12 @@ typedef struct Material {
 	float fhdist, xyfrict;
 	short dynamode, pad2;
 
+	float sss_radius[3], sss_col[3];
+	float sss_error, sss_scale, sss_ior;
+	float sss_colfac, sss_texfac;
+	float sss_front, sss_back;
+	short sss_flag, sss_preset;
+
 	/* yafray: absorption color, dispersion parameters and material preset menu */
 	float YF_ar, YF_ag, YF_ab, YF_dscale, YF_dpwr;
 	int YF_dsmp, YF_preset, YF_djit;
@@ -268,5 +274,8 @@ typedef struct Material {
 /* pr_back */
 #define MA_DARK			1
 
+/* sss_flag */
+#define MA_DIFF_SSS		1
+
 #endif
 
diff --git a/source/blender/render/extern/include/RE_pipeline.h b/source/blender/render/extern/include/RE_pipeline.h
index bd210501330..824050d802f 100644
--- a/source/blender/render/extern/include/RE_pipeline.h
+++ b/source/blender/render/extern/include/RE_pipeline.h
@@ -173,7 +173,7 @@ void RE_DataBase_ApplyWindow(struct Render *re);
 void RE_set_max_threads(int threads);
 
 /* the main processor, assumes all was set OK! */
-void RE_TileProcessor(struct Render *re, int firsttile);
+void RE_TileProcessor(struct Render *re, int firsttile, int threaded);
 
 /* only RE_NewRender() needed, main Blender render calls */
 void RE_BlenderFrame(struct Render *re, struct Scene *scene, int frame);
diff --git a/source/blender/render/intern/include/render_types.h b/source/blender/render/intern/include/render_types.h
index 609c391d599..45aff8e8b60 100644
--- a/source/blender/render/intern/include/render_types.h
+++ b/source/blender/render/intern/include/render_types.h
@@ -68,9 +68,12 @@ typedef struct RenderPart
 	/* result of part rendering */
 	RenderResult *result;
 	
-	int *rectp;			/* polygon index table */
+	int *rectp;						/* polygon index table */
 	int *rectz;						/* zbuffer */
 	long *rectdaps;					/* delta acum buffer for pixel structs */
+	int *rectbackp;					/* polygon index table for backside sss */
+	int *rectbackz;					/* zbuffer for backside sss */
+	long *rectall;					/* buffer for all faces for sss */
 	
 	rcti disprect;					/* part coordinates within total picture */
 	int rectx, recty;				/* the size */
@@ -169,6 +172,10 @@ struct Render
 	
 	struct GHash *orco_hash;
 
+	struct GHash *sss_hash;
+	ListBase *sss_points;
+	struct Material *sss_mat;
+
 	ListBase customdata_names;
 
 	/* arena for allocating data for use during render, for
@@ -262,7 +269,7 @@ typedef struct VlakRen {
 	unsigned int lay;
 	float n[3];
 	struct Material *mat;
-	char snproj, puno;
+	char noflag, puno;
 	char flag, ec;
 	RadFace *radface;
 	Object *ob;
@@ -393,7 +400,11 @@ typedef struct LampRen {
 /* vertex normals are tangent or view-corrected vector, for hair strands */
 #define R_TANGENT		128		
 
-
+/* vlakren->noflag (char) */
+#define R_SNPROJ_X		1
+#define R_SNPROJ_Y		2
+#define R_SNPROJ_Z		4
+#define R_FLIPPED_NO	8
 
 
 
diff --git a/source/blender/render/intern/include/rendercore.h b/source/blender/render/intern/include/rendercore.h
index 9f22480f0cb..651fd423e5d 100644
--- a/source/blender/render/intern/include/rendercore.h
+++ b/source/blender/render/intern/include/rendercore.h
@@ -87,6 +87,8 @@ void zbufshadeDA(void);	/* Delta Accum Pixel Struct */
 void zbufshade_tile(struct RenderPart *pa);
 void zbufshadeDA_tile(struct RenderPart *pa);
 
+void zbufshade_sss_tile(struct RenderPart *pa);
+
 /* -------- ray.c ------- */
 
 extern void freeoctree(Render *re);
diff --git a/source/blender/render/intern/include/sss.h b/source/blender/render/intern/include/sss.h
new file mode 100644
index 00000000000..660fb544227
--- /dev/null
+++ b/source/blender/render/intern/include/sss.h
@@ -0,0 +1,65 @@
+/**
+ * $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) 2007 Blender Foundation.
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): none yet.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+#ifndef SSS_H
+#define SSS_H
+
+/* Generic multiple scattering API */
+
+struct ScatterSettings;
+typedef struct ScatterSettings ScatterSettings;
+
+struct ScatterTree;
+typedef struct ScatterTree ScatterTree;
+
+ScatterSettings *scatter_settings_new(float refl, float radius, float ior,
+	float reflfac, float frontweight, float backweight);
+void scatter_settings_free(ScatterSettings *ss);
+
+ScatterTree *scatter_tree_new(ScatterSettings *ss[3], float scale, float error,
+	float (*co)[3], float (*color)[3], float *area, int totpoint);
+void scatter_tree_build(ScatterTree *tree);
+void scatter_tree_sample(ScatterTree *tree, float *co, float *color);
+void scatter_tree_free(ScatterTree *tree);
+
+/* Internal renderer API */
+
+struct Render;
+struct Material;
+struct VlakRen;
+
+void make_sss_tree(struct Render *re);
+void sss_add_points(Render *re, float (*co)[3], float (*color)[3], float *area, int totpoint);
+void free_sss(struct Render *re);
+
+int sample_sss(struct Render *re, struct Material *mat, float *co, float *col);
+int has_sss_tree(struct Render *re, struct Material *mat);
+
+#endif /*SSS_H*/
+
diff --git a/source/blender/render/intern/include/zbuf.h b/source/blender/render/intern/include/zbuf.h
index c308ae25bce..ed3f93adfdd 100644
--- a/source/blender/render/intern/include/zbuf.h
+++ b/source/blender/render/intern/include/zbuf.h
@@ -51,6 +51,7 @@ void zbuffer_shadow(struct Render *re, struct LampRen *lar, int *rectz, int size
 void zbuffer_solid(struct RenderPart *pa, unsigned int layer, short layflag);
 unsigned short *zbuffer_transp_shade(struct RenderPart *pa, struct RenderLayer *rl, float *pass);
 void convert_zbuf_to_distbuf(struct RenderPart *pa, struct RenderLayer *rl);
+void zbuffer_sss(RenderPart *pa, unsigned int lay, void *handle, void (*func)(void *, int, int, int, int));
 
 typedef struct APixstr {
     unsigned short mask[4]; /* jitter mask */
@@ -85,6 +86,9 @@ typedef struct ZSpan {
 	int polygon_offset;						/* offset in Z */
 	float shad_alpha;						/* copy from material, used by irregular shadbuf */
 	int mask, apsmcounter;					/* in use by apixbuf */
+
+	void *sss_handle;						/* used by sss */
+	void (*sss_func)(void *, int, int, int, int);
 	
 	void (*zbuffunc)(struct ZSpan *, int, float *, float *, float *, float *);
 	void (*zbuflinefunc)(struct ZSpan *, int, float *, float *);
diff --git a/source/blender/render/intern/source/convertblender.c b/source/blender/render/intern/source/convertblender.c
index e93a509efd1..afa9aa0b2e4 100644
--- a/source/blender/render/intern/source/convertblender.c
+++ b/source/blender/render/intern/source/convertblender.c
@@ -103,6 +103,7 @@
 #include "shadbuf.h"
 #include "shading.h"
 #include "texture.h"
+#include "sss.h"
 #include "zbuf.h"
 
 #ifndef DISABLE_YAFRAY /* disable yafray */
@@ -2990,6 +2991,8 @@ void RE_Database_Free(Render *re)
 	}
 	
 	if(re->r.mode & R_RAYTRACE) freeoctree(re);
+
+	free_sss(re);
 	
 	re->totvlak=re->totvert=re->totlamp=re->tothalo= 0;
 	re->i.convertdone= 0;
@@ -3480,6 +3483,11 @@ void RE_Database_FromScene(Render *re, Scene *scene, int use_camera_view)
 		
 		if(!re->test_break())
 			project_renderdata(re, projectverto, re->r.mode & R_PANORAMA, 0);
+
+		/* SSS */
+		if(!re->test_break())
+			if (re->r.renderer==R_INTERN)
+				make_sss_tree(re);
 	}
 	
 	if(re->test_break())
diff --git a/source/blender/render/intern/source/envmap.c b/source/blender/render/intern/source/envmap.c
index c1f756894fe..ec55bc0a5e2 100644
--- a/source/blender/render/intern/source/envmap.c
+++ b/source/blender/render/intern/source/envmap.c
@@ -410,7 +410,7 @@ static void render_envmap(Render *re, EnvMap *env)
 		env_set_imats(envre);
 				
 		if(re->test_break()==0) {
-			RE_TileProcessor(envre, 0);
+			RE_TileProcessor(envre, 0, 0);
 		}
 		
 		/* rotate back */
diff --git a/source/blender/render/intern/source/pipeline.c b/source/blender/render/intern/source/pipeline.c
index d93ee448f62..93c80749bca 100644
--- a/source/blender/render/intern/source/pipeline.c
+++ b/source/blender/render/intern/source/pipeline.c
@@ -1163,16 +1163,17 @@ static void *do_part_thread(void *pa_v)
 	
 	/* need to return nicely all parts on esc */
 	if(R.test_break()==0) {
-		
 		pa->result= new_render_result(&R, &pa->disprect, pa->crop, RR_USEMEM);
-		
-		if(R.osa)
+
+		if(R.sss_points)
+			zbufshade_sss_tile(pa);
+		else if(R.osa)
 			zbufshadeDA_tile(pa);
 		else
 			zbufshade_tile(pa);
 		
 		/* merge too on break! */
-		if(R.result->exrhandle)
+		if(!R.sss_points && R.result->exrhandle)
 			save_render_result_tile(&R, pa);
 		else
 			merge_render_result(R.result, pa->result);
@@ -1464,7 +1465,7 @@ static void threaded_tile_processor(Render *re)
 }
 
 /* currently only called by preview renders and envmap */
-void RE_TileProcessor(Render *re, int firsttile)
+void RE_TileProcessor(Render *re, int firsttile, int threaded)
 {
 	/* the partsdone variable has to be reset to firsttile, to survive esc before it was set to zero */
 	
@@ -1472,8 +1473,10 @@ void RE_TileProcessor(Render *re, int firsttile)
 
 	re->i.starttime= PIL_check_seconds_timer();
 
-	//	threaded_tile_processor(re);
-	render_tile_processor(re, firsttile);
+	if(threaded)
+		threaded_tile_processor(re);
+	else
+		render_tile_processor(re, firsttile);
 		
 	re->i.lastframetime= PIL_check_seconds_timer()- re->i.starttime;
 	re->stats_draw(&re->i);
diff --git a/source/blender/render/intern/source/rendercore.c b/source/blender/render/intern/source/rendercore.c
index 10201087f2d..1c5f56362bc 100644
--- a/source/blender/render/intern/source/rendercore.c
+++ b/source/blender/render/intern/source/rendercore.c
@@ -64,6 +64,7 @@
 #include "pixelshading.h"
 #include "shadbuf.h"
 #include "shading.h"
+#include "sss.h"
 #include "zbuf.h"
 
 #include "PIL_time.h"
@@ -1113,6 +1114,298 @@ void zbufshade_tile(RenderPart *pa)
 	MEM_freeN(pa->clipflag); pa->clipflag= NULL;
 }
 
+/* SSS preprocess tile render, fully threadable */
+typedef struct ZBufSSSHandle {
+	RenderPart *pa;
+	ListBase psmlist;
+	int totps;
+} ZBufSSSHandle;
+
+static void addps_sss(void *cb_handle, int facenr, int x, int y, int z)
+{
+	ZBufSSSHandle *handle = cb_handle;
+	RenderPart *pa= handle->pa;
+
+	if (pa->rectall) {
+		long *rs= pa->rectall + pa->rectx*y + x;
+
+		addps(&handle->psmlist, rs, facenr, z, 0);
+		handle->totps++;
+	}
+	if (pa->rectz) {
+		int *rz= pa->rectz + pa->rectx*y + x;
+		int *rp= pa->rectp + pa->rectx*y + x;
+
+		if (z < *rz) {
+			if(*rp == 0)
+				handle->totps++;
+			*rz= z;
+			*rp= facenr;
+		}
+	}
+	if (pa->rectbackz) {
+		int *rz= pa->rectbackz + pa->rectx*y + x;
+		int *rp= pa->rectbackp + pa->rectx*y + x;
+
+		if (z >= *rz) {
+			if(*rp == 0)
+				handle->totps++;
+			*rz= z;
+			*rp= facenr;
+		}
+	}
+}
+
+static void shade_sample_sss(ShadeSample *ssamp, Material *mat, VlakRen *vlr, int quad, int x, int y, float z, float *co, float *color, float *area)
+{
+	ShadeInput *shi= ssamp->shi;
+	ShadeResult shr;
+	float texfac, orthoarea, nor[3];
+
+	/* normal flipping must be disabled to make back scattering work, so that
+	   backside faces actually face any lighting from the back */
+	shi->puno= 0;
+	
+	/* cache for shadow */
+	shi->samplenr++;
+	
+	if(quad) 
+		shade_input_set_triangle_i(shi, vlr, 0, 2, 3);
+	else
+		shade_input_set_triangle_i(shi, vlr, 0, 1, 2);
+
+	/* we don't want flipped normals, they screw up back scattering */
+	if(vlr->noflag & R_FLIPPED_NO) {
+		shi->facenor[0]= -shi->facenor[0];
+		shi->facenor[1]= -shi->facenor[1];
+		shi->facenor[2]= -shi->facenor[2];
+	}
+
+	/* we estimate the area here using shi->dxco and shi->dyco. we need to
+	   enabled shi->osatex these are filled. we compute two areas, one with
+	   the normal pointed at the camera and one with the original normal, and
+	   then clamp to avoid a too large contribution from a single pixel */
+	shi->osatex= 1;
+
+	VECCOPY(nor, shi->facenor);
+	calc_view_vector(shi->facenor, x, y);
+	Normalize(shi->facenor);
+	shade_input_set_viewco(shi, x, y, z);
+	orthoarea= VecLength(shi->dxco)*VecLength(shi->dyco);
+
+	VECCOPY(shi->facenor, nor);
+	shade_input_set_viewco(shi, x, y, z);
+	*area= VecLength(shi->dxco)*VecLength(shi->dyco);
+	*area= MIN2(*area, 2.0f*orthoarea);
+
+	shi->osatex= 0;
+
+	shade_input_set_uv(shi);
+	shade_input_set_normals(shi);
+
+	/* if nodetree, use the material that we are currently preprocessing
+	   instead of the node material */
+	if(shi->mat->nodetree && shi->mat->use_nodes)
+		shi->mat= mat;
+
+	/* init material vars */
+	// note, keep this synced with render_types.h
+	memcpy(&shi->r, &shi->mat->r, 23*sizeof(float));
+	shi->har= shi->mat->har;
+	
+	/* render */
+	shade_input_set_shade_texco(shi);
+	
+	shade_samples_do_AO(ssamp);
+	shade_material_loop(shi, &shr);
+	
+	VECCOPY(co, shi->co);
+	VECCOPY(color, shr.combined);
+
+	/* texture blending */
+	texfac= shi->mat->sss_texfac;
+
+	if(texfac == 0.0f) {
+		if(shr.col[0]!=0.0f) color[0] /= shr.col[0];
+		if(shr.col[1]!=0.0f) color[1] /= shr.col[1];
+		if(shr.col[2]!=0.0f) color[2] /= shr.col[2];
+	}
+	else if(texfac != 1.0f) {
+		if(shr.col[0]!=0.0f) color[0] *= pow(shr.col[0], texfac)/shr.col[0];
+		if(shr.col[1]!=0.0f) color[1] *= pow(shr.col[1], texfac)/shr.col[1];
+		if(shr.col[2]!=0.0f) color[2] *= pow(shr.col[2], texfac)/shr.col[2];
+	}
+}
+
+void zbufshade_sss_tile(RenderPart *pa)
+{
+	ShadeSample ssamp;
+	ZBufSSSHandle handle;
+	RenderResult *rr= pa->result;
+	RenderLayer *rl= rr->layers.first;
+	VlakRen *vlr;
+	Material *mat= R.sss_mat;
+	float (*co)[3], (*color)[3], *area, *fcol= rl->rectf;
+	int x, y, seed, quad, totpoint;
+#if 0
+	PixStr *ps;
+	long *rs;
+	int z;
+#else
+	int *rz, *rp, *rbz, *rbp;
+#endif
+
+	set_part_zbuf_clipflag(pa);
+
+	/* setup pixelstr list and buffer for zbuffering */
+	handle.pa= pa;
+	handle.totps= 0;
+
+#if 0
+	handle.psmlist.first= handle.psmlist.last= NULL;
+	addpsmain(&handle.psmlist);
+
+	pa->rectall= MEM_callocN(sizeof(long)*pa->rectx*pa->recty+4, "rectall");
+#else
+	pa->rectp= MEM_mallocN(sizeof(int)*pa->rectx*pa->recty, "rectp");
+	pa->rectz= MEM_mallocN(sizeof(int)*pa->rectx*pa->recty, "rectz");
+	pa->rectbackp= MEM_mallocN(sizeof(int)*pa->rectx*pa->recty, "rectbackp");
+	pa->rectbackz= MEM_mallocN(sizeof(int)*pa->rectx*pa->recty, "rectbackz");
+#endif
+
+	/* create the pixelstrs to be used later */
+	zbuffer_sss(pa, rl->lay, &handle, addps_sss);
+
+	co= MEM_mallocN(sizeof(float)*3*handle.totps, "SSSCo");
+	color= MEM_mallocN(sizeof(float)*3*handle.totps, "SSSColor");
+	area= MEM_mallocN(sizeof(float)*handle.totps, "SSSArea");
+
+#if 0
+	/* create ISB (does not work currently!) */
+	if(R.r.mode & R_SHADOW)
+		ISB_create(pa, NULL);
+#endif
+
+	/* setup shade sample with correct passes */
+	shade_sample_initialize(&ssamp, pa, rl);
+	ssamp.shi[0].passflag= SCE_PASS_DIFFUSE|SCE_PASS_AO|SCE_PASS_RADIO;
+	ssamp.shi[0].passflag |= SCE_PASS_RGBA;
+	ssamp.shi[0].combinedflag= ~(SCE_PASS_SPEC);
+
+	/* initialize scanline updates for main thread */
+	rr->renrect.ymin= 0;
+	rr->renlay= rl;
+	
+	seed= pa->rectx*pa->disprect.ymin;
+#if 0
+	rs= pa->rectall;
+#else
+	rz= pa->rectz;
+	rp= pa->rectp;
+	rbz= pa->rectbackz;
+	rbp= pa->rectbackp;
+#endif
+	totpoint= 0;
+
+	for(y=pa->disprect.ymin; y<pa->disprect.ymax; y++, rr->renrect.ymax++) {
+		for(x=pa->disprect.xmin; x<pa->disprect.xmax; x++, fcol+=4) {
+			/* per pixel fixed seed */
+			BLI_thread_srandom(pa->thread, seed++);
+			
+#if 0
+			if(rs) {
+				/* for each sample in this pixel, shade it */
+				for(ps=(PixStr*)*rs; ps; ps=ps->next) {
+					vlr= RE_findOrAddVlak(&R, (ps->facenr-1) & RE_QUAD_MASK);
+					quad= (ps->facenr & RE_QUAD_OFFS);
+					z= ps->z;
+
+					shade_sample_sss(&ssamp, mat, vlr, quad, x, y, z,
+						co[totpoint], color[totpoint], &area[totpoint]);
+
+					totpoint++;
+
+					VECADD(fcol, fcol, color);
+					fcol[3]= 1.0f;
+				}
+
+				rs++;
+			}
+#else
+			if(rp) {
+				if(*rp != 0) {
+					/* shade front */
+					vlr= RE_findOrAddVlak(&R, (*rp-1) & RE_QUAD_MASK);
+					quad= ((*rp) & RE_QUAD_OFFS);
+
+					shade_sample_sss(&ssamp, mat, vlr, quad, x, y, *rz,
+						co[totpoint], color[totpoint], &area[totpoint]);
+
+					VECADD(fcol, fcol, color[totpoint]);
+					fcol[3]= 1.0f;
+					totpoint++;
+				}
+
+				rp++; rz++;
+			}
+
+			if(rbp) {
+				if(*rbp != 0 && *rbp != *(rp-1)) {
+					/* shade back */
+					vlr= RE_findOrAddVlak(&R, (*rbp-1) & RE_QUAD_MASK);
+					quad= ((*rbp) & RE_QUAD_OFFS);
+
+					shade_sample_sss(&ssamp, mat, vlr, quad, x, y, *rbz,
+						co[totpoint], color[totpoint], &area[totpoint]);
+					
+					/* to indicate this is a back sample */
+					area[totpoint]= -area[totpoint];
+
+					VECADD(fcol, fcol, color[totpoint]);
+					fcol[3]= 1.0f;
+					totpoint++;
+				}
+
+				rbz++; rbp++;
+			}
+#endif
+		}
+
+		if(y&1)
+			if(R.test_break()) break; 
+	}
+
+	/* note: after adding we do not free these arrays, sss keeps them */
+	if(totpoint > 0) {
+		sss_add_points(&R, co, color, area, totpoint);
+	}
+	else {
+		MEM_freeN(co);
+		MEM_freeN(color);
+		MEM_freeN(area);
+	}
+	
+#if 0
+	if(R.r.mode & R_SHADOW)
+		ISB_free(pa);
+#endif
+		
+	/* display active layer */
+	rr->renrect.ymin=rr->renrect.ymax= 0;
+	rr->renlay= render_get_active_layer(&R, rr);
+	
+	MEM_freeN(pa->clipflag); pa->clipflag= NULL;
+#if 0
+	MEM_freeN(pa->rectall); pa->rectall= NULL;
+	freeps(&handle.psmlist);
+#else
+	MEM_freeN(pa->rectz); pa->rectz= NULL;
+	MEM_freeN(pa->rectp); pa->rectp= NULL;
+	MEM_freeN(pa->rectbackz); pa->rectbackz= NULL;
+	MEM_freeN(pa->rectbackp); pa->rectbackp= NULL;
+#endif
+}
+
 /* ------------------------------------------------------------------------ */
 
 static void renderhalo_post(RenderResult *rr, float *rectf, HaloRen *har)	/* postprocess version */
diff --git a/source/blender/render/intern/source/renderdatabase.c b/source/blender/render/intern/source/renderdatabase.c
index 57f0f689871..6d9f0e4eb01 100644
--- a/source/blender/render/intern/source/renderdatabase.c
+++ b/source/blender/render/intern/source/renderdatabase.c
@@ -894,6 +894,8 @@ void set_normalflags(Render *re)
 		if((a1 & 255)==0) vlr= re->vlaknodes[a1>>8].vlak;
 		else vlr++;
 		
+		vlr->noflag= 0;
+
 		/* abuse of this flag... this is code that just sets face normal in direction of camera */
 		/* that convention we should get rid of */
 		if((vlr->flag & R_NOPUNOFLIP)==0) {
@@ -910,6 +912,7 @@ void set_normalflags(Render *re)
 				vlr->n[0]= -vlr->n[0];
 				vlr->n[1]= -vlr->n[1];
 				vlr->n[2]= -vlr->n[2];
+				vlr->noflag |= R_FLIPPED_NO;
 			}
 		}
 		
@@ -924,9 +927,9 @@ void set_normalflags(Render *re)
 		xn= fabs(vlr->n[0]);
 		yn= fabs(vlr->n[1]);
 		zn= fabs(vlr->n[2]);
-		if(zn>=xn && zn>=yn) vlr->snproj= 0;
-		else if(yn>=xn && yn>=zn) vlr->snproj= 1;
-		else vlr->snproj= 2;
+		if(zn>=xn && zn>=yn) vlr->noflag |= R_SNPROJ_X;
+		else if(yn>=xn && yn>=zn) vlr->noflag |= R_SNPROJ_Y;
+		else vlr->noflag |= R_SNPROJ_Z;
 
 	}
 }
diff --git a/source/blender/render/intern/source/shadeinput.c b/source/blender/render/intern/source/shadeinput.c
index d9126a3cdb2..51157cb83f4 100644
--- a/source/blender/render/intern/source/shadeinput.c
+++ b/source/blender/render/intern/source/shadeinput.c
@@ -429,11 +429,11 @@ void shade_input_set_uv(ShadeInput *shi)
 			/* most of this could become re-used for faces */
 			float detsh, t00, t10, t01, t11;
 			
-			if(vlr->snproj==0) {
+			if(vlr->noflag & R_SNPROJ_X) {
 				t00= v3[0]-v1[0]; t01= v3[1]-v1[1];
 				t10= v3[0]-v2[0]; t11= v3[1]-v2[1];
 			}
-			else if(vlr->snproj==1) {
+			else if(vlr->noflag & R_SNPROJ_Y) {
 				t00= v3[0]-v1[0]; t01= v3[2]-v1[2];
 				t10= v3[0]-v2[0]; t11= v3[2]-v2[2];
 			}
@@ -446,7 +446,7 @@ void shade_input_set_uv(ShadeInput *shi)
 			t00*= detsh; t01*=detsh; 
 			t10*=detsh; t11*=detsh;
 			
-			if(vlr->snproj==0) {
+			if(vlr->noflag & R_SNPROJ_X) {
 				shi->u= (shi->co[0]-v3[0])*t11-(shi->co[1]-v3[1])*t10;
 				shi->v= (shi->co[1]-v3[1])*t00-(shi->co[0]-v3[0])*t01;
 				if(shi->osatex) {
@@ -456,7 +456,7 @@ void shade_input_set_uv(ShadeInput *shi)
 					shi->dy_v=  shi->dyco[1]*t00- shi->dyco[0]*t01;
 				}
 			}
-			else if(vlr->snproj==1) {
+			else if(vlr->noflag & R_SNPROJ_Y) {
 				shi->u= (shi->co[0]-v3[0])*t11-(shi->co[2]-v3[2])*t10;
 				shi->v= (shi->co[2]-v3[2])*t00-(shi->co[0]-v3[0])*t01;
 				if(shi->osatex) {
diff --git a/source/blender/render/intern/source/shadeoutput.c b/source/blender/render/intern/source/shadeoutput.c
index c2aeee6aab6..7dd7383c60c 100644
--- a/source/blender/render/intern/source/shadeoutput.c
+++ b/source/blender/render/intern/source/shadeoutput.c
@@ -46,6 +46,7 @@
 #include "pixelblending.h"
 #include "rendercore.h"
 #include "shadbuf.h"
+#include "sss.h"
 #include "texture.h"
 
 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
@@ -1506,6 +1507,35 @@ void shade_lamp_loop(ShadeInput *shi, ShadeResult *shr)
 			/* accumulates in shr->diff and shr->spec and shr->shad (diffuse with shadow!) */
 			shade_one_light(lar, shi, shr, passflag);
 		}
+
+		if(ma->sss_flag & MA_DIFF_SSS) {
+			float sss[3], col[3], texfac= ma->sss_texfac;
+
+			/* this will return false in the preprocess stage */
+			if(sample_sss(&R, ma, shi->co, sss)) {
+				if(texfac==0.0f) {
+					VECCOPY(col, shr->col);
+				}
+				else if(texfac==1.0f) {
+					col[0]= col[1]= col[2]= 1.0f;
+				}
+				else {
+					col[0]= pow(shr->col[0], 1.0f-texfac);
+					col[1]= pow(shr->col[1], 1.0f-texfac);
+					col[2]= pow(shr->col[2], 1.0f-texfac);
+				}
+
+				shr->diff[0]= sss[0]*col[0];
+				shr->diff[1]= sss[1]*col[1];
+				shr->diff[2]= sss[2]*col[2];
+
+				if(shi->combinedflag & SCE_PASS_SHADOW)	{
+					shr->shad[0]= sss[0]*col[0];
+					shr->shad[1]= sss[1]*col[1];
+					shr->shad[2]= sss[2]*col[2];
+				}
+			}
+		}
 		
 		if(shi->combinedflag & SCE_PASS_SHADOW)	
 			VECCOPY(shr->combined, shr->shad) 	/* note, no ';' ! */
@@ -1544,23 +1574,26 @@ void shade_lamp_loop(ShadeInput *shi, ShadeResult *shr)
 	shr->alpha= shi->alpha;
 	
 	/* from now stuff everything in shr->combined: ambient, AO, radio, ramps, exposure */
-	shr->combined[0]+= shi->ambr + shi->r*shi->amb*shi->rad[0];
-	shr->combined[1]+= shi->ambg + shi->g*shi->amb*shi->rad[1];
-	shr->combined[2]+= shi->ambb + shi->b*shi->amb*shi->rad[2];
-	
-	/* add AO in combined? */
-	if(R.wrld.mode & WO_AMB_OCC) {
-		if(shi->combinedflag & SCE_PASS_AO) {
-			float aodiff[3];
-			ambient_occlusion_to_diffuse(shi, aodiff);
-			
-			shr->combined[0] += shi->r*aodiff[0];
-			shr->combined[1] += shi->g*aodiff[1];
-			shr->combined[2] += shi->b*aodiff[2];
+	if(!(ma->sss_flag & MA_DIFF_SSS) || !has_sss_tree(&R, ma)) {
+		shr->combined[0]+= shi->ambr + shi->r*shi->amb*shi->rad[0];
+		shr->combined[1]+= shi->ambg + shi->g*shi->amb*shi->rad[1];
+		shr->combined[2]+= shi->ambb + shi->b*shi->amb*shi->rad[2];
+		
+		/* add AO in combined? */
+		if(R.wrld.mode & WO_AMB_OCC) {
+			if(shi->combinedflag & SCE_PASS_AO) {
+				float aodiff[3];
+				ambient_occlusion_to_diffuse(shi, aodiff);
+				
+				shr->combined[0] += shi->r*aodiff[0];
+				shr->combined[1] += shi->g*aodiff[1];
+				shr->combined[2] += shi->b*aodiff[2];
+			}
 		}
+		
+		if(ma->mode & MA_RAMP_COL) ramp_diffuse_result(shr->combined, shi);
 	}
-	
-	if(ma->mode & MA_RAMP_COL) ramp_diffuse_result(shr->combined, shi);
+
 	if(ma->mode & MA_RAMP_SPEC) ramp_spec_result(shr->spec, shr->spec+1, shr->spec+2, shi);
 	
 	/* refcol is for envmap only */
diff --git a/source/blender/render/intern/source/sss.c b/source/blender/render/intern/source/sss.c
new file mode 100644
index 00000000000..26833562670
--- /dev/null
+++ b/source/blender/render/intern/source/sss.c
@@ -0,0 +1,984 @@
+/* 
+ * $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) 2007 Blender Foundation.
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): none yet.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+/* Possible Improvements:
+   - add fresnel terms
+   - adapt Rd table to scale, now with small scale there are a lot of misses?
+   - possible interesting method: perform sss on all samples in the tree,
+     and then use those values interpolated somehow later. can also do this
+	 filtering on demand for speed. since we are doing things in screen
+	 space now there is an exact correspondence
+   - avoid duplicate shading (filtering points in advance, irradiance cache
+     like lookup?)
+   - lower resolution samples
+*/
+
+#include <math.h>
+#include <string.h>
+#include <stdio.h>
+#include <string.h>
+
+/* external modules: */
+#include "MEM_guardedalloc.h"
+
+#include "BLI_arithb.h"
+#include "BLI_blenlib.h"
+#include "BLI_ghash.h"
+#include "BLI_memarena.h"
+#include "PIL_time.h"
+
+#include "DNA_material_types.h"
+
+#include "BKE_global.h"
+#include "BKE_main.h"
+#include "BKE_material.h"
+#include "BKE_node.h"
+#include "BKE_utildefines.h"
+
+/* this module */
+#include "render_types.h"
+#include "rendercore.h"
+#include "renderdatabase.h" 
+#include "shading.h"
+#include "sss.h"
+#include "zbuf.h"
+
+extern Render R; // meh
+
+/* Generic Multiple Scattering API */
+
+/* Relevant papers:
+	[1] A Practical Model for Subsurface Light Transport
+	[2] A Rapid Hierarchical Rendering Technique for Translucent Materials
+	[3] Efficient Rendering of Local Subsurface Scattering
+	[4] Implementing a skin BSSRDF (or several...)
+*/
+
+/* Defines */
+
+#define RD_TABLE_RANGE 		100.0f
+#define RD_TABLE_RANGE_2	10000.0f
+#define RD_TABLE_SIZE		10000
+
+#define MAX_OCTREE_NODE_POINTS	8
+#define MAX_OCTREE_DEPTH		15
+
+/* Struct Definitions */
+
+struct ScatterSettings {
+	float eta;		/* index of refraction */
+	float sigma_a;	/* absorption coefficient */
+	float sigma_s_; /* reduced scattering coefficient */
+	float sigma_t_; /* reduced extinction coefficient */
+	float sigma;	/* effective extinction coefficient */
+	float Fdr;		/* diffuse fresnel reflectance */
+	float D;		/* diffusion constant */
+	float A;
+	float alpha_;	/* reduced albedo */
+	float zr;		/* distance of virtual lightsource above surface */
+	float zv;		/* distance of virtual lightsource below surface */
+	float ld;		/* mean free path */
+	float ro;		/* diffuse reflectance */
+	float color;
+	float invsigma_t_;
+	float frontweight;
+	float backweight;
+
+	float *tableRd;  /* lookup table to avoid computing Rd */
+	float *tableRd2; /* lookup table to avoid computing Rd for bigger values */
+};
+
+typedef struct ScatterPoint {
+	float co[3];
+	float radiance[3];
+	float area;
+	int back;
+} ScatterPoint;
+
+typedef struct ScatterNode {
+	float co[3];
+	float radiance[3];
+	float backradiance[3];
+	float area;
+	int back;
+
+	int totpoint;
+	ScatterPoint *points;
+
+	float split[3];
+	struct ScatterNode *child[8];
+} ScatterNode;
+
+struct ScatterTree {
+	MemArena *arena;
+
+	ScatterSettings *ss[3];
+	float error, scale;
+
+	ScatterNode *root;
+	ScatterPoint *points;
+	ScatterPoint **refpoints;
+	ScatterPoint **tmppoints;
+	int totpoint;
+	float min[3], max[3];
+};
+
+typedef struct ScatterResult {
+	float radiance[3];
+	float backradiance[3];
+	float rdsum[3];
+} ScatterResult;
+
+/* Functions for BSSRDF reparametrization in to more intuitive parameters,
+   see [2] section 4 for more info. */
+
+static float f_Rd(float alpha_, float A, float ro)
+{
+	float sq;
+
+	sq= sqrt(3.0f*(1.0f - alpha_));
+	return (alpha_/2.0f)*(1.0f + exp((-4.0f/3.0f)*A*sq))*exp(-sq) - ro;
+}
+
+static float compute_reduced_albedo(ScatterSettings *ss)
+{
+	const float tolerance= 1e-8;
+	const int max_iteration_count= 20;
+	float d, fsub, xn_1= 0.0f , xn= 1.0f, fxn, fxn_1;
+	int i;
+
+	/* use secant method to compute reduced albedo using Rd function inverse
+	   with a given reflectance */
+	fxn= f_Rd(xn, ss->A, ss->ro);
+	fxn_1= f_Rd(xn_1, ss->A, ss->ro);
+
+	for(i= 0; i < max_iteration_count; i++) {
+		fsub= (fxn - fxn_1);
+		if(fabs(fsub) < tolerance)
+			break;
+		d= ((xn - xn_1)/fsub)*fxn;
+		if(fabs(d) < tolerance)
+			break;
+
+		xn_1= xn;
+		fxn_1= fxn;
+		xn= xn - d;
+
+		if(xn > 1.0f) xn= 1.0f;
+		if(xn_1 > 1.0f) xn_1= 1.0f;
+		
+		fxn= f_Rd(xn, ss->A, ss->ro);
+    }
+
+	/* avoid division by zero later */
+	if(xn <= 0.0f)
+		xn= 0.00001f;
+
+    return xn;
+}
+
+/* Exponential falloff functions */
+
+static float Rd_rsquare(ScatterSettings *ss, float rr)
+{
+	float sr, sv, Rdr, Rdv;
+
+	sr= sqrt(rr + ss->zr*ss->zr);
+	sv= sqrt(rr + ss->zv*ss->zv);
+
+	Rdr= ss->zr*(1.0f + ss->sigma*sr)*exp(-ss->sigma*sr)/(sr*sr*sr);
+	Rdv= ss->zv*(1.0f + ss->sigma*sv)*exp(-ss->sigma*sv)/(sv*sv*sv);
+
+	return /*ss->alpha_*/(1.0f/(4.0f*M_PI))*(Rdr + Rdv);
+}
+
+static float Rd(ScatterSettings *ss, float r)
+{
+	return Rd_rsquare(ss, r*r);
+}
+
+/* table lookups for Rd. this avoids expensive exp calls. we use two
+   separate tables as well for lower and higher numbers to improve
+   precision, since the number are poorly distributed because we do
+   a lookup with the squared distance for smaller distances, saving
+   another sqrt. */
+
+static void approximate_Rd_rgb(ScatterSettings **ss, float rr, float *rd)
+{
+	float indexf, t, idxf;
+	int index;
+
+	if(rr > (RD_TABLE_RANGE_2*RD_TABLE_RANGE_2)) {
+		rd[0]= Rd_rsquare(ss[0], rr);
+		rd[1]= Rd_rsquare(ss[1], rr);
+		rd[2]= Rd_rsquare(ss[2], rr);
+	}
+	else if(rr > RD_TABLE_RANGE) {
+		rr= sqrt(rr);
+		indexf= rr*(RD_TABLE_SIZE/RD_TABLE_RANGE_2);
+		index= (int)indexf;
+		idxf= (float)index;
+		t= indexf - idxf;
+
+		rd[0]= (ss[0]->tableRd2[index]*(1-t) + ss[0]->tableRd2[index+1]*t);
+		rd[1]= (ss[1]->tableRd2[index]*(1-t) + ss[1]->tableRd2[index+1]*t);
+		rd[2]= (ss[2]->tableRd2[index]*(1-t) + ss[2]->tableRd2[index+1]*t);
+	}
+	else {
+		indexf= rr*(RD_TABLE_SIZE/RD_TABLE_RANGE);
+		index= (int)indexf;
+		idxf= (float)index;
+		t= indexf - idxf;
+
+		rd[0]= (ss[0]->tableRd[index]*(1-t) + ss[0]->tableRd[index+1]*t);
+		rd[1]= (ss[1]->tableRd[index]*(1-t) + ss[1]->tableRd[index+1]*t);
+		rd[2]= (ss[2]->tableRd[index]*(1-t) + ss[2]->tableRd[index+1]*t);
+	}
+}
+
+static void build_Rd_table(ScatterSettings *ss)
+{
+	float r;
+	int i, size = RD_TABLE_SIZE+1;
+
+	ss->tableRd= MEM_mallocN(sizeof(float)*size, "scatterTableRd");
+	ss->tableRd2= MEM_mallocN(sizeof(float)*size, "scatterTableRd");
+
+	for(i= 0; i < size; i++) {
+		r= i*(RD_TABLE_RANGE/RD_TABLE_SIZE);
+		/*if(r < ss->invsigma_t_*ss->invsigma_t_)
+			r= ss->invsigma_t_*ss->invsigma_t_;*/
+		ss->tableRd[i]= Rd(ss, sqrt(r));
+
+		r= i*(RD_TABLE_RANGE_2/RD_TABLE_SIZE);
+		/*if(r < ss->invsigma_t_)
+			r= ss->invsigma_t_;*/
+		ss->tableRd2[i]= Rd(ss, r);
+	}
+}
+
+ScatterSettings *scatter_settings_new(float refl, float radius, float ior, float reflfac, float frontweight, float backweight)
+{
+	ScatterSettings *ss;
+	
+	ss= MEM_callocN(sizeof(ScatterSettings), "ScatterSettings");
+
+	/* see [1] and [3] for these formulas */
+	ss->eta= ior;
+	ss->Fdr= -1.440f/ior*ior + 0.710f/ior + 0.668f + 0.0636f*ior;
+	ss->A= (1.0f + ss->Fdr)/(1.0f - ss->Fdr);
+	ss->ld= radius;
+	ss->ro= MIN2(refl, 0.999f);
+	ss->color= ss->ro*reflfac + (1.0f-reflfac);
+
+	ss->alpha_= compute_reduced_albedo(ss);
+
+	ss->sigma= 1.0f/ss->ld;
+	ss->sigma_t_= ss->sigma/sqrt(3.0f*(1.0f - ss->alpha_));
+	ss->sigma_s_= ss->alpha_*ss->sigma_t_;
+	ss->sigma_a= ss->sigma_t_ - ss->sigma_s_;
+
+	ss->D= 1.0f/(3.0f*ss->sigma_t_);
+
+	ss->zr= 1.0f/ss->sigma_t_;
+	ss->zv= ss->zr + 4.0f*ss->A*ss->D;
+
+	ss->invsigma_t_= 1.0f/ss->sigma_t_;
+
+	ss->frontweight= frontweight;
+	ss->backweight= backweight;
+
+	/* precompute a table of Rd values for quick lookup */
+	build_Rd_table(ss);
+
+	return ss;
+}
+
+void scatter_settings_free(ScatterSettings *ss)
+{
+	MEM_freeN(ss->tableRd);
+	MEM_freeN(ss->tableRd2);
+	MEM_freeN(ss);
+}
+
+/* Hierarchical method as in [2]. */
+
+/* traversal */
+
+#define SUBNODE_INDEX(co, split) \
+	((co[0]>=split[0]) + (co[1]>=split[1])*2 + (co[2]>=split[2])*4)
+	
+static void add_radiance(ScatterTree *tree, float *frontrad, float *backrad, float area, float rr, ScatterResult *result)
+{
+	float rd[3];
+
+#if 0
+	rd[0]= Rd_rsquare(tree->ss[0], rr);
+	rd[1]= Rd_rsquare(tree->ss[1], rr);
+	rd[2]= Rd_rsquare(tree->ss[2], rr);
+#else
+	approximate_Rd_rgb(tree->ss, rr, rd);
+#endif
+
+	rd[0] *= area;
+	rd[1] *= area;
+	rd[2] *= area;
+
+	if(frontrad) {
+		result->radiance[0] += frontrad[0]*rd[0];
+		result->radiance[1] += frontrad[1]*rd[1];
+		result->radiance[2] += frontrad[2]*rd[2];
+	}
+	if(backrad) {
+		result->backradiance[0] += backrad[0]*rd[0];
+		result->backradiance[1] += backrad[1]*rd[1];
+		result->backradiance[2] += backrad[2]*rd[2];
+	}
+
+	result->rdsum[0] += rd[0];
+	result->rdsum[1] += rd[1];
+	result->rdsum[2] += rd[2];
+}
+
+static void traverse_octree(ScatterTree *tree, ScatterNode *node, float *co, int self, ScatterResult *result)
+{
+	float sub[3], dist;
+	int i, index = 0;
+
+	if(node->totpoint > 0) {
+		/* leaf - add radiance from all samples */
+		for(i=0; i<node->totpoint; i++) {
+			ScatterPoint *p= &node->points[i];
+
+			VECSUB(sub, co, p->co);
+			dist= INPR(sub, sub);
+
+			if(p->back)
+				add_radiance(tree, NULL, p->radiance, p->area, dist, result);
+			else
+				add_radiance(tree, p->radiance, NULL, p->area, dist, result);
+		}
+	}
+	else {
+		/* branch */
+		if (self)
+			index = SUBNODE_INDEX(co, node->split);
+
+		for(i=0; i<8; i++) {
+			if(node->child[i]) {
+				ScatterNode *subnode= node->child[i];
+
+				if(self && index == i) {
+					/* always traverse node containing the point */
+					traverse_octree(tree, subnode, co, 1, result);
+				}
+				else {
+					/* decide subnode traversal based on maximum solid angle */
+					VECSUB(sub, co, subnode->co);
+					dist= INPR(sub, sub);
+
+					/* actually area/dist > error, but this avoids division */
+					if(subnode->area>tree->error*dist) {
+						traverse_octree(tree, subnode, co, 0, result);
+					}
+					else {
+						add_radiance(tree, subnode->radiance,
+							subnode->backradiance, subnode->area, dist, result);
+					}
+				}
+			}
+		}
+	}
+}
+
+static void compute_radiance(ScatterTree *tree, float *co, float *radiance)
+{
+	ScatterResult result;
+	float rdsum[3];
+
+	memset(&result, 0, sizeof(result));
+
+	traverse_octree(tree, tree->root, co, 1, &result);
+
+	VecMulf(result.radiance, tree->ss[0]->frontweight);
+	VecMulf(result.backradiance, tree->ss[0]->backweight);
+
+	VECCOPY(radiance, result.backradiance);
+	VECADD(radiance, radiance, result.radiance);
+
+	VECCOPY(rdsum, result.rdsum);
+
+	/* the original paper doesn't do this, but we normalize over the
+	   sampled area and multiply with the reflectance. this is because
+	   our point samples are incomplete, there are no samples on parts
+	   of the mesh not visible from the camera. this can not only make
+	   it darker, but also lead to ugly color shifts */
+	if(rdsum[0] > 0.0f) radiance[0]= tree->ss[0]->color*radiance[0]/rdsum[0];
+	if(rdsum[1] > 0.0f) radiance[1]= tree->ss[1]->color*radiance[1]/rdsum[1];
+	if(rdsum[2] > 0.0f) radiance[2]= tree->ss[2]->color*radiance[2]/rdsum[2];
+}
+
+/* building */
+
+static void sum_leaf_radiance(ScatterTree *tree, ScatterNode *node)
+{
+	ScatterPoint *p;
+	float radiance, totradiance= 0.0f, inv;
+	int i;
+
+	node->co[0]= node->co[1]= node->co[2]= 0.0;
+	node->radiance[0]= node->radiance[1]= node->radiance[2]= 0.0;
+	node->backradiance[0]= node->backradiance[1]= node->backradiance[2]= 0.0;
+
+	/* compute total radiance, radiance weighted average position,
+	   and total area */
+	for(i=0; i<node->totpoint; i++) {
+		p= &node->points[i];
+
+		radiance= p->area*(p->radiance[0] + p->radiance[1] + p->radiance[2]);
+		totradiance += radiance;
+
+		node->co[0] += radiance*p->co[0];
+		node->co[1] += radiance*p->co[1];
+		node->co[2] += radiance*p->co[2];
+
+		if(p->back) {
+			node->backradiance[0] += p->radiance[0]*p->area;
+			node->backradiance[1] += p->radiance[1]*p->area;
+			node->backradiance[2] += p->radiance[2]*p->area;
+		}
+		else {
+			node->radiance[0] += p->radiance[0]*p->area;
+			node->radiance[1] += p->radiance[1]*p->area;
+			node->radiance[2] += p->radiance[2]*p->area;
+		}
+
+		node->area += p->area;
+	}
+
+	if(node->area > 0) {
+		inv= 1.0/node->area;
+		node->radiance[0] *= inv;
+		node->radiance[1] *= inv;
+		node->radiance[2] *= inv;
+		node->backradiance[0] *= inv;
+		node->backradiance[1] *= inv;
+		node->backradiance[2] *= inv;
+	}
+
+	if(totradiance > 0.0f) {
+		inv= 1.0/totradiance;
+		node->co[0] *= inv;
+		node->co[1] *= inv;
+		node->co[2] *= inv;
+	}
+	else {
+		/* make sure that if radiance is 0.0f, we still have these points in
+		   the tree at a good position, they count for rdsum too */
+		for(i=0; i<node->totpoint; i++) {
+			p= &node->points[i];
+
+			node->co[0] += p->co[0];
+			node->co[1] += p->co[1];
+			node->co[2] += p->co[2];
+		}
+
+		node->co[0] /= node->totpoint;
+		node->co[1] /= node->totpoint;
+		node->co[2] /= node->totpoint;
+	}
+}
+
+static void sum_branch_radiance(ScatterTree *tree, ScatterNode *node)
+{
+	ScatterNode *subnode;
+	float radiance, totradiance= 0.0f, inv;
+	int i, totnode;
+
+	node->co[0]= node->co[1]= node->co[2]= 0.0;
+	node->radiance[0]= node->radiance[1]= node->radiance[2]= 0.0;
+	node->backradiance[0]= node->backradiance[1]= node->backradiance[2]= 0.0;
+
+	/* compute total radiance, radiance weighted average position,
+	   and total area */
+	for(i=0; i<8; i++) {
+		if(node->child[i] == NULL)
+			continue;
+
+		subnode= node->child[i];
+
+		radiance= subnode->area*(subnode->radiance[0] + subnode->radiance[1] + subnode->radiance[2] + subnode->backradiance[0] + subnode->backradiance[1] + subnode->backradiance[2]);
+		totradiance += radiance;
+
+		node->co[0] += radiance*subnode->co[0];
+		node->co[1] += radiance*subnode->co[1];
+		node->co[2] += radiance*subnode->co[2];
+
+		node->radiance[0] += subnode->radiance[0]*subnode->area;
+		node->radiance[1] += subnode->radiance[1]*subnode->area;
+		node->radiance[2] += subnode->radiance[2]*subnode->area;
+
+		node->backradiance[0] += subnode->backradiance[0]*subnode->area;
+		node->backradiance[1] += subnode->backradiance[1]*subnode->area;
+		node->backradiance[2] += subnode->backradiance[2]*subnode->area;
+
+		node->area += subnode->area;
+	}
+
+	if(node->area > 0) {
+		inv= 1.0/node->area;
+		node->radiance[0] *= inv;
+		node->radiance[1] *= inv;
+		node->radiance[2] *= inv;
+		node->backradiance[0] *= inv;
+		node->backradiance[1] *= inv;
+		node->backradiance[2] *= inv;
+	}
+
+	if(totradiance > 0.0f) {
+		inv= 1.0/totradiance;
+		node->co[0] *= inv;
+		node->co[1] *= inv;
+		node->co[2] *= inv;
+	}
+	else {
+		/* make sure that if radiance is 0.0f, we still have these points in
+		   the tree at a good position, they count for rdsum too */
+		totnode= 0;
+
+		for(i=0; i<8; i++) {
+			if(node->child[i]) {
+				subnode= node->child[i];
+
+				node->co[0] += subnode->co[0];
+				node->co[1] += subnode->co[1];
+				node->co[2] += subnode->co[2];
+
+				totnode++;
+			}
+		}
+
+		node->co[0] /= totnode;
+		node->co[1] /= totnode;
+		node->co[2] /= totnode;
+	}
+}
+
+static void sum_radiance(ScatterTree *tree, ScatterNode *node)
+{
+	if(node->totpoint > 0) {
+		sum_leaf_radiance(tree, node);
+	}
+	else {
+		int i;
+
+		for(i=0; i<8; i++)
+			if(node->child[i])
+				sum_radiance(tree, node->child[i]);
+
+		sum_branch_radiance(tree, node);
+	}
+}
+
+static void subnode_middle(int i, float *mid, float *subsize, float *submid)
+{
+	int x= i & 1, y= i & 2, z= i & 4;
+
+	submid[0]= mid[0] + ((x)? subsize[0]: -subsize[0]);
+	submid[1]= mid[1] + ((y)? subsize[1]: -subsize[1]);
+	submid[2]= mid[2] + ((z)? subsize[2]: -subsize[2]);
+}
+
+static void create_octree_node(ScatterTree *tree, ScatterNode *node, float *mid, float *size, ScatterPoint **refpoints, int depth)
+{
+	ScatterNode *subnode;
+	ScatterPoint **subrefpoints, **tmppoints= tree->tmppoints;
+	int index, nsize[8], noffset[8], i, subco, usednodes, usedi;
+	float submid[3], subsize[3];
+
+	/* stopping condition */
+	if(node->totpoint <= MAX_OCTREE_NODE_POINTS || depth == MAX_OCTREE_DEPTH) {
+		for(i=0; i<node->totpoint; i++)
+			node->points[i]= *(refpoints[i]);
+
+		return;
+	}
+
+	subsize[0]= size[0]*0.5;
+	subsize[1]= size[1]*0.5;
+	subsize[2]= size[2]*0.5;
+
+	node->split[0]= mid[0];
+	node->split[1]= mid[1];
+	node->split[2]= mid[2];
+
+	memset(nsize, 0, sizeof(nsize));
+	memset(noffset, 0, sizeof(noffset));
+
+	/* count points in subnodes */
+	for(i=0; i<node->totpoint; i++) {
+		index= SUBNODE_INDEX(refpoints[i]->co, node->split);
+		tmppoints[i]= refpoints[i];
+		nsize[index]++;
+	}
+
+	/* here we check if only one subnode is used. if this is the case, we don't
+	   create a new node, but rather call this function again, with different 
+	   size and middle position for the same node. */
+	for(usedi=0, usednodes=0, i=0; i<8; i++) {
+		if(nsize[i]) {
+			usednodes++;
+			usedi = i;
+		}
+		if(i != 0)
+			noffset[i]= noffset[i-1]+nsize[i-1];
+	}
+	
+	if(usednodes<=1) {
+		subnode_middle(usedi, mid, subsize, submid);
+		create_octree_node(tree, node, submid, subsize, refpoints, depth+1);
+		return;
+	}
+
+	/* reorder refpoints by subnode */
+	for(i=0; i<node->totpoint; i++) {
+		index= SUBNODE_INDEX(tmppoints[i]->co, node->split);
+		refpoints[noffset[index]]= tmppoints[i];
+		noffset[index]++;
+	}
+
+	/* create subnodes */
+	for(subco=0, i=0; i<8; subco+=nsize[i], i++) {
+		if(nsize[i] > 0) {
+			subnode= BLI_memarena_alloc(tree->arena, sizeof(ScatterNode));
+			node->child[i]= subnode;
+			subnode->points= node->points + subco;
+			subnode->totpoint= nsize[i];
+			subrefpoints= refpoints + subco;
+
+			subnode_middle(i, mid, subsize, submid);
+
+			create_octree_node(tree, subnode, submid, subsize, subrefpoints,
+				depth+1);
+		}
+		else
+			node->child[i]= NULL;
+	}
+
+	node->points= NULL;
+	node->totpoint= 0;
+}
+
+/* public functions */
+
+ScatterTree *scatter_tree_new(ScatterSettings *ss[3], float scale, float error,
+	float (*co)[3], float (*color)[3], float *area, int totpoint)
+{
+	ScatterTree *tree;
+	ScatterPoint *points, **refpoints;
+	int i;
+
+	/* allocate tree */
+	tree= MEM_callocN(sizeof(ScatterTree), "ScatterTree");
+	tree->scale= scale;
+	tree->error= error;
+	tree->totpoint= totpoint;
+
+	tree->ss[0]= ss[0];
+	tree->ss[1]= ss[1];
+	tree->ss[2]= ss[2];
+
+	points= MEM_callocN(sizeof(ScatterPoint)*totpoint, "ScatterPoints");
+	refpoints= MEM_callocN(sizeof(ScatterPoint*)*totpoint, "ScatterRefPoints");
+
+	tree->points= points;
+	tree->refpoints= refpoints;
+
+	/* build points */
+	INIT_MINMAX(tree->min, tree->max);
+
+	for(i=0; i<totpoint; i++) {
+		VECCOPY(points[i].co, co[i]);
+		VECCOPY(points[i].radiance, color[i]);
+		points[i].area= fabs(area[i])/(tree->scale*tree->scale);
+		points[i].back= (area[i] < 0.0f);
+
+		VecMulf(points[i].co, 1.0f/tree->scale);
+		DO_MINMAX(points[i].co, tree->min, tree->max);
+
+		refpoints[i]= points + i;
+	}
+
+	return tree;
+}
+
+void scatter_tree_build(ScatterTree *tree)
+{
+	ScatterPoint *newpoints, **tmppoints;
+	float mid[3], size[3];
+	int totpoint= tree->totpoint;
+
+	newpoints= MEM_callocN(sizeof(ScatterPoint)*totpoint, "ScatterPoints");
+	tmppoints= MEM_callocN(sizeof(ScatterPoint*)*totpoint, "ScatterTmpPoints");
+	tree->tmppoints= tmppoints;
+
+	tree->arena= BLI_memarena_new(0x8000 * sizeof(ScatterNode));
+	BLI_memarena_use_calloc(tree->arena);
+
+	/* build tree */
+	tree->root= BLI_memarena_alloc(tree->arena, sizeof(ScatterNode));
+	tree->root->points= newpoints;
+	tree->root->totpoint= totpoint;
+
+	mid[0]= (tree->min[0]+tree->max[0])*0.5;
+	mid[1]= (tree->min[1]+tree->max[1])*0.5;
+	mid[2]= (tree->min[2]+tree->max[2])*0.5;
+
+	size[0]= (tree->max[0]-tree->min[0])*0.5;
+	size[1]= (tree->max[1]-tree->min[1])*0.5;
+	size[2]= (tree->max[2]-tree->min[2])*0.5;
+
+	create_octree_node(tree, tree->root, mid, size, tree->refpoints, 0);
+
+	MEM_freeN(tree->points);
+	MEM_freeN(tree->refpoints);
+	MEM_freeN(tree->tmppoints);
+	tree->refpoints= NULL;
+	tree->tmppoints= NULL;
+	tree->points= newpoints;
+	
+	/* sum radiance at nodes */
+	sum_radiance(tree, tree->root);
+}
+
+void scatter_tree_sample(ScatterTree *tree, float *co, float *color)
+{
+	float sco[3];
+
+	VECCOPY(sco, co);
+	VecMulf(sco, 1.0f/tree->scale);
+
+	compute_radiance(tree, sco, color);
+}
+
+void scatter_tree_free(ScatterTree *tree)
+{
+	if (tree->arena) BLI_memarena_free(tree->arena);
+	if (tree->points) MEM_freeN(tree->points);
+	if (tree->refpoints) MEM_freeN(tree->refpoints);
+		
+	MEM_freeN(tree);
+}
+
+/* Internal Renderer API */
+
+/* sss tree building */
+
+typedef struct SSSData {
+	ScatterTree *tree;
+	ScatterSettings *ss[3];
+} SSSData;
+
+typedef struct SSSPoints {
+	struct SSSPoints *next, *prev;
+
+	float (*co)[3];
+	float (*color)[3];
+	float *area;
+	int totpoint;
+} SSSPoints;
+
+static void sss_create_tree_mat(Render *re, Material *mat)
+{
+	SSSPoints *p;
+	ListBase layers, points;
+	float (*co)[3] = NULL, (*color)[3] = NULL, *area = NULL;
+	int totpoint = 0, osa;
+
+	if(re->test_break())
+		return;
+	
+	points.first= points.last= NULL;
+
+	/* do SSS preprocessing render */
+	layers= re->r.layers;
+	osa= re->osa;
+
+	re->r.layers.first= re->r.layers.last= NULL;
+	re->osa= 0;
+	re->sss_points= &points;
+	re->sss_mat= mat;
+
+	RE_TileProcessor(re, 0, 1);
+
+	re->sss_mat= NULL;
+	re->sss_points= NULL;
+	re->r.layers= layers;
+	re->osa= osa;
+
+	/* no points? no tree */
+	if(!points.first)
+		return;
+
+	/* merge points together into a single buffer */
+	if(!re->test_break()) {
+		for(totpoint=0, p=points.first; p; p=p->next)
+			totpoint += p->totpoint;
+		
+		co= MEM_mallocN(sizeof(*co)*totpoint, "SSSCo");
+		color= MEM_mallocN(sizeof(*color)*totpoint, "SSSColor");
+		area= MEM_mallocN(sizeof(*area)*totpoint, "SSSArea");
+
+		for(totpoint=0, p=points.first; p; p=p->next) {
+			memcpy(co+totpoint, p->co, sizeof(*co)*p->totpoint);
+			memcpy(color+totpoint, p->color, sizeof(*color)*p->totpoint);
+			memcpy(area+totpoint, p->area, sizeof(*area)*p->totpoint);
+			totpoint += p->totpoint;
+		}
+	}
+
+	/* free points */
+	for(p=points.first; p; p=p->next) {
+		MEM_freeN(p->co);
+		MEM_freeN(p->color);
+		MEM_freeN(p->area);
+	}
+	BLI_freelistN(&points);
+
+	/* build tree */
+	if(!re->test_break()) {
+		SSSData *sss= MEM_callocN(sizeof(*sss), "SSSData");
+		float ior= mat->sss_ior, cfac= mat->sss_colfac;
+		float *col= mat->sss_col, *radius= mat->sss_radius;
+		float fw= mat->sss_front, bw= mat->sss_back;
+
+		sss->ss[0]= scatter_settings_new(col[0], radius[0], ior, cfac, fw, bw);
+		sss->ss[1]= scatter_settings_new(col[1], radius[1], ior, cfac, fw, bw);
+		sss->ss[2]= scatter_settings_new(col[2], radius[2], ior, cfac, fw, bw);
+		sss->tree= scatter_tree_new(sss->ss, mat->sss_scale, mat->sss_error,
+			co, color, area, totpoint);
+
+		MEM_freeN(co);
+		MEM_freeN(color);
+		MEM_freeN(area);
+
+		scatter_tree_build(sss->tree);
+
+		BLI_ghash_insert(re->sss_hash, mat, sss);
+	}
+	else {
+		if (co) MEM_freeN(co);
+		if (color) MEM_freeN(color);
+		if (area) MEM_freeN(area);
+	}
+}
+
+void sss_add_points(Render *re, float (*co)[3], float (*color)[3], float *area, int totpoint)
+{
+	SSSPoints *p;
+	
+	if(totpoint > 0) {
+		p= MEM_callocN(sizeof(SSSPoints), "SSSPoints");
+
+		p->co= co;
+		p->color= color;
+		p->area= area;
+		p->totpoint= totpoint;
+
+		BLI_lock_thread(LOCK_CUSTOM1);
+		BLI_addtail(re->sss_points, p);
+		BLI_unlock_thread(LOCK_CUSTOM1);
+	}
+}
+
+static void sss_free_tree(SSSData *sss)
+{
+	scatter_tree_free(sss->tree);
+	scatter_settings_free(sss->ss[0]);
+	scatter_settings_free(sss->ss[1]);
+	scatter_settings_free(sss->ss[2]);
+	MEM_freeN(sss);
+}
+
+/* public functions */
+
+void make_sss_tree(Render *re)
+{
+	Material *mat;
+	
+	re->sss_hash= BLI_ghash_new(BLI_ghashutil_ptrhash, BLI_ghashutil_ptrcmp);
+
+	re->i.infostr= "SSS preprocessing";
+	re->stats_draw(&re->i);
+	
+	for(mat= G.main->mat.first; mat; mat= mat->id.next)
+		if(mat->id.us && (mat->sss_flag & MA_DIFF_SSS))
+			sss_create_tree_mat(re, mat);
+}
+
+void free_sss(Render *re)
+{
+	if(re->sss_hash) {
+		GHashIterator *it= BLI_ghashIterator_new(re->sss_hash);
+
+		while(!BLI_ghashIterator_isDone(it)) {
+			sss_free_tree(BLI_ghashIterator_getValue(it));
+			BLI_ghashIterator_step(it);
+		}
+
+		BLI_ghashIterator_free(it);
+		BLI_ghash_free(re->sss_hash, NULL, NULL);
+		re->sss_hash= NULL;
+	}
+}
+
+int sample_sss(Render *re, Material *mat, float *co, float *color)
+{
+	if(re->sss_hash) {
+		SSSData *sss= BLI_ghash_lookup(re->sss_hash, mat);
+
+		if(sss) {
+			scatter_tree_sample(sss->tree, co, color);
+			return 1;
+		}
+		else {
+			color[0]= 0.0f;
+			color[1]= 0.0f;
+			color[2]= 0.0f;
+		}
+	}
+
+	return 0;
+}
+
+int has_sss_tree(struct Render *re, struct Material *mat)
+{
+	return (re->sss_hash && BLI_ghash_lookup(re->sss_hash, mat));
+}
+
diff --git a/source/blender/render/intern/source/zbuf.c b/source/blender/render/intern/source/zbuf.c
index ee24b15029c..f919aa11b98 100644
--- a/source/blender/render/intern/source/zbuf.c
+++ b/source/blender/render/intern/source/zbuf.c
@@ -51,6 +51,7 @@
 #include "DNA_meshdata_types.h"
 
 #include "BKE_global.h"
+#include "BKE_material.h"
 #include "BKE_utildefines.h"
 
 #include "radio_types.h"
@@ -67,6 +68,7 @@
 #include "rendercore.h"
 #include "shadbuf.h"
 #include "shading.h"
+#include "sss.h"
 
 /* own includes */
 #include "zbuf.h"
@@ -2012,6 +2014,164 @@ void zbuffer_shadow(Render *re, LampRen *lar, int *rectz, int size, float jitx,
 	zbuf_free_span(&zspan);
 }
 
+static void zbuffill_sss(ZSpan *zspan, int zvlnr, float *v1, float *v2, float *v3, float *v4)
+{
+	double zxd, zyd, zy0, z;
+	float x0, y0, x1, y1, x2, y2, z0, z1, z2, xx1, *span1, *span2;
+	int x, y, sn1, sn2, rectx= zspan->rectx, my0, my2;
+	
+	/* init */
+	zbuf_init_span(zspan);
+	
+	/* set spans */
+	zbuf_add_to_span(zspan, v1, v2);
+	zbuf_add_to_span(zspan, v2, v3);
+	if(v4) {
+		zbuf_add_to_span(zspan, v3, v4);
+		zbuf_add_to_span(zspan, v4, v1);
+	}
+	else 
+		zbuf_add_to_span(zspan, v3, v1);
+	
+	/* clipped */
+	if(zspan->minp2==NULL || zspan->maxp2==NULL) return;
+	
+	if(zspan->miny1 < zspan->miny2) my0= zspan->miny2; else my0= zspan->miny1;
+	if(zspan->maxy1 > zspan->maxy2) my2= zspan->maxy2; else my2= zspan->maxy1;
+	
+	if(my2<my0) return;
+	
+	/* ZBUF DX DY, in floats still */
+	x1= v1[0]- v2[0];
+	x2= v2[0]- v3[0];
+	y1= v1[1]- v2[1];
+	y2= v2[1]- v3[1];
+	z1= v1[2]- v2[2];
+	z2= v2[2]- v3[2];
+	
+	x0= y1*z2-z1*y2;
+	y0= z1*x2-x1*z2;
+	z0= x1*y2-y1*x2;
+	
+	if(z0==0.0f) return;
+	
+	xx1= (x0*v1[0] + y0*v1[1])/z0 + v1[2];
+	zxd= -(double)x0/(double)z0;
+	zyd= -(double)y0/(double)z0;
+	zy0= ((double)my2)*zyd + (double)xx1;
+	
+	/* correct span */
+	sn1= (my0 + my2)/2;
+	if(zspan->span1[sn1] < zspan->span2[sn1]) {
+		span1= zspan->span1+my2;
+		span2= zspan->span2+my2;
+	}
+	else {
+		span1= zspan->span2+my2;
+		span2= zspan->span1+my2;
+	}
+	
+	for(y=my2; y>=my0; y--, span1--, span2--) {
+		sn1= floor(*span1);
+		sn2= floor(*span2);
+		sn1++; 
+		
+		if(sn2>=rectx) sn2= rectx-1;
+		if(sn1<0) sn1= 0;
+		
+		z= (double)sn1*zxd + zy0;
+		
+		for(x= sn1; x<=sn2; x++, z+=zxd)
+			zspan->sss_func(zspan->sss_handle, zvlnr, x, y, z);
+		
+		zy0 -= zyd;
+	}
+}
+
+void zbuffer_sss(RenderPart *pa, unsigned int lay, void *handle, void (*func)(void *, int, int, int, int))
+{
+	ZSpan zspan;
+	VlakRen *vlr= NULL;
+	VertRen *v1, *v2, *v3, *v4;
+	Material *ma=0, *sss_ma= R.sss_mat;
+	short nofill=0, env=0, wire=0;
+	char *clipflag= pa->clipflag;
+	int v, zvlnr;
+	
+	zbuf_alloc_span(&zspan, pa->rectx, pa->recty);
+
+	zspan.sss_handle= handle;
+	zspan.sss_func= func;
+	
+	/* needed for transform from hoco to zbuffer co */
+	zspan.zmulx=  ((float)R.winx)/2.0;
+	zspan.zmuly=  ((float)R.winy)/2.0;
+	
+	/* -0.5f to center the sample position */
+	zspan.zofsx= -pa->disprect.xmin - 0.5f;
+	zspan.zofsy= -pa->disprect.ymin - 0.5f;
+	
+	/* filling methods */
+	zspan.zbuffunc= zbuffill_sss;
+
+	/* fill front and back zbuffer */
+	if(pa->rectz) {
+		fillrect(pa->rectp, pa->rectx, pa->recty, 0); 
+		fillrect(pa->rectz, pa->rectx, pa->recty, 0x7FFFFFFF);
+	}
+	if(pa->rectbackz) {
+		fillrect(pa->rectbackp, pa->rectx, pa->recty, 0); 
+		fillrect(pa->rectbackz, pa->rectx, pa->recty, -0x7FFFFFFF);
+	}
+
+	for(v=0; v<R.totvlak; v++) {
+		if((v & 255)==0) vlr= R.vlaknodes[v>>8].vlak;
+		else vlr++;
+		
+		if((vlr->flag & R_VISIBLE) && material_in_material(vlr->mat, sss_ma)) {
+			/* three cases, visible for render, only z values and nothing */
+			if(vlr->lay & lay) {
+				if(vlr->mat!=ma) {
+					ma= vlr->mat;
+					nofill= ma->mode & MA_ONLYCAST;
+					env= (ma->mode & MA_ENV);
+					wire= (ma->mode & MA_WIRE);
+				}
+			}
+			else {
+				nofill= 1;
+				ma= NULL;	/* otherwise nofill can hang */
+			}
+			
+			if(nofill==0 && wire==0 && env==0) {
+				unsigned short partclip;
+				
+				v1= vlr->v1;
+				v2= vlr->v2;
+				v3= vlr->v3;
+				v4= vlr->v4;
+				
+				/* partclipping doesn't need viewplane clipping */
+				partclip= clipflag[v1->index] & clipflag[v2->index] & clipflag[v3->index];
+				if(v4)
+					partclip &= clipflag[v4->index];
+				
+
+				if(partclip==0) {
+					zvlnr= v+1;
+					zbufclip(&zspan, zvlnr, v1->ho, v2->ho, v3->ho, v1->clip, v2->clip, v3->clip);
+					if(v4) {
+						zvlnr+= RE_QUAD_OFFS;
+						zbufclip(&zspan, zvlnr, v1->ho, v3->ho, v4->ho, v1->clip, v3->clip, v4->clip);
+					}
+				}
+			}
+		}
+	}
+	
+	zbuf_free_span(&zspan);
+}
+
 /* ******************** VECBLUR ACCUM BUF ************************* */
 
 typedef struct DrawBufPixel {
diff --git a/source/blender/src/buttons_shading.c b/source/blender/src/buttons_shading.c
index bd8dae311e8..b290857c23d 100644
--- a/source/blender/src/buttons_shading.c
+++ b/source/blender/src/buttons_shading.c
@@ -3269,6 +3269,118 @@ static void material_panel_tramir_yafray(Material *ma)
 
 }
 
+static void material_sss_preset_cb(void *material_v, void *unused_v)
+{
+	static const float presets[11][7] = {
+	{0.909578, 0.905931, 0.665691, 6.961082,   6.400181, 1.895899, 1.300000},
+	{0.429632, 0.210025, 0.167767, 11.605177,  3.883766, 1.754386, 1.300000},
+	{0.439300, 0.216000, 0.141027, 9.435642,   3.347647, 1.790287, 1.300000},
+	{0.986552, 0.942637, 0.827285, 15.027623,  4.663968, 2.541380, 1.300000},
+	{0.221636, 0.007505, 0.002154, 4.761743,   0.574827, 0.394116, 1.300000},
+	{0.925008, 0.905025, 0.884275, 8.509412,   5.566180, 3.951266, 1.500000},
+	{0.855344, 0.740311, 0.291994, 14.266395,  7.227615, 2.036157, 1.300000},
+	{0.889319, 0.888034, 0.795811, 18.424364, 10.443473, 3.501882, 1.300000},
+	{0.573652, 0.312750, 0.174289, 3.673294,   1.366534, 0.682693, 1.300000},
+	{0.748679, 0.570766, 0.467133, 4.821475,   1.693699, 1.089971, 1.300000},
+	{0.947235, 0.931028, 0.851872, 10.898815,  6.575351, 2.508417, 1.300000}};
+
+	Material *ma= (Material*)material_v;
+
+	if(ma->sss_preset==0) return;
+
+	ma->sss_col[0]= presets[ma->sss_preset][0];
+	ma->sss_col[1]= presets[ma->sss_preset][1];
+	ma->sss_col[2]= presets[ma->sss_preset][2];
+	ma->sss_radius[0]= presets[ma->sss_preset][3];
+	ma->sss_radius[1]= presets[ma->sss_preset][4];
+	ma->sss_radius[2]= presets[ma->sss_preset][5];
+	ma->sss_ior= presets[ma->sss_preset][6];
+}
+
+static void material_sss_custom_set_cb(void *material_v, void *unused_v)
+{
+	Material *ma= (Material*)material_v;
+
+	ma->sss_preset= 0;
+	allqueue(REDRAWNODE, 0);
+}
+
+static void material_panel_sss(Material *ma)
+{
+	uiBlock *block;
+	uiBut *bt;
+	
+	block= uiNewBlock(&curarea->uiblocks, "material_panel_sss", UI_EMBOSS, UI_HELV, curarea->win);
+	uiNewPanelTabbed("Shaders", "Material");
+	if(uiNewPanel(curarea, block, "SSS", "Material", 640, 0, 318, 204)==0) return;
+
+	uiSetButLock(ma->id.lib!=NULL, ERROR_LIBDATA_MESSAGE);
+	
+	uiDefButBitS(block, TOG, MA_DIFF_SSS, B_MATPRV,"Subsurface Scattering",10,180,180,20, &(ma->sss_flag), 0, 0, 0, 0, "Enables diffuse subsurface scattering");
+
+	bt=uiDefButS(block, MENU, B_MATPRV, "Apple %x1|Chicken %x2|Cream %x3|Ketchup %x4|Marble %x5|Potato %x6|Skim Milk %x7|Skin 1 %x8|Skin 2 %x9|Whole Milk %x10|Custom %x0",
+		200,180,110,20, &ma->sss_preset, 0, 0, 0, 0, "");
+	uiButSetFunc(bt, material_sss_preset_cb, ma, NULL);
+
+	uiBlockBeginAlign(block);
+	uiDefButF(block, NUM, B_MATPRV, "Scale:", 10,150,145,20,
+			  &ma->sss_scale, 0.001, 1000, 1, 3, "Object scale");
+	bt=uiDefButF(block, NUM, B_MATPRV, "Radius R", 10,130,145,20,
+			  &ma->sss_radius[0], 0.00001, 10000, 0, 0,
+			  "Mean red scattering path length");
+	uiButSetFunc(bt, material_sss_custom_set_cb, ma, NULL);
+	bt=uiDefButF(block, NUM, B_MATPRV, "Radius G", 10,110,145,20,
+			  &ma->sss_radius[1], 0.00001, 10000, 0, 0,
+			  "Mean green scattering path length");
+	uiButSetFunc(bt, material_sss_custom_set_cb, ma, NULL);
+	bt=uiDefButF(block, NUM, B_MATPRV, "Radius B", 10,90,145,20,
+			  &ma->sss_radius[2], 0.00001, 10000, 0, 0,
+			  "Mean blue scattering path length");
+	uiButSetFunc(bt, material_sss_custom_set_cb, ma, NULL);
+	uiBlockEndAlign(block);
+
+	bt=uiDefButF(block, NUM, B_MATPRV, "IOR:", 10,60,145,20,
+			  &ma->sss_ior, 0.1, 2, 1, 3, "Index of refraction");
+	uiButSetFunc(bt, material_sss_custom_set_cb, ma, NULL);
+
+	uiBlockBeginAlign(block);
+	uiDefButF(block, NUM, B_MATPRV, "Error:", 10,30,145,20,
+			  &ma->sss_error, 0.0001, 10, 1, 3, "Error");
+	uiBlockEndAlign(block);
+
+	uiBlockBeginAlign(block);
+	bt=uiDefButF(block, NUMSLI, B_MATPRV, "R ", 165,150,145,20,
+			  &ma->sss_col[0], 0.0, 1.0, 0, 0,
+			  "Red scattering color");
+	uiButSetFunc(bt, material_sss_custom_set_cb, ma, NULL);
+	bt=uiDefButF(block, NUMSLI, B_MATPRV, "G ", 165,130,145,20,
+			  &ma->sss_col[1], 0.0, 1.0, 0, 0,
+			  "Green scattering color");
+	uiButSetFunc(bt, material_sss_custom_set_cb, ma, NULL);
+	bt=uiDefButF(block, NUMSLI, B_MATPRV, "B ", 165,110,145,20,
+			  &ma->sss_col[2], 0.0, 1.0, 0, 0,
+			  "Blue scattering color");
+	uiButSetFunc(bt, material_sss_custom_set_cb, ma, NULL);
+	uiBlockEndAlign(block);
+
+	uiBlockBeginAlign(block);
+	uiDefButF(block, NUMSLI, B_MATPRV, "Col ", 165,80,145,20,
+			  &ma->sss_colfac, 0.0, 1.0, 0, 0,
+			  "Blend factor for SSS colors");
+	uiDefButF(block, NUMSLI, B_MATPRV, "Tex ", 165,60,145,20,
+			  &ma->sss_texfac, 0.0, 1.0, 0, 0,
+			  "Texture scattering factor");
+	uiBlockEndAlign(block);
+
+	uiBlockBeginAlign(block);
+	uiDefButF(block, NUMSLI, B_MATPRV, "Front ", 165,30,145,20,
+			  &ma->sss_front, 0.0, 2.0, 0, 0,
+			  "Front scattering weight");
+	uiDefButF(block, NUMSLI, B_MATPRV, "Back ", 165,10,145,20,
+			  &ma->sss_back, 0.0, 10.0, 0, 0,
+			  "Back scattering weight");
+	uiBlockEndAlign(block);
+}
 
 static void material_panel_shading(Material *ma)
 {
@@ -3731,7 +3843,8 @@ void material_panels()
 				}
 				material_panel_tramir_yafray(ma);
 			}
-			
+
+			material_panel_sss(ma);
 			material_panel_texture(ma);
 			
 			mtex= ma->mtex[ ma->texact ];
diff --git a/source/blender/src/previewrender.c b/source/blender/src/previewrender.c
index d54336807f9..1526de1b74d 100644
--- a/source/blender/src/previewrender.c
+++ b/source/blender/src/previewrender.c
@@ -467,7 +467,7 @@ void BIF_previewrender(struct ID *id, struct RenderInfo *ri, struct ScrArea *are
 	/* entire cycle for render engine */
 	RE_SetCamera(re, sce->camera);
 	RE_Database_FromScene(re, sce, 1);
-	RE_TileProcessor(re, ri->curtile);	// actual render engine
+	RE_TileProcessor(re, ri->curtile, 0);	// actual render engine
 	RE_Database_Free(re);
 	
 	/* handle results */
@@ -842,7 +842,7 @@ void BIF_view3d_previewrender(ScrArea *sa)
 		/* OK, can we enter render code? */
 		if(ri->status==(PR_DISPRECT|PR_DBASE|PR_PROJECTED|PR_ROTATED)) {
 			//printf("curtile %d tottile %d\n", ri->curtile, ri->tottile);
-			RE_TileProcessor(ri->re, ri->curtile);
+			RE_TileProcessor(ri->re, ri->curtile, 0);
 	
 			if(ri->rect==NULL)
 				ri->rect= MEM_mallocN(sizeof(int)*ri->pr_rectx*ri->pr_recty, "preview view3d rect");