Implemented calculating Otsu value

src/c/Common/CHelpers.cpp

@@ -31,3 +31,30 @@ double* createEllipsoidKernel(Vec<size_t> radii, Vec<size_t>& kernelDims)
 	return kernel;
 }
 
+int calcOtsuThreshold(const double* normHistogram, int numBins)
+{
+	//code modified from http://www.dandiggins.co.uk/arlib-9.html
+	double totalMean = 0.0f;
+	for (int i=0; i<numBins; ++i)
+		totalMean += i*normHistogram[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];
+
+		temp1 = totalMean * class1Prob - class1Mean;
+		curThresh = (temp1*temp1) / (class1Prob*(1.0f-class1Prob));
+
+		if(curThresh>bestThresh)
+		{
+			bestThresh = curThresh;
+			bestIndex = i;
+		}
+	}
+
+	return bestIndex;
+}
\ No newline at end of file

src/c/Common/CHelpers.h

@@ -4,4 +4,6 @@
 #include <memory.h>
 #include <complex>
 
-double* createEllipsoidKernel(Vec<size_t> radii, Vec<size_t>& kernelDims);
\ No newline at end of file
+double* createEllipsoidKernel(Vec<size_t> radii, Vec<size_t>& kernelDims);
+
+int calcOtsuThreshold(const double* normHistogram, int numBins);

src/c/Common/CudaProcessBuffer.cu

@@ -614,6 +614,19 @@ void CudaProcessBuffer::otsuThresholdFilter(float alpha/*=1.0f*/)
 	throw std::logic_error("The method or operation is not implemented.");
 }
 
+
+double CudaProcessBuffer::otsuThresholdValue(const DevicePixelType* imageIn, Vec<size_t> dims)
+{
+	int arraySize;
+	double* hist = normalizeHistogram(imageIn,dims,arraySize);
+
+	double thrsh = calcOtsuThreshold(hist,arraySize);
+
+	delete[] hist;
+
+	return thrsh;
+}
+
 void CudaProcessBuffer::imagePow(int p)
 {
 	throw std::logic_error("The method or operation is not implemented.");

src/c/Common/CudaProcessBuffer.cuh

@@ -153,6 +153,8 @@ public:
 
 	void otsuThresholdFilter(float alpha=1.0f);
 
+	double otsuThresholdValue(const DevicePixelType* imageIn, Vec<size_t> dims);
+
 	/*
 	*	Raise each pixel to a power
 	*/

src/c/Common/CudaUtilities.cu

@@ -187,31 +187,3 @@ Vec<size_t> createGaussianKernel(Vec<float> sigma, float* kernel, Vec<int>& iter
 
 	return kernelDims;
 }
-
-int calcOtsuThreshold(const double* normHistogram, int numBins)
-{
-	//code modified from http://www.dandiggins.co.uk/arlib-9.html
-	double totalMean = 0.0f;
-	for (int i=0; i<numBins; ++i)
-		totalMean += i*normHistogram[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];
-
-		temp1 = totalMean * class1Prob - class1Mean;
-		curThresh = (temp1*temp1) / (class1Prob*(1.0f-class1Prob));
-
-		if(curThresh>bestThresh)
-		{
-			bestThresh = curThresh;
-			bestIndex = i;
-		}
-	}
-
-	return bestIndex;
-}

src/c/Common/CudaUtilities.cuh

@@ -112,5 +112,3 @@ struct Lock
 Vec<size_t> createGaussianKernel(Vec<float> sigma, float* kernel, int& iterations);
 
 Vec<size_t> createGaussianKernel(Vec<float> sigma, float* kernel, Vec<int>& iterations);
-
-int calcOtsuThreshold(const double* normHistogram, int numBins);

src/c/CudaMex/MexCommand.cpp

@@ -33,7 +33,7 @@ void MexCommand::init()
 // 	REGISTER_COMMAND(NormalizedCovariance);
  	REGISTER_COMMAND(NormalizedHistogram);
 // 	REGISTER_COMMAND(OtsuThresholdFilter);
-// 	REGISTER_COMMAND(OtsuThesholdValue);
+ 	REGISTER_COMMAND(OtsuThesholdValue);
  	REGISTER_COMMAND(ReduceImage);
 // 	REGISTER_COMMAND(SumArray);
 // 	REGISTER_COMMAND(ThresholdFilter);

src/c/CudaMex/MexCommand.h

@@ -319,19 +319,19 @@ public:
 // 	virtual std::string printUsage();
 //virtual std::string printHelp();
 // };
-// 
-// class OtsuThesholdValue : MexCommand
-// {
-// public:
-// 	OtsuThesholdValue(){}
-// 	virtual ~OtsuThesholdValue(){}
-// 
-// 	virtual void execute(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]);
-// 	virtual std::string check(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]);
-// 	virtual std::string printUsage();
-//virtual std::string printHelp();
-// };
-// 
+
+class OtsuThesholdValue : MexCommand
+{
+public:
+	OtsuThesholdValue(){}
+	virtual ~OtsuThesholdValue(){}
+
+	virtual void execute(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]);
+	virtual std::string check(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]);
+	virtual std::string printUsage();
+	virtual std::string printHelp();
+};
+
 // class ImagePow : MexCommand
 // {
 // public:

src/c/CudaMex/OtsuThesholdValue.cpp

-// #include "MexCommand.h"
-// 
-// void OtsuThesholdValue::execute( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[] )
-// {
-// 	Vec<size_t> imageDims;
-// 	ImageContainer* imageIn;
-// 	setupImagePointers(prhs[0],&imageIn);
-// 
-// 	double thresh;
-// 
-// 	thresh = (double)otsuThesholdValue(imageIn);
-// 
-// 	plhs[0] = mxCreateDoubleScalar(thresh);
-// 
-// 	delete imageIn;
-// }
-// 
-// std::string OtsuThesholdValue::check( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[] )
-// {
-// 	if (nrhs!=1)
-// 		return "Incorrect number of inputs!";
-// 
-// 	if (nlhs!=1)
-// 		return "Requires one output!";
-// 
-// 	if (!mxIsUint8(prhs[0]))
-// 		return "Image has to be formated as a uint8!";
-// 
-// 	size_t numDims = mxGetNumberOfDimensions(prhs[0]);
-// 	if (numDims>3 || numDims<2)
-// 		return "Image can only be either 2D or 3D!";
-// 
-// 	return "";
-// }
-// 
-// std::string OtsuThesholdValue::printUsage()
-// {
-// 	return "threshold = CudaMex('OtsuThesholdValue',imageIn)";
-// }
+#include "MexCommand.h"
+#include "CudaProcessBuffer.cuh"
+ 
+ void OtsuThesholdValue::execute( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[] )
+ {
+	 Vec<size_t> imageDims;
+	 HostPixelType* imageIn;
+	 CudaProcessBuffer cudaBuffer;
+	 setupImagePointers(prhs[0],&imageIn,&imageDims);
+ 
+ 	double thresh = cudaBuffer.otsuThresholdValue(imageIn,imageDims);
+ 
+ 	plhs[0] = mxCreateDoubleScalar(thresh);
+ }
+ 
+ std::string OtsuThesholdValue::check( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[] )
+ {
+ 	if (nrhs!=1)
+ 		return "Incorrect number of inputs!";
+ 
+ 	if (nlhs!=1)
+ 		return "Requires one output!";
+ 
+ 	if (!mxIsUint8(prhs[0]))
+ 		return "Image has to be formated as a uint8!";
+ 
+ 	size_t numDims = mxGetNumberOfDimensions(prhs[0]);
+ 	if (numDims>3 || numDims<2)
+ 		return "Image can only be either 2D or 3D!";
+ 
+ 	return "";
+ }
+ 
+ std::string OtsuThesholdValue::printUsage()
+ {
+ 	return "threshold = CudaMex('OtsuThesholdValue',imageIn)";
+ }
+
+ std::string OtsuThesholdValue::printHelp()
+ {
+	 std::string msg = "\tCalculates the optimal two class threshold using Otsu's method.\n";
+	 msg += "\n";
+	 return msg;
+ }
\ No newline at end of file