recursive_gaussian_convolution.cc 10 KB

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