阵列位置度的最佳拟合

今天跟朋友聊天，对于孔阵列的位置度问题。孔阵列不全基准，看基准约束的自由度。理论基础见下，聊着聊着也就勾起我想写出来的欲望。

# -*- coding: utf-8 -*-

import numpy as np

from scipy.optimize import minimize

# 实际测量位置

act = np.array([

    [-14.8454, 14.8453],

    [-14.8509, -14.8321],

    [14.8395, -14.8552],

    [14.8550, 14.8437]

])

# 名义位置

nom = np.array([

    [-14.8492, 14.8492],

    [-14.8492, -14.8492],

    [14.8492, -14.8492],

    [14.8492, 14.8492]

])

# 计算每个圆的位置度误差

def calc_position_errors(act, nom):

    errors = 2*np.sqrt(np.sum(np.square(act - nom), axis=1))

    return errors

# 计算并输出没有最佳拟合的位置度误差

position_errors = calc_position_errors(act, nom)

print("Initial position dev:")

for i, error in enumerate(position_errors):

    print(f"Position {i+1}: {round(error, 4)}")

# 定义旋转函数

def rotate(points, angle):

    rotation_matrix = np.array([[np.cos(angle), -np.sin(angle)],

                                [np.sin(angle), np.cos(angle)]])

    return np.dot(points, rotation_matrix.T)

# 定义目标函数（最小化位置度误差）

def objective(params, act, nom):

    angle, tx, ty = params  

    # 旋转

    transformed_positions = rotate(act, angle)

    # 平移

    transformed_positions[:, 0] += tx

    transformed_positions[:, 1] += ty

   

    # 计算位置度误差

    distances = calc_position_errors(transformed_positions, nom)

   

    # 返回误差平方和

    return np.sum(np.square(distances))

# 初始猜测

initial_guess = [0.0, 0.0, 0.0]

# 优化

result = minimize(objective, initial_guess, args=(act, nom), method='SLSQP')

# 输出优化结果

bestfit_angle, tx, ty = result.x

print(f"\nBest_fit angle: {round(bestfit_angle*180/np.pi,4)}")

print(f"tx, ty: ({round(tx,4)}, {round(ty,4)})")

print(f"Minimum error: {result.fun}")

# 应用最优旋转和平移

trans_positions = rotate(act, bestfit_angle)

trans_positions[:, 0] += tx

trans_positions[:, 1] += ty

# 计算旋转和平移后的误差

bestfit_position_dev = calc_position_errors(trans_positions, nom)

# 输出旋转和平移后的误差

print("\nTransformed position errors:")

for i, error in enumerate(bestfit_position_dev):

    print(f"Position {i+1}: {round(error, 4)}")

# 输出总误差（均方误差）

bestfit_mean_squared_error = np.mean(np.square(bestfit_position_dev))

print(f"Best_fit mean squared error: {bestfit_mean_squared_error}")

对比CALYPSO，相差结果也是满意的。

共勉吧！