OpenGL Texture Mapping [Equirectangular 2 Azimuthal Equidist

[felixrulz]

Hey, this is my first time here... was looking for some help from some cartography savy people

I have created an almost spherical disk by joining lots of triangular/rectangular polygons together and have a equidistant image of the earth loaded as a texture map. What I am having trouble doing is converting/mapping the azimuthal equidistant points on the disk to a latitude and longitude so I can determing the texture coordinate. Here is a snippet of code where the disk is created:

Code: Select all

void CreateDisk (double R, double x0, double y0, double z0)
{
	const double rad_inc = R/(disk_r_divs);
	const double theta_inc = 2*PI/(disk_az_divs);
	double theta, radius, b2[4], a2[4];
	double latitude, longitude;
	int n=0,i; // Current vertex

	for(radius=0; radius < (R+rad_inc/2); radius+=rad_inc)
	{
		b2[0] = radius;
		b2[1] = radius + rad_inc;
		b2[2] = radius;
		b2[3] = radius + rad_inc;

		for(theta=0; theta <= PI2; theta+=theta_inc)
		{
			a2[0] = theta;
			a2[1] = theta;
			a2[2] = theta + theta_inc;
			a2[3] = theta + theta_inc;

			for (i=0; i<4; i++)
			{
				// disk coordinates
				p_d[n].x = b2[i] * cos(a2[i]) + x0;
				p_d[n].y = b2[i] * sin(a2[i]) + y0;
				p_d[n].z = z0;

				// TEXTURE COORDINATES CALCULATED HERE

				n++;
			}
		}
	}
}

The use of b2[] and a2[] are to calculate the other points of the rectangle. And are not really relavent to the problem I am havving (I think) so this code can be interpreted like this for simplicity:

Code: Select all

void CreateDisk (double R, double x0, double y0, double z0)
{
	const double rad_inc = R/(disk_r_divs);
	const double theta_inc = 2*PI/(disk_az_divs);
	double theta, radius;
	double latitude, longitude; //used to calculate texture coordinates
	int n=0 // Current vertex

	for(radius=0; radius < R; radius+=rad_inc)
	{
		for(theta=0; theta <= PI2; theta+=theta_inc)
		{
				// disk coordinates
				p_d[n].x = b2[i] * cos(a2[i]) + x0;
				p_d[n].y = b2[i] * sin(a2[i]) + y0;
				p_d[n].z = z0;

				// TEXTURE COORDINATES CALCULATED HERE

				n++;
		}
	}
}

This code starts with a radius of 0 and then will sweep through theta (the azimuth of the map). Repeating with larger radii untill reaching R. Here is an example of the resulting disk:


The problem I am havving is converting the cartesian (or polar if this is easier) coordinates to a latitude and longitude for the equirectangular map (This needs to be inserted where indicated in the code sample). I have managed managed a polar aspect by simply using a ratio of the radius and azimuth. (Note that the +0.5 and /PI2 and /R are to convert/normalise the latitude and longitude into a 0:1 range which is required by OpenGL.

Code: Select all

		p_d[n].u = 0.5 + a2[i]/PI2; // + rotation
		p_d[n].v = (b2[i])/R;

This is the resulting image:


For the oblique projection I am struggeling. The latest equations I have been trying to use where taken from Map Projections – A Working Manual (pg 196) http://books.google.com.au/books/reader ... pg=GBS.PR1,
the equations are reproduced here at Wolfram http://mathworld.wolfram.com/AzimuthalE ... ction.html.

I implemented them in code as follows. b2[] can simply be considered a radius and a2[] an azimuth, p_d[n].y/x are the x/y coordinates:

Code: Select all

				if (radius == 0)
				{
					latitude = lat0;
					longitude = long0 + theta;
				}
				else
				{
					latitude = asin( cos(b2[i])*sin(lat0) + (p_d[n].y*sin(b2[i])*cos(lat0)/b2[i]) );

					if (lat0 == PI_2)
						longitude = long0 + atan( -p_d[n].x/p_d[n].y );
					else if (lat0 == -PI_2)
						longitude = long0 + atan( p_d[n].x/p_d[n].y );
					else
						longitude = long0 + atan( p_d[n].x*sin(b2[i])/( b2[i]*cos(lat0)*cos(b2[i])-p_d[n].y*sin(lat0)*sin(b2[i]) ) );
				}

				// convert lat/long to texture coordinates
				p_d[n].u = longitude/PI2;
				p_d[n].v = 0.5 + latitude/PI;

The results seem somewhat correct but I don't even get the entire azimuthal projection (for example with australia at the centre you can't see much else). I have been trying to fix this for near a week and thought I should ask for help.

Thanks.

[daan]

Hello, felixruiz, and welcome to the forums.

Your problem here is in your interpretation of the map projection formulae. You write, for example,

Code: Select all

latitude = asin( cos(b2[i])*sin(lat0) + (p_d[n].y*sin(b2[i])*cos(lat0)/b2[i]) );

but where you have

Code: Select all

cos(b2[i])*sin(lat0)

what you actually need is cos(c)*sin(lat0). c is computed as √(x²+y²). Likewise, you divide by

Code: Select all

b2[i]

but you need to divide by c.

Best of luck.
— daan

[felixrulz]

hey daan,

Thanks for your help but i am a little confused as to the problem you mentioned. In my code b2[] represents the radii for the current quadrilateral... since (by pythag) sqrt(x^2 + y^2) = R, in my code I just used R as the for-loops are generating polar coordinates.

[this is not my image - credit to 17calculus.com]

Am I overlooking something?

Thanks, Felix

[daan]

I see now what you intend by b2. Apologies for the confusion. I’m curious why your subscripts for b2 differ from your subscripts for p_d: You use i for b2 and n for p_d.

One problem for sure is that your atan function only yields a range of ±π/2 when you need a range of ±π. You should be using an atan2 implementation:

Code: Select all

  if (lat0 == PI_2)
    longitude = long0 + atan( -p_d[n].x/p_d[n].y );
  else if (lat0 == -PI_2)
    longitude = long0 + atan( p_d[n].x/p_d[n].y );
  else
    longitude = long0 + atan2( p_d[n].x*sin(b2[i]), ( b2[i]*cos(lat0)*cos(b2[i])-p_d[n].y*sin(lat0)*sin(b2[i]) ) );

Good luck!
— daan

[felixrulz]

I see now what you intend by b2. Apologies for the confusion.

I think some of the confusion can be attibuted to my bad variable names...

I’m curious why your subscripts for b2 differ from your subscripts for p_d: You use i for b2 and n for p_d.

I'll try to explain what I am doing. It probably isn't the best method as I am storing the same points several times...
[Note: I adapted this method from a tutorial written by swiftless]

As in the code segment below for each point on the disk, four points defining a quadilateral are calculated (using the offsets). So realy b2 and a2 are parallel arrays (or should be stored in the same array). These four points have been (hastily) drawn in the pic below:

Code: Select all

	
for(radius=0; radius <= R; radius+=rad_inc)
	{
		for(theta=0; theta <= PI2; theta+=theta_inc)
		{
			// Calculate polar coordinates for the QUAD
			b2[0] = radius;
			b2[1] = radius + rad_inc;
			b2[2] = radius;
			b2[3] = radius + rad_inc;

			a2[0] = theta;
			a2[1] = theta;
			a2[2] = theta + theta_inc;
			a2[3] = theta + theta_inc;

			for (i=0; i<4; i++)
			{
                               // STORE POINTS OF QUAD
				// ('i' used for both b2 and a2)
			}
		}
	}

the for loop using the 'i' variable is used for both b2 and a2 representing the polar coordinates. p_d (points of disk) then stores all the coordinates (and eventually the texture coordinates). 'n' is used as a counter for the number of points stored.
the p_d structure will then be used in a drawing function using :

Code: Select all

glBegin (GL_QUAD_STRIP);
	glVertex3f (p_d[].x, p_d[].y, p_d[].z);
glEnd();

-------------------------------------------------------

One problem for sure is that your atan function only yields a range of ±π/2 when you need a range of ±π.

Ahhh, i thought I may have had a problem with the signs in different quadrents. I wrote a matlab script and I was getting half the texture coordinates overlapping. I'll take a look at this, thanks.

[felixrulz]

Ok, so I think I have solved most of the mapping problems. I ended up using the method derived by T.I here: http://somedayspace.blogspot.com.au/201 ... n-and.html.

After testing formulas in MALTAB I got promising results:


I then converted to a function in C and tested. I seem to be very close now; there is just some detracting distortion where the image is joined along the date line:


I the problem (at least I think) is due to the wrapping of texture coordinates. I have tried to illustrate the problem in the figure below. The program is currently mapping from A-->B so the sector has a compressed version of nearly the entire map; what we want is to map from A to the right wrapping to B. In OpenGL I think this can be achieved by enabling GL_REPEAT and simple using that B'=B+1 to get the correct texture response.


I added a check:

Code: Select all

if (p_d[n].u < p_d[n-1].u)
        p_d[n].u += 1;

but this only fixed them problem from one of the polar aspects.

Let me know if you have any ideas on hot to fix this. Here is a link to the program I have so far: http://dl.dropbox.com/u/49415695/Deskto ... lose%5D.7z if it helps diagnose (or if you just want a look).
You will need the data folder with the image (and maybe some dll files I haven't included). The arrow keys control the centre of the projection, 't' toggles the texture, 'up-pg' 'dn-pg' zoom in and out.

[daan]

Well, at this point you are firmly out of the map projections realm and into the vagaries of OpenGL. My best guess is that you need to adjust for the central meridian by rotating the entire rectangular image. But, just a guess.

Best of luck.
— daan

[felixrulz]

Well, at this point you are firmly out of the map projections realm and into the vagaries of OpenGL.
— daan

I will go searching some OpenGL forums. Thanks for your help.
Should I countinue posting developments here or is that not advised?

[daan]

I would be interested to know what resolves your trouble.

Best of luck.
— daan

[felixrulz]

I have nearly fixed the problem by searching for jumps (above some threshold) in the texture coordinates. I have made this short video to illustrate the method I have used: [video]http://dl.dropbox.com/u/49415695/Mappin ... %20Fix.mp4[/video]

The only problem is that the problem still appears, but after you move the centre point around (using the arrow keys) the problem appears to vanish (it takes a little while). Even if I pause the program and reenter the original centre where I previously had distortion the problem dosn't appear. I'm a little baffled. I think i'll export to a .txt file and have a closer look in another program.
Here is the new version of my program: http://dl.dropbox.com/u/49415695/Deskto ... oser%5D.7z