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?
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.