function thmin = findrot(S, varargin)
%SAXS/FINDROT Find rotation to align a SAXS image with the major axes.
%   FINDROT(S) searches for and returns the degrees of rotation to maximize
%   symmetry about the major axes of SAXS image S.  If S has not been
%   centered or converted to polar coordinates, FINDROT will attempt to do
%   so before beginning the search.
%
%   FINDROT(S, 'debug') produces a polar plot of the masked image and
%   prints rotational values and corresponding objective function values.
%
%   See also:  SAXS, SAXS/CENTER, SAXS/POLAR.

debug = 0;  % no debugging output unless requested
if nargin > 1 & isequal(varargin{1}, 'debug')
    debug = 1;
end

switch S.type
    case 'xy'
        S = polar(S);
    case 'rth'
        % Good.
    otherwise
        error('SAXS object must be of type ''xy'' or ''rth''.')
end

% Determine pixel spacings.
rstep = (S.xaxis(2) - S.xaxis(1)) / (size(S.pixels, 2) - 1);
thstep = (S.yaxis(2) - S.yaxis(1)) / (size(S.pixels, 1) - 1);

% Remove far-radius columns containing NaNs.
% I hope no near-radius columns contain NaNs, or we're going to lop off
% most of our image.
nancols = sum(isnan(S.pixels));
lastcol = min(find(nancols)) - 1;
lastx = axiscoords(S, [lastcol 1]);
lastx = lastx(1);
S = crop(S, [S.xaxis(1), S.yaxis(1); lastx, S.yaxis(2)]);
S = history(S, 'Lop off far-radius NaN columns.');
    
% Zero out inner masked zone.
if S.cmask_radius > 0
    S.pixels(:, 1 : fix(S.cmask_radius / rstep)) = 0;
end

if debug
    image(S)
end

% Sample objective values.
step = 90 / 20;
th = -45 : step : +45;
for ith = 1 : length(th)
    objval(ith) = rotobs(th(ith), S, thstep, debug);
end

% Note:
% Since we shift rows in the polar image to rotate, our rotation angles are
% limited to the angular spacing in the polar image.

% Thinking:
% Images symmetric around only one axis will have one minimum objective
% value within a half-rotation.  Images symmetric around both perpendicular
% axes will have two minima 45 degrees apart.  We would like to use the
% angle that places more image intensity near the major axes, as opposed to
% the one that places more intensity between the axes.  One heuristic could
% be to use the smaller of the two angles.  Another could be to compare
% intensities along radii near the axes to intensities along radii near 45
% degrees.

% Search for local minimum around lowest sample.
% [Why not search the entire half-rotation interval?]
thlow = th(find(objval == min(objval(:))));
options = optimset('TolX', 0.01);
thmin = fminbnd(@rotobs, thlow - step, thlow + step, options, S, thstep, debug);

if debug
    if thmin>=0
        hline([thmin thmin+90 thmin+180 thmin+270]);
    else
        hline([thmin+90 thmin+180 thmin+270 thmin+360]);
    end
end


function obs = rotobs(th, S, thstep, debug)
% Calculates objective function for rotational angle
% The objective function is the sum of squared differences between top data
% and bottom data and between left data and right data at the given
% rotation.  Add the sums of squared differences of top-bottom and
% left-right.
%
% TH may be an array of angles to evaluate, whereupon ROTOBS returns a
% like-sized array of objective values in OBS.
%
% S is a SAXS object.

% Convert to an equivalent theta in the first quadrant.
if th <= 0
    th = th + 90;
end
if th > 90
    th = th - 90;
end

pix = S.pixels;
obs = 0;
for i = 1 : 2  % 1 -> x axis, 2 -> y axis
    th = th + (i - 1) * 90;  % rotate 90 degrees for y axis
    thpix = fix(th / thstep);
    th180pix = fix(180 / thstep);
    top = pix(thpix : thpix + th180pix - 1, :);  % top or left half of data
    if th == 1
        bot = flipud(pix(th180pix + 1 : end - 1, :));
    else
        bot = flipud([pix(thpix + th180pix : end - 1, :); pix(1 : thpix - 1, :)]);
    end  % bottom or right half of data at given rotation
    obs = obs + sum(sum((top - bot) .^ 2));
end
% How it was---
% for i = 1 : 2  % 1 -> x axis, 2 -> y axis
%     th = th + (i - 1) * 90;  % rotate 90 degrees for y axis
%     thpix = fix(th / thstep);
%     th180pix = fix(180 / thstep);
%     top = pix(thpix : thpix + th180pix - 1, :);  % top or left half of data
%     if th == 1
%         bot = flipud(pix(th180pix + 1 : end - 1, :));
%     else
%         bot = flipud([pix(thpix + th180pix : end - 1, :); pix(1 : thpix - 1, :)]);
%     end  % bottom or right half of data at given rotation
%     obs = obs + sum(sum((top - bot) .^ 2));
% end

if debug
    fprintf('%g, %g\n', th, obs);
end
