射线与球面的交集

Question

我最近一直在尝试让自己成为一个光线追踪器来计算给定位置的经度和纬度值。我的球体位于原点 {0, 0, 0} 并且我的相机始终指向原点。

#include <array>
#include <vector>
#include <iostream>

#include "glm/glm.hpp"
#include "glm/gtc/matrix_transform.hpp"

std::array<double, 3> ecef_to_geo(std::array<double, 3> ecef) {
    std::array<double, 3> geo{ 0, 0, 0 };   //Results go here (Lat, Lon, Altitude)
    double  a = 6378137.0;              //WGS-84 semi-major axis
    double e2 = 6.6943799901377997e-3;  //WGS-84 first eccentricity squared
    double a1 = 4.2697672707157535e+4;  //a1 = a*e2
    double a2 = 1.8230912546075455e+9;  //a2 = a1*a1
    double a3 = 1.4291722289812413e+2;  //a3 = a1*e2/2
    double a4 = 4.5577281365188637e+9;  //a4 = 2.5*a2
    double a5 = 4.2840589930055659e+4;  //a5 = a1+a3
    double a6 = 9.9330562000986220e-1;  //a6 = 1-e2
    double zp, w2, w, r2, r, s2, c2, s, c, ss;
    double g, rg, rf, u, v, m, f, p, x, y, z;
    double n, lat, lon, alt;

    x = ecef[0];
    y = ecef[1];
    z = ecef[2];
    zp = std::abs(z);
    w2 = x * x + y * y;
    w = std::sqrt(w2);
    r2 = w2 + z * z;
    r = std::sqrt(r2);
    geo[1] = std::atan2(y, x);       //Lon (final)
    s2 = z * z / r2;
    c2 = w2 / r2;
    u = a2 / r;
    v = a3 - a4 / r;
    if (c2 > 0.3) {
        s = (zp / r) * (1.0 + c2 * (a1 + u + s2 * v) / r);
        geo[0] = std::asin(s);      //Lat
        ss = s * s;
        c = std::sqrt(1.0 - ss);
    }
    else {
        c = (w / r) * (1.0 - s2 * (a5 - u - c2 * v) / r);
        geo[0] = std::acos(c);      //Lat
        ss = 1.0 - c * c;
        s = std::sqrt(ss);
    }
    g = 1.0 - e2 * ss;
    rg = a / std::sqrt(g);
    rf = a6 * rg;
    u = w - rg * c;
    v = zp - rf * s;
    f = c * u + s * v;
    m = c * v - s * u;
    p = m / (rf / g + f);
    geo[0] = geo[0] + p;      //Lat
    geo[2] = f + m * p / 2.0;     //Altitude
    if (z < 0.0) {
        geo[0] *= -1.0;     //Lat
    }
    geo[0] = geo[0] * 180 / glm::pi<double>();
    geo[1] = geo[1] * 180 / glm::pi<double>();
    return(geo);    //Return Lat, Lon, Altitude in that order
}

bool solve_quadratic(float a, float b, float c, float& t0, float& t1) {

    float discriminant = (b * b) - (4 * a * c);
    if (discriminant < 0) //no solution
        return false;

    if (discriminant == 0) { //1 solution
        t0 = -b / 2.0f * a;
    }
    else if (discriminant > 0) {//2 solutions;
        auto droot = std::sqrt(discriminant);
        t0 = (-b - droot) / (2.0f * a);
        t1 = (-b + droot) / (2.0f * a);
    }

    return true;
}

void sphere_intersection(glm::vec3 ray_origin, glm::vec3 ray_direction) {
    double r = 6378137.0;
    float a = glm::dot(ray_direction, ray_direction);
    float b = glm::dot(ray_origin, ray_direction) * 2.0f;
    float c = glm::dot(ray_origin, ray_origin) - (r * r);

    float t0{ 0 }, t1{ 0 };
    if (!solve_quadratic(a, b, c, t0, t1))
        return;

    glm::vec3 hit1 = ray_origin * ray_direction * t0;
    glm::vec3 hit2 = ray_origin * ray_direction * t1;
    auto r1 = ecef_to_geo({ hit1.x, hit1.y, hit1.z });
    auto r2 = ecef_to_geo({ hit2.x, hit2.y, hit2.z });

    std::cout << "------results for ray_direction(" << ray_direction.x << ", " << ray_direction.y << ", " << ray_direction.z << ") ---------- " << std::endl;
    std::cout << "t0 -> " << t0 << " | hit1 -> (" << r1[0] << ", " << r1[1] << ")" << std::endl;
    std::cout << "t1 -> " << t1 << " | hit2 -> (" << r2[0] << ", " << r2[1] << ")" << std::endl;
}


int main() {
//-x=anti meridian
//+x=meridian
//+y=india
//-y=panama
//-z=south pole
//+z=north pole

    std::vector<glm::vec3> camera_positions = { 
    glm::vec3{ 1000000, 0, 0 },
    glm::vec3{ 0, 1000000, 0 },
    glm::vec3{ 0, 0, 1000000 },
    glm::vec3{ -1000000, 0, 0 },
    glm::vec3{ -1000000, 1000000, 0 },
    glm::vec3{ -1000000, -1000000, 0 },
    glm::vec3{ -1000000, 0, -1000000 },
    glm::vec3{ -1000000, -1000000, -1000000 }
    };


    for (size_t i = 0; i < camera_positions.size(); i++) {
        sphere_intersection(camera_positions[i], camera_positions[i]);
        auto result1_e = ecef_to_geo(std::array<double, 3>{camera_positions[i].x, camera_positions[i].y, camera_positions[i].z});
        std::cout << "ecef_to_geo -> (" << result1_e[0] << ", " << result1_e[1] << ")" << std::endl << std::endl;
    }

}

正如我们在输出窗口中看到的；对于所有零/正

x

值，

t1

看起来像是正确的解决方案，对于负

x

值，

t0

看起来像是正确的解决方案（有时带有错误的符号）。我们可以确认，因为

ecef_to_geo

函数总是产生正确的纬度/经度值（球体/椭球体差异会产生轻微差异）。 我一直期望正

t

值是解决方案。我的错误可能在于射线方向，因为它与射线原点相同，但我不能肯定地说。