【图像分割】基于matlab蚁群优化模糊聚类图像分割【含Matlab源码 130期】

绝代码农

　　
一、蚁群算法简介
　　蚁群算法(ant colony optimization, ACO)，又称蚂蚁算法，是一种用来在图中寻找优化路径的机率型算法。它由Marco Dorigo于1992年在他的博士论文中提出，其灵感来源于蚂蚁在寻找食物过程中发现路径的行为。蚁群算法是一种模拟进化算法，初步的研究表明该算法具有许多优良的性质。针对PID控制器参数优化设计问题，将蚁群算法设计的结果与遗传算法设计的结果进行了比较，数值仿真结果表明，蚁群算法具有一种新的模拟进化优化方法的有效性和应用价值。
　　定义
各个蚂蚁在没有事先告诉他们食物在什么地方的前提下开始寻找食物。当一只找到食物以后，它会向环境释放一种挥发性分泌物pheromone (称为信息素,该物质随着时间的推移会逐渐挥发消失，信息素浓度的大小表征路径的远近)来实现的，吸引其他的蚂蚁过来，这样越来越多的蚂蚁会找到食物。
　　有些蚂蚁并没有像其它蚂蚁一样总重复同样的路，他们会另辟蹊径，如果另开辟的道路比原来的其他道路更短，那么，渐渐地，更多的蚂蚁被吸引到这条较短的路上来。
　　最后，经过一段时间运行，可能会出现一条最短的路径被大多数蚂蚁重复着。
　　蚁群算法是一种仿生学算法，是由自然界中蚂蚁觅食的行为而启发的。在自然界中，蚂蚁觅食过程中，蚁群总能够按照寻找到一条从蚁巢和食物源的最优路径。下图显示了这样一个觅食的过程。
　　在图（a）中，有一群蚂蚁，假如A是蚁巢，E是食物源（反之亦然）。
　　这群蚂蚁将沿着蚁巢和食物源之间的直线路径行驶。假如在A和E之间突然出现了一个障碍物（图（b）），那么，在B点（或D点）的蚂蚁将要做出决策，到底是向左行驶还是向右行驶？由于一开始路上没有前面蚂蚁留下的信息素（pheromone），蚂蚁朝着两个方向行进的概率是相等的。但是当有蚂蚁走过时，它将会在它行进的路上释放出信息素，并且这种信息素会议一定的速率散发掉。信息素是蚂蚁之间交流的工具之一。它后面的蚂蚁通过路上信息素的浓度，做出决策，往左还是往右。很明显，沿着短边的的路径上信息素将会越来越浓（图（c）），从而吸引了越来越多的蚂蚁沿着这条路径行驶。

蚁群算法原理
假如蚁群中所有蚂蚁的数量为m，所有城市之间的信息素用矩阵pheromone表示，最短路径为bestLength，最佳路径为bestTour
每只蚂蚁都有自己的内存，内存中用一个禁忌表（Tabu）来存储该蚂蚁已经访问过的城市，表示其在以后的搜索中将不能访问这些城市
一个允许访问的城市表（Allowed）来存储它还可以访问的城市
矩阵（Delta）来存储它在一个循环（或者迭代）中给所经过的路径释放的信息素
还有一些额外数据，例如一些控制参数(α，β，ρ，Q)
蚂蚁行走完全程的总成本或距离（tourLength）
假定算法总共运行MAX_GEN次，运行时间为t。
算法流程说明（结合流程图食用，效果更佳）
（1）初始化
　　设t=0，初始化bestLength为一个非常大的数（正无穷），bestTour为空。初始化所有的蚂蚁的Delt矩阵所有元素初始化为0，Tabu表清空，Allowed表中加入所有的城市节点。随机选择它们的起始位置（也可以人工指定）。在Tabu中加入起始节点，Allowed中去掉该起始节点。
　　（2）为每只蚂蚁选择下一个节点
　　为每只蚂蚁选择下一个节点，该节点只能从Allowed中以某种概率（公式1）搜索到，每搜到一个，就将该节点加入到Tabu中，并且从Allowed中删除该节点。该过程重复n-1次，直到所有的城市都遍历过一次。遍历完所有节点后，将起始节点加入到Tabu中。此时Tabu表元素数量为n+1（n为城市数量），Allowed元素数量为0。接下来按照（公式2）计算每个蚂蚁的Delta矩阵值。最后计算最佳路径，比较每个蚂蚁的路径成本，然后和bestLength比较，若它的路径成本比bestLength小，则将该值赋予bestLength，并且将其Tabu赋予BestTour。
　　（3）更新信息素矩阵Delta
（4）检查终止条件，是否到达MAX_GEN次
（5）输出最优值bestlength
　　算法流程图


二、源代码

clc
clear all
close all
%The input image should have square size
%All parameters are set to be exactly same as that of the paper
 
%Image Loading
%filename = 'flower';
I=imread('1.jpg');
%img = double(imread(['C:\Users\yxw\Desktop.jpg']));
img=double(I);
figure(1)
imshow(img/255);
[nrow, ncol, channel] = size(img);
R=img(:, :, 1);
G=img(:, :, 2);
B=img(:, :, 3);
 
%将彩色图像转化成灰度图像
intensity_img=zeros(nrow, ncol);
for rr = 1 : nrow
  for cc = 1 : ncol
    intensity_img(rr, cc)=(((R(rr, cc)).^2+(G(rr, cc)).^2+(B(rr, cc)).^2).^0.5)/(3.^0.5);
  end
end
 
figure(2)
imshow(intensity_img/255);
 
%使用canny算子进行边缘检测
edge_img = edge(intensity_img, 'canny');
figure(3)
imshow(edge_img);
 
%average = mean(mean(img))/255;
 
%参数设置
alpha=1;
beta=1;
num=0;%监督聚类中心的初始数目，初始化为0
statistic=60;
radius=50;% 聚类半径
lumda=0.40;
rho=0.95;
p1=1;
p2=1;
p3=1;
d=50;
 
%ACA图像分割程序
 
%初始化带有分类信息的图像矩阵
cluster_img = zeros(nrow, ncol, 4);
for rr=1:nrow
  for cc=1:ncol
    cluster_img(rr, cc, 1) = img(rr, cc, 1);
    cluster_img(rr, cc, 2) = img(rr, cc, 2);
    cluster_img(rr, cc, 3) = img(rr, cc, 3);
    cluster_img(rr, cc, 4) = 0;
  end
end
 
%初始化蚂蚁归类操作矩阵
ant_matrix=zeros(rr, cc, 1);
for rr=1:nrow
  for cc=1:ncol
    if ant_matrix(rr, cc, 1) == 0
      ant_matrix(rr, cc, 1) = 2;%ant_matrix(rr, cc, 1)代表该蚂蚁的状态，"0"表示未被归类，"1"表示已经被归类，"2"表示等待被归类, "3"表示在本次循环中成为类, "4"表示该点为边缘像素点
    end
  end
end
 
for rr = 1 : nrow
  for cc = 1 : ncol
    if edge_img(rr, cc)==1
      ant_matrix(rr, cc, 1) = 4;%划分为边缘像素点
      cluster_img(rr, cc, 1)=255;
      cluster_img(rr, cc, 2)=255;
      cluster_img(rr, cc, 3)=255;%边缘像素全部设为白色;
    end
  end
end
 
%初始化监督聚类中心
%对图像颜色进行统计
 
color_statistic = zeros(1331, 5);%color_statistic(i, 1)存储像素数目，
%color_statistic(i, 2)、color_statistic(i, 3)、color_statistic(i, 4)分别存储颜色的三个分量，color_statistic(i, 5)存储是否被作为监督聚类中心的信息
 
for rr = 1 : nrow
  for cc = 1 : ncol
    %对每个颜色分量进行分段处理
    if R(rr, cc) < 12.75
      x=1;
    elseif R(rr, cc) >= 12.75 && R(rr, cc) < 38.25
      x=2;
    elseif R(rr, cc) >= 38.25 && R(rr, cc) < 63.75
      x=3;
    elseif R(rr, cc) >= 63.75 && R(rr, cc) < 89.25
      x=4;
    elseif R(rr, cc) >= 89.25 && R(rr, cc) < 114.75
      x=5;
    elseif R(rr, cc) >= 114.75 && R(rr, cc) < 140.25
      x=6;
    elseif R(rr, cc) >= 140.25 && R(rr, cc) < 165.75
      x=7;
    elseif R(rr, cc) >= 165.75 && R(rr, cc) < 191.25
      x=8;
    elseif R(rr, cc) >= 191.25 && R(rr, cc) < 216.75
      x=9;
    elseif R(rr, cc) >= 216.75 && R(rr, cc) < 241.25
      x=10;
    elseif R(rr, cc) >= 241.25
      x=11;
    end
 
    if G(rr, cc) < 12.75
      y=1;
    elseif G(rr, cc) >= 12.75 && G(rr, cc) < 38.25
      y=2;
    elseif G(rr, cc) >= 38.25 && G(rr, cc) < 63.75
      y=3;
    elseif G(rr, cc) >= 63.75 && G(rr, cc) < 89.25
      y=4;
    elseif G(rr, cc) >= 89.25 && G(rr, cc) < 114.75
      y=5;
    elseif G(rr, cc) >= 114.75 && G(rr, cc) < 140.25
      y=6;
    elseif G(rr, cc) >= 140.25 && G(rr, cc) < 165.75
      y=7;
    elseif G(rr, cc) >= 165.75 && G(rr, cc) < 191.25
      y=8;
    elseif G(rr, cc) >= 191.25 && G(rr, cc) < 216.75
      y=9;
    elseif G(rr, cc) >= 216.75 && G(rr, cc) < 241.25
      y=10;
    elseif G(rr, cc) >= 241.25
      y=11;
    end
 
    if B(rr, cc) < 12.75
      z=1;
    elseif B(rr, cc) >= 12.75 && B(rr, cc) < 38.25
      z=2;
    elseif B(rr, cc) >= 38.25 && B(rr, cc) < 63.75
      z=3;
    elseif B(rr, cc) >= 63.75 && B(rr, cc) < 89.25
      z=4;
    elseif B(rr, cc) >= 89.25 && B(rr, cc) < 114.75
      z=5;
    elseif B(rr, cc) >= 114.75 && B(rr, cc) < 140.25
      z=6;
    elseif B(rr, cc) >= 140.25 && B(rr, cc) < 165.75
      z=7;
    elseif B(rr, cc) >= 165.75 && B(rr, cc) < 191.25
      z=8;
    elseif B(rr, cc) >= 191.25 && B(rr, cc) < 216.75
      z=9;
    elseif B(rr, cc) >= 216.75 && B(rr, cc) < 241.25
      z=10;
    elseif B(rr, cc) >= 241.25
      z=11;
    end
 
    %更新统计信息
    color_statistic(((x-1)*121+(y-1)*11+z), 2) = (color_statistic(((x-1)*121+(y-1)*11+z), 2) * color_statistic(((x-1)*121+(y-1)*11+z), 1) + R(rr, cc)) / (color_statistic(((x-1)*121+(y-1)*11+z), 1) + 1);
    color_statistic(((x-1)*121+(y-1)*11+z), 3) = (color_statistic(((x-1)*121+(y-1)*11+z), 3) * color_statistic(((x-1)*121+(y-1)*11+z), 1) + G(rr, cc)) / (color_statistic(((x-1)*121+(y-1)*11+z), 1) + 1);
    color_statistic(((x-1)*121+(y-1)*11+z), 4) = (color_statistic(((x-1)*121+(y-1)*11+z), 4) * color_statistic(((x-1)*121+(y-1)*11+z), 1) + B(rr, cc)) / (color_statistic(((x-1)*121+(y-1)*11+z), 1) + 1);
    color_statistic(((x-1)*121+(y-1)*11+z), 1) = color_statistic(((x-1)*121+(y-1)*11+z), 1) + 1;
 
  end
end
 
for i = 1 : 1331
  if color_statistic(i, 1) >= statistic
    num = num +1;%确定监督聚类中心的数目
    color_statistic(i, 5) = 1;
  end
end
 
center_ACA=zeros(num, 5);
for i=1 : num
  for j = 1 : 1331
    if (color_statistic(j, 1) >= statistic) && (color_statistic(j, 5)==1)
      center_ACA(i, 1) = color_statistic(j, 2);
      center_ACA(i, 2) = color_statistic(j, 3);
      center_ACA(i, 3) = color_statistic(j, 4);
      color_statistic(j, 5) = 0;%已被作为监督聚类中心
      center_ACA(i, 5)=1;%center_ACA(i, 4)用于储存该类具有的像素数，初始化为0,center_ACA(i, 5)用于储存该类具有的信息素浓度
      break
    end
  end
end
 
ant_info = zeros (num, 5);
%ant_info((r-1)*ncol+c, 1) 储存到该点的距离；ant_info((r-1)*ncol+c, 2) 储存该点的信息素；ant_info((r-1)*ncol+c,3) 储存到该点的启发函数
%ant_info((r-1)*ncol+c, 4) 储存到该点的概率；ant_info((r-1)*ncol+c, 5)储存类别
 
 
%主程序
sum_ph=0;%概率公式的分母
do=0;%判断主程序是否应继续循环
do2=1;%判断类合并程序是否应继续循环
judge=zeros(nrow, ncol);%判断该类是否与别的类可以合并的矩阵
n=0;%判断ACA主程序循环次数
m=0;%判断类合并程序循环次数
 
 
while (do<=500)
 
  %像素点聚类
  for rr = 1 : nrow
    for cc = 1 : ncol
      if ant_matrix(rr, cc, 1)==0
        ant_matrix(rr, cc, 1) = 2;%将前一个循环里未被归类的蚂蚁状态都改为待聚类
      end
    end
  end
 
  for rr = 1 : nrow
    for cc = 1 : ncol
      %计算像素(rr, cc)到各聚类中心的距离、信息素等多项信息
      if ant_matrix(rr, cc, 1)==2
        ant_info = zeros(num, 6);
        sum_ph=0;
        for i = 1 : num
 
          ant_info(i, 1) = sqrt(p1 * (R(rr, cc) - center_ACA(i, 1)).^2 + p2 * (G(rr, cc) - center_ACA(i, 2)).^2 + p3 * (B(rr, cc) - center_ACA(i, 3)).^2);
 
          ant_info(i, 2) = center_ACA(i, 5);
 
          ant_info(i, 3)=radius/(ant_info(i, 1)+0.00001);%保证距离为0的时候，启发函数很大但不为无穷。
 
          sum_ph = sum_ph + ant_info(i, 2).^alpha + ant_info(i, 3).^beta;
 
          ant_info(i, 5) = i;
 
        end
 
        %计算到各聚类中心的概率
        for i=1:num
          ant_info(i, 4) = (ant_info(i, 2).^alpha + ant_info(i, 3).^beta)/sum_ph;
        end
 
        rand('state', sum(100*clock));
        temp = find(cumsum(ant_info(:, 4)) >= rand(1), 1);
        %路径的概率选择计算
 
        if ant_info(temp, 4) >= lumda
 
          ant_matrix(rr, cc, 1) = 1;%该像素已经被归类
 
          cluster_img(rr, cc, 4) = temp;%记录该像素的类别
 
          center_ACA(temp, 4) = center_ACA(temp, 4) + 1;%该聚类中包含的像素数加1
 
          if ant_info(temp, 1) <= radius
            center_ACA(temp, 5) = center_ACA(temp, 5) + 1;%若距离小于radius，则信息素加1
          end
 
        else
          ant_matrix(rr, cc, 1) = 0;%像素未被归类，状态变为0
        end
      end
    end
  end
 
  %更新聚类中心信息
  for i = 1 : num
    if ~(center_ACA(i, 4)==0)
      sum1=0;
      sum2=0;
      sum3=0;
      for rr = 1 : nrow
        for cc = 1 : ncol
          if cluster_img(rr, cc, 4)==i
            sum1 = sum1 + cluster_img(rr, cc, 1);
            sum2 = sum2 + cluster_img(rr, cc, 2);
            sum3 = sum3 + cluster_img(rr, cc, 3);
          end
        end
      end
      center_ACA(i, 1) = sum1 / center_ACA(i, 4);
      center_ACA(i, 2) = sum2 / center_ACA(i, 4);
      center_ACA(i, 3) = sum3 / center_ACA(i, 4);
    end
  end
 
  %类间合并
  %初始化判断矩阵
  while(do2==1)
    judge=zeros(num, 1);
    for i = 1 : num
      if ~(center_ACA(i, 4)==0)
        judge(i, 1)=1;
      end
    end
 
    for i = 1 : num
      if ~(center_ACA(i, 4)==0)
        cluster_info = zeros(num, 2);%记录类间距离
        temp=[d; 0];%第一个记录上次的距离值，第二个记录类别
        for j = 1 : num
          if (~(j==i))&&(~(center_ACA(j, 4)==0))
            cluster_info(j, 1) = sqrt((center_ACA(i, 1) - center_ACA(j, 1)).^2 + (center_ACA(i, 2) - center_ACA(j, 2)).^2 + (center_ACA(i, 3) - center_ACA(j, 3)).^2);
            cluster_info(j, 2) = j;
          end
        end
        for j = 1 : num
          if cluster_info(j, 1)<temp(1,1) && (~(j==i)) && (~(center_ACA(j, 4)==0))
            temp(1, 1)=cluster_info(j, 1);
            temp(2, 1)=cluster_info(j, 2);
            %计算最相近的类
          end
        end
 
        if temp(1,1)<d
          for rr = 1 : nrow
            for cc = 1 : ncol
              if cluster_img(rr, cc, 4)==i;
                cluster_img(rr, cc, 4) = temp(2,1);
                %像素分类矩阵中，(rr, cc)点的像素被归类
              end
            end
          end
 
          center_ACA(temp(2, 1), 1) = (center_ACA(temp(2, 1), 1) * center_ACA(temp(2, 1), 4) + center_ACA(i, 1) * center_ACA(i, 4)) / (center_ACA(temp(2, 1), 4) + center_ACA(i, 4));
          center_ACA(temp(2, 1), 2) = (center_ACA(temp(2, 1), 2) * center_ACA(temp(2, 1), 4) + center_ACA(i, 2) * center_ACA(i, 4)) / (center_ACA(temp(2, 1), 4) + center_ACA(i, 4));
          center_ACA(temp(2, 1), 3) = (center_ACA(temp(2, 1), 3) * center_ACA(temp(2, 1), 4) + center_ACA(i, 3) * center_ACA(i, 4)) / (center_ACA(temp(2, 1), 4) + center_ACA(i, 4));
          center_ACA(temp(2, 1), 4) = center_ACA(temp(2, 1), 4) + center_ACA(i, 4);
          center_ACA(temp(2, 1), 5) = center_ACA(temp(2, 1), 5) + center_ACA(i, 5);%max(center_ACA(temp(2, 1), 5), center_ACA(i, 5)) + 0.3 * min(center_ACA(temp(2, 1), 5), center_ACA(i, 5));
          center_ACA(i, 4) = 0;
          center_ACA(i, 5) = 0;
 
          judge(i, 1)=0;
          judge(temp(2, 1), 1)=1;
 
        else
          judge(i, 1)=0;
        end
      end
    end
 
    do2=0;
 
    for i = 1 : num
      if judge(i, 1) == 1
        do2=1;
        break
      end
    end
  end
 
  %信息素挥发
  for i = 1 : num
    center_ACA(i, 5) = center_ACA(i, 5) * rho;
  end
  do = do + 1 %每循环一次，会显示一次do
end
 
for rr = 1 : nrow
  for cc = 1 : ncol
    if ant_matrix(rr, cc, 1)==1
      cluster_img(rr, cc, 1) = center_ACA(cluster_img(rr, cc, 4), 1);
      cluster_img(rr, cc, 2) = center_ACA(cluster_img(rr, cc, 4), 2);
      cluster_img(rr, cc, 3) = center_ACA(cluster_img(rr, cc, 4), 3);
    elseif ant_matrix(rr, cc, 1)==0
      cluster_img(rr, cc, 1) = 0;
      cluster_img(rr, cc, 2) = 0;
      cluster_img(rr, cc, 3) = 0;
    end
  end
end
%figure(),imshow(cluster_img(:, :, 1:3)./255);
imwrite(cluster_img(:, :, 1:3)./255, 'Result1.bmp', 'bmp');
 
 
%FCM主程序
C = num;%类的数目为C
m = 2;%类参数设置
e = nrow * ncol * C * (0.01).^2;
sum_d = e+1;
 
%初始化类矩阵
center_FCM = zeros(C, 3);
for i = 1 : C
  center_FCM(i, 1) = center_ACA(i, 1);
  center_FCM(i, 2) = center_ACA(i, 2);
  center_FCM(i, 3) = center_ACA(i, 3);
end
 
%初始化距离矩阵
distance=zeros(C, 1);
 
%利用ACA运行结果初始化隶属度矩阵
subjection = zeros(nrow, ncol, C);
subjection_temp = zeros(nrow, ncol, C);
for rr = 1 : nrow
  for cc = 1 : ncol
    if ~(ant_matrix(rr, cc, 1)==4)
      for i = 1 : C
        distance(i, 1) = sqrt((R(rr, cc)-center_FCM(i, 1)).^2 + (G(rr, cc)-center_FCM(i, 2)).^2 + (B(rr, cc)-center_FCM(i, 3)).^2);
      end
 
      do = 1;
      for i = 1 : C
        if distance(i, 1) == 0
          subjection(rr, cc, i) = 1;%若某个像素到聚类中心的距离为0，则其隶属度为1
          for j = 1 : C
            if ~(j==i)
              subjection(rr, cc, j) = 0;
              do = 0;
            end
          end
        end
        break
      end
 
      if do == 1
        normalize = 0;
        for i = 1 : C
          sum_distance = 0;
          for j = 1 : C
            sum_distance = sum_distance + (distance(i, 1)/distance(j, 1)).^(2/(m-1));
          end
          subjection(rr, cc, i) = 1 / sum_distance;
          normalize = normalize + subjection (rr, cc, i);
        end
        for i = 1 : C
          subjection(rr, cc, i) = (1 / normalize) * subjection(rr, cc, i);%隶属度归一化
        end
      end
    end
  end
end
 
 
end
figure
imshow(cluster_img2(:, :, 1:3)./255)
% center_FCM
%cluster_img2(:, :, 4)
imwrite(cluster_img2(:, :, 1:3)./255, 'Result2.bmp', 'bmp');

三、运行结果

四、备注
　　版本：2014a