结合RANSAC与最小二乘的直线拟合，python

实现效果：

 实现思路：

1、寻找内点集合：（1）从原始数据中随机选取两组点计算直线方程，（2）将其余点代入该方程并计算残差，取残差小于阈值T的点集作为“内点”，记录这些“内点”，（3）重复前两步，直到达到循环结束次数，得到很多组“内点”集合；

2、拟合直线：（1）找到上一步中“内点”数目最多的组合，利用这部分点进行最小二乘直线拟合

源代码&测试数据：

numpy == 1.21.4

matplotlib == 3.5.0

import numpy as np
import random
import matplotlib.pyplot as plt


# 用于正常显示中文
plt.rcParams['font.sans-serif'] = 'SimHei'
#用于正常显示符号
plt.rcParams['axes.unicode_minus'] = False
# 使用ggplot的绘图风格
plt.style.use('ggplot')


def LeastSquareWeightGetLine(x, y, w):
    """带有加权的直线拟合"""
    # 计算必要参数
    swxy = np.sum(w*x*y)
    swxx = np.sum(w*x*x)
    swx = np.sum(w*x)
    swy = np.sum(w*y)
    sw = np.sum(w)
    # 计算a, b
    a = (swxy - swx*swy/sw)/(swxx-swx*swx/sw)
    b = (swy-a*swx)/sw
    return [a, b]


def RANSACGetLine(x, y, T, miter=1000):
    """RANSAC方法拟合直线"""
    pnum = len(x)
    inpnums = []  # 记录内点数目
    inppos = []  # 记录内点索引
    for i in range(miter):
        # 随机选择两个点
        i1 = random.randint(0, pnum-1)
        i2 = random.randint(0, pnum-1)
        x1 = x[i1]
        y1 = y[i1]
        x2 = x[i2]
        y2 = y[i2]
        # 拟合直线
        k = (y2-y1) / (x2-x1+1e-14)
        b = y1 - k*x1
        # 利用直线计算y值
        yc = np.zeros(pnum, dtype=np.float32)
        for j in range(pnum):
            yc[j] = k*x[j] + b
        # 确定内点数目
        dy = np.abs(yc-y)
        label = np.zeros(pnum, dtype=np.int16)
        label[dy <= T] = 1
        inpnum = np.sum(label)
        # 记录内点数目
        inpnums.append(inpnum)
        # 记录内点索引
        inppos.append(np.where(label == 1))
    # 找到内点数目最多的组合，利用这部分内点进行最小二乘直线拟合
    mposindex = inpnums.index(max(inpnums))
    mpos = np.int32(inppos[mposindex][0])
    # print(mpos)
    inx = [x[i] for i in mpos]
    iny = [y[i] for i in mpos]
    w = np.ones(len(inx), dtype=np.float32)
    w = w/np.sum(w)
    return LeastSquareWeightGetLine(inx, iny, w)


if __name__ == '__main__':
    # 赋值
    # x = np.loadtxt('x.txt')
    # y = np.loadtxt('y.txt')
    x = np.array([0 ,1 ,2 ,3 ,4 ,5 ,6 ,7 ,8 ,9 ,10 ,11 ,12 ,13 ,14 ,15 ,16 ,17 ,
    18 ,19 ,20 ,21 ,22 ,23 ,24 ,25 ,26 ,27 ,28 ,29 ,30 ,31 ,32 ,33 ,34 ,
    35 ,36 ,37 ,38 ,39 ,40 ,41 ,42 ,43 ,44 ,45 ,46 ,47 ,48 ,49 ,50 ,51 ,
    52 ,53 ,54 ,55 ,56 ,57 ,58 ,59 ,60 ,61 ,62 ,63 ,64 ,65 ,66 ,67 ,68 ,
    69 ,70 ,71 ,72 ,73 ,74 ,75 ,76 ,77 ,78 ,79 ,80 ,81 ,82 ,83 ,84 ,85 ,
    86 ,87 ,88 ,89 ,90 ,91 ,92 ,93 ,94 ,95 ,96 ,97 ,98 ,99 ,100 ,101 ,
    102 ,103 ,104 ,105 ,106 ,107 ,108 ,109 ,110 ,111 ,112 ,113 ,114 ,
    115 ,116 ,117 ,118 ,119 ,120 ,121 ,122 ,123 ,124 ,125 ,126 ,127 ,
    128 ,129 ,130 ,131 ,132 ,133 ,134 ,135 ,136 ,137 ,138 ,139 ,140 ,
    141 ,142 ,143 ,144 ,145 ,146 ,147 ,148 ,149 ,150 ,151 ,152 ,153 ,
    154 ,155 ,156 ,157 ,158 ,159 ,160 ,161 ,162 ,163 ,164 ,165 ,166 ,
    167 ,168 ,169 ,170 ,171 ,172 ,173 ,174 ,175 ,176 ,177 ,178 ,179 ,
    180 ,181 ,182 ,183 ,184 ,185 ,186 ,187 ,188 ,189 ,190 ,191 ,192 ,
    193 ,194 ,195 ,196 ,197 ,198 ,199 ,200 ,201 ,202 ,203 ,204 ,205 ,
    206 ,207 ,208 ,209 ,210 ,211 ,212 ,213 ,214 ,215 ,216 ,217 ,218 ,
    219 ,220 ,221 ,222 ,223 ,224 ,225 ,226 ,227 ,228 ,229 ,230 ,231 ,
    232 ,233 ,234 ,235 ,236 ,237 ,238 ,239 ,240 ,241 ,242 ,243 ,244 ,
    245 ,246 ,247 ,248 ,249 ,250 ,251 ,252 ,253 ,254 ,255 ,256 ,257 ,
    258 ,259 ,260 ,261 ,262 ,263 ,264 ,265 ,266 ,267 ,268 ,269 ,270 ,
    271 ,272 ,273 ,274 ,275 ,276 ,277 ,278 ,279 ,280 ,281 ,282 ,283 ,
    284 ,285 ,286 ,287 ,288 ,289 ,290 ,291 ,292 ,293 ,294 ,295 ,296 ,
    297 ,298 ,299 ,300 ,301 ,302 ,303 ,304 ,305 ,306 ,307 ,308 ,309 ,
    310 ,311 ,312 ,313 ,314 ,315 ,316 ,317 ,318 ,319 ,320 ,321 ,322 ,
    323 ,324 ,325 ,326 ,327 ,328 ,329 ,330 ,331 ,332 ,333 ,334 ,335 ,
    336 ,337 ,338 ,339 ,340 ,341 ,342 ,343 ,344 ,345 ,346 ,347 ,348 ,
    349 ,350 ,351 ,352 ,353 ,354 ,355 ,356 ,357 ,358 ,359 ,360 ,361 ,
    362 ,363 ,364 ,365 ,366 ,367])
    y = np.array([-1.455 ,-2.342 ,-0.383 ,1.147 ,-1.323 ,-0.835 ,-0.889 ,-0.272 ,
    -2.109 ,-4.361 ,-1.635 ,0.122 ,-2.726 ,-3.524 ,-2.495 ,-0.765 ,-1.329 ,
    -3.054 ,-3.963 ,-4.808 ,-4.302 ,-3.015 ,-3.945 ,-2.923 ,-0.49 ,0.948 ,
    0.649 ,0.912 ,0.851 ,1.604 ,1.312 ,0.589 ,0.889 ,0.751 ,-0.3 ,0.397 ,
    -1.724 ,0.259 ,-0.491 ,0.558 ,-0.752 ,-0.821 ,-0.875 ,-1.378 ,-2.486 ,
    -1.431 ,-1.626 ,-1.472 ,-1.845 ,-3.174 ,-2.494 ,-2.342 ,-2.97 ,-1.11 ,
    -2.558 ,-3.27 ,-3.214 ,-2.799 ,-4.149 ,-4.085 ,-4.584 ,-4.504 ,-4.201 ,
    -4.602 ,-4.476 ,-5.419 ,-5.312 ,-5.468 ,-5.533 ,-6.138 ,-6.222 ,-6.223 ,
    -6.564 ,-6.662 ,-6.908 ,-6.845 ,-7.174 ,-7.236 ,-7.246 ,-7.734 ,-7.571 ,
    -7.955 ,-7.964 ,-8.267 ,-8.34 ,-8.249 ,-8.964 ,-9.053 ,-9.376 ,-9.554 ,
    -9.891 ,-9.638 ,-9.839 ,-10.143 ,-10.374 ,-10.226 ,-10.656 ,-10.84 ,
    -10.915 ,-11.29 ,-11.033 ,-11.74 ,-11.738 ,-11.715 ,-11.505 ,-12.305 ,
    -12.275 ,-12.704 ,-12.554 ,-12.984 ,-13.045 ,-13.362 ,-13.583 ,-13.749 ,
    -13.794 ,-13.94 ,-14.092 ,-14.55 ,-14.515 ,-14.937 ,-15.059 ,-14.667 ,
    -15.068 ,-15.631 ,-15.615 ,-16.129 ,-16.006 ,-16.381 ,-16.42 ,-16.567 ,
    -16.237 ,-16.798 ,-17.278 ,-17.227 ,-17.564 ,-17.867 ,-17.694 ,-18.194 ,
    -18.139 ,-18.286 ,-18.505 ,-18.657 ,-18.947 ,-19.112 ,-19.15 ,-19.554 ,
    -19.728 ,-19.728 ,-19.782 ,-20.159 ,-20.336 ,-20.637 ,-20.706 ,-21.061 ,
    -21.131 ,-21.235 ,-21.288 ,-21.769 ,-21.773 ,-21.959 ,-22.145 ,-22.14 ,
    -22.524 ,-22.686 ,-22.82 ,-23.053 ,-23.08 ,-23.389 ,-23.563 ,-23.794 ,
    -23.938 ,-24.077 ,-24.319 ,-24.453 ,-24.578 ,-24.815 ,-24.987 ,-25.083 ,
    -25.6 ,-25.557 ,-25.736 ,-25.936 ,-26.079 ,-26.273 ,-26.415 ,-26.56 ,
    -26.762 ,-27.027 ,-27.012 ,-27.374 ,-27.477 ,-27.767 ,-27.899 ,-27.875 ,
    -30.473 ,-28.705 ,-28.567 ,-28.734 ,-28.987 ,-29.154 ,-29.343 ,-29.342 ,
    -29.66 ,-29.869 ,-29.909 ,-30.211 ,-30.971 ,-30.611 ,-30.666 ,-30.89 ,
    -30.976 ,-31.371 ,-31.391 ,-31.628 ,-31.853 ,-32.141 ,-32.15 ,-32.329 ,
    -32.561 ,-32.822 ,-32.87 ,-33.14 ,-32.562 ,-33.482 ,-33.491 ,-33.882 ,
    -33.927 ,-34.322 ,-34.236 ,-34.581 ,-34.685 ,-34.616 ,-35.221 ,-35.178 ,
    -35.699 ,-36.29 ,-35.75 ,-35.896 ,-36.063 ,-36.282 ,-36.5 ,-36.655 ,
    -37.656 ,-37.058 ,-36.88 ,-37.363 ,-37.597 ,-37.723 ,-37.877 ,-38.277 ,
    -38.175 ,-38.411 ,-38.553 ,-38.81 ,-39.148 ,-39.09 ,-39.238 ,-39.367 ,
    -39.991 ,-39.938 ,-40.859 ,-40.273 ,-40.356 ,-40.768 ,-40.668 ,-41.778 ,
    -41.011 ,-41.328 ,-41.537 ,-41.578 ,-41.622 ,-42.055 ,-42.719 ,-42.473 ,
    -42.661 ,-42.678 ,-42.658 ,-43.005 ,-42.992 ,-43.81 ,-43.699 ,-43.72 ,
    -44.038 ,-44.71 ,-44.071 ,-44.75 ,-44.734 ,-44.872 ,-45.013 ,-45.275 ,
    -45.474 ,-45.808 ,-45.857 ,-46.005 ,-46.402 ,-45.615 ,-46.514 ,-47.449 ,
    -47.105 ,-47.033 ,-47.255 ,-47.47 ,-47.724 ,-47.794 ,-47.801 ,-47.77 ,
    -48.349 ,-48.36 ,-49.139 ,-47.765 ,-48.759 ,-49.704 ,-49.837 ,-49.323 ,
    -49.457 ,-49.728 ,-49.506 ,-50.965 ,-51.158 ,-50.685 ,-49.75 ,-51.065 ,
    -51.678 ,-51.579 ,-52.76 ,-50.725 ,-51.736 ,-51.861 ,-49.893 ,-52.627 ,
    -52.855 ,-52.91 ,-53.136 ,-53.266 ,-53.809 ,-53.565 ,-53.674 ,-54.797 ,
    -54.274 ,-55.306 ,-54.422 ,-53.082 ,-54.786 ,-55.174 ,-55.015 ,-58.145 ,
    -55.073 ,-52.608 ,-52.494 ,-49.882 ,-49.605 ,-50.287 ,-50.707 ,-50.961 ,
    -49.62 ,-48.789 ,-49.156 ,-49.695 ,-52.485 ,-51.045 ,-52.583 ,-54.236 ,
    -56.694 ,-58.868 ,-59.74 ,-57.928 ,-59.093 ,-58.384])
    y[0:int(len(y)/4)] = y[0:int(len(y)/4)] - 3.
    y[int(2*len(y)/4):int(3*len(y)/4)] = y[int(2*len(y)/4):int(3*len(y)/4)] + 3.
    # 拟合直线
    ab = RANSACGetLine(x, y, 0.5, len(x))
    # 显示
    showy = ab[0]*x + ab[1]
    plt.figure(figsize=(20, 5))
    plt.plot(x, showy, color='r', label='拟合直线')
    plt.plot(x, y, 'g.', label='散点')
    plt.title('RANSAC拟合直线测试')
    plt.legend()
    plt.show()