ReadImage/: ReadImage[filename_, n_] := 
   ReadList[filename, Table[Number, {n}]]
DispIt/: DispIt[t_] := ListPlot3D[t, Lighting -> False, 
    Mesh -> False, AspectRatio -> 3, 
    ViewPoint -> {0, 0, -1000}, Axes -> False]
DispIt/: DispIt[t_, a_] := ListPlot3D[t, Lighting -> False, 
    Mesh -> False, AspectRatio -> a, 
    ViewPoint -> {0, 0, -1000}, Axes -> False]
EpsOut/: EpsOut[which_] := Display["!psfix -epsf > out.ps", 
    which]
SetOptions[ListPlot, PlotJoined -> True];

(* synthesise a simple 3D image *)
nptsz = nptsy = nptsx = 32;
ims=Table[Table[Table[50 Random[],{nptsx}],{nptsy}],{nptsz}];
zstart=Floor[(nptsz/4)];
ystart=Floor[(nptsy/4)];
xstart=Floor[(nptsx/4)];
zw=Floor[(nptsz/4)];
yw=Floor[(nptsy/4)];
xw=Floor[(nptsx/4)];
zcen=Floor[((3 nptsz)/4)];
ycen=Floor[((3 nptsy)/4)];
xcen=Floor[((3 nptsx)/4)];
For[p=1,p≤nptsz,p++,For[m=1,m≤nptsy,m++,For[n=1,n≤nptsx,n++,ims[[m,n,p]]+=If[(((m-ycen))^2)+(((n-xcen))^2)+(((p-zcen))^2)<0.8 (xw^2),200,0];]]];
For[p=1,p≤zw,p++,For[m=1,m≤yw,m++,For[n=1,n≤xw,n++,ims[[ystart+m,xstart+n,zstart+p]]+=100]]];

(* Use synthetic image *)
im = ims;

{nptsx,nptsy,nptsz}
For[p=1,p≤nptsz,p++,
ListDensityPlot[im[[p]],Mesh->False,PlotRange->{0,350}];
]

(******************************************************
Now the Fourier Versions 
*******************************************************)(* this function rotates quadrants of a 2D image *)
RotImage2D[im_]:=Module[{nptsx,nptsy},nptsy=Length[im];nptsx=Length[im[[1]]];RotateLeft[Transpose[RotateLeft[Transpose[im],Floor[(nptsx/2)+1]]],Floor[(nptsy/2)+1]]]
(* inverse of the above *)
IRotImage2D[im_]:=Module[{nptsx,nptsy},nptsy=Length[im];nptsx=Length[im[[1]]];RotateRight[Transpose[RotateRight[Transpose[im],Floor[(nptsx/2)+1]]],Floor[(nptsy/2)+1]]]

Table[Mod[j+Floor[nptsz/2 ],nptsz]+1,{j,nptsz}]

RotImage3D[im_] := Table[RotImage2D[im[[Mod[j+Floor[nptsz/2 ],nptsz]+1]]],{j,nptsz}]

fim=Fourier[im];
(*
fim2=RotImage[fim];
ListDensityPlot[Re[Log[fim2]],Mesh->False];
ListDensityPlot[Im[Log[fim2]],Mesh->False];
*)

(* now form derivative of Gaussian edge filters *)

G[x_,s_]:=((E^(-((x^2)/(2 (s^2)))))/(Sqrt[2 π]))
Gx[x_,s_]:=(∂_xx)G[xx,s]/.xx->x
Gxx[x_,s_]:=(∂_({xx,2}))G[xx,s]/.xx->x

xmaskim=Table[Table[Table[0,{nptsx}],{nptsy},{nptsz}]];
zmaskim = ymaskim = xmaskim;
zzmaskim = yymaskim =xxmaskim = xmaskim;
zxmaskim = zymaskim = xymaskim = xmaskim;
ycen=(nptsy/2);
xcen=(nptsx/2);
zcen=(nptsz/2);
sx=1.5;
sy=1.5;
sz =1.5;
For[p=1,p≤nptsz,p++,For[m=1,m≤nptsy,m++,For[n=1,n≤nptsx,n++,xmaskim[[p,m,n]]=G[p-zcen,sz] G[m-ycen,sy] Gx[n-xcen,sx]]]];
ListDensityPlot[xmaskim[[xcen//Floor]],Mesh->False,PlotRange->All];

For[p=1,p≤nptsz,p++,For[m=1,m≤nptsy,m++,For[n=1,n≤nptsx,n++,ymaskim[[p,m,n]]=G[p-zcen,sz] Gx[m-ycen,sy] G[n-xcen,sx]]]];
ListDensityPlot[ymaskim[[ycen//Floor]],Mesh->False,PlotRange->All];

For[p=1,p≤nptsz,p++,For[m=1,m≤nptsy,m++,For[n=1,n≤nptsx,n++,zmaskim[[p,m,n]]=Gx[p-zcen,sz] G[m-ycen,sy] G[n-xcen,sx]]]];
ListDensityPlot[zmaskim[[zcen//Floor]],Mesh->False,PlotRange->All];


(* 2nd derivative images *)
For[p=1,p≤nptsz,p++,For[m=1,m≤nptsy,m++,For[n=1,n≤nptsx,n++,xxmaskim[[p,m,n]]=G[p-zcen,sz] G[m-ycen,sy] Gxx[n-xcen,sx]]]];
ListDensityPlot[xxmaskim[[xcen//Floor]],Mesh->False,PlotRange->All];

For[p=1,p≤nptsz,p++,For[m=1,m≤nptsy,m++,For[n=1,n≤nptsx,n++,yymaskim[[p,m,n]]=G[p-zcen,sz] Gxx[m-ycen,sy] G[n-xcen,sx]]]];
ListDensityPlot[yymaskim[[ycen//Floor]],Mesh->False,PlotRange->All];

For[p=1,p≤nptsz,p++,For[m=1,m≤nptsy,m++,For[n=1,n≤nptsx,n++,zzmaskim[[p,m,n]]=Gxx[p-zcen,sz] G[m-ycen,sy] G[n-xcen,sx]]]];
ListDensityPlot[zzmaskim[[zcen//Floor]],Mesh->False,PlotRange->All];

For[p=1,p≤nptsz,p++,For[m=1,m≤nptsy,m++,For[n=1,n≤nptsx,n++,xymaskim[[p,m,n]]=G[p-zcen,sz] Gx[m-ycen,sy] Gx[n-xcen,sx]]]];
ListDensityPlot[xymaskim[[zcen//Floor]],Mesh->False,PlotRange->All];

For[p=1,p≤nptsz,p++,For[m=1,m≤nptsy,m++,For[n=1,n≤nptsx,n++,zxmaskim[[p,m,n]]=Gx[p-zcen,sz] G[m-ycen,sy] Gx[n-xcen,sx]]]];
ListDensityPlot[zxmaskim[[zcen//Floor]],Mesh->False,PlotRange->All];

For[p=1,p≤nptsz,p++,For[m=1,m≤nptsy,m++,For[n=1,n≤nptsx,n++,zymaskim[[p,m,n]]=Gx[p-zcen,sz] Gx[m-ycen,sy] G[n-xcen,sx]]]];
ListDensityPlot[zymaskim[[zcen//Floor]],Mesh->False,PlotRange->All];

xmaskim[[1,1,1]]//N

dim = zymaskim;
max = -10^(6);
min  = -max;
(max =If[#> max,#,max])&/@Flatten[dim];
(min =If[#< min,#,min])&/@Flatten[dim];
{min,max}
For[p=1,p≤nptsz,p++,
ListDensityPlot[dim[[p]],Mesh->False,PlotRange->{min,max}];
]

(* 1st derivatives *)
fxmaskim=Fourier[N[xmaskim]];
fymaskim=Fourier[N[ymaskim]];
fzmaskim=Fourier[N[zmaskim]];
(* 2nd derivatives *)
fxxmaskim=Fourier[N[xxmaskim]];
fyymaskim=Fourier[N[yymaskim]];
fzzmaskim=Fourier[N[zzmaskim]];
fxymaskim=Fourier[N[xymaskim]];
fzxmaskim=Fourier[N[zxmaskim]];
fzymaskim=Fourier[N[zymaskim]];

xfiltimf=Chop[InverseFourier[N[fim fxmaskim]]]//RotImage3D;
ListDensityPlot[xfiltimf[[xcen]],Mesh->False];
yfiltimf=Chop[InverseFourier[N[fim fymaskim]]]//RotImage3D;
ListDensityPlot[yfiltimf[[ycen]],Mesh->False];
zfiltimf=Chop[InverseFourier[N[fim fzmaskim]]]//RotImage3D;
ListDensityPlot[zfiltimf[[zcen]],Mesh->False];
(* 2nd derivatives *)
xxfiltimf=Chop[InverseFourier[N[fim fxxmaskim]]]//RotImage3D;
ListDensityPlot[xxfiltimf[[xcen]],Mesh->False];
yyfiltimf=Chop[InverseFourier[N[fim fyymaskim]]]//RotImage3D;
ListDensityPlot[yyfiltimf[[ycen]],Mesh->False];
zzfiltimf=Chop[InverseFourier[N[fim fzzmaskim]]]//RotImage3D;
ListDensityPlot[zzfiltimf[[zcen]],Mesh->False];

xyfiltimf=Chop[InverseFourier[N[fim fxymaskim]]]//RotImage3D;
ListDensityPlot[xyfiltimf[[xcen]],Mesh->False];
zyfiltimf=Chop[InverseFourier[N[fim fzymaskim]]]//RotImage3D;
ListDensityPlot[zyfiltimf[[ycen]],Mesh->False];
zxfiltimf=Chop[InverseFourier[N[fim fzxmaskim]]]//RotImage3D;
ListDensityPlot[zxfiltimf[[zcen]],Mesh->False];

dim = zzfiltimf;
max = -10^(6);
min  = -max;
(max =If[#> max,#,max])&/@Flatten[dim];
(min =If[#< min,#,min])&/@Flatten[dim];
{min,max}
For[p=1,p≤nptsz,p++,
ListDensityPlot[dim[[p]],Mesh->False,PlotRange->{min,max}];
]

(* Form Gradient and Laplacian images *)
gradim = Sqrt[xfiltimf^2 + yfiltimf^2 + zfiltimf^2];
laplim = xxfiltimf^2 + yyfiltimf^2 + zzfiltimf^2;

max = -10^(6);
min  = -max;
(max =If[#> max,#,max])&/@Flatten[gradim];
(min =If[#< min,#,min])&/@Flatten[gradim];
{min,max}
For[p=1,p≤nptsz,p++,
ListDensityPlot[gradim[[p]],Mesh->False,PlotRange->{0,max}];
]

max = -10^(6);
min  = -max;
(max =If[#> max,#,max])&/@Flatten[laplim];
(min =If[#< min,#,min])&/@Flatten[laplim];
{min,max}
For[p=1,p≤nptsz,p++,
ListDensityPlot[laplim[[p]],Mesh->False,PlotRange->{min,max}];
]

(* Form the Canny image *)

cannyim = (xfiltimf^2 xxfiltimf +  yfiltimf^2 yyfiltimf+  zfiltimf^2 zzfiltimf+ 2 xfiltimf yfiltimf xyfiltimf  +2 xfiltimf zfiltimf zxfiltimf  + 2 zfiltimf yfiltimf zyfiltimf  )/(xfiltimf^2 + yfiltimf^2 + zfiltimf^2);

max = -10^(6);
min  = -max;
(max =If[#> max,#,max])&/@Flatten[cannyim];
(min =If[#< min,#,min])&/@Flatten[cannyim];
{min,max}

ListDensityPlot[cannyim[[zcen]],Mesh->False,PlotRange->{min,max}];

For[p=1,p≤nptsz,p++,
ListDensityPlot[cannyim[[p]],Mesh->False,PlotRange->{min,max}];
]