I've got a texture in openGL, which is basically a picture of a large
area of ground.

The coordinate system of the picture is a cartesian system
where the bottom left corner is (-1,-1) and the top right corner is
(1,1).  So the pixel at the center of the images is described by (0,0)
in texture coordinates.

I need a way of computing at runtime, a function G[x,y] = [lat,long]
where x and y describe any particular point on the picture.

The information I have to derive this function are several points on
the texture where I know both the texture coordinates and the lat,long
coordinates at that point.


My question is:

Assuming there are a minimally sufficient number of arbitrary points
available on the map ( I think you need a minimum of 4 ), how can I
derive the formula that lets me calulate (lat,long) from (x,y) for
any given map picture?

This is my first time using this service, and this seemed like a good
question to try it out with.

---

Clarification of Question by tsurugi-ga on 28 May 2005 20:17 PDT

The area of land is a map of a section of earth, so the latitude and
longitude lines will have curvature.  That's why you need (I think) a
minimum of four points.  The first attempt at solving the problem was
to create two polynomial functions:

lat = G1(x,y)
long = G2(x,y)

but that approach doesn't work because lines of constant latitude and
constant longitude aren't parallel to lines of constant x and constant
y, respectively.

So, the function to be generated will need two dependent variables and
two independent variables, like such:  [lat,long] = G(x,y)

ALSO:  I hadn't considered that particular sections of the globe (near
the poles)  would be much more complex to generate a function for, and
I'll be quite satisfied with an answer that assumes that the curvature
of longitude lines on the map sections are limited to second order. 
In other words, if it works for the United States, it works for me.  I
mean, a solution which took the poles into account would warrant a big
tip, but that's just because it'd be cool, not because I'd need it.

---

Request for Question Clarification by hedgie-ga on 31 May 2005 05:55 PDT

tsurugi-ga

   You are wrapping a section of a plane on a curved surface (sphere).
 This is inverse of a spherical projection (as discussed here)
http://astronomy.swin.edu.au/~pbourke/projection/spherical/

    There is always a distortion and solution is not unique 
- unless you specify which projection is used. 



so is  this the task?

given  table

  lat.i lon.i x.i y.i

for i= 1,2,..N >3

   what is an algorithm which will derive the formulas

lat(x,y)
lon(x,y)

  using stereographic (or just any reasonable) projection
  which will APROXIMATELY reproduce the table.

 Crux if the matter is this: You can minimize the error of the aproximation
 by using more complex projection, but you also would need larger N.

So - you have a choice here.


Hedgie

---

Clarification of Question by tsurugi-ga on 31 May 2005 14:00 PDT

I believe that you have the right idea. Unfortunately, we'll just get
bitmaps and reference points, we may not have any knowledge of the
projection that was used to generate the image. I understand that a
low number of reference points will yield a high uncertainty, but I
have to live with it. If we don't have enough reference points to get
any solution at all, then I can just reject the map. It would be nice
to have some idea of how many reference points are needed, but I
presume that would depend on the projection used. Clearly, we would
need at least 4, but getting 8 - 12 points is not unreasonable. I can
usually generate reference points by identifying features in the
bitmap, so I can control the point count and ensure non-orthogonality.

Exit criteria would be a function that will use all available
reference points to attempt to generate a mapping solution and an
estimate of the error, or determine that none is possible.

Thanks for your response and interest.

---

Clarification of Question by tsurugi-ga on 31 May 2005 14:10 PDT

I added another $10 for the error estimate

I take your clarification as an 'OK' to this formulation:

 given the table 4*N  of N rows, 4 colums (input table) 

x  y  lat lon

find the 'best' equation giving lat, lon for each x,y,
based on stereographic projection.

  Solution can be obtained by the Least Square algorithm,
  which will determine four parameters. 
Referring to Fig and terminology of this reference 

http://www.mathdaily.com/lessons/Stereographic_projection

 the equation given in the reference after lines

"Cartesian coordinates

There are a number of ways to perform stereographic projection onto a
sphere, based on your choice of where you put the plane and the
sphere. One can treat the plane as being complex numbers, and use the
following pair of transformations"

will give you your equation (the Z=u+iw, from  xsi eta zeta triple
corresponding point on the sphere, and so the lat and lon.

The unknown parameters can be found by least squares regression analysis
of which general overview is here:

http://www.statsoft.com/textbook/stpls.html

You need minimum four points, and you can use Stepwise regression to decide
what is a reasonable number.

A simple way of setting up your regression is to use LFT transform explained here:

http://www.mathpages.com/home/kmath464/kmath464.htm

You are seeking a,b,c,d, coefficients , so that 
Error  is minimized.   Errors is sum over the i of the
distance between projected (Calculated) and actual position on the plane.
For each i, Z is calculated form lan and lot in the table, using
projection eq. given above.

You compute transformed intput table

 X=x+iy  , Z= u+i*w    for each i

The problem can be linearized as follows 

 Z * (c *X +d) = (a *X +b) + error

where X = x+iy comes from your input table (actual point)
and Z is projection of the corresponding points (lat lon) on the same plane 

Given the quadratic expresson for the Error, software for determining
a b c d is widely available.

You add a constraint a*b +c*d =1 to limit the transform to rotations of the sphere
so you have only 4  independent parameters and so need at least N=4 equations.

This is just to start the ball rolling, and to see if I understand what you mean.

My understanding is:
>>>the large area is flat. It is not spherical.
>>>[lat,long] is [latitude,longitude].
>>>the coordinates of the point (x,y) in question are known. Only its
[lat,long] are to be solved.
>>>the other points [x's,y's] and [lats,longs] are known.

If that is correct, then you need only one other point to solve for
the [lat,long] of point (x,y).
Let the other point be (x1,y1),[lat1,long1].
And point (x,y) has the unknown [lat,long].

Then
lat = lat1 +(y -y1) 
lat = lat1 -y1 +y  ----***
Latitudes are horizontal lines, hence the y-coordinates.

long = long1 -x1 +x  ----***

---------
If a minimum of 3 other points (that means 4 points in all when (x,y)
is counted) are needed to solve for the [lat,long] of point (x,y),
then the coordinates of point (x,y) are not known?

Let us call point (x,y) as P.

P is a random point? Its coordinates, lat, and long are unknown? 

If yes, then 
>>>you must know the distance of P from each of the other 3 points to
solve for the coordinates of P. A formula can be derived if we know
the 3 distances.
>>>then, using any known (x1,y1)&[lat1,long1], and the just-solved
coordinates of P, you can get the [lat,long] of P.

So the area is spherical. It is a piece of the earth's surrface.

The points in the area are in cartesian coordinates and the area is
actually spherical.
Umm, too much computations. I pass.