通过QR分解从仿射变换矩阵提取旋转矩阵的异常问题排查

你的代码存在三个核心问题，导致旋转提取结果不符合预期：

• 角度未转弧度：math.cos和math.sin要求输入弧度值，但你直接传入了45度，生成的旋转矩阵完全错误。
• QR分解用法错误：你通过拼接叉乘构造3x3矩阵做QR分解，但这不是从2D仿射矩阵提取旋转分量的正确方式。仿射矩阵的线性变换部分是前两列组成的2x2矩阵，应该针对这个子矩阵做正交分解。
• 向量对比逻辑错误：你用提取的Q矩阵直接乘原始点src[0]和dst[0]，但这两个点是位置点而非方向向量，变换逻辑不符合仿射变换的推导逻辑。

修正后的完整代码 #

import numpy as np
import math

np.set_printoptions(suppress=True)

def compute_affine(ins, out): # 输入Nx3点集，输出Nx2点集
    l = len(ins)
    entry = lambda r, d: np.linalg.det(np.delete(np.vstack([r, ins.T, np.ones(l)]), d, axis=0))
    M = np.array([[(-1) ** i * entry(R, i) for R in out.T] for i in range(l + 1)])
    A, t = np.hsplit(M[1:].T / (-M[0])[:, None], [l - 1])
    t = t.flatten()
    # 验证变换正确性
    print("Affine transformation matrix:\n", A)
    print("Affine transformation translation vector:\n", t)
    print("TESTING:")
    for p, P in zip(ins, out):
        image_p = np.dot(A, p) + t
        result = "[OK]" if np.allclose(image_p, P) else "[ERROR]"
        print(p, " mapped to: ", image_p, " ; expected: ", P, result)
    return A, t

def dot_product_angle(v1, v2):
    norm1 = np.linalg.norm(v1)
    norm2 = np.linalg.norm(v2)
    if norm1 == 0 or norm2 == 0:
        print("Zero magnitude vector!")
        return 0
    # 用clip避免浮点误差导致arccos参数越界
    dot_ratio = np.clip(np.dot(v1, v2) / (norm1 * norm2), -1.0, 1.0)
    return np.degrees(np.arccos(dot_ratio))

if __name__=="__main__":
    src = np.array([
        [1, 0, 0],
        [0, 0, 1],
        [1, 0, 1],
        [1, 1, 1],
    ])

    uv = np.array([
        [0.1, 0],
        [0, -0.1],
        [0.1, -0.1],
        [1, 1]
    ])

    # 修正1：角度转为弧度
    angle = math.radians(45)
    rot_matrix = np.array([[math.cos(angle), -math.sin(angle), 0 ],
                           [math.sin(angle), math.cos(angle), 0],
                           [0, 0, 1]])

    dst = np.dot(src, rot_matrix)
    # 计算仿射变换矩阵
    src_A, src_t = compute_affine(src, uv)
    dst_A, dst_t = compute_affine(dst, uv)

    # 修正2：从仿射矩阵提取旋转分量的正确方法
    def extract_rotation(affine_mat):
        # 提取仿射矩阵的线性变换部分（前两列）
        linear_part = affine_mat[:, :2]
        # SVD分解得到正交矩阵U，对应旋转/正交变换
        U, _, Vt = np.linalg.svd(linear_part)
        # 确保行列式为1（纯旋转矩阵行列式为1，镜像为-1）
        if np.linalg.det(U) < 0:
            U[:, 1] *= -1
        return U

    src_rot = extract_rotation(src_A)
    dst_rot = extract_rotation(dst_A)

    # 修正3：使用方向向量而非位置点进行对比
    # 取src中的方向向量：两点之差
    src_dir = src[2] - src[0]
    # dst中对应的方向向量（经过旋转变换）
    dst_dir = dst[2] - dst[0]

    # 应用提取的旋转矩阵到方向向量
    transformed_src = src_rot @ src_dir[:2]
    transformed_dst = dst_rot @ dst_dir[:2]

    # 计算夹角
    angle_val = dot_product_angle(transformed_src, transformed_dst)
    print(f"夹角值: {angle_val:.4f} 度")

    # 考虑浮点误差，用isclose判断是否接近0
    if not np.isclose(angle_val, 0, atol=1e-3):
        raise Exception(f"夹角不符合预期，应为0，实际得到 {angle_val}")
    else:
        print("夹角符合预期（接近0）")

关键修复点说明 #

1. 弧度转换：将角度值转为弧度，保证旋转矩阵计算正确。
2. 旋转分量提取：使用SVD分解仿射矩阵的线性部分，得到正交旋转矩阵，同时修正行列式确保是纯旋转（无镜像），比QR分解更稳定可靠。
3. 向量选择：使用两点差作为方向向量，避免位置点的干扰，真正对比旋转后的方向一致性。
4. 数值稳定性：计算arccos时用np.clip处理浮点误差，避免参数超出[-1,1]范围导致报错。

内容的提问来源于stack exchange，提问作者иван зуйков