Nori pa4 分布式和Whitted-style光线追踪

Xial 发布于 2023-01-30 96 次阅读


1. Area lights

让任何模型都能成为光源。

Step 1

使用蒙特卡洛积分均匀采样来计算光照。

需要修改 Mesh 类来实现。

在光源模型的表面均匀采样,每次采样需要返回以下信息:

  1. 在模型表面获得的采样点的坐标 p
  2. 模型表面在 p 点的法线 n
    • 如果 mesh 提供了每个顶点的法线,使用这些信息来计算法线
    • 否则通过向量叉乘来计算。
  3. 本次采样的概率密度。由于采用均匀采样,每次的概率密度应该均为整个模型总表面积分之一。

DiscretePDF 类实现了管理所有三角面概率密度的方法。

在三角形内均匀取点可以通过下列给定的参数来计算:

$$
\begin{pmatrix}
\alpha\\
\beta
\end{pmatrix}
= \begin{pmatrix}
1 - \sqrt{1 - \xi_1}\\
\xi_2\, \sqrt{1 - \xi_1}
\end{pmatrix}
$$

其中 $\xi_1$ 和 $\xi_2$ 是 $[0,1]$ 上的均匀分布。

则使用三个顶点的线性组合计算得到三角形内的均匀采样,使用的系数为 $(\alpha, \beta, 1-\alpha-\beta)$

Mesh 类中添加声明和定义:

/// nori/mesh.h

class Mesh {
   public:
    // ...
    MeshSampleResult sampleTriangleUniform(Sampler* sampler) const;

   protected:
    // ...
    DiscretePDF m_dis_pdf;
};
/// src/mesh.cpp

void Mesh::activate() {
    if (!m_bsdf) {
        // ...
    }
    m_dis_pdf = DiscretePDF(getTriangleCount());
    for (uint32_t i = 0; i < getTriangleCount(); i++) {
        auto area = surfaceArea(i);
        m_dis_pdf.append(area);
    }
    m_dis_pdf.normalize();  // 归一化
}

MeshSampleResult Mesh::sampleTriangleUniform(Sampler* sampler) const {
    MeshSampleResult res;
    uint32_t idx = m_dis_pdf.sample(sampler->next1D());
    auto samp = sampler->next2D();
    auto alpha = 1 - sqrt(1 - samp.x());
    auto beta = samp.y() * sqrt(1 - samp.x());
    Point3f v[3];
    for (int i = 0; i < 3; i++) {
        v[i] = m_V.col(m_F(i, idx));
    }
    res.p = v[0] * alpha + v[1] * beta + v[2] * (1 - alpha - beta);
    if (m_N.size()) {   // 如果 mesh 提供了所有顶点处的法线
        Normal3f n[3];
        for (int i = 0; i < 3; i++) {
            n[i] = m_N.col(m_F(i, idx));
        }
        res.n = (n[0] * alpha + n[1] * beta + n[2] * (1 - alpha - beta)).normalized();
    } else {
        auto e1 = v[1] - v[0];
        auto e2 = v[2] - v[0];
        res.n = e1.cross(e2).normalized();
    }
    res.pdf = m_dis_pdf.getNormalization();
    return res;
}

Step 2

比着 BSDF 接口抄一份 Emitter 接口。

仿照 BSDFQueryRecord 写的 EmitterQueryRecord 里面放什么暂时先不管。

struct EmitterQueryRecord {
};

/**
 * \brief Superclass of all emitters
 */
class Emitter : public NoriObject {
   public:
    virtual Color3f sample(const Mesh* mesh, EmitterQueryRecord& eRec, Sampler* sampler) const = 0;
    virtual Color3f eval(const EmitterQueryRecord& eRec) const = 0;
    virtual float pdf(const EmitterQueryRecord& eRec) const = 0;

    /**
     * \brief Return the type of object (i.e. Mesh/Emitter/etc.)
     * provided by this instance
     * */
    EClassType getClassType() const { return EEmitter; }
};

然后在 area.cpp 里实现这个接口,记得添加到 CMakeLists.txt

#include <nori/emitter.h>
#include <nori/mesh.h>

NORI_NAMESPACE_BEGIN

class AreaLight : public Emitter {
   public:
    AreaLight(const PropertyList& prop) {
        m_radiance = prop.getColor("radiance");
    }

    Color3f sample(const Mesh* mesh, EmitterQueryRecord& eRec, Sampler* sampler) const {
        auto mRes = mesh->sampleTriangleUniform(sampler);
        eRec.n = mRes.n;
        eRec.p = mRes.p;
        eRec.px = (eRec.p - eRec.x).normalized();
        auto p = pdf(mesh, eRec);
        if (!isnan(p) && !isinf(p) && p > 0) {
            return eval(eRec) / p;
        }
        return Color3f(0.0f);
    }

    Color3f eval(const EmitterQueryRecord& eRec) const {
        return (eRec.n.dot(eRec.px) < 0.0f) ? m_radiance : Color3f(0.0f);
    }

    float pdf(const Mesh* mesh, const EmitterQueryRecord& eRec) const {
        auto cosTheta = eRec.n.dot(-eRec.px);
        if (cosTheta > 0.0f) {
            auto p = mesh->getDiscreatePDF().getNormalization();
            auto px = eRec.p - eRec.x;
            return p * px.dot(px) / cosTheta;
        }
        return 0;
    }

    std::string toString() const {
        return tfm::format(
            "AreaLight[\n"
            "  m_radiance =
            "]",
            m_radiance.toString());
    }

   private:
    Color3f m_radiance;
};

NORI_REGISTER_CLASS(AreaLight, "area")
NORI_NAMESPACE_END
struct EmitterQueryRecord {
    Point3f x;
    Point3f p;
    Normal3f n;
    Vector3f px;

    EmitterQueryRecord(const Point3f& _x)
        : x(_x) {}

    EmitterQueryRecord(const Point3f& _x, const Point3f& _p, const Normal3f& _n)
        : x(_x), p(_p), n(_n) {
        px = (_p - _x).normalized();
    }
};

2. Distribution Ray Tracing

渲染方程长这个样子。

$$


L_r(\vx,\omega_r) = \int_{\mathcal{H}^2} f_r (\vx,\omega_i,\omega_r)\,L_i (\vx,\omega_i) \cos\theta_i\, \mathrm{d}\omega_i.
$$

其中 $f_r$ 是材质函数,已经由 bsdf 实现好了。同时我们实现了 $L_i$ 的计算。

但是,这个任务中我们只计算直接光照,这就意味着,除了射线碰巧击中光源,其他地方返回都是 0,这样做的效率极低,只有一部分射线会碰到光源,其他射线会产生大量噪点。所以我们需要让射线更容易击中光源。

更好的办法是,直接在光源上采样并检查可见性,而不是在物体表面采样。这意味着我们需要将渲染方程从半球上改写到光源上。

$$

L_r(\vx,\omega_r) = \int_{\mathcal{L}} f_r (\vx,\vx\to\vy,\omega_r)\,L_e (\vy,\vy\to\vx) \, \mathrm{d} \vy?? \qquad(\text{warning: this is not (yet) correct})
$$

其中 $\mathbf{x}\to\mathbf{y}$ 表示从 $\mathbf{x}$ 到 $\mathbf{y}$ 的归一化方向, $L_e(\mathbf{x},\omega)$ 表示在 $\mathbf{x}$ 点向 $\omega$ 方向发射的 radiance

但是这个积分算式不正确,因为我们将积分变量从立体角改为了位置,应该有一个匹配的因子来修正这一点,这里我们把它称为几何项

$$

G(\vx\leftrightarrow\vy) :=V(\vx\leftrightarrow\vy)\frac{
|\vn_\vx \cdot(\vx\to\vy)|\,\cdot\,
|\vn_\vy \cdot(\vy\to\vx)|}{|\vx-\vy|^2}
$$

第一项 $V(\mathbf{x}\leftrightarrow\mathbf{y})$ 是可见性函数,0 表示不可见,1 表示可见。分子包含两个点乘的绝对值,相当于把光源面积投影到着色点单位圆上的那块面积,分母是我们在使用点光源渲染时已经观察到的平方反比衰减。最终的渲染方程是:

$$

L_r(\vx,\omega_r) = \int_{\mathcal{L}} f_r (\vx,\vx\to\vy,\omega_r)\,G(\vx\leftrightarrow\vy)\,L_e (\vy,\vy\to\vx)\, \mathrm{d} \vy
$$

只需要按照如下步骤操作:

  1. src/whitted.cpp 创建一个新的积分器
  2. 这个积分器应当先像 simple 一样先找到与摄像头相交的最近的一个点 $\mathbf{x}$ ,返回的值 L_i 应该由如下式子得出:
    $$

    L_i(\vc,\omega_c)= L_e(\vr(\vc, \omega_c), -\omega_c) + L_r(\vr(\vc, \omega_c), -\omega_c).
    $$
    其中 $L_e$ 是该物体本身发出的光(如果该物体是光源),另一部分则是下一步采样得出的。
  3. 对于每个采样点 $\mathbf{x}$ 都应该通过随机在光源上采样得到点 $\mathbf{y}$ 并计算:
    $$

    f_r (\vx,\vx\to\vy,\omega_r)\,G(\vx\leftrightarrow\vy)\,L_e (\vy,\vy\to\vx)
    $$
    将上式除以光源上当前采样点的概率密度。

最终得到主要代码:

class WhittedIntegrator : public Integrator {
   public:
    WhittedIntegrator(const PropertyList& props) {}

    Color3f Li(const Scene* scene, Sampler* sampler, const Ray3f& ray) const {
        Intersection its;
        if (!scene->rayIntersect(ray, its)) {
            return Color3f(0.0f);
        }
        Color3f L_emitter(0.0f);
        if (its.mesh->isEmitter()) {
            EmitterQueryRecord eRec(ray.o, its.p, its.shFrame.n);
            L_emitter = its.mesh->getEmitter()->eval(eRec);
        }
        auto light = scene->getRandomEmitter(sampler->next1D());
        EmitterQueryRecord eRec(its.p);
        auto Le = light->getEmitter()->sample(light, eRec, sampler);
        auto cosTheta = Frame::cosTheta(its.shFrame.toLocal(eRec.px));
        if (cosTheta < 0) {
            return L_emitter;
        }
        if (scene->rayIntersect(Ray3f(eRec.x, eRec.px, Epsilon, (eRec.p - eRec.x).norm() - Epsilon))) {
            return L_emitter;
        }
        BSDFQueryRecord bRec(its.toLocal(-ray.d), its.toLocal(eRec.px), ESolidAngle);
        auto fr = its.mesh->getBSDF()->eval(bRec);
        return L_emitter + Le * fr * cosTheta * scene->getEmitters().size();
    }

    std::string toString() const {
        return "WhittedIntegrator[]";
    }
};
class AreaLight : public Emitter {
   public:
    AreaLight(const PropertyList& prop) {
        m_radiance = prop.getColor("radiance");
    }

    Color3f sample(const Mesh* mesh, EmitterQueryRecord& eRec, Sampler* sampler) const {
        auto mRes = mesh->sampleTriangleUniform(sampler);
        eRec.n = mRes.n;
        eRec.p = mRes.p;
        eRec.px = (eRec.p - eRec.x).normalized();
        auto p = pdf(mesh, eRec);
        if (!isnan(p) && !isinf(p) && p > 0) {
            return eval(eRec) / p;
        }
        return Color3f(0.0f);
    }

    Color3f eval(const EmitterQueryRecord& eRec) const {
        return (eRec.n.dot(eRec.px) < 0.0f) ? m_radiance : Color3f(0.0f);
    }

    float pdf(const Mesh* mesh, const EmitterQueryRecord& eRec) const {
        auto cosTheta = eRec.n.dot(-eRec.px);
        if (cosTheta > 0.0f) {
            auto p = mesh->getDiscreatePDF().getNormalization();
            auto px = eRec.p - eRec.x;
            return p * px.dot(px) / cosTheta;
        }
        return 0;
    }

    std::string toString() const {
        return tfm::format(
            "AreaLight[\n"
            "  m_radiance =
            "]",
            m_radiance.toString());
    }

   private:
    Color3f m_radiance;
};

注意:其实还需要对八叉树进行一些修改使其支持多个 mesh,同时需要在场景中实现一个随机取光源的函数。

3. Dielectrics

src/mirror.cpp 里已经实现好了一个基于 Dirac delta function 的镜面反射的 mirror 材质。

并且还在 src/common.cpp 里实现好了一个 fresnel 函数用来计算介电材料的非偏振菲涅尔反射系数。

只需要实现一个电介质(绝缘体?晶体?)材质放在 src/dielectric.cpp 里面就好了。

只要根据入射角的菲涅尔结果选择光线是反射还是折射就好了。

主要代码如下:

static Vector3f refract(const Vector3f& wi, const Vector3f& n, float eta) {
    float cosThetaI = wi.dot(n);
    if (cosThetaI < 0)
        eta = 1 / eta;
    float cosThetaTSqr = 1 - (1 - cosThetaI * cosThetaI) * eta * eta;
    if (cosThetaTSqr <= 0)
        return Vector3f(0);
    float sign = (cosThetaI >= 0) ? 1 : -1;
    return n * (-cosThetaI * eta + sign * sqrt(cosThetaTSqr)) + wi * eta;
}

class Dielectric : public BSDF {
   public:
    // ...
    Color3f sample(BSDFQueryRecord& bRec, const Point2f& sample) const {
        bRec.measure = EDiscrete;
        float kr = fresnel(Frame::cosTheta(bRec.wi), m_extIOR, m_intIOR);
        if (sample.x() < kr) {
            bRec.wo = Vector3f(-bRec.wi.x(), -bRec.wi.y(), bRec.wi.z());
            bRec.eta = 1;
        } else {
            auto cosThetaI = Frame::cosTheta(bRec.wi);
            Vector3f n(0, 0, (cosThetaI < 0) ? -1 : 1);
            float factor = (cosThetaI < 0) ? (m_extIOR / m_intIOR) : (m_intIOR / m_extIOR);
            bRec.wo = refract(-bRec.wi, n, factor);
            bRec.eta = m_intIOR / m_extIOR;
        }
        return Color3f(1);
    }
    // ...
}

4. Whitted-style ray tracing

最后实现一个 Whitted-style 光线追踪就好了,只要根据采样点的材质选择光线的走向。

因为反射跟折射会让光线继续前进,直到射在漫反射(diffuse)材质上。使用如下递归方法来估计最终的光强:

$$

L_i(\vc, \omega_c) =
\begin{cases}
\frac{1}{0.95}c L_i(\vx, \omega_r),&\text{if $\xi < 0.95$}\\
0,&\text{otherwise}
\end{cases}
$$

其中 $\xi$ 是一个随机数。

WhittedIntegrator 修改成下面这样就好了:

Color3f Li(const Scene* scene, Sampler* sampler, const Ray3f& ray) const {
    Intersection its;
    if (!scene->rayIntersect(ray, its)) {
        return Color3f(0.0f);
    }
    if (its.mesh->getBSDF()->isDiffuse()) {
        Color3f L_emitter(0.0f);
        // ...
        return L_emitter + Le * fr * cosTheta * scene->getEmitters().size();
    } else {
        BSDFQueryRecord bRec(its.toLocal(-ray.d));
        Color3f refColor = its.mesh->getBSDF()->sample(bRec, sampler->next2D());
        if (sampler->next1D() < 0.95 && refColor.x() > 0) {
            return Li(scene, sampler, Ray3f(its.p, its.toWorld(bRec.wo))) / 0.95 * refColor;
        } else {
            return Color3f(0);
        }
    }
}
最后更新于 2023-02-05