Blender  V3.3
convolve.cc
Go to the documentation of this file.
1 // Copyright (c) 2007, 2008 libmv authors.
2 //
3 // Permission is hereby granted, free of charge, to any person obtaining a copy
4 // of this software and associated documentation files (the "Software"), to
5 // deal in the Software without restriction, including without limitation the
6 // rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
7 // sell copies of the Software, and to permit persons to whom the Software is
8 // furnished to do so, subject to the following conditions:
9 //
10 // The above copyright notice and this permission notice shall be included in
11 // all copies or substantial portions of the Software.
12 //
13 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
14 // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
15 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
16 // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
17 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
18 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
19 // IN THE SOFTWARE.
20 
21 #include "libmv/image/convolve.h"
22 
23 #include <cmath>
24 
25 #include "libmv/image/image.h"
26 
27 namespace libmv {
28 
29 // Compute a Gaussian kernel and derivative, such that you can take the
30 // derivative of an image by convolving with the kernel horizontally then the
31 // derivative vertically to get (eg) the y derivative.
32 void ComputeGaussianKernel(double sigma, Vec* kernel, Vec* derivative) {
33  assert(sigma >= 0.0);
34 
35  // 0.004 implies a 3 pixel kernel with 1 pixel sigma.
36  const float truncation_factor = 0.004f;
37 
38  // Calculate the kernel size based on sigma such that it is odd.
39  float precisehalfwidth = GaussianInversePositive(truncation_factor, sigma);
40  int width = lround(2 * precisehalfwidth);
41  if (width % 2 == 0) {
42  width++;
43  }
44  // Calculate the gaussian kernel and its derivative.
45  kernel->resize(width);
46  derivative->resize(width);
47  kernel->setZero();
48  derivative->setZero();
49  int halfwidth = width / 2;
50  for (int i = -halfwidth; i <= halfwidth; ++i) {
51  (*kernel)(i + halfwidth) = Gaussian(i, sigma);
52  (*derivative)(i + halfwidth) = GaussianDerivative(i, sigma);
53  }
54  // Since images should not get brighter or darker, normalize.
56 
57  // Normalize the derivative differently. See
58  // www.cs.duke.edu/courses/spring03/cps296.1/handouts/Image%20Processing.pdf
59  double factor = 0.;
60  for (int i = -halfwidth; i <= halfwidth; ++i) {
61  factor -= i * (*derivative)(i + halfwidth);
62  }
63  *derivative /= factor;
64 }
65 
66 template <int size, bool vertical>
67 void FastConvolve(const Vec& kernel,
68  int width,
69  int height,
70  const float* src,
71  int src_stride,
72  int src_line_stride,
73  float* dst,
74  int dst_stride) {
75  double coefficients[2 * size + 1];
76  for (int k = 0; k < 2 * size + 1; ++k) {
77  coefficients[k] = kernel(2 * size - k);
78  }
79  // Fast path: if the kernel has a certain size, use the constant sized loops.
80  for (int y = 0; y < height; ++y) {
81  for (int x = 0; x < width; ++x) {
82  double sum = 0;
83  for (int k = -size; k <= size; ++k) {
84  if (vertical) {
85  if (y + k >= 0 && y + k < height) {
86  sum += src[k * src_line_stride] * coefficients[k + size];
87  }
88  } else {
89  if (x + k >= 0 && x + k < width) {
90  sum += src[k * src_stride] * coefficients[k + size];
91  }
92  }
93  }
94  dst[0] = static_cast<float>(sum);
95  src += src_stride;
96  dst += dst_stride;
97  }
98  }
99 }
100 
101 template <bool vertical>
102 void Convolve(const Array3Df& in,
103  const Vec& kernel,
104  Array3Df* out_pointer,
105  int plane) {
106  int width = in.Width();
107  int height = in.Height();
108  Array3Df& out = *out_pointer;
109  if (plane == -1) {
110  out.ResizeLike(in);
111  plane = 0;
112  }
113 
114  assert(kernel.size() % 2 == 1);
115  assert(&in != out_pointer);
116 
117  int src_line_stride = in.Stride(0);
118  int src_stride = in.Stride(1);
119  int dst_stride = out.Stride(1);
120  const float* src = in.Data();
121  float* dst = out.Data() + plane;
122 
123  // Use a dispatch table to make most convolutions used in practice use the
124  // fast path.
125  int half_width = kernel.size() / 2;
126  switch (half_width) {
127 #define static_convolution(size) \
128  case size: \
129  FastConvolve<size, vertical>(kernel, \
130  width, \
131  height, \
132  src, \
133  src_stride, \
134  src_line_stride, \
135  dst, \
136  dst_stride); \
137  break;
141 #undef static_convolution
142  default : int dynamic_size = kernel.size() / 2;
143  for (int y = 0; y < height; ++y) {
144  for (int x = 0; x < width; ++x) {
145  double sum = 0;
146  // Slow path: this loop cannot be unrolled.
147  for (int k = -dynamic_size; k <= dynamic_size; ++k) {
148  if (vertical) {
149  if (y + k >= 0 && y + k < height) {
150  sum += src[k * src_line_stride] *
151  kernel(2 * dynamic_size - (k + dynamic_size));
152  }
153  } else {
154  if (x + k >= 0 && x + k < width) {
155  sum += src[k * src_stride] *
156  kernel(2 * dynamic_size - (k + dynamic_size));
157  }
158  }
159  }
160  dst[0] = static_cast<float>(sum);
161  src += src_stride;
162  dst += dst_stride;
163  }
164  }
165  }
166 }
167 
169  const Vec& kernel,
170  Array3Df* out_pointer,
171  int plane) {
172  Convolve<false>(in, kernel, out_pointer, plane);
173 }
174 
175 void ConvolveVertical(const Array3Df& in,
176  const Vec& kernel,
177  Array3Df* out_pointer,
178  int plane) {
179  Convolve<true>(in, kernel, out_pointer, plane);
180 }
181 
182 void ConvolveGaussian(const Array3Df& in, double sigma, Array3Df* out_pointer) {
183  Vec kernel, derivative;
184  ComputeGaussianKernel(sigma, &kernel, &derivative);
185 
186  Array3Df tmp;
187  ConvolveVertical(in, kernel, &tmp);
188  ConvolveHorizontal(tmp, kernel, out_pointer);
189 }
190 
191 void ImageDerivatives(const Array3Df& in,
192  double sigma,
193  Array3Df* gradient_x,
194  Array3Df* gradient_y) {
195  Vec kernel, derivative;
196  ComputeGaussianKernel(sigma, &kernel, &derivative);
197  Array3Df tmp;
198 
199  // Compute first derivative in x.
200  ConvolveVertical(in, kernel, &tmp);
201  ConvolveHorizontal(tmp, derivative, gradient_x);
202 
203  // Compute first derivative in y.
204  ConvolveHorizontal(in, kernel, &tmp);
205  ConvolveVertical(tmp, derivative, gradient_y);
206 }
207 
209  double sigma,
210  Array3Df* blurred_image,
211  Array3Df* gradient_x,
212  Array3Df* gradient_y) {
213  Vec kernel, derivative;
214  ComputeGaussianKernel(sigma, &kernel, &derivative);
215  Array3Df tmp;
216 
217  // Compute convolved image.
218  ConvolveVertical(in, kernel, &tmp);
219  ConvolveHorizontal(tmp, kernel, blurred_image);
220 
221  // Compute first derivative in x (reusing vertical convolution above).
222  ConvolveHorizontal(tmp, derivative, gradient_x);
223 
224  // Compute first derivative in y.
225  ConvolveHorizontal(in, kernel, &tmp);
226  ConvolveVertical(tmp, derivative, gradient_y);
227 }
228 
229 // Compute the gaussian blur of an image and the derivatives of the blurred
230 // image, and store the results in three channels. Since the blurred value and
231 // gradients are closer in memory, this leads to better performance if all
232 // three values are needed at the same time.
234  double sigma,
235  Array3Df* blurred_and_gradxy) {
236  assert(in.Depth() == 1);
237 
238  Vec kernel, derivative;
239  ComputeGaussianKernel(sigma, &kernel, &derivative);
240 
241  // Compute convolved image.
242  Array3Df tmp;
243  ConvolveVertical(in, kernel, &tmp);
244  blurred_and_gradxy->Resize(in.Height(), in.Width(), 3);
245  ConvolveHorizontal(tmp, kernel, blurred_and_gradxy, 0);
246 
247  // Compute first derivative in x.
248  ConvolveHorizontal(tmp, derivative, blurred_and_gradxy, 1);
249 
250  // Compute first derivative in y.
251  ConvolveHorizontal(in, kernel, &tmp);
252  ConvolveVertical(tmp, derivative, blurred_and_gradxy, 2);
253 }
254 
256  int window_size,
257  Array3Df* out_pointer) {
258  Array3Df& out = *out_pointer;
259  out.ResizeLike(in);
260  int half_width = (window_size - 1) / 2;
261 
262  for (int k = 0; k < in.Depth(); ++k) {
263  for (int i = 0; i < in.Height(); ++i) {
264  float sum = 0;
265  // Init sum.
266  for (int j = 0; j < half_width; ++j) {
267  sum += in(i, j, k);
268  }
269  // Fill left border.
270  for (int j = 0; j < half_width + 1; ++j) {
271  sum += in(i, j + half_width, k);
272  out(i, j, k) = sum;
273  }
274  // Fill interior.
275  for (int j = half_width + 1; j < in.Width() - half_width; ++j) {
276  sum -= in(i, j - half_width - 1, k);
277  sum += in(i, j + half_width, k);
278  out(i, j, k) = sum;
279  }
280  // Fill right border.
281  for (int j = in.Width() - half_width; j < in.Width(); ++j) {
282  sum -= in(i, j - half_width - 1, k);
283  out(i, j, k) = sum;
284  }
285  }
286  }
287 }
288 
290  int window_size,
291  Array3Df* out_pointer) {
292  Array3Df& out = *out_pointer;
293  out.ResizeLike(in);
294  int half_width = (window_size - 1) / 2;
295 
296  for (int k = 0; k < in.Depth(); ++k) {
297  for (int j = 0; j < in.Width(); ++j) {
298  float sum = 0;
299  // Init sum.
300  for (int i = 0; i < half_width; ++i) {
301  sum += in(i, j, k);
302  }
303  // Fill left border.
304  for (int i = 0; i < half_width + 1; ++i) {
305  sum += in(i + half_width, j, k);
306  out(i, j, k) = sum;
307  }
308  // Fill interior.
309  for (int i = half_width + 1; i < in.Height() - half_width; ++i) {
310  sum -= in(i - half_width - 1, j, k);
311  sum += in(i + half_width, j, k);
312  out(i, j, k) = sum;
313  }
314  // Fill right border.
315  for (int i = in.Height() - half_width; i < in.Height(); ++i) {
316  sum -= in(i - half_width - 1, j, k);
317  out(i, j, k) = sum;
318  }
319  }
320  }
321 }
322 
323 void BoxFilter(const Array3Df& in, int box_width, Array3Df* out) {
324  Array3Df tmp;
325  BoxFilterHorizontal(in, box_width, &tmp);
326  BoxFilterVertical(tmp, box_width, out);
327 }
328 
329 void LaplaceFilter(unsigned char* src,
330  unsigned char* dst,
331  int width,
332  int height,
333  int strength) {
334  for (int y = 1; y < height - 1; y++)
335  for (int x = 1; x < width - 1; x++) {
336  const unsigned char* s = &src[y * width + x];
337  int l = 128 + s[-width - 1] + s[-width] + s[-width + 1] + s[1] -
338  8 * s[0] + s[1] + s[width - 1] + s[width] + s[width + 1];
339  int d = ((256 - strength) * s[0] + strength * l) / 256;
340  if (d < 0)
341  d = 0;
342  if (d > 255)
343  d = 255;
344  dst[y * width + x] = d;
345  }
346 }
347 
348 } // namespace libmv
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei height
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint y
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei width
ATTR_WARN_UNUSED_RESULT const BMLoop * l
static DBVT_INLINE btScalar size(const btDbvtVolume &a)
Definition: btDbvt.cpp:52
static T sum(const btAlignedObjectArray< T > &items)
3D array (row, column, channel).
Definition: array_nd.h:325
void Resize(int height, int width, int depth=1)
Definition: array_nd.h:334
int Depth() const
Definition: array_nd.h:340
int Height() const
Definition: array_nd.h:338
int Width() const
Definition: array_nd.h:339
int Stride(int axis) const
Return the distance between neighboring elements along axis.
Definition: array_nd.h:172
T * Data()
Pointer to the first element of the array.
Definition: array_nd.h:186
#define static_convolution(size)
SyclQueue void void size_t num_bytes SyclQueue void const char void *memory_device_pointer KernelContext int kernel
SyclQueue void void * src
void BoxFilterVertical(const Array3Df &in, int window_size, Array3Df *out_pointer)
Definition: convolve.cc:289
Eigen::VectorXd Vec
Definition: numeric.h:61
void ImageDerivatives(const Array3Df &in, double sigma, Array3Df *gradient_x, Array3Df *gradient_y)
Definition: convolve.cc:191
void ConvolveHorizontal(const Array3Df &in, const Vec &kernel, Array3Df *out_pointer, int plane)
Definition: convolve.cc:168
void ConvolveVertical(const Array3Df &in, const Vec &kernel, Array3Df *out_pointer, int plane)
Definition: convolve.cc:175
void LaplaceFilter(unsigned char *src, unsigned char *dst, int width, int height, int strength)
Definition: convolve.cc:329
double GaussianDerivative(double x, double sigma)
Definition: convolve.h:41
double NormalizeL1(TVec *x)
Definition: numeric.h:221
void BoxFilterHorizontal(const Array3Df &in, int window_size, Array3Df *out_pointer)
Definition: convolve.cc:255
void BlurredImageAndDerivativesChannels(const Array3Df &in, double sigma, Array3Df *blurred_and_gradxy)
Definition: convolve.cc:233
double Gaussian(double x, double sigma)
Definition: convolve.h:32
void FastConvolve(const Vec &kernel, int width, int height, const float *src, int src_stride, int src_line_stride, float *dst, int dst_stride)
Definition: convolve.cc:67
double GaussianInversePositive(double y, double sigma)
Definition: convolve.h:45
void ConvolveGaussian(const Array3Df &in, double sigma, Array3Df *out_pointer)
Definition: convolve.cc:182
void Convolve(const Array3Df &in, const Vec &kernel, Array3Df *out_pointer, int plane)
Definition: convolve.cc:102
void ComputeGaussianKernel(double sigma, Vec *kernel, Vec *derivative)
Definition: convolve.cc:32
void BoxFilter(const Array3Df &in, int box_width, Array3Df *out)
Definition: convolve.cc:323
void BlurredImageAndDerivatives(const Array3Df &in, double sigma, Array3Df *blurred_image, Array3Df *gradient_x, Array3Df *gradient_y)
Definition: convolve.cc:208
static const pxr::TfToken out("out", pxr::TfToken::Immortal)