转自http://hi.baidu.com/yibobin/blog/item/3a84c8f7047dee2f720eec81.html
一维直方图熵阈值算法
% 1D entropy thresholding method
% Pun提出，Kapur对其阈值和熵进行改进
% 两类：object 和background
% P1 = sum(pi) i：1~T
% P2 = sum(pi) i:T+1~255
% HO = ln(P1) + H1/P1;
% HB = ln(P2) + H2/P2;
% H1 = -sum(pi*ln(pi)); i:1~T
% H2 = -sum(pi*ln(pi)); i:T+1~255
% H = HO + HB;
% T_best = argmax(H);
clc;
clear;
%I = imread(‘E:/test/chinalake.bmp’,’bmp’);
I = imread(‘E:/test/lena.png’,’png’);
I = double(I);
I = Medianfilter(I); % median filter
[height,width] = size(I);
Size = height * width; % the size of the image
h_T = sum(sum(I)); % the total gray value of the image
G_min = min(min(I)); % the min gray value of the image
G_max = max(max(I)); % the max gray value of the iamge
I_seg = zeros(height,width); % the array to store the segmented image
I_hist = zeros(1,256); % the array to store the hist of the image
thresh = 0; % the threshold
num1 = 0;
num2 = 0; % count the num of the pixel from the diffrient class
P1 = 0;
P2 = 0; % the probability of the different class
h_T1 = 0;
h_T2 = 0; % the total gray value of different class
max = 0;
H1 = 0;
H2 = 0; % the middle var
H_object = 0;
H_background = 0;
H_total = 0; % the total entropy
T_best = 0; % the best thresh
%%%%% 计算直方图 %%%%%%
for i=1:height % calculate the hist of the image
for j=1:width
I_hist(I(i,j)+1) = I_hist(I(i,j)+1) + 1;
end
end
for thresh=G_min:G_max % find the best threshold
H1 = 0;
h_T1 = 0;
H2 = 0;
for h=1:height
for w=1:width
if I(h,w) <= thresh
num1 = num1 + 1;
h_T1 = h_T1 + I(h,w);
end
end
end
num2 = Size – num1;
h_T2 = h_T – h_T1;
P1 = num1/Size;
P2 = num2/Size;
for i=1:thresh
px = I_hist(i+1)/Size;
H1 = H1 + (-px*ln(px));
end
for i=thresh+1:G_max
px = I_hist(i+1)/Size;
H2 = H2 + (-px*ln(px));
end
H_object = ln(P1) + H1/P1;
H_background = ln(P2) + H2/P2;
H_total = H_object + H_background;
if H_total > max
max = H_total;
T_best = thresh;
end
end
%%%%%%% Seg the image %%%%%%%%%
for i=1:height
for j=1:width
if I(i,j) > T_best
I_seg(i,j) = 255;
end
end
end
T_best
figure;
imshow(uint8(I_seg));
figure;
imhist(uint8(I));
**********************************************************
2维直方图熵阈值算法
% 二维直方图熵阈值法
% 参考 基于2D 熵阈值的铁谱磨粒图像分割方法，傅建平
%廖振强，张培林，汪传忠，（南京理工大学机械学院，南京 ），
%（军械工程学院，石家庄）
% pixel gray
% ^
% |
% | ==> 2D histgram
% |
% |
% |__________________________> local gray
clc;
clear;
%I = imread(‘E:/test/chinalake.bmp’,’bmp’);
I = imread(‘E:/test/lena.png’,’png’);
I = double(I);
[height,width] = size(I);
Size = height * width; % the size of the image
G_min = min(min(I)); % the min gray value of the image
G_max = max(max(I)); % the max gray value of the iamge
I_2Dhist = zeros(G_max+1,G_max+1); % the array to store the 2D hist of the image
I_mean = zeros(height,width); % the mean value of the local image
I_seg = zeros(height,width);
WS = 3; % mean filter’s window size 3*3
nr = floor(WS/2);
I_big = zeros(height+2*nr,width+2*nr); % the bigger array used to mean filter
I_big(nr+1:height+nr,nr+1:width+nr) = I; % copy data from the original image
%%%%%%%%%% mean filter %%%%%%%%%%%
%%%%%% 获取局部区域灰度信息 %%%%%
for i=1:height
for j=1:width
sum = 0;
num = 0;
for h=-nr:nr
for w=-nr:nr
sum = sum + I_big(i+h,j+w);
num = num + 1;
end
end
I_mean(i,j) = sum/num;
end
end
%%%% 构建2D直方图，横轴上以点象素灰度表示，纵轴上以局部区域灰度表示 %%%%
for i=1:height
for j=1:width
h = I(i,j)+1; % 横轴信息,避免0,所以加1，象素灰度
w = I_mean(i,j)+1; % 纵轴信息，避免0，所以加1，局部区域灰度
I_2Dhist(h,w) = I_2Dhist(h,w) + 1; % 统计灰度对<pixel,local>的出现次数，构建2D直方图
end
end
% find the best thresh : hor_thresh,and ver_thresh %
for ver_thresh=0:G_max
for hor_thresh=0:G_max
sum1 = 0;
sum2 = 0;
H1 = 0;
H2 = 0;
for i=0:ver_thresh
for j=0:hor_thresh
sum1 = sum1 + I_2Dhist(i+1,j+1);
end
end
for i=0:ver_thresh
for j=0:hor_thresh
P1 = I_2Dhist(i+1,j+1)/sum1;
H1 = H1 + P1*log(P1);
end
end
if i < G_max & j < G_max
for i=ver_thresh+1:G_max
for j=hor_thresh+1:G_max
sum2 = sum2 + I_2Dhist(i+1,j+1);
end
end
for i=ver_thresh+1:G_max
for j=hor_thresh+1:G_max
P2 = I_2Dhist(i+1,j+1)/sum2;
H2 = H2 +P2*log(P2);
end
end
end
H_total = H1 + H2;
if H_total > max
max = H_total;
ver_Tbest = ver_thresh;
hor_Tbest = hor_thresh;
end
end
end
%%%%%% 分割 %%%%%%%%
for i=1:height
for j=1:width
if I(i,j) < ver_Tbest & I_mean(i,j) < hor_Tbest
I_seg(i,j) = 0;
elseif I(i,j) >= ver_Tbest & I_mean(i,j) >= hor_Tbest
I_seg(i,j) = 255;
end
end
end
T_best
figure;
imshow(uint8(I_seg));
figure;
imhist(uint8(I));
