二维OTSU算法快速实现

        OTSU自适应阈值选取算法作为一种经典的阈值分割算法，在图像领域有广泛的应用。在一维OTSU算法上发展而来的二维OTSU算法因为计算时间复杂度高而制约了其应用。通过消除二维自适应阈值算法中的冗余计算，用迭代的方式得到查询表，从而大大的提高了二维阈值计算的速度。

        实现的代码：

// 二维OTSU算法实现.cpp : 定义控制台应用程序的入口点。

#include "stdafx.h"
#include <opencv2\opencv.hpp>
#include <highgui.h>
#include <cv.h>
#include <iostream> 
#include <stdio.h>
#include<opencv\highgui.h>
#include<opencv2\core\core.hpp>
#include "math.h"

int OTSU2d(IplImage * src)
{
	int height = src->height;
	int width = src->width;
	long pixel = height * width;
	int histogram[256][256];
	int i,j;

	for (i = 0;i < 256;i++)//初始化直方图
	{
		for (j = 0; j < 256;j++)
			histogram[i][j] = 0;
	}

	IplImage * temp = cvCreateImage(cvGetSize(src),8,1);
	cvSmooth(src,temp,CV_BLUR,3,0);
	
	for (i = 0;i < height;i++)//计算直方图
	{
		for (j = 0; j < width;j++)
		{
			int data1 = cvGetReal2D(src,i,j);
			int data2 = cvGetReal2D(temp,i,j);
			histogram[data1][data2]++;
		}
	}

	double p_histogram[256][256];
	for (i = 0; i < 256;i++)//直方图归一化
		for(j = 0; j < 256;j++)
			p_histogram[i][j] = (histogram[i][j]*1.0)/(pixel*1.0);

	double Pst_0[256][256];//Pst_0用来存储概率分布情况
	Pst_0[0][0] = p_histogram[0][0];
	for(i = 0;i < 256;i++)//计算概率分布情况
		for(j = 0;j < 256;j++)
		{
		   double temp = 0.0;
		   if(i-1 >= 0)
			   temp = temp + Pst_0[i-1][j];
		   if(j-1 >= 0)
			   temp = temp + Pst_0[i][j-1];
		   if(i-1 >= 0 && j-1 >= 0)
			   temp = temp - Pst_0[i-1][j-1];
		   temp = temp + p_histogram[i][j];
		   Pst_0[i][j] = temp;
		}
		
	double Xst_0[256][256];//存储x方向上的均值矢量
	Xst_0[0][0] = 0 * Pst_0[0][0];
	for(i = 0 ; i < 256;i++)//计算x方向上的均值矢量
		for(j = 0 ; j < 256;j++)
		{
		   double temp = 0.0;
		   if(i-1 >= 0)
			   temp = temp + Xst_0[i-1][j];
		   if(j-1 >= 0)
			   temp = temp + Xst_0[i][j-1];
		   if(i-1 >= 0 && j-1 >= 0)
			   temp = temp - Xst_0[i-1][j-1];
		   temp = temp + i * p_histogram[i][j];
		   Xst_0[i][j] = temp;
		}

	double Yst_0[256][256];//存储y方向上的均值矢量
	Yst_0[0][0] = 0 * Pst_0[0][0];
	for(i = 0 ; i < 256;i++)//计算y方向上的均值矢量
		for(j = 0 ; j < 256;j++)
		{
		   double temp = 0.0;
		   if(i-1 >= 0)
			   temp = temp + Yst_0[i-1][j];
		   if(j-1 >= 0)
			   temp = temp + Yst_0[i][j-1];
		   if(i-1 >= 0 && j-1 >= 0)
			   temp = temp - Yst_0[i-1][j-1];
		   temp = temp + j * p_histogram[i][j];
		   Yst_0[i][j] = temp;
		}
	 
	int threshold1;
	int threshold2;
	double variance = 0.0;
	double maxvariance = 0.0;
	for(i = 0;i < 256;i++)//计算类间离散测度
		for(j = 0;j < 256;j++)
		{
			 long double p0 = Pst_0[i][j];
			 long double v0 = pow(((Xst_0[i][j]/p0)-Xst_0[255][255]),2) + pow(((Yst_0[i][j]/p0)-Yst_0[255][255]),2);
			 long double p1 = Pst_0[255][255]-Pst_0[255][j]-Pst_0[i][255]+Pst_0[i][j];
			 long double vi = Xst_0[255][255]-Xst_0[255][j]-Xst_0[i][255]+Xst_0[i][j];
			 long double vj = Yst_0[255][255]-Yst_0[255][j]-Yst_0[i][255]+Yst_0[i][j];
			 long double v1 = pow(((vi/p1)-Xst_0[255][255]),2)+pow(((vj/p1)-Yst_0[255][255]),2);
			
			 variance = p0*v0+p1*v1;
		   
		   if(variance > maxvariance)
		   {
		      maxvariance = variance;
			  threshold1 = i;
			  threshold2 = j;
		   }
		}

	//printf("%d  %d",threshold1,threshold2);

	return (threshold1+threshold2)/2;

}

int _tmain(int argc, _TCHAR* argv[])
{
	const char * filename  = "rice.jpg";
	IplImage * img = cvLoadImage(filename,0);
	IplImage * dst = cvCreateImage(cvGetSize(img),8,1);
	int threshold = OTSU2d(img);
	cvThreshold(img,dst,threshold,255,CV_THRESH_BINARY);
	cvNamedWindow( "img", 1 );  
	cvShowImage("img", dst);  
	cvWaitKey(-1);
	return 0;
}

 

二值化前：
     

二值化后：