Hopefully would prove useful for someone!
Representation of rotation is angle only; and optimisation is steepest decent. metric is the sum of squared difference metric. Shared purely for pedagogical purposes. Python code is horribly unoptimised.
import numpy, cv2
import scipy
from scipy import ndimage
import matplotlib
from matplotlib import pyplot
#def rotate_correct(fixed, moving, init_angle=0):
DATA_DIR = 'your_data_path_here';
FILE_NAME = 'fixed.png';
ROT_FILE_NAME = 'moving.png';
INIT_ANGLE = 0;
#
fixed = cv2.imread( DATA_DIR + FILE_NAME, cv2.IMREAD_GRAYSCALE);
moving = cv2.imread(DATA_DIR + ROT_FILE_NAME, cv2.IMREAD_GRAYSCALE);
#
fixed_im = numpy.float32(fixed);
moving_im = numpy.float32(moving);
registered_im = numpy.copy(moving_im);
#
NUM_ITER = 25;
STEP_SIZE = 0.00001;
BORDER = 70;
N = registered_im[BORDER:-BORDER, BORDER:-BORDER].shape[0]*registered_im[BORDER:-BORDER, BORDER:-BORDER].shape[1]
INIT_SSD_COST = registered_im[BORDER:-BORDER, BORDER:-BORDER] - fixed_im[BORDER:-BORDER, BORDER:-BORDER];
INIT_SSD_COST = numpy.sum(INIT_SSD_COST**2);
print('Initial_SSD: {0}'.format(INIT_SSD_COST))
estimated_angle = INIT_ANGLE;
D = numpy.gradient(registered_im);
D0 = D[0];
D1 = D[1];
term1 = numpy.zeros_like(fixed_im);
#
prevSSD = INIT_SSD_COST;
prev_angle = INIT_ANGLE;
for iter in range(NUM_ITER):
term1[BORDER:-BORDER, BORDER:-BORDER] = registered_im[BORDER:-BORDER, BORDER:-BORDER] - fixed_im[BORDER:-BORDER, BORDER:-BORDER];
update = 0.0;
for x in range(BORDER,fixed_im.shape[1]-BORDER): #X direction
for y in range(BORDER,fixed_im.shape[0]-BORDER): #Y Direction
# rotation_matrix = numpy.array( [ [numpy.cos(angle), -numpy.sin(angle)], [numpy.sin(angle), numpy.cos(angle)] ]);
rotation_matrix_deriv = numpy.array( [ [-numpy.sin(estimated_angle), - numpy.cos(estimated_angle)], [numpy.cos(estimated_angle), -numpy.sin(estimated_angle)] ])
x1 = x*rotation_matrix_deriv[0][0] + y*rotation_matrix_deriv[0][1];
y1 = x*rotation_matrix_deriv[1][0] + y*rotation_matrix_deriv[1][1];
term2 = D0[y,x]*x1 + D1[y,x]*y1;
update = update + term1[y,x]*term2;
update = STEP_SIZE*update/N;
estimated_angle = estimated_angle - update;
D0 = scipy.ndimage.rotate(D[0], estimated_angle*180.0/numpy.pi,reshape=False);
D1 = scipy.ndimage.rotate(D[1], estimated_angle*180.0/numpy.pi,reshape=False);
registered_im = scipy.ndimage.rotate(moving_im, estimated_angle*180.0/numpy.pi,reshape=False);
SSD_COST = registered_im[BORDER:-BORDER, BORDER:-BORDER] - fixed_im[BORDER:-BORDER, BORDER:-BORDER];
SSD_COST = numpy.sum(SSD_COST**2);
currSSD = SSD_COST;
if currSSD > prevSSD:
estimated_angle = prev_angle;
print('Minima reached. Exiting loop');
break;
else:
prev_angle = estimated_angle;
prevSSD = currSSD;
print('Iter: {0}'.format(iter+1)),
print('Step Update: {0} Estimated Angle: {1}'.format(update*180.0/numpy.pi, estimated_angle*180.0/numpy.pi)),
print('SSD: {0}'.format(SSD_COST/INIT_SSD_COST))
if iter == NUM_ITER - 1:
print('Max iter reached. Exiting loop');
##
SSD_COST = registered_im[BORDER:-BORDER, BORDER:-BORDER] - fixed_im[BORDER:-BORDER, BORDER:-BORDER];
SSD_COST = numpy.sum(SSD_COST**2);
print('Final_SSD: {0}'.format(SSD_COST/INIT_SSD_COST))
print('Estimated_angle: {0}'.format(estimated_angle*180.0/numpy.pi))
##
DISPLAY = 0
BORDER = 1;
if DISPLAY == True:
#
pyplot.subplot(131)
pyplot.imshow(fixed_im[BORDER:-BORDER])
pyplot.gray()
pyplot.title('Before image difference')
#
pyplot.subplot(132)
pyplot.imshow(moving_im[BORDER:-BORDER])
pyplot.gray()
pyplot.title('After difference ')
#
pyplot.subplot(133)
pyplot.imshow(registered_im[BORDER:-BORDER])
pyplot.title('Registered image')
pyplot.gray()
#
pyplot.show()
