psolbg2.m

将PSO和LBG结合在一步迭代过程中

% =====================================================================
%                 PP algorithm 
%                 sub program for pp
%     Start Date: 15/07/2006
%   Last Changed: 21/08/2006
%  
%  Specification: The algorithm, proposed by Ji and Liao in 2006, incorporating
%                 both PSO algorithm and LBG to one iteration and using  
%                 particle-pair to explore in the problem space
%
% Copyright (c) Ji and Liao in 2006
% All rights Reserved
%==========================================================================

 function [gbest,psnr1,t2,gbest_value,tempz_cluster]=PSOLBG2(particle,pbest_value)

%Initialization of PSO parameters
itmax=15; %Maximum iteration number
c1=0.3;
c2=0.5;
GL=255;

for iter=1:itmax
W(iter)=0.1;
end


%**********************************************************
IndexIteration=1;
D=16;
num_cluster=256;
GenTrainingVector;
tic;
%**********************************************************

%**********************************************************
a=0;  
b=255;
num_particle=2; % number of particle, should be more than
m=0;
n=255;
V=round(m+(n-m)*rand(num_particle,num_cluster*D,1));% produce velocity randomly for all particles and all itmax

%**********************************************************
fitness=zeros(num_particle,itmax);% each particle's fitness value for all iterations
pbest=zeros(num_particle,num_cluster*D,itmax); % each particle's position for all iteration
%**********************************************************
fitness(:,1)=pbest_value;
pbest_value(:,1)=fitness(:,1);
[C,I]=min(abs(pbest_value(:,1)));
gbest_value(1)=C;
gbest(1,:,1)=particle(I,:,1);  % save the best particle 's position (in other words, the best codebook)
for i=1:num_particle; % for the 1st iteration, each particle's best position is themselves
       pbest(i,:,1)=particle(i,:,1);
end
% calculate psnr
Dis=gbest_value(1)/4096;
psnr = 10 * log10(255.^2*16 /Dis);

%******************************************************
while IndexIteration<=(itmax-1)
   IndexIteration=IndexIteration+1;

V(:,:,IndexIteration)=floor(W(IndexIteration)*V(:,:,IndexIteration-1)+c1*rand*(pbest(:,:,IndexIteration-1)...
    -particle(:,:,IndexIteration-1))+c2*rand*(repmat(gbest(:,:,IndexIteration-1),num_particle,1)-particle(:,:,IndexIteration-1)));
particle(:,:,IndexIteration)=particle(:,:,IndexIteration-1)+V(:,:,IndexIteration);

particle(:,:,IndexIteration)=min(particle(:,:,IndexIteration),255);

updateparticle_lhl;
calfitness;

%update pbest and pbest_value
pbest_value(:,IndexIteration)=min(pbest_value(:,IndexIteration-1),fitness(:,IndexIteration));
for (temp=1:num_particle)
    if pbest_value(temp,IndexIteration)== fitness(temp,IndexIteration)   
        pbest(temp,:,IndexIteration)=particle(temp,:,IndexIteration);
    else
        pbest(temp,:,IndexIteration)=pbest(temp,:,IndexIteration-1);
    end
end
% update gbest and gbest_value
[C,I]=min(abs(pbest_value(:,IndexIteration)));
gbest_value(IndexIteration)=min(C,gbest_value(IndexIteration-1));

    if gbest_value(IndexIteration)==C
        gbest(1,:,IndexIteration)=pbest(I,:,IndexIteration);  % save the best particle 's position (in other words, the best codebook)
    else
        gbest(1,:,IndexIteration)=gbest(1,:,IndexIteration-1);
    end
psnr(IndexIteration)=10 * log10(255.^2*16 /(gbest_value(1,IndexIteration)/4096));
% #########################################################################

end   % end for while
t2=toc;
tempz_cluster=reshape(gbest(1,:,IndexIteration),D,num_cluster)';
% deal with one point to one codeword, but the codeword is not the same as the point 
num_trainingvector=4096;
    for j=1:num_trainingvector 
        [C(j),I(j)]=min(sum((repmat(trainingvector(j,:),num_cluster,1)-tempz_cluster).^2,2));
       
    end
Dmin=C';
index=I';
size=zeros(1,256);
for i=1:4096
    temp=index(i);
    size(temp)=size(temp)+1;
end
[numofclustersize,geicodebook]=sort(size);
ntemp=find(numofclustersize==1);
numtotest=ntemp(end);
for i=1:numtotest
    temp=geicodebook(i);
    fortv=find(index==temp);
    if (trainingvector(fortv,1)~=tempz_cluster(temp,1))
        tempz_cluster(temp,:)=trainingvector(fortv,:);
    end
end
% reconstruct the image following
for j=1:num_trainingvector 
    [C,I]=min(sum((repmat(trainingvector(j,:),num_cluster,1)-tempz_cluster).^2,2));
    new_tv(j,:)=tempz_cluster(I,:);
end
recn_img=img_lena;
for i=0:Line/q-1
   for j=0:Col/q-1
       recn_img(i*q+1,[j*q+1:j*q+q])=new_tv(i*Col/q+j+1,[1:4]);
       recn_img(i*q+2,[j*q+1:j*q+q])=new_tv(i*Col/q+j+1,[5:8]);
       recn_img(i*q+3,[j*q+1:j*q+q])=new_tv(i*Col/q+j+1,[9:12]);
       recn_img(i*q+4,[j*q+1:j*q+q])=new_tv(i*Col/q+j+1,[13:16]);
    end
end
[C,I]=find(recn_img<0);
[C1,I1]=find(recn_img>255);
GL=255;
mse=mean(mean((recn_img-img_lena).^2));
psnr;
psnr1=10*log10(GL*GL/mse);
figure;
imshow(recn_img/255);