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