function Dj = interpj(Fp,Ip,Dj0,H9);
# This function interpolates using jinc function
# as a sparse matrix using Dj0 as pattern
#
# Usage: Dj = interpj(Fp,Ip,Dj0,H9)
# Fp: source points
# Ip: destination points
# Dj0: pattern matrix
# H9: adjacency matrix for Fp

# the magic constant is chosen that adjacent points are at 1st besselj zero
[i,j,v] = find(H9);
meandist = accumarray(i,abs(Fp(j)-Fp(i)),[],@mean);
magic91=3.8317; # 1st zero of Bessel J1
[i,j,v] = find(Dj0);
Dj = sparse(i,j,jinc(magic91*abs(Fp(i)-Ip(j))./meandist(i)),size(Dj0)(1),size(Dj0)(2));
Dj *= diag(1./sum(Dj));
endfunction
