skia/ext/recursive_gaussian_convolution_unittest.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 <functional>
   6 #include <numeric>
   7 #include <vector>
   8
   9 #include "base/basictypes.h"
  10 #include "base/files/file_path.h"
  11 #include "base/files/file_util.h"
  12 #include "base/logging.h"
  13 #include "base/time/time.h"
  14 #include "skia/ext/convolver.h"
  15 #include "skia/ext/recursive_gaussian_convolution.h"
  16 #include "testing/gtest/include/gtest/gtest.h"
  17 #include "third_party/skia/include/core/SkPoint.h"
  18 #include "third_party/skia/include/core/SkRect.h"
  19
  20 namespace {
  21
  22 int ComputeRowStride(int width, int channel_count, int stride_slack) {
  23   return width * channel_count + stride_slack;
  24 }
  25
  26 SkIPoint MakeImpulseImage(std::vector<unsigned char>* image,
  27                           int width,
  28                           int height,
  29                           int channel_index,
  30                           int channel_count,
  31                           int stride_slack) {
  32   const int src_row_stride = ComputeRowStride(
  33       width, channel_count, stride_slack);
  34   const int src_byte_count = src_row_stride * height;
  35   const int signal_x = width / 2;
  36   const int signal_y = height / 2;
  37
  38   image->resize(src_byte_count, 0);
  39   const int non_zero_pixel_index =
  40       signal_y * src_row_stride + signal_x * channel_count + channel_index;
  41   (*image)[non_zero_pixel_index] = 255;
  42   return SkIPoint::Make(signal_x, signal_y);
  43 }
  44
  45 SkIRect MakeBoxImage(std::vector<unsigned char>* image,
  46                      int width,
  47                      int height,
  48                      int channel_index,
  49                      int channel_count,
  50                      int stride_slack,
  51                      int box_width,
  52                      int box_height,
  53                      unsigned char value) {
  54   const int src_row_stride = ComputeRowStride(
  55       width, channel_count, stride_slack);
  56   const int src_byte_count = src_row_stride * height;
  57   const SkIRect box = SkIRect::MakeXYWH((width - box_width) / 2,
  58                                         (height - box_height) / 2,
  59                                         box_width, box_height);
  60
  61   image->resize(src_byte_count, 0);
  62   for (int y = box.top(); y < box.bottom(); ++y) {
  63     for (int x = box.left(); x < box.right(); ++x)
  64       (*image)[y * src_row_stride + x * channel_count + channel_index] = value;
  65   }
  66
  67   return box;
  68 }
  69
  70 int ComputeBoxSum(const std::vector<unsigned char>& image,
  71                   const SkIRect& box,
  72                   int image_width) {
  73   // Compute the sum of all pixels in the box. Assume byte stride 1 and row
  74   // stride same as image_width.
  75   int sum = 0;
  76   for (int y = box.top(); y < box.bottom(); ++y) {
  77     for (int x = box.left(); x < box.right(); ++x)
  78       sum += image[y * image_width + x];
  79   }
  80
  81   return sum;
  82 }
  83
  84 }  // namespace
  85
  86 namespace skia {
  87
  88 TEST(RecursiveGaussian, SmoothingMethodComparison) {
  89   static const int kImgWidth = 512;
  90   static const int kImgHeight = 220;
  91   static const int kChannelIndex = 3;
  92   static const int kChannelCount = 3;
  93   static const int kStrideSlack = 22;
  94
  95   std::vector<unsigned char> input;
  96   SkISize image_size = SkISize::Make(kImgWidth, kImgHeight);
  97   MakeImpulseImage(
  98       &input, kImgWidth, kImgHeight, kChannelIndex, kChannelCount,
  99       kStrideSlack);
 100
 101   // Destination will be a single channel image with stide matching width.
 102   const int dest_row_stride = kImgWidth;
 103   const int dest_byte_count = dest_row_stride * kImgHeight;
 104   std::vector<unsigned char> intermediate(dest_byte_count);
 105   std::vector<unsigned char> intermediate2(dest_byte_count);
 106   std::vector<unsigned char> control(dest_byte_count);
 107   std::vector<unsigned char> output(dest_byte_count);
 108
 109   const int src_row_stride = ComputeRowStride(
 110       kImgWidth, kChannelCount, kStrideSlack);
 111
 112   const float kernel_sigma = 2.5f;
 113   ConvolutionFilter1D filter;
 114   SetUpGaussianConvolutionKernel(&filter, kernel_sigma, false);
 115   // Process the control image.
 116   SingleChannelConvolveX1D(&input[0], src_row_stride,
 117                            kChannelIndex, kChannelCount,
 118                            filter, image_size,
 119                            &intermediate[0], dest_row_stride, 0, 1, false);
 120   SingleChannelConvolveY1D(&intermediate[0], dest_row_stride, 0, 1,
 121                            filter, image_size,
 122                            &control[0], dest_row_stride, 0, 1, false);
 123
 124   // Now try the same using the other method.
 125   RecursiveFilter recursive_filter(kernel_sigma, RecursiveFilter::FUNCTION);
 126   SingleChannelRecursiveGaussianY(&input[0], src_row_stride,
 127                                   kChannelIndex, kChannelCount,
 128                                   recursive_filter, image_size,
 129                                   &intermediate2[0], dest_row_stride,
 130                                   0, 1, false);
 131   SingleChannelRecursiveGaussianX(&intermediate2[0], dest_row_stride, 0, 1,
 132                                   recursive_filter, image_size,
 133                                   &output[0], dest_row_stride, 0, 1, false);
 134
 135   // We cannot expect the results to be really the same. In particular,
 136   // the standard implementation is computed in completely fixed-point, while
 137   // recursive is done in floating point and squeezed back into char*. On top
 138   // of that, its characteristics are a bit different (consult the paper).
 139   EXPECT_NEAR(std::accumulate(intermediate.begin(), intermediate.end(), 0),
 140               std::accumulate(intermediate2.begin(), intermediate2.end(), 0),
 141               50);
 142   int difference_count = 0;
 143   std::vector<unsigned char>::const_iterator i1, i2;
 144   for (i1 = control.begin(), i2 = output.begin();
 145        i1 != control.end(); ++i1, ++i2) {
 146     if ((*i1 != 0) != (*i2 != 0))
 147       difference_count++;
 148   }
 149
 150   EXPECT_LE(difference_count, 44);  // 44 is 2 * PI * r (r == 7, spot size).
 151 }
 152
 153 TEST(RecursiveGaussian, SmoothingImpulse) {
 154   static const int kImgWidth = 200;
 155   static const int kImgHeight = 300;
 156   static const int kChannelIndex = 3;
 157   static const int kChannelCount = 3;
 158   static const int kStrideSlack = 22;
 159
 160   std::vector<unsigned char> input;
 161   SkISize image_size = SkISize::Make(kImgWidth, kImgHeight);
 162   const SkIPoint centre_point = MakeImpulseImage(
 163       &input, kImgWidth, kImgHeight, kChannelIndex, kChannelCount,
 164       kStrideSlack);
 165
 166   // Destination will be a single channel image with stide matching width.
 167   const int dest_row_stride = kImgWidth;
 168   const int dest_byte_count = dest_row_stride * kImgHeight;
 169   std::vector<unsigned char> intermediate(dest_byte_count);
 170   std::vector<unsigned char> output(dest_byte_count);
 171
 172   const int src_row_stride = ComputeRowStride(
 173       kImgWidth, kChannelCount, kStrideSlack);
 174
 175   const float kernel_sigma = 5.0f;
 176   RecursiveFilter recursive_filter(kernel_sigma, RecursiveFilter::FUNCTION);
 177   SingleChannelRecursiveGaussianY(&input[0], src_row_stride,
 178                                   kChannelIndex, kChannelCount,
 179                                   recursive_filter, image_size,
 180                                   &intermediate[0], dest_row_stride,
 181                                   0, 1, false);
 182   SingleChannelRecursiveGaussianX(&intermediate[0], dest_row_stride, 0, 1,
 183                                   recursive_filter, image_size,
 184                                   &output[0], dest_row_stride, 0, 1, false);
 185
 186   // Check we got the expected impulse response.
 187   const int cx = centre_point.x();
 188   const int cy = centre_point.y();
 189   unsigned char value_x = output[dest_row_stride * cy + cx];
 190   unsigned char value_y = value_x;
 191   EXPECT_GT(value_x, 0);
 192   for (int offset = 0;
 193        offset < std::min(kImgWidth, kImgHeight) && (value_y > 0 || value_x > 0);
 194        ++offset) {
 195     // Symmetricity and monotonicity along X.
 196     EXPECT_EQ(output[dest_row_stride * cy + cx - offset],
 197               output[dest_row_stride * cy + cx + offset]);
 198     EXPECT_LE(output[dest_row_stride * cy + cx - offset], value_x);
 199     value_x = output[dest_row_stride * cy + cx - offset];
 200
 201     // Symmetricity and monotonicity along Y.
 202     EXPECT_EQ(output[dest_row_stride * (cy - offset) + cx],
 203               output[dest_row_stride * (cy + offset) + cx]);
 204     EXPECT_LE(output[dest_row_stride * (cy  - offset) + cx], value_y);
 205     value_y = output[dest_row_stride * (cy - offset) + cx];
 206
 207     // Symmetricity along X/Y (not really assured, but should be close).
 208     EXPECT_NEAR(value_x, value_y, 1);
 209   }
 210
 211   // Smooth the inverse now.
 212   std::vector<unsigned char> output2(dest_byte_count);
 213   std::transform(input.begin(), input.end(), input.begin(),
 214                  std::bind1st(std::minus<unsigned char>(), 255U));
 215   SingleChannelRecursiveGaussianY(&input[0], src_row_stride,
 216                                   kChannelIndex, kChannelCount,
 217                                   recursive_filter, image_size,
 218                                   &intermediate[0], dest_row_stride,
 219                                   0, 1, false);
 220   SingleChannelRecursiveGaussianX(&intermediate[0], dest_row_stride, 0, 1,
 221                                   recursive_filter, image_size,
 222                                   &output2[0], dest_row_stride, 0, 1, false);
 223   // The image should be the reverse of output, but permitting for rounding
 224   // we will only claim that wherever output is 0, output2 should be 255.
 225   // There still can be differences at the edges of the object.
 226   std::vector<unsigned char>::const_iterator i1, i2;
 227   int difference_count = 0;
 228   for (i1 = output.begin(), i2 = output2.begin();
 229        i1 != output.end(); ++i1, ++i2) {
 230     // The line below checks (*i1 == 0 <==> *i2 == 255).
 231     if ((*i1 != 0 && *i2 == 255) && ! (*i1 == 0 && *i2 != 255))
 232       ++difference_count;
 233   }
 234   EXPECT_LE(difference_count, 8);
 235 }
 236
 237 TEST(RecursiveGaussian, FirstDerivative) {
 238   static const int kImgWidth = 512;
 239   static const int kImgHeight = 1024;
 240   static const int kChannelIndex = 2;
 241   static const int kChannelCount = 4;
 242   static const int kStrideSlack = 22;
 243   static const int kBoxSize = 400;
 244
 245   std::vector<unsigned char> input;
 246   const SkISize image_size = SkISize::Make(kImgWidth, kImgHeight);
 247   const SkIRect box =  MakeBoxImage(
 248       &input, kImgWidth, kImgHeight, kChannelIndex, kChannelCount,
 249       kStrideSlack, kBoxSize, kBoxSize, 200);
 250
 251   // Destination will be a single channel image with stide matching width.
 252   const int dest_row_stride = kImgWidth;
 253   const int dest_byte_count = dest_row_stride * kImgHeight;
 254   std::vector<unsigned char> output_x(dest_byte_count);
 255   std::vector<unsigned char> output_y(dest_byte_count);
 256   std::vector<unsigned char> output(dest_byte_count);
 257
 258   const int src_row_stride = ComputeRowStride(
 259       kImgWidth, kChannelCount, kStrideSlack);
 260
 261   const float kernel_sigma = 3.0f;
 262   const int spread = 4 * kernel_sigma;
 263   RecursiveFilter recursive_filter(kernel_sigma,
 264                                    RecursiveFilter::FIRST_DERIVATIVE);
 265   SingleChannelRecursiveGaussianX(&input[0], src_row_stride,
 266                                   kChannelIndex, kChannelCount,
 267                                   recursive_filter, image_size,
 268                                   &output_x[0], dest_row_stride,
 269                                   0, 1, true);
 270   SingleChannelRecursiveGaussianY(&input[0], src_row_stride,
 271                                   kChannelIndex, kChannelCount,
 272                                   recursive_filter, image_size,
 273                                   &output_y[0], dest_row_stride,
 274                                   0, 1, true);
 275
 276   // In test code we can assume adding the two up should do fine.
 277   std::vector<unsigned char>::const_iterator ix, iy;
 278   std::vector<unsigned char>::iterator target;
 279   for (target = output.begin(), ix = output_x.begin(), iy = output_y.begin();
 280        target < output.end(); ++target, ++ix, ++iy) {
 281     *target = *ix + *iy;
 282   }
 283
 284   SkIRect inflated_rect(box);
 285   inflated_rect.outset(spread, spread);
 286   SkIRect deflated_rect(box);
 287   deflated_rect.inset(spread, spread);
 288
 289   int image_total = ComputeBoxSum(output,
 290                                   SkIRect::MakeWH(kImgWidth, kImgHeight),
 291                                   kImgWidth);
 292   int box_inflated = ComputeBoxSum(output, inflated_rect, kImgWidth);
 293   int box_deflated = ComputeBoxSum(output, deflated_rect, kImgWidth);
 294   EXPECT_EQ(box_deflated, 0);
 295   EXPECT_EQ(image_total, box_inflated);
 296
 297   // Try inverted image. Behaviour should be very similar (modulo rounding).
 298   std::transform(input.begin(), input.end(), input.begin(),
 299                  std::bind1st(std::minus<unsigned char>(), 255U));
 300   SingleChannelRecursiveGaussianX(&input[0], src_row_stride,
 301                                   kChannelIndex, kChannelCount,
 302                                   recursive_filter, image_size,
 303                                   &output_x[0], dest_row_stride,
 304                                   0, 1, true);
 305   SingleChannelRecursiveGaussianY(&input[0], src_row_stride,
 306                                   kChannelIndex, kChannelCount,
 307                                   recursive_filter, image_size,
 308                                   &output_y[0], dest_row_stride,
 309                                   0, 1, true);
 310
 311   for (target = output.begin(), ix = output_x.begin(), iy = output_y.begin();
 312        target < output.end(); ++target, ++ix, ++iy) {
 313     *target = *ix + *iy;
 314   }
 315
 316   image_total = ComputeBoxSum(output,
 317                               SkIRect::MakeWH(kImgWidth, kImgHeight),
 318                               kImgWidth);
 319   box_inflated = ComputeBoxSum(output, inflated_rect, kImgWidth);
 320   box_deflated = ComputeBoxSum(output, deflated_rect, kImgWidth);
 321
 322   EXPECT_EQ(box_deflated, 0);
 323   EXPECT_EQ(image_total, box_inflated);
 324 }
 325
 326 TEST(RecursiveGaussian, SecondDerivative) {
 327   static const int kImgWidth = 700;
 328   static const int kImgHeight = 500;
 329   static const int kChannelIndex = 0;
 330   static const int kChannelCount = 2;
 331   static const int kStrideSlack = 42;
 332   static const int kBoxSize = 200;
 333
 334   std::vector<unsigned char> input;
 335   SkISize image_size = SkISize::Make(kImgWidth, kImgHeight);
 336   const SkIRect box = MakeBoxImage(
 337       &input, kImgWidth, kImgHeight, kChannelIndex, kChannelCount,
 338       kStrideSlack, kBoxSize, kBoxSize, 200);
 339
 340   // Destination will be a single channel image with stide matching width.
 341   const int dest_row_stride = kImgWidth;
 342   const int dest_byte_count = dest_row_stride * kImgHeight;
 343   std::vector<unsigned char> output_x(dest_byte_count);
 344   std::vector<unsigned char> output_y(dest_byte_count);
 345   std::vector<unsigned char> output(dest_byte_count);
 346
 347   const int src_row_stride = ComputeRowStride(
 348       kImgWidth, kChannelCount, kStrideSlack);
 349
 350   const float kernel_sigma = 5.0f;
 351   const int spread = 8 * kernel_sigma;
 352   RecursiveFilter recursive_filter(kernel_sigma,
 353                                    RecursiveFilter::SECOND_DERIVATIVE);
 354   SingleChannelRecursiveGaussianX(&input[0], src_row_stride,
 355                                   kChannelIndex, kChannelCount,
 356                                   recursive_filter, image_size,
 357                                   &output_x[0], dest_row_stride,
 358                                   0, 1, true);
 359   SingleChannelRecursiveGaussianY(&input[0], src_row_stride,
 360                                   kChannelIndex, kChannelCount,
 361                                   recursive_filter, image_size,
 362                                   &output_y[0], dest_row_stride,
 363                                   0, 1, true);
 364
 365   // In test code we can assume adding the two up should do fine.
 366   std::vector<unsigned char>::const_iterator ix, iy;
 367   std::vector<unsigned char>::iterator target;
 368   for (target = output.begin(),ix = output_x.begin(), iy = output_y.begin();
 369        target < output.end(); ++target, ++ix, ++iy) {
 370     *target = *ix + *iy;
 371   }
 372
 373   int image_total = ComputeBoxSum(output,
 374                                   SkIRect::MakeWH(kImgWidth, kImgHeight),
 375                                   kImgWidth);
 376   int box_inflated = ComputeBoxSum(output,
 377                                    SkIRect::MakeLTRB(box.left() - spread,
 378                                                      box.top() - spread,
 379                                                      box.right() + spread,
 380                                                      box.bottom() + spread),
 381                                    kImgWidth);
 382   int box_deflated = ComputeBoxSum(output,
 383                                    SkIRect::MakeLTRB(box.left() + spread,
 384                                                      box.top() + spread,
 385                                                      box.right() - spread,
 386                                                      box.bottom() - spread),
 387                                    kImgWidth);
 388   // Since second derivative is not really used and implemented mostly
 389   // for the sake of completeness, we do not verify the detail (that dip
 390   // in the middle). But it is there.
 391   EXPECT_EQ(box_deflated, 0);
 392   EXPECT_EQ(image_total, box_inflated);
 393 }
 394
 395 }  // namespace skia