use crate::core::geometry::{
    Bounds3f, DirectionCone, Normal, Normal3f, Point2f, Point3f, Point3fi, Ray, Tuple, Vector3f,
    VectorLike, spherical_quad_area,
};
use crate::core::interaction::{Interaction, InteractionTrait, SurfaceInteraction};
use crate::core::pbrt::{Float, gamma};
use crate::core::shape::{Shape, ShapeIntersection, ShapeSample, ShapeSampleContext, ShapeTrait};
use crate::utils::Ptr;
use crate::utils::Transform;
use crate::utils::math::{SquareMatrix, clamp, difference_of_products, lerp, quadratic};
use crate::utils::mesh::DeviceBilinearPatchMesh;
use crate::utils::sampling::{
    bilinear_pdf, invert_spherical_rectangle_sample, sample_bilinear, sample_spherical_rectangle,
};
use core::ops::Add;

#[repr(C)]
#[derive(Debug, Copy, Clone)]
struct IntersectionData {
    pub t: Float,
    pub u: Float,
    pub v: Float,
}

#[repr(C)]
#[derive(Debug, Default, Copy, Clone)]
struct TextureDerivative {
    pub duds: Float,
    pub dvds: Float,
    pub dudt: Float,
    pub dvdt: Float,
}

#[repr(C)]
#[derive(Debug, Clone, Copy)]
pub struct BilinearIntersection {
    pub uv: Point2f,
    pub t: Float,
}

impl BilinearIntersection {
    pub fn new(uv: Point2f, t: Float) -> Self {
        Self { uv, t }
    }
}

#[repr(C)]
#[derive(Debug, Clone, Copy)]
pub struct BilinearPatchShape {
    pub mesh: Ptr<DeviceBilinearPatchMesh>,
    pub blp_index: i32,
    pub area: Float,
    pub rectangle: bool,
}

impl BilinearPatchShape {
    pub const MIN_SPHERICAL_SAMPLE_AREA: Float = 1e-4;
    fn mesh(&self) -> Ptr<DeviceBilinearPatchMesh> {
        self.mesh
    }

    #[inline(always)]
    fn get_vertex_indices(&self) -> [usize; 4] {
        unsafe {
            let base_ptr = self.mesh.vertex_indices.add((self.blp_index as usize) * 4);
            [
                *base_ptr.add(0) as usize,
                *base_ptr.add(1) as usize,
                *base_ptr.add(2) as usize,
                *base_ptr.add(3) as usize,
            ]
        }
    }

    #[inline(always)]
    fn get_points(&self) -> [Point3f; 4] {
        let [v0, v1, v2, v3] = self.get_vertex_indices();
        unsafe {
            [
                *self.mesh.p.add(v0),
                *self.mesh.p.add(v1),
                *self.mesh.p.add(v2),
                *self.mesh.p.add(v3),
            ]
        }
    }

    #[inline(always)]
    fn get_uvs(&self) -> Option<[Point2f; 4]> {
        if self.mesh.uv.is_null() {
            return None;
        }
        let [v0, v1, v2, v3] = self.get_vertex_indices();
        unsafe {
            Some([
                *self.mesh.uv.add(v0),
                *self.mesh.uv.add(v1),
                *self.mesh.uv.add(v2),
                *self.mesh.uv.add(v3),
            ])
        }
    }

    #[inline(always)]
    fn get_shading_normals(&self) -> Option<[Normal3f; 4]> {
        if self.mesh.n.is_null() {
            return None;
        }
        let [v0, v1, v2, v3] = self.get_vertex_indices();
        unsafe {
            Some([
                *self.mesh.n.add(v0),
                *self.mesh.n.add(v1),
                *self.mesh.n.add(v2),
                *self.mesh.n.add(v3),
            ])
        }
    }

    #[cfg(not(target_os = "cuda"))]
    pub fn new(mesh: Ptr<DeviceBilinearPatchMesh>, blp_index: i32) -> Self {
        let mut bp = BilinearPatchShape {
            mesh,
            blp_index,
            area: 0.,
            rectangle: false,
        };

        let [p00, p10, p01, p11] = bp.get_points();

        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 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,
        corners: &[Point3f; 4],
    ) -> Option<BilinearIntersection> {
        let &[p00, p01, p10, p11] = corners;
        let a = (p10 - p00).cross(p01 - p11).dot(ray.d);
        let c = (p00 - ray.o).cross(ray.d).dot(p01 - p00);
        let b = (p10 - ray.o).cross(ray.d).dot(p11 - 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()
                + p00.abs().max_component_value()
                + p10.abs().max_component_value()
                + p01.abs().max_component_value()
                + p11.abs().max_component_value());

        let hit1 = self.check_candidate(u1, ray, corners);
        let hit2 = if u1 != u2 {
            self.check_candidate(u2, ray, corners)
        } 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,
        corners: &[Point3f; 4],
    ) -> Option<IntersectionData> {
        if !(0.0..=1.0).contains(&u) {
            return None;
        }
        let &[p00, p01, p10, p11] = corners;
        let uo: Point3f = lerp(u, p00, p10);
        let ud: Point3f = Point3f::from(lerp(u, p01, 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,
        uv: Point2f,
        time: Float,
        wo: Vector3f,
    ) -> SurfaceInteraction {
        // Base geom and derivatives
        let corners = self.get_points();
        let [p00, p01, p10, p11] = corners;
        let p = lerp(uv[0], lerp(uv[1], p00, p01), lerp(uv[1], p10, p11));
        let mut dpdu = lerp(uv[1], p10, p11) - lerp(uv[1], p00, p01);
        let mut dpdv = lerp(uv[0], p01, p11) - lerp(uv[0], p00, p10);

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

        // Compute second moments
        let n = Normal3f::from(dpdu.cross(dpdv).normalize());
        let (mut dndu, mut dndv) = self.calculate_surface_curvature(&corners, &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 = p00.abs()
            + Vector3f::from(p01.abs())
            + Vector3f::from(p10.abs())
            + Vector3f::from(p11.abs());
        let p_error = gamma(6) * Vector3f::from(p_abs_sum);

        let flip_normal = self.mesh.reverse_orientation ^ self.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);

        self.apply_shading_normals(&mut isect, self.get_shading_normals(), uv, derivatives);
        isect
    }

    #[inline(always)]
    fn apply_texture_coordinates(
        &self,
        uv: Point2f,
        patch_uvs: Option<[Point2f; 4]>,
        dpdu: &mut Vector3f,
        dpdv: &mut Vector3f,
    ) -> (Point2f, Option<TextureDerivative>) {
        let Some(uvs) = patch_uvs else {
            return (uv, Some(TextureDerivative::default()));
        };
        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))
    }

    #[inline(always)]
    fn calculate_base_derivatives(
        &self,
        corners: &[Point3f; 4],
        uv: Point2f,
    ) -> (Point3f, Vector3f, Vector3f) {
        let (p00, p10, p01, p11) = (corners[0], corners[1], corners[2], corners[3]);

        let p = lerp(uv[0], lerp(uv[1], p00, p01), lerp(uv[1], p10, p11));
        let dpdu = lerp(uv[1], p10, p11) - lerp(uv[1], p00, p01);
        let dpdv = lerp(uv[0], p01, p11) - lerp(uv[0], p00, p10);
        (p, dpdu, dpdv)
    }

    #[inline(always)]
    fn calculate_surface_curvature(
        &self,
        corners: &[Point3f; 4],
        dpdu: &Vector3f,
        dpdv: &Vector3f,
        n: Normal3f,
    ) -> (Normal3f, Normal3f) {
        let (p00, p10, p01, p11) = (corners[0], corners[1], corners[2], corners[3]);
        let e = dpdu.dot(*dpdu);
        let f = dpdu.dot(*dpdv);
        let g = dpdv.dot(*dpdv);
        let d2pduv = (p00 - p01) + (p11 - 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,
        shading_normals: Option<[Normal3f; 4]>,
        uv: Point2f,
        derivatives: Option<TextureDerivative>,
    ) {
        let Some(normals) = shading_normals 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)?;
        ss.intr.get_common_mut().time = ctx.time;

        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, corners: &[Point3f; 4], u: Point2f) -> (Point2f, Float) {
        let (p00, p10, p01, p11) = (corners[0], corners[1], corners[2], corners[3]);
        if !self.mesh.image_distribution.is_null() {
            let image_distrib = self.mesh.image_distribution;
            let (uv, pdf, _) = image_distrib.sample(u);
            (uv, pdf)
        } else if !self.rectangle {
            let w = [
                (p10 - p00).cross(p01 - p00).norm(),
                (p10 - p00).cross(p11 - p10).norm(),
                (p01 - p00).cross(p11 - p01).norm(),
                (p11 - p10).cross(p11 - 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,
        corners: &[Point3f; 4],
        ctx: &ShapeSampleContext,
        u: Point2f,
        corner_dirs: &[Vector3f; 4],
    ) -> Option<ShapeSample> {
        let (p00, p10, p01, _p11) = (corners[0], corners[1], corners[2], corners[3]);
        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 = p10 - p00;
        let ev = p01 - p00;
        let (p, quad_pdf) = sample_spherical_rectangle(ctx.p(), p00, eu, ev, u);
        pdf *= quad_pdf?;

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

        let patch_uvs = self.get_uvs();
        let patch_normals = self.get_shading_normals();
        let n = self.compute_sampled_normal(patch_normals, &eu, &ev, uv);
        let st = patch_uvs.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: Interaction::Surface(intr),
            pdf,
        })
    }

    fn compute_sampled_normal(
        &self,
        patch_normals: Option<[Normal3f; 4]>,
        dpdu: &Vector3f,
        dpdv: &Vector3f,
        uv: Point2f,
    ) -> Normal3f {
        let mut n = Normal3f::from(dpdu.cross(*dpdv).normalize());

        if let Some(normals) = patch_normals {
            // 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);
        } else if self.mesh.reverse_orientation ^ self.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 [p00, p01, p10, p11] = self.get_points();
        let normals = self.get_shading_normals();
        if p00 == p10 || p10 == p11 || p11 == p01 || p01 == p00 {
            let dpdu = lerp(0.5, p10, p11) - lerp(0.5, p00, p01);
            let dpdv = lerp(0.5, p01, p11) - lerp(0.5, p00, p10);
            let mut n = Normal3f::from(dpdu.cross(dpdv).normalize());
            if let Some(normals) = normals {
                let interp_n = (normals[0] + normals[1] + normals[2] + normals[3]) / 4.;
                n = n.face_forward(interp_n);
            } else if self.mesh.reverse_orientation ^ self.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((p10 - p00).cross(p01 - p00).normalize());
        let mut n10 = Normal3f::from((p11 - p10).cross(p00 - p10).normalize());
        let mut n01 = Normal3f::from((p00 - p01).cross(p11 - p01).normalize());
        let mut n11 = Normal3f::from((p01 - p11).cross(p10 - p11).normalize());

        if let Some(normals) = normals {
            n00 = n00.face_forward(normals[0]);
            n10 = n10.face_forward(normals[1]);
            n01 = n01.face_forward(normals[2]);
            n11 = n11.face_forward(normals[3]);
        } else if self.mesh.reverse_orientation ^ self.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 [p00, p01, p10, p11] = self.get_points();
        Bounds3f::from_points(p00, p01).union(Bounds3f::from_points(p10, p11))
    }

    #[inline]
    fn intersect(&self, ray: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
        let t_max_val = t_max?;
        if let Some(bilinear_hit) =
            self.intersect_bilinear_patch(ray, t_max_val, &self.get_points())
        {
            let intr = self.interaction_from_intersection(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 corners = self.get_points();
        self.intersect_bilinear_patch(ray, t_max_val, &corners)
            .is_some()
    }

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

        // Compute surface normal for sampled bilinear patch (u, v)
        let patch_normals = self.get_shading_normals();
        let patch_uvs = self.get_uvs();
        let n = self.compute_sampled_normal(patch_normals, &dpdu, &dpdv, uv);
        let st = patch_uvs.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 = p00.abs()
            + Vector3f::from(p01.abs())
            + Vector3f::from(p10.abs())
            + Vector3f::from(p11.abs());
        let p_error = gamma(6) * Vector3f::from(p_abs_sum);
        let pi = Point3fi::new_with_error(p, p_error);
        Some(ShapeSample {
            intr: Interaction::Surface(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 corners = self.get_points();
        let [p00, p01, p10, p11] = corners;
        let v00 = (p00 - ctx.p()).normalize();
        let v10 = (p10 - ctx.p()).normalize();
        let v01 = (p01 - ctx.p()).normalize();
        let v11 = (p11 - ctx.p()).normalize();

        let use_area_sampling = self.rectangle
            || !self.mesh.image_distribution.is_null()
            || 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(&corners, ctx, u, &[v00, v10, v01, v11])
        }
    }

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

        let corners = self.get_points();
        let [p00, p01, p10, p11] = corners;
        let patch_uvs = self.get_uvs();
        let uv = if let Some(uvs) = patch_uvs {
            Point2f::invert_bilinear(si.common.uv, &uvs)
        } else {
            si.common.uv
        };

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

        let (_, dpdu, dpdv) = self.calculate_base_derivatives(&corners, 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 corners = self.get_points();
        let [p00, p01, p10, p11] = corners;

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

        let use_area_sampling = !self.rectangle
            || !self.mesh.image_distribution.is_null()
            || 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(),
                    p00,
                    p10 - p00,
                    p01 - p00,
                    isect.intr.p(),
                );
                pdf *= bilinear_pdf(u, &w);
            }
            pdf
        }
    }
}