Robust line model estimation using RANSACΒΆ

In this example we see how to robustly fit a line model to faulty data using the RANSAC algorithm.

../_images/plot_ransac_1.png

import numpy as np
from matplotlib import pyplot as plt

from skimage.measure import LineModel, ransac


np.random.seed(seed=1)

# generate coordinates of line
x = np.arange(-200, 200)
y = 0.2 * x + 20
data = np.column_stack([x, y])

# add faulty data
faulty = np.array(30 * [(180., -100)])
faulty += 5 * np.random.normal(size=faulty.shape)
data[:faulty.shape[0]] = faulty

# add gaussian noise to coordinates
noise = np.random.normal(size=data.shape)
data += 0.5 * noise
data[::2] += 5 * noise[::2]
data[::4] += 20 * noise[::4]

# fit line using all data
model = LineModel()
model.estimate(data)

# robustly fit line only using inlier data with RANSAC algorithm
model_robust, inliers = ransac(data, LineModel, min_samples=2,
                               residual_threshold=1, max_trials=1000)
outliers = inliers == False

# generate coordinates of estimated models
line_x = np.arange(-250, 250)
line_y = model.predict_y(line_x)
line_y_robust = model_robust.predict_y(line_x)

fig, ax = plt.subplots()
ax.plot(data[inliers, 0], data[inliers, 1], '.b', alpha=0.6,
        label='Inlier data')
ax.plot(data[outliers, 0], data[outliers, 1], '.r', alpha=0.6,
        label='Outlier data')
ax.plot(line_x, line_y, '-k', label='Line model from all data')
ax.plot(line_x, line_y_robust, '-b', label='Robust line model')
ax.legend(loc='lower left')
plt.show()

STDOUT


        

STDERR


        

Python source code: download (generated using skimage 0.11dev)

IPython Notebook: download (generated using skimage 0.11dev)

aW1wb3J0IG51bXB5IGFzIG5wCmZyb20gbWF0cGxvdGxpYiBpbXBvcnQgcHlwbG90IGFzIHBsdAoKZnJvbSBza2ltYWdlLm1lYXN1cmUgaW1wb3J0IExpbmVNb2RlbCwgcmFuc2FjCgoKbnAucmFuZG9tLnNlZWQoc2VlZD0xKQoKIyBnZW5lcmF0ZSBjb29yZGluYXRlcyBvZiBsaW5lCnggPSBucC5hcmFuZ2UoLTIwMCwgMjAwKQp5ID0gMC4yICogeCArIDIwCmRhdGEgPSBucC5jb2x1bW5fc3RhY2soW3gsIHldKQoKIyBhZGQgZmF1bHR5IGRhdGEKZmF1bHR5ID0gbnAuYXJyYXkoMzAgKiBbKDE4MC4sIC0xMDApXSkKZmF1bHR5ICs9IDUgKiBucC5yYW5kb20ubm9ybWFsKHNpemU9ZmF1bHR5LnNoYXBlKQpkYXRhWzpmYXVsdHkuc2hhcGVbMF1dID0gZmF1bHR5CgojIGFkZCBnYXVzc2lhbiBub2lzZSB0byBjb29yZGluYXRlcwpub2lzZSA9IG5wLnJhbmRvbS5ub3JtYWwoc2l6ZT1kYXRhLnNoYXBlKQpkYXRhICs9IDAuNSAqIG5vaXNlCmRhdGFbOjoyXSArPSA1ICogbm9pc2VbOjoyXQpkYXRhWzo6NF0gKz0gMjAgKiBub2lzZVs6OjRdCgojIGZpdCBsaW5lIHVzaW5nIGFsbCBkYXRhCm1vZGVsID0gTGluZU1vZGVsKCkKbW9kZWwuZXN0aW1hdGUoZGF0YSkKCiMgcm9idXN0bHkgZml0IGxpbmUgb25seSB1c2luZyBpbmxpZXIgZGF0YSB3aXRoIFJBTlNBQyBhbGdvcml0aG0KbW9kZWxfcm9idXN0LCBpbmxpZXJzID0gcmFuc2FjKGRhdGEsIExpbmVNb2RlbCwgbWluX3NhbXBsZXM9MiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHJlc2lkdWFsX3RocmVzaG9sZD0xLCBtYXhfdHJpYWxzPTEwMDApCm91dGxpZXJzID0gaW5saWVycyA9PSBGYWxzZQoKIyBnZW5lcmF0ZSBjb29yZGluYXRlcyBvZiBlc3RpbWF0ZWQgbW9kZWxzCmxpbmVfeCA9IG5wLmFyYW5nZSgtMjUwLCAyNTApCmxpbmVfeSA9IG1vZGVsLnByZWRpY3RfeShsaW5lX3gpCmxpbmVfeV9yb2J1c3QgPSBtb2RlbF9yb2J1c3QucHJlZGljdF95KGxpbmVfeCkKCmZpZywgYXggPSBwbHQuc3VicGxvdHMoKQpheC5wbG90KGRhdGFbaW5saWVycywgMF0sIGRhdGFbaW5saWVycywgMV0sICcuYicsIGFscGhhPTAuNiwKICAgICAgICBsYWJlbD0nSW5saWVyIGRhdGEnKQpheC5wbG90KGRhdGFbb3V0bGllcnMsIDBdLCBkYXRhW291dGxpZXJzLCAxXSwgJy5yJywgYWxwaGE9MC42LAogICAgICAgIGxhYmVsPSdPdXRsaWVyIGRhdGEnKQpheC5wbG90KGxpbmVfeCwgbGluZV95LCAnLWsnLCBsYWJlbD0nTGluZSBtb2RlbCBmcm9tIGFsbCBkYXRhJykKYXgucGxvdChsaW5lX3gsIGxpbmVfeV9yb2J1c3QsICctYicsIGxhYmVsPSdSb2J1c3QgbGluZSBtb2RlbCcpCmF4LmxlZ2VuZChsb2M9J2xvd2VyIGxlZnQnKQpwbHQuc2hvdygp