Python: Incorrect Rotation angle using scipy.optimize.minimize

Question

Two sets of x and y coordinates, one of which rotates around the origin (0,0) to become the other, are what I have. I want to know what the rotation angle is.

In any case, I'm seeing an incorrect angle and am unable to locate the script's mistake.

This script results in an approximate angle of -45°, although the correct angle for this example is 5°.

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

# function to calculate the error between both datasets
def error(rotation_angle, data1, data2):

    theta = np.radians(rotation_angle)
    rotation_matrix = np.array([[np.cos(theta), -np.sin(theta)],[np.sin(theta), np.cos(theta)]])
    
    # rotate data1
    rotated_data1 = np.matmul(data1, rotation_matrix)
    error = np.sum((rotated_data1 - data2) ** 2)

    return error

# calculated sample data for a rotation of 5°
data1 = np.array([[0, 0], [0, 1], [0, 2], [0, 3]])
data2 = np.array([[0, 0], [0.08715574, 0.9961947],[0.17431149, 1.9923894],[0.26146723, 2.98858409]])

# this value is low as expected (6.741376595818406e-17)
print(error(5, data1, data2))

initial_guess = 3

# minimize error to obtain the angle
result = minimize(error, initial_guess, args=(data1, data2), method='L-BFGS-B')

print(result)

fitted_rotation_angle = result.x[0]
print("rotation angle:", fitted_rotation_angle)

The result appears as follows:

       message: Optimization terminated successfully.
       success: True
        status: 0
           fun: 13.101511152940214
             x: [-4.500e+01]
           nit: 27
          nfev: 54
 final_simplex: (array([[-4.500e+01],
                       [-4.500e+01]]), array([ 1.310e+01,  1.310e+01]))
rotation angle: -45.0000000000001

On the other hand, an error-function plot appears to be in good shape:

angle = np.arange(4, 6., 0.01)
plt.plot(angle,[error(a, data1, data2) for a in angle])
plt.xlabel("rotation angle [°]")
plt.ylabel("error")
plt.show()

Error-Function over rotation angle.

What am I doing wrong?

Answer 1

The PROBLEM is that minimize is not sending a single number for rotation_angle. It's sending a one-element numpy array, which screws up your computations. Add:

def error(rotation_angle, data1, data2):
    if isinstance(rotation_angle,np.ndarray):
        rotation_angle = rotation_angle[0]

and things all work well.

Answer 2

This is an overly complex way of finding the rotation angle. There is a closed form solution.

An implementation of it can be found here: Calculate the angle between the rows of two matrices in numpy

Example, adapting to degrees:

def angle_rowwise(A, B):
    p1 = np.einsum('ij,ij->i',A,B)
    p2 = np.linalg.norm(A,axis=1)
    p3 = np.linalg.norm(B,axis=1)
    p4 = p1 / (p2*p3)
    return np.degrees(np.arccos(np.clip(p4,-1.0,1.0)))

angle_rowwise(data1, data2)

Prints:

array([       nan, 4.99999983, 5.00000012, 5.00000004])

This is about 30x faster than the optimizer.

Answer 3

You could simply find the polar angle of each vector and subtract them.

Note: rotations only defined up to modulo 360 degs; rotation of the origin is meaningless; no check on length of vectors; positive rotation is anticlockwise.

import numpy as np
from math import atan2, pi

def rotation_angle( a, b ):
    if sum( a * a ) < 1.0e-20: return 0.0
    degrees = ( atan2( b[1], b[0] ) - atan2( a[1], a[0] ) ) * 180.0 / pi
    return degrees

data1 = np.array( [ [0, 0], [0, 1], [0, 2], [0, 3] ] )
data2 = np.array( [ [0, 0], [0.08715574, 0.9961947], [0.17431149, 1.9923894], [0.26146723, 2.98858409] ] )
for a, b in zip( data1, data2 ): print( rotation_angle( a, b ) )

Output:

0.0
-4.999999833640675
-5.000000119029441
-5.000000040545375