Ant colony
brief introduction
Ant colony optimization (ACO), also known as ant algorithm, is a probabilistic algorithm used to find the optimal path in the graph. It was proposed by Marco Dorigo in his doctoral thesis in 1992. Its inspiration comes from the behavior of ants in finding paths in the process of looking for food. Ant colony algorithm is a simulated evolutionary algorithm. Preliminary research shows that the algorithm has many excellent properties. Aiming at the parameter optimization design of PID controller, the results of ant colony algorithm design are compared with those of genetic algorithm design. The numerical simulation results show that ant colony algorithm has the effectiveness and application value of a new simulated evolutionary optimization method.
definition
Each ant starts looking for food without telling them where the food is in advance. When an ant finds food, it will release a volatile secretion pheromone (called pheromone, which will gradually volatilize and disappear with the passage of time, and the pheromone concentration represents the distance of the path) to attract other ants, so that more and more ants will find food.
Some ants do not always repeat the same path as other ants. They will find another way. If the other path is shorter than the original one, then gradually, more ants are attracted to this shorter path.
Finally, after a period of operation, the shortest path may be repeated by most ants.
Ant colony algorithm is a bionic algorithm, which is inspired by the foraging behavior of ants in nature. In nature, in the process of ant foraging, ant colony can always find an optimal path from ant nest and food source. The following figure shows such a foraging process.
In figure (a), there is a group of ants. If a is the ant nest, E is the food source (and vice versa).
The ants will follow A straight path between the nest and the food source. If an obstacle suddenly appears between A and E (Figure (b)), then the ant at point B (or D) will have to make A decision, whether to drive left or right? Since there is no pheromone left by the ants in front at the beginning, the probability of ants moving in both directions is equal. But when an ant walks by, it will release pheromones on its way, and this pheromone will be released at A certain rate. Pheromone is one of the communication tools between ants. The ants behind it make decisions, left or right, based on the concentration of pheromones on the road. More and more ants will travel along the path c, which will attract more and more information.
Principle of ant colony algorithm
- If the number of all ants in the ant colony is m, the pheromone between all cities is represented by matrix pheromone. The shortest path is bestLength and the best path is bestTour
- Each ant has its own memory, in which a Tabu is used to store the cities that the ant has visited, indicating that it will not be able to visit these cities in the future search
- An Allowed city table to store the cities it can also access
- Matrix (Delta) to store the pheromone it releases to the path it passes through in a loop (or iteration)
- There are also some additional data, such as some control parameters( α，β，ρ， Q)
- The total cost or tourLength of the ant walking the whole journey
- Assume that the algorithm runs Max in total_ Gen times, running time t.
Description of algorithm flow (better effect if eaten in combination with flow chart)
(1) Initialize
Set t=0, initialize bestLength to a very large number (positive infinity), and bestTour is empty. Initialize the Delt matrix of all ants, initialize all elements to 0, clear the Tabu table, and add all city nodes to the Allowed table. Randomly select their starting positions (which can also be specified manually). Add a starting node in Tabu and remove it from Allowed.
(2) Select the next node for each ant
Select the next node for each ant. The node can only be searched with a certain probability (formula 1) from Allowed. For each one, add the node to tabu and delete the node from Allowed. The process is repeated n-1 times until all cities have been traversed once. After traversing all nodes, add the starting node to tabu. At this time, the number of elements in the tabu table is n+1 (n is the number of cities), and the number of Allowed elements is 0. Next, calculate the Delta matrix value of each ant according to (formula 2). Finally, calculate the best path, compare the path cost of each ant, and then compare it with bestLength. If its path cost is smaller than bestLength, give the value to bestLength and Its Tabu to BestTour.
(3) Update pheromone matrix Delta
(4) Check whether the termination condition reaches MAX_GEN times
(5) Output optimal value bestlength
algorithm flow chart
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);
%Convert color image into gray image
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);
%use canny Operator for edge detection
edge_img = edge(intensity_img, 'canny');
figure(3)
imshow(edge_img);
%average = mean(mean(img))/255;
%Parameter setting
alpha=1;
beta=1;
num=0;%The initial number of supervised clustering centers is initialized to 0
statistic=60;
radius=50;% Cluster radius
lumda=0.40;
rho=0.95;
p1=1;
p2=1;
p3=1;
d=50;
%ACA Image segmentation program
%Initialize the image matrix with classification information
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
%Initialize ant classification operation matrix
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)Represents the state of the ant,"0"Indicates not classified,"1"Indicates that it has been classified,"2"Indicates waiting to be classified, "3"Indicates that it becomes a class in this loop, "4"Indicates that the point is an edge pixel
end
end
end
for rr = 1 : nrow
for cc = 1 : ncol
if edge_img(rr, cc)==1
ant_matrix(rr, cc, 1) = 4;%Divided into edge pixels
cluster_img(rr, cc, 1)=255;
cluster_img(rr, cc, 2)=255;
cluster_img(rr, cc, 3)=255;%Set all edge pixels to white;
end
end
end
%Initialize supervised clustering center
%Statistics of image color
color_statistic = zeros(1331, 5);%color_statistic(i, 1)Number of stored pixels,
%color_statistic(i, 2),color_statistic(i, 3),color_statistic(i, 4)Store three components of color respectively, color_statistic(i, 5)Store information on whether it is used as a supervisory clustering center
for rr = 1 : nrow
for cc = 1 : ncol
%Segment each color component
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
%Update statistics
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;%Determine the number of supervised clustering centers
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;%Has been used as a supervised clustering center
center_ACA(i, 5)=1;%center_ACA(i, 4)Used to store the number of pixels of this class, initialized to 0,center_ACA(i, 5)Used to store the pheromone concentration of this class
break
end
end
end
ant_info = zeros (num, 5);
%ant_info((r-1)*ncol+c, 1) The distance stored to this point; ant_info((r-1)*ncol+c, 2) Store the pheromone of the point; ant_info((r-1)*ncol+c,3) Heuristic function stored at this point
%ant_info((r-1)*ncol+c, 4) The probability of storing to this point; ant_info((r-1)*ncol+c, 5)Storage category
%main program
sum_ph=0;%Denominator of probability formula
do=0;%Judge whether the main program should continue to cycle
do2=1;%Judge whether the class merging program should continue to cycle
judge=zeros(nrow, ncol);%A matrix that determines whether this class can be merged with other classes
n=0;%judge ACA Number of main program cycles
m=0;%Determine the number of cycles of class merging program
while (do<=500)
%Pixel clustering
for rr = 1 : nrow
for cc = 1 : ncol
if ant_matrix(rr, cc, 1)==0
ant_matrix(rr, cc, 1) = 2;%Change the status of unclassified ants in the previous cycle to be clustered
end
end
end
for rr = 1 : nrow
for cc = 1 : ncol
%Calculate pixel(rr, cc)Distance to each cluster center, pheromone and other information
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);%When the guaranteed distance is 0, the heuristic function is large but not infinite.
sum_ph = sum_ph + ant_info(i, 2).^alpha + ant_info(i, 3).^beta;
ant_info(i, 5) = i;
end
%Calculate the probability to each cluster center
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);
%Probability selection calculation of path
if ant_info(temp, 4) >= lumda
ant_matrix(rr, cc, 1) = 1;%The pixel has been classified
cluster_img(rr, cc, 4) = temp;%Record the category of the pixel
center_ACA(temp, 4) = center_ACA(temp, 4) + 1;%The number of pixels included in the cluster is increased by 1
if ant_info(temp, 1) <= radius
center_ACA(temp, 5) = center_ACA(temp, 5) + 1;%If the distance is less than radius，Then pheromone plus 1
end
else
ant_matrix(rr, cc, 1) = 0;%The pixel is not classified and the status changes to 0
end
end
end
end
%Update cluster center information
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
%Interclass merge
%Initialize judgment matrix
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);%Record the distance between classes
temp=[d; 0];%The first record is the last distance value and the second record category
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);
%Calculate the closest class
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);
%In the pixel classification matrix,(rr, cc)The pixels of the points are classified
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
%Pheromone volatilization
for i = 1 : num
center_ACA(i, 5) = center_ACA(i, 5) * rho;
end
do = do + 1 %Once per cycle 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 main program
C = num;%The number of classes is C
m = 2;%Class parameter setting
e = nrow * ncol * C * (0.01).^2;
sum_d = e+1;
%Initialization class matrix
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
%Initialize distance matrix
distance=zeros(C, 1);
%utilize ACA Initialization membership matrix of running results
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;%If the distance between a pixel and the cluster center is 0, its membership degree is 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);%Membership normalization
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');
Complete code or write on behalf of QQ1575304183