#include <math.h>
#include <core/ray.h>
#include <core/logger.h>
#include <core/transform.h>
#include <core/sampler.h>
#include <render/camera.h>
#include <render/scene.h>
#include <render/shape.h>
#include <algorithm>

Camera::Camera(const Properties &props)
    : width(props.get<int>("width")),
      height(props.get<int>("height")),
      m_fov(props.get<float>("fov")),
      cam_to_world(props.get<Matrix4x4>("to_world")),
      rect{props.get<int>("offset_x", 0),
           props.get<int>("offset_y", 0),
           props.get<int>("crop_width", 0),
           props.get<int>("crop_height", 0)}
{
    initRfilter(props.get<Properties>("rfilter", {}));
    configure();
}

Camera::Camera(int width, int height, const Matrix4x4 &cam_to_world, const Matrix3x3 &cam_to_ndc, int med_id, int type)
    : width(width), height(height), cam_to_world(cam_to_world), cam_to_ndc(cam_to_ndc), med_id(med_id)
{
    world_to_cam = cam_to_world.inverse();
    ndc_to_cam = cam_to_ndc.inverse();

    cpos = this->cam_to_world.block<3, 1>(0, 3);
    cframe.s = this->cam_to_world.block<3, 1>(0, 0);
    cframe.t = this->cam_to_world.block<3, 1>(0, 1);
    cframe.n = this->cam_to_world.block<3, 1>(0, 2);
    initRfilter(type);
    m_state = ESConfigured;
}

Camera::Camera(int width, int height, Float fov, const Matrix4x4 &to_world, int type)
    : width(width), height(height), m_fov(fov), cam_to_world(to_world)
{
    initRfilter(type);
    configure();
}

Camera::Camera(int width, int height, Float fov, const Vector &origin, const Vector &target, const Vector &up, int type)
    : width(width), height(height), m_fov(fov), cam_to_world(psdr::transform::look_at(origin, target, up))
{
    initRfilter(type);
    configure();
}

Camera::Camera(int width, int height, Float fov, const Matrix4x4 &cam_to_world, ReconstructionFilter *filter)
    : width(width), height(height), m_fov(fov), cam_to_world(cam_to_world), rfilter(filter), m_state(ESLoaded)
{
    sigX = Eigen::MatrixXd::Ones(height, width) * .5;
    sigY = Eigen::MatrixXd::Ones(height, width) * .5;
    configure();
}

Camera::Camera(const Camera &camera)
{
    copyFrom(camera);
}

Camera &Camera::operator=(const Camera &camera)
{
    copyFrom(camera);
    return *this;
}

inline void Camera::copyFrom(const Camera &camera)
{
    width = camera.width;
    height = camera.height;
    m_fov = camera.m_fov;
    cam_to_world = camera.cam_to_world;
    cam_to_ndc = camera.cam_to_ndc;
    med_id = camera.med_id;
    cam_to_world = camera.cam_to_world;
    world_to_cam = camera.world_to_cam;
    ndc_to_cam = camera.ndc_to_cam;
    cam_to_ndc = camera.cam_to_ndc;
    cpos = camera.cpos;
    cframe = camera.cframe;
    sigX = camera.sigX;
    sigY = camera.sigY;
    rfilter = camera.rfilter->clone();
    rect = camera.rect;
}
inline void Camera::initRfilter(const Properties &props)
{
    if (rfilter != nullptr)
        delete rfilter;
    std::string type = props.get<std::string>("type");
    if (type == "box")
    {
        rfilter = new BoxFilter();
    }
    else if (type == "tent")
    {
        rfilter = new TentFilter(props.get<float>("padding", 0.5));
    }
    else if (type == "aniso")
    {
        rfilter = new AnisotropicGaussianFilter();
        sigX = Eigen::MatrixXd::Ones(height, width) * .5;
        sigY = Eigen::MatrixXd::Ones(height, width) * .5;
    }
    else
    {
        assert(false);
    }
}

inline void Camera::initRfilter(int type)
{
    if (rfilter != nullptr)
        delete rfilter;
    switch (type)
    {
    case 0:
        rfilter = new TentFilter(0.5);
        break;
    case 1:
        rfilter = new BoxFilter();
        break;
    case 2:
        rfilter = new AnisotropicGaussianFilter();
        sigX = Eigen::MatrixXd::Ones(height, width) * .5;
        sigY = Eigen::MatrixXd::Ones(height, width) * .5;
        break;
    default:
        assert(false);
    }
}

void Camera::merge(const Camera &camera)
{
    cam_to_world += camera.cam_to_world;
    world_to_cam += camera.world_to_cam;
    ndc_to_cam += camera.ndc_to_cam;
    cam_to_ndc += camera.cam_to_ndc;
    cpos += camera.cpos;
    cframe.merge(camera.cframe);
}

void Camera::setZero()
{
    cpos.setZero();
    cam_to_ndc.setZero();
    ndc_to_cam.setZero();
    cam_to_world.setZero();
    world_to_cam.setZero();
    cframe.t.setZero();
    cframe.s.setZero();
    cframe.n.setZero();
}

/*  configure the camera given
int width, int height, Float fov,
const Matrix4x4 &cam_to_world,
Float near, Float far */
void Camera::configure()
{
    world_to_cam = cam_to_world.inverse();
    Float fov_factor = 1.0 / tan(m_fov * 0.5 * M_PI / 180.0);
    cam_to_ndc = Eigen::Scaling(fov_factor, fov_factor, 1.0);
    ndc_to_cam = cam_to_ndc.inverse();
    cpos = this->cam_to_world.block<3, 1>(0, 3);
    cframe.s = this->cam_to_world.block<3, 1>(0, 0);
    cframe.t = this->cam_to_world.block<3, 1>(0, 1);
    cframe.n = this->cam_to_world.block<3, 1>(0, 2);
    m_state = ESConfigured;
}

inline Float Camera::evalFilter(int ix, int iy, Float x, Float y) const
{
    assert(rfilter != nullptr);
    if (dynamic_cast<const AnisotropicGaussianFilter *>(rfilter))
    {
        Float sigs[2] = {sigX(iy, ix), sigY(iy, ix)};
        return rfilter->eval(Vector2(x, y), sigs);
    }
    else
        return rfilter->eval(Vector2(x, y));
}

inline Float Camera::_eval(int ix, int iy, const Vector &p, const Vector *n, bool eval_filter) const
{
    Vector dir = p - cpos;
    Float dist = dir.norm();
    dir /= dist;
    Float cosy = cframe.toLocal(dir).z();
    //! must use geonormal. its.wi.z() is shading normal.
    Float cosx = (n != nullptr ? n->dot(-dir) : 1.);
    Float G = cosx / (cosy * cosy * cosy * dist * dist);
    Vector2 uv = getPixelUV(dir);
    Vector2 offset = uv - Vector2(ix, iy);
    Float f = evalFilter(ix, iy, offset.x(), offset.y());
    Float baseVal = f * G;

    if (eval_filter)
    {
        Float pdf = detach(baseVal);
        return (pdf > FilterEpsilon ? baseVal / pdf : 1.);
    }
    else
    {
        auto fov_factor = cam_to_ndc(0, 0);
        auto aspect_ratio = Float(width) / Float(height);
        Float inv_pixel_area = 0.25 * fov_factor * fov_factor * aspect_ratio *
                               static_cast<Float>(width * height);
        return baseVal * inv_pixel_area;
    }
}

Float Camera::evalFilter(int ix, int iy, const Intersection &its) const
{
    return _eval(ix, iy, its.p, &its.geoFrame.n, true);
}

Float Camera::evalFilter(int ix, int iy, const Vector &p, const Vector &n) const
{
    return _eval(ix, iy, p, &n, true);
}

Float Camera::evalFilter(int ix, int iy, const Vector &p) const
{
    return _eval(ix, iy, p, nullptr, true);
}

Float Camera::eval(int ix, int iy, const Vector &p, const Vector &n) const
{
    return _eval(ix, iy, p, &n, false);
}

Float Camera::eval(int ix, int iy, const Vector &p) const
{
    return _eval(ix, iy, p, nullptr, false);
}

Vector2 Camera::getPixelUV(const Vector &dir) const
{
    Vector local = cframe.toLocal(dir);
    Vector ndc = cam_to_ndc * local;
    Float t = 1.0 / local.z();
    ndc *= t;
    auto aspect_ratio = Float(width) / Float(height);
    Vector2 ret;
    ret.x() = ((ndc.x() / 2.0 + 0.5)) * width;
    ret.y() = (ndc.y() * aspect_ratio / -2.0 + 0.5) * height;
    return ret;
}

Vector2 Camera::getPixelUV(const Ray &ray) const
{
    assert((ray.org - cpos).norm() < Epsilon);
    return getPixelUV(ray.dir);
}

inline void Camera::_samplePrimaryRayFromFilter(int ix, int iy, const Vector2 &rnd2, Float &x, Float &y) const
{
    assert(rfilter != nullptr);
    Vector2 p;
    if (dynamic_cast<const AnisotropicGaussianFilter *>(rfilter))
    {
        assert(sigX.rows() == height && sigX.cols() == width && sigY.rows() == height && sigY.cols() == width);
        assert(ix >= 0 && ix < width && iy >= 0 && iy < height);

        Float sigs[2] = {sigX(iy, ix), sigY(iy, ix)};
        assert(sigs[0] > Epsilon && sigs[1] > Epsilon);

        rfilter->sample(rnd2, p, sigs);
    }
    else
        rfilter->sample(rnd2, p);
    x = p[0];
    y = p[1];
}

void Camera::samplePrimaryRayFromFilter(int ix, int iy, const Vector2 &rnd2, Ray &ray) const
{
    Float x, y;
    _samplePrimaryRayFromFilter(ix, iy, rnd2, x, y);
    ray = samplePrimaryRay(ix + x, iy + y);
}

void Camera::samplePrimaryRayFromFilter(int ix, int iy, const Vector2 &rnd2, Ray &ray, Ray &dual) const
{
    Float x, y;
    _samplePrimaryRayFromFilter(ix, iy, rnd2, x, y);
    ray = samplePrimaryRay(ix + x, iy + y);

    Float x_dual = 1.0 - x,
          y_dual = 1.0 - y;
    dual = samplePrimaryRay(ix + x_dual, iy + y_dual);
}

void Camera::samplePrimaryRayFromFilter(int ix, int iy, const Vector2 &rnd2, Ray *rays) const
{
    Float x, y;
    _samplePrimaryRayFromFilter(ix, iy, rnd2, x, y);

    rays[0] = samplePrimaryRay(ix + x, iy + y);
    rays[1] = samplePrimaryRay(ix + 1.0 - x, iy + y);
    rays[2] = samplePrimaryRay(ix + x, iy + 1.0 - y);
    rays[3] = samplePrimaryRay(ix + 1.0 - x, iy + 1.0 - y);
}

Ray Camera::sampleDualRay(const Array2i &pixel_idx, Ray &ray) const
{
    Vector2 uv = getPixelUV(ray);
    uv = uv - pixel_idx.matrix().cast<Float>();
    assert(uv.x() >= -0.5 && uv.x() <= 1.5);
    assert(uv.y() >= -0.5 && uv.y() <= 1.5);
    Float x_dual = 1.0 - uv.x();
    Float y_dual = 1.0 - uv.y();
    return samplePrimaryRay(pixel_idx.x() + x_dual, pixel_idx.y() + y_dual);
}

Ray Camera::sampleDualRay(const Vector2i &pixel_idx, const Vector2 &pixel_uv) const
{
    int pixel_x = pixel_idx.x();
    int pixel_y = pixel_idx.y();
    Float dual_x = 1.0 - (pixel_uv.x() - pixel_x);
    Float dual_y = 1.0 - (pixel_uv.y() - pixel_y);
    assert(dual_x >= -a && dual_x <= a + 1.0);
    assert(dual_y >= -a && dual_y <= a + 1.0);
    return samplePrimaryRay(pixel_x + dual_x, pixel_y + dual_y);
}

Ray Camera::samplePrimaryRay(int pixel_x, int pixel_y,
                             const Vector2 &rnd2) const
{
    Vector2 screen_pos = Vector2((pixel_x + rnd2.x()) / Float(width),
                                 (pixel_y + rnd2.y()) / Float(height));
    auto aspect_ratio = Float(width) / Float(height);
    auto ndc = Vector3((screen_pos(0) - 0.5f) * 2.f,
                       (screen_pos(1) - 0.5f) * (-2.f) / aspect_ratio, Float(1));
    Vector3 dir = ndc_to_cam * ndc;
    dir.normalize();
    return Ray{xfm_point(cam_to_world, Vector3::Zero()),
               xfm_vector(cam_to_world, dir)};
}

Ray Camera::samplePrimaryRay(Float x, Float y) const
{
    x /= width;
    y /= height;
    auto aspect_ratio = Float(width) / Float(height);
    auto ndc =
        Vector3((x - 0.5f) * 2.f, (y - 0.5f) * (-2.f) / aspect_ratio, 1.0f);
    Vector3 dir = ndc_to_cam * ndc;
    dir.normalize();
    return Ray{xfm_point(cam_to_world, Vector3::Zero()),
               xfm_vector(cam_to_world, dir)};
}

std::tuple<Ray, Vector, Vector, Vector> Camera::samplePrimaryBoundaryRay(const Array2i &pixel_idx, Float rnd) const
{
    rnd *= 4.;
    Vector2 offset;
    Vector dir_visible;
    Vector dir_edge;
    if (rnd < 1)
    {
        offset = Vector2(0, rnd);
        dir_visible = Vector3(1, 0, 0);
        dir_edge = Vector3(0, 1, 0);
    }
    else if (rnd < 2)
    {
        offset = Vector2(rnd - 1, 1);
        dir_visible = Vector3(0, -1, 0);
        dir_edge = Vector3(1, 0, 0);
    }
    else if (rnd < 3)
    {
        offset = Vector2(1, 3 - rnd);
        dir_visible = Vector3(-1, 0, 0);
        dir_edge = Vector3(0, 1, 0);
    }
    else
    {
        offset = Vector2(4 - rnd, 0);
        dir_visible = Vector3(0, 1, 0);
        dir_edge = Vector3(1, 0, 0);
    }

    Vector2 screen_pos = Vector2((pixel_idx.x() + offset.x()) / Float(width),
                                 (pixel_idx.y() + offset.y()) / Float(height));
    auto aspect_ratio = Float(width) / Float(height);
    auto ndc = Vector3((screen_pos(0) - 0.5f) * 2.f,
                       (screen_pos(1) - 0.5f) * (-2.f) / aspect_ratio, Float(1));
    Vector3 dir_local = ndc_to_cam * ndc;
    dir_visible = ndc_to_cam * Vector3(dir_visible.x(), -dir_visible.y(), dir_visible.z());
    dir_edge = ndc_to_cam * Vector3(dir_edge.x(), -dir_edge.y(), dir_edge.z());
    dir_visible.normalize();
    dir_local.normalize();
    dir_edge.normalize();
    Ray ray{xfm_point(cam_to_world, Vector3::Zero()),
            xfm_vector(cam_to_world, dir_local)};
    return {ray, dir_visible, dir_local, dir_edge};
}

Float Camera::sampleDirect(const Vector &p, Vector2 &pixel_uv, Vector &dir) const
{
    Vector refP = xfm_point(world_to_cam, p);
    // if (refP.z() < clip_near) return 0.0;
    auto fov_factor = cam_to_ndc(0, 0);
    auto aspect_ratio = Float(width) / Float(height);
    //! might problematic. area = 2 x 2 x ndc_to_cam x ndc_to_cam * aspect_ratio
    Float inv_area = 0.25 * fov_factor * fov_factor * aspect_ratio;

    int xmin = 0, xmax = width;
    int ymin = 0, ymax = height;

    if (rect.isValid())
    {
        inv_area *= (Float)width / (Float)rect.crop_width * (Float)height /
                    rect.crop_height;
        xmin = rect.offset_x;
        xmax = rect.offset_x + rect.crop_width;
        ymin = rect.offset_y;
        ymax = rect.offset_y + rect.crop_height;
    }

    Vector pos_camera = cam_to_ndc * refP;
    pos_camera.x() /= pos_camera.z();
    pos_camera.y() /= pos_camera.z();
    Vector2 screen_pos =
        Vector2((pos_camera.x() * .5 + .5) * width,
                (-pos_camera.y() * .5 * aspect_ratio + .5) * height);
    if (screen_pos.x() >= xmin && screen_pos.x() <= xmax &&
        screen_pos.y() >= ymin && screen_pos.y() <= ymax)
    {
        pixel_uv.x() = screen_pos.x() - xmin;
        pixel_uv.y() = screen_pos.y() - ymin;
        Float dist = refP.norm(), inv_dist = 1. / dist;
        refP *= inv_dist;
        Float inv_cosTheta = 1. / refP.z();
        dir = (cpos - p) * inv_dist;
        Float inv_pixel_area = inv_area * (xmax - xmin) * (ymax - ymin);
        return inv_dist * inv_dist * inv_cosTheta * inv_cosTheta * inv_cosTheta *
               inv_pixel_area;
    }
    else
    {
        return 0.0;
    }
}

void Camera::sampleDirect(const Vector &p, Matrix2x4 &pixel_uvs, Array4 &weights, Vector &dir) const
{
    // assert(dynamic_cast<const TentFilter *>(rfilter));
    weights = Array4(0.0);
    Vector refP = xfm_point(world_to_cam, p);
    // if (refP.z() < clip_near) return;
    auto fov_factor = cam_to_ndc(0, 0);
    auto aspect_ratio = Float(width) / Float(height);
    Float inv_area = 0.25 * fov_factor * fov_factor * aspect_ratio;

    int xmin = 0, xmax = width;
    int ymin = 0, ymax = height;

    if (rect.isValid())
    {
        inv_area *= (Float)width / (Float)rect.crop_width * (Float)height /
                    rect.crop_height;
        xmin = rect.offset_x;
        xmax = rect.offset_x + rect.crop_width;
        ymin = rect.offset_y;
        ymax = rect.offset_y + rect.crop_height;
    }

    Vector pos_camera = cam_to_ndc * refP;
    pos_camera.x() /= pos_camera.z();
    pos_camera.y() /= pos_camera.z();
    Vector2 screen_pos =
        Vector2((pos_camera.x() * .5 + .5) * width,
                (-pos_camera.y() * .5 * aspect_ratio + .5) * height);
    if (screen_pos.x() >= xmin && screen_pos.x() <= xmax &&
        screen_pos.y() >= ymin && screen_pos.y() <= ymax)
    {
        Vector2 pixel_uv;
        pixel_uv.x() = screen_pos.x() - xmin;
        pixel_uv.y() = screen_pos.y() - ymin;
        Float dist = refP.norm(), inv_dist = 1. / dist;
        refP *= inv_dist;
        Float inv_cosTheta = 1. / refP.z();
        dir = (cpos - p) * inv_dist;
        Float inv_pixel_area = inv_area * (xmax - xmin) * (ymax - ymin);
        Float val0 = inv_dist * inv_dist * inv_cosTheta * inv_cosTheta *
                     inv_cosTheta * inv_pixel_area;

        // compute all involved neighboring pixels
        Vector2i pixel_index((int)pixel_uv.x(), (int)pixel_uv.y());
        Float x_offset = pixel_uv.x() - pixel_index.x();
        Float y_offset = pixel_uv.y() - pixel_index.y();
        int offsetX_min = (x_offset < a && pixel_index.x() > 0) ? -1 : 0;
        int offsetX_max =
            (x_offset > 1.0 - a && pixel_index.x() < xmax - 1) ? 1 : 0;
        int offsetY_min = (y_offset < a && pixel_index.y() > 0) ? -1 : 0;
        int offsetY_max =
            (y_offset > 1.0 - a && pixel_index.y() < ymax - 1) ? 1 : 0;

        int index = 0;
        for (int ox = offsetX_min; ox <= offsetX_max; ox++)
            for (int oy = offsetY_min; oy <= offsetY_max; oy++)
            {
                pixel_uvs(0, index) = pixel_uv.x() + ox;
                pixel_uvs(1, index) = pixel_uv.y() + oy;

                weights(index) = val0 * rfilter->eval(Vector2(x_offset - ox, y_offset - oy));
                index++;
            }
        return;
    }
    else
    {
        return;
    }
}

bool Camera::sampleDirect(const Vector &p, CameraDirectSamplingRecord &cRec) const
{
    Vector refP = xfm_point(world_to_cam, p);
    Float fov_factor = cam_to_ndc(0, 0);
    Float aspect_ratio = static_cast<Float>(width) / static_cast<Float>(height);
    Float inv_pixel_area = .25 * fov_factor * fov_factor * aspect_ratio
                           * static_cast<Float>(width * height);

    Vector pos_camera = cam_to_ndc * refP;
    pos_camera.x() /= pos_camera.z();
    pos_camera.y() /= pos_camera.z();

    Vector2 &pixel_uv = cRec.pixel_uv;
    pixel_uv.x() = (pos_camera.x() * .5 + .5) * static_cast<Float>(width);
    pixel_uv.y() = (-pos_camera.y() * .5 * aspect_ratio + .5) * static_cast<Float>(height);

    Vector2i &idx2 = cRec.pixel_idx;
    idx2 = pixel_uv.array().floor().matrix().template cast<int>();

    int padding = -1;
    if (dynamic_cast<const BoxFilter *>(rfilter))
        padding = 0;
    else if (dynamic_cast<const TentFilter *>(rfilter))
        padding = 1;

    auto [xmin, xmax, ymin, ymax] = getPixelAABB();
    if (padding < 0 || (idx2.x() >= xmin - padding && idx2.x() < xmax + padding &&
                        idx2.y() >= ymin - padding && idx2.y() < ymax + padding))
    {
        Float dist = refP.norm(), inv_dist = 1. / dist;
        refP *= inv_dist;
        Float inv_cosTheta = 1. / refP.z();
        cRec.dir = (cpos - p) * inv_dist;

        cRec.baseVal = inv_dist * inv_cosTheta;
        cRec.baseVal *= cRec.baseVal * inv_cosTheta * inv_pixel_area;
        return true;
    }
    return false;
}

std::pair<int, Float> Camera::sampleDirectPixel(const CameraDirectSamplingRecord &cRec, Float rnd) const
{
    auto [xmin, xmax, ymin, ymax] = getPixelAABB();
    const Vector2i &idx2 = cRec.pixel_idx;
    if (dynamic_cast<const BoxFilter *>(rfilter))
    {
        if (idx2.x() >= xmin && idx2.x() < xmax && idx2.y() >= ymin && idx2.y() < ymax)
            return {getPixelIndex1D(idx2), cRec.baseVal};
    }
    else if (dynamic_cast<const TentFilter *>(rfilter))
    {
        if (idx2.x() >= xmin - 1 && idx2.x() <= xmax && idx2.y() >= ymin - 1 && idx2.y() <= ymax)
        {
            Vector2 offset = cRec.pixel_uv - idx2.template cast<Float>();            
            const Float tent_a = dynamic_cast<const TentFilter *>(rfilter)->m_padding;

            int offsetX_min = offset.x() < tent_a       ? -1 : 0;
            int offsetX_max = offset.x() > 1.0 - tent_a ?  1 : 0;
            int offsetY_min = offset.y() < tent_a       ? -1 : 0;
            int offsetY_max = offset.y() > 1.0 - tent_a ?  1 : 0;

            Vector2i idx_offsets[4];
            Vector4 weights, cdf;

            int total = 0;
            for (int ox = offsetX_min; ox <= offsetX_max; ox++)
                for (int oy = offsetY_min; oy <= offsetY_max; oy++)
                {
                    Vector2i &cur_offset = idx_offsets[total];
                    cur_offset = Vector2i(ox, oy);
                    Vector2i cur_idx = idx2 + cur_offset;
                    if (cur_idx.x() >= xmin && cur_idx.x() < xmax && cur_idx.y() >= ymin && cur_idx.y() < ymax)
                    {
                        assert(total < 4);
                        weights[total] = rfilter->eval(offset - cur_offset.template cast<Float>());
                        cdf[total] = total ? cdf[total - 1] + weights[total] : weights[total];
                        total++;
                    }
                }

            if (total && cdf[total - 1] > Epsilon)
            {
                cdf /= static_cast<Float>(cdf[total - 1]);
                int i = static_cast<int>(std::lower_bound(&cdf[0], &cdf[0] + total, rnd) - &cdf[0]);
                i = std::clamp(i, 0, total - 1);
                Float pdf = i ? cdf[i] - cdf[i - 1] : cdf[i];
                return {getPixelIndex1D(idx2 + idx_offsets[i]), cRec.baseVal * weights[i] / pdf};
            }
        }
    }
    else
        PSDR_ASSERT_MSG(false, "Unsupported rfilter");

    return {-1, 0.};
}

void Camera::accumulateDirect(const CameraDirectSamplingRecord &cRec, const Spectrum &value, RndSampler *sampler, Spectrum *buffer) const
{
    auto [xmin, xmax, ymin, ymax] = getPixelAABB();
    const Vector2i &idx2 = cRec.pixel_idx;
    if (dynamic_cast<const BoxFilter *>(rfilter))
    {
        if (idx2.x() >= xmin && idx2.x() < xmax && idx2.y() >= ymin && idx2.y() < ymax) {
            int idx = getPixelIndex1D(idx2);
            buffer[idx] += value * cRec.baseVal;
        }
    }
    else if (dynamic_cast<const TentFilter *>(rfilter))
    {
        if (idx2.x() >= xmin - 1 && idx2.x() <= xmax && idx2.y() >= ymin - 1 && idx2.y() <= ymax)
        {
            Vector2 offset = cRec.pixel_uv - idx2.template cast<Float>();            
            const Float tent_a = dynamic_cast<const TentFilter *>(rfilter)->m_padding;

            int offsetX_min = offset.x() < tent_a       ? -1 : 0;
            int offsetX_max = offset.x() > 1.0 - tent_a ?  1 : 0;
            int offsetY_min = offset.y() < tent_a       ? -1 : 0;
            int offsetY_max = offset.y() > 1.0 - tent_a ?  1 : 0;

            for (int ox = offsetX_min; ox <= offsetX_max; ox++)
                for (int oy = offsetY_min; oy <= offsetY_max; oy++)
                {
                    Vector2i cur_offset = Vector2i(ox, oy), cur_idx = idx2 + cur_offset;
                    if (cur_idx.x() >= xmin && cur_idx.x() < xmax && cur_idx.y() >= ymin && cur_idx.y() < ymax)
                    {
                        int idx = getPixelIndex1D(cur_idx);
                        buffer[idx] += value * (cRec.baseVal * rfilter->eval(offset - cur_offset.template cast<Float>()));
                    }
                }
        }
    }
    else
        PSDR_ASSERT_MSG(false, "Unsupported rfilter");
}

Float Camera::sampleDirect(const Vector &p, Float rnd, int &pixel_idx, Vector &dir) const
{
    Float value = 0.;
    CameraDirectSamplingRecord cRec;
    if (sampleDirect(p, cRec) && cRec.baseVal > Epsilon) {
        dir = cRec.dir;
        std::tie(pixel_idx, value) = sampleDirectPixel(cRec, rnd);
    }
    return value;
}

Float Camera::geometric(const Vector &p, const Vector &n) const
{
    Vector dir = p - cpos;
    Float dist = dir.norm();
    dir /= dist;
    Float cosy = cframe.toLocal(dir).z();
    Float cosx = n.dot(-dir);
    return abs(cosx) / (cosy * cosy * cosy * dist * dist);
}

Float Camera::geometric(const Vector &p) const
{
    Vector dir = p - cpos;
    Float dist = dir.norm();
    dir /= dist;
    Float cosy = cframe.toLocal(dir).z();
    return 1. / (cosy * cosy * cosy * dist * dist);
}
