[Proj] Some ideas about the general inverse method

[Glynn Clements]

support.mn at elisanet.fi wrote:

> the basic idea might be that if the forward projection is available,
> it is always possible to find the inverse just by calling the forward one,
> by probing, until required amount of accuracy is achieved. 
> 
> The basic idea is the same as in binary search algotrithm.. or Newton
> iteration algotrithm.
> 
> 1) Assume that the requested point is exactly in the middle of the projection
> space (by very good luck). Calculate forward point.
> 
> 2) if not... in wihch part of the space it is.Here we compare the
> coordinates we have against the coordinates given by the user.
> Or use Newton methods, which are a bit similar.
> 
> 3) At this point we can exclude 75% of the remaining area since the
> point must be in some quarter space of the area and we know where.
> Divide lattitudes and longitudes so that the point is always in the
> halve which we are going to use in next iteration.
> 
> 4) Now repeat the assumption that the point must be in the middle
> of the remaining area.
> 
> 5) Loop until required amount of accuracy is achieved.

This method only works if the function is injective (two distinct
inputs never produce the same output, i.e. no "folding").

Although this is an important property for any cartographic
projection, it's only necessary for it to hold over the region for
which the projection is intended to be used. E.g. proj=alsk is a long
way from being injective once you get outside of Alaska.

In order for this method to work in general, you would need to be able
to reliably determine the initial rectangle; you can't just assume
180W-180E, 90S-90N.

Using Newton's method would require derivatives, which we don't have. 
And it would still require a reasonable starting point in order to
converge on the correct root.

> As we all know binary search is rather fast, since the area gets halved
> at every iteration. Especially if we can limit the area at the beginning
> strongly. We have for example some direct (hash) table, that already
> gives us a very small area to start with (not the entire world).
> 
> There might be several other ways to speed up that method...

The main one is to use proportional interpolation rather than a
midpoint. A couple of divisions is a trivial price to pay to reduce
the total number of pj_fwd() calls.

> The user could supply some limiting points where the following
> points must be. At that point the size of the tables and amount of
> search loops could be even smaller.

As mentioned above, this isn't just an optimisation. It's essential in
order to be able to find the correct root.

> Even if we first have just some version of the routine, even very
> slow and large... when we continue to improve it... finally we have
> something that is very fast, accurate and reliable.
> 
> The routine could use the parameter: "accepted amount of error"
> in meters for example. As soon as that is achieved, no more iterations
> is needed. People that need more accuracy, would also need more
> computer power (and time)...

Once you get to within a kilometre, the projection will be almost
affine, at which point proportional interpolation will converge very
quickly.

-- 
Glynn Clements <glynn at gclements.plus.com>