OpenCV-Otsu阈值分割

步骤：1、计算累加直方图；2、计算一阶累积矩；3、计算图像的总体灰度平均值；4、计算每个灰度值对应的方差；5、找到目标阈值。

C++：

#include<opencv2/opencv.hpp>
#include<iostream>
using namespace std;
using namespace cv;
//计算图像灰度直方图
Mat calcgrayhist(const Mat&image)
{
	Mat histogram = Mat::zeros(Size(256, 1), CV_32SC1);
	//图像宽高
	int rows = image.rows;
	int cols = image.cols;
	for (int i = 0; i < rows; i++)
	{
		for (int j = 0; j < rows; j++)
		{
			int index = int(image.at<uchar>(i, j));
			histogram.at<int>(0, index) += 1;
		}
	}
	return histogram;
};

//OTSU
int OTSU(const Mat&image,Mat &OTSU_image)
{	
	int rows = image.rows;
	int cols = image.cols;
	Mat histogram = calcgrayhist(image);
	//归一化直方图
	Mat normhist;
	histogram.convertTo(normhist, CV_32FC1, 1.0 / (rows*cols), 0.0);
	//计算累加直方图和一阶累积矩
	Mat zeroaccumulate = Mat::zeros(Size(256, 1), CV_32FC1);
	Mat oneaccumulate = Mat::zeros(Size(256, 1), CV_32FC1);
	for (int i = 0; i < 256; i++)
	{
		if (i == 0)
		{
			zeroaccumulate.at<float>(0, i) = normhist.at<float>(0, i);
			oneaccumulate.at<float>(0, i) =i * normhist.at<float>(0, i);
		}
		else
		{
			zeroaccumulate.at<float>(0, i) = zeroaccumulate.at<float>(0, i-1)
										+normhist.at<float>(0, i);
			oneaccumulate.at<float>(0, i) = oneaccumulate.at<float>(0, i - 1)
				+i* normhist.at<float>(0, i);
		}
	}
	//计算间类方差
	Mat variance= Mat::zeros(Size(256, 1), CV_32FC1);
	float mean = oneaccumulate.at<float>(0, 255);
	for (int i = 0; i < 255; i++)
	{
		if (zeroaccumulate.at<float>(0, i) == 0 || zeroaccumulate.at<float>(0, i) == 1)
			variance.at<float>(0, i) = 0;
		else
		{
			float cofficient = zeroaccumulate.at<float>(0, i)*(1.0 -
				zeroaccumulate.at<float>(0, i));
			variance.at<float>(0, i) = pow(mean*zeroaccumulate.at<float>(0, i)
				- oneaccumulate.at<float>(0, i), 2.0) / cofficient;
		}
	}
	//找到阈值；
	Point maxloc;
	//计算矩阵中最大值
	minMaxLoc(variance, NULL, NULL, NULL, &maxloc);
	int otsuthreshold = maxloc.x;
	threshold(image, OTSU_image, otsuthreshold, 255, THRESH_BINARY);
	return otsuthreshold;
};

int main()
{
	Mat img, gray_img,dst;
	img = imread("D:/testimage/orange.jpg");
	cvtColor(img, gray_img, COLOR_BGR2GRAY);
	OTSU(gray_img,dst);
	imshow("otsu image", dst);
	waitKey(0);
	return 0;
};

结果：
Python：

import cv2 as cv
import numpy as np
import math

def calc_grayhist(image):
    #图像宽高
    rows,cols=image.shape
    grayhist=np.zeros([256],np.uint64)
    for i in range(rows):
        for j in range(cols):
            grayhist[image[i][j]]+=1
    return grayhist

def OTSU(image):
    rows, cols = image.shape
    grayhist=calc_grayhist(image)
    #归一化直方图
    uniformgrayhist=grayhist/float(rows*cols)
    #计算零阶累积矩和一阶累积矩
    zeroaccumulat = np.zeros([256],np.float32)
    oneaccumulat = np.zeros([256], np.float32)
    for k in range(256):
        if k==0:
            zeroaccumulat[k]=uniformgrayhist[0]
            oneaccumulat[k]=k*uniformgrayhist[0]
        else:
            zeroaccumulat[k]=zeroaccumulat[k-1]+uniformgrayhist[k]
            oneaccumulat[k]=oneaccumulat[k-1]+k*uniformgrayhist[k]
    #计算间类方差
    variance=np.zeros([256],np.float32)
    for k in range(256):
        if zeroaccumulat[k]==0 or zeroaccumulat[k]==1:
            variance[k]=0
        else:
            variance[k]=math.pow(oneaccumulat[255]*zeroaccumulat[k]-
                  oneaccumulat[k],2)/(zeroaccumulat[k]*(1.0-zeroaccumulat[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 threshold

if __name__=="__main__":
    img=cv.imread("D:/testimage/cxk.jpg")
    gray_dst = cv.cvtColor(img, cv.COLOR_BGR2GRAY)
    Otsu_dst=OTSU(gray_dst)
    cv.imshow("gray dst",gray_dst)
    cv.imshow("Otsu dst", Otsu_dst)
    cv.imwrite("D:/testimage/result-cxk.jpg",Otsu_dst)
    cv.waitKey(0)
    cv.destroyAllWindows()