Merged changes in the trunk up to revision 46045.
[blender-staging.git] / source / blender / render / intern / source / sss.c
1 /*
2  *
3  * ***** BEGIN GPL LICENSE BLOCK *****
4  *
5  * This program is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU General Public License
7  * as published by the Free Software Foundation; either version 2
8  * of the License, or (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program; if not, write to the Free Software Foundation,
17  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18  *
19  * The Original Code is Copyright (C) 2007 Blender Foundation.
20  * All rights reserved.
21  *
22  * The Original Code is: all of this file.
23  *
24  * Contributor(s): none yet.
25  *
26  * ***** END GPL LICENSE BLOCK *****
27  */
28
29 /** \file blender/render/intern/source/sss.c
30  *  \ingroup render
31  */
32
33
34 /* Possible Improvements:
35  * - add fresnel terms
36  * - adapt Rd table to scale, now with small scale there are a lot of misses?
37  * - possible interesting method: perform sss on all samples in the tree,
38  *   and then use those values interpolated somehow later. can also do this
39  *   filtering on demand for speed. since we are doing things in screen
40  *   space now there is an exact correspondence
41  * - avoid duplicate shading (filtering points in advance, irradiance cache
42  *   like lookup?)
43  * - lower resolution samples
44  */
45
46 #include <math.h>
47 #include <string.h>
48 #include <stdio.h>
49 #include <string.h>
50
51 /* external modules: */
52 #include "MEM_guardedalloc.h"
53
54 #include "BLI_math.h"
55 #include "BLI_blenlib.h"
56 #include "BLI_utildefines.h"
57 #include "BLI_ghash.h"
58 #include "BLI_memarena.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         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= MIN2(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, float *co, 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, float *co, 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, usednodes, 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, usednodes=0, i=0; i<8; i++) {
690                 if (nsize[i]) {
691                         usednodes++;
692                         usedi = i;
693                 }
694                 if (i != 0)
695                         noffset[i]= noffset[i-1]+nsize[i-1];
696         }
697         
698         if (usednodes<=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                 DO_MINMAX(points[i].co, tree->min, tree->max);
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, float *co, float *color)
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= 0;
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_new(BLI_ghashutil_ptrhash, BLI_ghashutil_ptrcmp, "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, float *co, float *color)
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