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