1 // Copyright (c) 2013 The Chromium Authors. All rights reserved. 2 // Use of this source code is governed by a BSD-style license that can be 3 // found in the LICENSE file. 4 5 #include <algorithm> 6 #include <cmath> 7 #include <vector> 8 9 #include "base/logging.h" 10 #include "skia/ext/recursive_gaussian_convolution.h" 11 12 namespace skia { 13 14 namespace { 15 16 // Takes the value produced by accumulating element-wise product of image with 17 // a kernel and brings it back into range. 18 // All of the filter scaling factors are in fixed point with kShiftBits bits of 19 // fractional part. 20 template<bool take_absolute> 21 inline unsigned char FloatTo8(float f) { 22 int a = static_cast<int>(f + 0.5f); 23 if (take_absolute) 24 a = std::abs(a); 25 else if (a < 0) 26 return 0; 27 if (a < 256) 28 return a; 29 return 255; 30 } 31 32 template<RecursiveFilter::Order order> 33 inline float ForwardFilter(float in_n_1, 34 float in_n, 35 float in_n1, 36 const std::vector<float>& w, 37 int n, 38 const float* b) { 39 switch (order) { 40 case RecursiveFilter::FUNCTION: 41 return b[0] * in_n + b[1] * w[n-1] + b[2] * w[n-2] + b[3] * w[n-3]; 42 case RecursiveFilter::FIRST_DERIVATIVE: 43 return b[0] * 0.5f * (in_n1 - in_n_1) + 44 b[1] * w[n-1] + b[2] * w[n-2] + b[3] * w[n-3]; 45 case RecursiveFilter::SECOND_DERIVATIVE: 46 return b[0] * (in_n - in_n_1) + 47 b[1] * w[n-1] + b[2] * w[n-2] + b[3] * w[n-3]; 48 } 49 50 NOTREACHED(); 51 return 0.0f; 52 } 53 54 template<RecursiveFilter::Order order> 55 inline float BackwardFilter(const std::vector<float>& out, 56 int n, 57 float w_n, 58 float w_n1, 59 const float* b) { 60 switch (order) { 61 case RecursiveFilter::FUNCTION: 62 case RecursiveFilter::FIRST_DERIVATIVE: 63 return b[0] * w_n + 64 b[1] * out[n + 1] + b[2] * out[n + 2] + b[3] * out[n + 3]; 65 case RecursiveFilter::SECOND_DERIVATIVE: 66 return b[0] * (w_n1 - w_n) + 67 b[1] * out[n + 1] + b[2] * out[n + 2] + b[3] * out[n + 3]; 68 } 69 NOTREACHED(); 70 return 0.0f; 71 } 72 73 template<RecursiveFilter::Order order, bool absolute_values> 74 unsigned char SingleChannelRecursiveFilter( 75 const unsigned char* const source_data, 76 int source_pixel_stride, 77 int source_row_stride, 78 int row_width, 79 int row_count, 80 unsigned char* const output, 81 int output_pixel_stride, 82 int output_row_stride, 83 const float* b) { 84 const int intermediate_buffer_size = row_width + 6; 85 std::vector<float> w(intermediate_buffer_size); 86 const unsigned char* in = source_data; 87 unsigned char* out = output; 88 unsigned char max_output = 0; 89 for (int r = 0; r < row_count; 90 ++r, in += source_row_stride, out += output_row_stride) { 91 // Compute forward filter. 92 // First initialize start of the w (temporary) vector. 93 if (order == RecursiveFilter::FUNCTION) 94 w[0] = w[1] = w[2] = in[0]; 95 else 96 w[0] = w[1] = w[2] = 0.0f; 97 // Note that special-casing of w[3] is needed because of derivatives. 98 w[3] = ForwardFilter<order>( 99 in[0], in[0], in[source_pixel_stride], w, 3, b); 100 int n = 4; 101 int c = 1; 102 int byte_index = source_pixel_stride; 103 for (; c < row_width - 1; ++c, ++n, byte_index += source_pixel_stride) { 104 w[n] = ForwardFilter<order>(in[byte_index - source_pixel_stride], 105 in[byte_index], 106 in[byte_index + source_pixel_stride], 107 w, n, b); 108 } 109 110 // The value of w corresponding to the last image pixel needs to be computed 111 // separately, again because of derivatives. 112 w[n] = ForwardFilter<order>(in[byte_index - source_pixel_stride], 113 in[byte_index], 114 in[byte_index], 115 w, n, b); 116 // Now three trailing bytes set to the same value as current w[n]. 117 w[n + 1] = w[n]; 118 w[n + 2] = w[n]; 119 w[n + 3] = w[n]; 120 121 // Now apply the back filter. 122 float w_n1 = w[n + 1]; 123 int output_index = (row_width - 1) * output_pixel_stride; 124 for (; c >= 0; output_index -= output_pixel_stride, --c, --n) { 125 float w_n = BackwardFilter<order>(w, n, w[n], w_n1, b); 126 w_n1 = w[n]; 127 w[n] = w_n; 128 out[output_index] = FloatTo8<absolute_values>(w_n); 129 max_output = std::max(max_output, out[output_index]); 130 } 131 } 132 return max_output; 133 } 134 135 unsigned char SingleChannelRecursiveFilter( 136 const unsigned char* const source_data, 137 int source_pixel_stride, 138 int source_row_stride, 139 int row_width, 140 int row_count, 141 unsigned char* const output, 142 int output_pixel_stride, 143 int output_row_stride, 144 const float* b, 145 RecursiveFilter::Order order, 146 bool absolute_values) { 147 if (absolute_values) { 148 switch (order) { 149 case RecursiveFilter::FUNCTION: 150 return SingleChannelRecursiveFilter<RecursiveFilter::FUNCTION, true>( 151 source_data, source_pixel_stride, source_row_stride, 152 row_width, row_count, 153 output, output_pixel_stride, output_row_stride, b); 154 case RecursiveFilter::FIRST_DERIVATIVE: 155 return SingleChannelRecursiveFilter< 156 RecursiveFilter::FIRST_DERIVATIVE, true>( 157 source_data, source_pixel_stride, source_row_stride, 158 row_width, row_count, 159 output, output_pixel_stride, output_row_stride, b); 160 case RecursiveFilter::SECOND_DERIVATIVE: 161 return SingleChannelRecursiveFilter< 162 RecursiveFilter::SECOND_DERIVATIVE, true>( 163 source_data, source_pixel_stride, source_row_stride, 164 row_width, row_count, 165 output, output_pixel_stride, output_row_stride, b); 166 } 167 } else { 168 switch (order) { 169 case RecursiveFilter::FUNCTION: 170 return SingleChannelRecursiveFilter<RecursiveFilter::FUNCTION, false>( 171 source_data, source_pixel_stride, source_row_stride, 172 row_width, row_count, 173 output, output_pixel_stride, output_row_stride, b); 174 case RecursiveFilter::FIRST_DERIVATIVE: 175 return SingleChannelRecursiveFilter< 176 RecursiveFilter::FIRST_DERIVATIVE, false>( 177 source_data, source_pixel_stride, source_row_stride, 178 row_width, row_count, 179 output, output_pixel_stride, output_row_stride, b); 180 case RecursiveFilter::SECOND_DERIVATIVE: 181 return SingleChannelRecursiveFilter< 182 RecursiveFilter::SECOND_DERIVATIVE, false>( 183 source_data, source_pixel_stride, source_row_stride, 184 row_width, row_count, 185 output, output_pixel_stride, output_row_stride, b); 186 } 187 } 188 189 NOTREACHED(); 190 return 0; 191 } 192 193 } 194 195 float RecursiveFilter::qFromSigma(float sigma) { 196 DCHECK_GE(sigma, 0.5f); 197 if (sigma <= 2.5f) 198 return 3.97156f - 4.14554f * std::sqrt(1.0f - 0.26891f * sigma); 199 return 0.98711f * sigma - 0.96330f; 200 } 201 202 void RecursiveFilter::computeCoefficients(float q, float b[4]) { 203 b[0] = 1.57825f + 2.44413f * q + 1.4281f * q * q + 0.422205f * q * q * q; 204 b[1] = 2.4413f * q + 2.85619f * q * q + 1.26661f * q * q * q; 205 b[2] = - 1.4281f * q * q - 1.26661f * q * q * q; 206 b[3] = 0.422205f * q * q * q; 207 208 // The above is exactly like in the paper. To cut down on computations, 209 // we can fix up these numbers a bit now. 210 float b_norm = 1.0f - (b[1] + b[2] + b[3]) / b[0]; 211 b[1] /= b[0]; 212 b[2] /= b[0]; 213 b[3] /= b[0]; 214 b[0] = b_norm; 215 } 216 217 RecursiveFilter::RecursiveFilter(float sigma, Order order) 218 : order_(order), q_(qFromSigma(sigma)) { 219 computeCoefficients(q_, b_); 220 } 221 222 unsigned char SingleChannelRecursiveGaussianX(const unsigned char* source_data, 223 int source_byte_row_stride, 224 int input_channel_index, 225 int input_channel_count, 226 const RecursiveFilter& filter, 227 const SkISize& image_size, 228 unsigned char* output, 229 int output_byte_row_stride, 230 int output_channel_index, 231 int output_channel_count, 232 bool absolute_values) { 233 return SingleChannelRecursiveFilter(source_data + input_channel_index, 234 input_channel_count, 235 source_byte_row_stride, 236 image_size.width(), 237 image_size.height(), 238 output + output_channel_index, 239 output_channel_count, 240 output_byte_row_stride, 241 filter.b(), 242 filter.order(), 243 absolute_values); 244 } 245 246 unsigned char SingleChannelRecursiveGaussianY(const unsigned char* source_data, 247 int source_byte_row_stride, 248 int input_channel_index, 249 int input_channel_count, 250 const RecursiveFilter& filter, 251 const SkISize& image_size, 252 unsigned char* output, 253 int output_byte_row_stride, 254 int output_channel_index, 255 int output_channel_count, 256 bool absolute_values) { 257 return SingleChannelRecursiveFilter(source_data + input_channel_index, 258 source_byte_row_stride, 259 input_channel_count, 260 image_size.height(), 261 image_size.width(), 262 output + output_channel_index, 263 output_byte_row_stride, 264 output_channel_count, 265 filter.b(), 266 filter.order(), 267 absolute_values); 268 } 269 270 } // namespace skia 271
