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
20namespace {
21
22int ComputeRowStride(int width, int channel_count, int stride_slack) {
23  return width * channel_count + stride_slack;
24}
25
26SkIPoint 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
45SkIRect 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
70int 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
86namespace skia {
87
88TEST(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
153TEST(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
237TEST(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
326TEST(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
396