How to compute PSD of 2D signals

Question

I have a binary matrix $X$ that filtered by a 2D filter $H$. In MATLAB, the Output $Y$ is $Y=conv2(X,H)$. Now I want to see the effect of this filtering operation by comparing the Power Spectral Density of $X$ and $Y$. Each element $x_{ij}$ of $X$ is independent and identically distributed.

My question is :

1. How can I use fft2(.) in matlab to get PSD. This tutorial, http://faculty.olin.edu/bstorey/Notes/Fourier.pdf, explains computing PSD for 1D signals but I cant find anything similar for 2D. Can someone make a reference to a good resource.
2. How can I compute PSD for this typical input/output when H=[1 1;1 1].

MATLAB Code:

f = -2*(round(rand(N,N))-0.5);    
subplot(2,3,1)
imagesc(f);colormap(gray(256));axis equal;title('Input')

Af = fft2(f);    
pf = (abs(Af)).^2;
subplot(2,3,2)
imagesc(pf);colormap(gray(256));axis equal;title('PSD of Input')

Bf = fftshift(Af); % Center FFT
Bf = abs(Bf); % Get the magnitude
Bf = log(Bf+1); % Use log, for perceptual scaling, and +1 since log(0) is undefined
Bf = mat2gray(Bf); % Use mat2gray to scale the image between 0 and 1

subplot(2,3,3)
imagesc(Bf);colormap(gray(256));axis equal;title('Shifted PSD of Input')

%--------------------------------------------------------------------------

h = [1 1;1 1];
y= conv2(f,h);

subplot(2,3,4)
imagesc(y);colormap(gray(256));axis equal;title('Output')

Cf = fft2(y);    
yf = (abs(Cf)).^2;
subplot(2,3,5)
imagesc(yf);colormap(gray(256));axis equal;title('PSD of Output')

Df = fftshift(Cf); % Center FFT
Df = abs(Df); % Get the magnitude
Df = log(Df+1); % Use log, for perceptual scaling, and +1 since log(0) is undefined
Df = mat2gray(Df); % Use mat2gray to scale the image between 0 and 1

subplot(2,3,6)
imagesc(Df);colormap(gray(256));axis equal;title('Shifted PSD of Output')