Ransac算法例程

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

# 数据量。
SIZE = 50
# 产生数据。np.linspace 返回一个一维数组，SIZE指定数组长度。
# 数组最小值是0，最大值是10。所有元素间隔相等。
X = np.linspace(0, 10, SIZE)
Y = 3 * X + 10

fig = plt.figure()
# 画图区域分成1行1列。选择第一块区域。
ax1 = fig.add_subplot(1,1, 1)
# 标题
ax1.set_title("RANSAC")


# 让散点图的数据更加随机并且添加一些噪声。
random_x = []
random_y = []
# 添加直线随机噪声
for i in range(SIZE):
    random_x.append(X[i] + random.uniform(-0.5, 0.5)) 
    random_y.append(Y[i] + random.uniform(-0.5, 0.5)) 
# 添加随机噪声
for i in range(SIZE):
    random_x.append(random.uniform(0,10))
    random_y.append(random.uniform(10,40))
RANDOM_X = np.array(random_x) # 散点图的横轴。
RANDOM_Y = np.array(random_y) # 散点图的纵轴。

# 画散点图。
ax1.scatter(RANDOM_X, RANDOM_Y)
# 横轴名称。
ax1.set_xlabel("x")
# 纵轴名称。
ax1.set_ylabel("y")

# 使用RANSAC算法估算模型
# 迭代最大次数，每次得到更好的估计会优化iters的数值
iters = 100000
# 数据和模型之间可接受的差值
sigma = 0.25
# 最好模型的参数估计和内点数目
best_a = 0
best_b = 0
pretotal = 0
# 希望的得到正确模型的概率
P = 0.99
for i in range(iters):
    # 随机在数据中红选出两个点去求解模型
    sample_index = random.sample(range(SIZE * 2),2)
    x_1 = RANDOM_X[sample_index[0]]
    x_2 = RANDOM_X[sample_index[1]]
    y_1 = RANDOM_Y[sample_index[0]]
    y_2 = RANDOM_Y[sample_index[1]]

    # y = ax + b 求解出a，b
    a = (y_2 - y_1) / (x_2 - x_1)
    b = y_1 - a * x_1

    # 算出内点数目
    total_inlier = 0
    for index in range(SIZE * 2):
        y_estimate = a * RANDOM_X[index] + b
        if abs(y_estimate - RANDOM_Y[index]) < sigma:
            total_inlier = total_inlier + 1

    # 判断当前的模型是否比之前估算的模型好
    if total_inlier > pretotal:
        iters = math.log(1 - P) / math.log(1 - pow(total_inlier / (SIZE * 2), 2))
        pretotal = total_inlier
        best_a = a
        best_b = b

    # 判断是否当前模型已经符合超过一半的点
    if total_inlier > SIZE:
        break

# 用我们得到的最佳估计画图
Y = best_a * RANDOM_X + best_b

# 直线图
ax1.plot(RANDOM_X, Y)
text = "best_a = " + str(best_a) + "\nbest_b = " + str(best_b)
plt.text(5,10, text,
         fontdict={'size': 8, 'color': 'r'})
plt.show()
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

#include   
#include   
#include   
#include   
  
using namespace std;  
  
const int SIZE = 50;  
  
int main() {  
    // 产生数据  
    double X[SIZE];  
    double Y[SIZE];  
    for (int i = 0; i < SIZE; i++) {  
        X[i] = i * (10.0 / (SIZE - 1)); // 最小值为0，最大值为10，间隔相等  
        Y[i] = 3 * X[i] + 10;  
    }  
  
    // 添加随机噪声  
    srand(time(NULL)); // 初始化随机数种子  
    double random_x[SIZE * 2];  
    double random_y[SIZE * 2];  
    for (int i = 0; i < SIZE; i++) {  
        random_x[i] = X[i] + (rand() % 11 - 5) / 10.0; // 随机值在-0.5到0.5之间  
        random_y[i] = Y[i] + (rand() % 11 - 5) / 10.0;  
        random_x[i + SIZE] = (rand() % 11) / 10.0; // 随机值在0到1之间  
        random_y[i + SIZE] = (rand() % 31) / 10.0 + 1; // 随机值在1到4之间  
    }  
  
    // 使用RANSAC算法估算模型  
    int iters = 100000; // 迭代最大次数  
    double sigma = 0.25; // 数据和模型之间可接受的差值  
    double best_a = 0; // 最佳模型的参数估计和内点数目  
    double best_b = 0;  
    int pretotal = 0;  
    double P = 0.99; // 希望得到正确模型的概率  
    for (int i = 0; i < iters; i++) {  
        // 随机在数据中选出两个点去求解模型  
        int sample_index[2];  
        sample_index[0] = rand() % (SIZE * 2);  
        sample_index[1] = rand() % (SIZE * 2);  
        while (sample_index[1] == sample_index[0]) { // 保证两个点不同  
            sample_index[1] = rand() % (SIZE * 2);  
        }  
        double x_1 = random_x[sample_index[0]];  
        double x_2 = random_x[sample_index[1]];  
        double y_1 = random_y[sample_index[0]];  
        double y_2 = random_y[sample_index[1]];  
  
        // y = ax + b 求解出a，b  
        double a = (y_2 - y_1) / (x_2 - x_1);  
        double b = y_1 - a * x_1;  
  
        // 算出内点数目  
        int total_inlier = 0;  
        for (int index = 0; index < SIZE * 2; index++) {  
            double y_estimate = a * random_x[index] + b;  
            if (abs(y_estimate - random_y[index]) < sigma) {  
                total_inlier++;  
            }  
        }  
  
        // 判断当前的模型是否比之前估算的模型好  
        if (total_inlier > pretotal) {  
            iters = log(1 - P) / log(1 - pow(total_inlier / (SIZE * 2), 2));  
            pretotal = total_inlier;  
            best_a = a;  
            best_b = b;  
        }  
  
        // 判断是否当前模型已经符合超过一半的点  
        if (total_inlier > SIZE) {  
            break;  
        }  
    }  
  
    // 用我们得到的最佳估计画图  
    for (int i = 0; i < SIZE; i++) {  
        Y[i] = best_a * X[i] + best_b;  
    }  
  printf("best_a: %f,bast_b: %f",best_a,best_b);
}
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