add BLI_strcpy_rlen, replace strcat, which was used in misleading way.
[blender.git] / intern / cycles / render / bssrdf.cpp
1 /*
2  * Copyright 2011, Blender Foundation.
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
19 #include "bssrdf.h"
20
21 #include "util_algorithm.h"
22 #include "util_math.h"
23 #include "util_types.h"
24
25 #include "kernel_types.h"
26 #include "kernel_montecarlo.h"
27
28 #include "closure/bsdf_diffuse.h"
29 #include "closure/bssrdf.h"
30
31 CCL_NAMESPACE_BEGIN
32
33 /* Cumulative density function utilities */
34
35 static float cdf_lookup_inverse(const vector<float>& table, float2 range, float x)
36 {
37         int index = upper_bound(table.begin(), table.end(), x) - table.begin();
38
39         if(index == 0)
40                 return range[0];
41         else if(index == table.size())
42                 return range[1];
43         else
44                 index--;
45         
46         float t = (x - table[index])/(table[index+1] - table[index]);
47         float y = ((index + t)/(table.size() - 1));
48
49         return y*(range[1] - range[0]) + range[0];
50 }
51
52 static void cdf_invert(vector<float>& to, float2 to_range, const vector<float>& from, float2 from_range)
53 {
54         float step = 1.0f/(float)(to.size() - 1);
55
56         for(int i = 0; i < to.size(); i++) {
57                 float x = (i*step)*(from_range[1] - from_range[0]) + from_range[0];
58                 to[i] = cdf_lookup_inverse(from, to_range, x);
59         }
60 }
61
62 /* BSSRDF */
63
64 static float bssrdf_lookup_table_max_radius(const BSSRDFParams *ss)
65 {
66         /* todo: adjust when we use the real BSSRDF */
67         return ss->ld;
68 }
69
70 static void bssrdf_lookup_table_create(const BSSRDFParams *ss, vector<float>& sample_table, vector<float>& pdf_table)
71 {
72         const int size = BSSRDF_RADIUS_TABLE_SIZE;
73         vector<float> cdf(size);
74         vector<float> pdf(size);
75         float step = 1.0f/(float)(size - 1);
76         float max_radius = bssrdf_lookup_table_max_radius(ss);
77         float pdf_sum = 0.0f;
78
79         /* compute the probability density function */
80         for(int i = 0; i < pdf.size(); i++) {
81                 float x = (i*step)*max_radius;
82                 pdf[i] = bssrdf_cubic(ss->ld, x);
83                 pdf_sum += pdf[i];
84         }
85
86         /* adjust for area covered by each distance */
87         for(int i = 0; i < pdf.size(); i++) {
88                 float x = (i*step)*max_radius;
89                 pdf[i] *= M_2PI_F*x;
90         }
91
92         /* normalize pdf, we multiply in reflectance later */
93         if(pdf_sum > 0.0f)
94                 for(int i = 0; i < pdf.size(); i++)
95                         pdf[i] /= pdf_sum;
96
97         /* sum to account for sampling which uses overlapping sphere */
98         for(int i = pdf.size() - 2; i >= 0; i--)
99                 pdf[i] = pdf[i] + pdf[i+1];
100
101         /* compute the cumulative density function */
102         cdf[0] = 0.0f;
103
104         for(int i = 1; i < size; i++)
105                 cdf[i] = cdf[i-1] + 0.5f*(pdf[i-1] + pdf[i])*step*max_radius;
106         
107         /* invert cumulative density function for importance sampling */
108         float2 cdf_range = make_float2(0.0f, cdf[size - 1]);
109         float2 table_range = make_float2(0.0f, max_radius);
110
111         cdf_invert(sample_table, table_range, cdf, cdf_range);
112
113         /* copy pdf table */
114         for(int i = 0; i < pdf.size(); i++)
115                 pdf_table[i] = pdf[i];
116 }
117
118 void bssrdf_table_build(vector<float>& table)
119 {
120         vector<float> sample_table(BSSRDF_RADIUS_TABLE_SIZE);
121         vector<float> pdf_table(BSSRDF_RADIUS_TABLE_SIZE);
122
123         table.resize(BSSRDF_LOOKUP_TABLE_SIZE);
124
125         /* create a 2D lookup table, for reflection x sample radius */
126         for(int i = 0; i < BSSRDF_REFL_TABLE_SIZE; i++) {
127                 float refl = (float)i/(float)(BSSRDF_REFL_TABLE_SIZE-1);
128                 float ior = 1.3f;
129                 float radius = 1.0f;
130
131                 BSSRDFParams ss;
132                 bssrdf_setup_params(&ss, refl, radius, ior);
133                 bssrdf_lookup_table_create(&ss, sample_table, pdf_table);
134
135                 memcpy(&table[i*BSSRDF_RADIUS_TABLE_SIZE], &sample_table[0], BSSRDF_RADIUS_TABLE_SIZE*sizeof(float));
136                 memcpy(&table[BSSRDF_PDF_TABLE_OFFSET + i*BSSRDF_RADIUS_TABLE_SIZE], &pdf_table[0], BSSRDF_RADIUS_TABLE_SIZE*sizeof(float));
137         }
138 }
139
140 CCL_NAMESPACE_END
141