Jump to content

 
Photo

OpenGL Texture Mapping [Equirectangular to Azimuthal Equidistant]

- - - - -

  • Please log in to reply
4 replies to this topic

#1
felixrulz

felixrulz

    Newbie

  • Validated Member
  • Pip
  • 6 posts
  • Australia

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:

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:

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:
Posted Image

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.

p_d[n].u = 0.5 + a2[i]/PI2; // + rotation

		p_d[n].v = (b2[i])/R;
This is the resulting image:
Posted 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....d...&pg=GBS.PR1,
the equations are reproduced here at Wolfram http://mathworld.wol...Projection.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:
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.

#2
Michael Schmeling

Michael Schmeling

    Master Contributor

  • Validated Member
  • PipPipPipPip
  • 206 posts
  • Gender:Male
  • Location:Kassel, Germany
  • Germany

The PROJ.4 - Cartographic Projections Library (written in C) may help you with the calculations.
Michael Schmeling
Kassel, Germany
Arid Ocean Map Illustrations
http://maps.aridocean.com
Indie Cartographer
http://www.indiecartographer.com

#3
felixrulz

felixrulz

    Newbie

  • Validated Member
  • Pip
  • 6 posts
  • Australia

The PROJ.4 - Cartographic Projections Library (written in C) may help you with the calculations.


Thanks Michael, I'll take a look


Edit: Hopfully this works, seems that none of the posts I have made are showing up. See if it will accept this edit...


I think I have solved most of the mapping problems. I ended up using the method derived by T.I here: http://somedayspace....ection-and.html.

After testing formulas in MALTAB I got promising results:
Posted Image

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:
Posted Image

I think the problem problem 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 distorted sectors are displaying a compressed version of nearly the entire map. What should happen instead is to map from A wrapping to B (to the right). 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.
Posted Image

But when I added this check:
if (p_d[n].u < p_d[n-1].u)

		p_d[n].u += 1;
only the problem from one of the polar aspects was fixed (and not at all from oblique aspects).

Any ideas? Here is a link to the program I have so far if it will help with diagnoses: http://dl.dropbox.co.....ery.close].7z
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.


[EDIT]
The problem is no longer cartography related so I i'm not sure If I should keep posting here.
I have solved the majority of the problems using the method I have illustrated in this short video:
http://dl.dropbox.co...[EvenCloser].7z

Edited by felixrulz, 28 July 2012 - 02:30 AM.


#4
felixrulz

felixrulz

    Newbie

  • Validated Member
  • Pip
  • 6 posts
  • Australia

Thanks Michael, I'll take a look


Ok, so there is a function for Azimuthal Equidistant: "PJ_aeqd.c" it seems to include many variations I didn''t even know about...
This seems to be the section i am interested in:

INVERSE(s_inverse); /* spherical */

	double cosc, c_rh, sinc;



	if ((c_rh = hypot(xy.x, xy.y)) > PI) {		[color="#FF0000"][b]// for some reason the radius cannot be larger than PI[/b][/color]

		if (c_rh - EPS10 > PI) I_ERROR;   [color="#FF0000"] [b]// EPS10 is a small (1.e-10) used to prevent problems with rounding errors[/b][/color]

		c_rh = PI;

	} else if (c_rh < EPS10) {						   [color="#FF0000"] [b]// if at the centre of the projection just use centre lat/long[/b][/color]

		lp.phi = P->phi0;

		lp.lam = 0.;

		return (lp);

	}

	if (P->mode == OBLIQ || P->mode == EQUIT) 

	{

		sinc = sin(c_rh);								[color="#FF0000"][b]//precalculate for speed[/b][/color]

		cosc = cos(c_rh);							   [color="#FF0000"] [b]//precalculate for speed[/b][/color]

		if (P->mode == EQUIT) 

		{

 			[color="#FF0000"][b]//calculate latitude[/b][/color]

				lp.phi = aasin(P->ctx,	 xy.y * sinc / c_rh);

			xy.x *= sinc;

			xy.y = cosc * c_rh;

		}

		else 

		{

 			[color="#FF0000"][b]//calculate latitude

			I used: latitude = asin( cos(radius)*sin(lat0) + (y*sin(radius)*cos(lat0)/radius) );

			a bit different to what they have, I'll try to find where they define aasin (haven't been able to yet)

			[/b][/color]

			lp.phi = aasin(P->ctx,cosc * P->sinph0 + xy.y * sinc * P->cosph0

				c_rh);

			xy.y = (cosc - P->sinph0 * sin(lp.phi)) * c_rh;

			xy.x *= sinc * P->cosph0;

		}

[color="#FF0000"][b]//calculate longitude

			for the oblique case I used: 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]) ) );

			[/b][/color]

		lp.lam = xy.y == 0. ? 0. : atan2(xy.x, xy.y);

	} 

	else if (P->mode == N_POLE) 

	{

		lp.phi = HALFPI - c_rh;

		lp.lam = atan2(xy.x, -xy.y);

	} 

	else

	{

		lp.phi = c_rh - HALFPI;

		lp.lam = atan2(xy.x, xy.y);

	}



	return (lp);

I'll try and work out why their function is different.

Cheers

#5
felixrulz

felixrulz

    Newbie

  • Validated Member
  • Pip
  • 6 posts
  • Australia

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....ection-and.html.

After testing formulas in MALTAB I got promising results:
Posted Image

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:
Posted Image

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.
Posted Image

I added a check:
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: 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.



Problem solved:
See hereProblem Solved

Not sure why but new replies won't work so i added info in this one




0 user(s) are reading this topic

0 members, 0 guests, 0 anonymous users

-->