logsample(arr, rmin, rmax, xc, yc, nr, nw)

function logarr = logsample(arr, rmin, rmax, xc, yc, nr, nw)
% LOGSAMPLE Compute logpolar transform of image
% LOGARRAY = LOGSAMPLE(ARRAY, RMIN, RMAX, XC, YC, NR, NW) returns an
% array of samples on a logarithmic grid.
%
% ARRAY is the initial image array.
%
% RMIN and RMAX are the radii of the innermost and outermost rings of
% the logpolar sampling pattern, in terms of pixels in the original
% image. XC and YC are the position of the centre of the pattern in
% the original image, in terms of the array indices of ARRAY.
%
% NR and NW specify the number of rings and the number of wedges in the
% logsampling pattern.
%
% Any one of RMIN, RMAX, NR or NW may be given as the empty array. In
% this case, it will be calculated from the other three using the
% "circular samples" condition (see below).
%
% LOGARR(W+1, R+1) will contain the sample value for ring R and wedge
% W. Ring 0 lies at radius RMIN and ring (NR1) lies at radius RMAX in
% the original image. Wedge W lies in the direction of the positive
% xaxis, and W increases clockwise for an image in which the yaxis
% points down the screen (as is normal). The next section gives the
% exact relationship between ring and wedge indices and position in
% terms of the original image's x and y coordinates. The imtransform
% default of bilinear interpolation is adopted, but this could be
% changed later with a resampler structure.
%
% The logpolar formulae
% 
%
% For reference, the formulae relating positions in the logpolar array to
% positions in the original image are as follows. R and W are ring and
% wedge numbers in the logpolar array and X and Y are column and row
% numbers in the original array, all treated as if they could take
% noninteger values. For a sample at (X, Y):
%
% 2 2
% Radius of sample: P = sqrt( (X  XC) + (Y  YC) )
%
% Angle of sample: T = arctan( (Y  YC) / (X  XC) )
%
% Ring number: R = K * log( P / RMIN )
%
% where K = (NR  1) / log( RMAX / RMIN )
%
% Wedge number: W = NW * T / (2 * PI)
%
% The "circular samples" condition is
%
% RMAX = RMIN * exp( 2*pi*(NR1)/NW )
%
% If this is satisfied, the spatial separation of neighbouring pixels in
% the logpolar array is approximately the same along the wedges and round
% the rings.
%
% See also LOGSAMPBACK, LOGTFORM
% Copyright David Young 2010
t = logtform(rmin, rmax, nr, nw);
nr = t.tdata.nr; % Get computed values, in case default used
nw = t.tdata.nw;
[U, V, ~] = size(arr);
Udata = [1, V]  xc;
Vdata = [1, U]  yc;
Xdata = [0, nr1];
Ydata = [0, nw1];
Size = [nw, nr];
logarr = imtransform(arr, t, ...
'Udata', Udata, 'Vdata', Vdata, ...
'Xdata', Xdata, 'Ydata', Ydata, 'Size', Size);
end

