Zhiyuan LIU

2. Kernel algorithm

2.1 Start

The essential of RANSAC (Random SAmple Consensus) is to choose some samples randomly to construct a model. If the model satisfies the conditions, it will be kept, otherwise it will be dropped. This process will be repeated until the exit condition is satisfied. In this example, the model to construct is a plane and the samples are 3 points.

First, choose a batch of points randomly among all of the samples. This is to reduce the amount of calculation after in the iteration. A parameter should be given by user to define the number of the points in each batch. This parameter has influence on the time and the quality of the calculation.

// pick a batch of n vertices randomly from the vertex slice
vertexBatch := make([]plyfile.VertexMono, 0)
indexBatch := make([]int, 0)
for {
	indexRandom := r.Intn(len(vlist))

	// check if this vertex is already chosen. If not, append the vertex and its index to the batch
	if mymath.ExistIntList(indexBatch, indexRandom) {...}
	vertexBatch = append(vertexBatch, vlist[indexRandom])
	indexBatch = append(indexBatch, indexRandom)
	// end condition
	if len(vertexBatch) == volumeBatch {...}
}

2.2 Main iteration

Choose 3 points randomly to form a plane and compute the distance between all each point in the batch and this plane. If a point is enough near to the plane (controlled by a given parameter that defines the distance), it will be considered as an inline point. If the proportion of inline points in the batch surpasses a a given line, that plane will be a considered as validated and fed to the next step, otherwise this step will be repeated. After a given number of repetitions, the plane which has the most inline points will be taken into next step if none of them surpassed the minimal limit of inline points.

// pick 3 vertices randomly
	p1 := 0
	p2 := 0
	p3 := 0
	// check if the 3 vertices are different one from another
	for {...}

	// calculate the inlines for the plane formed by these 3 vertices (if the distance to the plane is inferior to the given threshold, this point would be regarded as inline) and take the score
	p := New_p_by_verticesMono(vlist[p1], vlist[p2], vlist[p3])
	for _, vertex := range vlist {
		if mymath.DistPointPlane(vertex.Ply_x, vertex.Ply_y, vertex.Ply_z, p.A, p.B, p.C, p.D) < minDistance {...}
	}
	NumIteration += 1

	// update the score (inline rate). If there is a higher score, note down the plane and 3 vertices which forms it
	score := InlineCount / float64(n)
	if score > highscore {...}

	// quit condition : the high score surpass the given parameter, or the number of iteration reach the maximal given in parameter too
	if NumIteration > minNforLoop && highscore > minScore || NumIteration >= maxNforLoop {...}
}

2.3 General regression

In this step, we compute the distance between each point among the entire sample group and the plane obtained in the last step and get a group of general inline points using the same method and parameter. The general inline points will be removed from the samples for the next iteration. This plane will be put into a pool to be tested its redundancy.

// Fit several points to x plane, using the least square method by solving a linear system Ax=b
. A : matrix of the size n rows and 3 columns. Rows are (xi, yi, 1) in which xi and yi are coordinates x and y of the points
A := mat.NewDense(n, 3, nil)
// B : matrix of the size n rows and 1 columns, (zi), the coordinates z of the points
b := mat.NewDense(n, 1, nil)
// initializing
for i := 0; i < n; i++ {
	A.Set(i, 0, vlist[i].Ply_x)
	A.Set(i, 1, vlist[i].Ply_y)
	A.Set(i, 2, 1)
	b.Set(i, 0, vlist[i].Ply_z)
}

//  compute : A * x = b --> ATA * x = AT * b --> x = (ATA)^-1 * AT * b

AT := A.T()
// multiply A by AT at the left side to prepare for the inverse
var ATA mat.Dense
ATA.Mul(AT, A)

// compute the inverse
var ATAInv mat.Dense
err := ATAInv.Inverse(&ATA)
if err != nil {
	log.Fatalf("A is not invertible: %v", err)
}

var ATAInvAT mat.Dense
ATAInvAT.Mul(&ATAInv, AT)
var X mat.Dense
X.Mul(&ATAInvAT, b)

P := New_p(X.At(0, 0), X.At(1, 0), -1, X.At(2, 0))
return P

2.4 Redundancy test

In this step we will check if there are 2 redundant planes in the pool. For all planes in the pool, we calculate the angle formed by 2 planes. If the angle is smaller than the given value, the 2 planes will be considered parallel. For every 2 parallel planes, we test some point-to-plane distance between them. If the results are small enough, the 2 planes will be considered redundant. For 2 redundant planes, we mix their inline points and use the least square method to generate one common plane.

// check that if the new plane belongs to the same planar facade as one of the planes already obtained. If it is the case, mix the vertices and re-adjust the two planes using least square method
if len(planes) > 0 {
	for index, planeObtained := range planes {
		// first check : if the angle formed by the two planes are inferior to the given value
		// second check : if the average distance from the points forming the old plane and the new plane is inferior to the given value
		if mymath.VectorsAngle32(planeObtained.A, planeObtained.B, planeObtained.C, newP.A, newP.B, newP.C) < maxAnglePlanes && planeObtained.DistAvrPointPlaneMono32(vetexInline) < 4 * maxDistance {...}
	}
}

// if vetexInline != nil (empty slice) means that no similar planes found, we add the new plane to the slice. Otherwise there was a similar plane
if vetexInline != nil {...}

// remove the inline points from the original slice, color these vertices and move them to a new slice to prepare for the painting
for i, index := range indexInline {
	vlist[index], vlist[len(vlist) - 1 - i] = vlist[len(vlist) - 1 - i], vlist[index]
}

vlist = append(vlist[:0], vlist[0:len(vlist) - inlineCount]...)

2.5 Exit condition

We iterate all these steps until the total times of iteration surpasses a given value, or the non-inline points left here are few enough. Then we will return the number of planes generated and their standard equations.

// decide whether the iteration continues
numIter++
if len(vlist) < maxVtoQuit || numIter > iterMax {
	ifContinue = false
}
return planes, len(planes), verticesOfPlanes

2.6 Parameters

The main parameters are here below. They need to be adjusted according to the particualt .ply file that to be analysd.

func PlaneMonoConsecRANSAC32(
	vlist []plyfile.VertexMono,	// the slice which contains all the points fetched from a .ply file
	maxDistance float32, 		// maximal distance between a inline point and its plane
	minScoreforRAN float32,		// the minimal proportion of inline points that a plane need to have in a batch to be considered as validated
	minVforPlane int,			// the minimal number of inline points that a plane need to have among all points to be considered as validated
	maxAnglePlanes float32,		// maximal angle between 2 planes that to be considered as parallel
	maxVtoQuit int,				// exit condition : the maximal number of points which are not considered inline fot any plane
	iterMax int,				// exit condition : the maximal number of iterations
	volumeBatch int 			// the size of batch
	) ([]Plane32, int, [][]plyfile.VertexMono)

Example of usage ：

Planes, numPlanes, _ := calculator.PlaneMonoConsecRANSAC32(vlist, 0.003,
																0.1,
																500,
																0.1,
																50,
																80,
																500)