ClangFormat: apply to source, most of intern
[blender.git] / intern / cycles / kernel / filter / filter_prefilter.h
1 /*
2  * Copyright 2011-2017 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 /* First step of the shadow prefiltering, performs the shadow division and stores all data
20  * in a nice and easy rectangular array that can be passed to the NLM filter.
21  *
22  * Calculates:
23  * unfiltered: Contains the two half images of the shadow feature pass
24  * sampleVariance: The sample-based variance calculated in the kernel. Note: This calculation is biased in general, and especially here since the variance of the ratio can only be approximated.
25  * sampleVarianceV: Variance of the sample variance estimation, quite noisy (since it's essentially the buffer variance of the two variance halves)
26  * bufferVariance: The buffer-based variance of the shadow feature. Unbiased, but quite noisy.
27  */
28 ccl_device void kernel_filter_divide_shadow(int sample,
29                                             CCL_FILTER_TILE_INFO,
30                                             int x,
31                                             int y,
32                                             ccl_global float *unfilteredA,
33                                             ccl_global float *unfilteredB,
34                                             ccl_global float *sampleVariance,
35                                             ccl_global float *sampleVarianceV,
36                                             ccl_global float *bufferVariance,
37                                             int4 rect,
38                                             int buffer_pass_stride,
39                                             int buffer_denoising_offset)
40 {
41   int xtile = (x < tile_info->x[1]) ? 0 : ((x < tile_info->x[2]) ? 1 : 2);
42   int ytile = (y < tile_info->y[1]) ? 0 : ((y < tile_info->y[2]) ? 1 : 2);
43   int tile = ytile * 3 + xtile;
44
45   int offset = tile_info->offsets[tile];
46   int stride = tile_info->strides[tile];
47   const ccl_global float *ccl_restrict center_buffer = (ccl_global float *)ccl_get_tile_buffer(
48       tile);
49   center_buffer += (y * stride + x + offset) * buffer_pass_stride;
50   center_buffer += buffer_denoising_offset + 14;
51
52   int buffer_w = align_up(rect.z - rect.x, 4);
53   int idx = (y - rect.y) * buffer_w + (x - rect.x);
54   unfilteredA[idx] = center_buffer[1] / max(center_buffer[0], 1e-7f);
55   unfilteredB[idx] = center_buffer[4] / max(center_buffer[3], 1e-7f);
56
57   float varA = center_buffer[2];
58   float varB = center_buffer[5];
59   int odd_sample = (sample + 1) / 2;
60   int even_sample = sample / 2;
61
62   /* Approximate variance as E[x^2] - 1/N * (E[x])^2, since online variance
63    * update does not work efficiently with atomics in the kernel. */
64   varA = max(0.0f, varA - unfilteredA[idx] * unfilteredA[idx] * odd_sample);
65   varB = max(0.0f, varB - unfilteredB[idx] * unfilteredB[idx] * even_sample);
66
67   varA /= max(odd_sample - 1, 1);
68   varB /= max(even_sample - 1, 1);
69
70   sampleVariance[idx] = 0.5f * (varA + varB) / sample;
71   sampleVarianceV[idx] = 0.5f * (varA - varB) * (varA - varB) / (sample * sample);
72   bufferVariance[idx] = 0.5f * (unfilteredA[idx] - unfilteredB[idx]) *
73                         (unfilteredA[idx] - unfilteredB[idx]);
74 }
75
76 /* Load a regular feature from the render buffers into the denoise buffer.
77  * Parameters:
78  * - sample: The sample amount in the buffer, used to normalize the buffer.
79  * - m_offset, v_offset: Render Buffer Pass offsets of mean and variance of the feature.
80  * - x, y: Current pixel
81  * - mean, variance: Target denoise buffers.
82  * - rect: The prefilter area (lower pixels inclusive, upper pixels exclusive).
83  */
84 ccl_device void kernel_filter_get_feature(int sample,
85                                           CCL_FILTER_TILE_INFO,
86                                           int m_offset,
87                                           int v_offset,
88                                           int x,
89                                           int y,
90                                           ccl_global float *mean,
91                                           ccl_global float *variance,
92                                           float scale,
93                                           int4 rect,
94                                           int buffer_pass_stride,
95                                           int buffer_denoising_offset)
96 {
97   int xtile = (x < tile_info->x[1]) ? 0 : ((x < tile_info->x[2]) ? 1 : 2);
98   int ytile = (y < tile_info->y[1]) ? 0 : ((y < tile_info->y[2]) ? 1 : 2);
99   int tile = ytile * 3 + xtile;
100   ccl_global float *center_buffer = ((ccl_global float *)ccl_get_tile_buffer(tile)) +
101                                     (tile_info->offsets[tile] + y * tile_info->strides[tile] + x) *
102                                         buffer_pass_stride +
103                                     buffer_denoising_offset;
104
105   int buffer_w = align_up(rect.z - rect.x, 4);
106   int idx = (y - rect.y) * buffer_w + (x - rect.x);
107
108   float val = scale * center_buffer[m_offset];
109   mean[idx] = val;
110
111   if (v_offset >= 0) {
112     if (sample > 1) {
113       /* Approximate variance as E[x^2] - 1/N * (E[x])^2, since online variance
114        * update does not work efficiently with atomics in the kernel. */
115       variance[idx] = max(
116           0.0f, (center_buffer[v_offset] - val * val * sample) / (sample * (sample - 1)));
117     }
118     else {
119       /* Can't compute variance with single sample, just set it very high. */
120       variance[idx] = 1e10f;
121     }
122   }
123 }
124
125 ccl_device void kernel_filter_write_feature(int sample,
126                                             int x,
127                                             int y,
128                                             int4 buffer_params,
129                                             ccl_global float *from,
130                                             ccl_global float *buffer,
131                                             int out_offset,
132                                             int4 rect)
133 {
134   ccl_global float *combined_buffer = buffer + (y * buffer_params.y + x + buffer_params.x) *
135                                                    buffer_params.z;
136
137   int buffer_w = align_up(rect.z - rect.x, 4);
138   int idx = (y - rect.y) * buffer_w + (x - rect.x);
139
140   combined_buffer[out_offset] = from[idx];
141 }
142
143 ccl_device void kernel_filter_detect_outliers(int x,
144                                               int y,
145                                               ccl_global float *image,
146                                               ccl_global float *variance,
147                                               ccl_global float *depth,
148                                               ccl_global float *out,
149                                               int4 rect,
150                                               int pass_stride)
151 {
152   int buffer_w = align_up(rect.z - rect.x, 4);
153
154   int n = 0;
155   float values[25];
156   float pixel_variance, max_variance = 0.0f;
157   for (int y1 = max(y - 2, rect.y); y1 < min(y + 3, rect.w); y1++) {
158     for (int x1 = max(x - 2, rect.x); x1 < min(x + 3, rect.z); x1++) {
159       int idx = (y1 - rect.y) * buffer_w + (x1 - rect.x);
160       float3 color = make_float3(
161           image[idx], image[idx + pass_stride], image[idx + 2 * pass_stride]);
162       color = max(color, make_float3(0.0f, 0.0f, 0.0f));
163       float L = average(color);
164
165       /* Find the position of L. */
166       int i;
167       for (i = 0; i < n; i++) {
168         if (values[i] > L)
169           break;
170       }
171       /* Make space for L by shifting all following values to the right. */
172       for (int j = n; j > i; j--) {
173         values[j] = values[j - 1];
174       }
175       /* Insert L. */
176       values[i] = L;
177       n++;
178
179       float3 pixel_var = make_float3(
180           variance[idx], variance[idx + pass_stride], variance[idx + 2 * pass_stride]);
181       float var = average(pixel_var);
182       if ((x1 == x) && (y1 == y)) {
183         pixel_variance = (pixel_var.x < 0.0f || pixel_var.y < 0.0f || pixel_var.z < 0.0f) ? -1.0f :
184                                                                                             var;
185       }
186       else {
187         max_variance = max(max_variance, var);
188       }
189     }
190   }
191
192   max_variance += 1e-4f;
193
194   int idx = (y - rect.y) * buffer_w + (x - rect.x);
195   float3 color = make_float3(image[idx], image[idx + pass_stride], image[idx + 2 * pass_stride]);
196   color = max(color, make_float3(0.0f, 0.0f, 0.0f));
197   float L = average(color);
198
199   float ref = 2.0f * values[(int)(n * 0.75f)];
200
201   /* Slightly offset values to avoid false positives in (almost) black areas. */
202   max_variance += 1e-5f;
203   ref -= 1e-5f;
204
205   if (L > ref) {
206     /* The pixel appears to be an outlier.
207      * However, it may just be a legitimate highlight. Therefore, it is checked how likely it is that the pixel
208      * should actually be at the reference value:
209      * If the reference is within the 3-sigma interval, the pixel is assumed to be a statistical outlier.
210      * Otherwise, it is very unlikely that the pixel should be darker, which indicates a legitimate highlight.
211      */
212
213     if (pixel_variance < 0.0f || pixel_variance > 9.0f * max_variance) {
214       depth[idx] = -depth[idx];
215       color *= ref / L;
216       variance[idx] = variance[idx + pass_stride] = variance[idx + 2 * pass_stride] = max_variance;
217     }
218     else {
219       float stddev = sqrtf(pixel_variance);
220       if (L - 3 * stddev < ref) {
221         /* The pixel is an outlier, so negate the depth value to mark it as one.
222         * Also, scale its brightness down to the outlier threshold to avoid trouble with the NLM weights. */
223         depth[idx] = -depth[idx];
224         float fac = ref / L;
225         color *= fac;
226         variance[idx] *= fac * fac;
227         variance[idx + pass_stride] *= fac * fac;
228         variance[idx + 2 * pass_stride] *= fac * fac;
229       }
230     }
231   }
232   out[idx] = color.x;
233   out[idx + pass_stride] = color.y;
234   out[idx + 2 * pass_stride] = color.z;
235 }
236
237 /* Combine A/B buffers.
238  * Calculates the combined mean and the buffer variance. */
239 ccl_device void kernel_filter_combine_halves(int x,
240                                              int y,
241                                              ccl_global float *mean,
242                                              ccl_global float *variance,
243                                              ccl_global float *a,
244                                              ccl_global float *b,
245                                              int4 rect,
246                                              int r)
247 {
248   int buffer_w = align_up(rect.z - rect.x, 4);
249   int idx = (y - rect.y) * buffer_w + (x - rect.x);
250
251   if (mean)
252     mean[idx] = 0.5f * (a[idx] + b[idx]);
253   if (variance) {
254     if (r == 0)
255       variance[idx] = 0.25f * (a[idx] - b[idx]) * (a[idx] - b[idx]);
256     else {
257       variance[idx] = 0.0f;
258       float values[25];
259       int numValues = 0;
260       for (int py = max(y - r, rect.y); py < min(y + r + 1, rect.w); py++) {
261         for (int px = max(x - r, rect.x); px < min(x + r + 1, rect.z); px++) {
262           int pidx = (py - rect.y) * buffer_w + (px - rect.x);
263           values[numValues++] = 0.25f * (a[pidx] - b[pidx]) * (a[pidx] - b[pidx]);
264         }
265       }
266       /* Insertion-sort the variances (fast enough for 25 elements). */
267       for (int i = 1; i < numValues; i++) {
268         float v = values[i];
269         int j;
270         for (j = i - 1; j >= 0 && values[j] > v; j--)
271           values[j + 1] = values[j];
272         values[j + 1] = v;
273       }
274       variance[idx] = values[(7 * numValues) / 8];
275     }
276   }
277 }
278
279 CCL_NAMESPACE_END