当前位置：首页 > news >正文

MATLAB实现图像分割：Otsu阈值法

news 2025/9/2 8:27:04

MATLAB实现图像分割：Otsu阈值法

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) == 3I_gray = rgb2gray(I);
elseI_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 == 0continue;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_varmax_var = between_var;threshold = T - 1; % 转换为0-255范围的阈值endend
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, []);
endfunction best_thresholds = recursive_otsu(p, start, finish, n, current)% 递归辅助函数if n == 0% 计算当前阈值组合的类间方差var = calculate_variance(p, current);best_thresholds = {current, var};return;endbest_var = -1;best_combination = [];for t = start:finish-nnew_current = [current, t];[combination, var] = recursive_otsu(p, t+1, finish, n-1, new_current);if var > best_varbest_var = var;best_combination = combination;endendbest_thresholds = best_combination;
endfunction var = calculate_variance(p, thresholds)% 计算给定阈值组合的类间方差thresholds = sort([0, thresholds, 256]);var = 0;for i = 1:length(thresholds)-1start = thresholds(i) + 1;finish = thresholds(i+1);if start > finishcontinue;endw = sum(p(start:finish));if w == 0continue;endmu = 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/contentcsf/23056.html

实际应用示例：细胞图像分割

% 细胞图像分割示例
function cell_segmentation_example()% 读取细胞图像I = imread('cell_image.png');% 转换为灰度图if size(I, 3) == 3I_gray = rgb2gray(I);elseI_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');endtitle(['检测到的细胞数量: ', 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