利用OpenCV 2.2 的Python接口实现Ostu（大津法）获取阈值

OpenCV（版本2.2）是啥不说了你懂的。其实本文只是使用OpenCV计算了一个直方图。完全可以写Python代码实现。

Python（版本是2.7）是啥也不说你也懂的。



什么是大津法？

请注意是大“津”法，不是大“律”法！我第一看到的时候阴差阳错的看成了大“律”法，然后还在百度和谷歌上溜了一圈，居然还找到了很多资料，后来在琢磨Ostu四个字母的时候发现，Ostu是日本的一个县，叫大津县，所以这里应当是“大津法”。



大津法，又叫最大类间方差法，简称Otsu，是由日本学者大津于1979年提出的，是一种自适应的阈值确定方法。按图像的灰度特性，将图像分成两个部分，称之为部分0和部分1（即通常意义的背景与目标），两个部分的类间方差越大，说明构成图像的两部分的差别越大，也就是说图片的背景和目标的差别越大。



对于图像f(x,y)，通过阈值t将图片分成两个部分，亮度小于t的称之为部分0，亮度大于等于t的称之为部分1。部分0像素点占整幅图像的比例为ω0，平均灰度μ0；部分1像素点占整幅图像的比例为ω1，平均灰度为μ1。图像的总平均灰度记为μ，类间方差记为g。

则有

μ = ω0μ0+ω1μ1

g = ω0(μ0-μ)2+ω1(μ1-μ)2



将第一个式子带入到第二个式子我们可以得到

g=ω0ω1(μ0-μ1)2

通过遍历所有的亮度（一般是[0,255]）求得所有的t值对应的g，最大的g出现的t就是Otsu算法得到的最佳阈值。



下面的Python程序展示了在OpenCV 2.2中使用Python 2.7计算Otsu最佳阈值，为了方便阅读，所有的ω使用w，μ使用u代替。

import cv

def OtsuGray( grayImage ,debug = 0):
    # 如果图片是Mat对象，则转换为Image对象
    if type(grayImage) == cv.cvmat:
        grayImage = cv.GetImage(grayImage)
        
    # 创建Hist
    hist = cv.CreateHist([256],cv.CV_HIST_ARRAY,[[0,256]])
    cv.ClearHist(hist)
    # 计算Hist
    cv.CalcHist([grayImage],hist)

    # 开始计算
    # 计算总亮度
    totalH = 0
    for h in range(0,256):
        v = cv.QueryHistValue_1D(hist,h)
        if v == 0 : continue
        totalH += v*h
        if debug > 3 : print "t=%d,%d,%d"%(h,totalH,v*h)
        

    width  = grayImage.width
    height = grayImage.height
    total  = width*height
    
    if debug > 1 : print "总像素:%d;总亮度:%d平均亮度:%0.2f"%(total,totalH,totalH/total)

    # t=0和t=255的时候无法构成分割，所以从t=1开始计算一致到t=255
    # 初始化v值
    v = 0

    gMax = 0.0
    tIndex = 0
    
    # temp
    n0Acc = 0
    n1Acc = 0
    n0H   = 0
    n1H   = 0
    for t in range(1,255):
        v = cv.QueryHistValue_1D(hist,t-1)
        if v == 0: continue
        
        n0Acc += v          #灰度小于t的像素的数目
        n1Acc = total - n0Acc #灰度大于等于t的像素的数目
        n0H += (t-1)*v          #灰度小于t的像素的总亮度
        n1H = totalH - n0H  #灰度大于等于t的像素的总亮度

        if n0Acc > 0 and n1Acc > 0:
            u0 = n0H/n0Acc # 灰阶小于t的平均灰度
            u1 = n1H/n1Acc # 灰阶大于等于t的平均灰度
            w0 = n0Acc/total # 灰阶小于t的像素比例
            w1 = 1.0-w0      # 灰阶大于等于t的像素的比例
            uD = u0-u1
            g = w0 * w1 * uD * uD

            if debug > 2: print "t=%3d; u0=%.2f,u1=%.2f,%.2f;n0H=%d,n1H=%d; g=%.2f"/
               %(t,u0,u1,u0*w0+u1*w1,n0H,n1H,g)

            if gMax < g:
                gMax   = g
                tIndex = t
        
    if debug >0 : print "gMaxValue=%.2f; t = %d ; t_inv = %d"/
               %(gMax,tIndex,255-tIndex)            

    return tIndex

展示一下二值化的效果：

原图如下：



t=102时g取得最大值



为了节省空间，只传这么一个例子吧~

各位老师，文章中若有错误之处，恳请您斧正。