Otsu fixed by implementing MATLAB's method

src/c/Common/CHelpers.cpp

 #include "CHelpers.h"
+#include <vector>
+#include "Defines.h"
 
 float* createEllipsoidKernel(Vec<size_t> radii, Vec<size_t>& kernelDims)
 {
@@ -31,28 +33,56 @@ float* createEllipsoidKernel(Vec<size_t> radii, Vec<size_t>& kernelDims)
 
 int calcOtsuThreshold(const double* normHistogram, int numBins)
 {
-	//code modified from http://www.dandiggins.co.uk/arlib-9.html
-	double totalMean = 0.0f;
+	double* omegas = new double[numBins];
+	double* mus = new double[numBins];
+	double* sigma_b_squared = new double[numBins];
+	double maxVal = 0.0;
+
+	std::vector<int> maxs;
+	maxs.reserve(numBins);
+
+	omegas[0] = normHistogram[0];
+	mus[0] = 0.0;
+
+	double lastOmega = normHistogram[0];
+	double lastMu = 0.0;
+	for (int i=1; i<numBins; ++i)
+	{
+		omegas[i] = lastOmega + normHistogram[i];
+		mus[i] = lastMu + (i*normHistogram[i]);
+
+		lastOmega = omegas[i];
+		lastMu = mus[i];
+	}
+
 	for (int i=0; i<numBins; ++i)
-		totalMean += i*normHistogram[i];
+	{
+		sigma_b_squared[i] = SQR(mus[numBins-1] * omegas[i] - mus[i]) / (omegas[i] * (1.0 - omegas[i]));
+		if (maxVal < sigma_b_squared[i])
+			maxVal = sigma_b_squared[i];
+	}
 
-	double class1Prob=0, class1Mean=0, temp1, curThresh;
-	double bestThresh = 0.0;
-	int bestIndex = 0;
 	for (int i=0; i<numBins; ++i)
 	{
-		class1Prob += normHistogram[i];
-		class1Mean += i * normHistogram[i];
+		if (sigma_b_squared[i] > maxVal*0.95)
+			maxs.push_back(i);
+	}
 
-		temp1 = totalMean * class1Prob - class1Mean;
-		curThresh = (temp1*temp1) / (class1Prob*(1.0f-class1Prob));
+	double mean = 0.0;
 
-		if(curThresh>bestThresh)
+	if (maxVal!=0 && !maxs.empty())
+	{
+		for (std::vector<int>::iterator it=maxs.begin(); it!=maxs.end(); ++it)
 		{
-			bestThresh = curThresh;
-			bestIndex = i;
+			mean += *it;
 		}
+
+		mean /= maxs.size();
 	}
 
-	return bestIndex;
+	delete[] omegas;
+	delete[] mus;
+	delete[] sigma_b_squared;
+
+	return floor(mean);
 }
\ No newline at end of file

src/c/Common/CudaHistogram.cuh

@@ -152,7 +152,7 @@ PixelType cOtsuThresholdValue(const PixelType* imageIn, Vec<size_t> dims, int de
 
 	delete[] hist;
 
-	double binSize = (maxVal-minVal)/arraySize;
+	double binSize = ((double)maxVal-minVal)/arraySize;
 
 	PixelType thrsh = (PixelType)(minVal + binSize*theshIdx);
 

src/c/Common/Defines.h

 #pragma once
 
-#define NUM_BINS (255)
+#define NUM_BINS (256)
 #define MAX_KERNEL_DIM (25)
 
 #define SQR(x) ((x)*(x))