彩色图像转灰度图_灰度图转二值图_双峰法自动阈值二值化[Matlab]

彩色图像转灰度图_灰度图转二值图_双峰法自动阈值二值化

彩色转灰度图

matlab提供了彩色图像转灰度图的函数——rgb2gray(colorImg)。该方法内部是通过Gray = 0.29900 * R + 0.58700 * G + 0.11400 * B该公式计算的。

这里采用了三种自定义函数实现彩色图到灰度图的转换

  • 采用三通道的均值作为灰度值:CvtGrayByMean

    % 彩色转灰度:均值
    function gray = CvtGrayByMean(image)
    image=double(image); %输入图转浮点
    gray= (image (:,:,1)+ image (:,:,2)+ image(:,:,3))/3; 
    gray=round(gray); %四舍五入
    gray=uint8(gray); %转整数
    end
  • 采用单通道,如绿色通道值,作为灰度值:CvtGrayByGreen

    % 彩色转灰度:green通道值
    function gray = CvtGrayByGreen(image)
    gray= image (:,:,2); 
    end
  • 采用灰度转换公式,Gray = 0.29900 * R + 0.58700 * G + 0.11400 * B:CvtGrayBy305911

    % 彩色转灰度:305911公式
    function gray = CvtGrayBy305911(image)
    image=double(image); %输入图转浮点
    gray=0.299* image (:,:,1)+0.587* image (:,:,2)+0.114* image(:,:,3); 
    gray=round(gray); %四舍五入
    gray=uint8(gray); %转整数
    end

完整代码如下

% input original image
input = imread('lena.jpg');
figure(1)
set(gcf,'Name','RGB2GRAY');

% show origianl image
subplot(2,3,1);
imshow(input);
title('input');

% show gray image derived from matlab 
gray = rgb2gray(input);
subplot(2,3,2);
imshow(gray);
title('gray by bgr2gray');

% 2.1 method 1: gray by mean 
gray1 = CvtGrayByMean(input);
subplot(2,3,4);
imshow(gray1);
title('gray by mean');

% 2.2 method 2: gray by Green channel
gray2 = CvtGrayByGreen(input);
subplot(2,3,5);
imshow(gray2);
title('gray by Green');

% 2.3 method 3: gray by 305911 formula
gray3 = CvtGrayBy305911(input);
subplot(2,3,6);
imshow(gray3);
title('gray by 305911');

结果:
rgb2gray


灰度图转二值图

matlab提供了阈值化函数,im2bw或者imbinarize,输入灰度图,然后根据阈值进行阈值话操作,通常会采用graythresh(gray)方法自动获取阈值,该方法内部采用的OTSU算法(大律算法、最大类间方差法)。

threshold = graythresh(gray);
bin = imbinarize(gray,threshold);
imshow(bin);

本文实现了几种常见的阈值化方法:

  • 指定阈值为127:CvtBinBy127

    % 灰度图转二值图:127为阈值
    function [binary] = CvtBinBy127(gray)
    binary = gray;
    low = gray<127;
    binary(low) = 0;
    high = gray>=127;
    binary(high) = 255;
    end
  • 灰度均值作为阈值:CvtBinByMean

    % 灰度图转二值图:均值为阈值
    function [binary] = CvtBinByMean(gray)
    m = mean(gray);
    binary = gray;
    low = gray<m;
    binary(low) = 0;
    high = gray>=m;
    binary(high) = 255;
    end
    
  • 双峰法,:CvtBinByHist(未实现自动阈值,是根据观测灰度直方图求波谷做为阈值,自动阈值求法见第三部分)

    % 灰度图转二值图:直方图双峰法阈值,难点为统计双峰之间的波谷,未实现,采用人工设定阈值
    function [binary] = CvtBinByHist(gray,threshold)
    binary = gray;
    low = gray<threshold;
    binary(low) = 0;
    high = gray>=threshold;
    binary(high) = 255;
    end

完整代码如下(灰度图来自前文):

figure(2)
set(gcf,'Name','GRAY2BIN_IMAGE');

% show origianl image
subplot(2,3,1);
imshow(gray);
title('original');

% show hist
subplot(2,3,2);
hist = imhist(gray, 256);
plot(hist);
% 设置x轴范围
xlim([0,260]);
title('Hist');

% show binary image derived from matlab, graythresh(OTSU) + imbinarize
subplot(2,3,3);
threshold = graythresh(gray);
bin = imbinarize(gray,threshold);
imshow(bin);
title('BinImgByOSTU');

% 3.1 method 1: bin by threshold of 127, 系统函数为 imbinarize(gray,0.5)
subplot(2,3,4);
bin1 = CvtBinBy127(gray);
imshow(bin1);
title('BinImgBy127');

% 3.2 method 1: bin by threshold of Mean
subplot(2,3,5);
bin2 = CvtBinByMean(gray);
imshow(bin2);
title('BinImgByMean');

% 3.3 method 3: bin by hist,双峰法【非自动双峰阈值】
% 注意:此处,没有实现双峰自动求取阈值,而是根据直方图绘制后人工得出观测输入图的两个峰值,再求取波谷阈值
subplot(2,3,6);
% 根据直方图手动获取阈值,根据图中双峰的位置,可知波谷应该在[130,155]之间产生
[var,index] = min(hist(130:155));
thresh = 130+index;
% 上式得出thresh = 130+index,var为统计的波谷灰度值像素数目
bin3 = CvtBinByHist(gray,thresh);
imshow(bin3);
title('BinImgByHist');

结果
gray2bin


双峰法自动求阈值

前文中的双峰法,难点在于求直方图中的两个最大波峰,一直没有头绪,只好采用观测的方法,人工观察出两个最高的波峰,然后求波峰间的最小值作为阈值,这种办法适应性很差,无法自动取阈值。

在网络上没有找到相应的解决办法,之后在吃饭中想到,这个问题不就是求局部最大(小)值点么,可以用导数,直方图hist实际是一个一维的数组,有很多梯度算子实际就是导数的思想,那么可以利用这个思想来求最大的两个波峰,然后取波峰间的波谷就可以了。即找到两个导数为0,且值为最大和第二大的点,那他们就是两个峰值,再找到它俩之间,导数为0的取值最小的点就是波谷,也就是我们要求的阈值。这里为了增加算法可靠度,采取以下两个策略:

  • 1、导数不一定是等于0才认为是局部最大值小值点,因为这是离散情况下的近似求法,所以采用一个阈值,小于这个值得导数都认为是局部最值点。
  • 【补充】在实验过程中,发现取导数阈值的做法不不科学。
    • 原因:若设置阈值,则会认为局部极大值附近的点都会满足条件,会被收入到备选点集合中,这样在寻找最大的两个峰值点的时候,找到的实际上是一个最大波峰上的最大的两个点,这显然是不对的。
    • 解决方案:不设置阈值,直接在导数数组中,找这样的点——该点前一个点导数>0,后一个点导数<0,那么就认为该点是极大值点,这样显然更可靠!
  • 2、求出两个波峰之后,直接在两个波峰之间求最小值作为最终得阈值,不再用求导得方法求最小值。
  • 【补充】在用hist求梯度前,一定要先平滑滤波,因为求导对噪声的干扰很大,必须先做出平滑曲线,平滑操作也决定了第一条中的阈值可以不用做了。显然平滑方案,比不平滑但同时取阈值,要科学。

步骤如下:

  • 读取图片,转换成灰度图
  • 求灰度图的直方图hist
  • 将hist进行平滑滤波
  • 将平滑后的hist进行梯度操作
  • 根据梯度结果求出所有局部极大值点。
  • 再从这些点中找到两个最大值作为两个波峰
  • 求两个波峰间的最小值,作为最终的二值化阈值。

存在的问题:当hist曲线的两端为最大值点的时候,上述做法会检测不到这个最值点(它不是最值点,但双峰法明显应该将它作为一个峰值)
解决方法:可以取直方图的左右两端各n各点(如n=3,可能n=1更好理解,但是因为有高斯平滑,这里可能会漏掉一些波峰,所以n最好大于1),然后在两端各求出一个极大值点,放到极大值点集合中,用于下一步的求最大的两个波峰。(本文未实现该部分代码)

结果图:

  • 图1:原灰度图,其灰度直方图,高斯blur之后的直方图曲线,根据直方图曲线求取的导数图

    1

  • 图2:根据导数图,求出的可能的极大值点(中间结果,用于观测代码合理性,在这里发现了相邻极大值点问题,下文简单采用只保留相邻点中的一个极大值点做法)

    2

  • 图3:上述删除相邻极值点之后的结果,该图中的所有点就是待选的极大值点。

    3

  • 图4:根据图3的结果求两个极大值点,然后求出它俩之间的最小值,作为波谷(阈值)

    4

  • 图5:根据上述阈值,得到的二值化图像

    5

代码如下:

% 双峰法求二值化得阈值【包含测试】

% 1.input image
input = imread('lena.jpg');
gray = rgb2gray(input);

% 2.convert to gray image
figure(1);
subplot(2,2,1);
imshow(gray);
title('gray');

% 3. calc hist of gray image 
subplot(2,2,2);
hist = imhist(gray, 256);
plot(hist);
grid on;
title('hist');

% 4. blur hist
% 梯度算子对噪声敏感,所以先对hist曲线进行高斯滤波器平滑,n为模板长度
n=5;  
%gaussian为1×n的模板,sigma = 1
gau = fspecial('gaussian', [1 n], 1);  
blurHist = conv(hist,gau);  
% 这里的模板操作后的长度范围没搞懂,没搞懂它的对应关系,反正模板操作会增加长度
% 不知道它的卷积的过程,不理解为何是舍弃最后几个元素
blurHist=blurHist(1:(max(size(blurHist))-max(size(gau))+1)); 
subplot(2,2,4);
plot(blurHist);
grid on;
title('gaus\_blur\_hist'); % 这里str自动支持latex,所以要转义

% 求blurHist一维数组的双峰,应该是可以用sobel 1*n的算子卷积计算的,没搞懂conv的机制,自己写算了
kernel = [-1 -4 -5 0 5 4 1]; % 梯度计算模板:1*7 sobel
len = max(size(blurHist));   % 数组元素个数
grad = zeros(1,len);        % 梯度计算结果

% 由于模板是1*7的,最开始3个元素和最后3个元素无法计算,全部取0
% 计算梯度
for i=4:len-3
    grad(i) = blurHist(i-3)*kernel(1) + blurHist(i-2)*kernel(2) + blurHist(i-3)*kernel(3) + blurHist(i+1)*kernel(5) + blurHist(i+2)*kernel(6) +blurHist(i+3)*kernel(7);
end
% 梯度图
subplot(2,2,3);
plot(grad);
grid on;
title('grad');

% 求波峰:遍历grad数组,寻找左右异号,且左大0,右小0的点
% result保存备选的局部极大值点
result = [];
for i = 2:(length(grad)-1)
   if(grad(i-1)*grad(i+1)<0 && grad(i-1)>0 && grad(i+1)<0 )
       result = [result,i];
   end
end    

% 重绘blurHist和grad图,标出初步的峰值点看看结果
% 发现图中有临近的两个点判定为峰值的情况,分析原因:x轴的取值间隔过大导致
% 解决办法,若result中有两个连续的点,则取值大的作为局部极大值,另一个舍弃
% TEST代码—观察效果
figure(2);
subplot(2,1,1);
plot(blurHist);
grid on;
title('gaus\_blur\_hist');
set(gca,'Xtick',result);

subplot(2,1,2);
plot(grad);
grid on;
title('grad');
set(gca,'Xtick',result);

% 通过上述观察,修正结果,删除相邻元素
% 在循环中直接删除会出问题,所以,先求出相邻元素的索引
indexs = [];
for i = 1:(length(result)-1)
    if(result(i)==result(i+1)-1)
        indexs = [indexs,i];
    end
end
result(indexs) = []; % 删除元素

% 再来看看结果:可以接受,得到了近似的极大值点,但存在误差,并不是完全准确
figure(3);
subplot(2,1,1);
plot(blurHist);
grid on;
title('gaus\_blur\_hist');
set(gca,'Xtick',result);

subplot(2,1,2);
plot(grad);
grid on;
title('grad');
set(gca,'Xtick',result);

% 寻找最大的两个波峰,索引为v1和v2
[~,maxp] = max(blurHist(result));
v1 = result(maxp);
% 删除第一个最大值,再重复操作,找到第二个最大值
result(maxp) = [];
[~,maxp] = max(blurHist(result));
v2 = result(maxp);

% 画图看结果
figure(4);
subplot(2,1,1);
plot(blurHist);
grid on;
title('gaus\_blur\_hist');
if v2>v1
    xtick = [v1 v2];
else 
    xtick = [v2 v1];
end
set(gca,'Xtick',xtick);

% 找v1 v2间的最小值,作为波谷,即最终的阈值(注意,阈值是直方图的x轴索引,而不是value)
threshold = xtick(1);
for i = xtick(1):xtick(2)
    if hist(i)<hist(threshold)
        threshold = i;
    end
end

% 波谷图
subplot(2,1,2);
plot(hist);
grid on;
title('min\_in\_hist');
set(gca,'Xtick',threshold);

% 又,hist的x轴是[2,256],故threshold再减一即可
threshold = threshold - 1;

% 二值化
figure(5);
bin = imbinarize(gray,threshold/255);
imshow(bin);
title('BinImgByHist');

核心代码【不包含中间过程的图像绘制】

% 双峰法求二值化得阈值

% 1.input image
input = imread('lena.jpg');
gray = rgb2gray(input);

% 2.convert to gray image
figure(1);
subplot(2,2,1);
imshow(gray);
title('gray');

% 3. calc hist of gray image 
subplot(2,2,2);
hist = imhist(gray, 256);
plot(hist);
grid on;
title('hist');

% 4. blur hist
% 梯度算子对噪声敏感,所以先对hist曲线进行高斯滤波器平滑,n为模板长度
n=5;  
%gaussian为1×n的模板,sigma = 1
gau = fspecial('gaussian', [1 n], 1);  
blurHist = conv(hist,gau);  
% 这里的模板操作后的长度范围没搞懂,没搞懂它的对应关系,反正模板操作会增加长度
% 不知道它的卷积的过程,不理解为何是舍弃最后几个元素
blurHist=blurHist(1:(max(size(blurHist))-max(size(gau))+1)); 
subplot(2,2,4);
plot(blurHist);
grid on;
title('gaus\_blur\_hist'); % 这里str自动支持latex,所以要转义

% 求blurHist一维数组的双峰,应该是可以用sobel 1*n的算子卷积计算的,没搞懂conv的机制,自己写算了
kernel = [-1 -4 -5 0 5 4 1]; % 梯度计算模板:1*7 sobel
len = max(size(blurHist));   % 数组元素个数
grad = zeros(1,len);        % 梯度计算结果

% 由于模板是1*7的,最开始3个元素和最后3个元素无法计算,全部取0
% 计算梯度
for i=4:len-3
    grad(i) = blurHist(i-3)*kernel(1) + blurHist(i-2)*kernel(2) + blurHist(i-3)*kernel(3) + blurHist(i+1)*kernel(5) + blurHist(i+2)*kernel(6) +blurHist(i+3)*kernel(7);
end

% 求波峰:遍历grad数组,寻找左右异号,且左大0,右小0的点
% result保存备选的局部极大值点
result = [];
for i = 2:(length(grad)-1)
   if(grad(i-1)*grad(i+1)<0 && grad(i-1)>0 && grad(i+1)<0 )
       result = [result,i];
   end
end    

% 修正结果,删除相邻元素
% 在循环中直接删除会出问题,所以,先求出相邻元素的索引
indexs = [];
for i = 1:(length(result)-1)
    if(result(i)==result(i+1)-1)
        indexs = [indexs,i];
    end
end
result(indexs) = []; % 删除元素

% 寻找最大的两个波峰,索引为v1和v2
[~,maxp] = max(blurHist(result));
v1 = result(maxp);
% 删除第一个最大值,再重复操作,找到第二个最大值
result(maxp) = [];
[~,maxp] = max(blurHist(result));
v2 = result(maxp);

% 找v1 v2间的最小值,作为波谷,即最终的阈值(注意,阈值是直方图的x轴索引,而不是value)
threshold = xtick(1);
for i = xtick(1):xtick(2)
    if hist(i)<hist(threshold)
        threshold = i;
    end
end
% 又,hist的x轴是[2,256],故threshold再减一即可
threshold = threshold - 1;

% 二值化
figure(5);
bin = imbinarize(gray,threshold/255);
imshow(bin);
title('BinImgByHist');
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值