Algorithm to compute average direction of undirected segments

Question

I have a set of 2D undirected segments, composed of two end-points. Statistically most of them are lying in more or less the same direction.

What I'd like to compute is the average direction of the segment set (for example if the set is globally N/S, it would return something ~ 0°, etc...). Note that I do not care which actual direction is returned (0° or 180° will equally do).

Clamping the direction of each segment in the [0..180°[ range and taking the average does not work (for example two segments, one 1° and the other -1°: the second will clamp to 179° and the average is wrong, here 90°, should be 0°).

I was also thinking of clustering the "normalized segments" end-points in two groups, and computing the direction of the segment composed of the 2 clusters mid-points, but that seems a bit complicated for the task. By "normalized segment" I mean the segment having both end-points on the unit-circle and the middle point at the origin.

Is there an known algorithm/formula for this?

Answer 1

Meta note: this answer calculates the “median” of the given lines. The other answer by MBo calculates the “mean” of the given lines.

---

Let us formalize the problem in the following way. We are given a collection of lines, and we want to find a line p such that the sum of angles between p and all given lines is the minimum possible. Here, an angle between two lines is the minimal of the angles at their intersection point, or 0 if they are parallel or coincide. Thus, an angle between two lines is always from 0 to 90 degrees.

To make things simpler to reason about, translate the lines so that they all pass through the origin. Obviously, this would not affect the answer.

---

To solve this, let us study the derivative of the said sum. Suppose we have an answer line p. Let there be x lines which are 0-90 degrees clockwise from p, and y lines which are 0-90 degrees counter-clockwise from p (x + y = n, the total number of lines given).

Now, rotate p by a small angle α clockwise. The answer will decrease by x * α and increase by y * α. So, if x > y, the answer will decrease, and if x < y, it will increase.

There are two cases where the quantities x and y change.

1. The line p coincides with one of the given lines.

2. The line q is orthogonal to one of the given lines.

Between any two such consecutive points on the circle, the derivative of our sum will be the constant x - y. So, the minimum will be at one of the “angles of interest”: either parallel or orthogonal to some of the given lines. This leads to O(n^2) algorithm: for each of the O(n) angles of interest, just compute the sum in O(n), and choose the angle which gives the minimum sum.

---

This can be accelerated further to O(n log n).

1. Generate the 2 n angles of interest in O(n).

2. Sort them in O(n log n).

3. Compute the answer, and also x and y, for the first such angle in O(n).

4. Move along the circle, maintaining the current answer and the values x and y. In each of the O(n) steps, calculations can be done in O(1).

Answer 2

As I understand the location of segments does not matter, only their direction.

So we can change the problem a bit: we have a set of vectors and we want to fit a line on them.

We can take different criteria for this. A commonly used one is least squares.

For this criteria the solution is:

double dvx=0,dvy=0;
for(const auto &direction:directions)
{
    dvx+=2*direction.dx*direction.dy;
    dvy+=squared(directions.dx)-squared(directions.dy);
}
return std::atan2(dvx,dvy)/2;//or may be +pi/2

Note: for this implementation directions will be weighted by their length, if you want to assign the same weight, direction vectors should be normalized.

This method is sometimes used in determining direction of lines in fingerprints recognition: http://jmit.us.edu.pl/cms/jmitjrn/22/28_Wieclaw_4.pdf

There are several ways to understand this method. One of them is geometrical:

We have a set of vectors with angle alpha[i] from X axis. We don't average these vectors. Instead we build vectors with angle 2*alpha[i], average them and take the half of the resulting angle. The trick is that if opposite directions differ by pi and after doubling they will differ by 2*pi which is no difference at all.

Answer 3

There is method to find mean of angles (in full circular range)

MeanAngle = ArcTan2(Sum{i=1..n}(Sin(Alpha[i])), Sum{i}(Cos(Alpha[i])))

It seems that for you case you can calculate mean of Cosines of direction vectors (because Cos(-alpha) = Cos(alpha)), and get ArcCos (in range 0..Pi)

MeanAngleWithoutDir = ArcCos(1/n * Sum{i=1..n}(Cos(Alpha[i])))

Probably angles should be clamped to (0..Pi) or (-Pi/2..Pi/2) range to avoid ambiguity.