Clear[Evaluate[Context[]<>"*"]];
s=20;
n[x_]:=Block[{ma,mi},
ma=Max@x;
mi=Min@x;
(x-mi)/(ma-mi)
];
data[a_,α_,b_]:=Table[1/(1+(i-7)^2+(j-5)^2)+1/(1+(i+2)^2+(j-2)^2),{i,-s,s},{j,-s,s}];
d=data[{0,15},π/7,10];
Print@Image@n@d;
mf=ImageData@Radon@Image@d;
Print@Image@n@mf;
fd=Fourier[d];
fd=fd*Table[(-1)^(i+j),{i,-s,s},{j,-s,s}];
Print@Image@n@Map[Abs,fd,{2}];
fdi[x_,y_]:=Block[{xi,yi,xip,yip,xf,yf,xin,xipn,yin,yipn},
xi=Floor@x;
yi=Floor@y;
xf=x-xi;
yf=y-yi;
xin=If[Abs@xi≤s,If[xi≥0,1,-1],0];
yin=If[Abs@yi≤s,If[yi≥0,1,-1],0];
xipn=If[Abs[xi+1]≤s,If[xi+1≥0,1,-1],0];
yipn=If[Abs[yi+1]≤s,If[yi+1≥0,1,-1],0];
xip=Mod[xi+2,2s+1,1];
yip=Mod[yi+2,2s+1,1];
xi=Mod[xi+1,2s+1,1];
yi=Mod[yi+1,2s+1,1];
xin*yin*fd[[yi,xi]]*(1-xf)*(1-yf)
+xin*yipn*fd[[yip,xi]]*(1-xf)*yf
+xipn*yin*fd[[yi,xip]]*xf*(1-yf)
+xipn*yipn*fd[[yip,xip]]*xf*yf
];
fr=Table[fdi[r*Cos@β,r*Sin@β],{r,-2s,2s},{β,0,π,0.02}];
β=-π/2.2;
Plot[{Abs@fdi[r Cos@β,r Sin@β],Re@fdi[r Cos@β,r Sin@β],Im@fdi[r Cos@β,r Sin@β]},{r,-10,10},PlotRange->All]
Image@n@Map[Abs,fr,{2}]
fr=Table[(-1)^i,{i,-2s,2s}]*fr;
r=Transpose@Map[InverseFourier,Transpose@fr];
Image@n@Map[Abs,r,{2}]