doorxp-blog

doorxp.com

射线与球的相交性检测

从图形来说

射线和圆相交, origin是射线起点, dir是射线的方向向量。p0,p1是两个交点,center为圆心,半径为R,d为圆心到射线的距离。

我们先以2D切面图来说明,当射线和圆相交的时候,可以看到,球心 center 到射线 ray 的距离 d <= R,这个即为相交的条件。那么射线与球相切就转化为了球心到射线的距离d的判断。先求出d:

  1. 设圆心在射线上的投影为c',则 origin,center, c' 形成了一个直角三角形。

  2. 获得射线起点到圆心的向量 Voc = Vcenter - Vorigin

  3. 在射线方向上的投影为: Poc = Voc·dir

  4. 勾股定理:d·d = Voc·Voc - Poc·Poc

可以求出d的数值,

  • d < R,射线穿过圆,与圆有两个交点。

  • d = R,射线与圆相切,有一个交点为切点。

  • d > R,射线在圆外,没有交点。

接下来求P0,P1:

  1. c',center,P0 or P1点构成直角三角形。

  2. P0 or P1到c'的距离 tca·tca = R·R - d·d;

  3. 有如下式子

  4. P0 = dir·( |Poc| - tca );
    P1 = dir·( |Poc| + tca );
要注意,没有交点的时候, tca·tca < 0 是没办法开平方的

推导三维情况可以照上面的去做,dot能保证投影点在同一个平面上的。

附代码

bool Intersect(const Ray& ray, const Sphere& sphere, float& t0, float& t1){
    Vector3 oc = sphere.GetCenter() - ray.GetOrigin();    float projoc = dot(ray.GetDirection(), oc);    if (projoc < 0)        return false;    float oc2 = dot(oc, oc);    float distance2 = oc2 - projoc * projoc;   //计算出的球心到射线的距离    if (distance2 > sphere.GetRadiusSquare())        return false;    float discriminant = sphere.GetRadiusSquare() - distance2;  //使用勾股定理,计算出另一条边的长度    if(discriminant < FLOAT_EPSILON)   //表明只有一个交点,射线与球相切
        t0 = t1 = projoc;    else
    {
        discriminant = sqrt(discriminant);
        t0 = projoc - discriminant;
        t1 = projoc + discriminant;        if (t0 < 0)
            t0 = t1;
    }    return true;
}

从方程角度来看

射线方程:ray : P(t) = O + D·t ( t >= 0 )
球的方程:sphere : sqr( P-C ) = R·R (sqr(x) = x^2 = x·x)

O=origin, D=direction, C=center, R=radius

射线方程表明的是如下一个点的集合P,当t从零增大时, D·t会沿着D向量的方向从零逐步变长,t 取值无限表示了射线单方向。从O点开始在D方向上无限个点构成了一条射线。
球的方程表明了任何点P,只要到C点的距离等于半径R,则表明点在球面上,这么一个球面上的点的集合。
因此当射线与球相交的时候,这个点既在射线上,又在球面上。等式射线的P(t) = 球的P成立。

联立两个方程,试着求解 t 有:


sqr( O + D·t - C ) = R·R

设 O-C=OC,有:


sqr( OC+D·t ) - R·R = 0
//展开得到如下式子
=> D·D·t·t + 2·OC·D·t + OC·OC - R·R = 0=> (D·D)·t·t + 2·(OC·D)·t + OC·OC - R·R = 0

因为 D 是单位向量有D·D = dot(D, D) = 1最后方程为:


t·t + 2·(OC·D)·t +  OC·OC - R·R = 0;

这是一个关于 t 的二次方程at^2 + bt + c = 0那么解就已经出来了:

  • t0 = -(b + √Δ) / 2a

  • t1 = -(b - √Δ) / 2a

    • a = D·D = dot(D, D) = 1;

    • b = 2·OC·D = 2·dot(OC, D);

    • c = OC·OC - R·R = dot(OC, OC) - R·R;

    • 判别式 Δ = sqr(b) - 4ac
      = 4·sqr( OC·D ) - 4·( OC·OC - R·R )
      = 4·( sqr( OC·D ) - OC·OC + R·R );
      如果判别式 Δ > 0,则表明球与射线相交。

根据以上方程,我们其中试着展开 t 的式子
t0 = -(b + √Δ) / 2a = -(b + √Δ) / 2·1
= -b/2 - √(Δ/4)
= -dot(OC, D) - √( sqr( dot(OC, D) ) - dot(OC, OC) + R·R )

求出 t 后可以根据P(t) = O + D * t 得到交点。

附chai3d中的计算代码


[html] view plain copy
  1. inline int cIntersectionSegmentSphere(const cVector3d& a_segmentPointA,  

  2.                                       const cVector3d& a_segmentPointB,  

  3.                                       const cVector3d& a_spherePos,  

  4.                                       const double& a_sphereRadius,  

  5.                                       cVector3d& a_collisionPoint0,  

  6.                                       cVector3d& a_collisionNormal0,  

  7.                                       cVector3d& a_collisionPoint1,  

  8.                                       cVector3d& a_collisionNormal1)  

  9. {  

  10.     // temp variables  

  11.     cVector3d AB, CA;  

  12.     a_segmentPointB.subr(a_segmentPointA, AB);  

  13.     a_segmentPointA.subr(a_spherePos, CA);  

  14.     double radiusSq = a_sphereRadius * a_sphereRadius;  

  15.   

  16.     double a = AB.lengthsq();  

  17.     double b = 2.0 * cDot(AB, CA);  

  18.     double c = CA.lengthsq() - radiusSq;  

  19.   

  20.     // invalid segment  

  21.     if (a == 0)  

  22.     {  

  23.         return (0);  

  24.     }  

  25.   

  26.     double d = b*b - 4*a*c;  

  27.   

  28.     // segment ray is located outside of sphere  

  29.     if (d < 0)  

  30.     {  

  31.         return (0);  

  32.     }  

  33.   

  34.     // segment ray intersects sphere  

  35.     d = sqrt(d);  

  36.     double e = 2.0 * a;  

  37.   

  38.     // compute both solutions  

  39.     double u0 = (-b + d) / e;  

  40.     double u1 = (-b - d) / e;  

  41.   

  42.     // check if the solutions are located along the segment AB  

  43.     bool valid_u0 = cContains(u0, 0.0, 1.0);  

  44.     bool valid_u1 = cContains(u1, 0.0, 1.0);  

  45.   

  46.     // two intersection points are located along segment AB  

  47.     if (valid_u0 && valid_u1)  

  48.     {  

  49.         if (u0 > u1) { cSwap(u0, u1); }  

  50.   

  51.         // compute point 0  

  52.         AB.mulr(u0, a_collisionPoint0);  

  53.         a_collisionPoint0.add(a_segmentPointA);  

  54.   

  55.         a_collisionPoint0.subr(a_spherePos, a_collisionNormal0);  

  56.         a_collisionNormal0.normalize();  

  57.   

  58.         // compute point 1  

  59.         AB.mulr(u1, a_collisionPoint1);  

  60.         a_collisionPoint1.add(a_segmentPointA);  

  61.   

  62.         a_collisionPoint1.subr(a_spherePos, a_collisionNormal1);  

  63.         a_collisionNormal1.normalize();  

  64.   

  65.         return (2);  

  66.     }  

  67.   

  68.     // one intersection point is located along segment AB  

  69.     else if (valid_u0)  

  70.     {  

  71.         // compute point 0  

  72.         AB.mulr(u0, a_collisionPoint0);  

  73.         a_collisionPoint0.add(a_segmentPointA);  

  74.   

  75.         a_collisionPoint0.subr(a_spherePos, a_collisionNormal0);  

  76.         a_collisionNormal0.normalize();  

  77.   

  78.         // check dot product to see if the intial segment point is located  

  79.         // inside the sphere.  

  80.         double dotProduct = cDot(AB, a_collisionNormal0);   //如果射线在球内部与球发生相交,就不能算作相交,保证射线运动的方向性,射线从球体出来时的交点不算作相交点  

  81.         if (dotProduct < 0.0)  

  82.         {  

  83.             return (1);  

  84.         }  

  85.         else  

  86.         {  

  87.             return (0);  

  88.         }  

  89.     }  

  90.   

  91.     // one intersection point is located along segment AB  

  92.     else if (valid_u1)  

  93.     {  

  94.         // compute point 0  

  95.         AB.mulr(u1, a_collisionPoint0);  

  96.         a_collisionPoint0.add(a_segmentPointA);  

  97.   

  98.         a_collisionPoint0.subr(a_spherePos, a_collisionNormal0);  

  99.         a_collisionNormal0.normalize();  

  100.   

  101.         // check dot product to see if the intial segment point is located  

  102.         // inside the sphere.  

  103.         double dotProduct = cDot(AB, a_collisionNormal0);  

  104.         if (dotProduct < 0.0)  

  105.         {  

  106.             return (1);  

  107.         }  

  108.         else  

  109.         {  

  110.             return (0);  

  111.         }  

  112.     }  

  113.   

  114.     // both points are located outside of the segment AB  

  115.     else  

  116.     {  

  117.         return (0);  

  118.     }  

  119. }</span>  



文/goteet(简书作者)
原文链接:http://www.jianshu.com/p/1b008ed86627
著作权归作者所有,转载请联系作者获得授权,并标注“简书作者”。


发表评论:

Powered By Z-BlogPHP 1.5.1 Zero

Copyright doorxp.com Rights Reserved.