如何基于Scipy实现所有点均在直线上方的约束直线拟合？

问题描述 #

现有一张包含白色点的图像，需要对这些白色点进行直线拟合，且要求所有点均位于拟合直线的上方。目前已编写无约束条件的拟合脚本，请求协助修改代码以实现该约束。

原无约束拟合代码 #

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

def fit_line(img):

    def func_linear(x, x0, y0, k):
        # y=k(x−xo)+yo
        f = lambda x:k*(x-x0)+y0
        return f(x)

    points = np.where(img>0)
    points = np.array([points[1], points[0]]).T
    
    x = points[:,0]
    y = points[:,1]

    p0 = [1, 1, 1]
    p , e = optimize.curve_fit(func_linear, x, y, p0)

    pt1 = (np.min(x).astype(int), func_linear(np.min(x), *p).astype(int))
    pt2 = (np.max(x).astype(int), func_linear(np.max(x), *p).astype(int))

    cv2.line(img, pt1, pt2, (255,0,0), 3)


img = cv2.imread("toy_2.png")

fit_line(img)

plt.imshow(img)

修改后的约束拟合代码 #

要实现所有点在直线上方的约束，我们改用带约束的优化方法，用scipy.optimize.minimize自定义损失函数并添加不等式约束：

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

def fit_line(img):
    # 简化直线方程为y = k*x + b，更便于约束定义
    def line_func(params, x):
        k, b = params
        return k * x + b

    # 损失函数：最小化直线与点的y方向误差平方和
    def loss_func(params, x, y):
        y_pred = line_func(params, x)
        return np.sum((y - y_pred)**2)

    # 约束条件：所有点的y值 >= 直线预测的y值（确保点在直线上方）
    def constraint_func(params, x, y):
        y_pred = line_func(params, x)
        return y - y_pred  # 要求返回值全部 >= 0

    # 提取白色点的(x,y)坐标
    points = np.where(img > 0)
    points = np.array([points[1], points[0]]).T
    x = points[:, 0]
    y = points[:, 1]

    # 用无约束拟合结果作为初始参数，加快收敛
    p0_unconstrained, _ = optimize.curve_fit(lambda x, k, b: k*x + b, x, y)
    initial_guess = p0_unconstrained

    # 设置不等式约束
    constraints = {'type': 'ineq', 'fun': constraint_func, 'args': (x, y)}

    # 执行带约束的最小化优化
    result = optimize.minimize(loss_func, initial_guess, args=(x, y), constraints=constraints)
    k_opt, b_opt = result.x

    # 计算直线端点用于绘制
    x_min, x_max = np.min(x), np.max(x)
    y_min = k_opt * x_min + b_opt
    y_max = k_opt * x_max + b_opt
    pt1 = (int(x_min), int(y_min))
    pt2 = (int(x_max), int(y_max))

    # 绘制约束拟合后的直线
    cv2.line(img, pt1, pt2, (255, 0, 0), 3)

# 加载图像并执行拟合
img = cv2.imread("toy_2.png")
fit_line(img)
plt.imshow(img)
plt.show()

代码说明 #

• 替换curve_fit为optimize.minimize，因为前者不支持直接添加不等式约束；
• 约束条件明确要求每个点的y坐标大于等于直线对应x位置的y值，确保所有点在直线上方；
• 用无约束拟合结果作为初始参数，提升优化的收敛效率；
• 损失函数选用y方向误差平方和，在满足约束的前提下让直线尽可能贴近点集。

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