function Edge4=Edge16IndexSimple(H, W, Tau);
%%%% 2D
N = H*W;
Edge4=zeros(N,(2*Tau+1)*(2*Tau+1));
Diff=zeros(1, (2*Tau+1)*(2*Tau+1));
in=1;
for i=-Tau:Tau,
    for j=-Tau:Tau,
        Diff(in)=i*H+j; 
        DiffCoord(:,in)=[i;j];
        in=in+1;
    end
end
Diff(find(Diff==0))=[];
DiffCoord(:,find(Diff==0))=[];
T=length(Diff);
DiffM=repmat(Diff, [N, 1]);
iss=[1:N;]';
ISM=repmat(iss, [1, T]);
ISM=ISM+DiffM;

ISMF=ISM;
for i=1:T,
    is0=ISM(:,i);    
    index0=find(is0<1 | is0>N);
    is1=ISM(index0,T-i+1);
    index1=find(is1>=1 | is1<=N);
    is0(index0(index1))=is1(index1);
    
    index00=find(is0<1 | is0>N);
    if ~isempty(index00)
        is0(index00)=iss(index00);
    end
    
    ISMF(:,i)=is0;
end
Edge4=[iss,ISMF];
return