skia/ext/recursive_gaussian_convolution.cc

   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