数字图像处理程序合集下(matlab实现)

---

前言

本篇文章什数字图像小程序合集的最后一部分，主要包含影像的特征提取和影像分类两块内容，这部分的内容较难，请大家细细品味，有问题可以在评论区提出，我会及时回答的！

---

1.影像特征提取

1.1 迭代阈值影像分割算法

function B=diedaiyuzhi(filename,T0,N,sita0)% 输入参数（文件名，初始阈值，迭代次数，终止阈值）
% 读入图像文件
A=double(imread(filename));
% 将图像矩阵灰度化
A=rgb2gray(A);
[row,col]=size(A);
% 将三维矩阵化为一维向量便于处理
A=A(:);
% list1，list2存储两个类别的数据
list1=[];
list2=[];
% 赋初始阈值
T1=T0;
% 进入循环
for j=1:N
    for i=1:length(A)
        if A(i)<=T1  % 根据阈值将灰度数据划分
            list1=[list1;A(i)];
        else
            list2=[list2;A(i)];
        end
    end
    % 计算两个类别的平均值，再重新计算阈值
    av1=mean(list1);
    av2=mean(list2);
    T0=T1;
    T1=(av1+av2)/2;
    % 判断迭代终止条件，两个阈值差小于一个给定的值，停止迭代
    if T1-T0<=sita0
        break
    end
end
% 将影像二值化
A(A<=T1)=0;
A(A>T1)=1;
% 重新转为三维矩阵便于展示
B=reshape(A,[row,col]);
% 设置参数并且展示分割后的影像
figure()
imshow(uint8(B))
title('阈值分割后的图像')
grid on
end

---

1.2 区域生长影像分割算法

function B=quyuzengzhang(filename,point,e) % 输入参数（文件名，生长种子的像素坐标，生长范围）
% 将影像读取并且灰度化
A=rgb2gray(double(imread(filename)));
[row,col]=size(A);
% 初始化标志矩阵
flag=zeros([row,col]);
list=[point(1),point(2)];% 生长集合中首先放入种子点坐标
% 得到生长点的8邻域
quyu=A(point(1)-1:point(1)+1,point(2)-1:point(2)+1);
% find函数查找邻域内与种子点灰度差<=e的点
[rows,cols]=find(quyu>=A(point(1),point(2)-e)&&quyu<=A(point(1),point(2)+e));
% 将此点所在位置的标识变为1
flag(quyu>=A(point(1),point(2)-e)&&quyu<=A(point(1),point(2)+e))=1;
% 将此点加入生长集合
list1=[list;rows+point(1)-2,cols+point(2)-2];
% 区域不封闭，继续迭代
while list1~=list
    % 种子点变为外侧新找到的点
    list=list1;
    % 生长范围更新为生长集合的平均值
    e=mean(A(list(:,1),list(:,2)));
    for i=1:length(list)
         % 重复上述过程生长区域不断延伸
         quyu=A(list(i,1)-1:list(i,1)+1,list(i,2)-1:list(i,2)+1);
         [rows,cols]=find(quyu>=A(list(i,1),list(i,2)-e)&&quyu<=A(list(i,1),liBst(i,2)+e));
         flag(quyu>=A(list(i,1),list(i,2)-e)&&quyu<=A(list(i,1),list(i,2)+e))=1;
         list1=[list;quyu(rows,cols)];
    end
end
% 最终的标志矩阵为一个二值矩阵
B=flag;
% 对分割结果进行展示
figure()
imshow(uint8(B))
title('区域增长后的图像')
grid on
end

---

2.影像分类

2.1 K-means聚类

function B = kmeans(A,K)% Nmeans聚类，A为待处理影像，K为自定的聚类数
[lines,samples,~] = size(A);%默认是三个波段
A = reshape(A(:,:,:),lines*samples,3);
Benter = zeros(K,3); % 最终的聚类中心
Benter1=zeros(K,3);
%初始化聚类中心
for N = 1:K
    Benter(N,:)= round(255*N/N);
end
Result = zeros(lines*samples,1); % 类别判定矩阵
Data=zeros(N,1);
while(Benter~=Benter1) % 迭代条件，两次的聚类中心不相等继续迭代
    Benter1 = Benter;
    %按最短距离聚类
    for i = 1:lines*samples
        for N = 1:K
            Data(N) = sqrt((A(i,1)-Benter(N,1))^2+(A(i,2)-Benter(N,2))^2+(A(i,3)-Benter(N,3))^2); % 计算每个像素点到各个聚类中心的欧式距离
        end
        [~,index] = min(Data); % 找到最小值的索引
        Result(i,1) = index; % 将此点划分为第i个类别
    end
    %重新计算聚类中心
    for N = 1:K
        % 找到所有属于第N个类别的像素点计算灰度平均值
        indexx = find(Result == N);
        Benter(N,3) = mean(A(indexx,1)+A(indexx,2)+A(indexx,3));  
    end
end
% 将原始图像与聚类分割后图像对比展示
figure(1);
subplot(1,2,1);
imshow(uint8(A));
title('原始图像');
B = reshape(Result,lines,samples);
subplot(1,2,2);
imshow(uint8(B));
title('K-means聚类分割后的影像');
end

---

2.2 ISODATA分类

这个分类较难，相当于是K-means分类的进阶版，包含了聚类的合并和分割，类别的不确定性以及自动更新是程序的难点，一般考试只作为了解即可以，但是原理一定要掌握！

function ISODATA(x,K,theta_N,theta_S,theta_c,L,I)
%% 输入参数
% x : data
% K : 预期的聚类中心数
% theta_N : 每一聚类中心中最少的样本数，少于此数就不作为一个独立的聚类
% theta_S ：一个聚类中样本距离分布的标准差
% theta_c : 两聚类中心之间的最小距离，如小于此数，两个聚类进行合并
% L : 在一次迭代运算中可以和并的聚类中心的最多对数
% I ：迭代运算的次数序号
%% 开始
%% step1：
% 获取类别数，并且得到各个类别的平均值
n = size(x,1);
N_c = K;
mean = cell(K,1);
for i=1:K
    mean{i} = x(i,:);
end
% 开始迭代，设置迭代标志
ite = 1;
while ite<I
    flag = 1;
    while flag
    %% step2
    class = cell(size(mean));
    for i=1:n
        num = Belong2(x(i,:),mean);
        class{num} =  [class{num};x(i,:)];
    end
    %% step3
    for i=1:N_c
        size_i = size(class{i},1);
        if size_i<theta_N
          class_i = class{i};
          mean = DeleteRow(mean,i);
          class = DeleteRow(class,i);
          N_c = N_c-1;
          for j=1:size_i
            class_ij = class_i(j,:);%the j'th row of class{i}
            num = Belong2(class_ij,mean);
            class{num} = [class{num};class_ij];
          end
        end
    end
 
    %% step4
    for i=1:N_c
        if ~isempty(mean{i})
            mean{i} = sum(class{i})./size(class{i},1);
        end
    end
    %% step5
    Dis = zeros(N_c,1);
    for i=1:N_c
        if ~isempty(class{i})
            N_i =size(class{i},1);
            tmp = bsxfun(@minus,class{i},mean{i});
            Dis(i) = sum(arrayfun(@(x)norm(tmp(x,:)),1:N_i))/N_i;
        end
    end
    %% step6
    D = 0;
    for i=1:N_c
        if ~isempty(class{i})
            N_i =size(class{i},1);
            D = D + N_i*Dis(i);
        end
    end
    D = D/n;
    %% step7
    flag = 0;
    if ite == I
        theta_c = 0;
        flag = 0;
    elseif ~(N_c > K/2)
        flag = 1;
    elseif mod(ite,2)==0 || ~(N_c<2*K)
        flag = 0;
    end
    %% 分裂处理
    %% step8
    if flag
        flag = 0;
        delta = cell(N_c,1);
        for i=1:N_c
            if ~isempty(class{i})
                 N_i =size(class{i},1);
                 tmp = bsxfun(@minus,class{i},mean{i});
                 delta{i} = arrayfun(@(x)norm(tmp(:,x)),1:size(tmp,2))/N_i;
            end
        end
 
    %% step9
    delta_max = cell(N_c,1);
    for i=1:N_c
        if ~isempty(class{i})
            max_i = max(delta{i});
            sub = find(delta{i}==max_i,1);
            delta_max{i} = [max_i,sub];
        end
    end
    %% step10  
    for i=1:N_c
        if delta_max{i}(1) > theta_S
            N_i =size(class{i},1);
            con1 = (Dis(i)>D && N_i>2*(theta_N + 1));
            con2 = ~(N_c>K/2);
            if con1 || con2
               %%%%这里分裂%%%%%
               flag = 1;%一旦发生分裂，那么分裂一次后就返回第二步；若没发生分裂，则直接进入合并处理步
               lamda = 0.5;
               max_sub = delta_max{i}(2);
               mean{i}(max_sub) = mean{i}(max_sub) + lamda * delta_max{i}(1);
               addOneMean =  mean{i};
               addOneMean(max_sub) = addOneMean(max_sub) - lamda * delta_max{i}(1);
               mean = [mean;addOneMean];
               N_c = N_c+1;
               break;
            end
        end
     end
 
    end
 
    end
    %% 合并处理
    if L
    %% step11
    Distance = zeros(N_c,N_c);
    for i=1:N_c-1
        for j=i:N_c
            Distance(i,j) = norm(mean{i}-mean{j});
        end
    end
    %% step12
    index = find(-Distance>theta_c);
    keepIndex = [Distance(index),index];
    [~, index] = sort(keepIndex(:,1));
    if size(index,1) > L
        index = index(1:L,:);
    end
    %% step13
    if size(index,1) ~= 0
        for id=1:size(index,1)
            [m_i, m_j]= seq2idx(index(id),N_c);
            %%%%%这里合并%%%%%
            N_mi = size(class{m_i},1);
            N_mj = size(class{m_j},1);
            mean{m_i} = (N_mi*mean{m_i} + N_mj*mean{m_j})/(N_mi+N_mj);
            mean = DeleteRow(mean,m_j);
            class{m_i} = [class{m_i};class{m_j}];
            class = DeleteRow(class,m_j);
        end  
    end
    end
    %% step14
    ite=ite+1;
end
   for  i=1:N_c
       fprintf('第%d类聚类中心为\n',i);
       disp(mean{i});
       fprintf('第%d类中元素为\n',i);
       disp(class{i});
   end
end
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function number = Belong2(x_i,means)
    INF = 10000;
    min = INF;
    kk = size(means,1);
    number = 1;
    for i=1:kk
        if ~isempty(means{i})
            if norm(x_i - means{i}) < min
                min = norm(x_i - means{i});
                number = i;
            end
        end
    end
end
 
function A_del = DeleteRow(A,r)
    n = size(A,1);
    if r == 1
        A_del = A(2:n,:);
    elseif r == n
        A_del = A(1:n-1,:);
    else
        A_del = [A(1:r-1,:);A(r+1:n,:)];
    end
end
 
 
function [row, col] = seq2idx(id,n)
    if mod(id,n)==0
        row = n;
        col = id/n;
    else
        row = mod(id,n);
        col = ceil(id/n);
    end
end

---

2.3 最短距离分类（马氏距离）

function B = zuiduan(A,yb) % 马氏最短距离分类，A是待分类的影像，yb是样本(本质为元胞数组)
[m,n,p]=size(A); %获取图像大小
A=double(A); %转double型数据
[k,~]=size(yb); %获取样本个数
a=reshape(A,m*n,p); %重组图像成为二维数组
B=zeros(m*n,1); %初始化
for i=1:m*n %遍历所有像元
    %首先计算初始值
    a1=a(i,:); 
    %计算第一类别的中心
    ave=mean(yb{1,:});
    %计算第一类别的协方差
    center=inv((cov(yb{1,:}))); 
    %计算马氏距离
    d1=sqrt((a1-ave)*center*(a1-ave)'); 
    p=1;  
     %遍历其他类别，找距离最近的类
    for j=2:k
        ave=mean(yb{j,:}); 
        %计算类别的协方差
        center=inv((cov(yb{j,:}))); 
        %计算第一类别的中心
        %计算马氏距离
        d=sqrt((a1-ave)*center*(a1-ave)'); 
        %找距离最近的类 如果先前样本距离较大，则替换成当前样本的分类
        if d1>d 
            p=j;
            d1=d;
        end      
    end   
    B(i,1)=p; %存放分类的标签结果
end
% 对分类结果进行展示
B=reshape(B,m,n);
figure()
imshow(unit8(B))
title('马氏距离分类结果')
end

---

2.4 最大似然分类

function B= zuidasiran(A,yb ) % 最大似然分类方法，A为待分类影像，yb为已处理的样本(本质为元胞数组)
%yb为已经选取好的分类样本,为n*1的元胞数组
[m,n,p]=size(A); 
[k,~]=size(yb); %从样本中获取类别数
a=reshape(A,m*n,p); %将三维矩阵化为二维数组，便于后续处理
% 初始化
B=zeros(m*n,1);% 初始化类别标签矩阵
junzhi=zeros(k,1);% 初始化每个类别的平均值 
xfc=zeros(k,1);% 初始化每个类别的协方差
for i=1:k 
    junzhi(i,1)=mean(yb{i,1}); %计算各个类别的平均值
    xfc(i,1)=xfc(yb{i,1});  %计算各个类别的协方差
end
for i=1:m*n 
    % 将第i个波段的数据单独提取出来
    a1=a(i,:); 
    %利用最大似然公式计算像元属于第一类别的概率
    f1=exp(-(a1-junzhi(1,1))^2/(2*xfc(1,1)))/sqrt(2*pi*xfc(1,1)); 
    p=1; % 标识像元属于第几个类别
    for j=2:k
        %利用最大似然公式计算像元属于其他类别的概率
        f2=exp(-(a1-junzhi(j,1))^2/(2*xfc(j,1)))/sqrt(2*pi*xfc(j,1)); 
        if f1<f2 %如果下一个分类概率大
            p=j; %将这个类别的类别号记录下来
            f1=f2; 
        end
    end    
    if f1<0.3 
        B(i,1)=k+1; % 如果概率小于0.3，这一个像元可能不属于任何类别
    else %如果概率大于0.3
        B(i,1)=p;  % 最大似然检验通过，将标签矩阵更新
    end
end
B=reshape(B,m,n); % 重新将标签矩阵转为三维，为了展示
Asc(B); %此处应用Asc以全色图的方式显示分类的结果
end

---

总结

总的来说，这个合集中包含了数字图像处理这门课程中所有的基础程序，当然还有些比较复杂的程序，例如影像特征提取的两个程序的变化性其实很大，大家可以在食用的基础做出自己的改进！创作不易，还请大家点点赞，让我们下篇文章再见！