MATLAB

#MATLAB| 来源: 网络整理| 查看: 265

「这是我参与11月更文挑战的第17天，活动详情查看：2021最后一次更文挑战」

前言

Hello！小伙伴！

非常感谢您阅读海轰的文章，倘若文中有错误的地方，欢迎您指出～

自我介绍 ଘ(੭ˊᵕˋ)੭

昵称：海轰

标签：程序猿｜C++选手｜学生

简介：因C语言结识编程，随后转入计算机专业，有幸拿过一些国奖、省奖...已保研。目前正在学习C++/Linux/Python

学习经验：扎实基础 + 多做笔记 + 多敲代码 + 多思考 + 学好英语！

概念

OTSU算法是由日本学者OTSU于1979年提出的一种对图像进行二值化的高效算法。（大津算法）

Otsu原理

对于图像 t (x,y)，前景(即目标)和背景的分割阈值记作 T，属于前景的像素点数占整幅图像的比例记为 ω0，平均灰度为 μ0；背景像素点数占整幅图像的比例为 ω1，平均灰度为 μ1；整幅图像的平均灰度记为μ，类间方差记为g。假设图像大小为M×N，图像中像素的灰度值小于阈值 T 的像素个数为 N0，像素灰度大于阈值T的像素个数为 N1 在这里插入图片描述注：（7）式是将（5）式代入（6）式得到的，我们的重点放在（7）式上。

Otsu用处

利用Otsu算法，我们可以得到一个阈值，利用该阈值对图像进行二值化等操作。相比于单阈值的固定阈值，otsu算法效果更好。

MATLAB中实现Otsu算法的是 garythresh()函数，一般都与im2bw()配套使用

例：

t=rgb2gray(imread('a1.jpg')); x=graythresh(t);%获取otsu算得的阈值 t=im2bw(t,x); graythresh()源码--MATLAB function [level em] = graythresh(I) %GRAYTHRESH Global image threshold using Otsu's method. % LEVEL = GRAYTHRESH(I) computes a global threshold (LEVEL) that can be % used to convert an intensity image to a binary image with IM2BW. LEVEL % is a normalized intensity value that lies in the range [0, 1]. % GRAYTHRESH uses Otsu's method, which chooses the threshold to minimize % the intraclass variance of the thresholded black and white pixels. % % [LEVEL EM] = GRAYTHRESH(I) returns effectiveness metric, EM, as the % second output argument. It indicates the effectiveness of thresholding % of the input image and it is in the range [0, 1]. The lower bound is % attainable only by images having a single gray level, and the upper % bound is attainable only by two-valued images. % % Class Support % ------------- % The input image I can be uint8, uint16, int16, single, or double, and it % must be nonsparse. LEVEL and EM are double scalars. % % Example % ------- % I = imread('coins.png'); % level = graythresh(I); % BW = im2bw(I,level); % figure, imshow(BW) % narginchk(1,1); validateattributes(I,{'uint8','uint16','double','single','int16'},{'nonsparse'}, ... mfilename,'I',1); if ~isempty(I) % Convert all N-D arrays into a single column. Convert to uint8 for % fastest histogram computation. I = im2uint8(I(:)); num_bins = 256; counts = imhist(I,num_bins); % Variables names are chosen to be similar to the formulas in % the Otsu paper. p = counts / sum(counts); omega = cumsum(p); mu = cumsum(p .* (1:num_bins)'); mu_t = mu(end); sigma_b_squared = (mu_t * omega - mu).^2 ./ (omega .* (1 - omega)); % Find the location of the maximum value of sigma_b_squared. % The maximum may extend over several bins, so average together the % locations. If maxval is NaN, meaning that sigma_b_squared is all NaN, % then return 0. maxval = max(sigma_b_squared); isfinite_maxval = isfinite(maxval); if isfinite_maxval idx = mean(find(sigma_b_squared == maxval)); % Normalize the threshold to the range [0, 1]. level = (idx - 1) / (num_bins - 1); else level = 0.0; end else level = 0.0; isfinite_maxval = false; end % compute the effectiveness metric if nargout > 1 if isfinite_maxval em = maxval/(sum(p.*((1:num_bins).^2)') - mu_t^2); else em = 0; end end Ostu算法（个人实现 MATLAB版） t=rgb2gray(imread('a1.jpg')); [m,n]=size(t); %counts为图片总像素个数值 counts=m*n; %count是一个256行一列的矩阵记录了每个灰度级上像素点的个数 count=imhist(t); %每个灰度级上像素点占总像素点数量的比例（概率） p=count/counts; %cumsum 计算累加概率 w0=cumsum(p); %u 计算的是平均灰度比如灰度级为0~125的平均灰度值 %(1:256)' 矩阵转置 1行256列转换为256行1列 u=cumsum(p.*(1:256)'); %u_end是全局平均灰度 u_end=u(end); %d 求方差 d是256行一列的矩阵 d=(w0*u_end-u).^2./(w0.*(1-w0)); %在d中寻找为最大方差的灰度值这里y是最大方差的灰度值 [x,y]=max(d); %为了和im2bw配合使用进行归一化 %x就是所得阈值 x=(y-1)/255

算法结果验证

t=rgb2gray(imread('a1.jpg')); [m,n]=size(t); counts=m*n; count=imhist(t); p=count/counts; w1=cumsum(p); u=cumsum(p.*(1:256)'); u_end=u(end);%u_end是全局平均灰度 d=(w1*u_end-u).^2./(w1.*(1-w1)); [x,y]=max(d); x=(y-1)/255 %x就是所得阈值 xx=graythresh(t)%graythresh计算出的阈值

结果　

%自己编写的代码计算出的阈值 x = 0.7569 %利用garythresh计算出的阈值 xx = 0.7569

结论该代码与gragthresh（）计算出的阈值大体上完全一致！

疑点解惑

在我们自己编写的代码中，计算方差用的是下面的公式，与原理中计算方差的公式（ g=ω0ω1(μ0-μ1)^ ）不一样难道是原理错了？哈哈不是，其实都是一个公式，我们在代码中的公式只是原理公式的变形罢了，这是为了减少变量，提高运算速度。

%d 求方差 d是256行一列的矩阵 d=(w1*u_end-u).^2./(w1.*(1-w1));

公式变形

在这里插入图片描述注意：上图中的u就是验证代码中的u_end,代表的是图像全局的平均灰度而u是代表灰度值 T前面的平均灰度。这里要注意图中的w0u0就是验证代码中的u，对每个灰度级数乘以对应的概率并累加，这是因为u0的计算公式是每个灰度级乘以对应的数量，再除以前面一共的像素数量，这里前面一共的像素数量可以用总量乘以前面占的概率w0。哈哈，对，就是这里，就可以把w0u0中的w0约掉，所以w0*u0==u（这里的u是验证代码中的u，代表灰度值乘以概率的累加值）

其实只需要知道g=w0w1(u1-u0)^2就行了，验证代码只是为了减少变量，做了一下变形，哎，数学还是要学好啊下面在给个版本的验证算法，原理差不多，只是个别地方做了精简（其实现在看来，这个算法反而搞复杂了）

C=imread('h.jpg'); C=rgb2gray(C); %读取图像 figure,imshow(C);title('原始灰度图像');%绘原图 count=imhist(C); %直方图统计 subplot(1,3,1); imhist(C); %绘制直方图 [r,t]=size(C); %图像矩阵大小 N=r*t; %图像像素个数 L=256; %制定凸显灰度级为256 count=count/N; %各级灰度出现概率 for i=2:L if count(i)~=0 ; st=i-1; break end end %以上循环语句实现寻找出现概率不为0的最小灰度值 for i=L:-1:1 if count(i)~=0; nd=i-1; break end end %实现找出出现概率不为0的最大灰度值 f=count(st+1:nd+1); p=st;q=nd-st; %p和q分别是灰度的起始和结束值 u=0; for i=1:q u=u+f(i)*(p+i-1); ua(i)=u; end %计算图像的平均灰度值 for i=1:q w(i)=sum(f(1:i)); end %计算出选择不同k的时候，A区域的概率 d=(u*w-ua).^2./(w.*(1-w)); %求出不同k值时类间方差 [y,tp]=max(d); %求出最大方差对应的灰度值 th=tp+p; if th=th) Y1(i,j)=X1(i,j); else Y1(i,j)=0; end end end %上面一段代码实现分割 figure,imshow(Y1);title('灰度门限分割图像');

【本文地址】

公司简介

联系我们