ransac

diff --git a/imago_pack/ransac.py b/imago_pack/ransac.py
new file mode 100644 (file)
index 0000000..3f8bb65
--- /dev/null
+++ b/imago_pack/ransac.py
@@ -0,0 +1,51 @@
+"""RANSAC estimation."""
+
+import random
+from math import sqrt
+import numpy as NP
+
+def initial_estimate(data):
+    return random.sample(data, 2)
+
+def points_to_line((x1, y1), (x2, y2)):
+    return (y2 - y1, x1 - x2, x2 * y1 - x1 * y2)
+
+def filter_near(data, line, distance):
+    a, b, c = line
+    dst = lambda (x, y): abs(a * x + b * y + c) / sqrt(a*a+b*b)
+    is_near = lambda p: dst(p) <= distance
+    return [p for p in data if is_near(p)]
+
+def least_squares(data):
+    x = NP.matrix([(a, 1) for (a, b) in data])
+    xt = NP.transpose(x)
+    y = NP.matrix([[b] for (a, b) in data])
+    [a,c] = NP.dot(NP.linalg.inv(NP.dot(xt, x)), xt).dot(y).flat
+    return (a, -1, c)
+
+def get_model(data):
+    if len(data) == 2:
+        return points_to_line(*data)
+    else:
+        return least_squares(data)
+
+def iterate(data, distance):
+    consensus = 0
+    consensual = initial_estimate(data)
+    while (len(consensual) > consensus):
+        consensus = len(consensual)
+        model = get_model(consensual)
+        consensual = filter_near(data, model, distance)
+    return consensus, model
+        
+def estimate(data, dist, k):
+    best = 0
+    model = None
+    for i in xrange(0, k):
+        new, new_model = iterate(data, dist)
+        if new > best:
+            best = new
+            model = new_model
+
+    return model
+