Fourier Image Analysis

Settings and Helper Functions

Avoid error messages for similar symbols (e.g. viewI, viewF, viewC, viewM)

Off[General::spell];Off[General::spell1];Off[General::munfl];

Re-align images

recenter[array_]:=RotateLeft[array,Floor[Dimensions[array]/2]]

Create {image,fourier} from source image or fourier-transform image

SetOptions[Fourier,FourierParameters->{1,1}];
SetOptions[InverseFourier,FourierParameters->{1,1}];
ItoIF[imageI_]:={imageI,Chop[Fourier[imageI]]};
FtoIF[imageF_]:={Chop[InverseFourier[imageF]],imageF}

Import color image and convert to {grayscale,fourier}

ImagetoIF[image_]:=ItoIF[recenter[ImageData[ColorSeparate[image,"I"]]]]

Given function of x and y, and range for x and y, create {grayscale,fourier}

createIF[fn_,range_]:=ItoIF[
recenter[Table[fn,{x,-range[[1]]/2,range[[1]]/2-1},{y,-range[[2]]/2,range[[2]]/2-1}]]];
createFI[fn_,range_]:=FtoIF[
recenter[Table[fn,{x,-range[[1]]/2,range[[1]]/2-1},{y,-range[[2]]/2,range[[2]]/2-1}]]];

Create {image,fourier} from function[x,y]

convolve[imageA_,imageB_] := FtoIF[imageA[[2]] imageB[[2]]];
multiply[imageA_,imageB_] := ItoIF[imageA[[1]] imageB[[1]]]

Functions for viewing image and Fourier transform data.
View image, with color values scaled to fill range. View fourier transform, scaled so mean+1/2 standard deviation is white (otherwise, DC terms swamps other terms and you just get an image with one visible dot in the middle)

viewI[im_]:=ImageAdjust[Image[recenter[Abs[im[[1]]]]]];
IScale[im_]:=Block[{a=Flatten[Abs[im[[2]]]]},Mean[a]+0.5StandardDeviation[a]];
viewF[im_,scale_:0]:=Block[
{s=scale},If[scale==0,s=IScale[im]];
ImageAdjust[Image[recenter[Abs[im[[2]]]]],1,{0,s}]];

Display image, transform and inverse(transform(image))

viewIF[im_]:=GraphicsRow[{viewI[im],viewF[im]}]

Layout with two source images pairs (image/fourier) and result image

view3[imageA_,imageB_,imageC_]:=Block[{sc=IScale[imageA]},GraphicsGrid[{{viewI[imageA],viewI[imageB],viewI[imageC]},{viewF[imageA,sc],viewF[imageB],viewF[imageC,sc]}},ImageSize->Full]]

Basic transforms

Manipulate[
Plot[{
c0 Cos[0 x+p0]+c1 Cos[1 x+p1]+c2 Cos[2 x+p2] +c3 Cos[3x+p3]+c4 Cos[4x+p4],
c0 Cos[0 x+p0],c1 Cos[1 x+p1],c2 Cos[2 x+p2] ,c3 Cos[3x+p3],c4 Cos[4x+p4] },
{x,-2π,2π},PlotPoints->100,PlotRange->{-1.5,1.5},ExclusionsStyle->Black],
{{c0,0},-1.5,1.5},{{p0,0},-π,π},
{{c1,0},-1.5,1.5},{{p1,0},-π,π},
{{c2,0},-1.5,1.5},{{p2,0},-π,π},
{{c3,0},-1.5,1.5},{{p3,0},-π,π},
{{c4,0},-1.5,1.5},{{p4,0},-π,π}
]

Manipulate[
Plot[{(π/4)Sign[Cos[x]],(∑_(i=0)^n)(((((-1))^i)Cos[(2i+1) x])/(2 i+1))},{x,-2π,2π},PlotPoints->100,ExclusionsStyle->Black],
{n,0,20,1}]

viewIF[createIF[Cos[((x+32)/2+y)/8*π],{64,64}]]

Manipulate[viewIF[createIF[Cos[π(x (2^(-xscale))+y (2^(-yscale))+offset)],{64,64}]],
{{xscale,5},0,5,1},{{yscale,3},0,5,1},{offset,0,64}]

Manipulate[viewIF[createFI[If[{x,y}=={x0,y0}||{x,y}==-{x0,y0},1,0],{64,64}]],
{{x0,5},0,32,1},{{y0,3},0,32,1}]

GraphicsRow[{Plot[#,{x,-2,2},ExclusionsStyle->Dashed],Plot[FourierTransform[#,x,t],{t,-30,30},PlotRange->Full,ExclusionsStyle->Dashed]}]&[UnitBox[x]]

GraphicsRow[{Plot[#,{x,-2,2}],Plot[FourierTransform[#,x,t],{t,-30,30},PlotRange->Full]}]&[UnitTriangle[x]]

GraphicsRow[{Plot[#,{x,-2,2}],Plot[FourierTransform[#,x,t],{t,-30,30},PlotRange->Full]}]&[1/(1+(x^2))]

image=ImagetoIF[Import[SystemDialogInput["FileOpen",NotebookDirectory[]]]];
viewIF[image]

square=createIF[If[Abs[x]<8&&Abs[y]<8,1/(17^2),0],Dimensions[image[[1]]]];
viewIF[square]

Convolution

DynamicModule[{f,g,c},Manipulate[
c[x_]:=Evaluate[Convolve[f[t],g[t],t,x]];
GraphicsColumn[
{Plot[{f[t],g[x-t],f[t]g[x-t]},{t,-2,2},PlotPoints->100,PlotRange->{-1,1},Filling->{3->Axis},ExclusionsStyle->Dashed],
Show[
Plot[{c[x]},{x,-2,2},PlotPoints->100,PlotRange->{-1,1}],
Graphics[{PointSize[Large],Red,Point[{x,c[x]}]}]]
}],
{f,{UnitBox->"Box",UnitTriangle->"Triangle",Sinc[8#]&->"Sinc",Exp[-50((#)^2)]&->"Narrow Gaussian",Exp[-2((#)^2)]&->"Wide Gaussian",(DiracDelta[#-.75]+DiracDelta[#+.75])&->"Samples"}},
{g,{UnitBox->"Box",UnitTriangle->"Triangle",Sinc[8#]&->"Sinc",Exp[-50((#)^2)]&->"Narrow Gaussian",Exp[-2((#)^2)]&->"Wide Gaussian",(DiracDelta[#-.75]+DiracDelta[#+.75])&->"Samples"}},
{{x,0},-2,2},SynchronousUpdating->False]]

DynamicModule[{f,g,c},Manipulate[
c[x_]:=Evaluate[Convolve[f[t],g[t],t,x]];
GraphicsColumn[
{Show[
Plot[{f[t],f[x]g[t-x]},{t,-2,2},PlotPoints->100,PlotRange->{-1,1},Filling->{2->Axis}],
Graphics[{PointSize[Large],Red,Point[{x,f[x]}]}]],
Plot[{c[t],f[x]g[t-x]},{t,-2,2},PlotPoints->100,PlotRange->{-1,1},Filling->{2->Axis}]
}],
{f,{UnitBox->"Box",UnitTriangle->"Triangle",Sinc[8#]&->"Sinc",Exp[-50((#)^2)]&->"Narrow Gaussian",Exp[-2((#)^2)]&->"Wide Gaussian",(DiracDelta[#-.75]+DiracDelta[#+.75])&->"Samples"}},
{g,{UnitBox->"Box",UnitTriangle->"Triangle",Sinc[8#]&->"Sinc",Exp[-50((#)^2)]&->"Narrow Gaussian",Exp[-2((#)^2)]&->"Wide Gaussian",(DiracDelta[#-.75]+DiracDelta[#+.75])&->"Samples"}},
{{x,0},-2,2},SynchronousUpdating->False]]

view3[image,square,convolve[image,square]]

circle=createIF[If[(x^2)+(y^2)<(8^2),(1/(π (8^2))),0],Dimensions[image[[1]]]];
viewIF[circle]

view3[image,circle,convolve[image,circle]]

fsquare=createFI[If[Max[Abs[x],Abs[y]]<20,1,0],Dimensions[image[[1]]]];
viewIF[fsquare]

view3[image,fsquare,convolve[image,fsquare]]

gaussian=createIF[N[(1/(16(Sqrt[2π])))Exp[-(((x^2)+(y^2))/(2*(16^2)))]],Dimensions[image[[1]]]];
viewIF[gaussian]

view3[image,gaussian,convolve[image,gaussian]]

Sampling

sample=createIF[If[Mod[x,16]==8 && Mod[y,16]==8,16,0],Dimensions[image[[1]]]];
viewIF[sample]

Display image, transform and inverse(transform(image))

sampled=multiply[image,sample];
view3[image,sample,sampled]

view3[sampled,square,convolve[sampled,square]]

view3[sampled,gaussian,convolve[sampled,gaussian]]

lowpass=createFI[If[Abs[x]<Dimensions[image][[2]]/32&&Abs[y]<Dimensions[image][[3]]/32,1,0],Dimensions[image[[1]]]];
viewIF[lowpass]

view3[sampled,lowpass,convolve[sampled,lowpass]]

lowpassImage=convolve[image,lowpass];
view3[image,lowpass,lowpassImage]

sampledLowpass=multiply[lowpassImage,sample];
view3[lowpassImage,sample,sampledLowpass]

view3[sampledLowpass,lowpass,convolve[sampledLowpass,lowpass]]