MAXIT=50; 
I=imread('D:\rsphoto\experiment.bmp');
I1=I(250:279,130:159);
imshow(I),figure,imshow(I1)
NC=30*30; 
n1=zeros(1,NC);
n=zeros(1,NC);
for i=1:NC
    a=ceil(i/30);
    b=mod(i,30);
    if b==0
        b=30;
    end
    n1(1,i)=I1(a,b);
    n(1,i)=double(n1(1,i));
end 
n2=zeros(1,NC);
tao=zeros(NC,NC);
rho=0.8; 
alpha=3; 
beta=1; 
Q=100; 
p0=0.1;
p=zeros(NC,NC);
r=15;
mant=NC; 
qifaxinxi=zeros(NC,NC);
fenmu=0;
vsum=0;
i=1;
s=1;
belong=zeros(1,NC);
v=n;
v1=zeros(1,NC*mant*MAXIT);
d=0;
v2=zeros(1,NC);
c=0;
for l=1:MAXIT
   for i=1:NC
    for k=1:mant
         if abs(n(k)-v(i))<=r
              tao(k,i)=1;
         else tao(k,i)=0;
         end
         j=1;
         for m=1;NC;
            if (abs(v(i)-n(m)))<=r
                belong(j)=n(m);
                j=j+1;
             end
          end
         l
          for m=1:j
              qifaxinxi(k,i)=1/(abs(n(k)-v(i))+1);
              fenmu=fenmu+((tao(m,i))^alpha)*((qifaxinxi(m,i))^beta);
          end
          p(k,i)=((tao(k,i))^alpha)*((qifaxinxi(k,i))^beta)/fenmu;
          fenmu=0;
          if p(k,i)>p0
              for m=1:j
                  vsum=vsum+belong(m);
               end
               v(i)=vsum/j;
               vsum=0;
               break;
           else tao(k,i)=rho*tao(k,i)+Q*qifaxinxi(k,i);
           end
       end
   end
end
v1 = flipud(sort(v));
v2(1)=v1(1);
for m=2:NC
    if v1(m)~=v1(m-1)
        c=c+1;
        v2(c)=v1(m);
    end
end
j=1;
fenlei=zeros(c,NC);
oushijuli=zeros(NC,c);
for m=1:NC
    for i=1:c
        oushijuli(m,i)=abs(n(m)-v2(i));
    end
end
for i=1:NC
    [x,y]=min(oushijuli(i,:));
    fenlei(y,j)=i;
    j=j+1;
end 
for i=1:NC
    [x,y]=min(oushijuli(i,:));
    fenlei(y,i)=i;
    if fenlei(1,i)~=0
        n2(1,i)=v(1);
    end
    if fenlei(2,i)~=0
        n2(2,i)=v(2);
    end
    if fenlei(3,i)~=0
         n2(3,i)=v(3);
     end
end
for i=1:NC
    a=ceil(i/30);
    b=mod(i,30);
    if b==0
        b=30;
    end
     if n2(1,i)~=0
        I2(a,b,1)=1;
        I2(a,b,2)=0;
        I2(a,b,3)=0;
     end
     if n2(2,i)~=0
        I2(a,b,1)=1;
        I2(a,b,2)=1;
        I2(a,b,3)=0;
     end
     if n2(3,i)~=0
        I2(a,b,1)=0;
        I2(a,b,2)=0;
        I2(a,b,3)=1;
     end
end
figure,imshow(I2)