|logsample(arr, rmin, rmax, xc, yc, nr, nw)
function logarr = logsample(arr, rmin, rmax, xc, yc, nr, nw)
% LOGSAMPLE Compute log-polar 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 log-polar 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
% log-sampling 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 (NR-1) lies at radius RMAX in
% the original image. Wedge W lies in the direction of the positive
% x-axis, and W increases clockwise for an image in which the y-axis
% 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 log-polar formulae
% For reference, the formulae relating positions in the log-polar array to
% positions in the original image are as follows. R and W are ring and
% wedge numbers in the log-polar array and X and Y are column and row
% numbers in the original array, all treated as if they could take
% non-integer 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*(NR-1)/NW )
% If this is satisfied, the spatial separation of neighbouring pixels in
% the log-polar 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, nr-1];
Ydata = [0, nw-1];
Size = [nw, nr];
logarr = imtransform(arr, t, ...
'Udata', Udata, 'Vdata', Vdata, ...
'Xdata', Xdata, 'Ydata', Ydata, 'Size', Size);