function [D] = process_color (size)

for i=1:size
    for j=1:size
        g(i,j)=exp(-((i-1)^2+(j-1)^2)/4);
        g1(i,j)=exp(-(i^2+j^2)/4);
    end% test different values of sigma
end
g=g+fliplr(g1)+flipud(g1)+fliplr(flipud(g1));
g=gauss2(3.5,size,3.5,size);
g2=gauss2(3.5,size,3.5,size);
[x y]=find(g==max(max(g)));
g=scroll2(g,x+1,y+1);
[x y]=find(g2==max(max(g2)));
g2=scroll2(g2,x+1,y+1);
g=g/sum(sum(g));
g2=g2/sum(sum(g2));
g=fft2(g)*2;
g2=fft2(g2)/2;

m = [1.7008 -1.1765 0.7961;1.9217 -2.0651 1.3798;2.6610 -0.3566 -2.1312];
m1 = [1.7008 -1.1765 0.7961];
m1 = [0.3195 0.5669 0.1607];
minv = inv(m);
template = rand(size,size,3)<0.2;

template = arrange(template,g,g2);
D=0.2*ones(size,size,3);
d=0.2-1/(size^2);
temp=template;
b1=size;
b2=size;
while d>0
    matv = reshape((reshape(temp,b1*b2,3))*m',b1,b2,3);
    filt(:,:,1)=real(ifft2(g.*fft2(matv(:,:,1))));
    filt(:,:,2)=real(ifft2(g2.*fft2(matv(:,:,2))));
    filt(:,:,3)=real(ifft2(g2.*fft2(matv(:,:,3))));
    %filt_vis = reshape(m1*(reshape(filt,b1*b2,3))',b1,b2);
    filt_vis = reshape((reshape(filt,b1*b2,3))*minv',b1,b2,3);
    [x,i]=max(filt_vis(:,:,1).*temp(:,:,1),[],1);
    [x,j]=max(x,[],2);i=i(j);
    D(i,j,1)=d;
    temp(i,j,1)=0;
    
    matv = reshape((reshape(temp,b1*b2,3))*m',b1,b2,3);
    filt(:,:,1)=real(ifft2(g.*fft2(matv(:,:,1))));
    filt(:,:,2)=real(ifft2(g2.*fft2(matv(:,:,2))));
    filt(:,:,3)=real(ifft2(g2.*fft2(matv(:,:,3))));
    %filt_vis = reshape(m1*(reshape(filt,b1*b2,3))',b1,b2);
    filt_vis = reshape((reshape(filt,b1*b2,3))*minv',b1,b2,3);
    [x,i]=max(filt_vis(:,:,2).*temp(:,:,2),[],1);
    [x,j]=max(x,[],2);i=i(j);
    D(i,j,2)=d;
    temp(i,j,2)=0;
    
    matv = reshape((reshape(temp,b1*b2,3))*m',b1,b2,3);
    filt(:,:,1)=real(ifft2(g.*fft2(matv(:,:,1))));
    filt(:,:,2)=real(ifft2(g2.*fft2(matv(:,:,2))));
    filt(:,:,3)=real(ifft2(g2.*fft2(matv(:,:,3))));
    %filt_vis = reshape(m1*(reshape(filt,b1*b2,3))',b1,b2);
    filt_vis = reshape((reshape(filt,b1*b2,3))*minv',b1,b2,3);
    [x,i]=max(filt_vis(:,:,3).*temp(:,:,3),[],1);
    [x,j]=max(x,[],2);i=i(j);
    D(i,j,3)=d;
    temp(i,j,3)=0;
    
    d=d-1/(size^2);
end
temp=template;
d=0.2+1/(size^2);
while d<1
    matv = reshape((reshape(temp,b1*b2,3))*m',b1,b2,3);
    filt(:,:,1)=real(ifft2(g.*fft2(matv(:,:,1))));
    filt(:,:,2)=real(ifft2(g2.*fft2(matv(:,:,2))));
    filt(:,:,3)=real(ifft2(g2.*fft2(matv(:,:,3))));
    %filt_vis = reshape(m1*(reshape(filt,b1*b2,3))',b1,b2);
    filt_vis = reshape((reshape(filt,b1*b2,3))*minv',b1,b2,3);
    [x,i]=min(filt_vis(:,:,1)+256*temp(:,:,1),[],1);
    [x,j]=min(x,[],2);i=i(j);
    D(i,j,1)=d;
    temp(i,j,1)=1;
    
    matv = reshape((reshape(temp,b1*b2,3))*m',b1,b2,3);
    filt(:,:,1)=real(ifft2(g.*fft2(matv(:,:,1))));
    filt(:,:,2)=real(ifft2(g2.*fft2(matv(:,:,2))));
    filt(:,:,3)=real(ifft2(g2.*fft2(matv(:,:,3))));
    %filt_vis = reshape(m1*(reshape(filt,b1*b2,3))',b1,b2);
    filt_vis = reshape((reshape(filt,b1*b2,3))*minv',b1,b2,3);
    [x,i]=min(filt_vis(:,:,2)+256*temp(:,:,2),[],1);
    [x,j]=min(x,[],2);i=i(j);
    D(i,j,2)=d;
    temp(i,j,2)=1;
    
    matv = reshape((reshape(temp,b1*b2,3))*m',b1,b2,3);
    filt(:,:,1)=real(ifft2(g.*fft2(matv(:,:,1))));
    filt(:,:,2)=real(ifft2(g2.*fft2(matv(:,:,2))));
    filt(:,:,3)=real(ifft2(g2.*fft2(matv(:,:,3))));
    %filt_vis = reshape(m1*(reshape(filt,b1*b2,3))',b1,b2);
    filt_vis = reshape((reshape(filt,b1*b2,3))*minv',b1,b2,3);
    [x,i]=min(filt_vis(:,:,3)+256*temp(:,:,3),[],1);
    [x,j]=min(x,[],2);i=i(j);
    
    D(i,j,3)=d;
    temp(i,j,3)=1;
    
    d=d+1/(size^2);
end


function mat = arrange(mat,g,g2)
[b1,b2,b3]=size(mat);
a=sum(sum(mat));
m = [1.7008 -1.1765 0.7961;1.9217 -2.0651 1.3798;2.6610 -0.3566 -2.1312];
%m=[1 0 0;0 1 0;0 0 1];
m1 = [1.7008 -1.1765 0.7961];
m1 = [0.3195 0.5669 0.1607];
minv = inv(m);
matv = reshape((reshape(mat,b1*b2,3))*m',b1,b2,3);
filt(:,:,1)=real(ifft2(g.*fft2(matv(:,:,1))));
filt(:,:,2)=real(ifft2(g2.*fft2(matv(:,:,2))));
filt(:,:,3)=real(ifft2(g2.*fft2(matv(:,:,3))));
%filt_vis = reshape(m1*(reshape(filt,b1*b2,3))',b1,b2);
filt_vis = reshape((reshape(filt,b1*b2,3))*minv',b1,b2,3);
[x,i]=max(filt_vis(:,:,1).*mat(:,:,1),[],1);
[x,j]=max(x,[],2);i=i(j);
[x,k]=min(filt_vis(:,:,1)+x*mat(:,:,1),[],1);;
[x,l]=min(x,[],2);k=k(l);
m0=mat(i,j,1);
mat(i,j,1)=mat(k,l,1);
mat(k,l,1)=m0;
last1 = [i j];

matv = reshape((reshape(mat,b1*b2,3))*m',b1,b2,3);
filt(:,:,1)=real(ifft2(g.*fft2(matv(:,:,1))));
filt(:,:,2)=real(ifft2(g2.*fft2(matv(:,:,2))));
filt(:,:,3)=real(ifft2(g2.*fft2(matv(:,:,3))));
%filt_vis = reshape(m1*(reshape(filt,b1*b2,3))',b1,b2);
filt_vis = reshape((reshape(filt,b1*b2,3))*minv',b1,b2,3);
[x,i]=max(filt_vis(:,:,2).*mat(:,:,2),[],1);
[x,j]=max(x,[],2);i=i(j);
[x,k]=min(filt_vis(:,:,2)+x*mat(:,:,2),[],1);;
[x,l]=min(x,[],2);k=k(l);
m0=mat(i,j,2);
mat(i,j,2)=mat(k,l,2);
mat(k,l,2)=m0;
last2 = [i j];

matv = reshape((reshape(mat,b1*b2,3))*m',b1,b2,3);
filt(:,:,1)=real(ifft2(g.*fft2(matv(:,:,1))));
filt(:,:,2)=real(ifft2(g2.*fft2(matv(:,:,2))));
filt(:,:,3)=real(ifft2(g2.*fft2(matv(:,:,3))));
%filt_vis = reshape(m1*(reshape(filt,b1*b2,3))',b1,b2);
filt_vis = reshape((reshape(filt,b1*b2,3))*minv',b1,b2,3);
[x,i]=max(filt_vis(:,:,3).*mat(:,:,3),[],1);
[x,j]=max(x,[],2);i=i(j);
[x,k]=min(filt_vis(:,:,3)+x*mat(:,:,3),[],1);;
[x,l]=min(x,[],2);k=k(l);
m0=mat(i,j,3);
mat(i,j,3)=mat(k,l,3);
mat(k,l,3)=m0;
last3=[i j];
fin=0;
fin1=0;
fin2=0;
fin3=0;

iterations=1;
while (~fin)
    if iterations>100
        a;
    end
    if ~fin1    
        matv = reshape((reshape(mat,b1*b2,3))*m',b1,b2,3);
        filt(:,:,1)=real(ifft2(g.*fft2(matv(:,:,1))));
        filt(:,:,2)=real(ifft2(g2.*fft2(matv(:,:,2))));
        filt(:,:,3)=real(ifft2(g2.*fft2(matv(:,:,3))));
        %filt_vis = reshape(m1*(reshape(filt,b1*b2,3))',b1,b2);
        filt_vis = reshape((reshape(filt,b1*b2,3))*minv',b1,b2,3);
        [x,i]=max(filt_vis(:,:,1).*mat(:,:,1),[],1);
        [x,j]=max(x,[],2);i=i(j);
        [x,k]=min(filt_vis(:,:,1)+x*mat(:,:,1),[],1);;
        [x,l]=min(x,[],2);k=k(l);
        m0=mat(i,j,1);
        mat(i,j,1)=mat(k,l,1);
        mat(k,l,1)=m0;
    end
    if sum((last1-[k l]).^2)<1
        fin1=1;
    end
    last1=[i j];
    if ~fin2    
        matv = reshape((reshape(mat,b1*b2,3))*m',b1,b2,3);
        filt(:,:,1)=real(ifft2(g.*fft2(matv(:,:,1))));
        filt(:,:,2)=real(ifft2(g2.*fft2(matv(:,:,2))));
        filt(:,:,3)=real(ifft2(g2.*fft2(matv(:,:,3))));
        %filt_vis = reshape(m1*(reshape(filt,b1*b2,3))',b1,b2);
        filt_vis = reshape((reshape(filt,b1*b2,3))*minv',b1,b2,3);
        [x,i]=max(filt_vis(:,:,2).*mat(:,:,2),[],1);
        [x,j]=max(x,[],2);i=i(j);
        [x,k]=min(filt_vis(:,:,2)+x*mat(:,:,2),[],1);;
        [x,l]=min(x,[],2);k=k(l);
        m0=mat(i,j,2);
        mat(i,j,2)=mat(k,l,2);
        mat(k,l,2)=m0;
    end
    if sum((last2-[k l]).^2)<1
        fin2=1;
    end
    last2=[i j];
    if ~fin3
        matv = reshape((reshape(mat,b1*b2,3))*m',b1,b2,3);
        filt(:,:,1)=real(ifft2(g.*fft2(matv(:,:,1))));
        filt(:,:,2)=real(ifft2(g2.*fft2(matv(:,:,2))));
        filt(:,:,3)=real(ifft2(g2.*fft2(matv(:,:,3))));
        %filt_vis = reshape(m1*(reshape(filt,b1*b2,3))',b1,b2);
        filt_vis = reshape((reshape(filt,b1*b2,3))*minv',b1,b2,3);
        [x,i]=max(filt_vis(:,:,3).*mat(:,:,3),[],1);
        [x,j]=max(x,[],2);i=i(j);
        [x,k]=min(filt_vis(:,:,3)+x*mat(:,:,3),[],1);    
        [x,l]=min(x,[],2);k=k(l);
        
        m0=mat(i,j,3);
        mat(i,j,3)=mat(k,l,3);
        mat(k,l,3)=m0;
    end
    if sum((last3-[k l]).^2)<1
        fin3=1;
    end
    last3=[i j];
    
    if fin1&fin2&fin3
        fin=1;
    end
    iterations=iterations+1
end
iterations