diff --git a/lib/libstereo/src/algorithms/clustersf.cu b/lib/libstereo/src/algorithms/clustersf.cu
index db488112343fb0c76b0d15af0487b266b1340586..cacbe1e40cbb7da52d90590cc249aef191971adb 100644
--- a/lib/libstereo/src/algorithms/clustersf.cu
+++ b/lib/libstereo/src/algorithms/clustersf.cu
@@ -3,6 +3,7 @@
 #include "../filters/salient_gradient.hpp"
 
 #include <opencv2/highgui.hpp>
+#include <opencv2/imgproc.hpp>
 
 struct StereoCSF::Impl {
 	Array2D<uchar> l;
@@ -10,10 +11,12 @@ struct StereoCSF::Impl {
 	Array2D<uchar> gl;
 	Array2D<uchar> gr;
 	Array2D<uchar> temp;
+	Bucket2D<ushort, 64> buckets_l;
 
 	Impl(int width, int height) :
 		l(width, height), r(width, height),
-		gl(width, height), gr(width, height), temp(width, height) {}
+		gl(width, height), gr(width, height), temp(width, height),
+		buckets_l(16, height) {}
 };
 
 StereoCSF::StereoCSF() : impl_(nullptr) {
@@ -32,13 +35,21 @@ void StereoCSF::compute(cv::InputArray l, cv::InputArray r, cv::OutputArray disp
 	mat2gray(l, impl_->l);
 	mat2gray(r, impl_->r);
 
-	SalientGradient sg = {impl_->l.data(), impl_->gl.data(), impl_->temp.data(), impl_->l.width, impl_->l.height};
-	parallel1DWarpSM(sg, l.rows(), l.cols());
+	SalientGradient sgl = {impl_->l.data(), impl_->gl.data(), impl_->temp.data(), impl_->buckets_l.data(), impl_->l.width, impl_->l.height};
+	parallel1DWarpSM(sgl, l.rows(), l.cols());
+	//SalientGradient sgr = {impl_->r.data(), impl_->gr.data(), impl_->temp.data(), impl_->r.width, impl_->r.height};
+	//parallel1DWarpSM(sgr, r.rows(), r.cols());
 
 	cv::Mat tmp;
-	impl_->gl.toGpuMat().download(tmp);
-	cv::resize(tmp,tmp, cv::Size(tmp.cols/2, tmp.rows/2));
-	cv::imshow("Gradients", tmp);
+	impl_->buckets_l.toGpuMat().download(tmp);
+	tmp.convertTo(tmp, CV_8UC1, 4.0);
+	cv::resize(tmp,tmp, cv::Size(tmp.cols*40, tmp.rows/2));
+	cv::applyColorMap(tmp, tmp, cv::COLORMAP_TURBO);
+	cv::imshow("Gradients Left", tmp);
+
+	//impl_->gr.toGpuMat().download(tmp);
+	//cv::resize(tmp,tmp, cv::Size(tmp.cols/2, tmp.rows/2));
+	//cv::imshow("Gradients Right", tmp);
 
 	cudaSafeCall(cudaDeviceSynchronize());
 	disparity.create(l.rows(), l.cols(), CV_32F);
diff --git a/lib/libstereo/src/bucket2d.hpp b/lib/libstereo/src/bucket2d.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..6eb71d42e3cd74a62fd38f05bd836c23a1f3cd56
--- /dev/null
+++ b/lib/libstereo/src/bucket2d.hpp
@@ -0,0 +1,136 @@
+#ifndef _FTL_LIBSTEREO_BUCKET2D_HPP_
+#define _FTL_LIBSTEREO_BUCKET2D_HPP_
+
+#include "memory.hpp"
+
+template<typename T, int DEPTH>
+class Bucket2D {
+public:
+	Bucket2D() : width(0), height(0), needs_free_(false) {
+		data_.buckets = nullptr;
+		data_.counters = nullptr;
+	}
+
+	Bucket2D(int w, int h) : width(w), height(h), needs_free_(true) {
+		data_.counters = allocateMemory2D<uint>(w, h, data_.counter_pitch);
+		data_.buckets = allocateMemory2D<T>(w*DEPTH, h, data_.bucket_pitch);
+	}
+
+	~Bucket2D() {
+		free();
+	}
+
+	void free() {
+		if (needs_free_ && data_.buckets) {
+			freeMemory(data_.counters);
+			freeMemory(data_.buckets);
+		}
+	}
+
+	Bucket2D<T, DEPTH> &operator=(const Bucket2D<T, DEPTH> &c) {
+		data_ = c.data_;
+		width = c.width;
+		height = c.height;
+		needs_free_ = false;
+		return *this;
+	}
+
+	struct Data {
+		__host__ __device__ inline uint& operator() (const int y, const int x) {
+			return counters[x + y*counter_pitch];
+		}
+
+		__host__ __device__ inline const uint& operator() (const int y, const int x) const {
+			return counters[x + y*counter_pitch];
+		}
+
+		__host__ __device__ inline const T& operator() (const int y, const int x, int d) const {
+			return buckets[d + x*DEPTH + y*bucket_pitch];
+		}
+
+		__host__ __device__ inline void add (const int y, const int x, T v) {
+			uint ix = atomicInc(&counters[x + y*counter_pitch], DEPTH);
+			buckets[ix + x*DEPTH + y*bucket_pitch] = v;
+		}
+
+		__host__ __device__ inline uint *ptr(const int y) { return &counters[y*counter_pitch]; }
+		__host__ __device__ inline const uint *ptr(const int y) const { return &counters[y*counter_pitch]; }
+		__host__ __device__ inline const T *ptr(const int y, const int x) const { return &buckets[x*DEPTH+y*bucket_pitch]; }
+
+		uint bucket_pitch;
+		uint counter_pitch;
+		T *buckets;
+		uint *counters;
+	};
+
+	void create(int w, int h) {
+		if (w == width && h == height) return;
+		width = w;
+		height = h;
+		free();
+		needs_free_ = true;
+		data_.counters = allocateMemory2D<uint>(w, h, data_.counter_pitch);
+		data_.buckets = allocateMemory2D<T>(w*DEPTH, h, data_.bucket_pitch);
+	}
+
+	inline Data &data() { return data_; }
+	inline const Data &data() const { return data_; }
+
+	void toMat(cv::Mat &m) {
+		#ifdef USE_GPU
+		cv::cuda::GpuMat gm;
+		toGpuMat(gm);
+		gm.download(m);
+		#else
+		m = cv::Mat(height, width, cv::traits::Type<int>::value, data_.counters);
+		#endif
+	}
+
+	cv::Mat toMat() {
+		#ifdef USE_GPU
+		cv::cuda::GpuMat gm;
+		toGpuMat(gm);
+		cv::Mat m;
+		gm.download(m);
+		return m;
+		#else
+		return cv::Mat(height, width, cv::traits::Type<int>::value, data_.counters);
+		#endif
+	}
+
+	const cv::Mat toMat() const {
+		#ifdef USE_GPU
+		cv::cuda::GpuMat gm(height, width, cv::traits::Type<int>::value, (void*)data_.counters, size_t(size_t(data_.counter_pitch)*sizeof(uint)));
+		cv::Mat m;
+		gm.download(m);
+		return m;
+		#else
+		return cv::Mat(height, width, cv::traits::Type<int>::value, data_.counters);
+		#endif
+	}
+
+	void toGpuMat(cv::cuda::GpuMat &m) {
+		#ifdef USE_GPU
+		m = cv::cuda::GpuMat(height, width, cv::traits::Type<int>::value, (void*)data_.counters, size_t(size_t(data_.counter_pitch)*sizeof(uint)));
+		#else
+		// TODO
+		#endif
+	}
+
+	cv::cuda::GpuMat toGpuMat() {
+		#ifdef USE_GPU
+		return cv::cuda::GpuMat(height, width, cv::traits::Type<int>::value, (void*)data_.counters, size_t(size_t(data_.counter_pitch)*sizeof(uint)));
+		#else
+		return cv::cuda::GpuMat(height, width, cv::traits::Type<int>::value);
+		#endif
+	}
+
+	int width;
+	int height;
+
+private:
+	Data data_;
+	bool needs_free_;
+};
+
+#endif
diff --git a/lib/libstereo/src/filters/salient_gradient.hpp b/lib/libstereo/src/filters/salient_gradient.hpp
index 093319518e75d51d5b220eb48709815be167199d..872e8c9a658b4511348952e4028a857180aff3e6 100644
--- a/lib/libstereo/src/filters/salient_gradient.hpp
+++ b/lib/libstereo/src/filters/salient_gradient.hpp
@@ -3,11 +3,13 @@
 
 #include "../util.hpp"
 #include "../array2d.hpp"
+#include "../bucket2d.hpp"
 
 struct SalientGradient {
 	Array2D<uchar>::Data image;
-	Array2D<uchar>::Data output;
+	Array2D<uchar>::Data angle;
 	Array2D<uchar>::Data magnitude;
+	Bucket2D<ushort, 64>::Data buckets;
 
 	int width, height;
 
@@ -48,31 +50,38 @@ struct SalientGradient {
 			// Reset histogram
 			//for (int i=thread.x; i < 32; i+=32) wsm.gradient_histogram[i] = 0;
 			wsm.gradient_histogram[thread.x] = 0;
+			//int maxmag = 0;
 
 			for (int x=thread.x; x<size.x; x+=stride.x) {
 				auto g = calculateGradient(x,y);
-				output(y,x) = uchar((g.y+PI2) / PI * 63.0f);
+				angle(y,x) = uchar(min(15,int((g.y+PI2) / PI * 16.0f)));
 				magnitude(y,x) = uchar(g.x);
 
-				//maxmag = warpMax(max(maxmag,int(g.x)));
+				//maxmag = max(maxmag,int(g.x));
 				atomicAdd(&wsm.gradient_histogram[min(31,int(g.x / 361.0f * 32.0f))], 1);
 				//atomicAdd(&wsm.gradient_histogram[0], 1);
 			}
 
-			__syncwarp();
+			//maxmag = int(float(warpMax(maxmag)) * 0.8f);
 
 			uchar gthresh = min(255, scan(wsm.gradient_histogram, thread.x, float(width)*0.95f) * (256/32));
 
 			// Apply threshold
 			for (int x=thread.x; x<size.x; x+=stride.x) {
 				uchar thresh = gthresh;
-				for (int u=-5; u<=5; ++u) {
-					thresh = (x+u >= 0 && x+u < width) ? max(thresh, magnitude(y,x+u)-abs(u)) : thresh;
+				for (int u=-3; u<=3; ++u) {
+					thresh = (x+u >= 0 && x+u < width) ? max(thresh, magnitude(y,x+u)) : thresh;
 				}
 
 				uchar m = magnitude(y,x);
-				//if (m < max(thresh, uchar(float(maxmag)*0.2f))) output(y,x) = 0;
-				output(y,x) = (m < thresh) ? 0 : 255;
+				//if (m < thresh) angle(y,x) = 0;
+				//output(y,x) = (m < thresh) ? 0 : 255;
+				if (m >= thresh) {
+					int a = angle(y,x);
+					buckets.add(y, a, ushort(x));
+					buckets.add(y, (a > 0) ? a-1 : 15, ushort(x));
+					buckets.add(y, (a < 15) ? a+1 : 0, ushort(x));
+				}
 			}
 
 			// Next step would be to bucket the results