将线段拟合到一组点

Question

我试图将一条线段拟合到一组点上，但我很难找到它的算法。我有一个二维线段L和一组二维点CL可以用任何合适的方式表示（我不在乎），比如支持向量和定义向量，两点，一个左右边界的线性方程。。。唯一重要的是直线有起点和终点，所以它不是无限的

我想将L拟合到C，这样c到L（其中c是C中的一个点）的所有距离之和都会最小化。这是一个最小二乘问题，但我（认为）不能使用多项式拟合，因为L只是一段。我在这方面的数学知识有点缺乏，因此如果有任何进一步阅读的提示，我也将不胜感激

以下是我的问题的一个例子：

橙色线应适合于蓝色点，以便每个点到该线的距离平方和最小。我不介意解决方案是否使用不同的语言或根本不使用代码，只要我可以从中提取算法

因为这更多的是一个数学问题，我不确定它是否适合SO，或者应该转移到交叉验证或数学交换

Answer 1

这里是python中的一个命题。点与线之间的距离根据此处提出的方法计算：Fit a line segment to a set of points

由于线段长度有限，因此必须使用min和max函数，或者if测试以确定我们是否必须使用垂直距离或到其中一个端点的距离，因此很难（不可能）得到解析解

因此，建议的解决方案将使用优化算法来接近最佳解决方案。它使用scipy.optimize.minimize，请参见：https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html

因为线段长度是固定的，所以我们只有三个自由度。在建议的解决方案中，我使用起始段点的x和y坐标以及段坡度作为自由参数。我使用getCoordinates函数从这3个参数和长度中获取段的起点和终点

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
import math as m
from scipy.spatial import distance

# Plot the points and the segment
def plotFunction(points,x1,x2):
    'Plotting function for plane and iterations'
    plt.plot(points[:,0],points[:,1],'ro')
    plt.plot([x1[0],x2[0]],[x1[1],x2[1]])
    plt.xlim(0, 1)
    plt.ylim(0, 1)
    plt.show()

# Get the sum of the distance between all the points and the segment
# The segment is defined by guess and length were:
# guess[0]=x coordinate of the starting point
# guess[1]=y coordinate of the starting point
# guess[2]=slope
# Since distance is always >0 no need to use root mean square values
def getDist(guess,points,length):
  start_pt=np.array([guess[0],guess[1]])
  slope=guess[2]
  [x1,x2]=getCoordinates(start_pt,slope,length)
  total_dist=0
  # Loop over each points to get the distance between the point and the segment
  for pt in points:
    total_dist+=minimum_distance(x1,x2,pt,length)

  return(total_dist)

# Return minimum distance between line segment x1-x2 and point pt
# Adapted from https://stackoverflow.com/questions/849211/shortest-distance-between-a-point-and-a-line-segment
def minimum_distance(x1, x2, pt,length):
  length2 = length**2  # i.e. |x1-x2|^2 - avoid a sqrt, we use length that we already know to avoid re-computation
  if length2 == 0.0:
    return distance.euclidean(p, v);
  # Consider the line extending the segment, parameterized as x1 + t (x2 - x1).
  # We find projection of point p onto the line.
  # It falls where t = [(pt-x1) . (x2-x1)] / |x2-x1|^2
  # We clamp t from [0,1] to handle points outside the segment vw.
  t = max(0, min(1, np.dot(pt - x1, x2 - x1) / length2));
  projection = x1 + t * (x2 - x1);  # Projection falls on the segment
  return distance.euclidean(pt, projection);


# Get coordinates of start and end point of the segment from start_pt,
# slope and length, obtained by solving slope=dy/dx, dx^2+dy^2=length
def getCoordinates(start_pt,slope,length):
    x1=start_pt
    dx=length/m.sqrt(slope**2+1)
    dy=slope*dx
    x2=start_pt+np.array([dx,dy])
    return [x1,x2]

if __name__ == '__main__':
    # Generate random points
    num_points=20
    points=np.random.rand(num_points,2)

    # Starting position
    length=0.5
    start_pt=np.array([0.25,0.5])
    slope=0

    #Use scipy.optimize, minimize to find the best start_pt and slope combination
    res = minimize(getDist, x0=[start_pt[0],start_pt[1],slope], args=(points,length), method="Nelder-Mead")

    # Retreive best parameters
    start_pt=np.array([res.x[0],res.x[1]])
    slope=res.x[2]
    [x1,x2]=getCoordinates(start_pt,slope,length)

    print("\n** The best segment found is defined by:")
    print("\t** start_pt:\t",x1)
    print("\t** end_pt:\t",x2)
    print("\t** slope:\t",slope)
    print("** The total distance is:",getDist([x1[0],x2[1],slope],points,length),"\n")

    # Plot results
    plotFunction(points,x1,x2)

Answer 2

这个解决方案与已经发布在这里的解决方案相对类似，但我认为它更高效、更优雅、更容易理解，这就是为什么尽管有相似之处，我还是发布了它

如前所述，min（max（…）公式很难解析地解决此问题，这就是为什么scipy.optimize很适合

该解决方案基于https://math.stackexchange.com/questions/330269/the-distance-from-a-point-to-a-line-segment中概述的点与有限线段之间距离的数学公式

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize, NonlinearConstraint


def calc_distance_from_point_set(v_):
    #v_ is accepted as 1d array to make easier with scipy.optimize
    #Reshape into two points
    v = (v_[:2].reshape(2, 1), v_[2:].reshape(2, 1))

    #Calculate t* for s(t*) = v_0 + t*(v_1-v_0), for the line segment w.r.t each point
    t_star_matrix = np.minimum(np.maximum(np.matmul(P-v[0].T, v[1]-v[0]) / np.linalg.norm(v[1]-v[0])**2, 0), 1)
    #Calculate s(t*)
    s_t_star_matrix = v[0]+((t_star_matrix.ravel())*(v[1]-v[0]))

    #Take distance between all points and respective point on segment
    distance_from_every_point = np.linalg.norm(P.T -s_t_star_matrix, axis=0)
    return np.sum(distance_from_every_point)

if __name__ == '__main__':

    #Random points from bounding box

    box_1 = np.random.uniform(-5, 5, 20)
    box_2 = np.random.uniform(-5, 5, 20)
    P = np.stack([box_1, box_2], axis=1)
    segment_length = 3
    segment_length_constraint = NonlinearConstraint(fun=lambda x: np.linalg.norm(np.array([x[0], x[1]]) - np.array([x[2] ,x[3]])), lb=[segment_length], ub=[segment_length])
    point = minimize(calc_distance_from_point_set, (0.0,-.0,1.0,1.0), options={'maxiter': 100, 'disp': True},constraints=segment_length_constraint).x
    plt.scatter(box_1, box_2)
    plt.plot([point[0], point[2]], [point[1], point[3]])

示例结果：