Otsu方法（大津法）是一种广泛使用的自动图像阈值分割技术，它通过最大化类间方差来确定最佳阈值。
Otsu方法原理
Otsu方法的基本思想是将图像像素分为前景和背景两类，通过寻找使两类间方差最大的阈值来实现最佳分割。
数学原理
设图像有L个灰度级（通常为0-255）
设阈值为T，将像素分为两类：
- C₀: 灰度值在[0, T]之间的像素
- C₁: 灰度值在[T+1, L-1]之间的像素
计算两类间的方差：
σ²_b(T) = ω₀(μ₀ - μ_T)² + ω₁(μ₁ - μ_T)²其中：
- ω₀, ω₁: 两类像素的概率
- μ₀, μ₁: 两类像素的平均灰度
- μ_T: 整幅图像的平均灰度
寻找使σ²_b(T)最大的T值作为最佳阈值
MATLAB中的Otsu阈值实现
MATLAB提供了内置函数graythresh来实现Otsu方法，该函数返回归一化的阈值（范围0-1）。
基本使用示例
% 读取图像
I = imread('coins.png');
figure;
subplot(2, 2, 1);
imshow(I);
title('原始图像');
% 转换为灰度图像（如果必要）
if size(I, 3) == 3
I_gray = rgb2gray(I);
else
I_gray = I;
end
% 使用graythresh计算Otsu阈值
level = graythresh(I_gray);
fprintf('Otsu阈值 (0-1范围): %.4f\n', level);
fprintf('Otsu阈值 (0-255范围): %.4f\n', level * 255);
% 应用阈值进行分割
BW = imbinarize(I_gray, level);
% 显示结果
subplot(2, 2, 2);
imshow(BW);
title('Otsu阈值分割结果');
% 显示灰度直方图和阈值位置
subplot(2, 2, 3);
imhist(I_gray);
hold on;
line([level*255, level*255], ylim, 'Color', 'r', 'LineWidth', 2);
title('灰度直方图与阈值');
hold off;
% 显示原始图像与分割结果的叠加
subplot(2, 2, 4);
imshow(I);
hold on;
% 创建分割边界覆盖图
boundaries = bwboundaries(BW);
for k = 1:length(boundaries)
boundary = boundaries{
k};
plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 1.5);
end
title('分割边界叠加');
hold off;
自定义Otsu阈值函数实现
为了更好地理解Otsu方法，我们可以自己实现一个Otsu阈值计算函数：
function [threshold, max_var] = otsu_threshold(image)
% OTSU_THRESHOLD 自定义Otsu阈值计算函数
% 输入: image - 灰度图像矩阵
% 输出: threshold - 计算得到的最佳阈值
% max_var - 最大类间方差
% 计算图像直方图
[counts, binLocations] = imhist(image);
% 归一化直方图，得到概率分布
p = counts / sum(counts);
% 初始化变量
max_var = 0;
threshold = 0;
total_mean = sum((0:255)' .* p); % 整体均值
% 遍历所有可能的阈值
for T = 1:256
% 计算两类概率
w0 = sum(p(1:T));
w1 = sum(p(T+1:end));
if w0 == 0 || w1 == 0
continue;
end
% 计算两类均值
mu0 = sum((0:T-1)' .* p(1:T)) / w0;
mu1 = sum((T:255)' .* p(T+1:end)) / w1;
% 计算类间方差
between_var = w0 * w1 * (mu0 - mu1)^2;
% 更新最大方差和阈值
if between_var > max_var
max_var = between_var;
threshold = T - 1; % 转换为0-255范围的阈值
end
end
end
使用自定义函数：
% 使用自定义Otsu函数
[custom_threshold, max_var] = otsu_threshold(I_gray);
fprintf('自定义Otsu阈值: %.4f\n', custom_threshold);
fprintf('最大类间方差: %.4f\n', max_var);
% 应用自定义阈值
BW_custom = I_gray > custom_threshold;
% 比较MATLAB内置函数和自定义函数的结果
figure;
subplot(1, 2, 1);
imshow(BW);
title(['MATLAB内置函数: ', num2str(level*255)]);
subplot(1, 2, 2);
imshow(BW_custom);
title(['自定义函数: ', num2str(custom_threshold)]);
多阈值Otsu方法
对于更复杂的图像，可能需要使用多阈值Otsu方法：
function thresholds = multi_otsu(image, n)
% MULTI_OTSU 多阈值Otsu方法
% 输入: image - 灰度图像
% n - 阈值数量（分割区域数为n+1）
% 输出: thresholds - 阈值向量
% 计算直方图
[counts, ~] = imhist(image);
p = counts / sum(counts);
% 初始化
thresholds = zeros(1, n);
max_var = 0;
% 生成所有可能的阈值组合（简化版本，实际应用中可能需要优化算法）
% 这里使用递归方法搜索最佳阈值组合
thresholds = recursive_otsu(p, 1, 256, n, []);
end
function best_thresholds = recursive_otsu(p, start, finish, n, current)
% 递归辅助函数
if n == 0
% 计算当前阈值组合的类间方差
var = calculate_variance(p, current);
best_thresholds = {
current, var};
return;
end
best_var = -1;
best_combination = [];
for t = start:finish-n
new_current = [current, t];
[combination, var] = recursive_otsu(p, t+1, finish, n-1, new_current);
if var > best_var
best_var = var;
best_combination = combination;
end
end
best_thresholds = best_combination;
end
function var = calculate_variance(p, thresholds)
% 计算给定阈值组合的类间方差
thresholds = sort([0, thresholds, 256]);
var = 0;
for i = 1:length(thresholds)-1
start = thresholds(i) + 1;
finish = thresholds(i+1);
if start > finish
continue;
end
w = sum(p(start:finish));
if w == 0
continue;
end
mu = sum((start-1:finish-1)' .* p(start:finish)) / w;
global_mu = sum((0:255)' .* p);
var = var + w * (mu - global_mu)^2;
end
end
参考 MATLAB实现图像分割otsuf 源程序代码 www.youwenfan.com/contentalf/23056.html
实际应用示例：细胞图像分割
% 细胞图像分割示例
function cell_segmentation_example()
% 读取细胞图像
I = imread('cell_image.png');
% 转换为灰度图
if size(I, 3) == 3
I_gray = rgb2gray(I);
else
I_gray = I;
end
% 应用高斯滤波去噪
I_filtered = imgaussfilt(I_gray, 1);
% 使用Otsu阈值分割
level = graythresh(I_filtered);
BW = imbinarize(I_filtered, level);
% 后处理：去除小对象和填充孔洞
BW_cleaned = bwareaopen(BW, 50); % 移除面积小于50像素的对象
BW_cleaned = imfill(BW_cleaned, 'holes'); % 填充孔洞
% 显示结果
figure;
subplot(2, 2, 1);
imshow(I);
title('原始细胞图像');
subplot(2, 2, 2);
imshow(I_filtered);
title('滤波后的图像');
subplot(2, 2, 3);
imshow(BW);
title('Otsu分割结果');
subplot(2, 2, 4);
imshow(BW_cleaned);
title('后处理结果');
% 计算并显示细胞统计信息
labeled_image = bwlabel(BW_cleaned);
stats = regionprops(labeled_image, 'Area', 'Centroid');
% 在图像上标记细胞
figure;
imshow(I);
hold on;
for i = 1:length(stats)
text(stats(i).Centroid(1), stats(i).Centroid(2), ...
num2str(i), 'Color', 'r', 'FontSize', 12, 'FontWeight', 'bold');
end
title(['检测到的细胞数量: ', num2str(length(stats))]);
hold off;
% 显示细胞面积分布
areas = [stats.Area];
figure;
histogram(areas, 20);
xlabel('细胞面积 (像素)');
ylabel('频率');
title('细胞面积分布');
end
Otsu方法的局限性及改进
虽然Otsu方法在许多情况下表现良好，但它也有一些局限性：
- 对噪声敏感：图像中的噪声会影响直方图形状，导致阈值选择不准确
- 双峰假设：Otsu方法假设直方图是双峰的，对于单峰或多峰直方图效果可能不佳
- 光照不均匀：在光照不均匀的图像中，全局阈值可能不适用
改进方法
% 自适应Otsu阈值处理
function adaptive_otsu_example()
I = imread('uneven_lighting.jpg');
I_gray = rgb2gray(I);
% 全局Otsu阈值
global_thresh = graythresh(I_gray);
BW_global = imbinarize(I_gray, global_thresh);
% 自适应Otsu阈值
BW_adaptive = imbinarize(I_gray, 'adaptive', ...
'ForegroundPolarity', 'dark', 'Sensitivity', 0.5);
% 显示结果比较
figure;
subplot(1, 3, 1);
imshow(I_gray);
title('原始图像');
subplot(1, 3, 2);
imshow(BW_global);
title('全局Otsu阈值');
subplot(1, 3, 3);
imshow(BW_adaptive);
title('自适应Otsu阈值');
end
总结
Otsu方法是一种简单而有效的自动图像阈值分割技术，特别适用于具有双峰直方图的图像。MATLAB提供了内置函数graythresh和imbinarize来简化Otsu阈值分割的实现。对于更复杂的情况，可以考虑使用自适应阈值方法或多阈值Otsu方法。
