use super::{
    BilinearIntersection, BilinearPatchShape, Bounds3f, DirectionCone, Interaction, Normal3f,
    Point2f, Point3f, Point3fi, Ray, ShapeIntersection, ShapeSample, ShapeSampleContext,
    ShapeTrait, SurfaceInteraction, Transform, Vector3f,
};
use crate::core::geometry::{Tuple, VectorLike, spherical_quad_area};
use crate::core::interaction::InteractionTrait;
use crate::core::pbrt::{Float, gamma};
use crate::utils::math::{SquareMatrix, clamp, difference_of_products, lerp, quadratic};
use crate::utils::mesh::BilinearPatchMesh;
use crate::utils::sampling::{
    bilinear_pdf, invert_spherical_rectangle_sample, sample_bilinear, sample_spherical_rectangle,
};
use std::sync::Arc;
use std::sync::OnceLock;

struct PatchData<'a> {
    mesh: &'a BilinearPatchMesh,
    p00: Point3f,
    p10: Point3f,
    p01: Point3f,
    p11: Point3f,
    n: Option<[Normal3f; 4]>,
    uv: Option<[Point2f; 4]>,
}

struct IntersectionData {
    t: Float,
    u: Float,
    v: Float,
}

struct TextureDerivative {
    duds: Float,
    dvds: Float,
    dudt: Float,
    dvdt: Float,
}

static BILINEAR_MESHES: OnceLock<Vec<Arc<BilinearPatchMesh>>> = OnceLock::new();
impl BilinearPatchShape {
    pub const MIN_SPHERICAL_SAMPLE_AREA: Float = 1e-4;

    pub fn new(_mesh: BilinearPatchMesh, mesh_index: usize, blp_index: usize) -> Self {
        let mut bp = BilinearPatchShape {
            mesh_index,
            blp_index,
            area: 0.,
            rectangle: false,
        };

        let (p00, p10, p01, p11) = {
            let data = bp.get_data();
            (data.p00, data.p10, data.p01, data.p11)
        };

        bp.rectangle = bp.is_rectangle(p00, p10, p01, p11);

        if bp.rectangle {
            bp.area = p00.distance(p01) + p00.distance(p10);
        } else {
            const NA: usize = 3;
            let mut p = [[Point3f::default(); NA + 1]; NA + 1];
            #[allow(clippy::needless_range_loop)]
            for i in 0..=NA {
                let u = i as Float / NA as Float;
                for j in 0..=NA {
                    let v = j as Float / NA as Float;
                    p[i][j] = lerp(u, lerp(v, p00, p01), lerp(v, p10, p11));
                }
            }

            let mut area = 0.;
            for i in 0..NA {
                for j in 0..NA {
                    let diag1 = p[i + 1][j + 1] - p[i][j];
                    let diag2 = p[i + 1][j] - p[i][j + 1];
                    let quad_area = 0.5 * diag1.cross(diag2).norm();
                    area += quad_area;
                }
            }
            bp.area = area;
        }
        bp
    }

    fn mesh(&self) -> &Arc<BilinearPatchMesh> {
        let meshes = BILINEAR_MESHES
            .get()
            .expect("Mesh has not been initialized");
        &meshes[self.mesh_index]
    }

    fn get_data(&self) -> PatchData<'_> {
        let mesh = self.mesh();
        let start_index = 4 * self.blp_index;
        let v = &mesh.vertex_indices[start_index..start_index + 4];
        let p00: Point3f = mesh.p[v[0]];
        let p10: Point3f = mesh.p[v[1]];
        let p01: Point3f = mesh.p[v[2]];
        let p11: Point3f = mesh.p[v[3]];
        let n = mesh
            .n
            .as_ref()
            .map(|normals| [normals[v[0]], normals[v[1]], normals[v[2]], normals[v[3]]]);
        let uv = mesh
            .uv
            .as_ref()
            .map(|uvs| [uvs[v[0]], uvs[v[1]], uvs[v[2]], uvs[v[3]]]);
        PatchData {
            mesh,
            p00,
            p10,
            p01,
            p11,
            n,
            uv,
        }
    }

    fn is_rectangle(&self, p00: Point3f, p10: Point3f, p01: Point3f, p11: Point3f) -> bool {
        if p00 == p01 || p01 == p11 || p11 == p10 || p10 == p00 {
            return false;
        }

        let n = Normal3f::from((p10 - p00).cross(p01 - p00).normalize());
        if (p11 - p00).normalize().dot(n.into()).abs() > 1e-5 {
            return false;
        }

        let p_center_vec = Vector3f::from(p00 + p01.into() + p10.into() + p11.into()) / 4.;
        let p_center = Point3f::from(p_center_vec);
        let d2 = [
            p00.distance_squared(p_center),
            p01.distance_squared(p_center),
            p10.distance_squared(p_center),
            p11.distance_squared(p_center),
        ];

        for i in 0..4 {
            if (d2[i] - d2[0]).abs() / d2[0] > 1e-4 {
                return false;
            }
        }
        true
    }

    fn intersect_bilinear_patch(
        &self,
        ray: &Ray,
        t_max: Float,
        data: &PatchData,
    ) -> Option<BilinearIntersection> {
        let a = (data.p10 - data.p00).cross(data.p01 - data.p11).dot(ray.d);
        let c = (data.p00 - ray.o).cross(ray.d).dot(data.p01 - data.p00);
        let b = (data.p10 - ray.o).cross(ray.d).dot(data.p11 - data.p10) - (a + c);

        let (u1, u2) = quadratic(a, b, c)?;

        let eps = gamma(10)
            * (ray.o.abs().max_component_value()
                + ray.d.abs().max_component_value()
                + data.p00.abs().max_component_value()
                + data.p10.abs().max_component_value()
                + data.p01.abs().max_component_value()
                + data.p11.abs().max_component_value());

        let hit1 = self.check_candidate(u1, ray, data);
        let hit2 = if u1 != u2 {
            self.check_candidate(u2, ray, data)
        } else {
            None
        };

        // Find best hit from candidates
        [hit1, hit2]
            .into_iter()
            .flatten()
            .filter(|hit| hit.t < t_max && hit.t > eps)
            .min_by(|a, b| a.t.partial_cmp(&b.t).unwrap())
            .map(|best_hit| {
                BilinearIntersection::new(Point2f::new(best_hit.u, best_hit.v), best_hit.t)
            })
    }

    fn check_candidate(&self, u: Float, ray: &Ray, data: &PatchData) -> Option<IntersectionData> {
        if !(0.0..=1.0).contains(&u) {
            return None;
        }
        let uo: Point3f = lerp(u, data.p00, data.p10);
        let ud: Point3f = Point3f::from(lerp(u, data.p01, data.p11) - uo);
        let deltao = uo - ray.o;
        let perp = ray.d.cross(ud.into());
        let p2 = perp.norm_squared();
        if p2 == 0. {
            return None;
        }
        let v_det = SquareMatrix::<Float, 3>::new([
            [deltao.x(), ray.d.x(), perp.x()],
            [deltao.y(), ray.d.y(), perp.y()],
            [deltao.z(), ray.d.z(), perp.z()],
        ])
        .determinant();

        if !(0.0..=p2).contains(&v_det) {
            return None;
        }

        let t_det = SquareMatrix::<Float, 3>::new([
            [deltao.x(), ud.x(), perp.x()],
            [deltao.y(), ud.y(), perp.y()],
            [deltao.z(), ud.z(), perp.z()],
        ])
        .determinant();

        Some(IntersectionData {
            t: t_det / p2,
            u,
            v: v_det / p2,
        })
    }

    fn interaction_from_intersection(
        &self,
        data: &PatchData,
        uv: Point2f,
        time: Float,
        wo: Vector3f,
    ) -> SurfaceInteraction {
        // Base geom and derivatives
        let p = lerp(
            uv[0],
            lerp(uv[1], data.p00, data.p01),
            lerp(uv[1], data.p10, data.p11),
        );
        let mut dpdu = lerp(uv[1], data.p10, data.p11) - lerp(uv[1], data.p00, data.p01);
        let mut dpdv = lerp(uv[0], data.p01, data.p11) - lerp(uv[0], data.p00, data.p10);

        // If textured, apply coordinates
        let (st, derivatives) = self.apply_texture_coordinates(data, uv, &mut dpdu, &mut dpdv);

        // Compute second moments
        let n = Normal3f::from(dpdu.cross(dpdv).normalize());
        let (mut dndu, mut dndv) = self.calculate_surface_curvature(data, &dpdu, &dpdv, n);
        if let Some(ref deriv) = derivatives {
            let dnds = Normal3f::from(dndu * deriv.duds + dndv * deriv.dvds);
            let dndt = Normal3f::from(dndu * deriv.dudt + dndv * deriv.dvdt);
            dndu = dnds;
            dndv = dndt;
        }

        let p_abs_sum = data.p00.abs()
            + Vector3f::from(data.p01.abs())
            + Vector3f::from(data.p10.abs())
            + Vector3f::from(data.p11.abs());
        let p_error = gamma(6) * Vector3f::from(p_abs_sum);

        let flip_normal = data.mesh.reverse_orientation ^ data.mesh.transform_swaps_handedness;
        let pi = Point3fi::new_with_error(p, p_error);
        let mut isect =
            SurfaceInteraction::new(pi, st, wo, dpdu, dpdv, dndu, dndv, time, flip_normal);
        if data.n.is_some() {
            self.apply_shading_normals(&mut isect, data, uv, derivatives);
        }
        isect
    }

    fn apply_texture_coordinates(
        &self,
        data: &PatchData,
        uv: Point2f,
        dpdu: &mut Vector3f,
        dpdv: &mut Vector3f,
    ) -> (Point2f, Option<TextureDerivative>) {
        let Some(uvs) = data.uv else {
            return (uv, None);
        };
        let uv00 = uvs[0];
        let uv01 = uvs[1];
        let uv10 = uvs[2];
        let uv11 = uvs[3];
        // Compute partial derivatives of (u, v) with respect to (s, t)
        let st = lerp(uv[0], lerp(uv[1], uv00, uv01), lerp(uv[1], uv10, uv11));
        let dstdu = lerp(uv[1], uv10, uv11) - lerp(uv[1], uv00, uv01);
        let dstdv = lerp(uv[0], uv01, uv11) - lerp(uv[0], uv00, uv10);

        let det = difference_of_products(dstdu[0], dstdv[1], dstdu[1], dstdv[0]);
        let inv_det = if det == 0.0 { 0.0 } else { 1.0 / det };
        let duds = dstdv[1] * inv_det;
        let dvds = -dstdv[0] * inv_det;
        let dudt = -dstdu[1] * inv_det;
        let dvdt = dstdu[0] * inv_det;

        let dpds = *dpdu * duds + *dpdv * dvds;
        let mut dpdt = *dpdu * dudt + *dpdv * dvdt;

        if dpdu.cross(*dpdv).dot(dpds.cross(dpdt)) < 0. {
            dpdt = -dpdt;
        }
        *dpdu = dpds;
        *dpdv = dpdt;

        let factors = TextureDerivative {
            duds,
            dvds,
            dudt,
            dvdt,
        };
        (st, Some(factors))
    }

    fn calculate_base_derivatives(
        &self,
        data: &PatchData,
        uv: Point2f,
    ) -> (Point3f, Vector3f, Vector3f) {
        let p = lerp(
            uv[0],
            lerp(uv[1], data.p00, data.p01),
            lerp(uv[1], data.p10, data.p11),
        );
        let dpdu = lerp(uv[1], data.p10, data.p11) - lerp(uv[1], data.p00, data.p01);
        let dpdv = lerp(uv[0], data.p01, data.p11) - lerp(uv[0], data.p00, data.p10);
        (p, dpdu, dpdv)
    }

    fn calculate_surface_curvature(
        &self,
        data: &PatchData,
        dpdu: &Vector3f,
        dpdv: &Vector3f,
        n: Normal3f,
    ) -> (Normal3f, Normal3f) {
        let e = dpdu.dot(*dpdu);
        let f = dpdu.dot(*dpdv);
        let g = dpdv.dot(*dpdv);
        let d2pduv = (data.p00 - data.p01) + (data.p11 - data.p10);
        let d2pduu = Vector3f::zero();
        let d2pdvv = Vector3f::zero();

        let e_min = n.dot(d2pduu.into());
        let f_min = n.dot(d2pduv.into());
        let g_min = n.dot(d2pdvv.into());

        let egf2 = difference_of_products(e, g, f, f);
        let inv_egf2 = if egf2 == 0. { 0. } else { 1. / egf2 };

        let dndu = Normal3f::from(
            (f_min * f - e_min * g) * inv_egf2 * *dpdu + (e_min * f - f_min * e) * inv_egf2 * *dpdv,
        );
        let dndv = Normal3f::from(
            (g_min * f - f_min * g) * inv_egf2 * *dpdu + (f_min * f - g_min * e) * inv_egf2 * *dpdv,
        );
        (dndu, dndv)
    }

    fn apply_shading_normals(
        &self,
        isect: &mut SurfaceInteraction,
        data: &PatchData,
        uv: Point2f,
        derivatives: Option<TextureDerivative>,
    ) {
        let Some(normals) = data.n else { return };
        let n00 = normals[1];
        let n10 = normals[1];
        let n01 = normals[2];
        let n11 = normals[3];
        let ns = lerp(uv[0], lerp(uv[1], n00, n01), lerp(uv[1], n10, n11)).normalize();
        let mut dndu = lerp(uv[1], n10, n11) - lerp(uv[1], n00, n01);
        let mut dndv = lerp(uv[0], n01, n11) - lerp(uv[0], n00, n10);

        if let Some(deriv) = derivatives {
            let dnds = dndu * deriv.duds + dndv * deriv.dvds;
            let dndt = dndu * deriv.dudt + dndv * deriv.dvdt;
            dndu = dnds;
            dndv = dndt;
        }
        let r = Transform::rotate_from_to(isect.n().normalize().into(), ns.into());

        isect.set_shading_geom(
            ns,
            r.apply_to_vector(isect.dpdu),
            r.apply_to_vector(isect.dpdv),
            dndu,
            r.apply_to_normal(dndv),
            true,
        );
    }

    fn sample_area_and_pdf(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option<ShapeSample> {
        let mut ss = self.sample(u)?;

        let mut intr_clone = (*ss.intr).clone();
        intr_clone.common.time = ctx.time;
        ss.intr = Arc::new(intr_clone);

        let mut wi = ss.intr.p() - ctx.p();
        let dist_sq = wi.norm_squared();
        if dist_sq == 0. {
            return None;
        }

        wi = wi.normalize();
        let abs_dot = Vector3f::from(ss.intr.n()).abs_dot(-wi);
        if abs_dot == 0. {
            return None;
        }

        ss.pdf *= dist_sq / abs_dot;

        if ss.pdf.is_infinite() { None } else { Some(ss) }
    }

    fn sample_parametric_coords(&self, data: &PatchData, u: Point2f) -> (Point2f, Float) {
        if let Some(image_distrib) = &data.mesh.image_distribution {
            let (uv, pdf, _) = image_distrib.sample(u);
            (uv, pdf)
        } else if !self.rectangle {
            let w = [
                (data.p10 - data.p00).cross(data.p01 - data.p00).norm(),
                (data.p10 - data.p00).cross(data.p11 - data.p10).norm(),
                (data.p01 - data.p00).cross(data.p11 - data.p01).norm(),
                (data.p11 - data.p10).cross(data.p11 - data.p01).norm(),
            ];
            let uv = sample_bilinear(u, &w);
            let pdf = bilinear_pdf(uv, &w);
            (uv, pdf)
        } else {
            (u, 1.0)
        }
    }

    fn sample_solid_angle(
        &self,
        data: &PatchData,
        ctx: &ShapeSampleContext,
        u: Point2f,
        corner_dirs: &[Vector3f; 4],
    ) -> Option<ShapeSample> {
        let mut pdf = 1.;
        if ctx.ns != Normal3f::zero() {
            let w = [
                0.01_f32.max(corner_dirs[0].dot(ctx.ns.into()).abs()),
                0.01_f32.max(corner_dirs[1].dot(ctx.ns.into()).abs()),
                0.01_f32.max(corner_dirs[2].dot(ctx.ns.into()).abs()),
                0.01_f32.max(corner_dirs[3].dot(ctx.ns.into()).abs()),
            ];
            let u = sample_bilinear(u, &w);
            pdf *= bilinear_pdf(u, &w);
        }

        let eu = data.p10 - data.p00;
        let ev = data.p01 - data.p00;
        let (p, quad_pdf) = sample_spherical_rectangle(ctx.p(), data.p00, eu, ev, u);
        pdf *= quad_pdf?;

        // Compute (u, v) and surface normal for sampled points on rectangle
        let uv = Point2f::new(
            (p - data.p00).dot(eu) / data.p10.distance_squared(data.p00),
            (p - data.p00).dot(ev) / data.p01.distance_squared(data.p00),
        );

        let n = self.compute_sampled_normal(data, &eu, &ev, uv);
        let st = data.uv.map_or(uv, |uvs| {
            lerp(
                uv[0],
                lerp(uv[1], uvs[0], uvs[1]),
                lerp(uv[1], uvs[2], uvs[3]),
            )
        });

        let pi = Point3fi::new_from_point(p);
        let mut intr = SurfaceInteraction::new_simple(pi, n, st);
        intr.common.time = ctx.time;
        Some(ShapeSample {
            intr: Arc::new(intr),
            pdf,
        })
    }

    fn compute_sampled_normal(
        &self,
        data: &PatchData,
        dpdu: &Vector3f,
        dpdv: &Vector3f,
        uv: Point2f,
    ) -> Normal3f {
        let mut n = Normal3f::from(dpdu.cross(*dpdv).normalize());

        if let Some(normals) = data.n {
            // Apply interpolated shading normal to orient the geometric normal
            let ns = lerp(
                uv[0],
                lerp(uv[1], normals[0], normals[2]),
                lerp(uv[1], normals[1], normals[3]),
            );
            n = n.face_forward(ns.into());
        } else if data.mesh.reverse_orientation ^ data.mesh.transform_swaps_handedness {
            n = -n;
        }
        n
    }
}

impl ShapeTrait for BilinearPatchShape {
    #[inline]
    fn area(&self) -> Float {
        self.area
    }

    #[inline]
    fn normal_bounds(&self) -> DirectionCone {
        let data = self.get_data();
        if data.p00 == data.p10
            || data.p10 == data.p11
            || data.p11 == data.p01
            || data.p01 == data.p00
        {
            let dpdu = lerp(0.5, data.p10, data.p11) - lerp(0.5, data.p00, data.p01);
            let dpdv = lerp(0.5, data.p01, data.p11) - lerp(0.5, data.p00, data.p10);
            let mut n = Normal3f::from(dpdu.cross(dpdv).normalize());
            if let Some(normals) = data.n {
                let interp_n = (normals[0] + normals[1] + normals[2] + normals[3]) / 4.;
                n = n.face_forward(interp_n.into());
            } else if data.mesh.reverse_orientation ^ data.mesh.transform_swaps_handedness {
                n *= -1.;
            }
            return DirectionCone::new_from_vector(Vector3f::from(n));
        }

        // Compute bilinear patch normals n10, n01, and n11
        let mut n00 = Normal3f::from((data.p10 - data.p00).cross(data.p01 - data.p00).normalize());
        let mut n10 = Normal3f::from((data.p11 - data.p10).cross(data.p00 - data.p10).normalize());
        let mut n01 = Normal3f::from((data.p00 - data.p01).cross(data.p11 - data.p01).normalize());
        let mut n11 = Normal3f::from((data.p01 - data.p11).cross(data.p10 - data.p11).normalize());

        if let Some(normals) = data.n {
            n00 = n00.face_forward(normals[0].into());
            n10 = n10.face_forward(normals[1].into());
            n01 = n01.face_forward(normals[2].into());
            n11 = n11.face_forward(normals[3].into());
        } else if data.mesh.reverse_orientation ^ data.mesh.transform_swaps_handedness {
            n00 = -n00;
            n10 = -n10;
            n01 = -n01;
            n11 = -n11;
        }

        // Compute average normal and return normal bounds for patch
        let n_avg = (n00 + n10 + n01 + n11).normalize();
        let cos_theta = n_avg
            .dot(n00)
            .min(n_avg.dot(n10))
            .min(n_avg.dot(n01))
            .min(n_avg.dot(n11));
        DirectionCone::new(n_avg.into(), clamp(cos_theta, -1., 1.))
    }

    #[inline]
    fn bounds(&self) -> Bounds3f {
        let data = self.get_data();
        Bounds3f::from_points(data.p00, data.p01).union(Bounds3f::from_points(data.p10, data.p11))
    }

    #[inline]
    fn intersect(&self, ray: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
        let t_max_val = t_max?;
        let data = self.get_data();
        if let Some(bilinear_hit) = self.intersect_bilinear_patch(ray, t_max_val, &data) {
            let intr = self.interaction_from_intersection(&data, bilinear_hit.uv, ray.time, -ray.d);

            Some(ShapeIntersection {
                intr,
                t_hit: bilinear_hit.t,
            })
        } else {
            None
        }
    }

    #[inline]
    fn intersect_p(&self, ray: &Ray, t_max: Option<Float>) -> bool {
        let t_max_val = t_max.unwrap_or(Float::INFINITY);
        let data = self.get_data();
        self.intersect_bilinear_patch(ray, t_max_val, &data)
            .is_some()
    }

    #[inline]
    fn sample(&self, u: Point2f) -> Option<ShapeSample> {
        let data = self.get_data();
        // Sample bilinear patch parametric coordinate (u, v)
        let (uv, pdf) = self.sample_parametric_coords(&data, u);
        // Compute bilinear patch geometric quantities at sampled (u, v)
        let (p, dpdu, dpdv) = self.calculate_base_derivatives(&data, uv);
        if dpdu.norm_squared() == 0. || dpdv.norm_squared() == 0. {
            return None;
        }

        // Compute surface normal for sampled bilinear patch (u, v)
        let n = self.compute_sampled_normal(&data, &dpdu, &dpdv, uv);
        let st = data.uv.map_or(uv, |patch_uvs| {
            lerp(
                uv[0],
                lerp(uv[1], patch_uvs[0], patch_uvs[1]),
                lerp(uv[1], patch_uvs[2], patch_uvs[3]),
            )
        });

        let p_abs_sum = data.p00.abs()
            + Vector3f::from(data.p01.abs())
            + Vector3f::from(data.p10.abs())
            + Vector3f::from(data.p11.abs());
        let p_error = gamma(6) * Vector3f::from(p_abs_sum);
        let pi = Point3fi::new_with_error(p, p_error);
        Some(ShapeSample {
            intr: Arc::new(SurfaceInteraction::new_simple(pi, n, st)),
            pdf: pdf / dpdu.cross(dpdv).norm(),
        })
    }

    #[inline]
    fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option<ShapeSample> {
        let data = self.get_data();
        let v00 = (data.p00 - ctx.p()).normalize();
        let v10 = (data.p10 - ctx.p()).normalize();
        let v01 = (data.p01 - ctx.p()).normalize();
        let v11 = (data.p11 - ctx.p()).normalize();

        let use_area_sampling = self.rectangle
            || data.mesh.image_distribution.is_some()
            || spherical_quad_area(v00, v10, v11, v01) <= Self::MIN_SPHERICAL_SAMPLE_AREA;
        if use_area_sampling {
            self.sample_area_and_pdf(ctx, u)
        } else {
            self.sample_solid_angle(&data, ctx, u, &[v00, v10, v01, v11])
        }
    }

    #[inline]
    fn pdf(&self, intr: &Interaction) -> Float {
        let Interaction::Surface(si) = intr else {
            return 0.0;
        };

        let data = self.get_data();
        let uv = if let Some(uvs) = &data.mesh.uv {
            Point2f::invert_bilinear(si.uv, uvs)
        } else {
            si.uv
        };

        let param_pdf = if let Some(image_distrib) = &data.mesh.image_distribution {
            image_distrib.pdf(uv)
        } else if self.rectangle {
            let w = [
                (data.p10 - data.p00).cross(data.p01 - data.p00).norm(),
                (data.p10 - data.p00).cross(data.p11 - data.p10).norm(),
                (data.p01 - data.p00).cross(data.p11 - data.p01).norm(),
                (data.p11 - data.p10).cross(data.p11 - data.p01).norm(),
            ];
            bilinear_pdf(uv, &w)
        } else {
            1.
        };

        let (_, dpdu, dpdv) = self.calculate_base_derivatives(&data, uv);
        let cross = dpdu.cross(dpdv).norm();
        if cross == 0. { 0. } else { param_pdf / cross }
    }

    #[inline]
    fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float {
        let ray = ctx.spawn_ray(wi);
        let Some(isect) = self.intersect(&ray, None) else {
            return 0.;
        };

        let data = self.get_data();

        let v00 = (data.p00 - ctx.p()).normalize();
        let v10 = (data.p10 - ctx.p()).normalize();
        let v01 = (data.p01 - ctx.p()).normalize();
        let v11 = (data.p11 - ctx.p()).normalize();

        let use_area_sampling = !self.rectangle
            || data.mesh.image_distribution.is_some()
            || spherical_quad_area(v00, v10, v01, v11) <= Self::MIN_SPHERICAL_SAMPLE_AREA;

        if use_area_sampling {
            let intr_wrapper = Interaction::Surface(isect.intr.clone());
            let isect_pdf = self.pdf(&intr_wrapper);
            let distsq = ctx.p().distance_squared(isect.intr.p());
            let absdot = Vector3f::from(isect.intr.n()).abs_dot(-wi);
            if absdot == 0. {
                return 0.;
            }
            let pdf = isect_pdf * distsq / absdot;
            if pdf.is_infinite() { 0. } else { pdf }
        } else {
            let mut pdf = 1. / spherical_quad_area(v00, v10, v01, v11);
            if ctx.ns != Normal3f::zero() {
                let w = [
                    0.01_f32.max(v00.dot(ctx.ns.into()).abs()),
                    0.01_f32.max(v10.dot(ctx.ns.into()).abs()),
                    0.01_f32.max(v01.dot(ctx.ns.into()).abs()),
                    0.01_f32.max(v11.dot(ctx.ns.into()).abs()),
                ];
                let u = invert_spherical_rectangle_sample(
                    ctx.p(),
                    data.p00,
                    data.p10 - data.p00,
                    data.p01 - data.p00,
                    isect.intr.p(),
                );
                pdf *= bilinear_pdf(u, &w);
            }
            pdf
        }
    }
}