function RD=compute_rd(X,Y)

% Assumption NxM grid
[N,M] = size(X);

% memory allocation
RD = NaN(size(X));

for row = 2:N-1
    for col=2:M-1
        nn=0;
        r1=sqrt((X(row+1,col)-X(row,col))^2+(Y(row+1,col)-Y(row,col))^2);if(~isnan(r1)); nn=nn+1; end
        r2=sqrt((X(row-1,col)-X(row,col))^2+(Y(row-1,col)-Y(row,col))^2);if(~isnan(r2)); nn=nn+1; end
        r3=sqrt((X(row,col+1)-X(row,col))^2+(Y(row,col+1)-Y(row,col))^2);if(~isnan(r3)); nn=nn+1; end
        r4=sqrt((X(row,col-1)-X(row,col))^2+(Y(row,col-1)-Y(row,col))^2);if(~isnan(r4)); nn=nn+1; end
        if isequal(nn,4)
            RD(row,col) = (r1+r2+r3+r4)/4;
        end
    end
end


end