otsu阈值分割

#otsu阈值分割
import sys,cv2,math
import numpy as np
import matplotlib.pyplot as plt
def calcGrayHist(image):
    rows,cols = image.shape
    grayHist = np.zeros([256],np.uint64)
    for r in range(rows):
        for c in range(cols):
            grayHist[image[r][c]] +=1#把图像灰度值作为索引
    return(grayHist)
def otsu(image):
    rows,cols = image.shape
    #灰度直方图
    grayHist = calcGrayHist(image)
    #归一化灰度直方图
    uniformGrayHist = grayHist/float(rows*cols)
    #计算零阶累积矩和一阶累积矩
    zeroCumuMoment = np.zeros([256],np.float32)
    oneCumuMoment = np.zeros([256],np.float32)
    for k in range(256):
        if k==0:
            zeroCumuMoment[k] = uniformGrayHist[0]
            oneCumuMoment[k] = (k)*uniformGrayHist[0]
        else:
            zeroCumuMoment[k] = zeroCumuMoment[k-1]+uniformGrayHist[k]
            oneCumuMoment[k] = oneCumuMoment[k-1]+k*uniformGrayHist[k]
    #计算类间方差
    variance = np.zeros([256],np.float32)
    for k in range(255):
        if zeroCumuMoment[k] == 0 or zeroCumuMoment[k] ==1:
            variance[k] = 0
        else:
            variance[k]=math.pow(oneCumuMoment[255]*zeroCumuMoment[k]-
           oneCumuMoment[k],2)/(zeroCumuMoment[k]*(1.0-zeroCumuMoment[k]))
    threshLoc = np.where(variance[0:255]==np.max(variance[0:255]))
    thresh = threshLoc[0][0]
    #
    threshold = np.copy(image)
    threshold[threshold>thresh]=255
    threshold[threshold<=thresh]=0
    return(thresh,threshold)
if __name__ =='__main__':
    src = cv2.imread('E:/sy2/6/img.jpg',cv2.IMREAD_GRAYSCALE)   
    re,ra = otsu(src)
    print(re)
    print(ra)
    cv2.imshow('ra',ra)
    cv2.waitKey(0)
    cv2.destroyAllWindows()