Unified Fourier Reconstruction (gridding)

This is some earlier work of mine. I was working on an algorithm to impelment Band-limited interpolation from unevenly spaced sampled data after Mehrdad Soumekh. It is also known as “gridding” in SAR and CAT communities (see The gridding method for image reconstruction by Fourier transformation Schomberg, H. and Timmer, J., for example).

The best way to delve into the secrets of UFR is probably to have a look at some code. This file accompanied the book “Fourier Array Imaging” by M. Soumekh. The book is out of print, but you might find it in a good technical university library.

Unfortunately, one cannot place m-files here. So I’ll just drop the code in the text. The convention will be that text in bold italics will mark separate Matlab scripts or functions. Normal text is the Matlab code itself. They should all run without a hitch – you’ll just need to fix some wrapped lines, I suspect. Also, the indentations disappeared, but I don’t have time to fix that now. If there are some other problems, please let me know.

The starting code from the file above is here. You can copy and paste it into Matlab editor as ufr.m, for instance.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INTERPOLATION VIA UFR %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This matlab program performs interpolation from nonuniform data stored
% in array g(alpha,beta). In its present form, the target function is a
% rectangular pulse. The transformation is POLAR.
%
% (x0,y0): target’s support
% (nx,ny): target’s image size
% mm: the number of samples within the Hamming window that cuts
% off the interpolating kernel
% (kx1,ky1): size of the Hamming window
% ma: the number of samples in the alpha domain
% mb: ” ” ” ” ” ” beta “
% (da,db): sample spacing in the (alpha,beta) domain
% a: alpha array
% b: beta array
% g: 2D array of avialble data in the (alpha,beta) domain
% kx: transformation in the kx domain; tkx: scalar form
% ky: transformation in the ky domain; tky: scalar form
% jac: Jacobian of the transformation
% ff: F.T. of the target function
% f: target function
%
% (ax,by): SIZE OF RECTANGULAR PULSE; REMOVE THIS LINE IN FUTURE USAGE
ax=1; by=.5;
%
% REPLACE THE FOLLOWING LINE WITH YOUR OWN PARAMETERS
nx=320; ny=160; x0=3; y0=1.5; ma=32; mb=50; da=pi/3; db=pi/(mb-1);
%
ff=0; f=0; kx=0; ky=0; a=0; b=0; g=0; mm=3;
nnx=nx/2+1; nny=ny/2+1; dx=(2*x0)/nx; dy=(2*y0)/ny; dkx=pi/x0; dky=pi/y0;
kx0=pi/dx; ky0=pi/dy; kx1=dkx*(mm+1); ky1=dky*(mm+1);
%
for i=1:nx; for j=1:ny; ff(i,j)=0; end; end;
for i=1:ma; a(i)=da*(i-1-ma/2); end;
for i=1:mb; b(i)=db*(i-1-mb/2); end;
%
for i=1:ma; for j=1:mb;
%
% MODIFY THE FOLLOWING LINE FOR OTHER TRANSFORMATIONS
tkx=a(i)*cos(b(j)); tky=a(i)*sin(b(j)); jac=abs(a(i)); % *** POLAR TRANS.
%
ij=(i-1)*mb+j; kx(ij)=tkx; ky(ij)=tky;
%
% REMOVE THE FOLLOWING 3 LINES; g(i,j) SHOLD BE INPUT FROM ANOTHER PROGRAM
if tkx == 0, sincx=1; else sincx=sin(tkx*ax)/(tkx*ax); end;
if tky == 0, sincy=1; else sincy=sin(tky*by)/(tky*by); end;
g(i,j)=sincx*sincy; % ***** INPUT
%
cp=g(i,j)*jac*da*db*(.97*2*x0)*(.97*2*y0);
%
% INTERPOLATION VIA UFR
%
fii=(tkx/dkx)+nnx;
fjj=(tky/dky)+nny;
inf=round(fii);
jnf=round(fjj);
for iii=-mm:mm;
iif=inf+iii;
if(iif >= 1 & iif <= nx),
cx=dkx*(iif-nnx);
d1=abs(cx-tkx);
if d1 <= kx1,
for jjj=-mm:mm;
jjf=jnf+jjj;
if(jjf >= 1 & jjf <= ny),
cy=dky*(jjf-nny);
d2=abs(cy-tky);
if d2 <= ky1,
w1=.54+.46*cos((pi*d1)/kx1);
w2=.54+.46*cos((pi*d2)/ky1);
if d1 == 0, sincx=1; else, sincx=sin(.97*x0*d1)/(.97*x0*d1); end;
if d2 == 0, sincy=1; else, sincy=sin(.97*y0*d2)/(.97*y0*d2); end;
ct=w1*w2*sincx*sincy*cp;
ff(iif,jjf)=ff(iif,jjf)+ct;
end;
end;
end;
end;
end;
end;
end; end;
%
% INVERSE TRANSFORM TO OBTAIN THE IMAGE
%
f = fftshift(ifft2(fftshift(ff)));
%
% FOLLOWING SHOWS YOU THE CONTOUR OF AVAILABLE DATA IN (Kx,Ky) DOMAIN
%
plot(kx,ky,’.');
figure(2)
mesh(abs(ff))
figure(3)
mesh(abs(f))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

My first attempt was to try reducing the number of loops this code had. Here’s the result:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function my_ufr_interp
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INTERPOLATION VIA UFR %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This matlab program performs interpolation from nonuniform data stored
% in array g(alpha,beta). In its present form, the target function is a
% rectangular pulse. The transformation is POLAR.
%
% (x0,y0): target’s support
% (nx,ny): target’s image size
% mm: the number of samples within the Hamming window that cuts
% off the interpolating kernel
% (kx1,ky1): size of the Hamming window
% ma: the number of samples in the alpha domain
% mb: ” ” ” ” ” ” beta “
% (da,db): sample spacing in the (alpha,beta) domain
% a: alpha array
% b: beta array
% g: 2D array of avialble data in the (alpha,beta) domain
% kx: transformation in the kx domain; tkx: scalar form
% ky: transformation in the ky domain; tky: scalar form
% jac: Jacobian of the transformation
% ff: F.T. of the target function
% f: target function
%
% (ax,by): SIZE OF RECTANGULAR PULSE; REMOVE THIS LINE IN FUTURE USAGE
ax=1; by=.5;
%
% REPLACE THE FOLLOWING LINE WITH YOUR OWN PARAMETERS
nx=320; ny=160; x0=3; y0=1.5; ma=32; mb=50; da=pi/3; db=pi/(mb-1);
%
ff=0; f=0; kx=0; ky=0; a=0; b=0; g=0; mm=3;
nnx=nx/2+1; nny=ny/2+1; dx=(2*x0)/nx; dy=(2*y0)/ny; dkx=pi/x0; dky=pi/y0;
kx0=pi/dx; ky0=pi/dy; kx1=dkx*(mm+1); ky1=dky*(mm+1);
%
ff=zeros(nx,ny);
a=da*(-ma/2:ma/2-1);
b=db*(-mb/2:mb/2-1);

[th,r]=meshgrid(b,a);
[tkx,tky]=pol2cart(th,r);

figure(1)
plot(tkx,tky,’.');

%tkx=round(tkx/dkx)*dkx;
%tky=round(tky/dky)*dky;

g=sparse(sinc(tkx*ax/pi).*sinc(tky*by/pi).*(abs(tkx)<=kx0).*(abs(tky)<=ky0));
jac=abs(r);
cp=g.*jac*da*db*(2*x0)*(2*y0);

kx=dkx*(-nx/2:nx/2-1);
ky=dky*(-ny/2:ny/2-1);

% for m=1:nx,
% d1=kx(m)-tkx;
% w1=abs(d1)<=dkx;
% % w1=abs(d1)<=kx1+dkx;
% for n=1:ny,
% d2=ky(n)-tky;
% w2=abs(d2)<=dky;
% % w2=abs(d2)<=ky1+dky;
% I=sparse(sinc(x0*sparse(d1.*w1)/pi).*sinc(y0*sparse(d2.*w2)/pi).*w1.*w2);
% F(m,n)=full(sum(sum(cp.*I)));
% end
% end

for m=1:nx,
d1=kx(m)-tkx;
w1=abs(d1)<=dkx;
% w1=abs(d1)<=kx1+dkx;
for n=1:ny,
d2=ky(n)-tky;
w2=abs(d2)<=dky;
% w2=abs(d2)<=ky1+dky;
I=sparse(sinc(sparse(x0*d1.*w1+y0*d2.*w2)/pi).*w1.*w2);
F(m,n)=full(sum(sum(cp.*I)));
end
end

% figure(2)
% imagesc(w1)
%
% figure(3)
% imagesc(w2)

figure(4)
mesh(F)

figure(5)
mesh(real(fftshift(ifft2(fftshift(F)))))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I wasn’t too happy with this – it was too slow. What I needed was an idea. In order to find it, I took the original code and reduced it to a one-dimensional version:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INTERPOLATION VIA UFR %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This matlab program performs interpolation from nonuniform data stored
% in array g(alpha,beta). In its present form, the target function is a
% rectangular pulse. The transformation is POLAR.
%
% (x0,y0): target’s support
% (nx,ny): target’s image size
% mm: the number of samples within the Hamming window that cuts
% off the interpolating kernel
% (kx1,ky1): size of the Hamming window
% ma: the number of samples in the alpha domain
% mb: ” ” ” ” ” ” beta “
% (da,db): sample spacing in the (alpha,beta) domain
% a: alpha array
% b: beta array
% g: 2D array of avialble data in the (alpha,beta) domain
% kx: transformation in the kx domain; tkx: scalar form
% ky: transformation in the ky domain; tky: scalar form
% jac: Jacobian of the transformation
% ff: F.T. of the target function
% f: target function
%
clear all
% (ax,by): SIZE OF RECTANGULAR PULSE; REMOVE THIS LINE IN FUTURE USAGE
ax=1;
%
% REPLACE THE FOLLOWING LINE WITH YOUR OWN PARAMETERS
nx=32; x0=3; ma=32; da=pi/3;
%
ff=0; f=0; kx=0; a=0; g=0; mm=3;
nnx=nx/2+1; dx=(2*x0)/nx; dkx=pi/x0;
kx0=pi/dx; kx1=dkx*(mm+1);
%
for i=1:nx; ff(i)=0; end;
for i=1:ma; a(i)=da*(i-1-ma/2); end;
%
for i=1:ma;
%
% MODIFY THE FOLLOWING LINE FOR OTHER TRANSFORMATIONS
tkx=a(i)^3/8; jac=abs(3*a(i)^2/8); % *** CUBIC TRANS.
%
kx(i)=tkx;
%
% REMOVE THE FOLLOWING 3 LINES; g(i,j) SHOLD BE INPUT FROM ANOTHER PROGRAM
if tkx == 0, sincx=1; else sincx=sin(tkx*ax)/(tkx*ax); end;
g(i)=sincx; % ***** INPUT
%
cp=g(i)*jac*da*(.97*2*x0);
%
% INTERPOLATION VIA UFR
%
fii=(tkx/dkx)+nnx;
inf=round(fii);
fiif(i)=inf;
for iii=-mm:mm;
iif=inf+iii;
if(iif >= 1 & iif <= nx),
cx=dkx*(iif-nnx);
d1=abs(cx-tkx);
if d1 <= kx1,
if d1 == 0, sincx=1; else, sincx=sin(.97*x0*d1)/(.97*x0*d1); end;
dd1(i)=d1;
ct=sincx*cp;
ff(iif)=ff(iif)+ct;
end;
end;
end;
end;
%
% INVERSE TRANSFORM TO OBTAIN THE IMAGE
%
f = fftshift(ifft2(fftshift(ff)));
%
% FOLLOWING SHOWS YOU THE CONTOUR OF AVAILABLE DATA IN (Kx,Ky) DOMAIN
%
figure(1)
stem(kx(1:22),dd1);
figure(2)
plot(kx,fiif,’.')
figure(3)
plot(kx,kx,’.')
kx(end)-kx(1)
size(dd1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

For whatever reason I transformed the code above into this:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INTERPOLATION VIA UFR %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This matlab program performs interpolation from nonuniform data stored
% in array g(alpha,beta). In its present form, the target function is a
% rectangular pulse. The transformation is POLAR.
%
% (x0,y0): target’s support
% (nx,ny): target’s image size
% mm: the number of samples within the Hamming window that cuts
% off the interpolating kernel
% (kx1,ky1): size of the Hamming window
% ma: the number of samples in the alpha domain
% mb: ” ” ” ” ” ” beta “
% (da,db): sample spacing in the (alpha,beta) domain
% a: alpha array
% b: beta array
% g: 2D array of avialble data in the (alpha,beta) domain
% kx: transformation in the kx domain; tkx: scalar form
% ky: transformation in the ky domain; tky: scalar form
% jac: Jacobian of the transformation
% ff: F.T. of the target function
% f: target function
%
clear all
% (ax,by): SIZE OF RECTANGULAR PULSE; REMOVE THIS LINE IN FUTURE USAGE
ax=1;
%
% REPLACE THE FOLLOWING LINE WITH YOUR OWN PARAMETERS
nx=32; x0=3; ma=32; da=pi/3;
%
ff=0; f=0; kx=0; a=0; g=0; mm=3;
nnx=nx/2+1; dx=(2*x0)/nx; dkx=pi/x0;
kx0=pi/dx; kx1=dkx*(mm+1);
%
for i=1:nx; ff(i)=0; end;
for i=1:ma; a(i)=da*(i-1-ma/2); end;
%
for i=1:ma;
%
% MODIFY THE FOLLOWING LINE FOR OTHER TRANSFORMATIONS
tkx=a(i)^3/8; jac=abs(3*a(i)^2/8); % *** CUBIC TRANS.
%
kx(i)=tkx;
%
% REMOVE THE FOLLOWING 3 LINES; g(i,j) SHOLD BE INPUT FROM ANOTHER PROGRAM
if tkx == 0, sincx=1; else sincx=sin(tkx*ax)/(tkx*ax); end;
g(i)=sincx; % ***** INPUT
%
cp=g(i)*jac*da*(.97*2*x0);
%
% INTERPOLATION VIA UFR
%
fii=(tkx/dkx)+nnx;
inf=round(fii);
% fiif(i)=inf;
for iii=-mm:mm;
iif=inf+iii;
if(iif >= 1 & iif <= nx),
cx=dkx*(iif-nnx);
d1=abs(cx-tkx);
if d1 <= kx1,
w=.54+.46*cos((pi*d1)/kx1);
if d1 == 0, sincx=1; else, sincx=sin(.97*x0*d1)/(.97*x0*d1); end;
dd1(i,iif)=d1;
gg(i,iif)=g(i);
ct=w*sincx*cp;
gg1(i,iif)=ct;
ff(iif)=ff(iif)+ct;
end;
end;
end;
end;
%
% INVERSE TRANSFORM TO OBTAIN THE IMAGE
%
f = fftshift(ifft2(fftshift(ff)));
%
% FOLLOWING SHOWS YOU THE CONTOUR OF AVAILABLE DATA IN (Kx,Ky) DOMAIN
%
figure(1)
plot(kx,g);
%pcolor(dd1)
figure(2)
plot(dx*(-nx/2:nx/2-1),f)
%pcolor(gg)
figure(3)
plot(dkx*(-nx/2:nx/2-1),ff)
%pcolor(gg1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

After some mental process, of which I have no recollection anymore, I figured out how to express this code without any loops in it and I applied the results to examples given in the original paper by Soumekh:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function art_ufr

%K=10;
%X0=10;

%alpha=linspace(-2,2,2*K+1);
%mu=linspace(S(-2),S(2),2*K+1);

X0=3;ax=1;
p=1; %auxiliary constant to regulate sampling
nx=32*p; ma=32*p; da=pi/3/p;
%
ff=0; f=0; kx=0; a=0; g=0; mm=3;
nnx=nx/2+1; dx=(2*X0)/nx; dkx=pi/X0;
kx0=pi/dx; kx1=dkx*(mm+1);

alpha=da*(-ma/2:ma/2-1);
mu=dkx*(-nx/2:nx/2-1);

F=2*ax*sinc(mu*ax/pi);
F=F/max(F); %normalize
G=2*ax*sinc(S(alpha)*ax/pi);

figure(1)
subplot(2,1,1)
stem(alpha,G)
xlabel(‘alpha’)
ylabel(‘G(alpha)’)
title(‘Sampled frequency spectrum of a rectangular pulse in alpha and mu domains.’)
subplot(2,1,2)
stem(S(alpha),G)
xlabel(‘mu’)
ylabel(‘G(S(alpha))=G(mu)’)

%Shannon’s reconstruction;
dalpha=da;
[mur alphar]=meshgrid(mu,alpha);
d=(T(mur)-alphar); %argument of the interp. kernel
w=abs(d)<T(kx1); %windowing mask
Fs=G*sparse(sinc(sparse(d.*w)/dalpha).*w);
Fs=Fs/max(Fs); %normalize

figure(2)
plot(mu,real(Fs))
title(‘Shannon”’’s reconstruction for a given span of mu’)
xlabel(‘mu’)
ylabel(‘F_s(mu)’)

%Modified Shannon’s reconstruction;
dalpha=da;
mur=meshgrid(mu);
Fms=Tp(mu+1e-12).*((G.*Sp(alpha))*sinc((T(mur)-alphar)/dalpha));
%infty=find(isinf(Fms)); %test for infinite results
%Fms(infty)=2*X0; %correct results
Fms=Fms/max(Fms); %normalize

figure(3)
plot(mu,real(Fms))
title(‘Modified Shannon”’’s reconstruction for a given span of mu’)
xlabel(‘mu’)
ylabel(‘F_{ms}(mu)’)

%Soumekh’s reconstruction
%dalpha=2*pi/K;
dalpha=da*pi;
[mur alphar]=meshgrid(mu,alpha);
d=mur-S(alphar); %argument of the interp. kernel
w=abs(d)<kx1; %windowing mask
Fsou=dalpha*(G.*Sp(alpha))*sparse(sinc(X0*sparse(d.*w)/pi).*w);
Fsou=Fsou/max(Fsou); %normalize

R=G./Tp(S(alpha));

figure(4)
plot(mu,real(Fsou))
title(‘Unified Fourier reconstruction for a given span of mu’)
xlabel(‘mu’)
ylabel(‘F_{ufr}(mu)’)

figure(5)
plot(mur(1,:),F,’+',mur(1,:),Fs,’:',mur(1,:),Fms,mur(1,:),Fsou,’-.’)
title(‘Comparison between interpolations and precise values’)
xlabel(‘mu’)
ylabel(‘F(mu)’)
legend(‘precise’,'Shannon’,'Modified Shannon’,'UFR’)
% plot(real(fftshift(fft(fftshift(R)))))
% plot(S(alpha),R)

function out=S(alpha)
%coordinate transformation
out=alpha.^3/8;

function out=T(mu)
%coordinate transformation
out=sign(mu).*abs(2*mu.^(1/3));

function out=Sp(alpha)
%analytically calculated Jacobian
out=abs(3/8*alpha.^2);

function out=Tp(mu)
%analytically calculated Jacobian
out=abs(2/3*mu.^(-2/3));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Somewhere around here, I had the idea to factor this algorithm using the k-nearest neighbor search. The results yielded an algorithm below and a paper with a nice trip to a conference.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function art_ufr
X0=3;ax=1;
p=5; %auxiliary constant to regulate sampling
nx=32*p; ma=32*p; da=pi/3/p;
%
ff=0; f=0; kx=0; a=0; g=0; mm=30;
nnx=nx/2+1; dx=(2*X0)/nx; dkx=pi/X0;
kx0=pi/dx; kx1=dkx*(mm+1);
%kx1=dkx*.1;

alpha=da*(-ma/2:ma/2-1);
mu=dkx*(-nx/2:nx/2-1);

F=2*ax*sinc(mu*ax/pi);
F=F/max(F); %normalize
G=2*ax*sinc(S(alpha)*ax/pi);

dalpha=da*pi;
const=kx1*100;

% %Soumekh with nearest neighbors engine
ind=interp1(T(mu),(1:length(mu)),alpha,’nearest’,'extrap’);
I=find(abs(mu(ind)-S(alpha))<const);
J=interp1(T(mu),(1:length(mu)),alpha(I),’nearest’,'extrap’);
sum=zeros(ma,nx);
ind2=sub2ind(size(sum),I,J);
sum(ind2)=sinc(X0*(mu(J)-S(alpha(I)))/pi);
Fsou_near=dalpha*(G.*Sp(alpha))*sparse(sum)
Fsou_near=Fsou_near/max(Fsou_near); %normalize

%Nearest neighbor interpolation
Fnear=interp1(alpha,G,mu,’nearest’,'extrap’);
Fnear=Fnear/max(Fnear); %normalize

% %Matrix-free version
ind=interp1(T(mu),(1:length(mu)),alpha,’nearest’,'extrap’);
I=find(abs(mu(ind)-S(alpha))<const); %here can be simply: I=1:ma;
J=interp1(T(mu),(1:length(mu)),alpha(I),’nearest’,'extrap’);
summ=sinc(X0*(mu(J)-S(alpha(I)))/pi);
Fsou_free(J)=dalpha*(G(I).*Sp(alpha(I))).*summ;
Fsou_free=Fsou_free/max(Fsou_free); %normalize

figure(1)
imagesc(sum)
figure(2)
%plot(mu,F,mu,Fsou_near,mu,Fnear,mu,Fsou_free)
plot(Fsou_free)
figure(3)
plot(abs(fftshift(fft(fftshift(Fsou_free)))))

function out=S(alpha)
%coordinate transformation
out=alpha.^3/8;

function out=T(mu)
%coordinate transformation
out=sign(mu).*abs(2*mu.^(1/3));

function out=Sp(alpha)
%analytically calculated Jacobian
out=abs(3/8*alpha.^2);

function out=Tp(mu)
%analytically calculated Jacobian
out=abs(2/3*mu.^(-2/3));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The fact that the algorithm can be factored into several steps based on the k-nearest neighbor search doesn’t mean much – it’s just another way to write this code. Using the full database of k-nearest neighbors, we would get exactly the same result as before. I tried to argue, however, that it will be enough to use only the nearest neighbors to get a sufficient precision. Here is some numerical analysis.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function casino

%K=10;
%X0=10;

%alpha=linspace(-2,2,2*K+1);
%mu=linspace(S(-2),S(2),2*K+1);

X0=3;ax=1;
p=5; %auxiliary constant to regulate sampling
nx=32*p; ma=32*p; da=pi/3/p;
%
ff=0; f=0; kx=0; a=0; g=0; mm=30;
nnx=nx/2+1; dx=(2*X0)/nx; dkx=pi/X0;
kx0=pi/dx; kx1=dkx*(mm+1);
%kx1=dkx*.1;

alpha=da*(-ma/2:ma/2-1);
mu=dkx*(-nx/2:nx/2-1);

F=2*ax*sinc(mu*ax/pi);
F=F/max(F); %normalize
G=2*ax*sinc(S(alpha)*ax/pi);

for dn=1:10,
%Soumekh’s reconstruction
%dalpha=2*pi/K;
dalpha=da*pi;
[mur alphar]=meshgrid(mu,alpha);
d=mur-S(alphar); %argument of the interp. kernel
w=abs(d)<kx1/40*dn; %windowing mask
Fsou=dalpha*(G.*Sp(alpha))*(sinc(X0*(d.*w)/pi).*w);
summm=(sinc(X0*(d.*w)/pi).*w);
Fsou=Fsou/max(Fsou); %normalize

error(dn)=sum((Fsou-F).^2);
operations(dn)=sum(sum(w));
end

dalpha=da*pi;
[mur alphar]=meshgrid(mu,alpha);
d=mur-S(alphar); %argument of the interp. kernel
w=abs(d)<kx1*1000; %windowing mask
Fsou=dalpha*(G.*Sp(alpha))*(sinc(X0*(d.*w)/pi).*w);
summm=(sinc(X0*(d.*w)/pi).*w);
Fsou=Fsou/max(Fsou); %normalize

e=sum((Fsou-F).^2)
op=sum(sum(w))

% R=G./Tp(S(alpha));

%Soumekh’s reconstruction w/ nearest neighbor search
mu=mu.’;
alpha=alpha.’;
m=length(mu);
n=length(alpha);
i=(1:m).’; %create index vector for x
I=i(:,ones(1,n)); %create index matrix
j=1:n; %create index vector for y
J=j(ones(1,m),:); %create index matrix
D=zeros(m,n);
D(:)=sqrt(sum(abs(mu(I(:),:)-S(alpha(J(:),:))).^2,2)); %calculate euclidean distance matrix

[s,I]=sort(D);
i=(1:m);

for kn=1:10,
W=zeros(m,n);
for k=1:kn, %k-nearest neighbors
j=I(k,:);
ind=sub2ind([m n],i,j);
W(ind)=1;
end

Fnn=dalpha*(G.*Sp(alpha.’))*(sinc(X0*(D.’.*W)/pi).*W);
summ=(sinc(X0*(D.’.*W)/pi).*W);
Fnn=Fnn/max(Fnn); %normalize

error_nn(kn)=sum((Fnn-F).^2);
operations_nn(kn)=sum(sum(W));
end

figure(1)
pcolor(abs(w))

figure(2)
pcolor(abs(W))

figure(3)
plot(mu,Fsou,mu,Fnn,mu,F)

figure(1)
plot(operations_nn,error_nn,’*')
title(‘Accuracy vs speed for 1 to 10 nearest neighbors’)
xlabel(‘number of calculations’)
ylabel(‘error’)

function out=S(alpha)
%coordinate transformation
out=alpha.^3/8;

function out=T(mu)
%coordinate transformation
out=sign(mu).*abs(2*mu.^(1/3));

function out=Sp(alpha)
%analytically calculated Jacobian
out=abs(3/8*alpha.^2);

function out=Tp(mu)
%analytically calculated Jacobian
out=abs(2/3*mu.^(-2/3));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

And that’s it! An open question now is how to do this in two, three, n-dimensions. That’s where the power of nearest neighbor search algorithms really pays of. There are two more pieces of code which might be of some help. I wrote them to test some examples in the book “Fourier Array Imaging”.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function jacob
%Soumekh, Fourier Array Imaging, pages 24 to 25

N=500;
X0=2*pi;
x=linspace(-X0,X0,N+1);
y=sin(x.^2);
xt=g(x);

gx=linspace(g(-X0),g(X0),N+1);

yt=interp1(xt,y,gx);
s=y.*J(x); %this line uses function J

%J=diff(xt)*N/(abs(x2)-abs(x(N)); %Jacobian calculated numerically
%J=[J 0];
%s=y.*J;

trapz(x,s)
trapz(gx,yt)

figure(1)
plot(x,x,x,xt,x,interp1(xt,x,xt))
figure(2)
plot(x,y,gx,yt)

dx=2*X0/N;
dxt=2*g(X0)/N;

dkx=2*pi/N/dx;
dkxt=2*pi/N/dxt;

kx=linspace(-N*dkx/2,N*dkx/2,N+1);
kxt=linspace(-N*dkxt/2,N*dkxt/2,N+1);

%kx=linspace(-4*pi,4*pi,N+1); %original freq. span for num. integration
Fn=trapz(x,meshgrid(s).*exp(-j*kxt’*g(x)),2);
Fnt=trapz(gx,meshgrid(yt).*exp(-j*kxt’*gx),2);

figure(3)%ideally, the spectra will be real only (see analytical expressions)
plot(kxt,real(Fn),kxt,real(Fnt),kxt,real(fftshift(fft(fftshift(yt))))*dxt)

%Shannon’s reconstruction; anti-aliasing filter
%gx=gx(1:end-1);
yt=yt(1:end-1);
kx=kx(1:end-1);
kxt=kxt(1:end-1);

n=-N/2:N/2;
n=n(1:end-1);
%n=1:N+1;
kxtt=g(kx);
kxttt=interp1(kxtt,kx,kx);
[kxttt n]=meshgrid(kxttt,n);
Fns=fftshift(fft(yt))*sinc((kxttt-n*dkxt)/dkxt);

%Soumekh’s reconstruction; anti-aliasing filter
%Fnss=dkxt*fftshift(fft(yt))*(sinc((kx-interp1(g((-N/2:N/2)*dkxt),(-N/2:N/2)*dkxt,n*dkxt))/pi));%./J(kx));
%Fnss=dkxt*fftshift(fft(yt))*(sinc((kx-g(n*dkxt))/pi)./J(n*dkxt));

figure(4)
%match vector lengths
x=x(1:end-1);
y=y(1:end-1);
plot(x,real(ifft(fftshift(Fns))),x,y)

function out=g(x)
%coordinate transformation
a=0.01;
out=x+a*x.^3;

function out=J(x)
%analytically calculated Jacobian
a=0.01;
out=abs(1+3*a*x.^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The second piece of code.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function jacob2
%Soumekh, Fourier Array Imaging, pages 26 to 27

global X0

N=100;
X0=50;
x=linspace(-X0,X0,N+1);
%dx=2*X0/N;
dx=2;
[x n]=meshgrid(x,(-N/2:N/2));
y=sum(dirac(x-n*dx))’;
figure(1)
stem(x(1,:),y)

yt=J(x(1,:)).*sum(dirac(g(x)-g(n*dx)));
figure(2)
stem(g(x(1,:)),yt)

function out=dirac(x)
global X0
out=abs(x==0);

function out=g(x)
%coordinate transformation
a=0.01;
out=x+a*x.^3;

function out=J(x)
%analytically calculated Jacobian
a=0.01;
out=abs(1+3*a*x.^2);