function [best_magnitude, conj_magnitude, best_angle] = cc_rotations(image,reference,angle_inc)
% function [ best_magnitude, conj_magnitude, best_angle] = cc_rotations(image,reference,angle_inc)
%
% conjugate_magnitude consists of the magnitude found at angles
% 180° opposed to the best_angle
%
% $Id: cc_rotations.m,v 1.7 2002/07/15 18:02:08 tobin Exp tobin $

image_f = fft2(image);
rotations = floor(180/angle_inc);
interpolation_algorithm = 'bicubic';

result_A = zeros(size(image));
result_B = zeros(size(image));

disp(sprintf('  There will be %d x 2 rotations',rotations));

start_time = cputime;

% I haven't yet figured out why this is necessary.
imsize = size(image);
[X,Y] = meshgrid(1:imsize(2),1:imsize(1));
X = (X ./ (imsize(1)/2)) - 1;
Y = (Y ./ (imsize(2)/2)) - 1;
phase_change_matrix = exp(2*pi*i*(X + Y));


for rotation=1:rotations
    angle = (rotation - 1) * angle_inc;
    disp(sprintf('    Processing rotation number %d (%03.2f° and %03.2f°)', ...
		  rotation,angle, angle + 180));
    
    % this code should probably be put into another function
    % since a lot is repeated here
    
    ref_rotated = imrotate(reference,angle,interpolation_algorithm);
    ref_rotated_big = pad(ref_rotated,size(image),0);
    ref_rotated_f = fft2(ref_rotated_big);
    result_f = image_f .* ref_rotated_f;
    result_A = real(fftshift(ifft2(result_f)));
    
    % rotations by 180 degrees can be done trivially, even in
    % sampled k-space (fourier-space), so we actually save some
    % time by doing it this way.  <-- unexpected bonus
    
    % Just to be sure, we'll rotate anyway... for now.
    
    ref_rotated = imrotate(reference,angle + 180,interpolation_algorithm);
    ref_rotated_big = pad(ref_rotated,size(image),0);
    ref_rotated_f = fft2(ref_rotated_big);
      
    % ref_rotated_f = rot90(ref_rotated_f,2) .* phase_change_matrix;
    
    result_f = image_f .* ref_rotated_f;
    result_B = real(fftshift(ifft2(result_f)));

    % now we need to combine the results.  the most obvious way to
    % do this is to try one result first and then do the other.

    if rotation == 1 
      best_angle = ones(size(result_A)) * angle;
      best_magnitude = result_A;
      conj_magnitude = result_B;
    end
    
    best_angle(result_A > best_magnitude) = angle;
    conj_magnitude(result_A > best_magnitude) = result_B(result_A > best_magnitude);
    best_magnitude = max(best_magnitude,result_A);

    best_angle(result_B > best_magnitude) = angle + 180;
    conj_magnitude(result_B > best_magnitude) = result_A(result_B > best_magnitude);
    best_magnitude = max(best_magnitude,result_B);
end


stop_time = cputime;
elapsed_time = stop_time - start_time;

disp(sprintf('Used %d seconds CPU time (%0.2f seconds per rotation)',elapsed_time,elapsed_time/rotations));


