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