在图像处理中Otsu方法，是以 Nobuyuki otsu 的名字命名的（日本人，大津展之），常用于基于图像分割的聚类。该算法的理论依据是：假定图像包含两类像素（前景像素和背景像素），直方图为双峰直方图，然后计算使得两类像素能分开的最佳阈值（类内方差），或等价的间类间方差最大。

Otsu算法原理：

对于图像 I(x,y)，前景(即目标)和背景的分割阈值记作 T，属于前景的像素点数占整幅图像的比例记为 ω0，平均灰度为 μ0；背景像素点数占整幅图像的比例为 ω1，平均灰度为 μ1；整幅图像的平均灰度记为μ，类间方差记为g。
假设图像大小为M×N，图像中像素的灰度值小于阈值 T 的像素个数为 N0，像素灰度大于阈值T的像素个数为 N1，那么：
　　　　　　ω0=N0/ M×N                (1)
　　　　　　ω1=N1/ M×N        (2)
　　　　　　N0+N1=M×N                      (3)
　　　　　　ω0+ω1=1　　　                 (4)
　　　　　　μ=ω0*μ0+ω1*μ1 (5)
　　　　　　g=ω0(μ0-μ)^2+ω1(μ1-μ)^2 (6)
　　　　　　g=ω0ω1(μ0-μ1)^2 　　    (7)　
采用遍历的方法使得类间方差g最大的阈值T，即为所求。Ostu方法可以形象地理解为：求取直方图有两个峰值的图像中那两个峰值之间的低谷值 T 。

Otsu算法实现：

matlab函数：

matlab中有现成的函数实现，函数名为： graythresh，该函数便是用Ostu方法求分割阈值Ｔ。用法如下：

imgScr=imread('..');
T = graythresh(imgScr);
BW = im2bw(imgScr,T);

我这里po下源码，方便理解：

function [level em] = graythresh(I)
%GRAYTHRESH Global image threshold using Otsu's method.
%   LEVEL = GRAYTHRESH(I) computes a global threshold (LEVEL) that can be
%   used to convert an intensity image to a binary image with IM2BW. LEVEL
%   is a normalized intensity value that lies in the range [0, 1].
%   GRAYTHRESH uses Otsu's method, which chooses the threshold to minimize
%   the intraclass variance of the thresholded black and white pixels.
%
%   [LEVEL EM] = GRAYTHRESH(I) returns effectiveness metric, EM, as the
%   second output argument. It indicates the effectiveness of thresholding
%   of the input image and it is in the range [0, 1]. The lower bound is
%   attainable only by images having a single gray level, and the upper
%   bound is attainable only by two-valued images.
%
%   Class Support
%   -------------
%   The input image I can be uint8, uint16, int16, single, or double, and it
%   must be nonsparse.  LEVEL and EM are double scalars.
%
%   Example
%   -------
%       I = imread('coins.png');
%       level = graythresh(I);
%       BW = im2bw(I,level);
%       figure, imshow(BW)
%
narginchk(1,1);
validateattributes(I,{'uint8','uint16','double','single','int16'},{'nonsparse'}, ...mfilename,'I',1);if ~isempty(I)% Convert all N-D arrays into a single column.  Convert to uint8 for% fastest histogram computation.I = im2uint8(I(:));num_bins = 256;counts = imhist(I,num_bins);% Variables names are chosen to be similar to the formulas in% the Otsu paper.p = counts / sum(counts);omega = cumsum(p);mu = cumsum(p .* (1:num_bins)');mu_t = mu(end);sigma_b_squared = (mu_t * omega - mu).^2 ./ (omega .* (1 - omega));% Find the location of the maximum value of sigma_b_squared.% The maximum may extend over several bins, so average together the% locations.  If maxval is NaN, meaning that sigma_b_squared is all NaN,% then return 0.maxval = max(sigma_b_squared);isfinite_maxval = isfinite(maxval);if isfinite_maxvalidx = mean(find(sigma_b_squared == maxval));% Normalize the threshold to the range [0, 1].level = (idx - 1) / (num_bins - 1);elselevel = 0.0;end
elselevel = 0.0;isfinite_maxval = false;
end% compute the effectiveness metric
if nargout > 1if isfinite_maxvalem = maxval/(sum(p.*((1:num_bins).^2)') - mu_t^2);elseem = 0;end
end

我试了两幅图像，效果如下：

opencv函数实现：

po下opencv的源码：

    getThreshVal_Otsu_8u( const Mat& _src )  {  Size size = _src.size();  if( _src.isContinuous() )  {  size.width *= size.height;  size.height = 1;  }  const int N = 256;  int i, j, h[N] = {0};  for( i = 0; i < size.height; i++ )  {  const uchar* src = _src.data + _src.step*i;  j = 0;  #if CV_ENABLE_UNROLLED  for( ; j <= size.width - 4; j += 4 )  {  int v0 = src[j], v1 = src[j+1];  h[v0]++; h[v1]++;  v0 = src[j+2]; v1 = src[j+3];  h[v0]++; h[v1]++;  }  #endif  for( ; j < size.width; j++ )  h[src[j]]++;  }  double mu = 0, scale = 1./(size.width*size.height);  for( i = 0; i < N; i++ )  mu += i*(double)h[i];  mu *= scale;  double mu1 = 0, q1 = 0;  double max_sigma = 0, max_val = 0;  for( i = 0; i < N; i++ )  {  double p_i, q2, mu2, sigma;  p_i = h[i]*scale;  mu1 *= q1;  q1 += p_i;  q2 = 1. - q1;  if( std::min(q1,q2) < FLT_EPSILON || std::max(q1,q2) > 1. - FLT_EPSILON )  continue;  mu1 = (mu1 + i*p_i)/q1;  mu2 = (mu - q1*mu1)/q2;  sigma = q1*q2*(mu1 - mu2)*(mu1 - mu2);  if( sigma > max_sigma )  {  max_sigma = sigma;  max_val = i;  }  }  return max_val;  }

上面的应该是opencv 2以后的版本，之前的版本更好理解一些，这里也po一下，

int otsuThreshold(IplImage *frame)
{const int GrayScale = 256;int width = frame->width;int height = frame->height;int pixelCount[GrayScale];float pixelPro[GrayScale];int i, j, pixelSum = width * height, threshold = 0;uchar* data = (uchar*)frame->imageData;for (i = 0; i < GrayScale; i++){pixelCount[i] = 0;pixelPro[i] = 0;}//统计灰度级中每个像素在整幅图像中的个数  for (i = 0; i < height; i++){for (j = 0; j < width; j++){pixelCount[(int)data[i * width + j]]++; }}//计算每个像素在整幅图像中的比例  float maxPro = 0.0;int kk = 0;for (i = 0; i < GrayScale; i++){pixelPro[i] = (float)pixelCount[i] / pixelSum;if (pixelPro[i] > maxPro){maxPro = pixelPro[i];kk = i;}}//遍历灰度级[0,255]， i作为阈值  float w0, w1, u0tmp, u1tmp, u0, u1, u, deltaTmp, deltaMax = 0;for (i = 0; i < GrayScale; i++)    {w0 = w1 = u0tmp = u1tmp = u0 = u1 = u = deltaTmp = 0;for (j = 0; j < GrayScale; j++){if (j <= i)   //背景部分  {w0 += pixelPro[j];u0tmp += j * pixelPro[j];}else        //前景部分  {w1 += pixelPro[j];u1tmp += j * pixelPro[j];}}u0 = u0tmp / w0;u1 = u1tmp / w1;u = u0tmp + u1tmp;deltaTmp = w0 * pow((u0 - u), 2) + w1 * pow((u1 - u), 2);if (deltaTmp > deltaMax){deltaMax = deltaTmp;threshold = i;}}return threshold;
}

参考：

https://en.wikipedia.org/wiki/Otsu%27s_method

http://www.cnblogs.com/ranjiewen/p/6385564.html

http://blog.csdn.net/glouds/article/details/38976573