最小二乘法求解3D-2D点投影矩阵遇简单平移难题

云上老周

2026-5-9

问题分析与解决方案

咱们来拆解你遇到的两个核心问题：为什么简单平移场景下非线性最小二乘估计投影矩阵失败，以及线性方法在分段变换场景下失效的原因，再给出对应的解决思路和代码示例。

一、简单平移场景：非线性最小二乘失效的核心原因

你的虚拟数据是把3D点的X坐标直接加3得到2D点，Y坐标完全复制3D点的Y，完全忽略了3D点的Z维度——但投影矩阵描述的是透视投影（或正交投影）变换，这类变换必然会引入Z维度的影响（透视投影中Z会影响X/Y的缩放）。你的数据和投影矩阵的模型本质不匹配，这才是优化无法收敛到正确结果的核心原因。

另外，你的初始猜测initial_guess_P是通过随机的内参K、外参R/T组合出来的，和真实需要的变换（X平移3，Y不变）差距极大，也会导致优化陷入局部最优或者无法收敛。

修正方案

如果一定要用投影矩阵来拟合，你需要先让数据符合投影模型，再进行优化：

用一个已知的投影矩阵生成真实的2D点（而不是简单平移）
调整初始猜测，让它更接近真实的投影矩阵

下面是修正后的代码示例：

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

# 构造一个真实的投影矩阵P_true
K = np.array([[500, 0, 320],
              [0, 500, 240],
              [0, 0, 1]])
R = np.eye(3)
T = np.array([0, 0, 1000])  # 相机在Z轴1000位置
P_true = np.hstack([K @ R, K @ T.reshape(-1,1)])

# 生成3D点并通过真实P生成2D点
points_3d = np.random.randint(100, 400, size=(50, 3))
points_3d_hom = np.hstack([points_3d, np.ones((50,1))])
points_2d_hom = P_true @ points_3d_hom.T
points_2d = (points_2d_hom[:2]/points_2d_hom[2]).T

# 投影函数（和你的保持一致）
def projection(P, points_3d):
    p_1 = P[0, :]
    p_2 = P[1, :]
    p_3 = P[2, :]
    projected_points_2d = []
    for n in range(points_3d.shape[0]):
        points = points_3d[n, :]
        if points.shape[0] == 3:
            points = np.concatenate((points, np.array([1])))
        x = np.sum(p_1*points) / np.sum(p_3*points)
        y = np.sum(p_2*points) / np.sum(p_3*points)
        projected_points_2d.append([x, y])
    return np.array(projected_points_2d)

# 目标函数
def objective_func(x, pts2d, pts3d):
    P = np.concatenate([x, [1]]).reshape(3,4)
    proj = projection(P, pts3d)
    return (proj - pts2d).flatten()

# 优化：用真实P的前11个元素作为初始猜测
initial_guess = P_true.flatten()[:11]
ls = least_squares(objective_func, initial_guess, args=(points_2d, points_3d), 
                   method='lm', verbose=2, max_nfev=50000)
P_est = np.concatenate([ls.x, [1]]).reshape(3,4)

# 评估结果
proj_est = projection(P_est, points_3d)
residual = np.sum(np.hypot(proj_est[:,0]-points_2d[:,0], proj_est[:,1]-points_2d[:,1]))
print(f"Residual: {residual:.2f}")

# 可视化
plt.scatter(points_2d[:,0], points_2d[:,1], c='red', label='Actual')
plt.scatter(proj_est[:,0], proj_est[:,1], c='green', marker='+', label='Estimated')
plt.legend()
plt.show()

这个示例中，数据是真实透视投影的结果，初始猜测接近真实解，优化就能收敛到正确的投影矩阵。

如果你的需求只是拟合X平移3、Y不变的简单变换，完全不需要用投影矩阵——直接用2D仿射变换（甚至更简单的平移）即可，参数更少，拟合更准确。

二、分段非线性变换场景：线性方法失效的原因

线性最小二乘（包括np.linalg.lstsq）只能估计全局线性变换，而你构造的是分段线性变换（X<250和X>250用不同的缩放平移），这本质是非线性模型，线性方法自然无法拟合。

解决方案

针对这类分段变换，有两种常用思路：

1. 分段拟合（最简单直接）

先根据X的阈值把数据分成两组，分别用线性方法拟合每组的变换，然后合并结果：

import numpy as np
import matplotlib.pyplot as plt

# 生成分段变换的数据
translate = np.random.randint(20,50)
translate2 = np.random.randint(0,10)
scale = np.random.random()
scale2 = np.random.random()*2
points_3d = np.random.randint(500, size=(50,3))

# 生成2D点
x_3d = points_3d[:,0]
y_3d = points_3d[:,1]
x_2d = np.where(x_3d <250, x_3d*scale + translate, x_3d*scale2 + translate2)
y_2d = y_3d*scale + translate
points_2d = np.column_stack([x_2d, y_2d])

# 分段拟合
mask = x_3d <250
# 第一组：X<250
X1 = np.hstack([points_3d[mask,:1], np.ones((mask.sum(),1))])
Y1 = points_2d[mask,:1]
A1, _, _, _ = np.linalg.lstsq(X1, Y1, rcond=None)
# 第二组：X>=250
X2 = np.hstack([points_3d[~mask,:1], np.ones((~mask.sum(),1))])
Y2 = points_2d[~mask,:1]
A2, _, _, _ = np.linalg.lstsq(X2, Y2, rcond=None)
# Y方向统一拟合
Xy = np.hstack([points_3d[:,1:2], np.ones((50,1))])
Yy = points_2d[:,1:2]
Ay, _, _, _ = np.linalg.lstsq(Xy, Yy, rcond=None)

# 预测函数
def predict_transform(points_3d):
    x_3d = points_3d[:,0]
    y_3d = points_3d[:,1]
    x_pred = np.where(x_3d <250, x_3d*A1[0,0]+A1[1,0], x_3d*A2[0,0]+A2[1,0])
    y_pred = y_3d*Ay[0,0]+Ay[1,0]
    return np.column_stack([x_pred, y_pred])

# 可视化
pred_pts = predict_transform(points_3d)
plt.scatter(points_2d[:,0], points_2d[:,1], c='red', label='Actual')
plt.scatter(pred_pts[:,0], pred_pts[:,1], c='green', marker='+', label='Predicted')
plt.legend()
plt.show()

2. 非线性全局优化

如果想把分段参数（包括阈值）作为变量一起优化，可以用scipy.optimize.least_squares定义分段的目标函数：

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

# 生成数据（和之前一致）
translate = np.random.randint(20,50)
translate2 = np.random.randint(0,10)
scale = np.random.random()
scale2 = np.random.random()*2
points_3d = np.random.randint(500, size=(50,3))
x_3d = points_3d[:,0]
y_3d = points_3d[:,1]
x_2d = np.where(x_3d <250, x_3d*scale + translate, x_3d*scale2 + translate2)
y_2d = y_3d*scale + translate
points_2d = np.column_stack([x_2d, y_2d])

# 目标函数：参数是[scale, translate, scale2, translate2, y_scale, y_translate, threshold]
def objective_func(params, pts3d, pts2d):
    s1, t1, s2, t2, sy, ty, thresh = params
    x_3d = pts3d[:,0]
    y_3d = pts3d[:,1]
    x_pred = np.where(x_3d < thresh, x_3d*s1 + t1, x_3d*s2 + t2)
    y_pred = y_3d*sy + ty
    diff = np.column_stack([x_pred, y_pred]) - pts2d
    return diff.flatten()

# 初始猜测
initial_guess = [scale+0.1, translate+5, scale2+0.1, translate2+2, scale+0.1, translate+5, 250]
ls = least_squares(objective_func, initial_guess, args=(points_3d, points_2d), 
                   method='lm', verbose=2)

# 预测
params = ls.x
pred_x = np.where(x_3d < params[6], x_3d*params[0]+params[1], x_3d*params[2]+params[3])
pred_y = y_3d*params[4]+params[5]
pred_pts = np.column_stack([pred_x, pred_y])

# 可视化
plt.scatter(points_2d[:,0], points_2d[:,1], c='red', label='Actual')
plt.scatter(pred_pts[:,0], pred_pts[:,1], c='green', marker='+', label='Predicted')
plt.legend()
plt.show()

这种方法可以连阈值一起优化，适合不知道确切分段阈值的场景。

内容的提问来源于stack exchange，提问作者Boris Mocialov

最小二乘法求解3D-2D点投影矩阵遇简单平移难题

一、简单平移场景：非线性最小二乘失效的核心原因 #

修正方案 #

二、分段非线性变换场景：线性方法失效的原因 #