一、Otsu 的基本思想（类间方差最大）

假设一幅灰度图中只有两类像素：

前景：目标（如文字、物体）
背景：不是目标的部分

我们要找一个灰度阈值 tt，把图像分成两类：

类 0（背景）：灰度 [0,t][0,t]
类 1（前景）：灰度 [t+1,L−1][t+1,L−1]，其中 LL 是灰度级数（一般 256）

Otsu 的思想：

找到那个使类间方差（两类之间差异）最大的阈值。

二、数学公式（单阈值 Otsu）

设灰度级为 0,1,…,L−10,1,…,L−1，直方图为 h(i)h(i)，总像素数为 NN。

概率分布

pi=h(i)N,i=0,…,L−1pi=Nh(i),i=0,…,L−1

给定阈值 tt，两类像素的概率（权重）：

ω0(t)=∑i=0tpi,ω1(t)=∑i=t+1L−1piω0(t)=i=0∑tpi,ω1(t)=i=t+1∑L−1pi

两类的平均灰度：

μ0(t)=1ω0(t)∑i=0ti pi,μ1(t)=1ω1(t)∑i=t+1L−1i piμ0(t)=ω0(t)1i=0∑tipi,μ1(t)=ω1(t)1i=t+1∑L−1ipi

图像整体平均灰度：

μT=∑i=0L−1i piμT=i=0∑L−1ipi

类间方差（between-class variance）定义为：

σB2(t)=ω0(t) (μ0(t)−μT)2+ω1(t) (μ1(t)−μT)2σB2(t)=ω0(t)(μ0(t)−μT)2+ω1(t)(μ1(t)−μT)2

常用等价形式：

σB2(t)=ω0(t) ω1(t)(μ0(t)−μ1(t))2σB2(t)=ω0(t)ω1(t)(μ0(t)−μ1(t))2

算法目标：

t=arg\max_tσ_B^2(t)

三、Otsu 算法实现步骤（单阈值）

计算图像灰度直方图 h(i)h(i)，并归一化得到 pipi
利用前缀和等方式快速计算：
- ω0(t)ω0(t)、ω1(t)ω1(t)
- μ0(t)μ0(t)、μ1(t)μ1(t) 或直接用累积和和 μTμT 计算
对每一个可能的阈值 tt（0 到 L-2，一般 0 到 254）：
- 计算对应的类间方差 σB2(t)σB2(t)
取使类间方差最大的 tt 作为 Otsu 阈值
按此阈值对图像进行二值化：
g(x,y) = \begin{cases} 0, & \text{if } f(x,y) \le t, \\[4pt] 255, & \text{if } f(x,y) > t . \end{cases}

四、Matlab 示例代码

1. 使用 Matlab 自带 Otsu（`graythresh` + `imbinarize`）

% 读入图像（灰度图）
I = imread('input.png');
if size(I, 3) == 3
    I = rgb2gray(I);   % 若为彩色图像，转为灰度
end

% Matlab 的 Otsu 阈值：graythresh 返回的是 [0,1] 之间的归一化阈值
level = graythresh(I);

% 二值化
BW = imbinarize(I, level);

% 显示结果
figure;
subplot(1,2,1); imshow(I);  title('Original');
subplot(1,2,2); imshow(BW); title(['Otsu (graythresh) level = ', num2str(level)]);

graythresh 的实现就是 Otsu 单阈值算法。注意它返回的是归一化阈值（相对于 0–1），实际灰度阈值为 T = level * 255。

2. 手写 Otsu 算法（示例实现）

下面这段代码不依赖 graythresh，按公式自己算一遍阈值。

% 读入图像
I = imread('input.png');
if size(I,3) == 3
    I = rgb2gray(I);
end

% 转 double 方便计算（也可以不转）
I = double(I);

% 参数
L = 256;              % 灰度级数 0~255
N = numel(I);         % 像素总数

% 1. 计算直方图 h(i)
[counts, ~] = imhist(uint8(I), L);   % counts: 256x1
p = counts / N;                      % 概率分布 p_i

% 2. 累计概率 ω(k) 和 累计均值 μ(k)
omega = cumsum(p);                   % ω_k
mu_k  = cumsum((0:L-1)' .* p);       % μ_k 的分子部分
mu_T  = mu_k(end);                   % 总均值 μ_T

% 3. 计算每个可能阈值 t 的类间方差 σ_B^2(t)
sigma_b2 = zeros(L, 1);
for t = 1:L
    if omega(t) == 0 || omega(t) == 1
        sigma_b2(t) = 0;
    else
        sigma_b2(t) = (mu_T * omega(t) - mu_k(t))^2 / (omega(t) * (1 - omega(t)));
    end
end

% 4. 找到最大类间方差对应的阈值 t*
[~, idx] = max(sigma_b2);
T = idx - 1;      % Matlab 索引从 1 开始，对应灰度 0~255，所以减 1

fprintf('Otsu threshold = %d\n', T);

% 5. 按 T 进行二值化
BW_manual = I > T;

% 显示结果
figure;
subplot(1,3,1); imshow(uint8(I));      title('Original');
subplot(1,3,2); imshow(BW_manual);     title(['Manual Otsu, T = ', num2str(T)]);
subplot(1,3,3); imhist(uint8(I));      hold on;
line([T T], ylim, 'LineWidth', 2);     title('Histogram with Otsu T');

五、Otsu 算法的特点

优点

无参数：不需要手工指定阈值，自动从数据中估计；
实现简单：只涉及直方图、累积和，计算高效；
对于前景和背景灰度分布差异较明显、近似双峰的图像效果很好；
被广泛集成到 OpenCV、Matlab 等库中，工程实践非常常见。

局限性

默认是单阈值二类分割，如果图像存在多类（多峰直方图）会不够理想；
假设背景和前景在直方图上区分明显，如果图像对比度很低或噪声很多，效果会下降；
标准 Otsu 是全局阈值，对光照不均匀的图像效果一般（可用局部/自适应方法改进）。

六、常见扩展与变体

多阈值 Otsu（Multi-level Otsu）
不只找一个阈值，而是找多个阈值 t1,t2,…t1,t2,…，将灰度分成多类，每个类之间的类间方差总和最大。
计算复杂度会显著增加（阈值组合数量多），常配合搜索优化或启发式算法。
二维 Otsu（2D Otsu）
不只考虑像素本身的灰度，还加入邻域平均、梯度等信息，用二维直方图计算最佳阈值，能更好对抗噪声，但运算量更大。
局部/块级 Otsu
将图像分成小块，对每块单独使用 Otsu 得到局部阈值，适用于光照不均、阴影较重的场景。
与其它方法组合
如先做滤波、直方图均衡或对数变换，再用 Otsu；或 Otsu 与形态学处理一起使用，提高分割质量。