22 #ifndef CVD_CONVOLUTION_H_
23 #define CVD_CONVOLUTION_H_
50 unsigned int i, argmax=0;
51 std::vector<double> kernel(k.size());
52 for (i=0;i<k.size();i++) {
53 double x = i +0.5 - k.size()/2.0;
54 sum += kernel[i] =
exp(-x*x/(2*stddev*stddev));
55 if (kernel[i] > kernel[argmax])
59 for (i=0;i<k.size();i++)
60 finalSum += k[i] = (T)(kernel[i]*maxval/kernel[argmax]);
71 template <
class S,
class T>
72 T
scaleKernel(
const std::vector<S>& k, std::vector<T>& scaled, T maxval)
74 unsigned int i,argmax=0;
75 for (i=1;i<k.size();i++)
78 scaled.resize(k.size());
80 for (i=0;i<k.size();i++)
81 sum += (scaled[i] = (T)((k[i]*maxval)/k[argmax]));
93 T* end = src + w*(h-4);
95 T sum= (T)(0.0544887*(src[0]+src[4*w])
96 + 0.2442010*(src[w]+src[3*w])
97 + 0.4026200*src[2*w]);
102 for (i=h-5;i>=0;i--) {
103 T* src = I.
data()+i*w;
106 T sum= (T)(0.0544887*(src[0]+src[4])
107 + 0.2442010*(src[1]+src[3])
115 namespace Exceptions {
119 namespace Convolution {
129 what =
"Incompatible image sizes in " +
function;
139 what =
"Odd sized kernel required in " +
function;
159 const double factor = 1.0/(win.
x*win.
y);
160 std::vector<sum_type> buffer(w*win.
y);
161 std::vector<sum_type> sums_v(w);
162 sum_type* sums = &sums_v[0];
163 sum_type* next_row = &buffer[0];
164 sum_type* oldest_row = &buffer[0];
166 const T* input = I.
data();
167 T* output = J[hwin.
y] - hwin.
x;
168 for (
int i=0; i<h; i++) {
169 sum_type hsum=sum_type();
170 const T* back = input;
172 for (j=0; j<win.
x-1; j++)
182 differences(sums+win.
x-1, oldest_row+win.
x-1, sums+win.
x-1, w-win.
x+1);
185 if (oldest_row == &buffer[0] + w*win.
y)
186 oldest_row = &buffer[0];
190 if (next_row == &buffer[0] + w*win.
y)
191 next_row = &buffer[0];
212 static const wider S = (A+B+C+B+A);
213 int width = I.
size().
x;
214 int height = I.
size().
y;
217 for (i=0; i<height; i++) {
222 p[0] = (T)(((c+c)*A+(b+b)*B + a*C) /S);
223 p[1] = (T)(((b+d)*A+(a+c)*B + b*C) /S);
224 for (j=0;j<width-4;j++,p++) {
226 p[2] = (T)(((a+e)*A + (b+d)*B + c*C)/S);
227 a = b; b = c; c = d; d = e;
229 p[2] = (T)(((a+c)*A + (b+d)*B + c*C) /S);
230 p[3] = (T)(((b+b)*A + (c+c)*B + d*C) /S);
233 for (j=0;j<width;j++) {
237 p[0] = (T)(((p[2*width]+p[2*width])*A+(b+b)*B + a*C) /S);
238 p[width] = (T)(((b+p[width*3])*A+(a+p[2*width])*B + b*C) /S);
239 for (i=0;i<height-4;i++) {
240 wider c = p[2*width];
241 p[2*width] = (T)(((a+p[4*width])*A + (b+p[3*width])*B + c*C)/S);
245 wider c = p[2*width];
246 p[2*width] = (T)(((a+c)*A + (b+p[width*3])*B + c*C) /S);
247 p[3*width] = (T)(((b+b)*A + (c+c)*B + p[width*3]*C) /S);
254 static const wider S = (A+B+C+D+C+B+A);
255 int width = I.
size().
x;
256 int height = I.
size().
y;
259 for (i=0; i<height; i++) {
264 p[0] = (T)(((d+d)*A + (c+c)*B + (b+b)*C + a*D)/S);
265 p[1] = (T)(((c+p[4])*A + (b+d)*B + (a+c)*C + b*D)/S);
266 p[2] = (T)(((b+p[5])*A + (a+p[4])*B + (b+d)*C + c*D)/S);
267 for (j=0;j<width-6;j++,p++) {
269 p[3] = (T)(((a+p[6])*A + (b+p[5])*B + (c+p[4])*C + d*D)/S);
274 p[3] = (T)(((a+e)*A + (b+p[5])*B + (c+e)*C + d*D)/S);
275 p[4] = (T)(((b+d)*A + (c+e)*B + (d+p[5])*C + e*D)/S);
276 p[5] = (T)(((c+c)*A + (d+d)*B + (e+e)*C + p[5]*D)/S);
279 for (j=0;j<width;j++) {
283 wider c = p[2*width];
284 wider d = p[3*width];
285 p[0] = (T)(((d+d)*A + (c+c)*B + (b+b)*C + a*D)/S);
286 p[width] = (T)(((c+p[4*width])*A + (b+d)*B + (a+c)*C + b*D)/S);
287 p[2*width] = (T)(((b+p[5*width])*A + (a+p[4*width])*B + (b+d)*C + c*D)/S);
288 for (i=0;i<height-6;i++) {
290 p[3*width] = (T)(((a+p[width*6])*A + (b+p[width*5])*B + (c+p[width*4])*C+d*D)/S);
295 wider e = p[4*width];
296 p[3*width] = (T)(((a+e)*A + (b+p[5*width])*B + (c+e)*C + d*D)/S);
297 p[4*width] = (T)(((b+d)*A + (c+e)*B + (d+p[5*width])*C + e*D)/S);
298 p[5*width] = (T)(((c+c)*A + (d+d)*B + (e+e)*C + p[5*width]*D)/S);
307 int r = (int)kernel.size()/2;
310 double factor = 1.0/divisor;
313 for (i=0; i<h-2*r; i++,src+=w) {
314 sum_type sum = src[r*w]*kernel[r], v;
316 sum += (src[m*w] + src[(2*r-m)*w]) * kernel[m];
317 *(src) = static_cast<T>(sum * factor);
320 int offset = r*w + r;
321 for (i=h-2*r-1;i>=0;i--) {
323 for (j=0;j<w-2*r;j++, src++) {
324 sum_type sum = src[r] * kernel[r], v;
326 sum += (src[m] + src[2*r-m])*kernel[m];
327 *(src+offset) = static_cast<T>(sum*factor);
334 static inline const B*
get(
const A* row,
int w, B* rowbuf) {
341 static inline const T*
get(
const T* row, int , T* ) {
351 static inline void cast_copy(
const T* from, S* to,
int count) {
352 for (
int i=0; i<count; i++)
353 to[i] = static_cast<S>(from[i]);
358 static inline void cast_copy(
const T* from, T* to,
int count) {
366 template <
class S>
static inline T
at(
const T* input,
const S& factor,
const S* kernel) {
return ConvolveMiddle<T,-1,C>::at(input,factor, kernel, N); }
370 template <
class S>
static inline T
at(
const T* input,
const S& factor,
const S* kernel) {
return ConvolveMiddle<T,N-1>::at(input,factor, kernel) + (input[-N]+input[N])*kernel[N-1]; }
374 template <
class S>
static inline T
at(
const T* input,
const S& factor,
const S* kernel,
int ksize) {
375 T hsum = *input * factor;
376 for (
int k=0; k<ksize; k++)
377 hsum += (input[-k-1] + input[k+1]) * kernel[k];
383 template <
class S>
static inline T
at(
const T* input,
const S& factor,
const S* kernel,
int ksize) {
384 T hsum = *input * factor;
385 for (
int k=0; k<ksize; k++)
386 hsum += (input[-k-1] + input[k+1]) * kernel[k];
392 template <
class S>
static inline T
at(
const T* input,
const S& factor,
const S* ) {
return *input * factor; }
395 template <
class T,
class S>
inline const T*
convolveMiddle(
const T* input,
const S& factor,
const S* kernel,
int ksize,
int n, T* output) {
396 #define CALL_CM(I) for (int j=0; j<n; ++j, ++input, ++output) { *output = ConvolveMiddle<T,I>::at(input, factor, kernel); } break
425 int ksize = (int)ceil(sigmas*sigma);
427 std::vector<sum_comp_type> kernel(ksize);
428 sum_comp_type ksum = sum_comp_type();
429 for (
int i=1; i<=ksize; i++)
430 ksum += (kernel[i-1] = static_cast<sum_comp_type>(
exp(-i*i/(2*sigma*sigma))));
431 for (
int i=0; i<ksize; i++)
432 kernel[i] /= (2*ksum+1);
433 double factor = 1.0/(2*ksum+1);
442 sum_type* rowbuf = aligned_rowbuf.
data();
443 sum_type* outbuf = aligned_outbuf.
data();
445 std::vector<sum_type*> rows(swin+1);
446 for (
int k=0;k<swin+1;k++)
447 rows[k] = buffer.
data() + k*w;
449 T* output = out.
data();
450 for (
int i=0; i<h; i++) {
451 sum_type* next_row = rows[swin];
454 for (
int j=0; j<ksize; j++) {
455 sum_type hsum =
static_cast<sum_type
>(input[j] * factor);
456 for (
int k=0; k<ksize; k++)
457 hsum += (input[std::max(j-k-1,0)] + input[j+k+1]) * kernel[k];
462 input = convolveMiddle<sum_type, sum_comp_type>(input,
static_cast<sum_comp_type
>(factor), &kernel.front(), ksize, w-swin, next_row+ksize);
464 for (
int j=w-ksize; j<w; j++, input++) {
465 sum_type hsum =
static_cast<sum_type
>(*input * factor);
466 const int room = w-j;
467 for (
int k=0; k<ksize; k++) {
468 hsum += (input[-k-1] + input[std::min(k+1,room-1)]) * kernel[k];
474 const sum_type* middle_row = rows[ksize];
476 for (
int k=0; k<ksize; k++) {
477 const sum_comp_type m = kernel[k];
478 const sum_type* row1 = rows[ksize-k-1];
479 const sum_type* row2 = rows[ksize+k+1];
485 for (
int r=0; r<ksize; r++) {
486 const sum_type* middle_row = rows[ksize+r+1];
488 for (
int k=0; k<ksize; k++) {
489 const sum_comp_type m = kernel[k];
490 const sum_type* row1 = rows[ksize+r-k];
491 const sum_type* row2 = rows[std::min(ksize+r+k+2, swin)];
498 }
else if (i == swin-1) {
499 for (
int r=0; r<ksize; r++) {
500 const sum_type* middle_row = rows[r+1];
502 for (
int k=0; k<ksize; k++) {
503 const sum_comp_type m = kernel[k];
504 const sum_type* row1 = rows[std::max(r-k-1,0)+1];
505 const sum_type* row2 = rows[r+k+2];
513 sum_type* tmp = rows[0];
514 for (
int r=0;r<swin; r++)
522 void van_vliet_blur(
const double b[],
const SubImage<float> in, SubImage<float> out);
524 void convolveGaussian(
const BasicImage<float>& I, BasicImage<float>& out,
double sigma,
double sigmas=3.0);
525 void convolveGaussian_fir(
const BasicImage<float>& I, BasicImage<float>& out,
double sigma,
double sigmas=3.0);
531 const int w = I.
size().
x;
532 O* o = out.
data()+w+1;
534 const double cross = k1*k2;
538 const double sum = k1*(a[0] + a[2] + a[w*2] + a[w*2+2]) + cross*(a[1] + a[w*2+1] + a[w] + a[w+2]) + k2*a[w+1];
539 *o++ = Pixel::scalar_convert<O,T,double>(sum);
static T at(const T *input, const S &factor, const S *kernel)
void cast_copy(const T *from, S *to, int count)
void convolveSeparableSymmetric(Image< T > &I, const std::vector< K > &kernel, K divisor)
int totalsize() const
What is the total number of elements in the image (i.e. size().x * size().y), including padding...
void convolveGaussian_fir(const BasicImage< float > &I, BasicImage< float > &out, double sigma, double sigmas=3.0)
T gaussianKernel(std::vector< T > &k, T maxval, double stddev)
void convolveGaussian(BasicImage< T > &I, double sigma, double sigmas=3.0)
const T * convolveMiddle(const T *input, const S &factor, const S *kernel, int ksize, int n, T *output)
static T at(const T *input, const S &factor, const S *kernel, int ksize)
IncompatibleImageSizes(const std::string &function)
static void cast_copy(const T *from, S *to, int count)
void copy(const BasicImage< S > &in, BasicImage< T > &out, ImageRef size=ImageRef(-1,-1), ImageRef begin=ImageRef(), ImageRef dst=ImageRef())
OddSizedKernelRequired(const std::string &function)
void convolveGaussian5_1(Image< byte > &I)
void zeroPixels(T *pixels, int count)
static void cast_copy(const T *from, T *to, int count)
Matrix< R, C, P > exp(const Matrix< R, C, P, B > &m)
ImageRef size() const
What is the size of this image?
static T at(const T *input, const S &factor, const S *kernel, int ksize)
const T * data() const
Returns the raw image data.
T scaleKernel(const std::vector< S > &k, std::vector< T > &scaled, T maxval)
void convolve_gaussian_3(const BasicImage< T > &I, BasicImage< O > &out, K k1, K k2)
void add_multiple_of_sum(const A *a, const A *b, const C &c, B *out, size_t count)
void convolveSymmetric(Image< T > &I)
void assign_multiple(const A *a, const B &c, C *out, size_t count)
void compute_triggs_M(const double b[], double M[][3])
static T at(const T *input, const S &factor, const S *kernel)
void convolveWithBox(const BasicImage< T > &I, BasicImage< T > &J, ImageRef hwin)
std::string what
The error message.
const B * getPixelRowTyped(const A *row, int n, B *rowbuf)
static T at(const T *input, const S &factor, const S *)
void differences(const A *a, const A *b, B *diff, size_t count)
void compute_van_vliet_b(double sigma, double b[])
static const B * get(const A *row, int w, B *rowbuf)
void van_vliet_blur(const double b[], const CVD::SubImage< float > in, CVD::SubImage< float > out)