\0;115;0cCycles: Cleanup, use ccl_restrict instead of ccl_restrict_ptr
[blender.git] / intern / cycles / kernel / filter / filter_transform.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 ccl_device void kernel_filter_construct_transform(const float *ccl_restrict buffer,
20                                                   int x, int y, int4 rect,
21                                                   int pass_stride,
22                                                   float *transform, int *rank,
23                                                   int radius, float pca_threshold)
24 {
25         int buffer_w = align_up(rect.z - rect.x, 4);
26
27         float features[DENOISE_FEATURES];
28
29         /* Temporary storage, used in different steps of the algorithm. */
30         float tempmatrix[DENOISE_FEATURES*DENOISE_FEATURES];
31         float tempvector[2*DENOISE_FEATURES];
32         const float *ccl_restrict pixel_buffer;
33         int2 pixel;
34
35
36
37
38         /* === Calculate denoising window. === */
39         int2 low  = make_int2(max(rect.x, x - radius),
40                               max(rect.y, y - radius));
41         int2 high = make_int2(min(rect.z, x + radius + 1),
42                               min(rect.w, y + radius + 1));
43
44
45
46
47         /* === Shift feature passes to have mean 0. === */
48         float feature_means[DENOISE_FEATURES];
49         math_vector_zero(feature_means, DENOISE_FEATURES);
50         FOR_PIXEL_WINDOW {
51                 filter_get_features(pixel, pixel_buffer, features, NULL, pass_stride);
52                 math_vector_add(feature_means, features, DENOISE_FEATURES);
53         } END_FOR_PIXEL_WINDOW
54
55         float pixel_scale = 1.0f / ((high.y - low.y) * (high.x - low.x));
56         math_vector_scale(feature_means, pixel_scale, DENOISE_FEATURES);
57
58         /* === Scale the shifted feature passes to a range of [-1; 1], will be baked into the transform later. === */
59         float *feature_scale = tempvector;
60         math_vector_zero(feature_scale, DENOISE_FEATURES);
61
62         FOR_PIXEL_WINDOW {
63                 filter_get_feature_scales(pixel, pixel_buffer, features, feature_means, pass_stride);
64                 math_vector_max(feature_scale, features, DENOISE_FEATURES);
65         } END_FOR_PIXEL_WINDOW
66
67         filter_calculate_scale(feature_scale);
68
69
70         /* === Generate the feature transformation. ===
71          * This transformation maps the DENOISE_FEATURES-dimentional feature space to a reduced feature (r-feature) space
72          * which generally has fewer dimensions. This mainly helps to prevent overfitting. */
73         float* feature_matrix = tempmatrix;
74         math_matrix_zero(feature_matrix, DENOISE_FEATURES);
75         FOR_PIXEL_WINDOW {
76                 filter_get_features(pixel, pixel_buffer, features, feature_means, pass_stride);
77                 math_vector_mul(features, feature_scale, DENOISE_FEATURES);
78                 math_matrix_add_gramian(feature_matrix, DENOISE_FEATURES, features, 1.0f);
79         } END_FOR_PIXEL_WINDOW
80
81         math_matrix_jacobi_eigendecomposition(feature_matrix, transform, DENOISE_FEATURES, 1);
82         *rank = 0;
83         if(pca_threshold < 0.0f) {
84                 float threshold_energy = 0.0f;
85                 for(int i = 0; i < DENOISE_FEATURES; i++) {
86                         threshold_energy += feature_matrix[i*DENOISE_FEATURES+i];
87                 }
88                 threshold_energy *= 1.0f - (-pca_threshold);
89
90                 float reduced_energy = 0.0f;
91                 for(int i = 0; i < DENOISE_FEATURES; i++, (*rank)++) {
92                         if(i >= 2 && reduced_energy >= threshold_energy)
93                                 break;
94                         float s = feature_matrix[i*DENOISE_FEATURES+i];
95                         reduced_energy += s;
96                 }
97         }
98         else {
99                 for(int i = 0; i < DENOISE_FEATURES; i++, (*rank)++) {
100                         float s = feature_matrix[i*DENOISE_FEATURES+i];
101                         if(i >= 2 && sqrtf(s) < pca_threshold)
102                                 break;
103                 }
104         }
105
106         /* Bake the feature scaling into the transformation matrix. */
107         for(int i = 0; i < (*rank); i++) {
108                 math_vector_mul(transform + i*DENOISE_FEATURES, feature_scale, DENOISE_FEATURES);
109         }
110         math_matrix_transpose(transform, DENOISE_FEATURES, 1);
111 }
112
113 CCL_NAMESPACE_END