边缘检测算法_【图像处理】基于蚁群算法的图像边缘检测

算法如下：

设置“信息素”的初始值

为常数，表示“信息素”；

是前期遗传算法获得的最优解。

具体步骤：

1.初始值构造

每一个个体有4个参数，每一个参数有16位二进制编码。所以初始共有64位编码，其中：

1——16 ：表示状态转移概率的“信息素”影响因素α的值；

17——32：表示启发式方向函数的影响因素β的值；(启发式策略如果感兴趣的你自己可以深入去了解)

33——48：表示“信息素”发散系数λ的值；

49——64：表示迭代次数N。

2.目标函数的设计

随机把m个人工蚂蚁放在n个像素中，评估通过每个个体和目标函数产生的边缘提取的图像。

目标函数：Objv = std2(u)

其中：u为迭代的边缘图像和原始图像的能量差，即u=I1-I2。

3.访问下一个像素值的选择

Gｋ是节点的集合，没访问之前的；

Sｋ是访问下一步节点的集合，根据概率转移公式得：

Cij为状态函数，从蚁群的转移公式获得；q0位常数，范围为0~1；q为随机数，q=(0，1]，if q<=q0，且选择下一个最小像素，若q>q0，则用下面公式：

其中：

α为“信息素”启发式因子，β为启发式方向函数的影响因子，pij(t)是蚂蚁k从像素i到j在时间t的概率。

4.更新每个像素的“信息素”浓度

设置“信息素”范围为[δmin，δmax]，防止其陷入局部最小值。

(i)全局“信息素”更新；

(ii)局部“信息素”更新。

(5)如果满足终止条件，则退出算法，否则进入(3)。

function edge_ACO%参考文献："An Ant Colony Optimization Algorithm For Image Edgeclose all; clear all; clc;% 读入图像img=imread('rice.png');img = double(img)./255;[nrow, ncol] = size(img);%公式(3.24.4)初始化 for nMethod = 1:4;  %四种不同的核函数, 参见式 (3.24.7)-(3.24.10)  %E: exponential; F: flat; G: gaussian; S:Sine; T:Turkey; W:Wave  fprintf('Welcome to demo program of image edge detection using ant colony.\nPlease wait......\n');    v = zeros(size(img));    v_norm = 0;    for rr =1:nrow        for cc=1:ncol            %定义像素团            temp1 = [rr-2 cc-1; rr-2 cc+1; rr-1 cc-2; rr-1 cc-1; rr-1 cc; rr-1 cc+1; rr-1 cc+2; rr cc-1];            temp2 = [rr+2 cc+1; rr+2 cc-1; rr+1 cc+2; rr+1 cc+1; rr+1 cc; rr+1 cc-1; rr+1 cc-2; rr cc+1];            temp0 = find(temp1(:,1)>=1 & temp1(:,1)<=nrow & temp1(:,2)>=1 & temp1(:,2)<=ncol & temp2(:,1)>=1 & temp2(:,1)<=nrow & temp2(:,2)>=1 & temp2(:,2)<=ncol);            temp11 = temp1(temp0, :);            temp22 = temp2(temp0, :);            temp00 = zeros(size(temp11,1));            for kk = 1:size(temp11,1)                temp00(kk) = abs(img(temp11(kk,1), temp11(kk,2))-img(temp22(kk,1), temp22(kk,2)));            end            if size(temp11,1) == 0                v(rr, cc) = 0;                v_norm = v_norm + v(rr, cc);            else                lambda = 10;                switch nMethod                    case 1%'F'                        temp00 = lambda .* temp00;                            case 2%'Q'                        temp00 = lambda .* temp00.^2;                           case 3%'S'                        temp00 = sin(pi .* temp00./2./lambda);                    case 4%'W'                    temp00 = sin(pi.*temp00./lambda).*pi.*temp00./lambda;                end                v(rr, cc) = sum(sum(temp00.^2));                v_norm = v_norm + v(rr, cc);            end        end    end % 归一化v = v./v_norm;      v = v.*100;    p = 0.0001 .* ones(size(img));     % 信息素函数初始化    %参数设置。alpha = 1;      %式(3.24.4)中的参数beta = 0.1;     %式(3.24.4)中的参数rho = 0.1;      %式(3.24.11)中的参数%式(3.24.12)中的参数    phi = 0.05;     %equation (12), i.e., (9) in IEEE-CIM-06    ant_total_num = round(sqrt(nrow*ncol)); % 记录蚂蚁的位置    ant_pos_idx = zeros(ant_total_num, 2);  % 初始化蚂蚁的位置    rand('state', sum(clock));    temp = rand(ant_total_num, 2);    ant_pos_idx(:,1) = round(1 + (nrow-1) * temp(:,1)); %行坐标   ant_pos_idx(:,2) = round(1 + (ncol-1) * temp(:,2)); %列坐标   search_clique_mode = '8';   %Figure 1   % 定义存储空间容量   if nrow*ncol == 128*128        A = 40;        memory_length = round(rand(1).*(1.15*A-0.85*A)+0.85*A);    elseif nrow*ncol == 256*256        A = 30;        memory_length = round(rand(1).*(1.15*A-0.85*A)+0.85*A);    elseif nrow*ncol == 512*512        A = 20;        memory_length = round(rand(1).*(1.15*A-0.85*A)+0.85*A);        end    ant_memory = zeros(ant_total_num, memory_length);   % 实施算法    if nrow*ncol == 128*128        % 迭代的次数total_step_num = 300;     elseif nrow*ncol == 256*256        total_step_num = 900;     elseif nrow*ncol == 512*512        total_step_num = 1500;     end     total_iteration_num = 3;   for iteration_idx = 1: total_iteration_num        delta_p = zeros(nrow, ncol);        for step_idx = 1: total_step_num           delta_p_current = zeros(nrow, ncol);             for ant_idx = 1:ant_total_num                ant_current_row_idx = ant_pos_idx(ant_idx,1);                ant_current_col_idx = ant_pos_idx(ant_idx,2);                % 找出当前位置的邻域                if search_clique_mode == '4'                    rr = ant_current_row_idx;                    cc = ant_current_col_idx;               ant_search_range_temp = [rr-1 cc; rr cc+1; rr+1 cc; rr cc-1];                elseif search_clique_mode == '8'                    rr = ant_current_row_idx;                    cc = ant_current_col_idx;                 ant_search_range_temp = [rr-1 cc-1; rr-1 cc; rr-1 cc+1; rr cc-1; rr cc+1; rr+1 cc-1; rr+1 cc; rr+1 cc+1];                end                 %移除图像外的位置                temp = find(ant_search_range_temp(:,1)>=1 & ant_search_range_temp(:,1)<=nrow & ant_search_range_temp(:,2)>=1 & ant_search_range_temp(:,2)<=ncol);                ant_search_range = ant_search_range_temp(temp, :);                 %计算概率转换函数                ant_transit_prob_v = zeros(size(ant_search_range,1),1);                ant_transit_prob_p = zeros(size(ant_search_range,1),1);                for kk = 1:size(ant_search_range,1)         temp = (ant_search_range(kk,1)-1)*ncol + ant_search_range(kk,2);                     if length(find(ant_memory(ant_idx,:)==temp))==0                         ant_transit_prob_v(kk) = v(ant_search_range(kk,1), ant_search_range(kk,2));                        ant_transit_prob_p(kk) = p(ant_search_range(kk,1), ant_search_range(kk,2));                    else                           ant_transit_prob_v(kk) = 0;                        ant_transit_prob_p(kk) = 0;                                        end                end                 if (sum(sum(ant_transit_prob_v))==0) | (sum(sum(ant_transit_prob_p))==0)                                    for kk = 1:size(ant_search_range,1)                        temp = (ant_search_range(kk,1)-1)*ncol + ant_search_range(kk,2);                        ant_transit_prob_v(kk) = v(ant_search_range(kk,1), ant_search_range(kk,2));                        ant_transit_prob_p(kk) = p(ant_search_range(kk,1), ant_search_range(kk,2));                    end                end                                         ant_transit_prob = (ant_transit_prob_v.^alpha) .* (ant_transit_prob_p.^beta) ./ (sum(sum((ant_transit_prob_v.^alpha) .* (ant_transit_prob_p.^beta))));                      % 产生一个随机数来确定下一个位置                rand('state', sum(100*clock));                     temp = find(cumsum(ant_transit_prob)>=rand(1), 1);                ant_next_row_idx = ant_search_range(temp,1);                ant_next_col_idx = ant_search_range(temp,2);               if length(ant_next_row_idx) == 0                    ant_next_row_idx = ant_current_row_idx;                    ant_next_col_idx = ant_current_col_idx;                end                ant_pos_idx(ant_idx,1) = ant_next_row_idx;                ant_pos_idx(ant_idx,2) = ant_next_col_idx;    delta_p_current(ant_pos_idx(ant_idx,1), ant_pos_idx(ant_idx,2)) = 1;                 if step_idx <= memory_length                   ant_memory(ant_idx,step_idx) = (ant_pos_idx(ant_idx,1)-1)*ncol + ant_pos_idx(ant_idx,2);                 elseif step_idx > memory_length                    ant_memory(ant_idx,:) = circshift(ant_memory(ant_idx,:),[0 -1]);                    ant_memory(ant_idx,end) = (ant_pos_idx(ant_idx,1)-1)*ncol + ant_pos_idx(ant_idx,2);                 end                 %更新信息素函数 p = ((1-rho).*p + rho.*delta_p_current.*v).*delta_p_current + p.*(abs(1-delta_p_current));             end             delta_p = (delta_p + (delta_p_current>0))>0;             p = (1-phi).*p;          end end     % 产生边缘图矩阵，运用信息素函数判断是否是边缘    % 调用子函数进行二值分割 T = func_seperate_two_class(p);     fprintf('Done!\n');    imwrite(uint8(abs((p>=T).*255-255)), gray(256), [filename '_edge_aco_' num2str(nMethod) '.jpg'], 'jpg');end %%%%%%%%子函数%%%%%%%function level = func_seperate_two_class(I)% 功能：进行二值分割I = I(:);% STEP 1: 通过直方图计算灰度均值, 设定T=mean(I)[counts, N]=hist(I,256);i=1;mu=cumsum(counts);T(i)=(sum(N.*counts))/mu(end);% STEP 2: 计算灰度值大于 T (MAT)像素的均值和灰度值小于 T(MBT)像素的均值% step 1mu2=cumsum(counts(N<=T(i)));MBT=sum(N(N<=T(i)).*counts(N<=T(i)))/mu2(end);mu3=cumsum(counts(N>T(i)));MAT=sum(N(N>T(i)).*counts(N>T(i)))/mu3(end);i=i+1;T(i)=(MAT+MBT)/2;% STEP 3 ：当T(i)~=T(i-1)，重复STEP 2Threshold=T(i);while abs(T(i)-T(i-1))>=1    mu2=cumsum(counts(N<=T(i)));    MBT=sum(N(N<=T(i)).*counts(N<=T(i)))/mu2(end);    mu3=cumsum(counts(N>T(i)));    MAT=sum(N(N>T(i)).*counts(N>T(i)))/mu3(end);    i=i+1;    T(i)=(MAT+MBT)/2;     Threshold=T(i);end% 归一化均值到 [i, 1]之间.level = Threshold;