[Bf-blender-cvs] [dba99c4] soc-2016-cycles_denoising: Cycles: Implement SSE3-optimized denoising kernel
Lukas Stockner
noreply at git.blender.org
Wed Aug 10 03:22:23 CEST 2016
Commit: dba99c499f55c73cb5a9e3c9c9c23f32ea599699
Author: Lukas Stockner
Date: Tue Aug 9 02:24:58 2016 +0200
Branches: soc-2016-cycles_denoising
https://developer.blender.org/rBdba99c499f55c73cb5a9e3c9c9c23f32ea599699
Cycles: Implement SSE3-optimized denoising kernel
===================================================================
M intern/cycles/kernel/kernel_filter.h
M intern/cycles/kernel/kernel_filter_util.h
M intern/cycles/util/util_math_matrix.h
M intern/cycles/util/util_types.h
===================================================================
diff --git a/intern/cycles/kernel/kernel_filter.h b/intern/cycles/kernel/kernel_filter.h
index 81d52d9..ffb7760 100644
--- a/intern/cycles/kernel/kernel_filter.h
+++ b/intern/cycles/kernel/kernel_filter.h
@@ -23,8 +23,254 @@ CCL_NAMESPACE_BEGIN
#define NORM_FEATURE_OFFSET 2
#define NORM_FEATURE_NUM 8
+#ifdef __KERNEL_SSE3__
+ccl_device void kernel_filter_estimate_params(KernelGlobals *kg, int sample, float *buffer, int x, int y, FilterStorage *storage, int4 rect)
+{
+ int buffer_w = align_up(rect.z - rect.x, 4);
+ int pass_stride = (rect.w - rect.y) * buffer_w;
+
+ __m128 features[DENOISE_FEATURES];
+ float *pixel_buffer;
+
+ int2 low = make_int2(max(rect.x, x - kernel_data.integrator.half_window),
+ max(rect.y, y - kernel_data.integrator.half_window));
+ int2 high = make_int2(min(rect.z, x + kernel_data.integrator.half_window + 1),
+ min(rect.w, y + kernel_data.integrator.half_window + 1));
+
+ __m128 feature_means[DENOISE_FEATURES] = {_mm_setzero_ps()};
+ FOR_PIXEL_WINDOW_SSE {
+ filter_get_features_sse(x4, y4, active_pixels, pixel_buffer, features, NULL, pass_stride);
+ math_add_vector_sse(feature_means, DENOISE_FEATURES, features);
+ } END_FOR_PIXEL_WINDOW_SSE
+
+ __m128 pixel_scale = _mm_set1_ps(1.0f / ((high.y - low.y) * (high.x - low.x)));
+ for(int i = 0; i < DENOISE_FEATURES; i++) {
+ feature_means[i] = _mm_mul_ps(_mm_hsum_ps(feature_means[i]), pixel_scale);
+ }
+
+ __m128 feature_scale[DENOISE_FEATURES] = {_mm_setzero_ps()};
+ FOR_PIXEL_WINDOW_SSE {
+ filter_get_features_sse(x4, y4, active_pixels, pixel_buffer, features, feature_means, pass_stride);
+ for(int i = 0; i < DENOISE_FEATURES; i++)
+ feature_scale[i] = _mm_max_ps(feature_scale[i], _mm_fabs_ps(features[i]));
+ } END_FOR_PIXEL_WINDOW_SSE
+ for(int i = 0; i < DENOISE_FEATURES; i++)
+ feature_scale[i] = _mm_rcp_ps(_mm_max_ps(_mm_hmax_ps(feature_scale[i]), _mm_set1_ps(0.01f)));
+
+ __m128 feature_matrix_sse[DENOISE_FEATURES*DENOISE_FEATURES];
+ __m128 feature_matrix_norm = _mm_setzero_ps();
+ math_matrix_zero_lower_sse(feature_matrix_sse, DENOISE_FEATURES);
+ FOR_PIXEL_WINDOW_SSE {
+ filter_get_features_sse(x4, y4, active_pixels, pixel_buffer, features, feature_means, pass_stride);
+ math_mul_vector_sse(features, DENOISE_FEATURES, feature_scale);
+ math_add_gramian_sse(feature_matrix_sse, DENOISE_FEATURES, features, _mm_set1_ps(1.0f));
+
+ filter_get_feature_variance_sse(x4, y4, active_pixels, pixel_buffer, features, feature_scale, pass_stride);
+ math_mul_vector_scalar_sse(features, DENOISE_FEATURES, _mm_set1_ps(kernel_data.integrator.filter_strength));
+ for(int i = 0; i < NORM_FEATURE_NUM; i++)
+ feature_matrix_norm = _mm_add_ps(feature_matrix_norm, features[i + NORM_FEATURE_OFFSET]);
+ } END_FOR_PIXEL_WINDOW_SSE
+ float feature_matrix[DENOISE_FEATURES*DENOISE_FEATURES];
+ math_hsum_matrix_lower(feature_matrix, DENOISE_FEATURES, feature_matrix_sse);
+
+ math_lower_tri_to_full(feature_matrix, DENOISE_FEATURES);
+
+ float *feature_transform = &storage->transform[0], singular[DENOISE_FEATURES];
+ __m128 feature_transform_sse[DENOISE_FEATURES*DENOISE_FEATURES];
+ int rank = svd(feature_matrix, feature_transform, singular, DENOISE_FEATURES);
+ float singular_threshold = 0.01f + 2.0f * (sqrtf(_mm_hsum_ss(feature_matrix_norm)) / (sqrtf(rank) * 0.5f));
+
+ rank = 0;
+ for(int i = 0; i < DENOISE_FEATURES; i++, rank++) {
+ float s = sqrtf(fabsf(singular[i]));
+ if(i >= 2 && s < singular_threshold)
+ break;
+ /* Bake the feature scaling into the transformation matrix. */
+ for(int j = 0; j < DENOISE_FEATURES; j++) {
+ feature_transform[rank*DENOISE_FEATURES + j] *= _mm_cvtss_f32(feature_scale[j]);
+ feature_transform_sse[rank*DENOISE_FEATURES + j] = _mm_set1_ps(feature_transform[rank*DENOISE_FEATURES + j]);
+ }
+ }
+
+ /* From here on, the mean of the features will be shifted to the central pixel's values. */
+ float feature_means_scalar[DENOISE_FEATURES];
+ float *center_buffer = buffer + (y - rect.y) * buffer_w + (x - rect.x);
+ filter_get_features(x, y, center_buffer, feature_means_scalar, NULL, pass_stride);
+ for(int i = 0; i < DENOISE_FEATURES; i++)
+ feature_means[i] = _mm_set1_ps(feature_means_scalar[i]);
+
+
+ /* === Estimate bandwidth for each r-feature dimension. ===
+ * To do so, the second derivative of the pixel color w.r.t. the particular r-feature
+ * has to be estimated. That is done by least-squares-fitting a model that includes
+ * both the r-feature vector z as well as z^T*z and using the resulting parameter for
+ * that dimension of the z^T*z vector times two as the derivative. */
+ int matrix_size = 2*rank + 1; /* Constant term (1 dim) + z (rank dims) + z^T*z (rank dims) */
+ __m128 XtX_sse[(2*DENOISE_FEATURES+1)*(2*DENOISE_FEATURES+1)], design_row[(2*DENOISE_FEATURES+1)];
+ float3 XtY[2*DENOISE_FEATURES+1];
+
+ math_matrix_zero_lower_sse(XtX_sse, matrix_size);
+ math_vec3_zero(XtY, matrix_size);
+ FOR_PIXEL_WINDOW_SSE {
+ filter_get_features_sse(x4, y4, active_pixels, pixel_buffer, features, feature_means, pass_stride);
+ __m128 weight = filter_fill_design_row_sse(features, active_pixels, rank, design_row, feature_transform_sse, NULL);
+ active_pixels = _mm_and_ps(active_pixels, _mm_cmpneq_ps(weight, _mm_setzero_ps()));
+
+ if(!_mm_movemask_ps(active_pixels)) continue;
+ weight = _mm_mul_ps(weight, _mm_rcp_ps(_mm_max_ps(_mm_set1_ps(1.0f), filter_get_pixel_variance_sse(pixel_buffer, active_pixels, pass_stride))));
+
+ math_add_gramian_sse(XtX_sse, matrix_size, design_row, weight);
+
+ __m128 color[3];
+ filter_get_pixel_color_sse(pixel_buffer, active_pixels, color, pass_stride);
+ math_mul_vector_scalar_sse(color, 3, weight);
+ for(int row = 0; row < matrix_size; row++) {
+ __m128 color_row[3] = {color[0], color[1], color[2]};
+ math_mul_vector_scalar_sse(color_row, 3, design_row[row]);
+ XtY[row] += math_sum_float3(color_row);
+ }
+ } END_FOR_PIXEL_WINDOW_SSE
+
+ float XtX[(2*DENOISE_FEATURES+1)*(2*DENOISE_FEATURES+1)];
+ math_hsum_matrix_lower(XtX, matrix_size, XtX_sse);
+
+ /* Solve the normal equation of the linear least squares system: Decompose A = X^T*X into L
+ * so that L is a lower triangular matrix and L*L^T = A. Then, solve
+ * A*x = (L*L^T)*x = L*(L^T*x) = X^T*y by first solving L*b = X^T*y and then L^T*x = b through
+ * forward- and backsubstitution. */
+ math_matrix_add_diagonal(XtX, matrix_size, 1e-3f); /* Improve the numerical stability. */
+ math_cholesky(XtX, matrix_size); /* Decompose A=X^T*x to L. */
+ math_substitute_forward_vec3(XtX, matrix_size, XtY); /* Solve L*b = X^T*y. */
+ math_substitute_back_vec3(XtX, matrix_size, XtY); /* Solve L^T*x = b. */
+
+ /* Calculate the inverse of the r-feature bandwidths. */
+ float *bandwidth_factor = &storage->bandwidth[0];
+ for(int i = 0; i < rank; i++)
+ bandwidth_factor[i] = sqrtf(2.0f * average(fabs(XtY[1+rank+i])) + 0.16f);
+ for(int i = rank; i < DENOISE_FEATURES; i++)
+ bandwidth_factor[i] = 0.0f;
+
+
+ float3 center_color = filter_get_pixel_color(center_buffer, pass_stride);
+ float sqrt_center_variance = sqrtf(filter_get_pixel_variance(center_buffer, pass_stride));
+ __m128 center_color_sse[3] = {_mm_set1_ps(center_color.x), _mm_set1_ps(center_color.y), _mm_set1_ps(center_color.z)};
+ __m128 sqrt_center_variance_sse = _mm_set1_ps(sqrt_center_variance);
+
+ const float candidate_bw[6] = {0.05f, 0.1f, 0.25f, 0.5f, 1.0f, 2.0f};
+ double lsq_bias[LSQ_SIZE], lsq_variance[LSQ_SIZE];
+ math_lsq_init(lsq_bias);
+ math_lsq_init(lsq_variance);
+ for(int g = 0; g < 6; g++) {
+ __m128 g_bandwidth_factor[DENOISE_FEATURES];
+ for(int i = 0; i < rank; i++)
+ /* Divide by the candidate bandwidth since the bandwidth_factor actually is the inverse of the bandwidth. */
+ g_bandwidth_factor[i] = _mm_set1_ps(bandwidth_factor[i]/candidate_bw[g]);
+
+ matrix_size = rank+1;
+ math_matrix_zero_lower_sse(XtX_sse, matrix_size);
+
+ FOR_PIXEL_WINDOW_SSE {
+ __m128 color[3];
+ filter_get_pixel_color_sse(pixel_buffer, active_pixels, color, pass_stride);
+ __m128 variance = filter_get_pixel_variance_sse(pixel_buffer, active_pixels, pass_stride);
+ active_pixels = _mm_and_ps(active_pixels, filter_firefly_rejection_sse(color, variance, center_color_sse, sqrt_center_variance_sse));
+ if(!_mm_movemask_ps(active_pixels)) continue;
+
+ filter_get_features_sse(x4, y4, active_pixels, pixel_buffer, features, feature_means, pass_stride);
+ __m128 weight = filter_fill_design_row_sse(features, active_pixels, rank, design_row, feature_transform_sse, g_bandwidth_factor);
+ active_pixels = _mm_and_ps(active_pixels, _mm_cmpneq_ps(weight, _mm_setzero_ps()));
+ if(!_mm_movemask_ps(active_pixels)) continue;
+
+ weight = _mm_mul_ps(weight, _mm_rcp_ps(_mm_max_ps(_mm_set1_ps(1.0f), variance)));
+
+ math_add_gramian_sse(XtX_sse, matrix_size, design_row, weight);
+ } END_FOR_PIXEL_WINDOW_SSE
+ math_hsum_matrix_lower(XtX, matrix_size, XtX_sse);
+
+ math_matrix_add_diagonal(XtX, matrix_size, 1e-4f); /* Improve the numerical stability. */
+ math_cholesky(XtX, matrix_size);
+ math_inverse_lower_tri_inplace(XtX, matrix_size);
+
+ float r_feature_weight_scalar[DENOISE_FEATURES+1];
+ math_vector_zero(r_feature_weight_scalar, matrix_size);
+ for(int col = 0; col < matrix_size; col++)
+ for(int row = col; row < matrix_size; row++)
+ r_feature_weight_scalar[col] += XtX[row]*XtX[col*matrix_size+row];
+ __m128 r_feature_weight[DENOISE_FEATURES+1];
+ for(int col = 0; col < matrix_size; col++)
+ r_feature_weight[col] = _mm_set1_ps(r_feature_weight_scalar[col]);
+
+ __m128 est_pos_color[3] = {_mm_setzero_ps()}, est_color[3] = {_mm_setzero_ps()};
+ __m128 est_variance = _mm_setzero_ps(), est_pos_variance = _mm_setzero_ps(), pos_weight_sse = _mm_setzero_ps();
+
+ FOR_PIXEL_WINDOW_SSE {
+ __m128 color[3];
+ filter_get_pixel_color_sse(pixel_buffer, active_pixels, color, pass_stride);
+ __m128 variance = filter_get_pixel_variance_sse(pixel_buffer, active_pixels, pass_stride);
+ active_pixels = _mm_and_ps(active_pixels, filter_firefly_rejection_sse(color, variance, center_color_sse, sqrt_center_variance_sse));
+
+ filter_get_features_sse(x4, y4, active_pixels, pixel_buffer, features, feature_means, pass_
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list