Code refactor: tweaks in SSS code to prepare for coming changes.
[blender-staging.git] / intern / cycles / kernel / kernel_subsurface.h
1 /*
2  * Copyright 2011-2013 Blender Foundation
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  * http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16
17 CCL_NAMESPACE_BEGIN
18
19 /* BSSRDF using disk based importance sampling.
20  *
21  * BSSRDF Importance Sampling, SIGGRAPH 2013
22  * http://library.imageworks.com/pdfs/imageworks-library-BSSRDF-sampling.pdf
23  *
24  */
25
26 ccl_device_inline float3 subsurface_scatter_eval(ShaderData *sd,
27                                                  const ShaderClosure *sc,
28                                                  float disk_r,
29                                                  float r,
30                                                  bool all)
31 {
32         /* this is the veach one-sample model with balance heuristic, some pdf
33          * factors drop out when using balance heuristic weighting */
34         float3 eval_sum = make_float3(0.0f, 0.0f, 0.0f);
35         float pdf_sum = 0.0f;
36         float sample_weight_inv = 0.0f;
37
38         if(!all) {
39                 float sample_weight_sum = 0.0f;
40
41                 for(int i = 0; i < sd->num_closure; i++) {
42                         sc = &sd->closure[i];
43
44                         if(CLOSURE_IS_BSSRDF(sc->type)) {
45                                 sample_weight_sum += sc->sample_weight;
46                         }
47                 }
48
49                 sample_weight_inv = 1.0f/sample_weight_sum;
50         }
51
52         for(int i = 0; i < sd->num_closure; i++) {
53                 sc = &sd->closure[i];
54                 
55                 if(CLOSURE_IS_BSSRDF(sc->type)) {
56                         /* in case of branched path integrate we sample all bssrdf's once,
57                          * for path trace we pick one, so adjust pdf for that */
58                         float sample_weight = (all)? 1.0f: sc->sample_weight * sample_weight_inv;
59
60                         /* compute pdf */
61                         float3 eval = bssrdf_eval(sc, r);
62                         float pdf = bssrdf_pdf(sc, disk_r);
63
64                         eval_sum += sc->weight * eval;
65                         pdf_sum += sample_weight * pdf;
66                 }
67         }
68
69         return (pdf_sum > 0.0f)? eval_sum / pdf_sum : make_float3(0.0f, 0.0f, 0.0f);
70 }
71
72 /* replace closures with a single diffuse bsdf closure after scatter step */
73 ccl_device void subsurface_scatter_setup_diffuse_bsdf(KernelGlobals *kg, ShaderData *sd, const ShaderClosure *sc, float3 weight, bool hit, float3 N)
74 {
75         sd->flag &= ~SD_CLOSURE_FLAGS;
76         sd->num_closure = 0;
77         sd->num_closure_left = kernel_data.integrator.max_closures;
78
79         if(hit) {
80                 Bssrdf *bssrdf = (Bssrdf *)sc;
81 #ifdef __PRINCIPLED__
82                 if(bssrdf->type == CLOSURE_BSSRDF_PRINCIPLED_ID) {
83                         PrincipledDiffuseBsdf *bsdf = (PrincipledDiffuseBsdf*)bsdf_alloc(sd, sizeof(PrincipledDiffuseBsdf), weight);
84
85                         if(bsdf) {
86                                 bsdf->N = N;
87                                 bsdf->roughness = bssrdf->roughness;
88                                 sd->flag |= bsdf_principled_diffuse_setup(bsdf);
89
90                                 /* replace CLOSURE_BSDF_PRINCIPLED_DIFFUSE_ID with this special ID so render passes
91                                  * can recognize it as not being a regular Disney principled diffuse closure */
92                                 bsdf->type = CLOSURE_BSDF_BSSRDF_PRINCIPLED_ID;
93                         }
94                 }
95                 else if(CLOSURE_IS_BSDF_BSSRDF(bssrdf->type) ||
96                         CLOSURE_IS_BSSRDF(bssrdf->type))
97 #endif  /* __PRINCIPLED__ */
98                 {
99                         DiffuseBsdf *bsdf = (DiffuseBsdf*)bsdf_alloc(sd, sizeof(DiffuseBsdf), weight);
100
101                         if(bsdf) {
102                                 bsdf->N = N;
103                                 sd->flag |= bsdf_diffuse_setup(bsdf);
104
105                                 /* replace CLOSURE_BSDF_DIFFUSE_ID with this special ID so render passes
106                                  * can recognize it as not being a regular diffuse closure */
107                                 bsdf->type = CLOSURE_BSDF_BSSRDF_ID;
108                         }
109                 }
110         }
111 }
112
113 /* optionally do blurring of color and/or bump mapping, at the cost of a shader evaluation */
114 ccl_device float3 subsurface_color_pow(float3 color, float exponent)
115 {
116         color = max(color, make_float3(0.0f, 0.0f, 0.0f));
117
118         if(exponent == 1.0f) {
119                 /* nothing to do */
120         }
121         else if(exponent == 0.5f) {
122                 color.x = sqrtf(color.x);
123                 color.y = sqrtf(color.y);
124                 color.z = sqrtf(color.z);
125         }
126         else {
127                 color.x = powf(color.x, exponent);
128                 color.y = powf(color.y, exponent);
129                 color.z = powf(color.z, exponent);
130         }
131
132         return color;
133 }
134
135 ccl_device void subsurface_color_bump_blur(KernelGlobals *kg,
136                                            ShaderData *sd,
137                                            ccl_addr_space PathState *state,
138                                            float3 *eval,
139                                            float3 *N)
140 {
141         /* average color and texture blur at outgoing point */
142         float texture_blur;
143         float3 out_color = shader_bssrdf_sum(sd, NULL, &texture_blur);
144
145         /* do we have bump mapping? */
146         bool bump = (sd->flag & SD_HAS_BSSRDF_BUMP) != 0;
147
148         if(bump || texture_blur > 0.0f) {
149                 /* average color and normal at incoming point */
150                 shader_eval_surface(kg, sd, state, state->flag, kernel_data.integrator.max_closures);
151                 float3 in_color = shader_bssrdf_sum(sd, (bump)? N: NULL, NULL);
152
153                 /* we simply divide out the average color and multiply with the average
154                  * of the other one. we could try to do this per closure but it's quite
155                  * tricky to match closures between shader evaluations, their number and
156                  * order may change, this is simpler */
157                 if(texture_blur > 0.0f) {
158                         out_color = subsurface_color_pow(out_color, texture_blur);
159                         in_color = subsurface_color_pow(in_color, texture_blur);
160
161                         *eval *= safe_divide_color(in_color, out_color);
162                 }
163         }
164 }
165
166 /* Subsurface scattering step, from a point on the surface to other
167  * nearby points on the same object.
168  */
169 ccl_device_inline int subsurface_scatter_multi_intersect(
170         KernelGlobals *kg,
171         LocalIntersection *ss_isect,
172         ShaderData *sd,
173         const ShaderClosure *sc,
174         uint *lcg_state,
175         float disk_u,
176         float disk_v,
177         bool all)
178 {
179         /* pick random axis in local frame and point on disk */
180         float3 disk_N, disk_T, disk_B;
181         float pick_pdf_N, pick_pdf_T, pick_pdf_B;
182
183         disk_N = sd->Ng;
184         make_orthonormals(disk_N, &disk_T, &disk_B);
185
186         if(disk_v < 0.5f) {
187                 pick_pdf_N = 0.5f;
188                 pick_pdf_T = 0.25f;
189                 pick_pdf_B = 0.25f;
190                 disk_v *= 2.0f;
191         }
192         else if(disk_v < 0.75f) {
193                 float3 tmp = disk_N;
194                 disk_N = disk_T;
195                 disk_T = tmp;
196                 pick_pdf_N = 0.25f;
197                 pick_pdf_T = 0.5f;
198                 pick_pdf_B = 0.25f;
199                 disk_v = (disk_v - 0.5f)*4.0f;
200         }
201         else {
202                 float3 tmp = disk_N;
203                 disk_N = disk_B;
204                 disk_B = tmp;
205                 pick_pdf_N = 0.25f;
206                 pick_pdf_T = 0.25f;
207                 pick_pdf_B = 0.5f;
208                 disk_v = (disk_v - 0.75f)*4.0f;
209         }
210
211         /* sample point on disk */
212         float phi = M_2PI_F * disk_v;
213         float disk_height, disk_r;
214
215         bssrdf_sample(sc, disk_u, &disk_r, &disk_height);
216
217         float3 disk_P = (disk_r*cosf(phi)) * disk_T + (disk_r*sinf(phi)) * disk_B;
218
219         /* create ray */
220 #ifdef __SPLIT_KERNEL__
221         Ray ray_object = ss_isect->ray;
222         Ray *ray = &ray_object;
223 #else
224         Ray *ray = &ss_isect->ray;
225 #endif
226         ray->P = sd->P + disk_N*disk_height + disk_P;
227         ray->D = -disk_N;
228         ray->t = 2.0f*disk_height;
229         ray->dP = sd->dP;
230         ray->dD = differential3_zero();
231         ray->time = sd->time;
232
233         /* intersect with the same object. if multiple intersections are found it
234          * will use at most BSSRDF_MAX_HITS hits, a random subset of all hits */
235         scene_intersect_local(kg,
236                               *ray,
237                               ss_isect,
238                               sd->object,
239                               lcg_state,
240                               BSSRDF_MAX_HITS);
241         int num_eval_hits = min(ss_isect->num_hits, BSSRDF_MAX_HITS);
242
243         for(int hit = 0; hit < num_eval_hits; hit++) {
244                 /* Quickly retrieve P and Ng without setting up ShaderData. */
245                 float3 hit_P;
246                 if(sd->type & PRIMITIVE_TRIANGLE) {
247                         hit_P = triangle_refine_local(kg,
248                                                       sd,
249                                                       &ss_isect->hits[hit],
250                                                       ray);
251                 }
252 #ifdef __OBJECT_MOTION__
253                 else  if(sd->type & PRIMITIVE_MOTION_TRIANGLE) {
254                         float3 verts[3];
255                         motion_triangle_vertices(
256                                 kg,
257                                 sd->object,
258                                 kernel_tex_fetch(__prim_index, ss_isect->hits[hit].prim),
259                                 sd->time,
260                                 verts);
261                         hit_P = motion_triangle_refine_local(kg,
262                                                              sd,
263                                                              &ss_isect->hits[hit],
264                                                              ray,
265                                                              verts);
266                 }
267 #endif  /* __OBJECT_MOTION__ */
268                 else {
269                         ss_isect->weight[hit] = make_float3(0.0f, 0.0f, 0.0f);
270                         continue;
271                 }
272
273                 float3 hit_Ng = ss_isect->Ng[hit];
274                 if(ss_isect->hits[hit].object != OBJECT_NONE) {
275                         object_normal_transform(kg, sd, &hit_Ng);
276                 }
277
278                 /* Probability densities for local frame axes. */
279                 float pdf_N = pick_pdf_N * fabsf(dot(disk_N, hit_Ng));
280                 float pdf_T = pick_pdf_T * fabsf(dot(disk_T, hit_Ng));
281                 float pdf_B = pick_pdf_B * fabsf(dot(disk_B, hit_Ng));
282
283                 /* Multiple importance sample between 3 axes, power heuristic
284                  * found to be slightly better than balance heuristic. pdf_N
285                  * in the MIS weight and denominator cancelled out. */
286                 float w = pdf_N / (sqr(pdf_N) + sqr(pdf_T) + sqr(pdf_B));
287                 if(ss_isect->num_hits > BSSRDF_MAX_HITS) {
288                         w *= ss_isect->num_hits/(float)BSSRDF_MAX_HITS;
289                 }
290
291                 /* Real distance to sampled point. */
292                 float r = len(hit_P - sd->P);
293
294                 /* Evaluate profiles. */
295                 float3 eval = subsurface_scatter_eval(sd, sc, disk_r, r, all) * w;
296
297                 ss_isect->weight[hit] = eval;
298         }
299
300 #ifdef __SPLIT_KERNEL__
301         ss_isect->ray = *ray;
302 #endif
303
304         return num_eval_hits;
305 }
306
307 ccl_device_noinline void subsurface_scatter_multi_setup(
308         KernelGlobals *kg,
309         LocalIntersection* ss_isect,
310         int hit,
311         ShaderData *sd,
312         ccl_addr_space PathState *state,
313         const ShaderClosure *sc)
314 {
315 #ifdef __SPLIT_KERNEL__
316         Ray ray_object = ss_isect->ray;
317         Ray *ray = &ray_object;
318 #else
319         Ray *ray = &ss_isect->ray;
320 #endif
321
322         /* Workaround for AMD GPU OpenCL compiler. Most probably cache bypass issue. */
323 #if defined(__SPLIT_KERNEL__) && defined(__KERNEL_OPENCL_AMD__) && defined(__KERNEL_GPU__)
324         kernel_split_params.dummy_sd_flag = sd->flag;
325 #endif
326
327         /* Setup new shading point. */
328         shader_setup_from_subsurface(kg, sd, &ss_isect->hits[hit], ray);
329
330         /* Optionally blur colors and bump mapping. */
331         float3 weight = ss_isect->weight[hit];
332         float3 N = sd->N;
333         subsurface_color_bump_blur(kg, sd, state, &weight, &N);
334
335         /* Setup diffuse BSDF. */
336         subsurface_scatter_setup_diffuse_bsdf(kg, sd, sc, weight, true, N);
337 }
338
339 /* subsurface scattering step, from a point on the surface to another nearby point on the same object */
340 ccl_device void subsurface_scatter_step(KernelGlobals *kg, ShaderData *sd, ccl_addr_space PathState *state,
341         const ShaderClosure *sc, uint *lcg_state, float disk_u, float disk_v, bool all)
342 {
343         float3 eval = make_float3(0.0f, 0.0f, 0.0f);
344
345         /* pick random axis in local frame and point on disk */
346         float3 disk_N, disk_T, disk_B;
347         float pick_pdf_N, pick_pdf_T, pick_pdf_B;
348
349         disk_N = sd->Ng;
350         make_orthonormals(disk_N, &disk_T, &disk_B);
351
352         if(disk_v < 0.5f) {
353                 pick_pdf_N = 0.5f;
354                 pick_pdf_T = 0.25f;
355                 pick_pdf_B = 0.25f;
356                 disk_v *= 2.0f;
357         }
358         else if(disk_v < 0.75f) {
359                 float3 tmp = disk_N;
360                 disk_N = disk_T;
361                 disk_T = tmp;
362                 pick_pdf_N = 0.25f;
363                 pick_pdf_T = 0.5f;
364                 pick_pdf_B = 0.25f;
365                 disk_v = (disk_v - 0.5f)*4.0f;
366         }
367         else {
368                 float3 tmp = disk_N;
369                 disk_N = disk_B;
370                 disk_B = tmp;
371                 pick_pdf_N = 0.25f;
372                 pick_pdf_T = 0.25f;
373                 pick_pdf_B = 0.5f;
374                 disk_v = (disk_v - 0.75f)*4.0f;
375         }
376
377         /* sample point on disk */
378         float phi = M_2PI_F * disk_v;
379         float disk_height, disk_r;
380
381         bssrdf_sample(sc, disk_u, &disk_r, &disk_height);
382
383         float3 disk_P = (disk_r*cosf(phi)) * disk_T + (disk_r*sinf(phi)) * disk_B;
384
385         /* create ray */
386         Ray ray;
387         ray.P = sd->P + disk_N*disk_height + disk_P;
388         ray.D = -disk_N;
389         ray.t = 2.0f*disk_height;
390         ray.dP = sd->dP;
391         ray.dD = differential3_zero();
392         ray.time = sd->time;
393
394         /* intersect with the same object. if multiple intersections are
395          * found it will randomly pick one of them */
396         LocalIntersection ss_isect;
397         scene_intersect_local(kg, ray, &ss_isect, sd->object, lcg_state, 1);
398
399         /* evaluate bssrdf */
400         if(ss_isect.num_hits > 0) {
401                 float3 origP = sd->P;
402
403                 /* Workaround for AMD GPU OpenCL compiler. Most probably cache bypass issue. */
404 #if defined(__SPLIT_KERNEL__) && defined(__KERNEL_OPENCL_AMD__) && defined(__KERNEL_GPU__)
405                 kernel_split_params.dummy_sd_flag = sd->flag;
406 #endif
407                 /* setup new shading point */
408                 shader_setup_from_subsurface(kg, sd, &ss_isect.hits[0], &ray);
409
410                 /* Probability densities for local frame axes. */
411                 float pdf_N = pick_pdf_N * fabsf(dot(disk_N, sd->Ng));
412                 float pdf_T = pick_pdf_T * fabsf(dot(disk_T, sd->Ng));
413                 float pdf_B = pick_pdf_B * fabsf(dot(disk_B, sd->Ng));
414
415                 /* Multiple importance sample between 3 axes, power heuristic
416                  * found to be slightly better than balance heuristic. pdf_N
417                  * in the MIS weight and denominator cancelled out. */
418                 float w = pdf_N / (sqr(pdf_N) + sqr(pdf_T) + sqr(pdf_B));
419                 w *= ss_isect.num_hits;
420
421                 /* Real distance to sampled point. */
422                 float r = len(sd->P - origP);
423
424                 /* Evaluate profiles. */
425                 eval = subsurface_scatter_eval(sd, sc, disk_r, r, all) * w;
426         }
427
428         /* optionally blur colors and bump mapping */
429         float3 N = sd->N;
430         subsurface_color_bump_blur(kg, sd, state, &eval, &N);
431
432         /* setup diffuse bsdf */
433         subsurface_scatter_setup_diffuse_bsdf(kg, sd, sc, eval, (ss_isect.num_hits > 0), N);
434 }
435
436 CCL_NAMESPACE_END
437