700字范文 > 数字图像处理11：阈值分割（基本全局阈值处理 Otsu 的最佳全局阈值图像平滑改善全

数字图像处理11：阈值分割（基本全局阈值处理 Otsu 的最佳全局阈值图像平滑改善全

时间：2020-02-27 03:05:22

阈值分割

从背景中提取物体的一种明显方法是，选择一个将这些模式分开的阈值 T。然后， f ( x , y ) > T f(x,y)>T f(x,y)>T的任何点 (x, y) 称为个对象点，否则该点称为背景点。

灰度阈值的成功与否直接关系到可区分直方图模式的波谷的宽度和深度。

而影响波谷特性的关键因素有以下几点：

(1)波峰间的间隔 (波峰离得越远，分离这些模式的机会越大);

(2)图像中的噪声内容 (模式随噪声的增加而展宽);

(3)物体和背景的相对尺寸;

(4)光源的均匀性;

(5)图像反射特性的均匀性。

基本全局阈值处理

当物体和背景像素的灰度分布十分明显时，可以用适用于整个图像的单个 (全局) 阈值。

实验思想

在大多数应用中，通常图像之间有较大变化，即使全局阈值是种合适的方法，也需要有能对每幅图像自动估计阈值的算法。下面的迭代算法可用于这一目的:

1.为全局阈值 T 选择一个初始估计值。

2.使用用初始的 T 分割该图像。这将产生两组像素:

G1 由灰度值大于 T 的所有像素组成，G2 由所有小于等于 T 的像素组成。

3.对 G1 和 G2 的像素分别计算平均灰度值 (均值)m1 和 m2。

4.计算一个新的阈值: T = (m1 + m2)/2

5.重复步骤 2 到步骤 4，直到连续迭代中的 T 值间的差小于个预定义的参数 dT 为止。

代码

主函数：

%% 运行 main_Global.m clc;clear;close all;%% 课本图 10.34im = imread('Fig1038(a)(noisy_fingerprint).tif'); % 原始图像 uint8% 得到图像的行数、列数、以及每像素的维数（防止出现RGB图像）[line,row,v]=size(im);im = im(:,:,1);t0 = sum(im(:))/(line*row);[im1,NK] = Global_Threshold(im,0,t0);x = 1:256;%% 将结果保存到当前目录下的result文件夹下imwrite(im, sprintf('result/%s.jpg','2_im1'));imwrite(im1, sprintf('result/%s.jpg','2_im2'));%% 显示图像figure(1); subplot(131); imshow(im); title('原图'); axis onsubplot(133); imshow(im1); title('全局阈值'); axis onsubplot(132); plot(x,NK); title('直方图'); axis on

功能函数：

function [im1,NK] = Global_Threshold(im,dT,T0)[h,w] = size(im);%% 得到直方图和累计直方图[NK,CH] = myHisteq(im);G1sum = 0;G2sum = 0;dTemp = 255;while(dTemp<dT)for i = 1:T0G1sum = G1sum + (i-1)*NK(i);endG1Ave = G1sum/CH(T0);for i = T0+1:256G2sum = G2sum + (i-1)*NK(i);endG2Ave = G2sum/(CH(256)-CH(T0));temp = round((G1Ave+G2Ave)/2);dTemp = abs(temp - T0);T0 = temp;endim1 = zeros(h,w);for i =1:hfor j = 1:wif(im(i,j)>=T0)im1(i,j) = 1;endendendim1 = logical(im1);

实验结果

用 Otsu 的最佳全局阈值

Otsu 方法是最佳的，因为它使得类间方差最大化。其基本思想是，适当的阈值化的类就其像素灰度值而言，应当是截然不同的，相反地，就其灰度值而言，给出最佳类间分离的阈值将是最佳的阈值。除了其最佳性之外，Otsu 方法还有一个重要的特性，即它完全以在一幅图像的直方图上执行计算为基础，直方图是很容易得到的一维阵列。

实验思想

Otsu 算法计算步骤如下:

1.计算输人图像的归一化直方图。使用 pi,i=0,1,2,.,L-1 表示该直方图的各个分量。

2. 对于 k=0, 1,2,.,L-1，计算累积和 P1(k) 和 P2(k), 计算公式为：

P 1 ( k ) = ∑ i = 0 k p i = 1 − P 2 ( k ) P1(k) = \sum_{i=0}^{k} pi = 1 - P2(k) P1(k)=i=0∑kpi=1−P2(k)

3.对于 k=0,1,2…,L-1, 计算累积均值 m1(k) 和 m2(k)。计算公式为：

m 1 ( k ) = 1 P 1 ( k ) ∗ ∑ i = 0 k ( i ∗ p i ) m1(k) = \frac{1}{P1(k)} * \sum_{i=0}^{k} (i*pi) m1(k)=P1(k)1∗i=0∑k(i∗pi) m 2 ( k ) = 1 P 2 ( k ) ∗ ∑ i = k + 1 L − 1 ( i ∗ p i ) m2(k) = \frac{1}{P2(k)} * \sum_{i=k+1}^{L-1} (i*pi) m2(k)=P2(k)1∗i=k+1∑L−1(i∗pi)

4.计算全局灰度均值 mg。计算公式为： m g = ∑ i = 0 L − 1 ( i ∗ p i ) mg = \sum_{i=0}^{L-1} (i*pi) mg=i=0∑L−1(i∗pi)

5.对于 k=0,1,2…,L-1, 计算类间方差 b 2 ( k ) b^2(k) b2(k)。计算公式：

b 2 ( k ) = P 1 ( k ) ∗ ( m 1 ( k ) − m g ) 2 + P 2 ( k ) ∗ ( m 2 ( k ) − m g ) 2 b^2(k) = P1(k)*(m1(k)-mg)^2 + P2(k)*(m2(k)-mg)^2 b2(k)=P1(k)∗(m1(k)−mg)2+P2(k)∗(m2(k)−mg)2

6.得到 Otsu 阈值 k，即使得 b2(k) 最大的 k 值。如果最大值不唯一，用相应检测到的各个值 k 的平均得到 k*。

代码

主函数：

%% 运行 main_Otsu.m clc;clear;close all;%% 课本图 10.39im = imread('Fig1039(a)(polymersomes).tif'); % 原始图像 uint8% 得到图像的行数、列数、以及每像素的维数（防止出现RGB图像）[line,row,v]=size(im);im = im(:,:,1);[im1,NK] = Global_Threshold(im,0,sum(im(:))/(line*row));[im2,NK] = my_Otsu(im);x = 1:256;%% 将结果保存到当前目录下的result文件夹下imwrite(im, sprintf('result/%s.jpg','3_im1'));imwrite(im1, sprintf('result/%s.jpg','3_im2'));imwrite(im2, sprintf('result/%s.jpg','3_im3'));%% 显示图像figure(1); subplot(221); imshow(im); title('原图'); axis onsubplot(222); plot(x,NK); title('直方图'); axis onsubplot(223); imshow(im1); title('全局阈值'); axis onsubplot(224); imshow(im2); title('otsu阈值'); axis on

功能函数：

function [im1,NK] = my_Otsu(im)[h,w] = size(im);%% 像素值从0到255，下标从1到256%% 得到直方图和累计直方图[NK,CH] = myHisteq(im);%归一化nK = NK/(h*w);cH = CH/(h*w);%% 计算P1k P2kP1 = cH;P2 = 1-cH;%% 累计均值的计算，仅有255个级别，比像素点取值总数小1m1 = zeros(255,1);for k = 1:255for j = 1:km1(k) = m1(k) + j*nK(j);endm1(k) = m1(k)/cH(k);endm2 = zeros(255,1);for k = 1:255for j = k+1:256m2(k) = m2(k) + j*nK(j);endm2(k) = m2(k)/(cH(256)-cH(k));end%% 全局灰度均值mG = 0;for j = 1:256mG = mG + j*nK(j);endmG = mG/cH(256);%% 计算类间方差phi2 = zeros(255,1);for k = 1: 255phi2(k) = P1(k)*(m1(k)-mG)^2 + P2(k)*(m2(k)-mG)^2;end%% 最优phi2值m = max(phi2);%% 防止存在相等的最优值num = 0;k = zeros(10,0);for i = 1:255if(phi2(i) == m) num = num+1;k(num) = i;endendBestk = sum(k)/num;%% 得到最优二值化图像im1 = zeros(h,w);for i =1:hfor j = 1:wif(im(i,j)>=Bestk)im1(i,j) = 1;endendendim1 = logical(im1);

实验结果

图像平滑改善全局阈值处理

这道题未在本次实验要求范围内，但给了这张图片，于是进行了下面的实验。

实验思想

噪声会将简单的阈值处理问题变得不可解决，当噪声不能在源头减少，并且阈值处理又是所以选择的分割方法时，那么通常能增强性能的一种技术是，在阈值处理之前平滑图像。

部分核心代码

主函数：

%% 运行 main_SmoothOtsu.m clc;clear;close all;%% 课本图 10.39im = imread('Fig1040(a)(large_septagon_gaussian_noise_mean_0_std_50_added).tif'); % 原始图像 uint8% 得到图像的行数、列数、以及每像素的维数（防止出现RGB图像）[line,row,v]=size(im);im = im(:,:,1);[im1,NK1] = my_Otsu(im);im2 = myAverage(im);[im3,NK2] = my_Otsu(im2);x = 1:256;%% 将结果保存到当前目录下的result文件夹下imwrite(im, sprintf('result/%s.jpg','4_im1'));imwrite(im1, sprintf('result/%s.jpg','4_im2'));imwrite(im2, sprintf('result/%s.jpg','4_im3'));imwrite(im3, sprintf('result/%s.jpg','4_im4'));%% 显示图像figure(1); subplot(231); imshow(im); title('原图'); axis onsubplot(232); plot(x,NK1); title('直方图'); axis onsubplot(233); imshow(im1); title('全局阈值'); axis onsubplot(234); imshow(im2); title('原图'); axis onsubplot(235); plot(x,NK2); title('直方图'); axis onsubplot(236); imshow(im3); title('全局阈值'); axis on

功能函数前面都已经给出过

function [img_2] = myAverage(img_1)size_1 = size(img_1);h = size_1(1);w = size_1(2);img_2 = zeros(h, w);%%边缘延拓四行四列a = img_1(2:-1:1,:);b = img_1(h:-1:h-1,:);img_1 = [a;img_1;b];c = img_1(:,2:-1:1);d = img_1(:,w:-1:w-1);img_1 = double([c,img_1,d]);%5X5均值模板L = 1/25*[1 1,1,1 1;1,1,1 1 1;1,1,1 1 1;1,1,1 1 1;1,1,1 1 1];for i= 1:hfor j = 1:wim = [img_1(i,j) img_1(i,j+1) img_1(i,j+2) img_1(i,j+3) img_1(i,j+4);...img_1(i+1,j) img_1(i+1,j+1) img_1(i+1,j+2) img_1(i+1,j+3) img_1(i+1,j+4);...img_1(i+2,j) img_1(i+2,j+1) img_1(i+2,j+2) img_1(i+2,j+3) img_1(i+2,j+4);...img_1(i+3,j) img_1(i+3,j+1) img_1(i+3,j+2) img_1(i+3,j+3) img_1(i+3,j+4);...img_1(i+4,j) img_1(i+4,j+1) img_1(i+4,j+2) img_1(i+4,j+3) img_1(i+4,j+4);...];img_2(i,j) = round(sum(sum(L.*im)));endendimg_2 =uint8(img_2);end

实验结果

图像分块的可变阈值

由于噪声或者光照不均匀等等一系列问题，导致无法直接进行阈值分割，图像平滑和边缘信息有益于阈值处理，但是在更经常的情况下，上述的两种做法效果不明显，只能使用可变阈值来进行解决。

实验思想

可变阈值处理最简单的方法之一是，把一幅图像分成不重叠的矩形。这种方法用于补偿光照和或反射的不均匀性。选择的矩形要足够小，以便每个矩形的光照都近似是均匀的。

代码

主函数：

%% 运行 main_divideOtsu.m clc;clear;close all;%% 课本图 10.39im = imread('Fig1046(a)(septagon_noisy_shaded).tif'); % 原始图像 uint8% 得到图像的行数、列数、以及每像素的维数（防止出现RGB图像）[line,row,v]=size(im);im = im(:,:,1);[im1,NK] = Global_Threshold(im,0,sum(im(:))/(line*row));[im2,NK] = my_Otsu(im);[im3,im4] = my_divideOtsu(im);x = 1:256;%% 将结果保存到当前目录下的result文件夹下imwrite(im, sprintf('result/%s.jpg','5_im1'));imwrite(im1, sprintf('result/%s.jpg','5_im2'));imwrite(im2, sprintf('result/%s.jpg','5_im3'));imwrite(im3, sprintf('result/%s.jpg','5_im4'));imwrite(im4, sprintf('result/%s.jpg','5_im5'));%% 显示图像figure(1); subplot(231); imshow(im); title('原图'); axis onsubplot(232); plot(x,NK); title('直方图'); axis onsubplot(233); imshow(im1); title('全局阈值'); axis onsubplot(234); imshow(im2); title('otsu阈值'); axis onsubplot(235); imshow(im3); title('划分区域'); axis onsubplot(236); imshow(im4); title('可变阈值'); axis on

功能函数：

function [imd,imv] = my_divideOtsu(im)[h,w] = size(im);h1 = floor(h/2);w1 = floor(w/3);im1 = im(1:h1,1:w1);im2 = im(1:h1,w1+1:2*w1);im3 = im(1:h1,2*w1+1:w);im4 = im(h1+1:h,1:w1);im5 = im(h1+1:h,w1+1:2*w1);im6 = im(h1+1:h,2*w1+1:w);imd = im;imd(h1,:) = 255;imd(:,w1) = 255;imd(:,2*w1) = 255;[im1,n] = my_Otsu(im1);[im2,n] = my_Otsu(im2);[im3,n] = my_Otsu(im3);[im4,n] = my_Otsu(im4);[im5,n] = my_Otsu(im5);[im6,n] = my_Otsu(im6);imv = [im1,im2,im3;im4,im5,im6];