Original subtitles

Welcome back to the lecture machine vision. Today we want to discuss chapter number four on curve fitting.

What does that mean?

Let us assume we make pictures like that. For instance in a factory we want to control the output of a machine.

So we make pictures of the objects that are produced.

Based on that we do some edge detection to determine the boundaries of the object.

As a result we get an image as the one shown on the right hand side.

An image that is black for all pixels except for the edge pixels.

Now we can represent this object as a set of edge pixels.

However for further processing doing all computations on the basis of individual pixels might be relatively slow.

And it might not allow to do some higher level reasoning on the shape of the object.

For that reason we might like to represent this object not just as a set of pixels but as a geometric shape.

For instance as a polygon and represent the polygon just by the coordinates of the vertices of the coordinates.

So we could reduce the amount of memory that is required for representing this shape by just six coordinates rather than hundreds or even thousands of pixels.

Now we want to discuss how we can get from this edge image to the geometric representation.

For that we want to start with the simplest elements that we can use namely lines and line segments.

But before we do that let's repeat a little bit how we can represent two dimensional geometric objects.

So this should be just a repetition of what you should already know from your lectures in mathematics.

The first thing that we need is a dot product.

The dot product of two vectors say p and q is defined as a sum over the pairwise products of the elements of the two vectors.

So if p and q are two dimensional vectors and the elements of p are p1 and p2 and the elements of q are q1 and q2,

then the dot product is just p1 times q1 plus p2 times q2.

Within this lecture I will always use the notation with the special brackets that you find here for the dot product.

I could also write p transpose times q that leads to the same result,

but to make clear that we are faced with a dot product I use this special notation with the brackets.

The dot product has some mathematical properties. For instance it is linear in both of its arguments.

That means if we apply a linear combination of vectors p and r with prefactors alpha and beta as first argument,

and another linear combination of vectors q and s with scalar prefactors gamma and delta to a second argument of the dot product,

then we can multiply out the dot product and we get four terms.

Each of those terms contains just a pairwise dot product of two of the basic vectors that are used in the linear combinations

and all the prefactors are multiplied out, are put in front of the dot product.

Another important property which leads to a geometric interpretation of the dot product is that the dot product of p and q is equal to

the Euclidean length of the vector p times the Euclidean length of vector q times the cosine of the angle between vector p and q.

That means there is a relationship between the angle between two vectors and their dot product.

And from that it follows that the dot product of a vector with itself is nothing else than the squared length of the vector

because the cosine of zero degree is one and that the dot product of two vectors which are orthogonal is always zero

because the cosine of 90 degree is zero.

So that means if you have two vectors p and q and p is a unit vector that means it has length one.

And we calculate the dot product of p and q and multiply it with this unit vector p.

What happens? Well, what we get is a kind of projection of the vector q onto the direction that is given by vector p.

So the vector dot product of p and q times p is given as the red vector in the diagram.

So the dot product of p and q provides the length of this vector because the vector p has unit length.

So we can see that the dot product can be used to calculate some orthogonal projection of a vector onto a certain direction.

Based on this we can define lines and line segments.

So the first definition of a line segment with two end points p and q is given as a set of points x

which can be composed as one minus tau times p plus tau times q.

Vector is a real number between zero and one.

That means the definition of a line segment is given by a set of points and all these points are represented in this way.

We can extend that line segment to a full line, that means a line that has an infinite extension,

by relaxing the values of tau to the real numbers, so we do not restrict them anymore to the interval between zero and one,

but to all real numbers. And then we get a representation of the full straight line.

For all values of tau between zero and one we get the line segment between p and q.

For all values of tau less than zero we get those elements of the line which are outside of the line segment between p and q,

but closer to p than q. And for values of tau greater than one we get all those points which are outside of the line segment between p and q,

but which are closer to q than to p.

For some purposes it is useful to be able to calculate the perpendicular point of a certain point onto such a line.

So we assume that a certain point r is given. This point might be part of the line, but it might also be not an element of the line.

And we want to calculate that point on the line that is closest to r. And that point is the point where the angle between the vector that connects r and l

and the vector that connects p and q is equal to 90°. That means this perpendicular point, say l, must be part of the line.

So l must have a representation that is given by one minus tau times p plus tau times q for a certain value of tau.

Furthermore, the vector from r to l should be orthogonal to the vector of p to q. That means the dot product of l minus r and p minus q must be equal to zero.

Now we can substitute in the dot product the right hand side of the first equation. That means instead of l we write one minus tau times p plus tau times q.

Then we multiply out the dot product and we get one equation that depends only on tau.

We resolve this equation with respect to tau and get as result that tau is equal to the dot product of p minus r and p minus q over the dot product of p minus q and p minus q.

So in the denominator you clearly see we are faced with a squared Euclidean length of the vector between p and q.

Beside the representation of a line that we have derived on the slide before, there is a second way to represent a line in the two dimensional space called the normal form.

For that we need a vector n that is orthogonal to the line and has a length one. So it is an orthogonal unit vector.

For that vector we know it is orthogonal, so the dot product of n and the vector q minus p is equal to zero.

We observe that the dot product of n and an arbitrary point x on the line is equal to the dot product of n and one minus tau times p plus tau times q minus one.

So we just substituted x by the definition of such a point that we have derived on the slide before.

Now we can multiply out the dot product and get the dot product of n and p plus tau times the dot product of n and q minus p.

The dot product of n and q minus p is equal to zero, as we have seen, because n is orthogonal to that vector.

So what remains is just the dot product of n and p.

If we reformulate that, then we get that for each point x on the line it holds that the dot product of n and x minus the dot product of n and p is equal to zero.

If we substitute minus the dot product of n and p by a new constant, say c, then we get the standard representation, normal form representation of a line, zero equals the dot product of n and x plus c.

At the bottom right we see a sketch of the situation.

We have the line, we have the two points p and q on the line that were used on the slide before to represent the line,

and we have the unit normal vector n that points into a direction that is orthogonal to the line.

The representation in normal form has some advantages over the representation that we have seen before,

because it allows us to calculate the distance of an arbitrary point to the line in a very easy way.

Let us assume we want to calculate the distance of a point capital R to that line.

With the standard definition of a line, we first would have to calculate the perpendicular point L,

and afterwards calculate the vector that connects L and R and evaluate the Euclidean length of this vector.

With the units with a normal form, we can simplify that and get as result the dot product of n and R.

R is a vector that starts in the origin of the coordinate system and points to the point capital R.

The dot product of n and R plus c.

This is a sign number. The sign tells us whether we are on the left hand side or on the right hand side of the line.

If you are just interested in the distance, we just take the absolute value of this number.

How can we derive that from the plot on the right hand side?

On the right hand side we see the point capital R and the point capital L, which is the perpendicular point to R.

We also see the small vector R that connects the origin of the coordinate system with the point capital R.

We see the vector n, the normal vector.

We see the line that starts in the origin and follows the direction of n.

We know that calculating the dot product of a vector and a unit vector is the same as projecting that vector onto the direction that is given by the unit vector.

So calculating the dot product of n and the vector R is nothing else than calculating the length between the origin of the coordinate system and the point R'.

Furthermore, we know that we are interested in the distance between the points R and L, which is exactly the same as the distance between the points R' and L'.

So what we need is to know how large is the distance between the origin and the point L'.

However, we know that L' is an element of the line. So we know that it has exactly the same orthogonal distance as the point P or the point Q.

That means if we project the point P or Q or L or L' or any other point of the line onto the direction that is given by n, we will always get the same number, the same length.

So this is the same as the dot product of n and P or the dot product of n and L or the dot product of n and Q.

And we know that this is already -c. And therefore we can derive the formula that is shown here.

Furthermore, you might sometimes be interested in getting an arbitrary point on the line if you know the line just represented in normal form.

How do we get that? Well, what could we do? We could pick an arbitrary point - doesn't matter where - in space and then project that point onto the line.

The simplest way is to take the origin as such a point and project the origin onto the line.

That means we are searching a point that has a representation tau times n, where tau might be a real number.

And we are searching for such a value of tau that the dot product of tau times n and n is equal to -c.

And when we resolve this equation we get that the value of tau needs to be chosen as -c.

This yields then one point, namely the point -c times n. And we know that this is element of the line.

And if we look more carefully on the drawing we will find that this point -c times n is exactly the point L' in the drawing.

To finalize this brief repetition of the basics of two-dimensional geometry, let us have a look at vectors.

If we are faced with an arbitrary two-dimensional vector V with elements V1 and V2, we want to calculate a unit normal vector for that vector.

How can we do that? The first step is to normalize the vector. That means to make it a vector of length 1.

This can be done by dividing the vector by its vector length. That means calculating 1 over the Euclidean length of V times the vector.

The second step is to rate-rotate it by 90° either to the left or to the right. In the two-dimensional case this is rather easy.

If we have a vector V1, V2 and we want to rotate it by 90°, we get the vector -V2, V1.

If we want to rotate it by -90°, we get the vector V2, -V1.

That means a unit normal vector for vector V has the shape 1 over the length of vector V times -V2, V1.

Every vector of course has a polar representation in the two-dimensional space.

Each vector V can be represented as R times cos(phi) sin(phi).

That means the vector cos(phi) sin(phi) is a unit vector that points into that direction of that vector V and has unit length.

The scalar R is a non-negative number and it tells us how far away from the origin that point is.

If we are given a vector V with elements entries V1 and V2, we might be interested in calculating its polar representation.

Determining its distance from the origin is rather simple, so we just need to calculate the Euclidean length of vector V.

But how do we get the angle phi?

We have to consider trigonometric calculations here.

The first way that you would find is to somehow use the Arcus Tangent to calculate that vector.

That's possible, but it requires some case distinction to get the angle phi directly.

The simpler way is to use a specific function that is known as I_A_10_2, the Arcus Tangent with two arguments,

that is defined to take the two components of the two-dimensional vector V1 and V2 and that provides this angle phi into which this vector V1 and V2 points.

The function I_A_10_2 is a standard function, so you can use it in mathematical formula as well as in programming languages.

That means every time you need to know the angle of a two-dimensional vector, use this function I_A_10_2,

rather than somehow using the Arcus Tangent of some ratio of some numbers and then adding some case distinction.

Now we have all the prerequisites for our further processing of the edge bitmaps.

We assume that we started with a gray value image and extracted the gray value edges.

Now we know all the pixels that are gray level edge pixels, but we want to get to a shape description.

The first way that we want to consider is to use a technique called the Huff Transform,

which transforms the content of the edge image into a different space and extracts the dominant lines in that space.

From the results of previous slides we know that we can represent each line in the two-dimensional space in its normal form.

If we represent the normal vector as the vector with elements cosine phi and sine phi and the position of a point as the vector x, y,

then we can represent each line as x times cosine phi plus y times sine phi plus z equals zero.

That's just the normal form of a line.

Furthermore, we can limit ourselves to angles phi between 0 and 180° and the values of z from the real numbers.

Why don't we need to consider larger angles between 180° and 360°?

Well, the answer is easy, because each line has two representations, one with a normal vector pointing to the left and one with a normal vector pointing to the right.

Here we restrict ourselves to just one of these representations by limiting the angles between 0 and 180°.

But still we can represent each line in such a representation.

That means we can represent all the lines in the two-dimensional space by two numbers, the pair of phi and z.

Now let us plot this parameter space.

On the left-hand side we see an image space, the original image space in which we detect edge pixels.

It has two coordinate axes x and y as usual in an image.

On the right-hand side we find the parameter space. It has the two axes phi and z.

And each line in the image is represented by one point in that parameter space.

If we consider one point in the image, indicated by the red circle, then we can find several lines that are passing through this point.

For instance this vertical line. It has a certain representation in parameter space.

For instance this line has an angle of phi equal to 0 and a certain offset z.

And that pair of parameters is represented in the parameter space at a certain position that is indicated by the blue circle.

But that's not the only line that passes through the red point.

We might also consider that line that has a different angle and a different offset.

So it has a different representation in the parameter space, which maybe is located at that position.

Like that we go on and consider further lines in the image space, which refer to further points in the parameter space.

For instance the horizontal line in the image space has a representation in parameter space at that point.

And that diagonal line has a representation in parameter space at another position.

And even those four lines are not all. There are infinitely many lines that are passing through the red point of interest.

And all of these have a representation in parameter space.

If we calculate those points we find that those are located on a sine curve.

A sine curve with a certain amplitude and a certain phase shift.

Now let's consider another point in the image. Another edge point that we have found.

For instance the point on the left hand side was indicated by a red circle.

Again we can consider the different lines that are passing through that point.

And we get another sine curve in the parameter space as a splice of all the lines that represent that point.

All the lines that intersect at that point.

And we observe that the two sine curves that we obtained in the parameter space are intersecting at a single parameter set.

And that's obviously the parameter set that represents the line that connects the two red points in the image space.

So what did we find?

If we observe lines in the image space, those lines are represented in the parameter space as points.

If we consider points in the image space we can represent them in the parameter space as sine curves.

And the sine curves implicitly represent all the lines that intersect in the point of interest in the image space.

And if we consider several points on the same line in the image space,

then what we can find in the parameter space are several sine curves that are all intersecting in the same point in the parameter space.

Namely the point that refers to the line in the image space on which all the original points are located.

And this gives the idea of the half transform.

For each point, for each edge point in the image we create the sine curve in the parameter space.

And then we consider the points of intersections of the sine curves.

And those points where many sine curves intersect are likely to represent lines that are very dominant in the image.

So here's the basic procedure.

We first take all the edge pixels and for each edge pixel we calculate the sine curve and we might draw it or somehow represent it, that sine curve.

So for instance here we have this parameter space, this phi and c, and the blue line, the blue curve represents one edge pixel.

The green represents another edge pixel, the red one represents another edge pixel and so on and so forth.

After we have done that we check for points of intersection in the parameter space.

So in this case we might have found one very dominant point of intersection.

That means many sine curves are intersecting in that point for a certain parameter value of phi and c.

And then this parameter value reveals the line on which all those points are located in the image.

In practice, however, things turn out not to be that easy.

Because we often do not get just a unique point of intersection, but we get an area where many points of intersections can be found.

But due to numerical imprecision in determining the point position of edge pixels, etc., they do not intersect in a single point.

But as we can see here, all these points of intersection are very close by, but they are not exactly the same.

Furthermore, it might happen that the image contains several dominant lines, so that we get several of these dominant areas of points of intersection.

And then we would need to extract those.

Furthermore, we get points of intersection which are not showing dominant lines.

So where only, for instance, two sine curves are intersecting, which are not supported by more than two points,

and therefore cannot be considered to be really dominant, apparent lines in the image.

So how can we find those areas with a lot of points of intersection close by?

The idea is to use something that is called an accumulator array.

So for that purpose, we split up the whole parameter space into small cells.

So for instance, we might consider the width of each cell to be, for instance, 1 degree or 10 degrees or 0.1 degree.

And in the vertical direction, we might consider the spacing to be 1 pixel, for instance, height, or 10 pixels height or 0.1 pixels height.

It depends on which precision we want to get.

So now we have a lot of cells, and we partition the space, the parameter space, in many cells.

For each cell, we implement a counter.

So what we do then is that we count how many of the sine curves are passing through each of these cells.

And then afterwards, we can check the counter values and search for the local maxima of the counter values.

And those local maxima reveal the dominant areas of intersection of sine curves.

So in this case, we would have one global maximum where five of the sine curves intersect.

And we might find a second local maximum, which is still dominant.

Where three sine curves intersect.

And then apart from that, we might also find some other areas where only two sine curves intersect.

So now we pick the most dominant local maxima.

That means those maxima, which are supported by very many sine curves.

So in practice, very many means 50 or 100 or 200 or 500 maybe.

And treat those as the dominant lines in the image.

So this yields the Haftrand form.

So in total, summing up, first we create an accumulator array of adequate precision.

Adequate means, well, the more cells we have, the more computational effort we have, the more accurate will be the result.

The smaller the number of cells is, the less computational effort we get.

And however, the precision of course suffers from smaller resolution.

And then we initialize the counters for each of the cells with zero.

In step number two, we go through all the edge pixels.

For each edge pixel, we calculate the sine curve that refers to it in parameter space.

And then we check which accumulator cells are, through which accumulator cells the sine curve passes through.

And increment the counter for those accumulator cells.

And at the end, we search the accumulator array for the dominant local maxima.

And these provide the dominant lines in the image.

So this procedure is known as Haftransform in computer vision.

And sometimes in a more generic measurement system sense, it is also known as a Radon transform.

But for us, the word Haftransform is sufficient to know.

So with the Haftransform, we can determine lines or straight lines.

The Haftransform does not provide endpoints of those lines or line segments for us.

So that means, if we are interested in line segments rather than lines, we need to do some further processing,

which is part of one of the assignments.

So what we need to do for that purpose is somehow determine all the points that refer to the line, that support the line.

Then check where do we find those points in a dense way, where do we have gaps.

And based on that, do some reasoning to determine the endpoints of the line segment.

Let's have a look at some examples.

Again, this image with a cup on a mousepad.

When we use the edge extraction with a KANI operator, we get an edge image that looks like that.

Then we transform everything in the Hafts parameter space.

And visualizing this accumulator array as an image yields a result that looks like that.

On the horizontal axis, the horizontal axis provides the angles between 0 and a little bit more than 180 degrees.

And the vertical axis are the C values.

The brightness of the pixels, so each pixel in this visualization, refers to one accumulator cell.

And the brightness refers to the counter number of the accumulator cells.

That means the brighter the pixel is, the larger the counter value is.

And we clearly see there are four dominant local maxima.

And those four dominant local maxima refer to four lines in the image, namely those that are colored here.

That means the blue, the cyan, the magenta and the green line segment here.

Each of those four refers to one of the four dominant maxima in the image.

Beside that, of course, we also see that there are some other local maxima, but they are not that dominant as those four.

And therefore, if you want to represent this scene with lines, then these four lines are the most important ones.

Let's talk about the properties of the Huff Transform.

As I've already said, the result depends on the size and the precision of the accumulator array.

The more cells we have, the smaller each cell is, the more accurate the result will be.

However, the more computational effort we need.

Determining significant peaks in the accumulator array might be difficult in practice.

I found that there are many bad implementations of the Huff Transform, which you find in computer vision libraries.

And only a few really good implementations of the Huff Transform.

Therefore, if you are working with the bad ones, you might be frustrated very easily.

If you are using better ones, the results are nicer and the quality of the results are better.

So check which kind of implementation you are using.

On the assignments, we will provide an own implementation of the Huff Transform,

which somehow circumvents some of these problems and yields, I think, nice results.

Maybe it's not the quickest implementation, but at least you do not get frustrated that easily.

For instance, if you are trying in comparison the MATLAB provided implementation,

you will easily start to become frustrated from the Huff Transform.

But that's not a matter of the method itself, but more of the implementation.

One property of the Huff Transform is that it ignores all gradient information.

So it adds pixels individually and it does not consider any knowledge about direction of an edge that passes through that point.

However, especially if you are dealing with edge pixels, usually we know the orientation of the edge,

because from the gradient at that position we can derive at least get a rough approximation of the edge orientation.

But this information is completely ignored by the Huff Transform, at least in its standard implementation.

And in natural scenes, if you would go outside and take a picture from some natural scenes,

then it seems that the accumulator arrays easily is flooded and the apparent lines or the dominant lines are somehow not that easily found.

So the Huff Transform works well on technical images where the relevant objects are clearly visible

and the background is not that dominant. If you are going to natural scenes, Huff Transform doesn't work that well in practice.

There are some extensions of the Huff Transform. So the standard form is implemented for lines, detecting lines,

but there are also extensions to work to extract for instance circles or ellipses.

They are not discussed here in the lecture, but the extension is rather easy.

So it is just a small and easy transformation of knowledge that is needed to derive those methods.

For circles, the Huff Transform is used rather often. For ellipses, things become more difficult.

The parameter space gets five-dimensional.

For ellipses, the Huff Transform is not used that often.

Then there is another extension called the randomized Huff Transform.

It avoids some of the problems of the standard Huff Transform, but I don't introduce it in the lecture due to time reasons.

Furthermore, there is one version called generalized Huff Transform that allows to extract any kind of shape.

You provide a certain template of a shape and then the generalized Huff Transform searches the image for those kind of shapes.

Besides the way to use the Huff Transform to extract lines from an image, there are also alternative ways to get lines and other structures from an image.

The second approach that we want to discuss is based out of several steps.

The first step is to somehow transform the representation of edge pixels from an edge bitmap to a list of edge pixels.

That means we somehow collect all the pixels in such a way that we have those organized in lists of adjacent edge pixels.

Based on that representation, we do some polyline segmentation. That means we somehow segment these lists of edge pixels into sublists,

where each sublist represents one polyline segment.

The last step is to fit the lines as best as possible to the edge pixels.

Let's follow this path. First, discuss the edge following, then the polyline segmentation, and at the end the line fitting step.

Edge following. The basic idea of edge following is that we travel through the image,

we collect all the edge pixels which we have found in the image, and then arrange those in lists of pixel coordinates,

where we arrange those lists in such a way that neighboring elements in the list refer to neighboring pixels in the image.

As a result, we get a set of lists of edge pixels. Each list of edge pixels represents one contour,

and the whole set of all these lists represents all the contours which we have found.

On the right hand side, for instance, one contour that will be represented afterwards as a list of edge pixels is the large contour around the mouse pad.

This contour is shown in the image at the bottom with a dark blue color.

All these pixels finally are arranged in one list, and the list is ordered in the sense that neighboring elements in the list refer to neighboring elements in the image.

Beside that large contour, there are other contours that we can find in the edge image, and each of these contours becomes one list of edge pixels.

For instance, the upper bound of the cup is composed out of two different contours, and each of these two contours becomes one edge list.

In the image below, those are represented with two different shades of magenta color.

Besides that, some more contours are found, and those will become individual additional lists of pixels.

The procedure is so simple that I even don't write down an algorithm here.

The idea is we travel through the image from top to bottom, from left to right.

Once we find an edge pixel, we start a new list of edge pixels.

Then we use the gradient information to follow the edge, and we collect all the edge pixels, which we find on that list, until there are no more edge pixels that we can add.

As a result, we get these lists of edge pixels.

If we apply this method to a little bit more complicated scenes, we get the result as shown here.

In these two images, black pixels means there is no edge at those positions, while all the other colors refer to edge pixels.

Each contour that is found, subsumed, or collected in one list of edge pixels is colored with the same color.

On the right hand side, it is very clear that we get several lists of edge pixels.

For instance, one that is colorized on the image with dark green pixels.

And others, which are colorized with black pixels.

At that point, it is worth mentioning that these lists of edge pixels are not yet fixed to a certain object.

Or it might happen that several edges, which belong to different objects or to different areas, are somehow connected with each other.

One example can be seen here.

The list of pixels starts here, then it goes here, here, here, here, here, downward here, and here again.

As we can see, these edges that are collected here belong to different areas.

For instance, some belong to the boundary between the white field marking and the robot,

while others belong to the boundary between the field marking and the soccer field.

There is no guarantee that this list contains some semantic information.

Since we collected all the edges of the pixels in this same object, we need to split them up.

As we can see in this example, as the lists are not only representing lines,

we at least get some information about which pixels belong together or might belong together.

Or several parts, which are put together.

To represent them as polylines, we need to split those lists of pixels up into several parts,

where we can represent each part, each segment, as a line segment.

That's our task. We want to subdivide the pixel list in such a way that the sub-lists can be represented by individual line segments.

There are several algorithms to do that.

Here we only want to introduce one, the so-called Ramadaglas-Poika algorithm.

How does it proceed? It's rather straightforward, I would say.

The assumption of the Ramadaglas-Poika algorithm is that in one list of pixels it represents one line.

That means when we connect the two endpoints in the list,

then the line that we get should somehow fit to all the points in the list.

And if that's not fulfilled, then we somehow split up this list into parts,

so that those parts more likely can be represented by individual line segments.

So, how does it proceed?

Well, as I said, we first connect the first point of the list with the last point of the list,

indicated here by the sine line.

Then we calculate the orthogonal distance of all pixels in that list to this line that we just estimated.

Some of these distances will be small and support the assumption that those pixels belong to that line that connects the first and the last pixel.

Other distances drawn in orange in this image,

they are too large to support this hypothesis.

That means if we find one point with a distance that is larger than a certain acceptance threshold,

then we argue that we somehow need to split up this list of pixels into parts.

So, for that purpose, we select the pixel that has the largest distance from the line between the endpoints,

and we split the list of pixels at that point into two sub-lists.

And then we iterate this algorithm on each of the sub-lists.

So, for instance, after this partitioning step, we repeat the first step and create the line between the first point and the corner point that we just selected.

This yields the first sub-list.

Again, we check the distance between the points in that sub-list and the estimated line.

And in this case, we find that the distance is small enough, that all the distances are below an acceptance threshold.

In that case, we stop the splitting process.

For the second sub-list that was generated, that starts at the new corner point and ends at the endpoint of the list of edge pixels,

we estimate the line by connecting the corner point with the endpoint.

Then we calculate again the distances, and we find that some of the distances are too large.

Again, we select the point that maximizes the distance to the estimated line, which is shown with a red ring.

And at that point, we split up this list into another two sub-lists.

And then we again start this algorithm on each of the two sub-lists that were generated.

That means we create a line that connects the two corner points that define that sub-list, namely the first and the last point of the sub-list.

Then we check whether the distances of all intermediate points are small enough, below an acceptance threshold, and if so, we stop the splitting process.

If not, then we select the point that has the largest distance from that line and split up this list into two sub-lists.

And like that we go on until we cannot split anything anymore.

So if we apply this algorithm on our example image, we will find a result that looks as the right image.

So the left image is the input, that each color refers to one list of edge pixels, and the right picture is the output of the Ramadaklas-Polka algorithm.

So again, each color refers to one sub-list that is a result of that splitting approach.

As we can see, the round shapes are split up into a lot of small segments.

Each of those segments can be represented with a line segment.

And the long straight edges are preserved.

Another example for the other two example images are shown here.

So now each sub-list of edge pixels that were generated is shown by a red line, and the end points are shown by yellow pixels.

So the yellow pixels are the end points, and the red lines represent the list of edge pixels that were represented by this line segment or by this pixel list.

As we can see, the complicated structures are now split into several line segments.

If we zoom in at the university main building, we see how many small segments are created.

We can see that, for instance, at the windows, we get individual line segments for the structures that we find here.

Okay, so now we did the edge following, and we did the polyline segmentation.

And now, as a last step, we do some line fitting.

With the polyline segmentation, we obtained some list of edge pixels.

And we obtained some first initial guess of a line that represents those pixels, namely the line that connects the two end points of the list of edge pixels.

However, this estimate might be suboptimal. It is acceptable somehow, but there might be other lines which fit better to these pixels.

As well, the result of the Huff Transform might also be suboptimal, especially if we use a cell partitioning of the Accumulator array that has two large cells.

On the right-hand side, we see an example. So here, the black points are the edge pixels.

The red line shows which of those are neighboring the list of edge pixels.

And the side line shows the line that connects the two end points.

And maybe this is not really optimal. Especially on the right-hand side, it is far away from most of the pixels.

And we might consider that another line, say that one here in yellow, might fit better than the side line did.

But how can we estimate such an optimal line?

Well, we know that we have to represent the line with a normal vector, n, and an offset parameter, c.

So what we need to do is, we need to somehow determine optimal values for n and c.

But what means optimal in this case?

Well, let us consider again a set of edge pixels and a certain line.

Say this is our first guess. It does not matter how we created this first guess.

So just a first guess, which we want to evaluate. How good is it?

What we can do is, that we can calculate the orthogonal distance of each of those points from the list.

Represented by the yellow dashed lines.

And of course, for a good line, those distances would be small.

In the best case, they would all be zero.

But in practice, it happens that you never find a line that makes all those distances zero.

But somehow, you have to find a compromise.

So now we have to rotate this line a little bit and shift it in such a way that we find the best compromise among the distances that we get.

So we want to find parameters n and c that minimize those distances.

Since we cannot minimize all the distances at the same time, we need to formulate a compromise.

And this can be done by summing up the square distances, d i^2, over all the pixels.

And stating that we want to find those normal vector n and offset parameter c, such that this sum becomes minimal.

So we might treat that as a kind of search problem.

We are varying the normal vector n and the parameter c over a large possible number of values.

And for each of those values, we check the distances of each pixel and then we select the solution with the smallest sum of square distances.

If we do that, we have to consider that n needs to be a unit vector.

So a vector of length 1.

And this needs to be stated as an additional constraint in this optimization problem.

So typically we add that as a constraint subject to the dot product of n with itself must be equal to 1.

So the Euclidean length of n must be equal to 1.

So this states a mathematical optimization problem with constraints.

It states what we want to have.

It does not necessarily or directly state how we get that.

But it already states what is our goal of calculation.

This method is known as the total least sum of squares approach.

So sum of squares, so sum of square distances.

And we want to minimize it, so we want to find the least sum of squares.

And the word total comes from the fact that we minimize the orthogonal distances.

So the next question is how can we solve this optimization problem?

That means how can we find the values of n and c without searching blindly values?

So the mathematical theory from optimization theory states that we have to determine or that we have to build the so called Lagrange function.

The Lagrange function is a function that is composed out of the original function that we want to minimize.

That means the sum over the di-squares.

And adds an additional term for each constraint.

First we have to reformulate the constraints in such a way that they have the form something, some function that depends on n and c equals 0.

We can get that by subtracting 1 in the constraint.

And then we add that constraint to the Lagrange function and multiply it with an additional variable, which is called the Lagrange variable.

In this case it is called lambda. So for each constraint we get one of these Lagrange multipliers.

So the Lagrange function in this case then becomes the sum of the di-squared minus lambda times the dot product of n with itself minus 1.

So if we substitute di by the way we can evaluate di, namely as the dot product of n and xi plus c, then we get the result in the line below.

Now the optimization theory states that the solution of our optimization problem can be found by zeroing the Lagrange function with respect to the original variables of the problem.

That means we have to calculate the partial derivative of the Lagrange function with respect to the parameter c.

And with respect to the normal vector n.

And for each of these derivatives we must zero them and then if we are lucky we can resolve the system of equations with respect to the unknown variables n and c.

So let's do that. First let's start with the partial derivative of the Lagrange function with respect to c.

We get just using the basic mathematical rules for calculating derivatives we get that this derivative is 2 times the sum over of the dot product of n and xi summing up over all pixels i plus 2 times the number of pixels capital C n times c.

And this should be equal to zero.

So if we zero it and resolve the resulting equation with respect to the unknown variable c we get that c must be equal to minus 1 over n times the dot product of times the sum of the dot product of n and xi.

Or if we like to rewrite it we can also write it's minus 1 over n times the dot product of n and the sum of all xi.

And if we want to rewrite that again it's nothing else than minus the dot product of n times the average position of the pixels. So 1 over n times the sum of xi which yields just the average position.

In the next step we calculate the partial derivatives of the Lagrange function with respect to the two elements of the vector n.

That means with respect to n1 and n2. Again just using the basic rules for calculating derivatives we get the two equations here.

We both set those terms to zero.

After that we substitute c by our previous result. That means by minus the dot product of n and the average x position or the average pixel position.

So 1 over n times the sum of xi. And afterwards we somehow rearrange all the terms in such a way that we have terms like something times n1 plus something else times n2 equals something that does only depend on lambda.

And by doing that we get the terms at the bottom of the slide. So as we can see those prefactors which are abbreviated with alpha, beta and gamma they only depend on the pixel coordinates.

They only depend on the xi values. Please note that here for the xi values the second index refers to the first or second entry of the xi vector.

So i refers to the index of the pixel that we address and 1 and 2 as the second index refers to the position in the vector.

So the first element in the vector or the second element in the vector.

So these terms alpha, beta and gamma only depend on the pixel coordinates that we know that we can derive or that we can calculate from the list of pixel positions that we have.

So at the end if we simplify everything we get the equation alpha, beta, a row vector alpha, beta times n equals lambda times the first element of n and a row vector of beta and gamma times the vector n is equal to lambda times n2.

So by rewriting that as a matrix vector equation we get that the matrix, the symmetric matrix that contains the entries alpha, beta, beta and gamma times the normal vector n is equal to lambda times n.

So we have to search for at appropriate vectors n an appropriate value lambda.

And if you go back to your math class you should have learned that these kinds of problems are known as eigenvalue eigenvector problems in mathematics.

So n is known as a so called eigenvector and lambda is called as a so called eigenvalue.

So that means what we have to do is we have to solve this eigenvector eigenvalue problem and find an appropriate eigenvalue for lambda and an appropriate eigenvector n for n.

So from the mathematical theory of eigenvalues and eigenvector it turns out that if we found an eigenvector n for this matrix then we can also scale this eigenvector as we like.

So if n is an eigenvector then 2 times n is also an eigenvector and 3 times n and 4 times n and 0.1 times n as well.

Since we know that n in our case must be a unit vector we can choose an eigenvector that has unit length.

That means we search an eigenvector and if it doesn't have unit length then we scale it in such a way that the scaled vector has unit length.

However still there are two solutions. For these kind of matrices as we find here we always find two real eigenvalues.

So two numbers lambda1 and lambda2 and appropriate eigenvectors for these two values which solve the system of equation.

For these kind of matrices it turns out that all the eigenvalues are real valued and that all the eigenvalues are non-negative.

Though they might be zero, they might be positive but they are never negative.

That is a result from mathematics which I don't like to derive here which is maybe also a little bit beyond the scope of a basic mathematics class but it's still true.

So if you are interested in the result then consider the literature for it to check how this result is derived.

But the fact is these kind of matrices always have as many eigenvalues as there are rows and columns and they are all non-negative real numbers.

So which of those eigenvalues should we choose?

Well a little bit of geometric analysis turns out that the smaller eigenvalue, so in this case lambda2, minimizes the distances.

While lambda1, the larger eigenvalue, maximizes the distances somehow.

And we don't want to maximize the distances, we want to minimize them.

So we choose the smaller eigenvalue lambda2 and for that eigenvalue lambda2 we calculate the corresponding eigenvector with unit length which provides the vector n.

So calculating eigenvectors and eigenvalues is not that difficult for a 2x2 matrix.

If you are interested, look in the mathematics literature for it, then you will find that you just have to solve a quadratic problem.

And that's always possible and therefore it is not that difficult to solve that here.

If you don't want to derive it yourself, you find software libraries that are doing this job for you in more or less all software packages on linear algebra.

Of course including Matlab which also provides a method to do that.

So now let's summarize this total d-squares method to estimate the best line that fits a certain set of points.

The first step is that we need to calculate the pixel positions of all the pixels that we want to consider.

So what we have to do is we have to calculate some sums.

The sum over all x positions, the sum over all y positions, the sum over the squares of x positions, the sum over the squares of y positions and the sum of the mixed x, y product.

And these sums are needed to establish these entries alpha, beta and gamma for the matrix that we need to create.

Then in the second step we create this matrix, alpha, beta, beta, gamma and calculate the eigenvalues lambda and the eigenvectors that are corresponding to these eigenvalues.

We select the smaller among the two eigenvalues and pick the appropriate eigenvector for that eigenvalue and normalize it to unit length.

Once we have done that we calculate the offset value c based on n.

So just by taking c equal to minus the dot product of n and the average pixel position.

So and then you are done. Then you know the line that represents a certain set of pixels.

But still you don't know where the endpoints are if you are interested in line segment and not an infinitely extended line.

If you want to do that then you need additional techniques which you will find in the assignments and the exercise classes.

Now we analyze again, we analyze somehow where the pixels are located along the line and then we select which points are the endpoints.

So have a look at the result of this calculation.

So here is a small clipping of this example image of the main building of KIT.

The left upper image shows the result of the line segmentation, of the polyline segmentation step of the Ramad-Douglas-Poiker algorithm.

And for instance at that position here we find that the line that is actually creating these edges is vertical.

But these lines that connect endpoints are not vertical but they are skewed a little bit.

After we have fitted all those lines again with the total least squares approach we find that the situation has improved.

And that now the lines that we have estimated are really aligned with the structure that they should represent.

So it is a small improvement but it is an improvement.

To summarize this part of total least squares let us make a comparison.

You might know already least squares methods from other lectures.

For instance from statistics lectures in which as well least squares methods are often used.

Those are typically so called ordinary least squares method.

And the question is what makes the difference between the total least squares method that we consider here and these ordinary least squares method that you might know from other lectures.

Well, the task itself is already a little bit different in the two methods.

In the total least squares approach we deal with points in the two dimensional plane and we don't care about the orientation of the coordinate system in that plane.

So we don't treat the x-axis in a different way than the y-axis.

The result should be independent of the choice of the coordinate axis.

While in the ordinary least squares we are not fitting a line in a two dimensional plane but we fit a linear function to a set of points.

And here x and y have different meanings.

So x is the independent variable, y is the dependent variable.

And we want to find the relationship in which way y depends on x.

So we treat y as a function of x.

So that's something different.

We are not interested in the line as such.

We are interested in a linear function and the line is just the graph that visualizes the linear function.

As a consequence in the ordinary least squares approach, the graph that we get, the function that we get, of course depends on the choice of the coordinate system.

So that's relevant.

So the method is anisotropic if you went to state it in one word.

While in the total least squares the result is independent of the choice of coordinate systems in the two dimensional plane.

So this result is isotropic.

Yeah, and what do we minimize?

In the total least squares as we have seen we minimize the orthogonal distance of the points to the line.

In the ordinary least squares we do not minimize an orthogonal distance.

What we minimize is the distance between the point and the linear function just considering the difference in y values, not considering differences in x values.

And that, of course, is an important difference so that the result that you would get with ordinary least squares looks different than what you get with total least squares.

And furthermore with ordinary least squares you will not be able to fit something like vertical lines.

With total least squares no problem.

So if you are dealing with least squares methods always check which case you are facing.

Are you facing a case where the coordinate system doesn't matter?

Then choose total least squares.

If you are treating a case where you have a functional relationship between y and x then choose ordinary least squares.

The total least squares approach has many advantages.

It can be calculated very efficiently.

It yields nice results.

However, it suffers from one problem because it is not robust with respect to outliers.

What is an outlier?

Well, an outlier is a point that deviates very much from the line.

That means if we consider the points which you find below, total least squares performs well and there is no problem.

But there is also not a single outlier there.

If we vary one point to that position this point becomes an outlier.

It is somehow an errorless point if you like what I like to say.

And total least square is not good in dealing with these points because the estimate that it would yield looks more like that.

And as you can see the line is queued and does not really fit to the majority of points anymore.

So what we need is a method that is more robust with respect to these outliers.

So that a single outlier does not have such a strong impact on the solution.

There are several ways to implement these ideas.

One idea is to reduce the influence of outliers.

That means those methods still consider outliers and outliers still have an influence, but not that strong as with total least squares.

One of those methods is called M estimators and it is discussed in the lecture.

Another way is to ignore outliers completely.

There are as well several methods.

One of that is called RANSAK which will be discussed in the lecture.

So first of all let's have a look at M estimators.

When we consider total least squares we minimize the sum over the squared distances.

So the squared distance is plotted on the right hand side with a red curve.

This one here.

And obviously the larger the distance is, so that means the farther away from zero on the horizontal axis we are, the stronger increases di^2.

That means for small deviations di^2 still is small.

For large distances di, di^2 becomes huge.

Now we might think why do we use the square function here?

Can't we use functions that do not grow that quickly as the square function?

And this gives the idea for M estimators.

So instead of using the square function we use another function, say rho, which is often called the loss function of the estimator, and apply di to rho.

And for rho we use functions that do not grow as quickly as the square function.

But still in the vicinity of zero we want to preserve the nice properties of the square function,

so we search for those functions which have similar behavior for small values and which do not grow that quickly for large values of di.

On the right hand side we see several possible options and the formula.

So the standard way of total e^2 would be the square function, which is shown here.

So the farther away from zero we get the stronger the square function looks like.

The first way to reduce this influence is to use the so-called Huber function.

The Huber function is drawn in blue.

The Huber function is exactly the same as the squared error term in the vicinity of zero.

But from a certain point on, from here on, it does not grow radically, but it only grows linearly.

So this can be also seen here in the formula.

For in the vicinity of zero, in the first of these two cases, the function grows with a square of d.

And in the second case, if we are away from zero, then the growth depends linearly on d.

The parameter k, that defines the switching point between quadratic and linear behavior,

has to be chosen manually according to the problem.

So k has to be chosen in such a way that typical inliers are smaller than k, have distances smaller than k,

and typical outliers have distances larger than k.

And obviously the growth of the Huber function is smaller than the growth of the square function.

That means points with large distances, outliers, have smaller impact on the solution than for the squared loss function.

We can even further decrease the influence of large outliers using other error terms.

For instance, one choice would be the Cauchy loss function that is plotted in green.

In the vicinity of zero, it approximates the squared error term.

But for larger values of d, it does not grow that much anymore.

So the growth itself, so the slope of the function decreases for larger values of d.

That means here, cross-outliers have even less influence on the result than for the Huber loss function.

And another choice, which is even more strict than Cauchy, is the Tucky loss function plotted in violet,

where points which are far away do not have a stronger impact than points, for instance, at this position.

So here the influence of cross-outliers is limited, even completely limited.

All these functions have some parameters, so for instance the kappa in the case of the Cauchy function,

the k in case of the Huber function, the a in case of the Tucky function,

and those values have to be chosen appropriately depending on the problem.

That was the M estimator. The second option is to ignore outliers completely.

And how can we implement that? Well, we could also say it's a kind of M estimator,

but this M estimator ignores points completely, which are far away than a certain threshold.

That means we could say we want to minimize the sum over sigma of di,

where sigma is a function that is equal to zero if the distance is smaller than a threshold,

or it's one if the distance is larger than a threshold.

That means what we do with this kind of loss function is we just count the outliers.

We just count those points which are far away from the line than a certain threshold theta is.

Theta again is a parameter that we have to choose application dependent in an appropriate manner.

So here the distance itself doesn't matter at all.

It just matters whether it's larger or smaller than the threshold.

So for this case here, that we can see here at the bottom,

we can argue that for this choice of line, if we assume a certain threshold size theta,

which is shown here by the dashed lines, so the dashed lines have a distance from the solid line of exactly theta,

we just count the number of inliers and the number of outliers.

So in this case the number of outliers is two. So the error term will yield an error of two.

For that choice that we can see here, well if we count how many points are outside of this band between the two dashed lines,

we find it's eight points. So the error term is eight. That's a bad fit.

This idea is known in literature as the RANSAC approach.

RANSAC stands for random sample consensus. So why random sample consensus?

Because the algorithm that implements this idea uses some randomization.

It randomly picks two points from the set of edge pixels.

Then it creates a line that passes through these two points.

So afterwards it checks how many outliers do we have for that line.

And this process is repeated again and again. And since we always pick randomly two points,

the solutions that we get in each trial will be different.

And at the end we check which of the trials yielded the best solution, the smallest number of outliers.

And that's the solution of the RANSAC algorithm.

Of course, if we like, we can add one additional step.

Once we found the best solution, we know which points are outliers and which points are inliers.

Based on that, we could do a total least sum of squares fit only on the inliers to obtain an even better solution.

Here is a procedure shown as a kind of processing plot.

So we start in the first trial, pick randomly two points, estimate the line parameters, evaluate the error and obtain the number of outliers.

Then we repeat that in a second trial, independent of the first trial and so on and so on.

Up to the case trial.

Usually based on an assumption of how many outliers we typically expect to have,

we determine the number of trials in such a way that the probability of randomly picking two outliers all the time in all the trials is so small that we can neglect it.

After we have done that, we pick the trial with the smallest number of outliers and that yields the line parameters.

Here is an example again, a sequence of points.

In the first trial, we might pick the two red points and fit a line to them and then check how many outliers do we get.

In this case, for that choice of theta, we get six outliers.

So six points outlier of this band that is bounded by the dashed lines.

Independent of this first choice, we repeat the process with the second trial and in the second trial we pick randomly another point.

So these might be the points which are now shown in red.

Again, we connect these two points, fit a line to them and count the number of outliers.

In that case, we find seven points outside of the band which is bounded by the two dashed lines.

In a third trial, we might pick these two points shown in red, estimate a line and count the outliers equal to three.

And in a fourth trial, we might pick those two points which are unfortunately one of them is an outlier.

We again estimate the line that connects both, count the number of outliers and get ten.

So at the end, we select the trial with the smallest number of outliers and we find, yeah, it was the trial number three.

That means that becomes the overall solution of the RANSAC estimate.

Let us compare the properties of the two estimators for robust estimation, the M estimator and the RANSAC.

The idea in the M estimator is to use a loss function that puts less influence on points that are far away on the outliers.

So the outliers still have an influence on the solution, but the solution is smaller than for the sum of squares approach.

In RANSAC, we completely ignore the outliers.

We just consider inliers and ignore the outliers completely.

So it doesn't matter how far a point an outlier is away from the line.

The influence does not depend on that. The influence is just zero.

The parameters which we need to determine to use that method is, well, for the M estimator, we have to determine the error terms or the loss function

and maybe the width parameter of the loss function.

For the RANSAC, we need to choose an acceptance threshold and we need to determine a number of trials that we run RANSAC.

And the algorithm, well, which algorithm do we use for the M estimator that becomes a numerical optimizer that minimizes a kind of weighted least squares.

So the implementation of the M estimator is traced back to the total least squares method, but with some additional weights for all the points.

And that yields a numeric iterative solution algorithm that calculates the solution in an iterative manner.

For RANSAC, we make repeated guesses from pairs of points and then check the number of outliers.

So, I've implemented or we have implemented together with Ahivi a demo tool that allows to play around with these methods and to check how the solutions change.

Unfortunately, in this slideshow, it is hard to integrate this tool.

Therefore, I will create an extra video with just this demo tool and put it on Ilya. So please check this extra video with a demo tool.

Now we have seen how we can fit lines and how we can determine polylines from edge images.

However, lines are not the only structures that are relevant in practice.

Often we want to estimate circles or pieces of circles and ellipses to structures.

So the question is how can we fit those structures to edge lists?

Well, have a look at those on the subsequent slides.

First of all, we need to know how we can determine or how can we represent circles and ellipses.

So for lines we already know that there is a normal form which is very suitable for the estimation task.

For circles, the representation must be different.

In a classical way, a circle is represented as a set of points that all have the same distance from a center.

So what we need to know is the center point and the radius of the circle.

So this classical approximation would look like that.

So if we want to consider all points x, y, which are on a certain circle, which is defined by its midpoint, m1, m2, and its radius r,

then we would say the set of all points for which x minus m1 squared plus y minus m2 squared minus r squared equals 0 describes all those points.

Or in other words, the set of all points x, y, which have a Euclidean distance equal to r, describe that circle.

If we want to determine the distance of a point from the circle, what do we do?

Say we want to consider this point vector x and calculate the smallest distance of that point from the circle.

What we do is we calculate the distance of that point from the center point of the circle and subtract the radius.

That means the Euclidean distance of a point x, y from the circle becomes the square root of x minus m1 squared plus y minus m2 squared minus r.

This is a signed number, of course, so the sign tells us whether we are inside or outside of the circle.

If you are not interested in the sign itself, we take the absolute value of this number.

So the first term with the square root is just the Euclidean distance between the point x and the center point, and this distance is compared with the radius of the circle.

This is known as the Euclidean distance of a point from the circle.

However, in practice it is useful to consider also another kind of distance measure,

which is not a Euclidean distance for sure, but a related distance, which is also known as the algebraic distance.

In that distance we do not calculate the Euclidean distance between x and m, between the point x and the center point, but the squared Euclidean distance,

and subtract the squared radius.

Of course, for all points on the circle this distance also is equal to zero.

And for all the points outside of the circle, which are not exactly on the circle, this distance grows with increasing distance from the circle.

That makes sense.

But indeed it is not the same as the Euclidean distance.

But the advantage is that we get rid of the square root, and square roots are hard to deal with in efficient calculations.

If we compare these two distance measurements, we find that the distances that are calculated for certain points look like that.

Here on the plot we see on the horizontal axis the distance or the position of a real point.

Zero would mean the point is located in the center, at the center point of the circle.

R means it is located on the circle, so at a distance of R from the center.

And all the points that are larger than R means the point is outside of the circle.

And the vertical axis, or the two vertical axes, better to say, represent the distance measure that is given by either the Euclidean distance in blue,

or the ultra-biotic distance in red.

For the Euclidean distance we obtain a result that is intuitive in what we expect.

So the error increases in a linear way as far away from the circle we get the larger the Euclidean distance is.

Compared to that, the algebraic distance increases faster for points outside of the circle and increases not that fast for points inside of the circle.

However, for points which are close to the circle, both distances are very similar.

So both distance measures do not deviate considerably as long as we are close to the circle itself.

But for points which are far away from the circle, both distances differ considerably.

Now let's use again some method to estimate the circle parameters from a set of points, from a set of noisy points which are scattered around a circle.

So here we find the points in black, and as we can see they deviate rather strong from the real circle.

So the impact of the algebraic distance here becomes apparent.

And if we start to minimize the squared distance, the squared sum of distance of all the points,

in order to obtain the parameters of the circle, we get two different solutions.

If we use the Euclidean distance, we get a radius that is smaller than what we get when we use the algebraic distance to solve the minimization task.

Of course, this is what we expect, because in the algebraic distance points that are outside of the circle have more impact on the solution,

because the algebraic distance grows more rapidly for points outside of the circle than for points inside of the circle.

And then we get a deviation like that, as we can see it here.

However, if the points are not scattered around in such a strong manner as we can see it here in the example,

but if the points are just scattered a little bit around the real circle, then the influence of the algebraic distance remains rather small,

and the results are almost the same.

That means if the points are spread a lot, then it is better to use the Euclidean distance, although the computation is more difficult.

If the points are not spread very much around the circle, then the algebraic distance is very suitable.

And the advantage of the algebraic distance is that the computation of the circle is much easier.

So, if you want to minimize the square Euclidean distance for a circle, then we cannot solve that analytically.

We do not find a closed-form solution, but we need a numerical optimization.

That means an iterative approach that needs several iterations to find the solution.

For the algebraic distance, that is not the case, so we find a closed-form solution.

For that purpose, we rewrite the algebraic distance a little bit.

Instead of writing x minus m1 squared plus y minus m2 squared minus r3 squared,

we can multiply these brackets and then rearrange all the terms.

Then we get equal representation, first terms that do not depend on the parameters of the circle itself, that is x squared plus y squared,

then parameters or some terms which do not depend on x and y, this is m1 squared plus m2 squared minus r squared,

then some terms that depend on x, which is minus 2m1 times x, and some terms that depend on y, which is minus 2m2 times y.

Now, if we rename minus 2 times m1 as a, minus 2 times m2 as b, and m1 squared plus m2 squared minus r squared as c,

then we get an algebraic representation that has a form ax plus by plus c plus x squared plus y squared.

We find that there is a unique mapping between the old parameters of the circle m1, m2, r,

and the new parameters with which we represent the circle a, b, and c.

So a is just minus 2 times m1, b is equal to minus 2 times m2, and c is equal to m1 squared plus m2 squared minus r squared.

So there is a unique mapping between these two parameter sets, so that we can represent a circle either using m1, m2, and r, or a, b, and c.

So now what we need to do is we want to minimize the sum of the squared algebraic distance of all the points that we consider for the estimation of the circle.

So we sum up over all the points, so i is the index for each point, and then we take the squared euclidean distance of that point.

Now we can again calculate the partial derivatives of this term with respect to a, b, and c,

zero those partial derivatives, rewrite everything and arrange it as a matrix vector multiplication.

And after some steps of calculation, which are not really difficult, but which are a little bit too lengthy to show them here on the slide,

we get as a result a matrix on the left hand side, a 3 by 3 matrix, times the vector of unknown parameters a, b, and c, equals some vector on the right hand side.

And we see that the matrix on the left hand side and the vector on the right hand side are built out of the coordinates of the pixels that we want to consider.

So that means once we know the pixels that we want to use to estimate the circle parameters, we can first create this matrix and this vector on the right hand side.

Then we get a system of linear equations and as long as this matrix has full rank, which is typically fulfilled if we consider at least three points,

then we can resolve the system of equations and get solutions for the unknown parameters a, b, and c.

Once we have found a, b, and c in such a way, we can calculate m1, m2, and r in a straightforward manner, as it is shown here.

So that means we find a closed form solution in terms of a linear system of equations to calculate the unknown parameters of the circle.

And that is computationally rather efficient.

Let's have a look at how this kind of circle estimation works in practice.

And for that purpose, let's consider some work that has been done in our group many years ago by Christoph Speck, a former PhD student.

He actually was working in the field of forensic computer vision and one of his tasks was to analyze bullet casings.

So bullet casings are say parts of bullets of guns and once you shoot with a gun, the bullet cases get some typical structures impressed by the gun.

So that means if someone was murdered and police comes, they are typically interested in analyzing these bullet casings

because the structures that you find on the bullet cases give some hint or give some cues on the gun with which the bullet casing was shot.

That means once you find a gun, you can use that gun in an experimental shooting and then collect the bullet casings there as well

and then compare whether the structures that you find are the same.

And then you know whether the bullet was shot from that gun that you found or whether that gun was not the one that was used for shooting.

So in this case, to give you an idea, mainly those structures here in the center matter.

The structures which you see outside here are not that relevant, but the structures inside.

Because they are created by the gun when the bullet is shot.

And as you can see, we find a lot of circular structures here.

And if you have a considerable microscope or a certain camera with suitable magnification, we get pictures like that.

And one of the tasks to analyze these structures is to estimate the circles.

On the right hand side, we get the circles that are estimated using first some kind of half transform for circles to get some basic idea where circles are located

and afterwards use some circle fitting with the algebraic distance to get the final fit.

And we clearly see that those circles are nicely fit by this method.

Now let's extend our findings to ellipses.

Ellipses can be seen as generalized circles.

They also have a center like circles.

They have radii, but ellipses have two radii, two different radii, not just one like circles.

And they have a turning angle that describes the orientation of the ellipse.

The figure on the right hand side illustrates these properties of ellipses.

So here we have the center point of the ellipse.

Then we have two main axes.

The axis with a larger radius R1 and the orthogonal axis with a smaller radius R2.

Furthermore, these two axes do not necessarily have the same orientation as the coordinate axis of that plane,

but they might be turned with respect to that coordinate axis.

And this turning between the coordinate axis and the main axis of the ellipse is given by this turning angle theta.

That means in total we need five parameters R1, R2, the two parameters of the vector M,

and the angle theta to describe an ellipse.

However, in practice it is sometimes convenient to use other parameterizations of ellipses.

So as we have seen it for circles, where the usage of the radius and the center is just one way to represent an ellipse,

there are also other ways to represent an ellipse.

And the general way to do it is to use six parameters A, H, B, G, F and C.

And to state that all points that meet this equation A times x squared plus H times xy plus B times y squared plus G times x plus F times y plus C equals 0

are those points which are located on the ellipse.

We need one additional constraint, namely that 4 times AB minus H squared is greater than 0.

Only if this condition is met we are really facing an ellipse.

If this constraint is not met, then we might be faced with for instance a hyperbola, which we are not interested in here.

So let's again stick to this representation with six parameters.

However, we have seen already that five parameters are actually enough to represent an ellipse.

But in this modified representation we have six parameters.

So there is one degree of freedom in the choice of parameters.

That means if we find a suitable parameter set A, H, B, G, F, C,

we could multiply all these parameters with the same positive constant

and we would still get a representation of the same ellipse where the representation is different.

So one ellipse has several representations of this form.

This is not really nice if we want to find the best parameter set,

because it introduces some problems in the optimization algorithm that should find the best fit.

Therefore we want to somehow restrict these parameters in such a way

that we get rid of this additional degree of freedom.

There are several possibilities how we could do that.

One way would be to state that A must be equal to 1, or that A plus B must be equal to 1,

or that the Euclidean length of the parameter vector is equal to 1.

In some papers you even find the possibility to set C equal to 1.

However, the last version is a little bit tricky,

because that solution is not invariant with respect to translation.

So there are some ellipses which cannot be represented with this choice of C.

There are several algorithms to solve this problem, how we can fit an ellipse to a set of points.

Meanwhile the most popular one is the approach of Fitzgibbon, Pilou and Fischer.

I don't want to go into details, because the mathematical details are a little bit more difficult.

It is based on a generalist eigenvalue problem,

so it becomes a little bit more difficult than what we have seen in the line fit problem,

where we had a normal eigenvalue problem, now we have a generalist eigenvalue problem,

and that is then more difficult to be solved numerically.

I don't want to go into these details, but I just want to mention there is this approach,

and of course you also find implementations of it in most software libraries for computer vision.

So I recommend to use this method if you need it.

It is based on the idea to use the constraint 4ab minus h squared equal to 1,

to get rid of this degree of freedom, and based on that it is possible to get a nice mathematical modeling.

When we use this method on the example picture with a cup on the mousepad,

and apply it to those structures which are really elliptical,

we get the results that are shown on the picture at the bottom.

So the colored ellipses are those that are fitted to the edge lists that we find.

And if we for instance consider the violet, or the green, or also the cyan ellipse that we estimated,

we can see that they really fit well to what we expect.

For the red ellipse that's a little bit maybe unexpected, because it also adapts to non-elliptical parts of the cup.

But still it's not really bad what we get.

So that's all we want to discuss about line fitting, or about curve fitting in general.

So let's summarize our findings.

First we introduce the half transform as a method to extract lines from an image.

The idea is to consider all the edge pixels for each edge pixel,

create all the lines that pass through these edge pixels,

and then use the parameter space and increment all the cells in this accumulator array in the parameter space,

which represent lines that pass through this point.

And afterwards we are searching for local maxima in this accumulator array,

and those local maxima provide to us some dominant lines in the image.

So a very general approach.

It does not use gradient information, it is somehow reasonable, it yields reasonable results.

However, it is not considering any local information of a pixel and its vicinity.

It's just considering there are pixels, but it does not consider neighborhood of pixels,

it does not consider gradient orientation or something like that.

The other way was to collect edge pixels in a list of edge pixels,

using the neighborhood and using the gradient information to collect edge pixels, which are somehow belonging to the same edge.

First doing some edge following and then applying the Ramada-Klas-Polka algorithm

to split up those edge element lists into smaller pieces,

where we can assume that each of the smaller sub-lists can be represented by an individual line.

And then we were discussing line estimation techniques, starting with the total least squares approach,

that just minimizes the sum of squared distances of points from the line.

This yields somehow very good results, but it suffers from outliers.

So if we are faced with outliers in the data, then this approach does not perform as good as we would expect it to be.

To avoid these problems with outliers, we get to know M-estimators,

which are somehow changing the error term that is used, the loss function, to reduce the influence of outliers.

And we discussed the RANSAC, which completely ignores outliers and just counts outliers,

but does not distinguish between an outlier that is far away and that is very far away,

but only considers there is an outlier or there isn't.

Finally, we discussed the case of circles and ellipses.

We have introduced several ways to represent circles and ellipses as parametric curves,

and we have then derived two different distance measures, the Euclidean distance and the algebraic distance.

We have found for circles that it is easier, computationally easier, to deal with the algebraic distance,

but the results are more accurate if we would use the Euclidean distance.

But at least if the points that we consider for circle estimation are not scattered too much,

as long as they are still very close to the circle,

the algebraic distance approximates the Euclidean distance sufficiently well.

While the points are scattered very much, then the algebraic distance overestimates the radius of the circle that we get.

For ellipses, we just introduced or mentioned the algorithm with which we can estimate the ellipses without going into details.