"RANSAC, a simple but highly effective algorithm to deal with outliers"
First things first, we use opencv library to graphically show the iteration processes. For the implementation of these codes, please look at
- C++
- Python
#include <iostream>
#include "opencv4/opencv2/opencv.hpp"
#include <math.h>
using namespace std;
using namespace cv;
import cv2
from random import randint
import numpy as np
from math import atan, sqrt
Following variables are used in the algorithm. Of course, they can be adjusted based on the desired output or result. Feel free to play with them. The "canvas.png" can be downloaded here.
- C++
- Python
const int nSample = 100;
const int maxSize = 1000;
int particles_x[nSample];
int particles_y[nSample];
Mat canvas_in = imread("canvas.png");
Mat canvas;
resize(canvas_in, canvas_in, Size(maxSize, maxSize));
int max_iteration = 1000;
int th_inliers = 100;
int max_inliers = -1, best_x1, best_x2, best_y1, best_y2; float best_theta;
nSample = 100
maxSize = 1000
particles_x = np.array([])
particles_y = np.array([])
canvas_in = cv2.imread("canvas.png")
canvas_in = cv2.resize(canvas_in, (maxSize, maxSize))
max_iteration = 1000
th_inliers = 100
max_inliers = -1
best_x1 = 0; best_x2 = 0; best_y1 = 0; best_y2 = 0; best_theta = 0.0
Before diving into the algorithm, random samples or particles need to be genarated to feed the algorithm. This step can be substituted to a specific problem to solve with RANSAC.
- C++
- Python
srand((unsigned int)time(NULL));
for (int i=0; i<nSample; i++) {
int rnd_x = round((float)rand()/RAND_MAX*maxSize), rnd_y = rnd_x;
while (rnd_x == rnd_y) {
rnd_y = 300 + round((float)rand()/RAND_MAX*400);
}
particles_x[i] = rnd_x;
particles_y[i] = rnd_y;
circle(canvas_in, Point(particles_x[i], particles_y[i]), 10, Scalar(128,128,128), -1);
}
for i in range(nSample):
rnd_x = randint(0, maxSize)
rnd_y = rnd_x
while rnd_x == rnd_y:
rnd_y = randint(300, 700)
particles_x = np.append(particles_x, rnd_x)
particles_y = np.append(particles_y, rnd_y)
canvas_in = cv2.circle(canvas_in, (int(particles_x[i]), int(particles_y[i])), 10, (128,128,128), -1)
Now, let's continue to the main RANSAC algorithm. There are three main steps: random sampling, scoring, and evaluation. All of them are iterated in a loop to get the best score of each iteration. The linear equations used in the algorithm are as follows.
y = mx + b
m = (y2 - y1) / (x2 - x1)
θ = atan (m)
Where, m is a slope of the line, b is a y-intercept, x and y are the values we want to solve. θ is the line angle. In the algorithm, the score is calculated based on the number of inliers.
- C++
- Python
for (int i=1; i<=max_iteration; i++) {
//step 1: random sample
int rnd_sample1 = round((float)rand()/RAND_MAX*nSample), rnd_sample2 = rnd_sample1;
while (rnd_sample1 == rnd_sample2) {
rnd_sample2 = round((float)rand()/RAND_MAX*nSample);
}
int sample1_posX = particles_x[rnd_sample1];
int sample1_posY = particles_y[rnd_sample1];
int sample2_posX = particles_x[rnd_sample2];
int sample2_posY = particles_y[rnd_sample2];
//step 2: compute the sample [y = mx+b]
float m = (float)(sample2_posY - sample1_posY)/(sample2_posX - sample1_posX);
float b = (float)sample1_posY - (m * sample1_posX);
float theta = -atan(m) * 180 / M_PI;
int new_x, new_y, cnt_inliers = 0;
canvas = canvas_in.clone();
for (int j=0; j<nSample; j++) {
if (fabs(theta) <= 45.0) { //horizontal line, find y
new_x = particles_x[j];
new_y = round((m * new_x) + b);
}
else { //vertical line, find x
new_y = particles_y[j];
new_x = round((new_y - b) / m);
}
int distance = (int)round(sqrt(((new_x-particles_x[j]) * (new_x-particles_x[j])) + ((new_y-particles_y[j]) * (new_y-particles_y[j]))));
if (distance <= th_inliers) {
circle(canvas, Point(particles_x[j], particles_y[j]), 10, Scalar(0,255,0), -1);
cnt_inliers++;
}
//line(canvas, Point(new_x, new_y), Point(particles_x[j], particles_y[j]), Scalar(0,0,0), 1);
}
//step 3: save the best score
int min_x = 0;
int max_x = maxSize;
int min_y = round((m * min_x) + b);
int max_y = round((m * max_x) + b);
if (cnt_inliers > max_inliers) {
max_inliers = cnt_inliers;
best_x1 = min_x; best_y1 = min_y;
best_x2 = max_x; best_y2 = max_y;
best_theta = theta;
}
//plot the results
line(canvas, Point(min_x, min_y), Point(max_x, max_y), Scalar(255,0,0), 5);
circle(canvas, Point(sample1_posX, sample1_posY), 10, Scalar(0,0,0), -1);
circle(canvas, Point(sample2_posX, sample2_posY), 10, Scalar(0,0,0), -1);
string theta_str = to_string((int)round(theta));
putText(canvas, "Iteration: "+to_string(i)+" score: "+to_string(cnt_inliers)+" theta: "+theta_str+"deg.", Point(10,maxSize-10), FONT_HERSHEY_COMPLEX, 1, Scalar(255,0,0), 2);
putText(canvas, "The best score: "+to_string(max_inliers), Point(10,maxSize-50), FONT_HERSHEY_COMPLEX, 1, Scalar(0,0,255), 2);
cout << "i:" << i << " rnd:[" << rnd_sample1 << " " << rnd_sample2 << "] m:" << m << " b:" << b << " thera:" << theta;
cout << " score:" << cnt_inliers << endl;
imshow("RANSAC @The JPID Coder", canvas);
waitKey(1);
}
for i in range(max_iteration):
#step 1: random sample
rnd_sample1 = randint(0, nSample-1)
rnd_sample2 = rnd_sample1
while rnd_sample1 == rnd_sample2:
rnd_sample2 = randint(0, nSample-1)
sample1_posX = particles_x[rnd_sample1]
sample1_posY = particles_y[rnd_sample1]
sample2_posX = particles_x[rnd_sample2]
sample2_posY = particles_y[rnd_sample2]
while sample1_posX == sample2_posX: #avoid m is divided by 0
rnd_sample2 = randint(0, nSample-1)
sample2_posX = particles_x[rnd_sample2]
sample2_posY = particles_y[rnd_sample2]
#step 2: compute the sample [y = mx+b]
m = (sample2_posY - sample1_posY)/(sample2_posX - sample1_posX)
b = sample1_posY - (m * sample1_posX)
theta = -atan(m) * 180 / 3.142857
new_x = 0; new_y = 0; cnt_inliers = 0
canvas = canvas_in.copy()
for j in range(nSample):
if abs(theta) <= 45.0: #horizontal line, find y
new_x = round(particles_x[j])
new_y = round((m * new_x) + b)
else: #vertical line, find x
new_y = round(particles_y[j])
new_x = round((new_y - b) / m)
distance = round(sqrt(((new_x-particles_x[j]) * (new_x-particles_x[j])) + ((new_y-particles_y[j]) * (new_y-particles_y[j]))))
if distance <= th_inliers:
canvas = cv2.circle(canvas, (int(particles_x[j]), int(particles_y[j])), 10, (0,255,0), -1)
cnt_inliers += 1
#canvas = cv2.line(canvas, (new_x,new_y), (int(particles_x[j]),int(particles_y[j])), (0,0,0), 1)
#step 3: save the best score
min_x = 0
max_x = maxSize
min_y = round((m * min_x) + b)
max_y = round((m * max_x) + b)
if cnt_inliers > max_inliers:
max_inliers = cnt_inliers
best_x1 = min_x; best_y1 = min_y
best_x2 = max_x; best_y2 = max_y
best_theta = theta
#plot the results
canvas = cv2.line(canvas, (min_x, min_y), (max_x, max_y), (255,0,0), 5)
canvas = cv2.circle(canvas, (int(sample1_posX), int(sample1_posY)), 10, (0,0,0), -1)
canvas = cv2.circle(canvas, (int(sample2_posX), int(sample2_posY)), 10, (0,0,0), -1)
canvas = cv2.putText(canvas, "Iteration: "+str(i)+" score: "+str(cnt_inliers)+" theta: "+str(round(theta))+"deg.", (10,maxSize-10), cv2.FONT_HERSHEY_SIMPLEX, 1, (255,0,0), 2)
canvas = cv2.putText(canvas, "The best score: "+str(max_inliers), (10,maxSize-50), cv2.FONT_HERSHEY_SIMPLEX, 1, (0,0,255), 2)
print("i:", i, " rnd:[", rnd_sample1, " ", rnd_sample2, "] m:", m, " b:", b, " theta:", theta, " score:", cnt_inliers)
cv2.imshow("RANSAC @The JPID Coder", canvas)
cv2.waitKey(1)
Finally, plot the final result after the iteration is complete. It is similar to Least Squares-Based Linear Regression in terms of results. However, each of them has its own advantages and disadvantages.
- C++
- Python
canvas = canvas_in.clone();
cout << endl << "The best score:" << max_inliers << " theta:" << best_theta << endl;
line(canvas, Point(best_x1, best_y1), Point(best_x2, best_y2), Scalar(0,0,255), 5);
string theta_str = to_string((int)round(best_theta));
putText(canvas, "The best score: "+to_string(max_inliers)+" theta: "+theta_str+"deg.", Point(10,maxSize-10), FONT_HERSHEY_COMPLEX, 1, Scalar(0,0,255), 2);
while (1) {
imshow("RANSAC @The JPID Coder", canvas);
waitKey(1);
}
canvas = canvas_in.copy()
print("\nThe best score:", max_inliers, " theta:", best_theta)
canvas = cv2.line(canvas, (best_x1, best_y1), (best_x2, best_y2), (0,0,255), 5)
canvas = cv2.putText(canvas, "The best score: "+str(max_inliers)+" theta: "+str(round(best_theta))+"deg.", (10,maxSize-10), cv2.FONT_HERSHEY_SIMPLEX, 1, (0,0,255), 2)
while True:
cv2.imshow("RANSAC @The JPID Coder", canvas)
cv2.waitKey(1)
Please look at for the implementation of this algorithm.
References
[1] M. A. Fischler and R. C. Bolles (1981), "Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography", doi:10.1145/358669.358692.
[2] Wikipedia, Random Sample Consensus.
Home Algorithms Applications Live Simulations YouTube About us
Comments
Post a Comment