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