function FINAL=calcpsd(MATRIX,HEIGHT,WIDTH);

% Computes the DSP calculation and returns the matrix of the DSP.
%
% EXAMPLES:
%
% >> calcpsd('cell.jpg')
%    where cell.jpg is the cellular pattern image
% OR
% >> calcpsd(MATRIX)
%    where MATRIX is the double-type matrix of
%    the cellular pattern image, i.e. MATRIX=double(imread('cell.jpg')).
%
% You can also specify the dimensions of the image/matrix
% (they are the same), but it's not compulsory.
% In this case, type:
% >> calcpsd(...,HEIGHT,WIDTH)
%    where HEIGHT and WIDTH are the number of rows and columns.
%
% If you want to get the PSD matrix and call it, say, FINAL, type:
% >> FINAL=calcpsd(...);
%
% Don't forget the semi-colon, or the whole matrix will be displayed!


if isstr(MATRIX)
   MATRIX=double(imread(MATRIX)); % Turns the image into a matrix
end
% The previous loop allows us to input the image either as a function argument
% (then CELLIMAGE is a string, for instance cell.jpg),
% or the image matrix (then CELLIMAGE is a matrix).

if nargin==1 % In case the dimensions of the matrix were not input as arguments.
   [HEIGHT WIDTH]=size(MATRIX);
end

if nargin~=1 & nargin~=3
   error (' Number of expected input arguments: 1 or 3');
end


WIND=MATRIX.*fspecial('gaussian',[HEIGHT WIDTH],0.2*min(HEIGHT,WIDTH)); % Windows the image with a gaussian filter
WIND=grayscale(WIND); % Rescales the entries between 0 and 255.

FOURTRANS=fft2(WIND); % Calculates the 2D Fourier transform
ABSOLU=abs(FOURTRANS); % Calculates the modulus of the Fourier coefficients
POWSPECDEN=ABSOLU.^2;
LOGA=log(POWSPECDEN);
FINAL=fftshift(LOGA); % Brings the main frequency peak in the center od the image
FINAL=grayscale(FINAL);