Douglas-Peucker - Shortest arc from a point to a c

Question

I have seen many examples in various programming languages that are using the Douglas-Peucker polyline simplification algorithm to produce a GPolyline to be used on Google Maps. The algorithm, when expressed for polylines on a plan, involves the calculation of the distance between a point and a line (passing through two other points).

Now all the examples I have seen so far are applying the algorithm in a very naïve way, simply by replace x and y by the latitude and longitude. This may produce acceptable results as long as the polyline is very localized, not too close to a pole, and does not cross the 180° meridian, but I would like to implement a more general version of the algorithm.

So, If I am not mistaken, I would need to compute the length of the shortest arc on the surface of a sphere, from a point to the circle passing through two other points of the surface of the sphere, the center of which coinciding with the center of the sphere (the earth).

Does anyone know the formula that computes this length?

Thanks in advance

Answer 1

I'll try to express everything in terms of unit vectors p, q, and r, which can be thought of as points on a unit sphere Σ centered at the origin 0. You can convert that to terrestrial quantities by scaling up by the radius of the earth. There is some background material here.

We want to find the great circle distance d from p to the great circle C going through q and r. C is the intersection of a plane P and the sphere Σ, where P is the plane that passes through q, r, and the origin 0. d is simply the angle θ (expressed in radians) between p and P. The normal vector for P is the normalized cross product q×r/sinφ, where φ is the angle between q and r.

We end up with

θ = arcsin(p⋅(q×r)/sinφ)

As I said, everything here gets scaled up by the radius R of the earth. So the three points are *R***p**, *R***q**, *R***r**, and the distance is Rθ.

However, if all you want is to find a point/line combo with the shortest distance, you can omit multiplying by R. In fact you can omit the arcsin() and just look at the relative sizes of the p⋅(q×r)/sinφ.