Probably wrong bilinear patch intersection

ray-tracing-pbrt-scene/src/shape.rs

@@ -58,6 +58,103 @@ pub(crate) enum ShapeType {
     },
 }
 
+fn quadratic(a: Float, b: Float, c: Float) -> (Option<Float>, Option<Float>) {
+    if a == 0.0 {
+        if b == 0.0 {
+            return (None, None);
+        }
+        return (Some(-c / b), None);
+    }
+
+    let discrim = b * b - 4.0 * a * c;
+
+    if discrim < 0.0 {
+        return (None, None);
+    }
+
+    let root_discrim = Float::sqrt(discrim);
+
+    let q = if b < 0.0 {
+        -0.5 * (b - root_discrim)
+    } else {
+        -0.5 * (b + root_discrim)
+    };
+
+    let mut t0 = q / a;
+    let mut t1 = c / q;
+
+    if t1 < t0 {
+        std::mem::swap(&mut t1, &mut t0);
+    }
+
+    (Some(t0), Some(t1))
+}
+
+fn bilinear_intersection(
+    ray: Ray,
+    min: Float,
+    max: Float,
+    patch: [Pos3; 4],
+) -> Option<(Float, Float, Float)> {
+    let a = Dir3::dot(
+        Dir3::cross(patch[2] - patch[0], patch[1] - patch[3]),
+        ray.dir(),
+    );
+    let c = Dir3::dot(
+        Dir3::cross(patch[0] - ray.start(), ray.dir()),
+        patch[1] - patch[0],
+    );
+    let b = Dir3::dot(
+        Dir3::cross(patch[2] - ray.start(), ray.dir()),
+        patch[3] - patch[2],
+    ) - (a + c);
+
+    let (u0, u1) = quadratic(a, b, c);
+
+    [u0, u1]
+        .into_iter()
+        .flatten()
+        .filter(|&u| {
+            if 0.0 <= u && u <= 1.0 {
+                return true;
+            }
+            false
+        })
+        .filter_map(|u| {
+            let uo = patch[0].as_dir() * u + patch[2].as_dir() * (1.0 - u);
+            let ud = patch[1].as_dir() * u + patch[3].as_dir() * (1.0 - u) - uo;
+
+            let deltao = uo - ray.start().as_dir();
+
+            let perp = Dir3::cross(ray.dir(), ud);
+
+            let p2 = perp.length_squared();
+
+            let v1 = deltao.x() * ray.dir().y()
+                + perp.z()
+                + ray.dir().x() * perp.y() * deltao.z()
+                + perp.x() * deltao.y() * ray.dir().z()
+                - deltao.x() * perp.y() * deltao.z()
+                - ray.dir().x() * deltao.y() * perp.z()
+                - perp.x() * ray.dir().y() * deltao.z();
+
+            let t1 = deltao.x() * ud.y()
+                + perp.z()
+                + ud.x() * perp.y() * deltao.z()
+                + perp.x() * deltao.y() * ud.z()
+                - deltao.x() * perp.y() * deltao.z()
+                - ud.x() * deltao.y() * perp.z()
+                - perp.x() * ud.y() * deltao.z();
+
+            if t1 > 0.0 && 0.0 < v1 && v1 <= p2 {
+                Some((t1 / p2, u, v1 / p2))
+            } else {
+                None
+            }
+        })
+        .next()
+}
+
 impl Shape {
     pub(crate) fn intersect<R: Rng>(
         &self,
@@ -100,23 +197,19 @@ impl Shape {
                         std::mem::swap(&mut t1, &mut t0);
                     }
 
-                    if let Some(t) = [t0, t1]
-                        .into_iter()
-                        .filter(|&t| {
-                            if min <= t && t <= max {
-                                let p = ray.at(t0);
-                                let mut phi = Float::atan2(p.y(), p.x());
-                                if phi < 0.0 {
-                                    phi += 2.0 * FloatConsts::PI;
-                                }
-                                if zmin <= p.z() && p.z() <= zmax && phi <= phimax {
-                                    return true;
-                                }
+                    if let Some(t) = [t0, t1].into_iter().find(|&t| {
+                        if min <= t && t <= max {
+                            let p = ray.at(t0);
+                            let mut phi = Float::atan2(p.y(), p.x());
+                            if phi < 0.0 {
+                                phi += 2.0 * FloatConsts::PI;
                             }
-                            false
-                        })
-                        .next()
-                    {
+                            if zmin <= p.z() && p.z() <= zmax && phi <= phimax {
+                                return true;
+                            }
+                        }
+                        false
+                    }) {
                         return Some(Intersection::new(
                             t,
                             (ray.at(t) - Pos3::zero()).normalize(),
@@ -134,7 +227,14 @@ impl Shape {
                 s,
                 uv,
             } => (),
-            ShapeType::BilinearMesh { indices, p, n, uv } => (),
+            ShapeType::BilinearMesh { indices, p, n, uv } => {
+                if indices.is_empty()
+                    && let Some((t, _u, _v)) =
+                        bilinear_intersection(ray, min, max, [p[0], p[1], p[2], p[3]])
+                {
+                    return Some(Intersection::new(t, Dir3::up(), None, None, 0.0));
+                }
+            }
             ShapeType::LoopSubDiv { levels, indices, p } => todo!(),
             ShapeType::Disk {
                 height,