[patch] add curvature based distance calculation to LatLon class, add way segment lengths between nodes to InspectorDialog

[cmuelle8]

An elliptical function for length calculation has not yet been added.

AUTO value of DistanceTypes in LatLon class should later on decide sensible what distance Calculation to use - i.e. in very high zoom with close LatValues, curvature based algorithm can be used.

Additionally the ellipsoid and datum are now gettable from josm's main projection instance..

Additional info on the net can be found here:

​http://clynchg3c.com/Technote/geodesy/radiigeo.pdf

and here:

​https://en.wikipedia.org/wiki/Earth_radius#Radii_of_curvature

Greetings

ps: please review and test - this is core functionality, so please do not apply without having tested it a second/third time - that said, the main distance function is still greatCircleDistance - performance-wise I guess we could use the curvature based function as well, and it promises to be (a little) more accurate when editing locally - comparing the functions at runtime can be done by drawing a way and then using Ctrl+I to get the Inspector..

[cmuelle8]

​http://clynchg3c.com/Technote/geodesy/radiigeo.pdf seems to have a more recent copy

[cmuelle8]

euclidianDistance still needs rework, it's currently not used actively except in Inspector summary..

[stoecker]

The patch of MapStatus.java is broken.

[cmuelle8]

MapStatus.java patch has been broken due to a collision in mid air..

The euclidianDistance() function interprets the output of distance() as an angular distance d. distance() itself calculates length of a line in a Cartesian coordinate system, with a value range from -180 to 180 on the x-axis and -90 to 90 on the y-axis.

d is a fraction of a full circle (360 degrees or 2*pi), so euclidianDistance() scales that fraction to the same great circle eventually used by greatCircleDistance() - given by earth's mean radius.

[cmuelle8]

josm-add-length-calculation-based-on-curvature_enhanceWayInfo.patch

[cmuelle8]

updated javadoc in code, stress the error conditions introduced using distance() as a basis for calculation of euclidianDistance(), the most important ones being:

• cannot cross projection bounds to eventually find a shorter distance on the great circle
• error introduced by stretching longitudinal degree to the same unit distance for every latitude

[bastiK]

If you want to be precise, there is a problem with the curvatureDistance algorithm:

If the line is divided into segments, to get better accuracy, the intermediate points should lie on that line. But LatLon.getCenter() is very approximate and will not give a point on the shortest path.

Have you checked ​https://en.wikipedia.org/wiki/Geographical_distance and ​https://en.wikipedia.org/wiki/Vincenty%27s_formulae ?

[stoecker]

Replying to cmuelle8:

MapStatus.java patch has been broken due to a collision in mid air..

No. You can't replace a configuration based variable check with an static if. The necessary minimum is to provide also the additional display as configuration option.

[cmuelle8]

Replying to stoecker:

Replying to cmuelle8:

MapStatus.java patch has been broken due to a collision in mid air..

No. You can't replace a configuration based variable check with an static if. The necessary minimum is to provide also the additional display as configuration option.

The patch does not replace the configuration based variable! It only adds a decimal place IF distance values are smaller than 100 - that way space is optimally used. If we have large distances, we usually are not interested in many decimal places (and hence depend on default or configuration setting, like it has been).

It's called value range dependent decimal formatting. If it's done right it obsoletes the need to have a setting variable at all.

Actually, what you want to do for a fixed width field is show a fixed amount of significant digits (and not a fixed amount of decimals in general).

For small numbers you would only limit the amount of decimal numbers shown, because the application precision is limited (by ortho or gps precision, or by the method of distance calculation used). But then, even if underlying imagery is not very precise you can get relations to other ground objects right (i.e. you may fit smaller objects you know about to bigger ones identifiable by the imagery - the overall result will have correct proportions and be only wrong in an /overall/ offset to reality).

For instance: if you know that a parking space on ground is exactly 2,90 m wide you would want to use that exact width to map it, even though we may not be able to place the overall object using a 0,01m resolution.

[cmuelle8]

If I wanted to configure that with the current state of the code, I would get two decimal points for large distances as well, which is neither wanted nor needed.

An alternative might be to just add another decimal to whatever is configured, but a better alternative is to just use a full range dependent decimal formatting, i.e. one that always shows at least 3 or 4 digits, but only 2 or 3 decimal points depending on the value range.

[cmuelle8]

Replying to bastiK:

If you want to be precise, there is a problem with the curvatureDistance algorithm:

If the line is divided into segments, to get better accuracy, the intermediate points should lie on that line. But LatLon.getCenter() is very approximate and will not give a point on the shortest path.

Have you checked ​https://en.wikipedia.org/wiki/Geographical_distance and ​https://en.wikipedia.org/wiki/Vincenty%27s_formulae ?

Yes I have, I want to include a function computing the definite integral of the elliptic curve eventually later on. But first things first.

You're correct that LatLon.getCenter() will not /always/ give a point on the shortest path, but please explain why you think that it is very approximate when it does?

I've already written in the ticket description that I do not intend to use curvatureDistance() as a general replacement, but for high zoom situations only (in which case LatLon.getCenter() /is/ guaranteed to deliver a point on the shortest path, at least for ​merc projections (if the prime meridian is not greenwich, we're just talking about an offset fix)).

EDIT1: To use it generally (low zoom), the shortest path problem could IMO be fixed by simply using the antipodal point to m if this.distance(other) > 180 degrees.

EDIT2: It is true, that a.getCenter(b) will not produce a point at exactly half of true distance, maybe this is what you're referring to, but please explain why it should not lie exactly on that line!? For the algorithm to work correctly we solely need it to be on that line, it does not need to hit the middle of the line.

[bastiK]

Okay, but it should be more clearly documented that this method can only be used for distances < 10 km.

In principle it is strange to interpolate two LatLon values for precise distance measurement. This introduces errors that depend on the orientation of the segment and the latitude. Near the poles, it goes crazy.

Why not simply use greatCircleDistance and multiply the mean radius?

There is a nice online tool for cross-reference: ​http://geo.javawa.nl/coordcalc/

[cmuelle8]

Replying to bastiK:

Okay, but it should be more clearly documented that this method can only be used for distances < 10 km.

Should not can. Should because the recursion involved for large distances makes it potentially more expensive than using greatCircleDistance, not because it gets erroneous (once >180 degree case is fixed, of course).

In principle it is strange to interpolate two LatLon values for precise distance measurement. This introduces errors that depend on the orientation of the segment and the latitude. Near the poles, it goes crazy.

Nope, I do not get your "point".
See source:/trunk/src/org/openstreetmap/josm/data/coor/LatLon.java#L313​

greatCircleDistance does exactly the same: other.lat() - this.lat()) / 2 (but for means of direct calculation, not for means of testing the distance result) There are much greater errors introduced at other points in the calculation ( toradians, sin, cos, you name it ).

Why not simply use greatCircleDistance and multiply the mean radius?

Because the method I use is documented elsewhere (see javadoc) and curvature radii functions in class Ellipsoid are there for a reason :)

EDIT:
The problem interpolating LatLons you refer to (esp. at the poles) appears if you measure their distance on the mercator projection directly - this is what euclidianDistance does, and you can see the effect easily by drawing a straight line from (80,-170) to (80,170) and then hitting Ctrl+I.

curvatureDistance does not do this. It simply uses the interpolation to get another point on exactly the same line. All distances measured use a very similar method to the one in greatCircleDistance, with the exception that for the east-west bound integration another "great circle" radius is used than for north-south bound integration. I've explained this thoroughly in the javadoc to the function, just like I've explained the "mistakes" euclidianDistance deliberately makes.

[bastiK]

Take for example the two points

N 81.03862° W 79.10156°
N 80.17871° E 20.74219°

Your code reports

      euclidian distance: 11114967.6 m
      great circle distance: 1598085.6 m
      curvature distance: 1816146.3 m

The javawa tool reports

Distance:	1603.316 km

Assuming the javawa tool is correct, this is a 0.4% error for great circle and 13.3% error for curvature distance.

[cmuelle8]

However, given that getCenter() is not guaranteed to be in the exact middle, we must think how it effects

return eps > Math.abs(_ret/2-ret) ? ret : r(a, m, ret)+r(m, b, ret);

If javawa tool is correct which it seems to be, I've rechecked using ​http://www.ngs.noaa.gov/cgi-bin/Inv_Fwd/inverse2.prl, then curvatureDistance() is not yet - most probably because of the above line. The algorithm presumably can't use eps > Math.abs(_ret/2-ret) as a recursion stop condition, if the distance results of r(a, m, ret) and r(m, b, ret) are not equally distributed. We need to replace _ret/2 by the correct fraction calculated by each of the former.

[cmuelle8]

Complete output of ​http://www.ngs.noaa.gov/cgi-bin/Inv_Fwd/inverse2.prl for reference

Ellipsoid : GRS80 / WGS84  (NAD83)        
  Equatorial axis,    a   =    6378137.0000
  Polar axis,         b   =    6356752.3141
  Inverse flattening, 1/f =  298.25722210088
   
  First  Station : A                             
   ---------------- 
    LAT =  81  2 19.03200 North (N81.03862)
    LON =  79  6  5.61600 West  (W79.10156)
   
  Second Station : B                             
   ---------------- 
    LAT =  80 10 43.35600 North (N80.17871)
    LON =  20 44 31.88400 East  (E20.74219)
   
  Forward azimuth        FAZ =  42 40 26.9709 From North
  Back azimuth           BAZ = 321 45 24.9975 From North
  Ellipsoidal distance     S =   1603316.3096 m

[Don-vip]

Can you please add in the patch a unit test that checks this particular distance? Thanks.

[bastiK]

Taking two points of the same latitude, the curvature distance algorithm travels along a line of constant latitude. Lines of constant latitude are not great circles, so the actual center would be more towards the pole.

Btw., I'm grateful you brought this up, because our best distance function (greatCircleDistance) is indeed too simple and it would be good to have a more precise method.

[cmuelle8]

Replying to bastiK:

Taking two points of the same latitude, the curvature distance algorithm travels along a line of constant latitude. Lines of constant latitude are not great circles, so the actual center would be more towards the pole.

Yep, it needs a better "getCenter()"..

EDIT: It would only work on true-direction projections, but mercator is true-angle only. If recursion is done, it does not travel along great circle distances in most cases where a difference in longitude is involved. So while we can measure an angle distance between two points using mercator, it is not guaranteed that this is an angle distance on a great circle. Having pulled in "getCenter()" was a mistake, as it traces the mercator angle distance, which is not the shortest angle distance on the ellipsoid/sphere.

[bastiK]

You can try, but this may be difficult. Geotools uses GeographicLib for distance measurement, so I assume this is well tested: ​net.sf.geographiclib.Geodesic (line 526)

I haven't looked at it in detail, so it may or may not be suitable to be included in JOSM.

[cmuelle8]

I've come to the conclusion that it is fixable.

• in 3D space there exist planes intersecting both points (a) and (b) on the ellipsoid's surface
• using a third point (c), let it be the center of the ellipsoid, a definite plane P out of the former is chosen

• all of the following are guaranteed to lie in plane P
  • a non-spherical, possibly ​non-isosceles triangle given by (a), (b), (c)
  • a "great ellipse" for (a) and (b) which we do not know an ​eccentricity for
  • a ​chord of that ellipse given by triangle side (a b)

• "​great circles" have a zero eccentricity and constant circumference for arbitrary choices of (a) and (b)
• eccentricity and circumference of the "​great ellipse" are a function of (a) and (b) and thus varying

For the recursion of curvatureDistance to work flawlessly, partial sums must calculate lengths of
"great ellipse"-segments which lie in plane P. But actual ratio of partial lengths is not important.

A bisecting ray br from (c) intersecting the chord will also intersect the "great ellipse".

br will not intersect at the middle of the chord if sides (a c) and (b c), the "radii" at (a) and (b),
are of different length (mostly if lat(a)!=lat(b) and (lat(a)+lat(b))/2!=0). Similar applies
to the intersection point of br with the "great ellipse".

Let the recursion step be return thres > ret ? ret : r(a, m, ret)+r(m, b, ret);,
then using br to simply find any point guaranteed to lie on the "greate ellipse" will suffice.

thres is not eps - it needs to be chosen in a different way, reasonably large,
acknowledging that for results below a certain threshold distance, all true distances
of arbitrarily chosen (a) and (b) do not differ from thres distance more than eps.

[cmuelle8]

For someone potentially jumping in on this, ​Latitude and in particular ​#Latitude_and_coordinate_systems might be a good read.

Note that OSM generally deals with geodetic latitude and geodetic points that do not record h component,
which is the distance of points p on earth's surface and pr on the reference ellipsoid's surface,
such that vector p-pr is a ​surface normal on the ellipsoid.

Prolonged surface normals of geodetic points do not generally intersect with the ellipsoid's center (origin),
exceptions being points with latitudes of exactly 0°, 90° or -90°. However, all surface normals will intersect

• the axis running through the poles (minor axis), with the intersection point possibly lying "outside" of the ellipsoid of revolution
• a circle in equatorial plane (rotated major axis), centered at the origin, "drawn" by the focal points of the rotated ellipse

For points with fixed longitude and a variable latitude in ]-90°,0°[ or ]0°,90°[ the ellipse normals intersect the major axis in ]0,f1[ with f1 being one of the focal points. Repeating this example with an offset of 180° to the longitude, normals will interset ]f2,0[ with f2 being the other focal point of the same ellipse. Both focal points are equidistant wrt center at [0, 0].

Geodetic latitudes measure the angle at the intersection points of the ellipsoid normals with said circle in equatorial plane.

Different reference ellipsoids for different sections of the geoid may displace their origin from the geoid center. To decode LatLon data the ​geodetic datum or system that was used for encoding has to be known. A geodetic datum not only defines a reference ellipsoid, but adds at least two reference points to define origin and global rotation of a coordinate system "wrapping" its surface.

The ​geodetic datum WGS84 defines an ellipsoid globally distributed by the prominence of GPS and the data output of its consumer receivers. The reference points used to define the coordinate system were the poles and Greenwich for historic reasons: Origin was/is defined to be the "middle" of an arc running from pole to pole through Greenwich.

It is noteworthy that any LatLon coordinate system uses two-dimensional coordinates only. Coordinates of the third dimension are indirectly given by a function of LatLon. To reconstruct a three dimensional point, h over pr needs to be known in addition to the geodetic datum with its associated reference ellipsoid.

Geocentric latitudes are a component of spherical polar coordinates and measure the angle between equatorial plane and a line from p to geoid center for any p, "bypassing" any ellipsoid definition.

These lines are always surface normals to a sphere around geoid center with a radius equal to euclidian distance r of that line, but only in rare cases surface normals to the geoid or one of its reference ellipsoids.

A geocentric coordinate of p missing the r component is almost useless to reconstruct p: Using earth's mean radius will preserve direction but produces a point above or below p, unless r matched mean radius. Interpreting the coordinate using a reference ellipsoid skews the ray, it will not preserve original direction from geoid center to p, but produce a point near the original of p on earth's surface.

GPS receivers convert a measured geocentric coordinate of p to a geodetic point with WGS84 datum before output.

ps: Note that ellipsoids of revolution are also called spheroids (which I personally find misleading)

pps: ​The GPSTK user reference elaborates on geocentric and geodetic coordinates and the conversion process, pp.23-26,69

[bastiK]

Replying to cmuelle8:

For the recursion of curvatureDistance to work flawlessly, partial sums must calculate lengths of
"great ellipse"-segments which lie in plane P. But actual ratio of partial lengths is not important.

No, for example for two points on the equator, the shortest path is not to run along the equator! Due to the flattening of the earth, you can take a shortcut by making a slight curve towards one of the poles. :)

See ​Geodesics on an ellipsoid.

Btw., your method will work reasonably well, but why put so much effort into something when the result is still not perfect. The problem is already fully solved in theory and there are even well tested implementations we can copy.

[cmuelle8]

Replying to bastiK:

Replying to cmuelle8:
See ​Geodesics on an ellipsoid.

Btw., your method will work reasonably well, but why put so much effort into something when the result is still not perfect. The problem is already fully solved in theory and there are even well tested implementations we can copy.

Because copying alone is not understanding. Thanks for the link above, I've read fragments of it from time to time, but obviously not grasped it to the full extent yet - maybe I should take the time to fully read it on a quiet occasion. :)

EDIT: That said, if you need a copy of an implementation in LatLon now, noone will stop you importing an existing solution..

[cmuelle8]

Another reason to explore this in great detail could be the performance issue - the more accurate we can be with a fast approximation, the less cpu cycles might be wasted when a lot of primitives are selected, or ways have many segments..

Also think of the next step: Calculating areas. The formula for calculating polygons is easy - just add up the area of the minimal bbox around every segment and divide by two. Add outer polygons, subtract inner ones. Transported to bboxes on curved surfaces the same can be done - but solving lots of definite elliptic integrals for every bbox might cost a lot of cpu time.

So yes, in theory the exact case is solved, but ways to balance between those expensive methods and good approximations that reduce costs in practical applications I find to be sparse. We have already "highlight" on objects you hover the mouse over, if we find methods of inexpensive, good approximations we can enable distance and area calculation for each object hovered over (ideally caching the result, like you cache styles, or like the multipolygon-cache caches areas)..

The search is (in the very long run), one for iterative/recursive methods that can be interrupted at any time, such that the calculation is not lost but gives a usable approximation. E.g. for the area calculations it'd be nice to have an algorithm in the end that starts on all segments at the same time and if interrupted early, outputs the best approximation that could be made in the time it was permitted to run. Of course, it'd be nice to know (and eventually deliver that to the user as well) the error class the result is in..

[bastiK]

Fair enough, but be aware that ​Lambert's formula is rather simple and only a minor extension of the current greatCircleDistance() function. It has a global error of about 10 m and your algorithms should beat that in accuracy or performance, otherwise I don't see the point.

I'm not familiar with area-calculation on an ellipsoid, but given how complicated a simple distance is, this might be hellishly difficult and really a job for an external library.

[cmuelle8]

Replying to bastiK:

Replying to cmuelle8:
No, for example for two points on the equator, the shortest path is not to run along the equator! Due to the flattening of the earth, you can take a shortcut by making a slight curve towards one of the poles. :)

This is correct, but not fully relevant to our use case: You will find two disjunct geodesics on the ellipsoid for any arc segment of a great ellipse through points a and b, not only if they both lie on equator, but also if they are equally displaced for instance.

However, we do not define an osm way between a and b by the possibly two shortest arcs between them, but by a single arc. That single arc is defined by the shortest path on the great ellipse, not by any of the possibly two true shortest arcs. For some subset of a and b the different arcs coincide.

So, for your example, with respect to our application we measure a distance for, it would be wrong to associate the true shortest distance between a and b on ground with the distances of the ways drawn between a and b in JOSM. We indeed need to run along the equator to measure a meaningful distance of the line shown, defining a single arc on ground only.

a and b are associated a single "great ellipse" in any case, but the geodesic is not disamiguous for all cases.

[cmuelle8]

Visualizing the two true (or more in some cases) shortest arcs between two points on ground, we'd actually have an application for the true shortest distance. Such a visualization will need to paint arcs in the default mapview.

The "great ellipse" arc will deviate from the "great circle" arc as well, but applying this correction to the lines drawn in the mapview, which are projections of "great circle" arcs, will presumably not be noticeable.

Btw, a nice understanding of the problem can be obtained by exploring the intersection point IP of the surface normals. Combining comment:20 and comment:21, it can be seen that if IP does not exist (the normals do not intersect), then there is a maximum of three planes P1, P2 and P3 defined by points (a), (b) and

• the center of the polar axis, denoted as (c) in comment:20
• the intersection of the ellipsoid surface normal in (a) with the polar axis
• the intersection of the normal in (b) with the polar axis

The intersection of the planes with the ellipsoid yield a maximum of three different arcs.

If the arcs obtained by intersecting the last two planes are unequal, I suggest these to enclose the geodesic.

If the normals do intersect, and IP is neither in equatorial plane nor on polar axis, all three arcs will be equal and be the true shortest arc between the points. The same applies if IP lies uncentered in equatorial plane, but at least one of the normals lies not.
(The great ellipse through the points coincides with the rotated ellipse defining the ellipsoid in these cases.)

If both normals lie in equatorial plane, IP coincides with the center of the ellipsoid. Other two points, displaced from the center along the polar axis in proportion to the eccentricity need to define planes P2, P3 in this case. There will be two geodesics and two true shortest arcs left and right to equator in this case. Similar applies if IP lies on polar axis.

Put unspecific, first is the general case, second for points with no difference in longitude and third for those with no difference in latitude.

[cmuelle8]

Replying to bastiK:

You can try, but this may be difficult. Geotools uses GeographicLib for distance measurement, so I assume this is well tested: ​net.sf.geographiclib.Geodesic (line 526)

I'm currently reading up on the ​paper to GeographicLib by Karney which seems to be a more exhaustive resource than ​Geodesics on an ellipsoid regarding the math involved. It claims an accuracy for geodesic calculations on ellipsoids close to machine precision at 15 nm using double floats.

Unfortunately the introductory sentence might be misinterpreted, if a reader is biased on distance between two points:

A geodesic is the natural “straight line”, defined as the line of minimum curvature, for the surface of the earth (Hilbert and Cohn-Vossen, 1952, pp. 220–222).

A geodesic is a line of minimum curvature would fit the case for ellipsoids much better. Having two disjunct geodesics between two points is among the common cases. Having infinite geodesics is true for the poles.

It also has a nice Appendix, pp.24 continuing. E.g. Appendix C: Area of a spherical polygon. Or Appendix A: Equations for a geodesic, quote:

Laplace (1829,Book 1,§8) shows that the path of a geodesic on a surface is the same as the motion of a particle constrained to the surface but subject to no external forces.

---

I'm +1 on having Lambert's formula in JOSM code, but IMO it won't replace the need for something iterative/recursive that operates on ever reduced partial segments. E.g. if we had a function unlike getCenter(), one that indeed does return a center (or near center) LatLon coordinate lying on a "great ellipse" (if true distance matters: on a geodesic), we can dream about approximating geoid distances as well (by combining it with the srtm data, or by points having 'ele' tags in the db): We can't do such things using "plain" distance functions, I suppose. Of course, this won't happen overnight and unit tests along that path, may we choose it, should not be forgotten.

Lastly, lets take note of one very interesting thing about geodesics on the ellipsoid: If there is more than one geodesic connecting two points they are very likely to have different distances on the geoid, i.e. using the geoid we'll find out which of the geodesic arcs is the true shortest of shortest arcs. This won't matter at all for aerial or marinal navigation though.. It's more food for thought.

[cmuelle8]

Replying to Don-vip:

Can you please add in the patch a unit test that checks this particular distance? Thanks.

I will once we have a basic understanding of what to measure how in which cases :). There are quite some disambiguities.

[bastiK]

Replying to cmuelle8:

Unfortunately the introductory sentence might be misinterpreted, if a reader is biased on distance between two points:

A geodesic is the natural “straight line”, defined as the line of minimum curvature, for the surface of the earth (Hilbert and Cohn-Vossen, 1952, pp. 220–222).

A geodesic is a line of minimum curvature would fit the case for ellipsoids much better. Having two disjunct geodesics between two points is among the common cases.

I don't think it is very common. There is generally only one geodesic connecting two points, except when both points are exactly on the equator and possibly for pairs of antipodal points.

Replying to cmuelle8:

You will find two disjunct geodesics on the ellipsoid for any arc segment of a great ellipse through points a and b, not only if they both lie on equator, but also if they are equally displaced for instance.

I don't think so.

However, we do not define an osm way between a and b by the possibly two shortest arcs between them, but by a single arc. That single arc is defined by the shortest path on the great ellipse, not by any of the possibly two true shortest arcs.

As far as I know, OSM has no specification of how exactly to connect two points in a way. I think with typical segment lengths of < 1 km and unrecorded elevation differences, such a specification would be somewhat pointless.

a and b are associated a single "great ellipse" in any case, but the geodesic is not disamiguous for all cases.

Neither is the great ellipse: two paths for antipodal points and infinitely many for the poles.

---

I'm +1 on having Lambert's formula in JOSM code, but IMO it won't replace the need for something iterative/recursive that operates on ever reduced partial segments. E.g. if we had a function unlike getCenter(), one that indeed does return a center (or near center) LatLon coordinate lying on a "great ellipse" (if true distance matters: on a geodesic), we can dream about approximating geoid distances as well (by combining it with the srtm data, or by points having 'ele' tags in the db): We can't do such things using "plain" distance functions, I suppose. Of course, this won't happen overnight and unit tests along that path, may we choose it, should not be forgotten.

It would be fun to have a draw mode that paints the connection between two points not as a straight line, but interpolates it to a geodesic or great ellipse curve. :)

[cmuelle8]

Replying to bastiK

Neither is the great ellipse: two paths for antipodal points and infinitely many for the poles.

Still, a lot less special cases than depending on geodesics. Simply excluding way segments terminating at antipodal points will do. You will still be able to define ways to those by simply using two segments, which also solves the pole case.

Ruling out segments on equator seems strange.

​Pittmans method deserves to be refered to as well. The paper in the link has a very nice figure including non-intersecting normals and the geodesic between two points on page 2. It also explains intersecting planes of interest, among them the ones containing the normals in the terminating points I had described. It relates their intersection curves to the geodesic on page 3.

[bastiK]

Replying to cmuelle8:

Replying to bastiK

Neither is the great ellipse: two paths for antipodal points and infinitely many for the poles.

Still, a lot less special cases than depending on geodesics. Simply excluding way segments terminating at antipodal points will do. You will still be able to define ways to those by simply using two segments, which also solves the pole case.

Ruling out segments on equator seems strange.

You don't rule them out, but have to break the symmetry by an arbitrary choice (e.g. always go south). Geodesics are a very fundamental and natural concept, whereas other intersection curves are artificial and ad hoc choices with no useful property other than "slightly larger than the geodesic" and "possibly easier to compute".

​Pittmans method deserves to be refered to as well. The paper in the link has a very nice figure including non-intersecting normals and the geodesic between two points on page 2. It also explains intersecting planes of interest, among them the ones containing the normals in the terminating points I had described. It relates their intersection curves to the geodesic on page 3.

Neat!

[cmuelle8]

Replying to bastiK:

Replying to cmuelle8:

You will find two disjunct geodesics on the ellipsoid for any arc segment of a great ellipse through points a and b, not only if they both lie on equator, but also if they are equally displaced for instance.

I don't think so.

Imo you should think so: Given two non-antipodal, disjunct points on equator, let A be the area between the two possible geodesics, such that the shortest arc on equator is contained in A and halving A. Since the ellipsoid surface is continuous, A must shrink to zero first if both points are displaced gradually into different "hemiellipsoids" (defined by equator), for the geodesic to eventually become unique.

EDIT: You can also find an access to it like this:

Let A, B be non-antipodal, disjunct points on equator.
Let N, S be disjunct intersection points of polar axis and ellipsoid.

Then the geodesic on the northern hemiellipsoid between A, B defines an ellipsoidal triangle, together with the meridian arcs from N to A and N to B.

Another, non-intersecting ellipsoidal triangle with disjunct sides is defined by the geodesic on the southern hemiellipsoid between A, B, together with meridian arcs from S to A and S to B.

Let 2G be the bi-gon between the two geodesics.

We have essentially tri-sected the area between the two meridians defined by N-A-S and N-B-S. What will happen to 2G, if A and B are moved along their meridians, in opposite direction?

[cmuelle8]

Replying to bastiK:

You don't rule them out, but have to break the symmetry by an arbitrary choice (e.g. always go south). Geodesics are a very fundamental and natural concept, whereas other intersection curves are artificial and ad hoc choices with no useful property other than "slightly larger than the geodesic" and "possibly easier to compute".

What about "natural straight line" and guaranteed to lie between geodesics?
I think this would be a very useful property.

Using the surface normals once again we can find a definition of that path without presuming an intersecting plane. Something like: 'A natural straight line is given by an arc of maximum curvature on an ellipsoid.' For a bi-gon bound by two arcs with minimum curvature, I suppose this yields an great ellipse arc terminating at the gons. And because such an arc will repeat after a full revolution, I think it is a very natural concept nontheless, aside from uniqueness.

[bastiK]

Replying to cmuelle8:

Replying to bastiK:

Replying to cmuelle8:

You will find two disjunct geodesics on the ellipsoid for any arc segment of a great ellipse through points a and b, not only if they both lie on equator, but also if they are equally displaced for instance.

I don't think so.

Imo you should think so: Given two non-antipodal, disjunct points on equator, let A be the area between the two possible geodesics, such that the shortest arc on equator is contained in A and halving A. Since the ellipsoid surface is continuous, A must shrink to zero first if both points are displaced gradually into different "hemiellipsoids" (defined by equator), for the geodesic to eventually become unique.

You are right, I just considered displacement to the same side.

Replying to cmuelle8:

Replying to bastiK:

You don't rule them out, but have to break the symmetry by an arbitrary choice (e.g. always go south). Geodesics are a very fundamental and natural concept, whereas other intersection curves are artificial and ad hoc choices with no useful property other than "slightly larger than the geodesic" and "possibly easier to compute".

What about "natural straight line" and guaranteed to lie between geodesics?
I think this would be a very useful property.

Using the surface normals once again we can find a definition of that path without presuming an intersecting plane. Something like: 'A natural straight line is given by an arc of maximum curvature on an ellipsoid.' For a bi-gon bound by two arcs with minimum curvature, I suppose this yields an great ellipse arc terminating at the gons.

I don't get what you mean by maximum curvature. The curvature is different at each point along the curve. You could take some kind of average along the line, but then it is far from obvious why this would give great ellipses.

I was exaggerating a bit by saying that curves other than geodesics are useless. We can have multiple "distance" functions side by side. But so far I'm not convinced, that great ellipses are particularly special in the large selection of possible choices one can pick for connecting two points on a spheroid.

[cmuelle8]

Replying to bastiK:

You are right, I just considered displacement to the same side.

Btw, if abs(lon(a)-lon(b)) <= (1−f)*pi && lat(a)==0 && lat(b)==0
then equator is a shortest geodesic. I believe ​the wiki article to be true in that respect.

If a and b are moved ever closer together on equator, there must be a time when the geodesic north, the one south and equator all have minimum curvature (when the bi-gon surface collapses to an arc), and hence the same length. Beyond this step the area would grow in size and hence is not of interest anymore. It's like climbing a hill, once you reach a certain height, it's shorter crossing the top to descend to the same height you would otherwise reach on equidistant paths around it.

I don't get what you mean by maximum curvature. The curvature is different at each point along the curve. You could take some kind of average along the line, but then it is far from obvious why this would give great ellipses.

Not average, but sum. Take a look at the first quote in comment:30. (Of course, substituting this with maximum curvature, you'll end up with an path of infinite length, so a further constraint is needed. Like 'all surface normals of the path must lie in one plane'.)

[bastiK]

Replying to cmuelle8:

Replying to bastiK:

You are right, I just considered displacement to the same side.

Btw, if abs(lon(a)-lon(b)) <= (1−f)*pi && lat(a)==0 && lat(b)==0
then equator is a shortest geodesic. I believe ​the wiki article to be true in that respect.

So if you want to get from one point of the equator to another point on the equator, you should follow the equator line, exept if the other point is on the opposite side of the earth or less than 67 km apart from the opposite point. In that case you have two optimal routes, one in the northern and one in the southern hemisphere.

I will remember that for my travels. :)

I don't get what you mean by maximum curvature. The curvature is different at each point along the curve. You could take some kind of average along the line, but then it is far from obvious why this would give great ellipses.

Not average, but sum. Take a look at the first quote in comment:30. (Of course, substituting this with maximum curvature, you'll end up with an path of infinite length, so a further constraint is needed. Like 'all surface normals of the path must lie in one plane'.)

Alright, but I don't see where this is headed.

[bastiK]

Closing due to inactivity. As stated above, the ellipsoidal distance function can be improved, but the current version of the patch is not the right direction IMO.