www.pudn.com > ACA_FCM.rar > ACA_FCM.m
function ACA_FCM
%The input image should have square size
%All parameters are set to be exactly same as that of the paper
%Image Loading
filename = 'flower';
img = double(imread([filename '.jpg']));
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=0.40;
beta=3;
num=0;%监督聚类中心的初始数目，初始化为0
statistic=60;
radius=90;% 聚类半径
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)=e)
for i = 1 : C
sum_subjection1 = 0;
sum_subjection2 = 0;
sum_subjection3 = 0;
sum_subjection=0;
for rr = 1 : nrow
for cc = 1 : ncol
if ~(ant_matrix(rr, cc, 1)==4)
sum_subjection1 = sum_subjection1 + (subjection(rr, cc, i).^m) * R(rr, cc);
sum_subjection2 = sum_subjection2 + (subjection(rr, cc, i).^m) * G(rr, cc);
sum_subjection3 = sum_subjection3 + (subjection(rr, cc, i).^m) * B(rr, cc);
sum_subjection = sum_subjection + subjection(rr, cc, i).^m;
end
end
end
center_FCM(i, 1) = sum_subjection1/sum_subjection;
center_FCM(i, 2) = sum_subjection2/sum_subjection;
center_FCM(i, 3) = sum_subjection3/sum_subjection;
end
for rr = 1 : nrow
for cc = 1 : ncol
for i = 1 : C
subjection_temp(rr, cc, i) = subjection(rr, cc, i);
end
end
end
%隶属度矩阵计算
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
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;
end
end
end
end
end
%距离和计算
sum_d = 0;
for rr = 1 : nrow
for cc = 1 : ncol
for i = 1 : C
sum_d = sum_d + (subjection(rr, cc, i) - subjection_temp(rr, cc, i)).^2;
end
end
end
sum_d
end
%带有分类信息的图像矩阵
cluster_img2 = zeros(nrow, ncol, 4);
for rr = 1 : nrow
for cc = 1 : ncol
if ~(ant_matrix(rr, cc, 1)==4)
temp = max(subjection(rr, cc, :));
for i = 1 : C
if subjection(rr, cc, i) == temp
cluster_img2(rr, cc, 1) = center_FCM(i, 1);
cluster_img2(rr, cc, 2) = center_FCM(i, 2);
cluster_img2(rr, cc, 3) = center_FCM(i, 3);
cluster_img2(rr, cc, 4) = i;
end
end
elseif ant_matrix(rr, cc, 1)==4
cluster_img2(rr, cc, 1) = 255;
cluster_img2(rr, cc, 2) = 255;
cluster_img2(rr, cc, 3) = 255;
end
end
end
center_FCM
%cluster_img2(:, :, 4)
imwrite(cluster_img2(:, :, 1:3)./255, 'Result2.bmp', 'bmp');