Skip to content

Commit

Permalink
An implementation of recursive gaussian filters.
Browse files Browse the repository at this point in the history 
Based on the Young and Vliet paper 'Recursive implementation of the Gaussian filter'.

BUG=155269

Review URL: https://chromiumcodereview.appspot.com/14740020

git-svn-id: svn://svn.chromium.org/chrome/trunk/src@199358 0039d316-1c4b-4281-b951-d872f2087c98
  • Loading branch information
motek@chromium.org committed May 10, 2013
1 parent f4cb55d commit e8a35d1
Show file tree
Hide file tree
Showing 9 changed files with 888 additions and 73 deletions.
173 changes: 103 additions & 70 deletions chrome/browser/thumbnails/content_analysis.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,12 +14,15 @@

#include "base/logging.h"
#include "skia/ext/convolver.h"
#include "skia/ext/recursive_gaussian_convolution.h"
#include "third_party/skia/include/core/SkBitmap.h"
#include "third_party/skia/include/core/SkSize.h"
#include "ui/gfx/color_analysis.h"

namespace {

const float kSigmaThresholdForRecursive = 1.5f;

template<class InputIterator, class OutputIterator, class Compare>
void SlidingWindowMinMax(InputIterator first,
InputIterator last,
Expand Down Expand Up @@ -75,37 +78,6 @@ void ApplyGaussianGradientMagnitudeFilter(SkBitmap* input_bitmap,
DCHECK(input_bitmap->getPixels());
DCHECK_EQ(SkBitmap::kA8_Config, input_bitmap->config());

const int tail_length = static_cast<int>(4.0f * kernel_sigma + 0.5f);
const int kernel_size = tail_length * 2 + 1;
const float sigmasq = kernel_sigma * kernel_sigma;
std::vector<float> smoother_weights(kernel_size, 0.0);
float kernel_sum = 1.0f;

smoother_weights[tail_length] = 1.0f;

for (int ii = 1; ii <= tail_length; ++ii) {
float v = std::exp(-0.5f * ii * ii / sigmasq);
smoother_weights[tail_length + ii] = v;
smoother_weights[tail_length - ii] = v;
kernel_sum += 2.0f * v;
}

for (int i = 0; i < kernel_size; ++i)
smoother_weights[i] /= kernel_sum;

std::vector<float> gradient_weights(kernel_size, 0.0);
gradient_weights[tail_length] = 0.0;
for (int ii = 1; ii <= tail_length; ++ii) {
float v = sigmasq * smoother_weights[tail_length + ii] / ii;
gradient_weights[tail_length + ii] = v;
gradient_weights[tail_length - ii] = -v;
}

skia::ConvolutionFilter1D smoothing_filter;
skia::ConvolutionFilter1D gradient_filter;
smoothing_filter.AddFilter(0, &smoother_weights[0], smoother_weights.size());
gradient_filter.AddFilter(0, &gradient_weights[0], gradient_weights.size());

// To perform computations we will need one intermediate buffer. It can
// very well be just another bitmap.
const SkISize image_size = SkISize::Make(input_bitmap->width(),
Expand All @@ -115,49 +87,110 @@ void ApplyGaussianGradientMagnitudeFilter(SkBitmap* input_bitmap,
input_bitmap->config(), image_size.width(), image_size.height());
intermediate.allocPixels();

skia::SingleChannelConvolveX1D(
input_bitmap->getAddr8(0, 0),
static_cast<int>(input_bitmap->rowBytes()),
0, input_bitmap->bytesPerPixel(),
smoothing_filter,
image_size,
intermediate.getAddr8(0, 0),
static_cast<int>(intermediate.rowBytes()),
0, intermediate.bytesPerPixel(), false);
skia::SingleChannelConvolveY1D(
intermediate.getAddr8(0, 0),
static_cast<int>(intermediate.rowBytes()),
0, intermediate.bytesPerPixel(),
smoothing_filter,
image_size,
input_bitmap->getAddr8(0, 0),
static_cast<int>(input_bitmap->rowBytes()),
0, input_bitmap->bytesPerPixel(), false);

// Now the gradient operator (we will need two buffers).
SkBitmap intermediate2;
intermediate2.setConfig(
input_bitmap->config(), image_size.width(), image_size.height());
intermediate2.allocPixels();

skia::SingleChannelConvolveX1D(
input_bitmap->getAddr8(0, 0),
static_cast<int>(input_bitmap->rowBytes()),
0, input_bitmap->bytesPerPixel(),
gradient_filter,
image_size,
intermediate.getAddr8(0, 0),
static_cast<int>(intermediate.rowBytes()),
0, intermediate.bytesPerPixel(), true);
skia::SingleChannelConvolveY1D(
input_bitmap->getAddr8(0, 0),
static_cast<int>(input_bitmap->rowBytes()),
0, input_bitmap->bytesPerPixel(),
gradient_filter,
image_size,
intermediate2.getAddr8(0, 0),
static_cast<int>(intermediate2.rowBytes()),
0, intermediate2.bytesPerPixel(), true);

if (kernel_sigma <= kSigmaThresholdForRecursive) {
// For small kernels classic implementation is faster.
skia::ConvolutionFilter1D smoothing_filter;
skia::SetUpGaussianConvolutionKernel(
&smoothing_filter, kernel_sigma, false);
skia::SingleChannelConvolveX1D(
input_bitmap->getAddr8(0, 0),
static_cast<int>(input_bitmap->rowBytes()),
0, input_bitmap->bytesPerPixel(),
smoothing_filter,
image_size,
intermediate.getAddr8(0, 0),
static_cast<int>(intermediate.rowBytes()),
0, intermediate.bytesPerPixel(), false);
skia::SingleChannelConvolveY1D(
intermediate.getAddr8(0, 0),
static_cast<int>(intermediate.rowBytes()),
0, intermediate.bytesPerPixel(),
smoothing_filter,
image_size,
input_bitmap->getAddr8(0, 0),
static_cast<int>(input_bitmap->rowBytes()),
0, input_bitmap->bytesPerPixel(), false);

skia::ConvolutionFilter1D gradient_filter;
skia::SetUpGaussianConvolutionKernel(&gradient_filter, kernel_sigma, true);
skia::SingleChannelConvolveX1D(
input_bitmap->getAddr8(0, 0),
static_cast<int>(input_bitmap->rowBytes()),
0, input_bitmap->bytesPerPixel(),
gradient_filter,
image_size,
intermediate.getAddr8(0, 0),
static_cast<int>(intermediate.rowBytes()),
0, intermediate.bytesPerPixel(), true);
skia::SingleChannelConvolveY1D(
input_bitmap->getAddr8(0, 0),
static_cast<int>(input_bitmap->rowBytes()),
0, input_bitmap->bytesPerPixel(),
gradient_filter,
image_size,
intermediate2.getAddr8(0, 0),
static_cast<int>(intermediate2.rowBytes()),
0, intermediate2.bytesPerPixel(), true);
} else {
// For larger sigma values use the recursive filter.
skia::RecursiveFilter smoothing_filter(kernel_sigma,
skia::RecursiveFilter::FUNCTION);
skia::SingleChannelRecursiveGaussianX(
input_bitmap->getAddr8(0, 0),
static_cast<int>(input_bitmap->rowBytes()),
0, input_bitmap->bytesPerPixel(),
smoothing_filter,
image_size,
intermediate.getAddr8(0, 0),
static_cast<int>(intermediate.rowBytes()),
0, intermediate.bytesPerPixel(), false);
unsigned char smoothed_max = skia::SingleChannelRecursiveGaussianY(
intermediate.getAddr8(0, 0),
static_cast<int>(intermediate.rowBytes()),
0, intermediate.bytesPerPixel(),
smoothing_filter,
image_size,
input_bitmap->getAddr8(0, 0),
static_cast<int>(input_bitmap->rowBytes()),
0, input_bitmap->bytesPerPixel(), false);
if (smoothed_max < 127) {
int bit_shift = 8 - static_cast<int>(
std::log10(static_cast<float>(smoothed_max)) / std::log10(2.0f));
for (int r = 0; r < image_size.height(); ++r) {
uint8* row = input_bitmap->getAddr8(0, r);
for (int c = 0; c < image_size.width(); ++c, ++row) {
*row <<= bit_shift;
}
}
}

skia::RecursiveFilter gradient_filter(
kernel_sigma, skia::RecursiveFilter::FIRST_DERIVATIVE);
skia::SingleChannelRecursiveGaussianX(
input_bitmap->getAddr8(0, 0),
static_cast<int>(input_bitmap->rowBytes()),
0, input_bitmap->bytesPerPixel(),
gradient_filter,
image_size,
intermediate.getAddr8(0, 0),
static_cast<int>(intermediate.rowBytes()),
0, intermediate.bytesPerPixel(), true);
skia::SingleChannelRecursiveGaussianY(
input_bitmap->getAddr8(0, 0),
static_cast<int>(input_bitmap->rowBytes()),
0, input_bitmap->bytesPerPixel(),
gradient_filter,
image_size,
intermediate2.getAddr8(0, 0),
static_cast<int>(intermediate2.rowBytes()),
0, intermediate2.bytesPerPixel(), true);
}

unsigned grad_max = 0;
for (int r = 0; r < image_size.height(); ++r) {
Expand All @@ -175,7 +208,7 @@ void ApplyGaussianGradientMagnitudeFilter(SkBitmap* input_bitmap,
bit_shift = static_cast<int>(
std::log10(static_cast<float>(grad_max)) / std::log10(2.0f)) - 7;
for (int r = 0; r < image_size.height(); ++r) {
const uint8* grad_x_row =intermediate.getAddr8(0, r);
const uint8* grad_x_row = intermediate.getAddr8(0, r);
const uint8* grad_y_row = intermediate2.getAddr8(0, r);
uint8* target_row = input_bitmap->getAddr8(0, r);
for (int c = 0; c < image_size.width(); ++c) {
Expand Down
1 change: 1 addition & 0 deletions chrome/chrome_tests_unit.gypi
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1767,6 +1767,7 @@
'../skia/ext/convolver_unittest.cc',
'../skia/ext/image_operations_unittest.cc',
'../skia/ext/platform_canvas_unittest.cc',
'../skia/ext/recursive_gaussian_convolution_unittest.cc',
'../skia/ext/refptr_unittest.cc',
'../skia/ext/skia_utils_ios_unittest.mm',
'../skia/ext/skia_utils_mac_unittest.mm',
Expand Down
35 changes: 35 additions & 0 deletions skia/ext/convolver.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -672,4 +672,39 @@ void SingleChannelConvolveY1D(const unsigned char* source_data,
}
}

void SetUpGaussianConvolutionKernel(ConvolutionFilter1D* filter,
float kernel_sigma,
bool derivative) {
DCHECK(filter != NULL);
DCHECK_GT(kernel_sigma, 0.0);
const int tail_length = static_cast<int>(4.0f * kernel_sigma + 0.5f);
const int kernel_size = tail_length * 2 + 1;
const float sigmasq = kernel_sigma * kernel_sigma;
std::vector<float> kernel_weights(kernel_size, 0.0);
float kernel_sum = 1.0f;

kernel_weights[tail_length] = 1.0f;

for (int ii = 1; ii <= tail_length; ++ii) {
float v = std::exp(-0.5f * ii * ii / sigmasq);
kernel_weights[tail_length + ii] = v;
kernel_weights[tail_length - ii] = v;
kernel_sum += 2.0f * v;
}

for (int i = 0; i < kernel_size; ++i)
kernel_weights[i] /= kernel_sum;

if (derivative) {
kernel_weights[tail_length] = 0.0;
for (int ii = 1; ii <= tail_length; ++ii) {
float v = sigmasq * kernel_weights[tail_length + ii] / ii;
kernel_weights[tail_length + ii] = v;
kernel_weights[tail_length - ii] = -v;
}
}

filter->AddFilter(0, &kernel_weights[0], kernel_weights.size());
}

} // namespace skia
13 changes: 10 additions & 3 deletions skia/ext/convolver.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,9 +115,9 @@ class ConvolutionFilter1D {
// There will be |filter_length| values in the return array.
// Returns NULL if the filter is 0-length (for instance when all floating
// point values passed to AddFilter were clipped to 0).
const Fixed* GetSingleFilter(int* specified_filter_length,
int* filter_offset,
int* filter_length) const;
SK_API const Fixed* GetSingleFilter(int* specified_filter_length,
int* filter_offset,
int* filter_length) const;

inline void PaddingForSIMD() {
// Padding |padding_count| of more dummy coefficients after the coefficients
Expand Down Expand Up @@ -219,6 +219,13 @@ SK_API void SingleChannelConvolveY1D(const unsigned char* source_data,
int output_channel_count,
bool absolute_values);

// Set up the |filter| instance with a gaussian kernel. |kernel_sigma| is the
// parameter of gaussian. If |derivative| is true, the kernel will be that of
// the first derivative. Intended for use with the two routines above.
SK_API void SetUpGaussianConvolutionKernel(ConvolutionFilter1D* filter,
float kernel_sigma,
bool derivative);

} // namespace skia

#endif // SKIA_EXT_CONVOLVER_H_
53 changes: 53 additions & 0 deletions skia/ext/convolver_unittest.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,8 @@

#include <string.h>
#include <time.h>
#include <algorithm>
#include <numeric>
#include <vector>

#include "base/basictypes.h"
Expand Down Expand Up @@ -473,4 +475,55 @@ TEST(Convolver, SeparableSingleConvolutionEdges) {
EXPECT_NEAR(output[test_column + dest_row_stride * (kImgHeight - 4)], 40, 1);
}

TEST(Convolver, SetUpGaussianConvolutionFilter) {
ConvolutionFilter1D smoothing_filter;
ConvolutionFilter1D gradient_filter;
SetUpGaussianConvolutionKernel(&smoothing_filter, 4.5f, false);
SetUpGaussianConvolutionKernel(&gradient_filter, 3.0f, true);

int specified_filter_length;
int filter_offset;
int filter_length;

const ConvolutionFilter1D::Fixed* smoothing_kernel =
smoothing_filter.GetSingleFilter(
&specified_filter_length, &filter_offset, &filter_length);
EXPECT_TRUE(smoothing_kernel);
std::vector<float> fp_smoothing_kernel(filter_length);
std::transform(smoothing_kernel,
smoothing_kernel + filter_length,
fp_smoothing_kernel.begin(),
ConvolutionFilter1D::FixedToFloat);
// Should sum-up to 1 (nearly), and all values whould be in ]0, 1[.
EXPECT_NEAR(std::accumulate(
fp_smoothing_kernel.begin(), fp_smoothing_kernel.end(), 0.0f),
1.0f, 0.01f);
EXPECT_GT(*std::min_element(fp_smoothing_kernel.begin(),
fp_smoothing_kernel.end()), 0.0f);
EXPECT_LT(*std::max_element(fp_smoothing_kernel.begin(),
fp_smoothing_kernel.end()), 1.0f);

const ConvolutionFilter1D::Fixed* gradient_kernel =
gradient_filter.GetSingleFilter(
&specified_filter_length, &filter_offset, &filter_length);
EXPECT_TRUE(gradient_kernel);
std::vector<float> fp_gradient_kernel(filter_length);
std::transform(gradient_kernel,
gradient_kernel + filter_length,
fp_gradient_kernel.begin(),
ConvolutionFilter1D::FixedToFloat);
// Should sum-up to 0, and all values whould be in ]-1.5, 1.5[.
EXPECT_NEAR(std::accumulate(
fp_gradient_kernel.begin(), fp_gradient_kernel.end(), 0.0f),
0.0f, 0.01f);
EXPECT_GT(*std::min_element(fp_gradient_kernel.begin(),
fp_gradient_kernel.end()), -1.5f);
EXPECT_LT(*std::min_element(fp_gradient_kernel.begin(),
fp_gradient_kernel.end()), 0.0f);
EXPECT_LT(*std::max_element(fp_gradient_kernel.begin(),
fp_gradient_kernel.end()), 1.5f);
EXPECT_GT(*std::max_element(fp_gradient_kernel.begin(),
fp_gradient_kernel.end()), 0.0f);
}

} // namespace skia
Loading

0 comments on commit e8a35d1

Please sign in to comment.