blog/A-IT/00-CS基础/50-algorithm/计算几何

// 3D vector  
class Vector3  
{  
public:  
    Vector3(float fx, float fy, float fz)  
        :x(fx), y(fy), z(fz)  
    {  
    }  
  
    // Subtract  
    Vector3 operator - (const Vector3& v) const  
    {  
        return Vector3(x - v.x, y - v.y, z - v.z) ;  
    }  
  
    // Dot product  
    float Dot(const Vector3& v) const  
    {  
        return x * v.x + y * v.y + z * v.z ;  
    }  
  
    // Cross product  
    Vector3 Cross(const Vector3& v) const  
    {  
        return Vector3(  
            y * v.z - z * v.y,  
            z * v.x - x * v.z,  
            x * v.y - y * v.x ) ;  
    }  
  
public:  
    float x, y, z ;  
};  
  
// Determine whether two vectors v1 and v2 point to the same direction
// v1 = Cross(AB, AC)
// v2 = Cross(AB, AP)  
bool SameSide(Vector3 A, Vector3 B, Vector3 C, Vector3 P)  
{  
    Vector3 AB = B - A ;  
    Vector3 AC = C - A ;  
    Vector3 AP = P - A ;  
  
    Vector3 v1 = AB.Cross(AC) ;  
    Vector3 v2 = AB.Cross(AP) ;  
  
    // v1 and v2 should point to the same direction  
    return v1.Dot(v2) >= 0 ;  
}  
  
// Same side method// Determine whether point P in triangle ABC  
bool PointinTriangle1(Vector3 A, Vector3 B, Vector3 C, Vector3 P)  
{  
    return SameSide(A, B, C, P) &&  
        SameSide(B, C, A, P) &&  
        SameSide(C, A, B, P) ;  
}

// Determine whether point P in triangle ABC  
bool PointinTriangle(Vector3 A, Vector3 B, Vector3 C, Vector3 P)  
{  
    Vector3 v0 = C - A ;  
    Vector3 v1 = B - A ;  
    Vector3 v2 = P - A ;  
  
    float dot00 = v0.Dot(v0) ;  
    float dot01 = v0.Dot(v1) ;  
    float dot02 = v0.Dot(v2) ;  
    float dot11 = v1.Dot(v1) ;  
    float dot12 = v1.Dot(v2) ;  
  
    float inverDeno = 1 / (dot00 * dot11 - dot01 * dot01) ;  
  
    float u = (dot11 * dot02 - dot01 * dot12) * inverDeno ;  
    if (u < 0 || u > 1) // if u out of range, return directly  
    {  
        return false ;  
    }  
  
    float v = (dot00 * dot12 - dot01 * dot02) * inverDeno ;  
    if (v < 0 || v > 1) // if v out of range, return directly  
    {  
        return false ;  
    }  
  
    return u + v <= 1 ;  
}

#include <glm/glm.hpp>
#include <glm/gtc/matrix_transform.hpp>
#include <glm/gtc/type_ptr.hpp>
#include <glm/gtx/intersect.hpp>
#include <glm/gtx/transform.hpp>
#include <iostream>

// 检查射线与球体的交点
bool intersectRaySphere(
    const glm::vec3& rayOrigin, 
    const glm::vec3& rayDirection, 
    const glm::vec3& sphereCenter, 
    float sphereRadius, 
    glm::vec3& intersectionPoint) 
{
    float t0, t1;  // 解

    glm::vec3 L = sphereCenter - rayOrigin;
    float tca = glm::dot(L, rayDirection);
    if (tca < 0) return false;
    float d2 = glm::dot(L, L) - tca * tca;
    float radius2 = sphereRadius * sphereRadius;
    if (d2 > radius2) return false;
    float thc = sqrt(radius2 - d2);
    t0 = tca - thc;
    t1 = tca + thc;

    if (t0 > t1) std::swap(t0, t1);

    if (t0 < 0) {
        t0 = t1;  // 如果 t0 是负数，则使用 t1
        if (t0 < 0) return false;  // 如果 t1 也小于零，则射线在球体之后
    }

    intersectionPoint = rayOrigin + t0 * rayDirection;
    return true;
}

int main() {
    // 示例参数
    glm::vec3 rayOrigin(0.0f, 0.0f, 5.0f);
    glm::vec3 rayDirection(0.0f, 0.0f, -1.0f);
    glm::vec3 sphereCenter(0.0f, 0.0f, 0.0f);
    float sphereRadius = 1.0f;

    // 检查射线与球体的交点
    glm::vec3 intersectionPoint;
    if (intersectRaySphere(rayOrigin, rayDirection, sphereCenter, sphereRadius, intersectionPoint)) {
        std::cout << "Intersection Point: " << glm::to_string(intersectionPoint) << std::endl;
    } else {
        std::cout << "No intersection" << std::endl;
    }

    return 0;
}