# Random Sample Consensus (RANSAC)

## 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.

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 $k>\frac{log(1-z)}{log(1-b)}$
• $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

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

Scene 2

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

schneider1 = np.stack((_data['x'], _data['y'], _data['z']), axis=1)
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, …