%Lecture_050712.m
%
%2-D FFTs

%% Making images using meshgrid

nPix = 200;  %number of pixels in each row and column
wDeg = 10;  %degrees of visual angle

[x,y] = meshgrid( linspace(-wDeg/2,wDeg/2,nPix));

figure(1)
clf
subplot(1,2,1)
imagesc(x)
colormap(gray)
colorbar
axis equal
axis tight

subplot(1,2,2)
imagesc(y)
colormap(gray)
colorbar
axis equal
axis tight

%% 'radius' image

img = sqrt(x.^2+y.^2);
figure(2)
clf
showImage(img,x,y)

%% grating
sf = .5; %cycles/deg
img = sin(2*pi*y*sf);
figure(2)
clf
showImage(img,x,y)

%% oriented grating

theta = 0;  %degrees
phase = 90;
ramp = cos(theta*pi/180)*x + sin(theta*pi/180)*y;
img = sin(2*pi*ramp*sf-phase*pi/180);
showImage(img,x,y);

%% Gaussian

center = [-2,-2];
sigma  = 1;

Gauss = exp( -((x-center(1)).^2 + (y-center(2)).^2)/(2*sigma^2));
showImage(Gauss,x,y);

%% Gabor

Gabor = Gauss.*img;
showImage(Gabor,x,y);

%% Circle

circle = (x-center(1)).^2 + (y-center(2)).^2 <= sigma^2;

showImage(circle,x,y);

%% plaid

plaid =  sin(2*pi*x*sf).*sin(2*pi*y*sf);

showImage(plaid,x,y);

%% checkerboard

checkerboard = sign(plaid);
showImage(checkerboard,x,y);

%% radial grating

ang = atan2(y,x);
rad = sqrt(x.^2+y.^2);

img = sin(10*ang);
showImage(img,x,y);

%% 2D FFT of a sinusoid

img = sin(2*pi*rad*sf);
showImage(img,x,y);

sf= 2;
theta = 45;  %degrees
phase = 90;
ramp = cos(theta*pi/180)*x + sin(theta*pi/180)*y;
img = sin(2*pi*ramp*sf-phase*pi/180);

img = img +1;

figure(2)
clf
subplot(1,2,1)
imagesc(img)
colormap(gray)
axis equal
axis tight


F= fftshift(fft2(img));

Fimg = abs(F);
subplot(1,2,2)
imagesc(Fimg)
colormap(gray)
axis equal
axis tight

%% FF2 of a Gabor

center = [-2,-2];
sigma  = 2;
sf = .5; %cycles/degree
theta = 45;  %angle (degrees)
phase = 90; %degrees
k = 5;  %zero padding factor

ramp = cos(theta*pi/180)*x + sin(theta*pi/180)*y;
grating = sin(2*pi*ramp*sf-phase*pi/180);

Gauss = exp( -((x-center(1)).^2 + (y-center(2)).^2)/(2*sigma^2));

Gabor = Gauss.*grating;

plotFFT2(Gabor,x,y,5,1);

%% Noise

img = randn(nPix,nPix);
plotFFT2(img,x,y);

%% Filtering by modulating amplitudes in Fourier space

%low-pass filtered noise
 Y = myfft2(img,x,y);

%cutoff all amplitudes for spatial frequencies above 2 c/deg

%Y.amp(Y.sf>2) = 0;  %sharp cutoff
 sigsf = 1;
Y.amp = Y.amp.* exp(-Y.sf.^2./(2*sigsf^2));

img2 = myifft2(Y);
plotFFT2(img2,x,y,k);

%% bandpass spatial frequency filtering

sigsf = 1;
sigCenter = 2; %c/deg

Y.amp = Y.amp.* exp(-(Y.sf-sigCenter).^2./(2*sigsf^2));
img2 = myifft2(Y);
plotFFT2(img2,x,y,k);

%% filtering angles

img = randn(nPix,nPix);
Y=myfft2(img,x,y);

%vonMises function
sigmaAng =20;  
centerAng = 0;

VonMises = exp(-sigmaAng*cos(pi*(2*(Y.angle-centerAng))/180));
Y.amp = Y.amp.*VonMises;
angleImg = myifft2(Y);

plotFFT2(angleImg,x,y,k);

%% restricted spatial frequencies and orientation

img = randn(nPix,nPix);
Y=myfft2(img,x,y);

%vonMises function
sigmaAng =20;  
centerAng = 0;

VonMises = exp(-sigmaAng*cos(pi*(2*(Y.angle-centerAng))/180));
Y.amp = Y.amp.*VonMises;

sigsf = .5;
sigCenter = 2; %c/deg

Y.amp = Y.amp.* exp(-(Y.sf-sigCenter).^2./(2*sigsf^2));


img2 = myifft2(Y);
plotFFT2(img2,x,y,k);


 




