recursive_gaussian_convolution_unittest.cc

// Copyright (c) 2013 The Chromium Authors. All rights reserved.
// Use of this source code is governed by a BSD-style license that can be
// found in the LICENSE file.

#include <functional>
#include <numeric>
#include <vector>

#include "base/basictypes.h"
#include "base/file_util.h"
#include "base/files/file_path.h"
#include "base/logging.h"
#include "base/time/time.h"
#include "skia/ext/convolver.h"
#include "skia/ext/recursive_gaussian_convolution.h"
#include "testing/gtest/include/gtest/gtest.h"
#include "third_party/skia/include/core/SkPoint.h"
#include "third_party/skia/include/core/SkRect.h"

namespace {

int ComputeRowStride(int width, int channel_count, int stride_slack) {
  return width * channel_count + stride_slack;
}

SkIPoint MakeImpulseImage(std::vector<unsigned char>* image,
                          int width,
                          int height,
                          int channel_index,
                          int channel_count,
                          int stride_slack) {
  const int src_row_stride = ComputeRowStride(
      width, channel_count, stride_slack);
  const int src_byte_count = src_row_stride * height;
  const int signal_x = width / 2;
  const int signal_y = height / 2;

  image->resize(src_byte_count, 0);
  const int non_zero_pixel_index =
      signal_y * src_row_stride + signal_x * channel_count + channel_index;
  (*image)[non_zero_pixel_index] = 255;
  return SkIPoint::Make(signal_x, signal_y);
}

SkIRect MakeBoxImage(std::vector<unsigned char>* image,
                     int width,
                     int height,
                     int channel_index,
                     int channel_count,
                     int stride_slack,
                     int box_width,
                     int box_height,
                     unsigned char value) {
  const int src_row_stride = ComputeRowStride(
      width, channel_count, stride_slack);
  const int src_byte_count = src_row_stride * height;
  const SkIRect box = SkIRect::MakeXYWH((width - box_width) / 2,
                                        (height - box_height) / 2,
                                        box_width, box_height);

  image->resize(src_byte_count, 0);
  for (int y = box.top(); y < box.bottom(); ++y) {
    for (int x = box.left(); x < box.right(); ++x)
      (*image)[y * src_row_stride + x * channel_count + channel_index] = value;
  }

  return box;
}

int ComputeBoxSum(const std::vector<unsigned char>& image,
                  const SkIRect& box,
                  int image_width) {
  // Compute the sum of all pixels in the box. Assume byte stride 1 and row
  // stride same as image_width.
  int sum = 0;
  for (int y = box.top(); y < box.bottom(); ++y) {
    for (int x = box.left(); x < box.right(); ++x)
      sum += image[y * image_width + x];
  }

  return sum;
}

}  // namespace

namespace skia {

TEST(RecursiveGaussian, SmoothingMethodComparison) {
  static const int kImgWidth = 512;
  static const int kImgHeight = 220;
  static const int kChannelIndex = 3;
  static const int kChannelCount = 3;
  static const int kStrideSlack = 22;

  std::vector<unsigned char> input;
  SkISize image_size = SkISize::Make(kImgWidth, kImgHeight);
  MakeImpulseImage(
      &input, kImgWidth, kImgHeight, kChannelIndex, kChannelCount,
      kStrideSlack);

  // Destination will be a single channel image with stide matching width.
  const int dest_row_stride = kImgWidth;
  const int dest_byte_count = dest_row_stride * kImgHeight;
  std::vector<unsigned char> intermediate(dest_byte_count);
  std::vector<unsigned char> intermediate2(dest_byte_count);
  std::vector<unsigned char> control(dest_byte_count);
  std::vector<unsigned char> output(dest_byte_count);

  const int src_row_stride = ComputeRowStride(
      kImgWidth, kChannelCount, kStrideSlack);

  const float kernel_sigma = 2.5f;
  ConvolutionFilter1D filter;
  SetUpGaussianConvolutionKernel(&filter, kernel_sigma, false);
  // Process the control image.
  SingleChannelConvolveX1D(&input[0], src_row_stride,
                           kChannelIndex, kChannelCount,
                           filter, image_size,
                           &intermediate[0], dest_row_stride, 0, 1, false);
  SingleChannelConvolveY1D(&intermediate[0], dest_row_stride, 0, 1,
                           filter, image_size,
                           &control[0], dest_row_stride, 0, 1, false);

  // Now try the same using the other method.
  RecursiveFilter recursive_filter(kernel_sigma, RecursiveFilter::FUNCTION);
  SingleChannelRecursiveGaussianY(&input[0], src_row_stride,
                                  kChannelIndex, kChannelCount,
                                  recursive_filter, image_size,
                                  &intermediate2[0], dest_row_stride,
                                  0, 1, false);
  SingleChannelRecursiveGaussianX(&intermediate2[0], dest_row_stride, 0, 1,
                                  recursive_filter, image_size,
                                  &output[0], dest_row_stride, 0, 1, false);

  // We cannot expect the results to be really the same. In particular,
  // the standard implementation is computed in completely fixed-point, while
  // recursive is done in floating point and squeezed back into char*. On top
  // of that, its characteristics are a bit different (consult the paper).
  EXPECT_NEAR(std::accumulate(intermediate.begin(), intermediate.end(), 0),
              std::accumulate(intermediate2.begin(), intermediate2.end(), 0),
              50);
  int difference_count = 0;
  std::vector<unsigned char>::const_iterator i1, i2;
  for (i1 = control.begin(), i2 = output.begin();
       i1 != control.end(); ++i1, ++i2) {
    if ((*i1 != 0) != (*i2 != 0))
      difference_count++;
  }

  EXPECT_LE(difference_count, 44);  // 44 is 2 * PI * r (r == 7, spot size).
}

TEST(RecursiveGaussian, SmoothingImpulse) {
  static const int kImgWidth = 200;
  static const int kImgHeight = 300;
  static const int kChannelIndex = 3;
  static const int kChannelCount = 3;
  static const int kStrideSlack = 22;

  std::vector<unsigned char> input;
  SkISize image_size = SkISize::Make(kImgWidth, kImgHeight);
  const SkIPoint centre_point = MakeImpulseImage(
      &input, kImgWidth, kImgHeight, kChannelIndex, kChannelCount,
      kStrideSlack);

  // Destination will be a single channel image with stide matching width.
  const int dest_row_stride = kImgWidth;
  const int dest_byte_count = dest_row_stride * kImgHeight;
  std::vector<unsigned char> intermediate(dest_byte_count);
  std::vector<unsigned char> output(dest_byte_count);

  const int src_row_stride = ComputeRowStride(
      kImgWidth, kChannelCount, kStrideSlack);

  const float kernel_sigma = 5.0f;
  RecursiveFilter recursive_filter(kernel_sigma, RecursiveFilter::FUNCTION);
  SingleChannelRecursiveGaussianY(&input[0], src_row_stride,
                                  kChannelIndex, kChannelCount,
                                  recursive_filter, image_size,
                                  &intermediate[0], dest_row_stride,
                                  0, 1, false);
  SingleChannelRecursiveGaussianX(&intermediate[0], dest_row_stride, 0, 1,
                                  recursive_filter, image_size,
                                  &output[0], dest_row_stride, 0, 1, false);

  // Check we got the expected impulse response.
  const int cx = centre_point.x();
  const int cy = centre_point.y();
  unsigned char value_x = output[dest_row_stride * cy + cx];
  unsigned char value_y = value_x;
  EXPECT_GT(value_x, 0);
  for (int offset = 0;
       offset < std::min(kImgWidth, kImgHeight) && (value_y > 0 || value_x > 0);
       ++offset) {
    // Symmetricity and monotonicity along X.
    EXPECT_EQ(output[dest_row_stride * cy + cx - offset],
              output[dest_row_stride * cy + cx + offset]);
    EXPECT_LE(output[dest_row_stride * cy + cx - offset], value_x);
    value_x = output[dest_row_stride * cy + cx - offset];

    // Symmetricity and monotonicity along Y.
    EXPECT_EQ(output[dest_row_stride * (cy - offset) + cx],
              output[dest_row_stride * (cy + offset) + cx]);
    EXPECT_LE(output[dest_row_stride * (cy  - offset) + cx], value_y);
    value_y = output[dest_row_stride * (cy - offset) + cx];

    // Symmetricity along X/Y (not really assured, but should be close).
    EXPECT_NEAR(value_x, value_y, 1);
  }

  // Smooth the inverse now.
  std::vector<unsigned char> output2(dest_byte_count);
  std::transform(input.begin(), input.end(), input.begin(),
                 std::bind1st(std::minus<unsigned char>(), 255U));
  SingleChannelRecursiveGaussianY(&input[0], src_row_stride,
                                  kChannelIndex, kChannelCount,
                                  recursive_filter, image_size,
                                  &intermediate[0], dest_row_stride,
                                  0, 1, false);
  SingleChannelRecursiveGaussianX(&intermediate[0], dest_row_stride, 0, 1,
                                  recursive_filter, image_size,
                                  &output2[0], dest_row_stride, 0, 1, false);
  // The image should be the reverse of output, but permitting for rounding
  // we will only claim that wherever output is 0, output2 should be 255.
  // There still can be differences at the edges of the object.
  std::vector<unsigned char>::const_iterator i1, i2;
  int difference_count = 0;
  for (i1 = output.begin(), i2 = output2.begin();
       i1 != output.end(); ++i1, ++i2) {
    // The line below checks (*i1 == 0 <==> *i2 == 255).
    if ((*i1 != 0 && *i2 == 255) && ! (*i1 == 0 && *i2 != 255))
      ++difference_count;
  }
  EXPECT_LE(difference_count, 8);
}

TEST(RecursiveGaussian, FirstDerivative) {
  static const int kImgWidth = 512;
  static const int kImgHeight = 1024;
  static const int kChannelIndex = 2;
  static const int kChannelCount = 4;
  static const int kStrideSlack = 22;
  static const int kBoxSize = 400;

  std::vector<unsigned char> input;
  const SkISize image_size = SkISize::Make(kImgWidth, kImgHeight);
  const SkIRect box =  MakeBoxImage(
      &input, kImgWidth, kImgHeight, kChannelIndex, kChannelCount,
      kStrideSlack, kBoxSize, kBoxSize, 200);

  // Destination will be a single channel image with stide matching width.
  const int dest_row_stride = kImgWidth;
  const int dest_byte_count = dest_row_stride * kImgHeight;
  std::vector<unsigned char> output_x(dest_byte_count);
  std::vector<unsigned char> output_y(dest_byte_count);
  std::vector<unsigned char> output(dest_byte_count);

  const int src_row_stride = ComputeRowStride(
      kImgWidth, kChannelCount, kStrideSlack);

  const float kernel_sigma = 3.0f;
  const int spread = 4 * kernel_sigma;
  RecursiveFilter recursive_filter(kernel_sigma,
                                   RecursiveFilter::FIRST_DERIVATIVE);
  SingleChannelRecursiveGaussianX(&input[0], src_row_stride,
                                  kChannelIndex, kChannelCount,
                                  recursive_filter, image_size,
                                  &output_x[0], dest_row_stride,
                                  0, 1, true);
  SingleChannelRecursiveGaussianY(&input[0], src_row_stride,
                                  kChannelIndex, kChannelCount,
                                  recursive_filter, image_size,
                                  &output_y[0], dest_row_stride,
                                  0, 1, true);

  // In test code we can assume adding the two up should do fine.
  std::vector<unsigned char>::const_iterator ix, iy;
  std::vector<unsigned char>::iterator target;
  for (target = output.begin(), ix = output_x.begin(), iy = output_y.begin();
       target < output.end(); ++target, ++ix, ++iy) {
    *target = *ix + *iy;
  }

  SkIRect inflated_rect(box);
  inflated_rect.outset(spread, spread);
  SkIRect deflated_rect(box);
  deflated_rect.inset(spread, spread);

  int image_total = ComputeBoxSum(output,
                                  SkIRect::MakeWH(kImgWidth, kImgHeight),
                                  kImgWidth);
  int box_inflated = ComputeBoxSum(output, inflated_rect, kImgWidth);
  int box_deflated = ComputeBoxSum(output, deflated_rect, kImgWidth);
  EXPECT_EQ(box_deflated, 0);
  EXPECT_EQ(image_total, box_inflated);

  // Try inverted image. Behaviour should be very similar (modulo rounding).
  std::transform(input.begin(), input.end(), input.begin(),
                 std::bind1st(std::minus<unsigned char>(), 255U));
  SingleChannelRecursiveGaussianX(&input[0], src_row_stride,
                                  kChannelIndex, kChannelCount,
                                  recursive_filter, image_size,
                                  &output_x[0], dest_row_stride,
                                  0, 1, true);
  SingleChannelRecursiveGaussianY(&input[0], src_row_stride,
                                  kChannelIndex, kChannelCount,
                                  recursive_filter, image_size,
                                  &output_y[0], dest_row_stride,
                                  0, 1, true);

  for (target = output.begin(), ix = output_x.begin(), iy = output_y.begin();
       target < output.end(); ++target, ++ix, ++iy) {
    *target = *ix + *iy;
  }

  image_total = ComputeBoxSum(output,
                              SkIRect::MakeWH(kImgWidth, kImgHeight),
                              kImgWidth);
  box_inflated = ComputeBoxSum(output, inflated_rect, kImgWidth);
  box_deflated = ComputeBoxSum(output, deflated_rect, kImgWidth);

  EXPECT_EQ(box_deflated, 0);
  EXPECT_EQ(image_total, box_inflated);
}

TEST(RecursiveGaussian, SecondDerivative) {
  static const int kImgWidth = 700;
  static const int kImgHeight = 500;
  static const int kChannelIndex = 0;
  static const int kChannelCount = 2;
  static const int kStrideSlack = 42;
  static const int kBoxSize = 200;

  std::vector<unsigned char> input;
  SkISize image_size = SkISize::Make(kImgWidth, kImgHeight);
  const SkIRect box = MakeBoxImage(
      &input, kImgWidth, kImgHeight, kChannelIndex, kChannelCount,
      kStrideSlack, kBoxSize, kBoxSize, 200);

  // Destination will be a single channel image with stide matching width.
  const int dest_row_stride = kImgWidth;
  const int dest_byte_count = dest_row_stride * kImgHeight;
  std::vector<unsigned char> output_x(dest_byte_count);
  std::vector<unsigned char> output_y(dest_byte_count);
  std::vector<unsigned char> output(dest_byte_count);

  const int src_row_stride = ComputeRowStride(
      kImgWidth, kChannelCount, kStrideSlack);

  const float kernel_sigma = 5.0f;
  const int spread = 8 * kernel_sigma;
  RecursiveFilter recursive_filter(kernel_sigma,
                                   RecursiveFilter::SECOND_DERIVATIVE);
  SingleChannelRecursiveGaussianX(&input[0], src_row_stride,
                                  kChannelIndex, kChannelCount,
                                  recursive_filter, image_size,
                                  &output_x[0], dest_row_stride,
                                  0, 1, true);
  SingleChannelRecursiveGaussianY(&input[0], src_row_stride,
                                  kChannelIndex, kChannelCount,
                                  recursive_filter, image_size,
                                  &output_y[0], dest_row_stride,
                                  0, 1, true);

  // In test code we can assume adding the two up should do fine.
  std::vector<unsigned char>::const_iterator ix, iy;
  std::vector<unsigned char>::iterator target;
  for (target = output.begin(),ix = output_x.begin(), iy = output_y.begin();
       target < output.end(); ++target, ++ix, ++iy) {
    *target = *ix + *iy;
  }

  int image_total = ComputeBoxSum(output,
                                  SkIRect::MakeWH(kImgWidth, kImgHeight),
                                  kImgWidth);
  int box_inflated = ComputeBoxSum(output,
                                   SkIRect::MakeLTRB(box.left() - spread,
                                                     box.top() - spread,
                                                     box.right() + spread,
                                                     box.bottom() + spread),
                                   kImgWidth);
  int box_deflated = ComputeBoxSum(output,
                                   SkIRect::MakeLTRB(box.left() + spread,
                                                     box.top() + spread,
                                                     box.right() - spread,
                                                     box.bottom() - spread),
                                   kImgWidth);
  // Since second derivative is not really used and implemented mostly
  // for the sake of completeness, we do not verify the detail (that dip
  // in the middle). But it is there.
  EXPECT_EQ(box_deflated, 0);
  EXPECT_EQ(image_total, box_inflated);
}

}  // namespace skia