最小误差阈值分割和最大类间方差阈值分割python 最小误差阈值原理步骤

基于最小误差法的胸片分割系统–matlab深度学习实战GUI项目

最小误差阈值分割法

“最小误差法所谓图像分割是指根据灰度、彩色、空间纹理、几何形状等特征把图像划分成若干个互不相交的区域,使得这些特征在同一区域内,表现出一致性或相似性,而在不同区域间表现出明显的不同。
根据图像中背景和目标像素的概率分布密度来实现的，其思想是找到一个阈值，并根据该阈值进行划分，计算出目标点误分为背景的概率和背景点误分为目标点的概率，得出总的误差划分概率。当总的误差划分概率最小时，便得到所需要的最佳阈值。

图像预处理

图像增强结果

最小误差法分割结果

区域生长算法的原理就是求图像中相似的像素的最大连通集合。把相似灰度值的像素点集合，然后依次检测该像素点周围的8个邻域像素点，把灰度值相近的点加入到之前的集合中，如果不相似则跳过该点。然后依次浏览该图像中的其他的像素点，直到图像中所有的像素点都被检查完毕为止，这样获得的区域就是按照区域生长算法得到的新区域。

参考来源: 区域生长图像处理

区域生长处理源码

主函数

clc; clear all; close all;
I = imread(fullfile(pwd, 'images/test.jpg'));
X = imadjust(I, [0.2 0.8], [0 1]);
% 阈值分割
bw = im2bw(X, graythresh(X));
[r, c] = find(bw);
rect = [min(c) min(r) max(c)-min(c) max(r)-min(r)];
Xt = imcrop(X, rect);
% 自动获取种子点
seed_point = round([size(Xt, 2)*0.15+rect(2) size(Xt, 1)*0.4+rect(1)]);
% 区域生长分割
X = im2double(im2uint8(mat2gray(X)));
X(1:rect(2), :) = 0;
X(:, 1:rect(1)) = 0;
X(rect(2)+rect(4):end, :) = 0;
X(:, rect(1)+rect(3):end) = 0;
[J, seed_point, ts] = Regiongrowing(X, seed_point);
figure(1);
subplot(1, 2, 1); imshow(I, []);
hold on;
plot(seed_point(1), seed_point(2), 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');
title('种子点自动选择');
hold off;
subplot(1, 2, 2); imshow(J, []); title('区域生长图像');
% 形态学后处理
bw = imfill(J, 'holes');
bw = imopen(bw, strel('disk', 5));
% 提取边缘
ed = bwboundaries(bw);
figure;
subplot(1, 2, 1); imshow(bw, []); title('形态学后处理图像');
subplot(1, 2, 2); imshow(I); 
hold on;
for k = 1 : length(ed)
    % 边缘
    boundary = ed{k};
    plot(boundary(:,2), boundary(:,1), 'g', 'LineWidth', 2);
end
hold off;
title('边缘标记图像');

自定义选取种子点位置
主函数：
main2.m

clc; clear all; close all;
I = imread(fullfile(pwd, 'images/test.jpg'));
X = imadjust(I, [0.2 0.8], [0 1]);
% 区域生长分割
X = im2double(im2uint8(mat2gray(X)));
[J, seed_point, ts] = Regiongrowing(X);
figure(1);
subplot(1, 2, 1); imshow(I, []);
hold on;
plot(seed_point(1), seed_point(2), 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');
title('种子点选择');
hold off;
subplot(1, 2, 2); imshow(J, []); title('区域生长分割结果');

形态学去噪处理。

提取边缘特征

主函数:

main.m

clc; clear all; close all;
warning off all;
% 读取图像
filename = fullfile(pwd, 'images/test.jpg');
Img = imread(filename);
% 灰度化
if ndims(Img) == 3
    I = rgb2gray(Img);
else
    I = Img;
end
% 直接二值化
bw_direct = im2bw(I);
figure; imshow(bw_direct); title('直接二值化分割');
% 圈选胃区域空气
c = [1524 1390 1454 1548 1652 1738 1725 1673 1524];
r = [1756 1909 2037 2055 1997 1863 1824 1787 1756];
bw_poly = roipoly(bw_direct, c, r);
figure;
imshow(I, []);
hold on;
plot(c, r, 'r-', 'LineWidth', 2);
hold off;
title('胃区域空气选择');
% 设置胃内空气为255
J = I;
J(bw_poly) = 255;
% 图像增强
J = mat2gray(J);
J = imadjust(J, [0.532 0.72], [0 1]);
J = im2uint8(mat2gray(J));
figure; imshow(J, []); title('图像增强处理');
% 直方图统计
[counts, gray_style] = imhist(J);
% 亮度级别
gray_level = length(gray_style);
% 计算各灰度概率
gray_probability  = counts ./ sum(counts);
% 统计像素均值
gray_mean = gray_style' * gray_probability;
% 初始化
gray_vector = zeros(gray_level, 1);
w = gray_probability(1);
mean_k = 0;
gray_vector(1) = realmax;
ks = gray_level-1;
for k = 1 : ks
    % 迭代计算
    w = w + gray_probability(k+1);
    mean_k = mean_k + k * gray_probability(k+1);
    % 判断是否收敛
    if (w < eps) || (w > 1-eps)
        gray_vector(k+1) = realmax;
    else
        % 计算均值
        mean_k1 = mean_k / w;
        mean_k2 = (gray_mean-mean_k) / (1-w);
        % 计算方差
        var_k1 = (((0 : k)'-mean_k1).^2)' * gray_probability(1 : k+1);
        var_k1 = var_k1 / w;
        var_k2 = (((k+1 : ks)'-mean_k2).^2)' * gray_probability(k+2 : ks+1);
        var_k2 = var_k2 / (1-w);
        % 计算目标函数
        if var_k1 > eps && var_k2 > eps
            gray_vector(k+1) = 1+w * log(var_k1)+(1-w) * log(var_k2)-2*w*log(w)-2*(1-w)*log(1-w);
        else
            gray_vector(k+1) = realmax;
        end
    end
end
% 极值统计
min_gray_index = find(gray_vector == min(gray_vector));
min_gray_index = mean(min_gray_index);
% 计算阈值
threshold_kittler = (min_gray_index-1)/ks;
% 阈值分割
bw__kittler = im2bw(J, threshold_kittler);
% 显示
figure; imshow(bw__kittler, []); title('最小误差法分割');
% 形态学后处理
bw_temp = bw__kittler;
% 反色
bw_temp = ~bw_temp;
% 填充孔洞
bw_temp = imfill(bw_temp, 'holes');
% 去噪
bw_temp = imclose(bw_temp, strel('disk', 5));
bw_temp = imclearborder(bw_temp);
% 区域标记
[L, num] = bwlabel(bw_temp);
% 区域属性
stats = regionprops(L);
Ar = cat(1, stats.Area);
% 提取目标并清理
[Ar, ind] = sort(Ar, 'descend');
bw_temp(L ~= ind(1) & L ~= ind(2)) = 0;
% 去噪
bw_temp = imclose(bw_temp, strel('disk',20));
bw_temp = imfill(bw_temp, 'holes');
figure;
subplot(1, 2, 1); imshow(bw__kittler, []); title('待处理二值图像');
subplot(1, 2, 2); imshow(bw_temp, []); title('形态学后处理图像');
% 提取肺边缘
ed = bwboundaries(bw_temp);
% 显示
figure;
subplot(2, 2, 1); imshow(I, []); title('原图像');
subplot(2, 2, 2); imshow(J, []); title('增强图像');
subplot(2, 2, 3); imshow(bw_temp, []); title('二值化图像');
subplot(2, 2, 4); imshow(I, []); hold on;
for k = 1 : length(ed)
    % 边缘
    boundary = ed{k};
    plot(boundary(:,2), boundary(:,1), 'g', 'LineWidth', 2);
end
title('肺边缘显示标记');
figure;
subplot(1, 2, 1); imshow(bw_temp, []); title('二值图像');
subplot(1, 2, 2); imshow(I, []); hold on;
for k = 1 : length(ed)
    % 边缘
    boundary = ed{k};
    plot(boundary(:,2), boundary(:,1), 'g', 'LineWidth', 2);
end
title('肺边缘显示标记');

引用说明

% MAINFORM MATLAB code for MainForm.fig
% MAINFORM, by itself, creates a new MAINFORM or raises the existing
% singleton*.
%
% H = MAINFORM returns the handle to a new MAINFORM or the handle to
% the existing singleton*.
%
% MAINFORM(‘CALLBACK’,hObject,eventData,handles,…) calls the local
% function named CALLBACK in MAINFORM.M with the given input arguments.
%
% MAINFORM(‘Property’,‘Value’,…) creates a new MAINFORM or raises the
% existing singleton*. Starting from the left, property value pairs are
% applied to the GUI before MainForm_OpeningFcn gets called. An
% unrecognized property name or invalid value makes property application
% stop. All inputs are passed to MainForm_OpeningFcn via varargin.
%
% *See GUI Options on GUIDE’s Tools menu. Choose “GUI allows only one
% instance to run (singleton)”.
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help MainForm

阈值分割测试函数：

clc; clear all; close all;
I = imread('IMG00016.jpg');
bw = im2bw(I, graythresh(I));
figure; 
subplot(1, 2, 1); imshow(I, []); title('原图像');
subplot(1, 2, 2); imshow(bw, []); title('阈值分割图像');

自定义区域生长函数:

function [J, seed_point, ts] = Regiongrowing(I, seed_point)
% 统计耗时
t1 = cputime;
% 参数检测
if nargin < 2
    % 显示并选择种子点
    figure; imshow(I,[]);  hold on;
    seed_point = ginput(1);   
    plot(seed_point(1), seed_point(2), 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');
    title('种子点选择');
    hold off;
end
% 变量初始化
seed_point = round(seed_point);
x = seed_point(2);
y = seed_point(1);
I = double(I);
rc = size(I);
J = zeros(rc(1), rc(2));
% 参数初始化
seed_pixel = I(x,y);
seed_count = 1;
pixel_free = rc(1)*rc(2); 
pixel_index = 0;
pixel_list = zeros(pixel_free, 3);
pixel_similarity_min = 0;
pixel_similarity_limit = 0.1;
% 邻域
neighbor_index = [-1 0;
        1 0;
        0 -1;
        0 1];
% 循环处理
while pixel_similarity_min < pixel_similarity_limit && seed_count < rc(1)*rc(2)   
    % 增加邻域点
    for k = 1 : size(neighbor_index, 1)
        % 计算相邻位置
        xk = x + neighbor_index(k, 1); 
        yk = y + neighbor_index(k, 2);
        % 区域生长
        if xk>=1 && yk>=1 && xk<=rc(1) && yk<=rc(2) && J(xk,yk) == 0
            % 满足条件
            pixel_index = pixel_index+1;
            pixel_list(pixel_index,:) = [xk yk I(xk,yk)]; 
            % 更新状态
            J(xk, yk) = 1;
        end
    end
    % 更新空间
    if pixel_index+10 > pixel_free 
        pixel_free = pixel_free+pixel_free;
        pixel_list(pixel_index+1:pixel_free,:) = 0;
    end
    % 统计迭代
    pixel_similarity = abs(pixel_list(1:pixel_index,3) - seed_pixel);
    [pixel_similarity_min, index] = min(pixel_similarity);
    % 更新状态
    J(x,y) = 1; 
    seed_count = seed_count+1;
    seed_pixel = (seed_pixel*seed_count + pixel_list(index,3))/(seed_count+1);
    % 存储位置
    x = pixel_list(index,1); 
    y = pixel_list(index,2);
    pixel_list(index,:) = pixel_list(pixel_index,:);
    pixel_index = pixel_index-1;
end
% 返回结果
J = mat2gray(J);
J = im2bw(J, graythresh(J));
% 统计耗时
t2 = cputime;
ts = t2 - t1;