skia/ext/recursive_gaussian_convolution_unittest.cc - chromium/src - Git at Google

 // Copyright 2013 The Chromium Authors
 // 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/files/file_path.h"
 #include "base/files/file_util.h"
 #include "base/ranges/algorithm.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);
   base::ranges::transform(input, input.begin(),
                           [](unsigned char c) { return 255U - c; });
   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).
   base::ranges::transform(input, input.begin(),
                           [](unsigned char c) { return 255U - c; });
   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
	// Copyright 2013 The Chromium Authors
	// 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/files/file_path.h"
	#include "base/files/file_util.h"
	#include "base/ranges/algorithm.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);
	base::ranges::transform(input, input.begin(),
	[](unsigned char c) { return 255U - c; });
	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).
	base::ranges::transform(input, input.begin(),
	[](unsigned char c) { return 255U - c; });
	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