遗传学算法概述
从之前转载的博客《非常好的理解遗传算法的例子》中可以知道,遗传学算法主要有6个步骤:
1. 个体编码
2. 初始群体
3. 适应度计算
4. 选择运算
5. 交叉运算
6. 变异运算
这是一个仿生的过程,模仿生物进化和自然选择。在该算法中,个体编码就相当于生物最基本的组成——基因,初始群体就是刚开始那些个原始的生物体。
在恶劣的环境中,适者生存的自然法则将让适应能力更好的生物继续存活繁衍下去,而适应能力差的生物将会被淘汰。因此遗传算法通过计算适应度来模拟这个自然选择的过程,用于筛选优良的个体保留下来。光这么讲也许太抽象,可能读者也无法知道具体如何实现的,也不知道这个适应度是如何去评价去个体的优劣程度的。下面举例个简单的例子来说明,如下图:假如现在要求一条曲线方程的最大值,并且已知有若干个点。那么曲线的值也就是求最大值过程中的适应度,因为值越大就越接近最大值,其适应度就越好。因此在这个求最大值的过程中,曲线对应的值可以直接作为适应度。在下一步选择运算过程中,这些适应能力好的个体也有更高的概率被保留下来。
选择运算就是自然选择的过程。在适应度计算过程中,各个个体的适应能力已经知道,因此部分适应能力差的个体就将会被淘汰。就想自然选择一样,大自然不会对生物体适应能力进行排序,然后一刀切全部淘汰弱者。在选择运算中,算法也是根据适应度来确定一个选择概率,按照概率来确定哪些个体允许保留下来。自然地,适应度高的,被保留下来的概率也更高,而适应度低的保留下来的可能性就越低。
进行选择之后,保留下来的个体就像生物一样,需要在自己的生命周期内繁衍后代来延续自己的物种。然后,就是脸红心跳的配对阶段了,两个个体之间需要贡献各自的染色体进行交换。这一过程中新的个体经过交叉运算获得由父母双方部分基因的组合,新的个体中基因来自于父母双方。当然了,算法中交差运算没有男女之分,交叉运算后所有的个体都是交叉运算产生的新个体,原本的个体被新个体代替不进入下一次迭代。如《非常好的理解遗传算法的例子》中过程5的描述:
交叉的时候进行交换基因的位置也是随机产生的,交换个体同样按照各自适应度以一定概率进行选择。因此适应度高的个体进行交差运算的概率也更高,这个在动物界比较普遍,即强壮的首领拥有交配权,以便优良基因更多机会保留下去。并且交叉本身也有一定发生概率,如果不进行交叉,那么交叉的结果就是新个体与原个体相同。这个理解起来也方便,就是说,要是进行交叉了,新个体为交叉结果;那要是不孕育新个体呢,那抱歉,你们两个个体还要经过下一次迭代。
最后,就是生物进化的过程了,这是一个小概率事件,但却是造就如今这样美丽而伟大的生命世界的关键。变异是随机的,因为随机,因此绝大部分变异都是恶性的,但也正是这种不确定性才导致了生物的多样性,让我们的世界变得美妙多彩。因为有自然选择的过程,变异后恶化的基因会被抛弃,因此,整个生物界依靠这几个关键步骤保持着不断进化,不断适应的趋势。但是,同样地,我们可以预料到,在变异运算这一步中,变异出现的概率会非常小,否则整个算法也就随机了。说了那么多,归结起来,变异在算法中的体现说白了就是按照一个小概率去挑选一个个体,然后将其某个编码后的码字取个反,当然,这是在二进制编码的情况下来说啦。
本身遗传学算法也是一门很复杂的学问,研究到现在,求解过程越来越像一个生命体的进化过程,但是基本思想还是如上述所说。同样地,下面我也将上述最基础的遗传学算法用matlab进行实现,各位读者也可以参照着代码来理解算法。
Matlab算法实现
下面演示,利用遗传学算法来求解最佳分割阈值。这个方法参考论文《Infrared Image Segmentation Based on Otsu and Genetic Algorithm》,该论文基于大津法和遗传算法来分割红外图像。不说算法有多好,其实也差不多,还不一定有Otsu来得直接(抱歉,我喜欢喷一喷(*^__^*) ),但是通过这个案例能够有效地描述遗传算法的使用。这里要感谢一下我师兄将该文献翻译的成中文,让我看起来方便多了。
大津法大家应该都听说过,目的是求最优的前景和背景的分界面,在图像分割算法中,大津法也是最有名的算法,个人觉得,没有之一。这里实际是使用遗传算法来求这个最优的分界面,就如上述我举例说求曲线的最大值一样,只是这里面最适应度的计算不是那么直接。下面从总体到各部分细节一一介绍。
(一)分割算法
这里将文中的算法命名为OtsuWithGenetic,算法不需要人为添加任何参数来辅助运算。从算法的名称来看,就可以知道,算法结合了大津法和遗传算法两者的特征。因此,在进行遗传算法的迭代流程之前,需要计算图像的直方图。而进行遗传学迭代运算的过程,在本案例中,我将其封装成一个名为BestThresh的函数,它的输入同样仅需要图像的直方图就可以了。在求取最佳阈值之后,就对图像进行二值化,得到最终结果图。
所以分割算法整体分为三个步骤:
1. 计算图像直方图
2. 利用遗传学算法求最佳阈值
3. 图像分割,进行二值化
function ImOut = OtsuWithGenetic(ImOriginal)
% 大津分割算法,结合遗传算法
% 在性能方面暂未进行相应优化,只对求解方程进行了优化
% 部分代码结构比较凌乱 ImGray = double(rgb2gray(ImOriginal)); %转换为双精度灰度图 [m,n] = size(ImGray); %计算图片的尺寸
ImHist = zeros(1,256); %灰度直方图
ImOut = zeros(m,n); %分割结果图像 % 统计各灰度级的像素点个数,计算灰度直方图
for i=1:m
for j=1:n
ImHist(ImGray(i,j)+1) = ImHist(ImGray(i,j)+1)+1;
end
end th = BestThresh(ImHist); % 分割图像,高于阈值为前景255,低于阈值为背景0
for i = 1:m
for j = 1:n
if (ImGray(i,j)+1>th)
ImOut(i,j) = 255;
else
ImOut(i,j) = 0;
end
end
end ImOut=uint8(ImOut); %读入读出变换
end
(二)遗传算法求解大津最佳阈值
这一步实际上就是BestThresh函数的描述。从整个函数的代码中我们可以很清楚地看到遗传学算法描述中的6个步骤。
function bestindividual = BestThresh(ImHist)
popsize=40; %群体大小
chromlength=8; %字符串长度(个体长度)
pc=0.9; %交叉概率
pm=0.001; %变异概率 pop = InitEntity(popsize,chromlength); %随机产生初始群体
for i = 1:50 %50为迭代次数
[objvalue] = FitnessValue(pop, ImHist); %计算目标函数
fitvalue = CalFitValue(objvalue); %计算群体中每个个体的适应度
[newpop] = Selection(pop,fitvalue); %选择
[newpop] = CrossOver(newpop,pc); %交叉
[newpop] = Mutation(newpop,pm); %变异
pop = newpop;
end
[bestindividual,bestfit] = BestEntity(pop,fitvalue); %求出群体中适应值最大的个体及其适应值
bestindividual = DecodeChrom(bestindividual,1,chromlength);
end
这里在一开始定义了四个参数:popsize, chromlength, pc, pm。具体含义见注释,这里重点介绍pc和pm两个参数。pc是交叉概率,pm是变异概率,在上述描述中我们说到,交叉和变异都是按照一定概率进行的,从这里也可以看到,两者的概率之间的差别。这里我默认将两者的概率定为了0.9和0.001,同样,迭代次数和群体的大小也定死了,实际上,运算次数再多,群体再大,也是浪费运算资源。
从代码中看到,在初始化完个体后,就按顺序进入迭代运算,迭代经过适应度的计算、选择运算、交叉运算和变异运算。其中计算目标函数和个体适应度计算实际是一个内容,CalFitValue实质是进行了矫正防止出现负数,而FitnessValue是适应度计算函数。
迭代结束以后,将适应度最大的作为最终的阈值,并且经过解码之后转换为实际的灰度值。
1. 初始化群体+个体编码
在本案例中,需要计算最佳的阈值,也就是计算作为阈值的那个最优的灰度值。对灰度值进行二进制编码获得随机的40个群体代码如下:
% 1 初始化群体(编码)
% popsize表示群体的大小,chromlength表示染色体的长度(二值数的长度)
function entities = InitEntity(entitySize,chromLength)
% rand随机产生随机数矩阵,round将其二值化
entities = round(rand(entitySize, chromLength));
end
该函数生成一个40行的矩阵,代码40个不同的初始个体;每个个体进行8位二进制编码,因此像素灰度本身也仅有8位深度。
2. 计算适应度
适应度的计算参考论文当中适应度函数
这里这里,是适应度函数,是目标图像像素数目,是背景图像像素数目;是目标图像像素总和,是背景图像像素总和。
% 2 计算适应函数
% 实现目标函数的计算
function [objvalue] = FitnessValue(pop, Imhist)
graysum = Imhist.*(0:1:255); % 各灰度值像素总和
x = DecodeChrom(pop,1,8);% 将pop每行转化成十进制数
loop = size(pop,1);
objvalue = zeros(loop,1);
for i = 1 : loop
P1 = sum(graysum(x(i)+1 : 255));
P2 = sum(graysum(1 : uint8(x(i))));
I = sum(Imhist(x(i)+1 : 255));
J = sum(Imhist(1 : uint8(x(i))));
objvalue(i) = I * J * ((P1-P2)^2) * 2 ;%计算目标函数值
end
end
function fitValue = CalFitValue(objValue)
Cmin = 0;
num = size(objValue,1);
fitValue = zeros(1,num);
for i = 1 : num
if objValue(i) + Cmin > 0
temp = Cmin + objValue(i);
else
temp = 0.0;
end
fitValue(i) = temp;
end
end
3. 选择运算
通过计算得到的适应度就可以知道各个体的适应能力,将这些适应度进行累加,各个体就对应于一个区间。将累加结果进行归一化,然后用随机函数来进行选择,落入哪个区间就选择区间对应的个体。
% 3 选择复制
function [selectedEntities] = Selection(entities,fitValue)
totalFit = sum(fitValue); %求适应值之和
fitValue = fitValue / totalFit; %单个个体被选择的概率
fitValue = cumsum(fitValue); %如 fitValue=[1 2 3 4],则 cumsum(fitValue)=[1 3 6 10]
[num,cols] = size(entities);
ms = rand(num,1); %产生用于选择的随机数 selectedEntities = zeros(num,cols); %初始化选择后的个体大小
for i = 1 : num
for j = 1 : num
if ms(i) < fitValue(j)
selectedEntities(i,:) = entities(j,:);
break;
end
end
end
end
4. 交叉运算
交叉运算利用一对个体产生新个体,函数如下:
% 4 交叉
% 群体中两个个体进行一定的概率 pc 的交叉,交换基因起始位置cpoint随机
function newEntities = CrossOver(entities,pc)
[num, cols]= size(entities);
newEntities = ones(size(entities));
for i = 1 : 2 : num-1
if(rand < pc)
cpoint = round(rand*cols); % 随机产生交换起始位置
newEntities(i,:) = [entities(i,1:cpoint), entities(i+1,cpoint+1:cols)];
newEntities(i+1,:) = [entities(i+1,1:cpoint), entities(i,cpoint+1:cols)];
else
newEntities(i,:) = entities(i,:);
newEntities(i+1,:) = entities(i+1,:);
end
end
end
5. 变异运算
在经过交叉运算以后,还要进行编译运算,由于是二进制编码,变异运算就可以直接将变异位置的二进制数取反。
% 5 变异
function [newEntities] = Mutation(entities,pm)
[num, cols] = size(entities);
newEntities = ones(size(entities));
for i = 1 : num
if(rand < pm)
mpoint = round(rand*cols); % 计算变异位置
if mpoint <= 0
mpoint = 1;
end
newEntities(i) = entities(i);
if any(newEntities(i,mpoint)) == 0 % 进行二进制取反
newEntities(i,mpoint) = 1;
else
newEntities(i,mpoint) = 0;
end
else
newEntities(i,:) = entities(i,:);
end
end
end
6. 计算最优值
在迭代结束后,所有的个体基本上呈现向最优值靠拢的趋势,这时候需要通过比较来取其中的最大值,作为最优解。这一步是一个比较过程,原理很简单:
% 6 求群体中对应的最大适应值
function [bestEntity,bestFit] = BestEntity(entities, fitValue)
num = size(entities, 1);
bestEntity = entities(1,:);
bestFit = fitValue(1);
for i = 2 : num
if fitValue(i) > bestFit
bestEntity = entities(i,:);
bestFit = fitValue(i);
end
end
end
(三)二进制转换函数
在运算中,由于使用数组来表示二进制编码数据,因此在运算过程中经常需要进行二进制与十进制之间的转换,转换函数如下:
% 将二进制编码转换成十进制
function decimals = DecodeChrom(entities,spoint,length)
temp = entities(:, spoint:spoint+length-1);
decimals = DecodeBinary(temp);
end % 将二进制转化为十进制
function decimals = DecodeBinary(pop)
[num,cols] = size(pop);
temp = zeros(num,cols);
for i = 1 : cols
temp(:,i) = 2.^(cols-i) .* pop(:,i);
end
decimals = sum(temp,2);
end