Routes on a sphere surface - Find geodesic?

Question

I'm working with some friends on a browser based game where people can move on a 2D map. It's been almost 7 years and still people play this game so we are thinking of a way to give them something new. Since then the game map was a limited plane and people could move from (0, 0) to (MAX_X, MAX_Y) in quantized X and Y increments (just imagine it as a big chessboard).
We believe it's time to give it another dimension so, just a couple of weeks ago, we began to wonder how the game could look with other mappings:

• Unlimited plane with continous movement: this could be a step forward but still i'm not convinced.
• Toroidal World (continous or quantized movement): sincerely I worked with torus before but this time I want something more...
• Spherical world with continous movement: this would be great!

What we want Users browsers are given a list of coordinates like (latitude, longitude) for each object on the spherical surface map; browsers must then show this in user's screen rendering them inside a web element (canvas maybe? this is not a problem). When people click on the plane we convert the (mouseX, mouseY) to (lat, lng) and send it to the server which has to compute a route between current user's position to the clicked point.

What we have We began writing a Java library with many useful maths to work with Rotation Matrices, Quaternions, Euler Angles, Translations, etc. We put it all together and created a program that generates sphere points, renders them and show them to the user inside a JPanel. We managed to catch clicks and translate them to spherical coords and to provide some other useful features like view rotation, scale, translation etc. What we have now is like a little (very little indeed) engine that simulates client and server interaction. Client side shows points on the screen and catches other interactions, server side renders the view and does other calculus like interpolating the route between current position and clicked point.

Where is the problem? Obviously we want to have the shortest path to interpolate between the two route points. We use quaternions to interpolate between two points on the surface of the sphere and this seemed to work fine until i noticed that we weren't getting the shortest path on the sphere surface:

We though the problem was that the route is calculated as the sum of two rotations about X and Y axis. So we changed the way we calculate the destination quaternion: We get the third angle (the first is latitude, the second is longitude, the third is the rotation about the vector which points toward our current position) which we called orientation. Now that we have the "orientation" angle we rotate Z axis and then use the result vector as the rotation axis for the destination quaternion (you can see the rotation axis in grey):

What we got is the correct route (you can see it lays on a great circle), but we get to this ONLY if the starting route point is at latitude, longitude (0, 0) which means the starting vector is (sphereRadius, 0, 0). With the previous version (image 1) we don't get a good result even when startin point is 0, 0, so i think we're moving towards a solution, but the procedure we follow to get this route is a little "strange" maybe?

In the following image you get a view of the problem we get when starting point is not (0, 0), as you can see starting point is not the (sphereRadius, 0, 0) vector, and as you can see the destination point (which is correctly drawn!) is not on the route.

The magenta point (the one which lays on the route) is the route's ending point rotated about the center of the sphere of (-startLatitude, 0, -startLongitude). This means that if i calculate a rotation matrix and apply it to every point on the route maybe i'll get the real route, but I start to think that there's a better way to do this.

Maybe I should try to get the plane through the center of the sphere and the route points, intersect it with the sphere and get the geodesic? But how?

Sorry for being way too verbose and maybe for incorrect English but this thing is blowing my mind!

EDIT: Code below works just fine! Thanks to everybody:

public void setRouteStart(double srcLat, double srcLng, double destLat, destLng) {
    //all angles are in radians
    u = Choords.sphericalToNormalized3D(srcLat, srcLng);
    v = Choords.sphericalToNormalized3D(destLat, destLng);
    double cos = u.dotProduct(v);
    angle = Math.acos(cos);
    if (Math.abs(cos) >= 0.999999) {
        u = new V3D(Math.cos(srcLat), -Math.sin(srcLng), 0);
    } else {
        v.subtract(u.scale(cos));
        v.normalize();
    }
}

public static V3D sphericalToNormalized3D( double radLat, double radLng) {
    //angles in radians
    V3D p = new V3D();
    double cosLat = Math.cos(radLat);
    p.x = cosLat*Math.cos(radLng);
    p.y = cosLat*Math.sin(radLng);
    p.z = Math.sin(radLat);
    return p;
}

public void setRouteDest(double lat, double lng) {
    EulerAngles tmp = new AngoliEulero(
       Math.toRadians(lat), 0, -Math.toRadians(lng));
    qtEnd.setInertialToObject(tmp);
    //do other stuff like drawing dest point...
}

public V3D interpolate(double totalTime, double t) {
    double _t = angle * t/totalTime;
    double cosA = Math.cos(_t);
    double sinA = Math.sin(_t);
    V3D pR = u.scale(cosA);
    pR.sum(
       v.scale(sinA)
    );
    return pR;
}

Answer 1

Your problem is purely two-dimensional, in the plane formed by the sphere centre and your source and destination points. Using quaternions is actually making things more complex, because in addition to a position on a 3D sphere, a quaternion encodes an orientation.

You may already have something to interpolate on a circle, but just in case, here is some code that should work.

V3D u, v;
double angle;

public V3D geographicTo3D(double lat, double long)
{
    return V3D(sin(Math.toRadians(long)) * cos(Math.toRadians(lat)),
               cos(Math.toRadians(long)) * cos(Math.toRadians(lat)),
               sin(Math.toRadians(lat)));
}

public V3D setSourceAndDest(double srcLat, double srcLong,
                            double dstLat, double dstLong)
{
    u = geographicTo3D(srcLat, srcLong);
    V3D tmp = geographicTo3D(dstLat, dstLong);
    angle = acos(dot(u, tmp));
    /* If there are an infinite number of routes, choose
     * one arbitrarily. */
    if (abs(dot(u, tmp)) >= 0.999999)
        v = V3D(cos(srcLong), -sin(srcLong), 0);
    else
        v = normalize(tmp - dot(u, tmp) * u);
}

public V3D interpolate(double totalTime, double t)
{
    double a = t / totalTime * angle;
    return cos(a) * u + sin(a) * v;
}

Answer 2

Make sure both quaternions are on the same hemisphere on the hypersphere. If their dot-product is less than 0, then they are not. In that case negate one of them (negate each of the its numbers), so they are on the same hemisphere and will give you the shortest path. Pseudocode:

quaternion from, to;

// is "from" and "to" on the same hemisphere=
if(dot(from, to) < 0.0)
{
    // put "from" to the other hemisphere, so its on the same as "to"
    from.x = -from.x;
    from.y = -from.y;
    from.z = -from.z;
    from.w = -from.w;
}

// now simply slerp them

My answer here explains in detail what negating each term of the quaternion does and why it is still the same orientation, just on the other side of the hypersphere.

---

EDIT the interpolation function should look like this:

public V3D interpolate(double totalTime, double t) {
    double _t = t/totalTime;    
    Quaternion tmp;
    if(dot(qtStart, qtEnd) < 0.0)
    {
        tmp.x = -qtEnd.x;
        tmp.y = -qtEnd.y;
        tmp.z = -qtEnd.z;
        tmp.w = -qtEnd.w;
    }
    else
    {
        tmp = qtEnd;
    }
    Quaternion q = Quaternion.Slerp(qtStart, tmp, _t);
    RotationMatrix.inertialQuatToIObject(q);
    V3D p = matInt.inertialToObject(V3D.Xaxis.scale(sphereRadius));
    //other stuff, like drawing point ...
    return p;
}

Answer 3

Since you want a V3D back from your interpolator, the simplest approach is to skip the quaternions entirely. Convert the start and end points to V3D and slerp between them.

If you insist on using quaternions then the quaternion representing the rotation from P to Q has direction P x Q and w of P . Q.