gauss2.cuh
1.9 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
#ifndef STIM_CUDA_GAUSS2_H
#define STIM_CUDA_GAUSS2_H
#include <stim/image/image.h>
#include <stim/math/filters/sepconv2.cuh>
#include <stim/math/constants.h>
namespace stim {
template<typename T>
T gauss1d(T x, T mu, T std) {
return ((T)1 / (T)sqrt(stim::TAU * std * std) * exp(-(x - mu)*(x - mu) / (2 * std*std)));
}
///Perform a 2D gaussian convolution on an input image.
///@param in is a pointer to the input image
///@param stdx is the standard deviation (in pixels) along the x axis
///@param stdy is the standard deviation (in pixels) along the y axis
///@param nstds specifies the number of standard deviations of the Gaussian that will be kept in the kernel
template<typename T, typename K>
stim::image<T> cpu_gauss2(const stim::image<T>& in, K stdx, K stdy, size_t nstds = 3) {
size_t kx = stdx * nstds * 2 + 1; //calculate the kernel sizes
size_t ky = stdy * nstds * 2 + 1;
size_t X = in.width() - kx + 1; //calculate the size of the output image
size_t Y = in.height() - ky + 1;
stim::image<T> r(X, Y, in.channels()); //create an output image
K* gx = (K*)malloc(kx * sizeof(K)); //allocate space for the gaussian kernels
K* gy = (K*)malloc(ky * sizeof(K));
K mux = (K)kx / (K)2; //calculate the mean of the gaussian (center of the kernel)
K muy = (K)ky / (K)2;
for (size_t xi = 0; xi < kx; xi++) //evaluate the kernels
gx[xi] = gauss1d((K)xi, mux, stdx);
for (size_t yi = 0; yi < ky; yi++)
gy[yi] = gauss1d((K)yi, muy, stdy);
std::vector< stim::image<T> > clist = in.split(); //split the input image into channels
std::vector< stim::image<T> > R = r.split(); //split the output image into channels
for (size_t c = 0; c < in.channels(); c++) //for each channel
cpu_sepconv2(R[c].data(), clist[c].data(), gx, gy, clist[c].width(), clist[c].height(), kx, ky);
r.merge(R); //merge the blurred channels into the final image
return r;
}
}
#endif