# Quaternions and Lat/Lon

## Recommended Posts

raider8654    100
Hi all,

For the past few days I have been struggling with a problem concerning a graphics program I am writing. It displays a 3d sphere with lat/lon coordinates for the user to see. An overlay of the US is also shown, with true 3d mapping. Everything works up to this point.

Here is what I have:

* The globe can be zoomed in and out by adjusting the horizontal field of view angle
* The globe can also be rotated
* The latitude and longitude are shown for the point that the screen is centered on (this works correctly)

Now, here is what I want to have:

* Show the latitude and longitude for the point where the mouse cursor is centered on (having issues with this)

Now, this is not like in other posts where people ask, how do I unproject a 2d point (the mouse coordinates) to 3d. No. This is different. Instead, I am asking about retrieving the latitude and longitude at the mouse coordinates, and this is feasible since we know the size of the sphere in its 3d space as well as the x,y coordinates and field of view.

So, here is what I have done so far.

I have managed to correctly calculate the change in horizontal angle and change in vertical angle from the center. This is accomplished by way of using hx=tan(fov/2) * d, where d is the distance to the center of the sphere. hx gives me the horizontal length at the current field of view (which necessary since, again, the sphere can be zoomed in and out). And then taking arcsin(((cursor.x-400)*(hx/400))/SPHERE_RADIUS) to get the change in horizontal angle based on horizontal change from the center of the screen (i.e., mouse cursor x movement from center).

So you might think, okay, he's got the latitude and longitude of what the screen is centered on, as well as the change in vertical and horizontal angles. So what's the problem? Can't he just subtract the horizontal angle change from the longitude, and similar for latitude?

As it turns out, it's not that simple. For example, when the user is looking directly at latitude 90 degrees (North pole), then any mouse deviation from the center of the screen necessitates a change in latitude even if only horizontal mouse movement has taken place. But then, when looking at latitude 0 degrees (the equator), then horizontal mouse movement will only affect longitude.

Herein lies the problem: what to do when the user is looking at between those two extremes? Say, latitude 45 degrees? I had managed to use cos(latitude)*(change in vertical) and also add to that the sin(latitude)*(combination of latitude and longitude) to mix the two.. And eventually got it working at lat=0 and lat=90. But never for in-between.

Going back to how I arrived at displaying the centered point latitude and longitude in the first place. I use quaternion for the rotation:

[code] // Get orig quaternion
LatLonToQuaternion(g_angRotateX, g_angRotateY, x, y, z, w);[/code]

That gives me a quaternion from my centered position's latitude and longitude. Yes the lat/lon are flipped in the function's input. But the output quaternion is correct, I have verified.

So then, now that I have a quaternion... I should be able to multiply that by my change in latitude and longitude, yes?

[code] // Rotate by lonX
RotationXYZToQuaternion(0.0f, lonX, 0.0f, tx, ty, tz, tw);
QuaternionMultiply(tx, ty, tz, tw, x, y, z, w, ux, uy, uz, uw);
x=ux;
y=uy;
z=uz;
w=uw;

// Rotate by latX
RotationXYZToQuaternion(latX, 0.0f, 0.0f, tx, ty, tz, tw);
QuaternionMultiply(tx, ty, tz, tw, x, y, z, w, ux, uy, uz, uw);
QuaternionToLatLon(ux, uy, uz, uw, lonF, latF); [/code]

The above code uses a function called RotationXYZToQuaternion to give me a quaternion from a given X,Y,Z rotation. So I first get a quaternion for the change in horizontal angle. Then I multiply that quaternion by the original quaternion that I derived from the lat/lon at the center of the screen. Next, I do the same thing for the change in vertical angle, and finally convert the result back into a latitude and longitude.

In this case the above code works at lat=0, but not at lat=90 or anywhere else in that I do not get back the expected latitude and longitude for where the user is pointing the mouse cursor.

So, this is where I am after about four days of struggling with it -- I have managed to get this far. But am really stuck.

Basically it boils down to this problem:

[b]Given:[/b]

x=longitude at screen center
y=latitude at screen center
qx=change in the sphere's horizontal angle, as calculated from the deviation of the mouse cursor's x position from screen center (*not the same thing as change in longitude, unless lat=0)
qy=change in the sphere's vertical angle, as calculated from the deviation of the mouse cursor's y position from screen center (*also not the same as change in latitude, unless lat=0)

[b]Determine[/b] the new latitude and longitude.

So I ask, if anyone can help me on this, how to solve the problem. What do I do when the user is not centered at latitude=0? How is the change in vertical and horizontal angles applied when latitude=30 degrees, for example?

Any help would be [i]immensely[/i] appreciated!

Thanks.

##### Share on other sites
jyk    2094
I admit that I'd have to read that again (and maybe again after that ;) in order to understand what it is you're doing, but meanwhile, let me ask this. If the objective is to determine the latitude and longitude corresponding to the cursor position, why not simply raycast against the sphere and compute the latitude and longitude from the intersection point?

##### Share on other sites
haegarr    7372
[quote name='Jack Smith' timestamp='1305949219' post='4813747']
...
Now, this is not like in other posts where people ask, how do I unproject a 2d point (the mouse coordinates) to 3d. No. This is different. Instead, I am asking about retrieving the latitude and longitude at the mouse coordinates, and this is feasible since we know the size of the sphere in its 3d space as well as the x,y coordinates and field of view.
...
[/quote]
Well, it is AFAIS not really different; it just needs 2 more steps. So I second jyk's implicit suggestion:

1. Compute a ray in global space that describes all locations that are projected onto the pixel covered by the mouse pointer's hot spot.
2. Transform the ray using the inverse of the local-to-global transformation of the sphere.
3. Compute the nearest intersection point of the ray with the sphere. This gives you a Cartesian co-ordinate.
4. Convert the Cartesian co-odinate into a spherical one.
5. Adapt the azimuthal angle to yield in the latitude.

The above should be sufficient as long as a sphere is used. When using an ellipsoid, a non-uniform scaling incorporated into the object's transformation should do the trick, although the term longitude / latitude becomes versatile; you have to decide which kind of longitude / latitude to use then. (But if you want to use the spherical harmonics series to approximate the real shape of the earth ...)

##### Share on other sites
raider8654    100
haegarr, jyk:

Thank you for your replies. From what I understand, ray intersection functions provided by Direct3D (the graphics library I am using) work with polygons. So for example, my sphere would need to be made out of (say, triangles) to compute the intersection from the ray to the sphere.

Is this correct?

See, my "sphere" is actually just the shape to which my 3d map of line segments conforms. In other words, I do not have any triangles that would intersect against the ray. Even if I did, I think I would lose accuracy because only if my sphere were made of trillions of triangles would I get any accuracy as to which lat/lon the mouse cursor has hit on.

So, if I do go the route of ray intersection.. It would be nice, I admit, since maybe that will get me what I want. But I need some help on implementing it.

What I gather at this point is that I can get a 3d ray, and would need to come up with the Cartesian coordinate at which it intersects with the sphere. So if I know the radius of the sphere in its 3d space, and am given a 3d ray, how would I get these coordinates of intersection?

Thanks again.

##### Share on other sites
jyk    2094
[quote name='Jack Smith' timestamp='1305994440' post='4813879']
From what I understand, ray intersection functions provided by Direct3D (the graphics library I am using) work with polygons. So for example, my sphere would need to be made out of (say, triangles) to compute the intersection from the ray to the sphere.

Is this correct?[/quote]
Nope; you can compute the intersection with the sphere directly (you don't have to use a mesh representation).

[quote]So if I know the radius of the sphere in its 3d space, and am given a 3d ray, how would I get these coordinates of intersection?[/quote]
Google/search for (e.g.) 'ray sphere intersection', 'line sphere intersection', 'sphere raytrace', or 'sphere raycast'. (The algorithm is straightforward and is well documented online.)

##### Share on other sites
haegarr    7372
Hidden
[quote name='Jack Smith' timestamp='1305994440' post='4813879']

Thank you for your replies. From what I understand, ray intersection functions provided by Direct3D (the graphics library I am using) work with polygons. So for example, my sphere would need to be made out of (say, triangles) to compute the intersection from the ray to the sphere.

Is this correct?
[/quote]

No, it need not, and IMHO it even shall not.

Instead, using the implicit formula for a sphere with radius R, centered around the origin, i.e.
x[sup]2[/sup] + y[sup]2[/sup] + z[sup]2[/sup] == R[sup]2[/sup]
so that any point [ x y z ][sup]t[/sup] in space belongs to sphere's surface if it fulfills the above condition.

A ray, in general, can be defined by any point [b]r[/b][sub]0[/sub] on the ray (lets name it its origin) and the direction [b]r[/b][sub]d[/sub] of the ray (let's call it its track), the latter one scaled to an arbitrary length by a factor k, so that
[b]r[/b]( k ) := [b]r[/b][sub]0[/sub] + k * [b]r[/b][sub]d[/sub]

For any k set into the above equation, we get a point on the ray.

Now let us define the transformation [b]M[/b] that maps the sphere from is model (or local) space into the global space. This transformation is a matrix that is typically constructed using a rotation and a translation, as usual. D3D supports this kind of transformation matrices, as you already know. You use this transformation to place the model of the sphere in the world.

There is further a view transformation matrix [b]V[/b], i.e. the inverse of the camera transformation. This transformation is also supported by D3D already.

The current mouse co-ordinate [ m[sub]h[/sub] m[sub]v[/sub] ][sup]t[/sup] in screen space has 1st to be transformed into view space. This gives
[b]m[/b] = [ m[sub]h[/sub]/p[sub]h[/sub]*v[sub]w[/sub] m[sub]v[/sub]/p[sub]v[/sub]*v[sub]h[/sub] v[sub]d[/sub] ][sup]t[/sup]
where p[sub]h[/sub] and p[sub]v[/sub] are the screen's/window's pixel amount in horizontal and vertical dimension, resp., v[sub]w[/sub] and v[sub]h[/sub] are the world sizes of the view plane in horizontal and vertical dimensions, resp., and v[sub]d[/sub] ist the distance of the view plane from the eye point.

With this point the ray given in view space can be computed as either
a) for orthogonal projection
[b]r[/b]( k ) := [b]m[/b] - v[sub]d [/sub]* [ 0 0 1 ][sup]t[/sup] + k * [ 0 0 1 ][sup]t[/sup]
assuming that the view plane lies in positive z direction, or else
b) for perspective projection
[b]r[/b]( k ) := k * [b]r[/b][sub]d[/sub] with [b]r[/b][sub]d[/sub] := [b]m[/b] / | [b]m[/b] |

This is in fact an un-project step. However, now there is a ray, given by its origin (a point vector) and its track (a direction vector), both transformable by a matrix. We transform the ray now into model local space of the sphere. We do so by 1st transforming into the global space and afterwards into the model space. Using row vector notation (as you've spoken of D3D), this looks like
[b]r[/b]'( k ) = [b]r[/b]( k ) * [b]V[/b][sup]-1[/sup] * [b]M[/b][sup]-1[/sup]

Now the ray is in the same co-ordinate system as the implicit form of the sphere is. Hence we can set in the one into the other, i.e. we compute all the points that are both onto the surface of the sphere as well as onto the ray (in other word, the intersection points). That looks like

x[sup]2[/sup] + y[sup]2[/sup] + z[sup]2[/sup] - R[sup]2[/sup] == 0
with
x := r'[sub]0x[/sub] - k * r'[sub]dx[/sub]
[sup][size=3]
y := r'[sub]0y[/sub] - k * r'[sub]dy[/sub]
[sub][size=3]
z := r'[sub]0z[/sub] - k * r'[sub]dz[/sub]
what give you a quadratic equation with k as variable. Solve this by using the famous p,q-formula, and you'll get 0, 1, or 2 solution for k. 0 solutions means that the ray passes the sphere because the mouse is outside the projection. 1 solution means that the ray tangents the sphere. 2 solution means that the ray enters the sphere in front and exits it in the back; you're then interested in the lower k only, of course.

Now set in k (if any is found) into [b]r[/b]'( k ) to compute the cartesian co-ordinates of the intersection of interest. Use this to compute the spherical co-ordinates as usual. However, spherical co-ordinates are usually defined with an angle in [0,pi] reaching from pole to pole, but you're interested in the latitude within [-pi,pi], so make a simple adjustment.

Well, double check my writings as always, because it is much stuff written here from the top of my head, and it happens easily to make a mistake.[/size][/sub][/size][/sup]

raider8654    100
[quote name='jyk' timestamp='1305998501' post='4813903']
[quote name='Jack Smith' timestamp='1305994440' post='4813879']
From what I understand, ray intersection functions provided by Direct3D (the graphics library I am using) work with polygons. So for example, my sphere would need to be made out of (say, triangles) to compute the intersection from the ray to the sphere.

Is this correct?[/quote]
Nope; you can compute the intersection with the sphere directly (you don't have to use a mesh representation).

[quote]So if I know the radius of the sphere in its 3d space, and am given a 3d ray, how would I get these coordinates of intersection?[/quote]
Google/search for (e.g.) 'ray sphere intersection', 'line sphere intersection', 'sphere raytrace', or 'sphere raycast'. (The algorithm is straightforward and is well documented online.)
[/quote]

Thank you so much! I have now come up with a function that gives me the Latitude and Longitude of where the user is pointing the cursor on the sphere. In case anyone is interested:

[code]
BOOL FindIntersection(POINT &curPos, D3DVIEWPORT9 &vpMap, D3DXMATRIX &matProjection, D3DXMATRIX &matView,
D3DXMATRIX &matWorld, FLOAT &lat, FLOAT &lon)
{
D3DXMATRIX m;
D3DXVECTOR3 rayOrigin;
D3DXVECTOR3 rayDir;
D3DXVECTOR3 p1;
D3DXVECTOR3 p2;
D3DXVECTOR3 v1;
D3DXVECTOR3 v2;
D3DXVECTOR3 v;
FLOAT determinant;
FLOAT theta;
FLOAT phi;
FLOAT rho;
FLOAT S;
FLOAT a;
FLOAT b;
FLOAT c;
FLOAT t1;
FLOAT t2;
bool bDoesIntersect;

vpMap.Width=792;
vpMap.Height=546;

D3DXVECTOR3 inP1(curPos.x, curPos.y, 0.0f);
D3DXVec3Unproject(&rayOrigin, &inP1, &vpMap, &matProjection, &matView, &matWorld);
D3DXVECTOR3 inP2(curPos.x, curPos.y, 1.0f);
D3DXVec3Unproject(&v2, &inP2, &vpMap, &matProjection, &matView, &matWorld);

D3DXVec3Normalize(&rayDir, D3DXVec3Subtract(&v, &rayOrigin, &v2));

// p = td + p0
// a = d*d
// b = 2d*(p0-pc)
// c = (p0-pc)*(p0-pc)-r^2

a = D3DXVec3Dot(&rayDir, &rayDir);
b = D3DXVec3Dot(D3DXVec3Scale(&v, &rayDir, 2), &rayOrigin);
c = D3DXVec3Dot(&rayOrigin, &rayOrigin) - pow(SPHERE_SIZE, 2);

// Calculate determinant
determinant = pow(b, 2) - 4*a*c;

if (determinant >= 0)
{
// There is at least one point of intersection
t1 = (-b + sqrt(determinant)) / (2*a);
t2 = (-b - sqrt(determinant)) / (2*a);

// Plug into p = td + p0 and solve for p

a = sqrt(pow(p1.x-rayOrigin.x,2)+pow(p1.y-rayOrigin.y,2)+pow(p1.z-rayOrigin.z,2));
b = sqrt(pow(p2.x-rayOrigin.x,2)+pow(p2.y-rayOrigin.y,2)+pow(p2.z-rayOrigin.z,2));

if (a > b)
{
// p1 is closest to viewer
v.x = p1.x;
v.y = p1.y;
v.z = p1.z;
}
else
{
// p2 is closest to viewer
v.x = p2.x;
v.y = p2.y;
v.z = p2.z;
}

// Now convert from Cartesian to spherical coordinates
rho = sqrt(pow(v.x,2)+pow(v.z,2)+pow(v.y,2));
S = sqrt(pow(v.x,2)+pow(v.z,2));
phi = acos(v.y/rho);
if (v.x >= 0)
{
theta = asin(v.z/S);
}
else
{
theta = PI - asin(v.z/S);
}

if (theta >= PI)
{
theta -= 2*PI;
}

bDoesIntersect = true;
lat = phi;
lon = theta;
}
else
{
// Cursor does not intersect with sphere
bDoesIntersect = false;
}

return bDoesIntersect;
}
[/code]

Note that I have the vector 'v' has the y and z flipped. That is just how it works in my 3d world.

Once I get the lat/lon, which are returned as radians, I convert them to degrees and add the appropriate offset (in my case I subtract 90 degrees from the latitude, and also if longitude exceeds 180 I subtract 360 from it).

In any case the above code is specific to my software, but could be adapted to someone else's application if they needed it.

So now, regardless of zoom level and rotation, the user can see the latitude and longitude that they have the cursor pointed at. "no intersection" is displayed whenever the cursor hovers outside the sphere. So it works perfectly.

Thanks again!!!