Random Sample Consensus (RANSAC)

This page is available for practice as an interactive jupyter notebook.

Plane determination using 3 points

If 3 points are given, plane parameters can be computed using the Hessian Normal Form (HNF).

HNF: $d=-\hat{n}p$ with $d\in\mathbb{R}$ a scalar representing the distance to the coordinate system origin, $\hat{n}\in\mathbb{R}^3$ normlized normal vector and $p\in\mathbb{R}^3$ some 3D point.

HNF can be determined from thre points $p_1,p_2,p_3 \in \mathbb{R}^3$ as follows:

  • Normal vector: $n = (p_2 - p_1)\times(p_3 - p_1)$

  • Normalization: $\hat{n}=\frac{n}{| n |}$

  • Calculat distance to origin: $d = - \hat{n}p_1$

Segmentation with normal vectors/planes

Local normal vectors and plane parameters are usefull features for such tasks as segmentation, classification, visualization and many more.

nv.png

This image shows how different regions like ground and vegetation can be distinguished by different normal vector distributions.

Segmentation with RANSAC

In the following task we will determine plane parameters from a cloud of points employing random sample consensus (RANSAC). The main idea of RANSAC is to estimate multiple sets of plane parametersbased on randomly selected sub sets and evaluate each plane on the whole set of points. The best matching plane is then kept.

Algorithm

  • Compute min. number of iterations
    • $w$: Prob. that a drawn observation belongs to the model
    • $n$: Number of necessary observations (here: 3)
    • $z$: Desired probability z with which the model shall be found
  • Randomly pick 3 points
  • Compute plane parameters from these three points
  • Count consensus
    • Define max. distance point-plane
  • Repeat k times
  • Choose plane with max. consensus set
  • Compute final plane parameters using all plane points

Task

You are given two point sets schneiderberg15cm_1.ply and schneiderberg15cm_2.ply. Goal is to fit one plane to the first scene and two planes to the second scene.

  • Detect the correct planes in scene 1 and two planes in scene 2 by using RANSAC. You have to estimate the required number of draws. You have to define the probabilities w, that a drawn observation belongs to a model, and z with which the model shall be found.

  • Choose a sufficient value for the max. point distance to a wall max_diff. Justify your decision in the report.

  • Plot the planes and outliers of both scenes

Scene 1 scene1.png

Scene 2

scene2.png

import random
#Usefull for the selection of random points:
print(random.sample(range(1, 100), 4))

[26, 92, 37, 50]



def find_random_points(samplesize, points):
    """Finds a set of unique random points from a point cloud
    :params:
        samplesize - size of the generated set
        points - set of points from which random points are selected
    
    :return:
        list of random points
    """
    ### BEGIN SOLUTION
    random_indx = random.sample(range(1, points.shape[0]), samplesize)
    random_points = []
    for ind in random_indx:
        random_points.append(points[ind,:])
    ### END SOLUTION
    return random_points

# use HNF here generate plane from 3 points
def plane_from_3_points(points):
    """Determines the normal vector of a plane given by 3 points
    
    :params:
        points - list of 3 points
    :return:
        (n,d) - Hesse normal from of a plane
                n - normalized normal vector
                d - distance from origin
    """
    ### BEGIN SOLUTION
    a = points[0]
    b = points[1]
    c = points[2]
    
    v_ab = b - a 
    v_ca = a - c
    n = np.cross(v_ab, v_ca)
    n = n/np.linalg.norm(n)
    
    d = np.inner(n, a)
    ### END SOLUTION
    return n, d
    

def count_consensus(n,d,max_diff,points):
    """Counts the consensus
       counts the consensus and gives back the
       line numbers of the corresponding points
        
        :params:
            n - normal vector of the plane
            d - distance of the plane from origin
            max_diff - max distance of a point to plane
            points - array of points
        
        :return:
            new_points - list - Consensual points
            outliers - list - Outlier points
    """
    ### BEGIN SOLUTION
    new_points = []
    outliers = []
    for p in points:
        diff = np.inner(n,p) - d
        if abs(diff) <= max_diff:
            new_points.append(p)
        else:
            outliers.append(p)
    ### END SOLUTION
    return new_points, outliers

import math

# implement RANSAC 
def ransac(points, z, w, max_dif):
    '''Finds the normal vectors of the planes in a scene
    
    :params: 
        z: likelihood, that plane can be found
        w: outliers in percent
        t: treshold for accepted observations
    :return:
        n: - normalized normal vector
        d: - distance from origin
        new_points: - list - Consensual points
        outliers: - list - Outlier points
    '''
    ### BEGIN SOLUTION
    samplesize = 3
    b = w**samplesize
    k_max = math.ceil(np.log10(1 - z)/np.log10(1-b))
    max_points = []
    for k in range(k_max):
        p_random = find_random_points(samplesize,points)
        
        n, d = plane_from_3_points(p_random)
        
        #[m_list(k)] = count_consensus(n_list(k,:),d(k),max_diff,points); 
        new_points, new_outliers = count_consensus(n,d,max_diff,points)
        if  len(new_points) > len(max_points):
            max_points = new_points
            outliers = new_outliers
            n_max = n
            d_max = d
    ### END SOLUTION
    return n_max, d_max, max_points, outliers    

from plyfile import PlyData
import numpy as np

# load data

_data = PlyData.read('schneiderberg15cm_1.ply').elements[0].data
schneider1 = np.stack((_data['x'], _data['y'], _data['z']), axis=1) 
_data = PlyData.read('schneiderberg15cm_2.ply').elements[0].data
schneider2 = np.stack((_data['x'], _data['y'], _data['z']), axis=1) 


print(schneider1)
print(schneider2)

[[ -0.29499999   3.61199999  -0.14300001]
 [-13.90999985   8.6079998    4.2670002 ]
 [-14.03199959   8.54300022   4.19799995]
 ..., 
 [-10.16800022   9.54500008  13.10000038]
 [-10.04800034   9.20199966  12.78800011]
 [  0.19599999   1.83099997   8.70300007]]
[[ -0.29499999   3.61199999  -0.14300001]
 [-13.90999985   8.6079998    4.2670002 ]
 [-14.03199959   8.54300022   4.19799995]
 ..., 
 [ -3.56599998  11.59200001   8.68299961]
 [ -1.09500003   3.76200008  16.45499992]
 [ -3.63199997  11.74100018  16.63400078]]

#Apply RANSAC to determine planes

# Scene 1
##
z = 0.99 # Wahrscheinlichkeit zum Finden der Ebene
w = 0.55 # Anzahl der Ebenenpunkte
max_diff = 0.1 # Maximaler Abstand eines Punktes zur Ebene
n, d, new_points, outliers = ransac(schneider1, z, w, max_diff)
print("Result: {}".format((len(new_points),len(outliers))))

Result: (3929, 1842)

# plot results
# you may use following example code or ploting framework you prefer
import ipyvolume.pylab as ply

ply.scatter(schneider2[:,0], schneider2[:,1], schneider2[:,2], size=0.5, color='red')
ply.scatter(schneider1[:,0], schneider1[:,1], schneider1[:,2], size=0.5, color='green')
ply.show()

VBox(children=(Figure(camera=PerspectiveCamera(fov=46.0, position=(0.0, 0.0, 2.0), quaternion=(0.0, 0.0, 0.0, …
Author: Artem Leichter
Last modified: 2019-04-30