function [srad, sang, S] = specxture(f)
%SPECTURE Computes spectral texture of an image.
%   [SRAD, SANG, S] = SPECXTURE(F) computes SRAD, the spectral energy
%   distribution as a function of radius from the center of the
%   spectrum, SANG, the spectral energy distibution as a function of
%   angle for 0 to 180 degrees in increments of 1 degree, and S = 
%   log(1 + spectrum of f), normalized to the range [0, 1]. The
%   maximum value of radios is min(M,N), where M and N are the number
%   of rows and columns of image (region) f. Thus, SRAD is a row
%   vector of lenght = (min(M, N)/2) -1); and SANG is a row vector of
%   lenght 180.
% Obtain the centerd spectrum, S, of f. The variables of S are% (u,v),
% running from 1:M and 1:N, with the center (zero frequency).
% at [M/2+1,N/2+1]
S = fftshift(fft2(f));
S = abs(S);
[M,N]= size(S);
x0 = M/2 +1;
y0 = N/2 +1;

% Maximum radius that guarantees a circle centered at (x0, y0) that
% does not exceed the boundaries of S.
rmax = min(M,N)/2 -1;

% compute srad;
srad = zeros (1,rmax);
srad(1) = S(x0, y0);
for r = 2:rmax;
    [xc,yc]= halfcircle(r, x0, y0);
    srad(r) = sum(S(sub2ind(size(S), xc, yc)));
end

% Compute sang.
[xc, yc] = halfcircle(rmax, x0, y0);
sang= zeros(1,length(xc));
for a = 1: length(xc)
    [xr,yr]= radial(x0,y0, xc(a), yc(a));
    sang(a)= sum(S(sub2ind(size(S), xr, yr)));
end

% Output the log of the spectrum for easier viewing, scaled to the
% range [0,1].
S = mat2gray(log(1+S));

%--------------------------------------------------------------------%
function [xc, yc] = halfcircle(r, x0, y0)
%   Computes the integer coordinates of a half circle of radius r and
%   center at (x0, y0) using one degree increments.
%
%   Goes from 91 to 270 because we want that the half circle to be in the
%   region defined by tp right and top left quadramts , in the
%   standart image coordinates.

theta=91:270;
theta = theta*pi/180;
[xc,yc]= pol2cart(theta, r);
xc = round(xc)'+ x0; %Column vector
yc = round(yc)'+y0;

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

function [xr, yr]= radial(x0, y0, x, y);
%   Computes the coordinates of a straight line segment extending
%   from (x0, y0) to (x, y).
%
%   Based on function intline.m. xr and yr are
%   returned as column vectors.

[xr, yr] = intline(x0, x, y0, y);

