射线到线段最短距离, 以及最近点的计算方式

[tiye]

在 Soluble 当中经常遇到, 虽然是类似 SDF 计算每个点的着色, 但是 Soluble 简化了图案, 比较简单.
这样, 就很容易遇到了一个计算射线距离线段最短的距离, 比如这个图案当中,

https://webgpu.art/soluble/?tab=parallel-mirror

为了计算光晕的强度, 就需要知道射线到线段的距离, 有两种可能,

• 两者所在直线的距离, 就是最短的距离. 这个比较好计算, 两者方向向量计算叉积, 得到垂直方向, 然后选个线段, 计算点积就得到了长度.
• 或者, 线段端点到直线的距离是最近的, 点到直线的距离, 任意选取点, 计算线段和射线方向的点积, 可以求出点到直线最近的点, 然后也得到了长度. 两个端点分别求出来, 取小的.
• 实际当中还有另一种可能, 因为是射线, 最近距离可能是在背面... 这个情况我这里只是粗略估算了. 按照两个端点看, 只要有一个在背面, 就不管了.

实际当中取的是直线间最短, 还是端点到直线最短, 是比较麻烦的, 这个情况就需要计算最近点是不是在线段内,
或者说, 最近的点, 到两个端点的向量, 是不是反向的?
这个也就涉及到了两直线距离最近点的具体位置了. 此外, 也可以推算直线射线到最短点, 经过了多长距离.

我这里没有完全准确得计算方案, 大致有思路, 基于 Gemini 的提示再处理成方程算了一下.
大致得到一个结果, 效果看开头给的图像... 这里记录一下大致步骤.

/// 返回的结构, 射线距离线段
struct RayReachSegment {
  /// 距离
  distance: f32,
  /// 是否在射线前方
  positive_side: bool,
  /// 射线经过多少距离, 达到两者最近位置
  traveled: f32,
}

/// 线段
struct Segment {
  start: vec3<f32>,
  end: vec3<f32>,
}

/// 计算射线到线段最短位置
fn ray_nearest_point_to_line(
  /** 射线起点, 也就是观察者 */ viewer_position: vec3f,
  /** 射线单位向量 */ray_unit: vec3f,
  /** 线段 */ s: Segment) -> RayReachSegment {
  /** 线段 A 端距离观察者 */
  let a = s.start - viewer_position;
  /** 线段 B 端距离观察者 */
  let b = s.end - viewer_position;

  /** 找到一个单位向量, 跟射线垂直, 也跟线段垂直(包含长度信息) *//
  let n = cross(b - a, ray_unit);

  // find projection of a of segment on ray direction, and use the Pythagorean theorem for another distance
  /** 计算观察者距离 A 端点在射线方向长度 */
  let a_proj = dot(ray_unit, a);
  /** 从最近点到 A 点方向 */
  let shadow_a = a - ray_unit * a_proj;

  /** 计算观察者距离 A 端点在射线方向长度 */
  let b_proj = dot(ray_unit, b);
  /** 从最近点到 A 点方向 */
  let shadow_b = b - ray_unit * b_proj;

  // 这里注意 n 方向垂直两条线, 也就是两者最近点之间的方向向量
  /** 叉积计算基于最近点往 A 跟最近方向的方向关系 */
  let direct_an = cross(shadow_a, n);
  /** 叉积计算基于最近点往 B 跟最近方向的方向关系 */
  let direct_bn = cross(shadow_b, n);

  /** A B 在最近点同侧的情况下, 叉积产生的正负号正负号应该是一直的, 乘积为正 */
  let same_side = dot(direct_an, direct_bn) >= 0.0;

  if same_side {
   // 那么意味着, 两直线最近位置, 没有落在线段上,, 那取端点近的一个就好了
   
   /** 直接勾股定理得到长度 */
    let a_distance_min = sqrt(dot(a, a) - a_proj * a_proj);
    let b_distance_min = sqrt(dot(b, b) - b_proj * b_proj);
    if a_distance_min < b_distance_min {
      return RayReachSegment(a_distance_min, a_proj > 0., a_proj);
    } else {
      return RayReachSegment(b_distance_min, b_proj > 0., b_proj);
    }
  } else {
    // 最近位置落在线段上, 这时要计算直线之间具体的最近位置了

    let n0 = normalize(n);
    /// smaller distance to segments ends
    /**  任取一点, 计算跟两直线垂直的方向得点积, 就是距离了, 这个没什么障碍. 最多正负号麻烦点 */
    let d_min = abs(dot(n0, a));

    /** 线段取一个单位方向 */
    let ab_unit = normalize(s.end - s.start);
    /** 其实还是线段到观察者 */
    let ac = (viewer_position - s.start);

    /** 换算最近距离的线段信息, 进度和方向 */
    let perp_reach: vec3f = dot(ac, n0) * n0;

    /** 这部分比较绕了, 方程当中我设定了一个 k 表示射线到最近点的距离
     上边有两者距离的线段的信息, 那么, 比如 AB 线段, 经过平移, 理论上跟 CD 会相交的, 我认为在 D 点相交(对应后面图片)
      https://cos-sh.tiye.me/cos-up/342e0f6b7e3a7b1b4d18dad16cd84f79/IMG_20240601_113231.jpg
      于是就可以联立一个公式, A 平移过来的 A'', 跟 D 点组成线段, 会跟 AB 本身是平行的, 叉积为 0.
      通过两边点积的方式消元, 可以得到一个 k, 我的步骤不清楚, 反而不如 Gemini 给出来的.
      https://g.co/gemini/share/dbf766e185a7
      最后就这么个式子(不严谨) */
    let k = dot(
      cross(s.start + perp_reach - viewer_position, ab_unit),
      cross(ray_unit, ab_unit)
    ) + .0;
    
    /** 判断线段在前在后, 比较粗略了 */
    // let front = dot(a, ray_unit) >= 0.0 && dot(b, ray_unit) >= 0.0;
    let front = k > 0.;
    return RayReachSegment(d_min, front, k);
  };
}

仔细看的话, 这里计算 k 的过程, 第一步的平移就有问题, 平移的方向是不确定的, 所以正负号不知道.
所以这个计算步骤, 只是实际使用当中没有直接撞上问题, 但并不准确.

由于实现这个之前, 网上到处找公式比较费劲, 所以把已知的计算步骤记录一下. 方便查阅.
也希望有人能整理更清楚的计算过程, 还有补充图示.
万一以后某天我真的很有时间了, 那是还没有人补的话, 那也只能我自己补一个了.