Rotas na superfície de uma esfera - Encontre geodésico?

Estou trabalhando com alguns amigos em um jogo baseado em navegador, onde as pessoas podem se mover em um mapa 2D. Já faz quase 7 anos e as pessoas ainda jogam esse jogo, então estamos pensando em uma maneira de lhes dar algo novo. Desde então, o mapa do jogo era um plano limitado e as pessoas podiam passar de (0, 0) para (MAX_X, MAX_Y) em incrementos quantificados de X e Y (imagine-o como um grande tabuleiro de xadrez).
Acreditamos que é hora de dar uma outra dimensão, então, apenas algumas semanas atrás, começamos a nos perguntar como o jogo poderia parecer com outros mapeamentos:

• Avião ilimitado com movimento contínuo: isso pode ser um passo à frente, mas ainda não estou convencido.
• Mundo Toroidal (movimento contínuo ou quantizado): sinceramente já trabalhei com o toro, mas desta vez quero algo mais ...
• Mundo esférico com movimento contínuo: isso seria ótimo!

O que queremos Os navegadores de usuários recebem uma lista de coordenadas como (latitude, longitude) para cada objeto no mapa de superfície esférico; os navegadores devem mostrar isso na tela do usuário, renderizando-os dentro de um elemento da web (tela talvez? isso não é um problema). Quando as pessoas clicam no plano, convertemos o (mouseX, mouseY) em (lat, lng) e o enviamos ao servidor que precisa calcular uma rota entre a posição atual do usuário e o ponto clicado.

O que temos Começamos a escrever uma biblioteca Java com muitas matemáticas úteis para trabalhar com matrizes de rotação, quaterniões, ângulos de Euler, traduções etc. Reunimos tudo isso e criamos um programa que gera pontos de esfera, os renderiza e os mostra ao usuário dentro de um JPanel. Conseguimos capturar cliques e convertê-los em cordas esféricas e fornecer alguns outros recursos úteis, como rotação de exibição, escala, tradução etc. O que temos agora é como um pequeno (muito pouco) mecanismo que simula a interação do cliente e do servidor. O lado do cliente mostra pontos na tela e captura outras interações, o lado do servidor renderiza a visualização e faz outros cálculos, como interpolar a rota entre a posição atual e o ponto clicado.

Onde está o problema? Obviamente, queremos ter o caminho mais curto para interpolar entre os dois pontos da rota . Usamos quaternions para interpolar entre dois pontos na superfície da esfera e isso pareceu funcionar bem até que notei que não estávamos obtendo o caminho mais curto na superfície da esfera:

Achamos que o problema era que a rota é calculada como a soma de duas rotações em torno dos eixos X e Y. Então, mudamos a maneira como calculamos o quaternion de destino: obtemos o terceiro ângulo (o primeiro é latitude, o segundo é longitude, o terceiro é a rotação sobre o vetor que aponta para a nossa posição atual) que chamamos de orientação. Agora que temos o ângulo "orientação", giramos o eixo Z e, em seguida, usamos o vetor de resultado como eixo de rotação para o quaternion de destino (você pode ver o eixo de rotação em cinza):

O que obtivemos é a rota correta (você pode vê-la em um grande círculo), mas chegaremos a isso SOMENTE se o ponto de rota inicial estiver em latitude, longitude (0, 0), o que significa que o vetor inicial é (sphereRadius, 0 0). Com a versão anterior (imagem 1), não obtemos um bom resultado, mesmo quando o ponto inicial é 0, 0, então acho que estamos caminhando em direção a uma solução, mas o procedimento a seguir para obter essa rota é um pouco "estranho" " talvez?

Na imagem a seguir, você tem uma visão do problema que ocorre quando o ponto inicial não é (0, 0), como você pode ver, o ponto inicial não é o vetor (sphereRadius, 0, 0) e como você pode ver o ponto de destino (que está corretamente desenhado!) não está na rota.

O ponto magenta (aquele que fica na rota) é o ponto final da rota girado sobre o centro da esfera de (-startLatitude, 0, -startLongitude). Isso significa que, se eu calcular uma matriz de rotação e aplicá-la a todos os pontos da rota, talvez eu obtenha a rota real, mas começo a pensar que há uma maneira melhor de fazer isso.

Talvez eu deva tentar fazer o avião atravessar o centro da esfera e os pontos da rota, interceptá-lo com a esfera e obter a geodésica? Mas como?

Desculpe por ser muito detalhado e talvez por um inglês incorreto, mas essa coisa está me surpreendendo!

EDIT: Código abaixo funciona muito bem! Obrigado a todos:

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;
}

Seu problema é puramente bidimensional, no plano formado pelo centro da esfera e seus pontos de origem e destino. Usar quaternions é realmente tornar as coisas mais complexas, porque além de uma posição em uma esfera 3D, um quaternion codifica uma orientação.

Você já pode ter algo para interpolar em um círculo, mas por precaução, aqui está um código que deve funcionar.

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;
}

Certifique-se de que ambos os quaternions estejam no mesmo hemisfério na hiperesfera. Se o produto escalar for menor que 0, eles não serão. Nesse caso, negue um deles (negue cada um de seus números), para que eles estejam no mesmo hemisfério e fornecerão o caminho mais curto. Pseudo-código:

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

Minha resposta aqui explica em detalhes o que nega cada termo do quaternion e por que ainda é a mesma orientação, apenas do outro lado da hiperesfera.

EDIT a função de interpolação deve ficar assim:

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;
}

Como você deseja ter o V3Dseu interpolador, a abordagem mais simples é pular os quaternions completamente. Converta os pontos inicial e final em V3De faça slerp entre eles.

Se você insiste em usar quaternions, o quaternion que representa a rotação de Ppara Qtem direção P x Qe wde P . Q.