Subsurface scattering:
[blender.git] / source / blender / render / intern / source / sss.c
1 /* 
2  * $Id$
3  *
4  * ***** BEGIN GPL LICENSE BLOCK *****
5  *
6  * This program is free software; you can redistribute it and/or
7  * modify it under the terms of the GNU General Public License
8  * as published by the Free Software Foundation; either version 2
9  * of the License, or (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, write to the Free Software Foundation,
18  * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
19  *
20  * The Original Code is Copyright (C) 2007 Blender Foundation.
21  * All rights reserved.
22  *
23  * The Original Code is: all of this file.
24  *
25  * Contributor(s): none yet.
26  *
27  * ***** END GPL LICENSE BLOCK *****
28  */
29
30 /* Possible Improvements:
31    - add fresnel terms
32    - adapt Rd table to scale, now with small scale there are a lot of misses?
33    - possible interesting method: perform sss on all samples in the tree,
34      and then use those values interpolated somehow later. can also do this
35          filtering on demand for speed. since we are doing things in screen
36          space now there is an exact correspondence
37    - avoid duplicate shading (filtering points in advance, irradiance cache
38      like lookup?)
39    - lower resolution samples
40 */
41
42 #include <math.h>
43 #include <string.h>
44 #include <stdio.h>
45 #include <string.h>
46
47 /* external modules: */
48 #include "MEM_guardedalloc.h"
49
50 #include "BLI_arithb.h"
51 #include "BLI_blenlib.h"
52 #include "BLI_ghash.h"
53 #include "BLI_memarena.h"
54 #include "PIL_time.h"
55
56 #include "DNA_material_types.h"
57
58 #include "BKE_global.h"
59 #include "BKE_main.h"
60 #include "BKE_material.h"
61 #include "BKE_node.h"
62 #include "BKE_utildefines.h"
63
64 /* this module */
65 #include "render_types.h"
66 #include "rendercore.h"
67 #include "renderdatabase.h" 
68 #include "shading.h"
69 #include "sss.h"
70 #include "zbuf.h"
71
72 extern Render R; // meh
73
74 /* Generic Multiple Scattering API */
75
76 /* Relevant papers:
77         [1] A Practical Model for Subsurface Light Transport
78         [2] A Rapid Hierarchical Rendering Technique for Translucent Materials
79         [3] Efficient Rendering of Local Subsurface Scattering
80         [4] Implementing a skin BSSRDF (or several...)
81 */
82
83 /* Defines */
84
85 #define RD_TABLE_RANGE          100.0f
86 #define RD_TABLE_RANGE_2        10000.0f
87 #define RD_TABLE_SIZE           10000
88
89 #define MAX_OCTREE_NODE_POINTS  8
90 #define MAX_OCTREE_DEPTH                15
91
92 /* Struct Definitions */
93
94 struct ScatterSettings {
95         float eta;              /* index of refraction */
96         float sigma_a;  /* absorption coefficient */
97         float sigma_s_; /* reduced scattering coefficient */
98         float sigma_t_; /* reduced extinction coefficient */
99         float sigma;    /* effective extinction coefficient */
100         float Fdr;              /* diffuse fresnel reflectance */
101         float D;                /* diffusion constant */
102         float A;
103         float alpha_;   /* reduced albedo */
104         float zr;               /* distance of virtual lightsource above surface */
105         float zv;               /* distance of virtual lightsource below surface */
106         float ld;               /* mean free path */
107         float ro;               /* diffuse reflectance */
108         float color;
109         float invsigma_t_;
110         float frontweight;
111         float backweight;
112
113         float *tableRd;  /* lookup table to avoid computing Rd */
114         float *tableRd2; /* lookup table to avoid computing Rd for bigger values */
115 };
116
117 typedef struct ScatterPoint {
118         float co[3];
119         float radiance[3];
120         float area;
121         int back;
122 } ScatterPoint;
123
124 typedef struct ScatterNode {
125         float co[3];
126         float radiance[3];
127         float backradiance[3];
128         float area;
129         int back;
130
131         int totpoint;
132         ScatterPoint *points;
133
134         float split[3];
135         struct ScatterNode *child[8];
136 } ScatterNode;
137
138 struct ScatterTree {
139         MemArena *arena;
140
141         ScatterSettings *ss[3];
142         float error, scale;
143
144         ScatterNode *root;
145         ScatterPoint *points;
146         ScatterPoint **refpoints;
147         ScatterPoint **tmppoints;
148         int totpoint;
149         float min[3], max[3];
150 };
151
152 typedef struct ScatterResult {
153         float radiance[3];
154         float backradiance[3];
155         float rdsum[3];
156 } ScatterResult;
157
158 /* Functions for BSSRDF reparametrization in to more intuitive parameters,
159    see [2] section 4 for more info. */
160
161 static float f_Rd(float alpha_, float A, float ro)
162 {
163         float sq;
164
165         sq= sqrt(3.0f*(1.0f - alpha_));
166         return (alpha_/2.0f)*(1.0f + exp((-4.0f/3.0f)*A*sq))*exp(-sq) - ro;
167 }
168
169 static float compute_reduced_albedo(ScatterSettings *ss)
170 {
171         const float tolerance= 1e-8;
172         const int max_iteration_count= 20;
173         float d, fsub, xn_1= 0.0f , xn= 1.0f, fxn, fxn_1;
174         int i;
175
176         /* use secant method to compute reduced albedo using Rd function inverse
177            with a given reflectance */
178         fxn= f_Rd(xn, ss->A, ss->ro);
179         fxn_1= f_Rd(xn_1, ss->A, ss->ro);
180
181         for(i= 0; i < max_iteration_count; i++) {
182                 fsub= (fxn - fxn_1);
183                 if(fabs(fsub) < tolerance)
184                         break;
185                 d= ((xn - xn_1)/fsub)*fxn;
186                 if(fabs(d) < tolerance)
187                         break;
188
189                 xn_1= xn;
190                 fxn_1= fxn;
191                 xn= xn - d;
192
193                 if(xn > 1.0f) xn= 1.0f;
194                 if(xn_1 > 1.0f) xn_1= 1.0f;
195                 
196                 fxn= f_Rd(xn, ss->A, ss->ro);
197     }
198
199         /* avoid division by zero later */
200         if(xn <= 0.0f)
201                 xn= 0.00001f;
202
203     return xn;
204 }
205
206 /* Exponential falloff functions */
207
208 static float Rd_rsquare(ScatterSettings *ss, float rr)
209 {
210         float sr, sv, Rdr, Rdv;
211
212         sr= sqrt(rr + ss->zr*ss->zr);
213         sv= sqrt(rr + ss->zv*ss->zv);
214
215         Rdr= ss->zr*(1.0f + ss->sigma*sr)*exp(-ss->sigma*sr)/(sr*sr*sr);
216         Rdv= ss->zv*(1.0f + ss->sigma*sv)*exp(-ss->sigma*sv)/(sv*sv*sv);
217
218         return /*ss->alpha_*/(1.0f/(4.0f*M_PI))*(Rdr + Rdv);
219 }
220
221 static float Rd(ScatterSettings *ss, float r)
222 {
223         return Rd_rsquare(ss, r*r);
224 }
225
226 /* table lookups for Rd. this avoids expensive exp calls. we use two
227    separate tables as well for lower and higher numbers to improve
228    precision, since the number are poorly distributed because we do
229    a lookup with the squared distance for smaller distances, saving
230    another sqrt. */
231
232 static void approximate_Rd_rgb(ScatterSettings **ss, float rr, float *rd)
233 {
234         float indexf, t, idxf;
235         int index;
236
237         if(rr > (RD_TABLE_RANGE_2*RD_TABLE_RANGE_2)) {
238                 rd[0]= Rd_rsquare(ss[0], rr);
239                 rd[1]= Rd_rsquare(ss[1], rr);
240                 rd[2]= Rd_rsquare(ss[2], rr);
241         }
242         else if(rr > RD_TABLE_RANGE) {
243                 rr= sqrt(rr);
244                 indexf= rr*(RD_TABLE_SIZE/RD_TABLE_RANGE_2);
245                 index= (int)indexf;
246                 idxf= (float)index;
247                 t= indexf - idxf;
248
249                 rd[0]= (ss[0]->tableRd2[index]*(1-t) + ss[0]->tableRd2[index+1]*t);
250                 rd[1]= (ss[1]->tableRd2[index]*(1-t) + ss[1]->tableRd2[index+1]*t);
251                 rd[2]= (ss[2]->tableRd2[index]*(1-t) + ss[2]->tableRd2[index+1]*t);
252         }
253         else {
254                 indexf= rr*(RD_TABLE_SIZE/RD_TABLE_RANGE);
255                 index= (int)indexf;
256                 idxf= (float)index;
257                 t= indexf - idxf;
258
259                 rd[0]= (ss[0]->tableRd[index]*(1-t) + ss[0]->tableRd[index+1]*t);
260                 rd[1]= (ss[1]->tableRd[index]*(1-t) + ss[1]->tableRd[index+1]*t);
261                 rd[2]= (ss[2]->tableRd[index]*(1-t) + ss[2]->tableRd[index+1]*t);
262         }
263 }
264
265 static void build_Rd_table(ScatterSettings *ss)
266 {
267         float r;
268         int i, size = RD_TABLE_SIZE+1;
269
270         ss->tableRd= MEM_mallocN(sizeof(float)*size, "scatterTableRd");
271         ss->tableRd2= MEM_mallocN(sizeof(float)*size, "scatterTableRd");
272
273         for(i= 0; i < size; i++) {
274                 r= i*(RD_TABLE_RANGE/RD_TABLE_SIZE);
275                 /*if(r < ss->invsigma_t_*ss->invsigma_t_)
276                         r= ss->invsigma_t_*ss->invsigma_t_;*/
277                 ss->tableRd[i]= Rd(ss, sqrt(r));
278
279                 r= i*(RD_TABLE_RANGE_2/RD_TABLE_SIZE);
280                 /*if(r < ss->invsigma_t_)
281                         r= ss->invsigma_t_;*/
282                 ss->tableRd2[i]= Rd(ss, r);
283         }
284 }
285
286 ScatterSettings *scatter_settings_new(float refl, float radius, float ior, float reflfac, float frontweight, float backweight)
287 {
288         ScatterSettings *ss;
289         
290         ss= MEM_callocN(sizeof(ScatterSettings), "ScatterSettings");
291
292         /* see [1] and [3] for these formulas */
293         ss->eta= ior;
294         ss->Fdr= -1.440f/ior*ior + 0.710f/ior + 0.668f + 0.0636f*ior;
295         ss->A= (1.0f + ss->Fdr)/(1.0f - ss->Fdr);
296         ss->ld= radius;
297         ss->ro= MIN2(refl, 0.999f);
298         ss->color= ss->ro*reflfac + (1.0f-reflfac);
299
300         ss->alpha_= compute_reduced_albedo(ss);
301
302         ss->sigma= 1.0f/ss->ld;
303         ss->sigma_t_= ss->sigma/sqrt(3.0f*(1.0f - ss->alpha_));
304         ss->sigma_s_= ss->alpha_*ss->sigma_t_;
305         ss->sigma_a= ss->sigma_t_ - ss->sigma_s_;
306
307         ss->D= 1.0f/(3.0f*ss->sigma_t_);
308
309         ss->zr= 1.0f/ss->sigma_t_;
310         ss->zv= ss->zr + 4.0f*ss->A*ss->D;
311
312         ss->invsigma_t_= 1.0f/ss->sigma_t_;
313
314         ss->frontweight= frontweight;
315         ss->backweight= backweight;
316
317         /* precompute a table of Rd values for quick lookup */
318         build_Rd_table(ss);
319
320         return ss;
321 }
322
323 void scatter_settings_free(ScatterSettings *ss)
324 {
325         MEM_freeN(ss->tableRd);
326         MEM_freeN(ss->tableRd2);
327         MEM_freeN(ss);
328 }
329
330 /* Hierarchical method as in [2]. */
331
332 /* traversal */
333
334 #define SUBNODE_INDEX(co, split) \
335         ((co[0]>=split[0]) + (co[1]>=split[1])*2 + (co[2]>=split[2])*4)
336         
337 static void add_radiance(ScatterTree *tree, float *frontrad, float *backrad, float area, float rr, ScatterResult *result)
338 {
339         float rd[3];
340
341 #if 0
342         rd[0]= Rd_rsquare(tree->ss[0], rr);
343         rd[1]= Rd_rsquare(tree->ss[1], rr);
344         rd[2]= Rd_rsquare(tree->ss[2], rr);
345 #else
346         approximate_Rd_rgb(tree->ss, rr, rd);
347 #endif
348
349         rd[0] *= area;
350         rd[1] *= area;
351         rd[2] *= area;
352
353         if(frontrad) {
354                 result->radiance[0] += frontrad[0]*rd[0];
355                 result->radiance[1] += frontrad[1]*rd[1];
356                 result->radiance[2] += frontrad[2]*rd[2];
357         }
358         if(backrad) {
359                 result->backradiance[0] += backrad[0]*rd[0];
360                 result->backradiance[1] += backrad[1]*rd[1];
361                 result->backradiance[2] += backrad[2]*rd[2];
362         }
363
364         result->rdsum[0] += rd[0];
365         result->rdsum[1] += rd[1];
366         result->rdsum[2] += rd[2];
367 }
368
369 static void traverse_octree(ScatterTree *tree, ScatterNode *node, float *co, int self, ScatterResult *result)
370 {
371         float sub[3], dist;
372         int i, index = 0;
373
374         if(node->totpoint > 0) {
375                 /* leaf - add radiance from all samples */
376                 for(i=0; i<node->totpoint; i++) {
377                         ScatterPoint *p= &node->points[i];
378
379                         VECSUB(sub, co, p->co);
380                         dist= INPR(sub, sub);
381
382                         if(p->back)
383                                 add_radiance(tree, NULL, p->radiance, p->area, dist, result);
384                         else
385                                 add_radiance(tree, p->radiance, NULL, p->area, dist, result);
386                 }
387         }
388         else {
389                 /* branch */
390                 if (self)
391                         index = SUBNODE_INDEX(co, node->split);
392
393                 for(i=0; i<8; i++) {
394                         if(node->child[i]) {
395                                 ScatterNode *subnode= node->child[i];
396
397                                 if(self && index == i) {
398                                         /* always traverse node containing the point */
399                                         traverse_octree(tree, subnode, co, 1, result);
400                                 }
401                                 else {
402                                         /* decide subnode traversal based on maximum solid angle */
403                                         VECSUB(sub, co, subnode->co);
404                                         dist= INPR(sub, sub);
405
406                                         /* actually area/dist > error, but this avoids division */
407                                         if(subnode->area>tree->error*dist) {
408                                                 traverse_octree(tree, subnode, co, 0, result);
409                                         }
410                                         else {
411                                                 add_radiance(tree, subnode->radiance,
412                                                         subnode->backradiance, subnode->area, dist, result);
413                                         }
414                                 }
415                         }
416                 }
417         }
418 }
419
420 static void compute_radiance(ScatterTree *tree, float *co, float *radiance)
421 {
422         ScatterResult result;
423         float rdsum[3];
424
425         memset(&result, 0, sizeof(result));
426
427         traverse_octree(tree, tree->root, co, 1, &result);
428
429         VecMulf(result.radiance, tree->ss[0]->frontweight);
430         VecMulf(result.backradiance, tree->ss[0]->backweight);
431
432         VECCOPY(radiance, result.backradiance);
433         VECADD(radiance, radiance, result.radiance);
434
435         VECCOPY(rdsum, result.rdsum);
436
437         /* the original paper doesn't do this, but we normalize over the
438            sampled area and multiply with the reflectance. this is because
439            our point samples are incomplete, there are no samples on parts
440            of the mesh not visible from the camera. this can not only make
441            it darker, but also lead to ugly color shifts */
442         if(rdsum[0] > 0.0f) radiance[0]= tree->ss[0]->color*radiance[0]/rdsum[0];
443         if(rdsum[1] > 0.0f) radiance[1]= tree->ss[1]->color*radiance[1]/rdsum[1];
444         if(rdsum[2] > 0.0f) radiance[2]= tree->ss[2]->color*radiance[2]/rdsum[2];
445 }
446
447 /* building */
448
449 static void sum_leaf_radiance(ScatterTree *tree, ScatterNode *node)
450 {
451         ScatterPoint *p;
452         float radiance, totradiance= 0.0f, inv;
453         int i;
454
455         node->co[0]= node->co[1]= node->co[2]= 0.0;
456         node->radiance[0]= node->radiance[1]= node->radiance[2]= 0.0;
457         node->backradiance[0]= node->backradiance[1]= node->backradiance[2]= 0.0;
458
459         /* compute total radiance, radiance weighted average position,
460            and total area */
461         for(i=0; i<node->totpoint; i++) {
462                 p= &node->points[i];
463
464                 radiance= p->area*(p->radiance[0] + p->radiance[1] + p->radiance[2]);
465                 totradiance += radiance;
466
467                 node->co[0] += radiance*p->co[0];
468                 node->co[1] += radiance*p->co[1];
469                 node->co[2] += radiance*p->co[2];
470
471                 if(p->back) {
472                         node->backradiance[0] += p->radiance[0]*p->area;
473                         node->backradiance[1] += p->radiance[1]*p->area;
474                         node->backradiance[2] += p->radiance[2]*p->area;
475                 }
476                 else {
477                         node->radiance[0] += p->radiance[0]*p->area;
478                         node->radiance[1] += p->radiance[1]*p->area;
479                         node->radiance[2] += p->radiance[2]*p->area;
480                 }
481
482                 node->area += p->area;
483         }
484
485         if(node->area > 0) {
486                 inv= 1.0/node->area;
487                 node->radiance[0] *= inv;
488                 node->radiance[1] *= inv;
489                 node->radiance[2] *= inv;
490                 node->backradiance[0] *= inv;
491                 node->backradiance[1] *= inv;
492                 node->backradiance[2] *= inv;
493         }
494
495         if(totradiance > 0.0f) {
496                 inv= 1.0/totradiance;
497                 node->co[0] *= inv;
498                 node->co[1] *= inv;
499                 node->co[2] *= inv;
500         }
501         else {
502                 /* make sure that if radiance is 0.0f, we still have these points in
503                    the tree at a good position, they count for rdsum too */
504                 for(i=0; i<node->totpoint; i++) {
505                         p= &node->points[i];
506
507                         node->co[0] += p->co[0];
508                         node->co[1] += p->co[1];
509                         node->co[2] += p->co[2];
510                 }
511
512                 node->co[0] /= node->totpoint;
513                 node->co[1] /= node->totpoint;
514                 node->co[2] /= node->totpoint;
515         }
516 }
517
518 static void sum_branch_radiance(ScatterTree *tree, ScatterNode *node)
519 {
520         ScatterNode *subnode;
521         float radiance, totradiance= 0.0f, inv;
522         int i, totnode;
523
524         node->co[0]= node->co[1]= node->co[2]= 0.0;
525         node->radiance[0]= node->radiance[1]= node->radiance[2]= 0.0;
526         node->backradiance[0]= node->backradiance[1]= node->backradiance[2]= 0.0;
527
528         /* compute total radiance, radiance weighted average position,
529            and total area */
530         for(i=0; i<8; i++) {
531                 if(node->child[i] == NULL)
532                         continue;
533
534                 subnode= node->child[i];
535
536                 radiance= subnode->area*(subnode->radiance[0] + subnode->radiance[1] + subnode->radiance[2] + subnode->backradiance[0] + subnode->backradiance[1] + subnode->backradiance[2]);
537                 totradiance += radiance;
538
539                 node->co[0] += radiance*subnode->co[0];
540                 node->co[1] += radiance*subnode->co[1];
541                 node->co[2] += radiance*subnode->co[2];
542
543                 node->radiance[0] += subnode->radiance[0]*subnode->area;
544                 node->radiance[1] += subnode->radiance[1]*subnode->area;
545                 node->radiance[2] += subnode->radiance[2]*subnode->area;
546
547                 node->backradiance[0] += subnode->backradiance[0]*subnode->area;
548                 node->backradiance[1] += subnode->backradiance[1]*subnode->area;
549                 node->backradiance[2] += subnode->backradiance[2]*subnode->area;
550
551                 node->area += subnode->area;
552         }
553
554         if(node->area > 0) {
555                 inv= 1.0/node->area;
556                 node->radiance[0] *= inv;
557                 node->radiance[1] *= inv;
558                 node->radiance[2] *= inv;
559                 node->backradiance[0] *= inv;
560                 node->backradiance[1] *= inv;
561                 node->backradiance[2] *= inv;
562         }
563
564         if(totradiance > 0.0f) {
565                 inv= 1.0/totradiance;
566                 node->co[0] *= inv;
567                 node->co[1] *= inv;
568                 node->co[2] *= inv;
569         }
570         else {
571                 /* make sure that if radiance is 0.0f, we still have these points in
572                    the tree at a good position, they count for rdsum too */
573                 totnode= 0;
574
575                 for(i=0; i<8; i++) {
576                         if(node->child[i]) {
577                                 subnode= node->child[i];
578
579                                 node->co[0] += subnode->co[0];
580                                 node->co[1] += subnode->co[1];
581                                 node->co[2] += subnode->co[2];
582
583                                 totnode++;
584                         }
585                 }
586
587                 node->co[0] /= totnode;
588                 node->co[1] /= totnode;
589                 node->co[2] /= totnode;
590         }
591 }
592
593 static void sum_radiance(ScatterTree *tree, ScatterNode *node)
594 {
595         if(node->totpoint > 0) {
596                 sum_leaf_radiance(tree, node);
597         }
598         else {
599                 int i;
600
601                 for(i=0; i<8; i++)
602                         if(node->child[i])
603                                 sum_radiance(tree, node->child[i]);
604
605                 sum_branch_radiance(tree, node);
606         }
607 }
608
609 static void subnode_middle(int i, float *mid, float *subsize, float *submid)
610 {
611         int x= i & 1, y= i & 2, z= i & 4;
612
613         submid[0]= mid[0] + ((x)? subsize[0]: -subsize[0]);
614         submid[1]= mid[1] + ((y)? subsize[1]: -subsize[1]);
615         submid[2]= mid[2] + ((z)? subsize[2]: -subsize[2]);
616 }
617
618 static void create_octree_node(ScatterTree *tree, ScatterNode *node, float *mid, float *size, ScatterPoint **refpoints, int depth)
619 {
620         ScatterNode *subnode;
621         ScatterPoint **subrefpoints, **tmppoints= tree->tmppoints;
622         int index, nsize[8], noffset[8], i, subco, usednodes, usedi;
623         float submid[3], subsize[3];
624
625         /* stopping condition */
626         if(node->totpoint <= MAX_OCTREE_NODE_POINTS || depth == MAX_OCTREE_DEPTH) {
627                 for(i=0; i<node->totpoint; i++)
628                         node->points[i]= *(refpoints[i]);
629
630                 return;
631         }
632
633         subsize[0]= size[0]*0.5;
634         subsize[1]= size[1]*0.5;
635         subsize[2]= size[2]*0.5;
636
637         node->split[0]= mid[0];
638         node->split[1]= mid[1];
639         node->split[2]= mid[2];
640
641         memset(nsize, 0, sizeof(nsize));
642         memset(noffset, 0, sizeof(noffset));
643
644         /* count points in subnodes */
645         for(i=0; i<node->totpoint; i++) {
646                 index= SUBNODE_INDEX(refpoints[i]->co, node->split);
647                 tmppoints[i]= refpoints[i];
648                 nsize[index]++;
649         }
650
651         /* here we check if only one subnode is used. if this is the case, we don't
652            create a new node, but rather call this function again, with different 
653            size and middle position for the same node. */
654         for(usedi=0, usednodes=0, i=0; i<8; i++) {
655                 if(nsize[i]) {
656                         usednodes++;
657                         usedi = i;
658                 }
659                 if(i != 0)
660                         noffset[i]= noffset[i-1]+nsize[i-1];
661         }
662         
663         if(usednodes<=1) {
664                 subnode_middle(usedi, mid, subsize, submid);
665                 create_octree_node(tree, node, submid, subsize, refpoints, depth+1);
666                 return;
667         }
668
669         /* reorder refpoints by subnode */
670         for(i=0; i<node->totpoint; i++) {
671                 index= SUBNODE_INDEX(tmppoints[i]->co, node->split);
672                 refpoints[noffset[index]]= tmppoints[i];
673                 noffset[index]++;
674         }
675
676         /* create subnodes */
677         for(subco=0, i=0; i<8; subco+=nsize[i], i++) {
678                 if(nsize[i] > 0) {
679                         subnode= BLI_memarena_alloc(tree->arena, sizeof(ScatterNode));
680                         node->child[i]= subnode;
681                         subnode->points= node->points + subco;
682                         subnode->totpoint= nsize[i];
683                         subrefpoints= refpoints + subco;
684
685                         subnode_middle(i, mid, subsize, submid);
686
687                         create_octree_node(tree, subnode, submid, subsize, subrefpoints,
688                                 depth+1);
689                 }
690                 else
691                         node->child[i]= NULL;
692         }
693
694         node->points= NULL;
695         node->totpoint= 0;
696 }
697
698 /* public functions */
699
700 ScatterTree *scatter_tree_new(ScatterSettings *ss[3], float scale, float error,
701         float (*co)[3], float (*color)[3], float *area, int totpoint)
702 {
703         ScatterTree *tree;
704         ScatterPoint *points, **refpoints;
705         int i;
706
707         /* allocate tree */
708         tree= MEM_callocN(sizeof(ScatterTree), "ScatterTree");
709         tree->scale= scale;
710         tree->error= error;
711         tree->totpoint= totpoint;
712
713         tree->ss[0]= ss[0];
714         tree->ss[1]= ss[1];
715         tree->ss[2]= ss[2];
716
717         points= MEM_callocN(sizeof(ScatterPoint)*totpoint, "ScatterPoints");
718         refpoints= MEM_callocN(sizeof(ScatterPoint*)*totpoint, "ScatterRefPoints");
719
720         tree->points= points;
721         tree->refpoints= refpoints;
722
723         /* build points */
724         INIT_MINMAX(tree->min, tree->max);
725
726         for(i=0; i<totpoint; i++) {
727                 VECCOPY(points[i].co, co[i]);
728                 VECCOPY(points[i].radiance, color[i]);
729                 points[i].area= fabs(area[i])/(tree->scale*tree->scale);
730                 points[i].back= (area[i] < 0.0f);
731
732                 VecMulf(points[i].co, 1.0f/tree->scale);
733                 DO_MINMAX(points[i].co, tree->min, tree->max);
734
735                 refpoints[i]= points + i;
736         }
737
738         return tree;
739 }
740
741 void scatter_tree_build(ScatterTree *tree)
742 {
743         ScatterPoint *newpoints, **tmppoints;
744         float mid[3], size[3];
745         int totpoint= tree->totpoint;
746
747         newpoints= MEM_callocN(sizeof(ScatterPoint)*totpoint, "ScatterPoints");
748         tmppoints= MEM_callocN(sizeof(ScatterPoint*)*totpoint, "ScatterTmpPoints");
749         tree->tmppoints= tmppoints;
750
751         tree->arena= BLI_memarena_new(0x8000 * sizeof(ScatterNode));
752         BLI_memarena_use_calloc(tree->arena);
753
754         /* build tree */
755         tree->root= BLI_memarena_alloc(tree->arena, sizeof(ScatterNode));
756         tree->root->points= newpoints;
757         tree->root->totpoint= totpoint;
758
759         mid[0]= (tree->min[0]+tree->max[0])*0.5;
760         mid[1]= (tree->min[1]+tree->max[1])*0.5;
761         mid[2]= (tree->min[2]+tree->max[2])*0.5;
762
763         size[0]= (tree->max[0]-tree->min[0])*0.5;
764         size[1]= (tree->max[1]-tree->min[1])*0.5;
765         size[2]= (tree->max[2]-tree->min[2])*0.5;
766
767         create_octree_node(tree, tree->root, mid, size, tree->refpoints, 0);
768
769         MEM_freeN(tree->points);
770         MEM_freeN(tree->refpoints);
771         MEM_freeN(tree->tmppoints);
772         tree->refpoints= NULL;
773         tree->tmppoints= NULL;
774         tree->points= newpoints;
775         
776         /* sum radiance at nodes */
777         sum_radiance(tree, tree->root);
778 }
779
780 void scatter_tree_sample(ScatterTree *tree, float *co, float *color)
781 {
782         float sco[3];
783
784         VECCOPY(sco, co);
785         VecMulf(sco, 1.0f/tree->scale);
786
787         compute_radiance(tree, sco, color);
788 }
789
790 void scatter_tree_free(ScatterTree *tree)
791 {
792         if (tree->arena) BLI_memarena_free(tree->arena);
793         if (tree->points) MEM_freeN(tree->points);
794         if (tree->refpoints) MEM_freeN(tree->refpoints);
795                 
796         MEM_freeN(tree);
797 }
798
799 /* Internal Renderer API */
800
801 /* sss tree building */
802
803 typedef struct SSSData {
804         ScatterTree *tree;
805         ScatterSettings *ss[3];
806 } SSSData;
807
808 typedef struct SSSPoints {
809         struct SSSPoints *next, *prev;
810
811         float (*co)[3];
812         float (*color)[3];
813         float *area;
814         int totpoint;
815 } SSSPoints;
816
817 static void sss_create_tree_mat(Render *re, Material *mat)
818 {
819         SSSPoints *p;
820         ListBase layers, points;
821         float (*co)[3] = NULL, (*color)[3] = NULL, *area = NULL;
822         int totpoint = 0, osa;
823
824         if(re->test_break())
825                 return;
826         
827         points.first= points.last= NULL;
828
829         /* do SSS preprocessing render */
830         layers= re->r.layers;
831         osa= re->osa;
832
833         re->r.layers.first= re->r.layers.last= NULL;
834         re->osa= 0;
835         re->sss_points= &points;
836         re->sss_mat= mat;
837
838         RE_TileProcessor(re, 0, 1);
839
840         re->sss_mat= NULL;
841         re->sss_points= NULL;
842         re->r.layers= layers;
843         re->osa= osa;
844
845         /* no points? no tree */
846         if(!points.first)
847                 return;
848
849         /* merge points together into a single buffer */
850         if(!re->test_break()) {
851                 for(totpoint=0, p=points.first; p; p=p->next)
852                         totpoint += p->totpoint;
853                 
854                 co= MEM_mallocN(sizeof(*co)*totpoint, "SSSCo");
855                 color= MEM_mallocN(sizeof(*color)*totpoint, "SSSColor");
856                 area= MEM_mallocN(sizeof(*area)*totpoint, "SSSArea");
857
858                 for(totpoint=0, p=points.first; p; p=p->next) {
859                         memcpy(co+totpoint, p->co, sizeof(*co)*p->totpoint);
860                         memcpy(color+totpoint, p->color, sizeof(*color)*p->totpoint);
861                         memcpy(area+totpoint, p->area, sizeof(*area)*p->totpoint);
862                         totpoint += p->totpoint;
863                 }
864         }
865
866         /* free points */
867         for(p=points.first; p; p=p->next) {
868                 MEM_freeN(p->co);
869                 MEM_freeN(p->color);
870                 MEM_freeN(p->area);
871         }
872         BLI_freelistN(&points);
873
874         /* build tree */
875         if(!re->test_break()) {
876                 SSSData *sss= MEM_callocN(sizeof(*sss), "SSSData");
877                 float ior= mat->sss_ior, cfac= mat->sss_colfac;
878                 float *col= mat->sss_col, *radius= mat->sss_radius;
879                 float fw= mat->sss_front, bw= mat->sss_back;
880
881                 sss->ss[0]= scatter_settings_new(col[0], radius[0], ior, cfac, fw, bw);
882                 sss->ss[1]= scatter_settings_new(col[1], radius[1], ior, cfac, fw, bw);
883                 sss->ss[2]= scatter_settings_new(col[2], radius[2], ior, cfac, fw, bw);
884                 sss->tree= scatter_tree_new(sss->ss, mat->sss_scale, mat->sss_error,
885                         co, color, area, totpoint);
886
887                 MEM_freeN(co);
888                 MEM_freeN(color);
889                 MEM_freeN(area);
890
891                 scatter_tree_build(sss->tree);
892
893                 BLI_ghash_insert(re->sss_hash, mat, sss);
894         }
895         else {
896                 if (co) MEM_freeN(co);
897                 if (color) MEM_freeN(color);
898                 if (area) MEM_freeN(area);
899         }
900 }
901
902 void sss_add_points(Render *re, float (*co)[3], float (*color)[3], float *area, int totpoint)
903 {
904         SSSPoints *p;
905         
906         if(totpoint > 0) {
907                 p= MEM_callocN(sizeof(SSSPoints), "SSSPoints");
908
909                 p->co= co;
910                 p->color= color;
911                 p->area= area;
912                 p->totpoint= totpoint;
913
914                 BLI_lock_thread(LOCK_CUSTOM1);
915                 BLI_addtail(re->sss_points, p);
916                 BLI_unlock_thread(LOCK_CUSTOM1);
917         }
918 }
919
920 static void sss_free_tree(SSSData *sss)
921 {
922         scatter_tree_free(sss->tree);
923         scatter_settings_free(sss->ss[0]);
924         scatter_settings_free(sss->ss[1]);
925         scatter_settings_free(sss->ss[2]);
926         MEM_freeN(sss);
927 }
928
929 /* public functions */
930
931 void make_sss_tree(Render *re)
932 {
933         Material *mat;
934         
935         re->sss_hash= BLI_ghash_new(BLI_ghashutil_ptrhash, BLI_ghashutil_ptrcmp);
936
937         re->i.infostr= "SSS preprocessing";
938         re->stats_draw(&re->i);
939         
940         for(mat= G.main->mat.first; mat; mat= mat->id.next)
941                 if(mat->id.us && (mat->sss_flag & MA_DIFF_SSS))
942                         sss_create_tree_mat(re, mat);
943 }
944
945 void free_sss(Render *re)
946 {
947         if(re->sss_hash) {
948                 GHashIterator *it= BLI_ghashIterator_new(re->sss_hash);
949
950                 while(!BLI_ghashIterator_isDone(it)) {
951                         sss_free_tree(BLI_ghashIterator_getValue(it));
952                         BLI_ghashIterator_step(it);
953                 }
954
955                 BLI_ghashIterator_free(it);
956                 BLI_ghash_free(re->sss_hash, NULL, NULL);
957                 re->sss_hash= NULL;
958         }
959 }
960
961 int sample_sss(Render *re, Material *mat, float *co, float *color)
962 {
963         if(re->sss_hash) {
964                 SSSData *sss= BLI_ghash_lookup(re->sss_hash, mat);
965
966                 if(sss) {
967                         scatter_tree_sample(sss->tree, co, color);
968                         return 1;
969                 }
970                 else {
971                         color[0]= 0.0f;
972                         color[1]= 0.0f;
973                         color[2]= 0.0f;
974                 }
975         }
976
977         return 0;
978 }
979
980 int has_sss_tree(struct Render *re, struct Material *mat)
981 {
982         return (re->sss_hash && BLI_ghash_lookup(re->sss_hash, mat));
983 }
984