added 2D gaussian filter (CPU implementation), and added the ability to split an…

…d merge stim::image by channel

added 2D gaussian filter (CPU implementation), and added the ability to split an…
…d merge stim::image by channel
David Mayerich
1 parent dbeb83f2
Showing 2 changed files with 67 additions and 6 deletions Show diff stats
stim/image/image.h
stim/math/filters/gauss2.h
@@ -54,7 +54,7 @@ class image{
 	size_t bytes(){ return size() * sizeof(T); }
-	inline size_t idx(size_t x, size_t y, size_t c = 0){
+	inline size_t idx(size_t x, size_t y, size_t c = 0) const {
 		return y * R[0] * R[1] + x * R[0] + c;
 	}
@@ -125,6 +125,11 @@ public:
 		free(img);
 	}
+	///Resize an image
+	void resize(size_t x, size_t y, size_t c = 1) {
+		allocate(x, y, c);
+	}
+
 	stim::image<T>& operator=(const stim::image<T>& I){
 		init();
 		if(&I == this)									//handle self-assignment
@@ -340,7 +345,7 @@ public:
 	/// Return an image representing a specified channel
 	/// @param c is the channel to be returned
-	image<T> channel(size_t c){		
+	image<T> channel(size_t c) const {		
 		image<T> r(X(), Y(), 1);				//create a new image
 		for(size_t x = 0; x < X(); x++){
 			for(size_t y = 0; y < Y(); y++){
@@ -351,7 +356,7 @@ public:
 	}
 	/// Returns an std::vector containing each channel as a separate image
-	std::vector<image<T>> split() {
+	std::vector<image<T>> split() const {
 		std::vector<image<T>> r;			//create an image array
 		r.resize(C());						//create images for each channel
@@ -361,6 +366,15 @@ public:
 		return r;
 	}
+	/// Merge a series of single-channel images into a multi-channel image
+	void merge(std::vector<image<T>>& list) {
+		size_t x = list[0].width();				//calculate the size of the image
+		size_t y = list[0].height();
+		allocate(x, y, list.size());			//re-allocate the image
+		for (size_t c = 0; c < list.size(); c++)		//for each channel
+			set_channel(list[c].channel(0).data(), c);	//insert the channel into the output image
+	}
+
 	T& operator()(size_t x, size_t y, size_t c = 0){
 		return img[idx(x, y, c)];
 	}
@@ -394,15 +408,15 @@ public:
 		}
 	}
-	size_t channels(){
+	size_t channels() const{
 		return C();
 	}
-	size_t width(){
+	size_t width() const{
 		return X();
 	}
-	size_t height(){
+	size_t height() const{
 		return Y();
 	}
+#ifndef STIM_CUDA_GAUSS2_H
+#define STIM_CUDA_GAUSS2_H
+
+#include <stim/image/image.h>
+#include <stim/math/filters/sepconv2.h>
+#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 = T>
+	stim::image<T> cpu_gauss2(const stim::image<T>& in, K stdx, K stdy, size_t nstds = 3) {
+		size_t kx = stdx * nstds * 2;					//calculate the kernel sizes
+		size_t ky = stdy * nstds * 2;
+		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>> IN = 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(), IN[c].data(), gx, gy, IN[c].width(), IN[c].height(), kx, ky);
+
+		r.merge(R);										//merge the blurred channels into the final image
+		return r;
+	}
+}
+
+
+#endif
 \ No newline at end of file