【图像处理】:Otsu算法最大类间方差法（大津算法：附源码）

这里写自定义目录标题

数学原理
算法评价
参考链接

数学原理

以灰度图像为例，对于图像M×N大小的矩阵，即图像中的像素，每一个值即为像素值，其中灰度图像像素值在(0~255)之间。
主要实现前景(即目标)和背景的分割：
主要公式：
前景的像素点数占整幅图像的比例记为ω0，前景平均灰度记为μ0
背景像素点数占整幅图像的比例记为ω1，其平均灰度记为μ1
图像的总平均灰度记为μ，类间方差记为maximum。
假设图像的背景较暗，并且图像的大小为M×N，图像中像素的灰度值小于阈值optimal threshold的像素个数记作N0
，像素灰度大于等于阈值optimalthreshold 的像素个数记作N1，
则有：
　　　　　　 ω0 = N0 / ( M × N ) 　　　　　　 (1) 　　
　　　　　　 ω1 = N1 / ( M × N ) 　　　　　　 (2) 　　　　　　
　　　　 N0 + N1 = M × N 　　　　　　 (3) 　　
　　　　　　 1 = ω 0 + ω 1 　　　　　　 (4) 　　　　　
　　　　　　 μ = ω0 × μ0 + ω1 × μ1 　　　　 (5) 　　
　　　　maximum = ω0 × ( μ0 − μ ) 2 + ω1 × ( μ1 − μ ) 2 (6) 　
将式(5)代入式(6)，得到等价公式(7)：
　　　 maximum = ω0 × ω1 × (μ0 − μ1 ) 2 (7) 　　　
采用遍历的方法得到使类间方差maximum最大的阈值optimal threshold
实现过程：

#!/usr/bin/env python
# -*- coding:utf-8 -*-
# author:longc
# datetime:2023/11/16 10:30
# software: PyCharm
# function: 图像处理逻辑

import cv2
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image


# otsu算法
def otsu(gray):
    pixel_number = gray.shape[0] * gray.shape[1]
    mean_weigth = 1.0 / pixel_number
    # #统计各灰度级的像素个数，灰度级分为256级
    # bins必须写到257，否则255这个值只能分到[254,255)区间
    his, bins = np.histogram(gray, np.arange(0, 257))  # 计算灰度的直方图，计数统计区间为0-257
    print("bins", bins)
    print("his", his)
    # 绘制直方图
    plt.figure(figsize=(12, 8))
    # plt.hist(gray, 256, [0, 256], label='灰度级直方图')  # 运行比较慢，如果电脑卡顿，可以将本行代码注释掉
    plt.show()

    final_thresh = -1
    final_value = -1
    intensity_arr = np.arange(256)  # 灰度分为256级，0级到255级

    # ************************************************************ 采用遍历的方法得到类间方差最大的阈值
    
    for t in bins[1:-1]:  # 遍历1到254级 (一定不能有超出范围的值)
        pcb = np.sum(his[:t])  # 小于当前灰度对应的所有像素点计数
        pcf = np.sum(his[t:])  # 大于当前灰度对应的所有像素点计数
        Wb = pcb * mean_weigth  # 像素被分类为背景的概率
        Wf = pcf * mean_weigth  # 像素被分类为目标的概率
        # if t == 100:
        #     print("1>>>", intensity_arr[:t])
        #     print("2>>>", his[:t])
        #     print("3>>>", np.sum(intensity_arr[:t] * his[:t]))
        #     print("4>>>", float(pcb))
        #     print("5>>>", np.sum(intensity_arr[:t] * his[:t]) / float(pcb))

        mub = np.sum(intensity_arr[:t] * his[:t]) / float(pcb)  # 分类为背景的像素均值
        muf = np.sum(intensity_arr[t:] * his[t:]) / float(pcf)  # 分类为目标的像素均值

        # print mub, muf
        value = Wb * Wf * (mub - muf) ** 2  # 计算目标和背景类间方差
        # 采用遍历的方法得到使类间方差value最大的阈值final_value和二值化对应最大的final_thresh
        
        if value > final_value:
            final_thresh = t  # 进行二值化的操作值
            final_value = value
    print("final_thresh>>>", final_thresh)
    print("final_value>>>", final_value)
    
    # 二值化操作处理
    # final_img = gray.copy()
    # print(final_thresh)
    # final_img[gray > final_thresh] = 255
    # final_img[gray < final_thresh] = 0
    # cv2.imwrite("final_img.jpg", final_img)
    plt.imshow(gray)
    plt.show()

    # 二值化图像（多种方法对比）
    ret, binary_image = cv2.threshold(gray, final_thresh-15, 255, cv2.THRESH_BINARY)
    plt.imshow(binary_image, cmap='gray')
    plt.show()

    # ret, binary_image1 = cv2.threshold(gray, final_thresh, 255, cv2.THRESH_TRUNC)
    # plt.imshow(binary_image1)
    # plt.show()
    #
    # ret, binary_image2 = cv2.threshold(gray, final_thresh, 255, cv2.THRESH_TOZERO)
    # plt.imshow(binary_image2)
    # plt.show()
    #
    # ret, binary_image3 = cv2.threshold(gray, final_thresh, 255, cv2.THRESH_TOZERO_INV)
    # plt.imshow(binary_image3)
    # plt.show()


imggray = cv2.imread("IMG_0004_3.jpg", 0)
plt.title("imggray")
plt.imshow(imggray, cmap='gray')
plt.show()

# 进行OSTU运算
otsu(imggray)