Nori pa3 蒙特卡洛采样和环境光遮蔽

Xial 发布于 2023-01-21 39 次阅读

Part 1: Monte Carlo Sampling

蒙特卡洛采样

前置知识

蒙特卡洛积分

求解积分 $I=\displaystyle \int_\Omega{f(x)dx}$
构造概率密度函数 $p(x)$

$$ I=\displaystyle \int_{\Omega}{\frac{f(x)}{p(x)}p(x)dx} = E(\frac{f(x)}{p(x)})\approx \frac{1}{N}\sum_{i=1}^N{\frac{f(x_i)}{p(x_i)}} $$

问题转化为求解 $\frac{f(x)}{p(x)}$ 的期望，通过在定义域上采样可以近似计算
当 $N \rightarrow \infty $ 的时候，可以得到积分的准确解 (无偏性)

该方法非常适合计算机实现

举例：$f(x)=3x^2$, 计算其在区间 $[a, b] (a=1,b=3)$ 的积分值

蒙特卡洛积分求解：

使用区间 $[a, b]$ 的均匀分布进行随机采样

$$ I = \frac{1}{N}\sum_{i=1}^N{\frac{f(x_i)}{p(x_i)}}\ ,\quad p(x)=\frac{1}{b-a} $$

  • 若采样样本是 {2} , $I=\frac{b-a}{N}\sum_{i=1}^N{\frac{f(x_i)}{p(x_i)}} = 24$;
  • 若采样样本是 {1,2,3} , $I=\frac{b-a}{N}\sum_{i=1}^N{\frac{f(x_i)}{p(x_i)}} = 28$;
  • 若采样样本是 {1.0, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3.0} , $I=\frac{b-a}{N}\sum_{i=1}^N{\frac{f(x_i)}{p(x_i)}} = 26.5$;
  • ...

收敛速度取决于被采样函数的采样方差

方差缩减：

  • 重要性采样
  • 分层采样
  • 拟蒙特卡洛采样
  • ......

重要性采样

概率密度函数需要满足

$$ \displaystyle \int_{\Omega}{p(x)dx} = 1 $$

求解右边函数 $f(x)$ 在定义域上的积分

如何设计概率密度函数在采样次数尽可能少的情况下得到准确结果？

$$ \frac{1}{N}\sum_N{\frac{f(x_i)}{p(x_i)}} $$

$\Rightarrow$ 采样得到的若干样本方差尽可能小

$$ \Rightarrow p(x) \propto f(x) $$

最理想情况

$$ p(x) = \frac{f(x)}{\int_{\Omega}{f(x)dx}}$$

此时方差为 0 ！

逆变换采样方法

已知概率密度函数 $p(x)$，如何按照概率密度函数进行采样？

求累积分布函数 CDF：$P(x)=\displaystyle \int_0^x{p(t)dt} $

对 CDF 求逆 $ P^{-1}(x) $

给定 $[0, 1]$ 均匀分布随机变量 $\epsilon$，生成符合 $p(x)$ 分布的随机变量 $ X=P^{-1}(\epsilon) $

例子：
$$
\begin{align}
p(x) &=\frac{8x}{\pi^2} \
P(x) &= \displaystyle \int_0^x{p(x)dx} = \frac{4x^2}{\pi^2} \
P^{-1}(t) &= \frac{\pi}{2}\sqrt{t} \
X_i &= \frac{\pi}{2}\sqrt{\xi}
\end{align}
$$

实操

Part 1 的任务是把 $[0,1]\times[0,1]$ 区域的均匀采样转化成不同区域的均匀采样。

其中 squareTo...Pdf 函数的作用是检验分布是否正确。

因此如果两个函数一起错了而且匹配是不会报错的。

Tent

帐篷分布？让我看看谁搭帐篷了

作业要求上甚至把解法都写明白了。

概率密度分布长这样：

$$
p(x,y) = p_1(x)p_1(y)\ ,\quad p_1(t) =
\begin{cases}
1-|t| , &-1 \leq t \leq 1\
0 , &otherwise
\end{cases}
$$

步骤甚至也写明白了：

  1. 计算 $p_1(t)$ 的 CDF $P_1(t)$
  2. 获得反函数 $P^{-1}(t)$
  3. 把 $[0,1]$ 均匀分布的随机变量 $\xi$ 通过 $P^{-1}(t)$ 映射到 $P^{-1}(\xi)$

于是可以得到代码：

Point2f Warp::squareToTent(const Point2f& sample) {
    auto pdf1 = [](float t) {
        return t < 0.5 ? sqrt(2 * t) - 1 : 1 - sqrt(2 - 2 * t);
    };
    return Point2f(pdf1(sample.x()), pdf1(sample.y()));
}

float Warp::squareToTentPdf(const Point2f& p) {
    auto pdf = [](float t) {
        return abs(t) <= 1 ? 1 - abs(t) : 0;
    };
    return pdf(p.x()) * pdf(p.y());
}

Uniform Disk

单位圆盘上的均匀分布

有两个 $[0,1]$ 的随机变量，考虑使用极坐标表示，一个随机变量用来表示旋转角度，范围在 $[0,2\pi]$ ,另一个再经过计算得到距离远点距离的表达。

考虑距离原点相同距离的位置的周长，$f(r) = 2\pi r$

因为要满足 $p(x) \propto f(x)$

设 $p(x) = af(x)$

$\displaystyle \int_0^1{p(x)dx} = 1$

得到 $p(x) = 2x$

求 CDF 得：$P(x) = \displaystyle \int_0^x{p(t)dt} = x^2$

求反函数：$P^{-1}(t) = \sqrt{x}$

得到代码：

Point2f Warp::squareToUniformDisk(const Point2f& sample) {
    auto theta = sample.x() * M_PI * 2;
    auto r = sqrt(sample.y());
    return Point2f(r * cos(theta), r * sin(theta));
}

float Warp::squareToUniformDiskPdf(const Point2f& p) {
    return sqrt(p.x() * p.x() + p.y() * p.y()) <= 1 ? INV_PI : 0;
}

Uniform Sphere

单位球面上的均匀分布

两个变量一个用来表示旋转角 $\phi$，另一个用来表达仰角 $\theta$

旋转角是均匀的，不需要另外计算。

对于仰角，特定仰角的周长为 $f(x) = 2\pi \sin{x}$

积分计算系数，得：$p(x) = \frac{\sin{x}}{2}$

积分得 CDF：$P(x) = \frac{1-\cos{x}}{2}$

求反函数得：$P^{-1}(t) = \arccos{(1-2t)}$

得到代码：

Vector3f Warp::squareToUniformSphere(const Point2f& sample) {
    auto theta = acos(1 - 2 * sample.x());
    auto phi = sample.y() * M_PI * 2;
    return Vector3f(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
}

float Warp::squareToUniformSpherePdf(const Vector3f& v) {
    return INV_FOURPI;
}

Uniform Hemisphere

均匀分布的半球，方向为 $(0,0,1)$

同上，将仰角的取值缩小到 $[0,\frac{\pi}{2}]$

直接得到代码：

Vector3f Warp::squareToUniformHemisphere(const Point2f& sample) {
    auto theta = acos(1 - sample.x());
    auto phi = sample.y() * M_PI * 2;
    return Vector3f(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
}

float Warp::squareToUniformHemispherePdf(const Vector3f& v) {
    return v.z() < 0 ? 0 : INV_TWOPI;
}

Cosine Hemisphere

有一个余弦作为密度的半球。

根据仰角的密度函数为 $p(\theta) = \frac{\cos{\theta}}{\pi}$

同仰角的圆周上：$f(x) = 2\pi \sin{x} \frac{\cos{x}}{\pi} = \sin{(2x)}$

求系数得：$p(x) = \sin{(2x)}$

积分得 CDF：$P(x) = \frac{1-\cos{(2x)}}{2}$

求反函数：$P^{-1}(t) = \frac{\arccos{(1-2t)}}{2}$

得到代码：

Vector3f Warp::squareToCosineHemisphere(const Point2f& sample) {
    auto theta = acos(1 - 2 * sample.x()) * 0.5F;
    auto phi = sample.y() * M_PI * 2;
    return Vector3f(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
}

float Warp::squareToCosineHemispherePdf(const Vector3f& v) {
    return v.z() < 0 ? 0 : (v.z() * INV_PI);
}

Beckmann

不知道是个什么分布，还是在上半球。

概率密度分布公式给出了：

$$ D(\theta,\phi) = \frac{1}{2\pi} \cdot \frac{2e^{\frac{-\tan^2{\theta}}{\alpha^2}}}{\alpha^2\cos^3{\theta}} $$

虽然写的挺迷惑，参数既带着 $\theta$ 又带着 $\phi$ 但其实这个公式只跟 $\theta$ 有关。

比着上一个再算一遍。

同仰角的圆周上：

$$ f(x) = 2\pi \sin{x}\ D(x,\phi) = \sin{x}\ \frac{2e^{\frac{-\tan^2{x}}{\alpha^2}}}{\alpha^2\cos^3{x}} $$

还给出了 $D(\theta,\phi)$ 满足：

$$ \displaystyle \int_0^{2\pi} \int_0^{\frac{\pi}{2}} D(\theta,\phi) \sin{\theta}\ d\theta\ d\phi = 1 $$

求系数得：$p(x) = f(x)$

甚至还给了提示：

使用变量替换：$x = \cos{x} \quad \tan^2{\theta} = \frac{1-x^2}{x^2}$

以及一个公式：

$$ \displaystyle \int f'(x)\ e^{f(x)}\ dx = e^{f(x)}+C,\ where\ C \in \R $$

积分求 CDF：

$$ P(x) = 1-e^{-\frac{\tan^2{x}}{\alpha^2}} $$

求反函数得：

$$ P^{-1}(t) = \arctan{\sqrt{-\alpha^2\ \ln{(1-t)}}} $$

得到代码：

Vector3f Warp::squareToBeckmann(const Point2f& sample, float alpha) {
    auto phi = sample.x() * 2 * M_PI;
    auto theta = atan(sqrt(-alpha * alpha * log(1.0F - sample.y())));
    return Vector3f(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
}

float Warp::squareToBeckmannPdf(const Vector3f& m, float alpha) {
    if (m.z() <= 0)
        return 0;
    float alpha2 = alpha * alpha;
    float tan_theta2 = (m.x() * m.x() + m.y() * m.y()) / (m.z() * m.z());
    float cos_theta3 = m.z() * m.z() * m.z();
    return INV_PI * exp(-tan_theta2 / alpha2) / (alpha2 * cos_theta3);
}

Part 2: Two simple rendering algorithms

进入正题

Part 2.1: Point lights

点光源。

ajax-simple.xml 文件中多出了如下属性：

<!-- An excerpt from scenes/pa3/ajax-simple.xml: -->
<integrator type="simple">
    <point name="position" value="-20, 40, 20"/>
    <color name="energy" value="3.76e4, 3.76e4, 3.76e4"/>
</integrator>

需要像 normals.cpp 一样新建一个积分器来处理这种类型的东西。

新建了一个 simple.cpp 在 normals.cpp 同目录，再把它也放到 CMakeLists.txt 中去。

相比 normals ，simple 标签中多了一个 point 和一个 color，属性。

在构造函数中解析这两个属性。

class SimpleIntegrator : public Integrator {
   public:
    SimpleIntegrator(const PropertyList& props) {
        m_light_pos = props.getPoint("position");
        m_light_energy = props.getColor("energy");
    }

    // ...

   private:
    Point3f m_light_pos;
    Color3f m_light_energy;
};

NORI_REGISTER_CLASS(SimpleIntegrator, "simple");

然后根据公式补全整个函数即可：

class SimpleIntegrator : public Integrator {
   public:
    SimpleIntegrator(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);
        auto p = m_light_pos;
        auto x = its.p;
        auto xp = p - x;
        if (scene->rayIntersect(Ray3f(x + xp * Epsilon, xp))) {
            return Color3f(0.0f);
        }
        auto phi = m_light_energy;
        auto cos_theta = its.shFrame.n.dot(xp.normalized());
        return phi * INV_FOURPI * INV_PI * std::max(0.0F, cos_theta) / xp.dot(xp);
    }

    std::string toString() const {
        return tfm::format(
            "SimpleIntegrator[\n"
            "  position =
            "  energy =
            "]",
            m_light_pos.toString(), m_light_energy.toString());
    }

   private:
    Point3f m_light_pos;
    Color3f m_light_energy;
};

用时大概 3.4 s

Part 2.2: Ambient occlusion

环境光遮蔽：

说白了就是在 lightbox 中放置的物体。四周光照都是均匀的，明暗区别是由物体自身互相遮挡导致的。

光强公式给出了：

$$ L(x) = \displaystyle \int_{H^2(x)}{V(x,x+\alpha \omega_i)\ \frac{\cos{\theta}}{\pi}}\ d\omega_i $$

其中 V 函数为可见性函数。

这个积分在 x 点所在的半球面上取得。

可以使用蒙特卡罗积分来做。

需要在半球面上采样。由于 V 函数的取值只有 0 和 1，因此 $p(x) \propto f(x)$ 可以取 $p(x) = \frac{\cos{x}}{\pi}$ 。正是在 Part1 中求过的 Cosine Hemisphere 。

这样变成求 $E(V(x,x+\alpha \omega_i))$

代码（需要 #include <nori/sampler.h> 和 #include <nori/warp.h>）：

注意：需要使用 toWorld 将切线空间中的坐标转换到世界空间。

class AoIntegrator : public Integrator {
   public:
    AoIntegrator(const PropertyList& props) {
        /* No parameters this time */
    }

    Color3f Li(const Scene* scene, Sampler* sampler, const Ray3f& ray) const {
        Intersection its;
        if (!scene->rayIntersect(ray, its))
            return Color3f(0.0f);

        auto samp = Warp::squareToCosineHemisphere(sampler->next2D());
        auto poi = its.shFrame.toWorld(samp).normalized();
        return Color3f(!scene->rayIntersect(Ray3f(its.p + poi * Epsilon, poi)));
    }

    std::string toString() const {
        return "AoIntegrator[]";
    }
};

NORI_REGISTER_CLASS(AoIntegrator, "ao");

用时大概 58.0 s （八叉树可能不大够用了）

最后更新于 2023-02-01
啥也没有呀