PDA

View Full Version : OpenGL Texture Mapping [Equirectangular 2 Azimuthal Equidistant]

felixrulz
07-18-2012, 10:00 PM
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 theta_inc = 2*PI/(disk_az_divs);
double latitude, longitude;
int n=0,i; // Current vertex

{

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 theta_inc = 2*PI/(disk_az_divs);
double latitude, longitude; //used to calculate texture coordinates
int n=0 // Current vertex

{
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:
http://dl.dropbox.com/u/49415695/Disk.png

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:
http://dl.dropbox.com/u/49415695/Polar%20aspect.png

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?id=nPdOAAAAMAAJ&printsec=frontcover&output=reader&pg=GBS.PR1,
the equations are reproduced here at Wolfram http://mathworld.wolfram.com/AzimuthalEquidistantProjection.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:

{
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.

felixrulz
07-25-2012, 09:36 AM
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/2012/03/azimuthal-equidistant-projection-and.html.

After testing formulas in MALTAB I got promising results:
http://dl.dropbox.com/u/49415695/NearlingWorkingMappingMatlab.PNG

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:
http://dl.dropbox.com/u/49415695/distortion.png

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.
http://dl.dropbox.com/u/49415695/TextureWrappingProb.PNG

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/Desktop5%5BAzimuthalEqui-very.close%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.

Redrobes
07-26-2012, 05:03 AM
With a quick look its hard for me or anyone to be too precise about a solution but from a glance here I would say that n your first post you have atan( ... ) where the ... is some func of x/y but usually you use atan with some func of y/x. Whether that's right or not in this case tho I don't know.

From the second post you can use GL_REPEAT with OpenGL or you can mod the texture coord yourself. I think in your case changing the coord would be better. To do this for getting into the principle range from a -ve value you normally add on a repeat value and then call fmod( x ) where x is the repeating amount.

But looking at the image I dont think that is your issue. What it looks like to me is the gourand shader for opengl where it will linearly ramp the texture coord from one vertex to another. If you have texture coords that run from say 0 to 1 around the globe and your last face has vertexes with texture coord of 0.95 on the left and 0 on the right then that last face is going to shade the whole texture map around. What you need is probably a duplicated set of vertices on the seam with a texture coord of 1.0 so it maps back to the 0. I.e. you now have 0.95 on the left and 1.0 on the right.

If you set the normals to be radial then it should blend perfectly. If you use a texture map with the extra border pixel duplicated then it will antialias correctly as well.

Well thats my my best guess from looking at it.

felixrulz
07-27-2012, 03:49 AM
But looking at the image I dont think that is your issue. What it looks like to me is the gourand shader for opengl where it will linearly ramp the texture coord from one vertex to another. If you have texture coords that run from say 0 to 1 around the globe and your last face has vertexes with texture coord of 0.95 on the left and 0 on the right then that last face is going to shade the whole texture map around. What you need is probably a duplicated set of vertices on the seam with a texture coord of 1.0 so it maps back to the 0. I.e. you now have 0.95 on the left and 1.0 on the right.

Thanks for that Redrobes, I looked ito the problem and it is exactly that. Here are some pictures I made with MATLAB. I have put a screenshot my program from the same aspect next to it. The blue quadrilaterals are where the texture coordinates where seperated by more than some threshold (i.e. they have wrapped around).
http://dl.dropbox.com/u/49415695/1.PNG
http://dl.dropbox.com/u/49415695/2.PNG
http://dl.dropbox.com/u/49415695/3.PNG

I made an attempt at correcting by adding 1 to the smaller of the coordinates but this didn't seem to correct the problem, but I still need to try a few things.

Redrobes
07-27-2012, 05:03 AM
If your using the more normal and newer vertex arrays, then you have faces which share the vertices at the ends of the globe. Texture coords are mapped to vertices and the texture lookup engine interpolates between them. You cant just add or subtract from the values at the vertices because at the two different times you need different values. So you need another set of vertices with the same position coords but with different texture coords.

If your using old school OpenGL with immediate mode geometry vertices then you will be entering the texture coord on the fly. In that case you need to ensure you don't wrap the texture coord value back around to 0 or else it will reverse the interpolation but instead use GL_REPEAT and let the texture engine repeat mode do it for you.

felixrulz
07-28-2012, 03:24 AM
Hey Redrobes. I'm very new to OpenGL (this is my first program) so I'm not 100% sure what method I am using. As in the first code sample I posted, I populate an array of structures with the x-y and u-v coordinates. I then loop through them in a display function:

LOOP THROUGH POINTS
{
glTexCoord2f (p_d[index].u, p_d[index].v);
glVertex3f (p_d[index].x, p_d[index].y, p_d[index].z);
}
glEnd();

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:
https://www.dropbox.com/s/si9e0pc7dycr4mg/MappingWrapFix.mp4

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. Here is the new version of my program: http://dl.dropbox.com/u/49415695/Desktop6%5BEvenCloser%5D.7z (LINK)

The problem is no longer cartography related however... perhaps i should not continue posting here.

Redrobes
07-28-2012, 08:12 PM
OK thats old school GL. Now when you come to putting in the last row of texture coords then you need to ensure that the values your using in that last row are not the same as the first row even tho the vertex positions will be the same. I.e. if you have 11 points numbered in 0 to 10 inclusively then you would have pos[0] as 0, pos[1] as 1, pos[9] as 9 and then pos[10] as 0 so that it wraps around. But you need your texture coords to be tex[0] = 0, tex[1] = 1, tex[9] = 9, tex[10] = 10 and use GL_REPEAT.

Maybe stack overflow would be better suited now. Its pretty hard to debug apps like this without being able to step it through. Now that you have a good idea what is going on it should be quite quick to see where its occurring and then to fix. If you get stuck then print out the position and the tex coord in a table and see if when it wraps your not using a wrapped tex coord too so that the interpolator thinks ok ill just run right back through the whole texture map back to 0 in that last face.

felixrulz
07-29-2012, 10:06 PM
Now when you come to putting in the last row of texture coords then you need to ensure that the values your using in that last row are not the same as the first row even tho the vertex positions will be the same.

I don't think its quite that simple except in the polar aspects (image on left). From the image below you can see that the distortion is not always in the same location (image on right).
http://dl.dropbox.com/u/49415695/distortion.png
This is why I used the method described in the video above. (Which only partially worked).

Maybe stack overflow would be better suited now. Its pretty hard to debug apps like this without being able to step it through.

I'm not sure what you mean by "stack overflow"; I thought this was a bad thing when you ran out of space on the stack?
I've been writhing the algorithm section of my program in MATLAB first to give step by step visual cues. I was going to try outputing values to a file where I could analyse like you sugested

FILE *fp= fopen("test.txt","w");

Redrobes
07-30-2012, 04:11 AM
Yeah a stack overflow is exactly that but there is a site (part of google I believe) called stack overflow where you can get help on coding.

http://stackoverflow.com/

But if you do get a table of values printed out then post the table as well and I'll glance my eye at it. There's quite a few good programmers around here as well making maps. someone will be able to see it if not I.

felixrulz
07-31-2012, 09:44 PM
if you do get a table of values printed out then post the table as well and I'll glance my eye at it. There's quite a few good programmers around here as well making maps. someone will be able to see it if not I.

I created some tables and analysed them in MATLAB... and i'm glad I did. it led to quite an embarrassing bug in my code:

* sorry about indenting couldn't get it to work

for (i=0; i<4; i++) { // FOR THE FOUR VERTICES OF QUAD.

// disk coordinates
p_d[n].x = b2[i] * cos(a2[i]) + x0;
p_d[n].y = b2[i] * sin(a2[i]) + y0;

// calcuate texture coordinates
AziCalc (lat_p, long_p, b2[i], a2[i], p_d[n].x, R);
p_d[n].u = 0.5 - longitude/PI2;
p_d[n].v = latitude/PI;

// I HAD THE TEXTURE WRAPPING CODE HERE!

n++;
}

// THE CODE SHOULD HAVE BEEN HERE!
getmax2(n-4, n-1, &max, &index); //get max value for last quad
getmin2(n-4, n-1, &min, &index); //get min value for last quad
while (max - min > 0.5){ // if texture coordinates too far apart
p_d[index].u += 1; //wrap coordinates
getmax2(n-4, n-1, &min, &index); //get max value for last quad
getmin2(n-4, n-1, &min, &index); //get min value for last quad
}
}

The program now works well. Here is the fixed version: Fixed version (http://dl.dropbox.com/u/49415695/Desktop7%5BWorking%5D.7z)
(I compiled this version as 'Release' as this should hopfully get rid of some dll dependencies - not sure though)

Thank you very much Redrobes, you were very helpfull. I also have been reading some of your thread on bitmaps and am enjoying it :)

I have uploaded some of the code here (http://dl.dropbox.com/u/49415695/Source1.c) if anyone wants to look. I

was also wndering if there was a way to add '[SOLVED]' to thread title?