【图像分割】基于Matlab Tsallis熵算法灰度图像分割【含Matlab源码 715期】

一、简介

在这里插入图片描述
二、源代码

%*********************************************
%Image thresholding using Tsallis entropy

%2021-3-26 
%*********************************************
clear all;
close all;
clc

format short g;
% disp(strcat('***************apply on whiteRose.jpg'));

% I=imread('whiteRose.jpg');
I=imread('C:\Users\lenovo\Desktop\113172244Tsallis-entropy\Tsallis entropy/cell.jpg');
% I=imread('ndmz.png');
% I=imread('star.jpg');
% I=imread('cell2.jpg');
imgray = rgb2gray(I);%灰度转化
thresh=graythresh(imgray);%求阈值
I2=im2bw(imgray,thresh);%二值图转化
subplot(2,3,1)
imshow(I);title('原图');
subplot(2,3,2)
imshow(imgray);title('灰度图');
subplot(2,3,3)
imhist(imgray);title('灰度直方图');
subplot(2,3,4)
omega=num2str(uint8(round(thresh*255)));
imshow(I2);title(['OTSU 阈值：',omega]);
disp(strcat('OTSU 计算灰度阈值：',num2str(uint8(round(thresh*255)))));
thresh1=thresh
[m,n]=size(imgray);
N=m*n ; %灰度值的总数
%****************************************
%KAPUR_entroy最大熵值法求阈值 figure(1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[count,x]=imhist(imgray,256);
A=[x,count];
p=count/N;%求每个灰度值出现的概率
%用'1985-histo'文献里提到的新算法计算阈值
a=min(find(p~=0.0));
P=p(p~=0);%提取p不为零的项，剔除奇点
p=p(a:256);%把前几项为零的奇点剔除
b=length(p);
ps=cumsum(p);%概率1-T求和
PS=cumsum(P);
Hs=-cumsum(P.*log(P));
HA=log(PS)+Hs./PS;
HB=log(1-PS)+(Hs(end)-Hs)./(1-PS); 
L2=length(HA);%剔除1-PS为零的点
HA=HA(1:(L2)-1);
HB=HB(1:(L2)-1);
Ts=HA+HB;
s=mean(find(Ts==max(Ts))+a-2);
disp(strcat('KAPUR_entroy计算灰度阈值：',num2str(s)));
thresh2=double(s)/255
I3=im2bw(imgray,double(s)/255);   
subplot(2,3,5);imshow(I3);title(['KAPUR 阈值：',num2str(s)]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Image thresholding using Tsallis entropy
q=0.5;
for T=1:1:b
    sumA=0.0;sumB=0.0;
      for i=1:1:T
       sumA=sumA+(p(i)/ps(T))^q;
      end
        for j=T+1:1:b
        sumB=sumB+(p(j)/(1-ps(T)))^q;
        end
    c(T)=sumA;
    d(T)=sumB;
end
L=min(find(d==0));%剔除d=0的奇点
c=c(1:L-1);
d=d(1:L-1);
 g=[c',d'];
    Sqa=(1-c)./(q-1);
    Sqb=(1-d)./(q-1);
    g=[Sqa',Sqb'];
    ST=Sqa+Sqb+(1-q).*Sqa.*Sqb;%Tsallis entropy定义式
    t=mean(find(ST==max(ST))+a-2);
disp(strcat('Tsallis entropy计算灰度阈值：',num2str(t)));
thresh3=double(t)/255
I3=im2bw(imgray,thresh3);   
subplot(2,3,6);imshow(I3);title(['Tsallis 阈值：',num2str(t)]);
%*****************************************
figure(2)                            %%%*************不同q值下的情况*******
%*****************************************

% disp(strcat('***************apply on star.jpg'));
% I=imread('whiteRose.jpg');
%  I=imread('ROSE2.jpg');
% I=imread('ndmz.png');
% I=imread('star.jpg');
% I=imread('inftest.jpg');
% I=imread('kejian.jpg');
I=imread('fengz.jpg');
% I=imread('cell.jpg');
%  I=imread('xingzuo.jpg');
imgray=rgb2gray(I);%灰度转化
thresh=graythresh(imgray);%求阈值
I2=im2bw(imgray,thresh);%二值图转化
subplot(2,3,1)
imshow(I);title('原图');
subplot(2,3,2)
imshow(imgray);title('灰度图');
subplot(2,3,3)
imhist(imgray);title('灰度直方图');
subplot(2,3,4)
omega=num2str(uint8(round(thresh*255)));
imshow(I2);title(['OTSU 阈值：',omega]);
disp(strcat('OTSU 计算灰度阈值：',num2str(uint8(round(thresh*255)))));
thresh1=thresh
[m,n]=size(imgray);
N=m*n ; %灰度值的总数
%****************************************
%KAPUR_entroy最大熵值法求阈值 figure(2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[count,x]=imhist(imgray,256);
A=[x,count];
p=count/N;%求每个灰度值出现的概率
%用'1985-histo'文献里提到的新算法计算阈值
a=min(find(p~=0.0));
P=p(p~=0);%提取不为零的项
p=p(a:256);
b=length(p);
ps=cumsum(p);%概率1-T求和
PS=cumsum(P);
Hs=-cumsum(P.*log(P));
HA=log(PS)+Hs./PS;
HB=log(1-PS)+(Hs(end)-Hs)./(1-PS);
L2=length(HA);%剔除1-PS为零的奇点
HA=HA(1:(L2)-1);
HB=HB(1:(L2)-1);
Ts=HA+HB;
s=mean(find(Ts==max(Ts))+a-2);
disp(strcat('KAPUR_entroy计算灰度阈值：',num2str(s)));
thresh2=double(s)/255
I3=im2bw(imgray,double(s)/255);   
subplot(2,3,5);imshow(I3);title(['KAPUR 阈值：',num2str(s)]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Image thresholding using Tsallis entropy ****figure(2)
q=[0.1,0.5,1,1.5,2];
for n=1:1:5
    if q(n)~=1%计算q/=1时候的情况Tsallis entropy
for T=1:1:b
    sumA=0.0;sumB=0.0;
    for i=1:1:T 
        sumA=sumA+(p(i)/ps(T))^q(n);
    end
    for j=T+1:1:b
        sumB=sumB+(p(j)/(1-ps(T)))^q(n);
    end
   
    c(T)=sumA;
    d(T)=sumB;
end
L=min(find(d==0));%处理d=0的奇点
c=c(1:L-1);
d=d(1:L-1);
 f=[c',d'];
    Sqa=(1-c)./(q(n)-1);
    Sqb=(1-d)./(q(n)-1);
    h=[Sqa',Sqb'];
    ST=Sqa+Sqb+(1-q(n)).*Sqa.*Sqb;
    ss=[f,h,ST'];
    find(ST==max(ST));
    t(n)=mean(find(ST==max(ST))+a-2);
    else
三、运行结果

在这里插入图片描述
四、备注

版本：2014a