Gabor patch visibility in CPD using fourier analysis - @binocularity Feb 2024
​
https://cdn.sinauer.com/wolfe4e/wa03.03.html
https://mathematica.stackexchange.com/questions/88168/averaging-over-a-circle-from-fft-to-plot-the-amplitude-vs-wavelength
https://mathematica.stackexchange.com/questions/87214/calculation-of-2d-fft-for-an-image/87240#87240
​
In[]:=
SetDirectory[NotebookDirectory[]];​​​​(*Constants*)​​maxBinormal=Max[Table[PDF[BinormalDistribution[0.0],{x,y}],{x,-Pi,Pi},{y,-Pi,Pi}]];​​scaleBinormal=1/maxBinormal;
In[]:=
(*Functions*)​​gaborFunction[f_,x_,y_,vector_,contrast_]:=Module[{vec=Normalize[vector],c=Clip[contrast,{0,1}],u,v},​​u=vec[[1]];​​v=vec[[2]];​​c*Sin[fxu+fyv]*(scaleBinormal*PDF[BinormalDistribution[0.0],{x,y}])]​​​​gaborPatch2D[f_,vector_,contrast_,res_]:=Module[{},​​Rasterize[DensityPlot[((gaborFunction[f,x,y,vector,contrast])+1)/2,​​{x,-Pi,Pi},{y,-Pi,Pi},​​MaxRecursion->1,​​ColorFunction->GrayLevel,ColorFunctionScaling->False,​​Frame->False,​​PlotPoints->res,PlotRangePadding->0,PlotRange->All],​​ImageSize->res/2,RasterSize->res]​​]
In[]:=
contrast=1;​​frequency=10;​​vector={1.0,0.0};​​resolution=260;​​img=gaborPatch2D[frequency,vector,contrast,resolution]
Out[]=
In[]:=
(*Fourieranalysistocalculatefrequencypowerspectrum*)​​imgRaw=ImageData[RemoveAlphaChannel[ColorConvert[img,"GrayScale"]]];​​imDat={imgRaw};​​ftImg=Fourier[imDat[[1]]];​​absFtImg=Abs[ftImg];​​fft1=RotateLeft[absFtImg,Floor[Dimensions[absFtImg]/2]];​​Image[fft1,ImageSize->resolution/2]
Out[]=
In[]:=
radius=resolution/2;​​fftPolar=ImageTransformation[Image[Abs[fft1]],{Cos[#[[1]]],Sin[#[[1]]]}#[[2]]&,Round[{π*radius,radius}],PlotRange->{{0,2π},{0,1}},DataRange->{{-1,1},{-1,1}}]​​
Out[]=
In[]:=
(*ListPlot[Mean/@ImageData[fftPolar,DataReversed->True],DataRange->{0,radius},Joined->True,PlotRange->{Automatic,{0,10}}]*)
In[]:=
x=Mean/@ImageData[fftPolar,DataReversed->True];​​peaks=FindPeaks[x];​​peaks=Reverse[#]&/@peaks;​​(*Selectpeaksover1hz(notDC)andwithnon-negligablepower*)​​keyPeaks=Select[peaks,(#[[2]]>1&&#[[1]]>1/130)&];​​powerMax=Max[keyPeaks[[All,1]]];​​plotMax=Ceiling[powerMax];​​freqMax=Select[keyPeaks,(#[[1]]>=powerMax)&][[1]][[2]];​​y=Range[0,radius-1,1];​​freqPowerSpectrum=Partition[Riffle[x,y],2];​​​​ListPlot[{freqPowerSpectrum,peaks},PlotStyle->{Automatic,PointSize[0.01]},Joined->{True,False},PlotRange->{{0,plotMax},{0,130}},AxesLabel->{"Power","Cycles across image"},PlotLabel->StringJoin["Frequency Power Spectrum,\n Max power at ",ToString[freqMax],"hz (cycles across image)"]]​​
Out[]=
In[]:=
​
What is freqMax in CPD (cycles per degree of visual angle)?
It depends on the display, assume we have:
14inch HD (1920x1080 pixels, 16:9 aspect ratio) laptop display viewed from 300mm
displaying a 260x260 Gabor patch with no magnification factor
​
14inch is 35.56cm diagonal display or 31cmx17.4cm dimensions
This is a pixel pitch of 0.162mm or a dpi of 157 https://pixelcalculator.com/en
​
Working in mm
In[]:=
viewingDistance=300;​​displayDiagonal=14*25.4;​​displayWidth=displayDiagonal/Sqrt[16^2+9^2]*16;​​displayPixelPitch=displayWidth/1920
Out[]=
0.161423
In[]:=
gaborPatchSize=resolution*displayPixelPitch
Out[]=
41.97
In[]:=
visualAngle=ArcTan[(gaborPatchSize/2)/viewingDistance]/Degree
Out[]=
4.00132
In[]:=
CPD=freqMax/visualAngle
Out[]=
2.49917
Typical limit of CPD for human vision at high contrast is 10 cpd.