# MATLAB--数字图像处理 Otsu算法(MATLAB原理验证)

+关注继续查看

## 概念

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

## 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
%   -------
%       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

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

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<=160
th=tp+p;
else
th=160
end                            %根据具体情况适当修正门限

Y1=zeros(r,t);
for i=1:r

for j=1:t
X1(i,j)=double(C(i,j));

end
end
for i=1:r
for j=1:t
if (X1(i,j)>=th)
Y1(i,j)=X1(i,j);
else Y1(i,j)=0;
end
end
end
%上面一段代码实现分割
figure,imshow(Y1);title('灰度门限分割图像');

|
5天前
|

27 0
|
2月前
|

22 0
|
3月前
|

32 0
|
6月前
|

85 0
|
8月前
|

758 1
|
9月前
|

44 0
|
9月前
|

60 0
|

ROS功能包:livox_camera_lidar_calibration提供了一个手动校准Livox雷达和相机之间外参的方法，已经在Mid-40，Horizon和Tele-15上进行了验证。其中包含了计算相机内参，获得标定数据，优化计算外参和雷达相机融合应用相关的代码。本方案中使用了标定板角点作为标定目标物，由于Livox雷达非重复性扫描的特点，点云的密度较大，比较易于找到雷达点云中角点的准确位置。相机雷达的标定和融合也可以得到不错的结果。 在前几篇中介绍了livox_camera_lidar_calibration功能包.以及在gazebo中搭建了标定场景.并进行外参标定,进行了简单的验
394 0
|
canal 算法
leetcode算法125.验证回文串

80 0
|

197 0