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