function [D] = process (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);
[x y]=find(g==max(max(g)));
g=scroll2(g,x+1,y+1);
g=g/sum(sum(g));
g=fft2(g);

template = rand(size)<0.2;
template = arrange(template,g);
D=0.2*ones(size);
d=0.2-1/(size^2);
temp=template;

while d>0
    filt=real(ifft2(g.*fft2(temp)));
    [x,i]=max(filt.*temp,[],1);
    [x,j]=max(x,[],2);i=i(j);
    D(i,j)=d;
    temp(i,j)=0;
    d=d-1/(size^2);
end
temp=template;
d=0.2+1/(size^2);
while d<1
    filt=real(ifft2(g.*fft2(temp)));
    [x,i]=min(filt+temp,[],1);
    [x,j]=min(x,[],2);i=i(j);
    D(i,j)=d;
    temp(i,j)=1;
    d=d+1/(size^2);
end


function mat = arrange(mat,g)
[b1,b2]=size(mat);
a=sum(sum(mat));
filt=real(ifft2(g.*fft2(mat)));
[x,i]=max(filt.*mat,[],1);
[x,j]=max(x,[],2);i=i(j);
[x,k]=min(filt+x*mat,[],1);;
[x,l]=min(x,[],2);k=k(l);
m=mat(i,j);
mat(i,j)=mat(k,l);
mat(k,l)=m;
fin=0;
iterations=1;
while (~fin)
    last=[i j];
    filt=abs(ifft2(g.*fft2(mat)));
    [x,i]=max(filt.*mat,[],1);
    [x,j]=max(x,[],2);i=i(j);
    [x,k]=min(filt+x*mat,[],1);
    [x,l]=min(x,[],2);k=k(l);
    mat(i,j)=m;
    mat(i,j)=mat(k,l);
    mat(k,l)=m;
    if sum((last-[k l]).^2)<1
        fin=1;
    end
    iterations=iterations+1
end
iterations