３次元空間、複数三角形内に均一に、点をばらまく

目的

乱数を散らすとき、複数の三角形にまたがり均一に分布させたい時があります。
例えばパーティクルをメッシュの形状に合わせて発生させたり、
レイトレーシングにおける光源の直接サンプリングといったものが、今回想定する用途です。

乱数

何はともあれ、ソースになる連続した乱数が必要です。
ただ、当然ながら、一般に実数の表現もビットで構成されている以上、離散的な数値になってしまいますが。
あまり踏み込むと脱線しすぎるので、ここについては、コードを張るだけにします。
以下のような乱数モジュールを使います。
エンジンをテンプレートで挿げ替えることだけは想定しています。
まあ、状況によっては必要ないかもしれませんが。

    struct Xor {
        Xor() {

        }
        Xor(uint32_t seed) {
            _y = std::max(seed, 1u);
        }
        uint32_t generate() {
            _y = _y ^ (_y << 13); _y = _y ^ (_y >> 17);
            return _y = _y ^ (_y << 5);
        }
    private:
        uint32_t _y = 2463534242;
    };
    template <class Generator>
    struct RandomEngine : public Generator {
    public:
        RandomEngine() :Generator() {
        }
        RandomEngine(uint32_t seed) :Generator(seed) {
        }

        void discard(int n) {
            for (int i = 0; i < n; ++i) {
                this->generate();
            }
        }

        double continuous() {
            uint32_t uniform = this->generate();
            constexpr double c = 1.0 / static_cast<double>(0xffffffffLL + 1);
            return static_cast<double>(uniform) * c;
        }
        // 0 <= x < 1
        double continuous(double a, double b) {
            return a + continuous() * (b - a);
        }
        // a <= x < b
        double continuous(double a, double b) {
            return a + continuous() * (b - a);
        }
        // その他のメソッド省略
    };

たとえばクライアント側はこうなりますね

    lc::RandomEngine<lc::Xor> e;
    e.discard(100);
    float value = e.continuous();

ちなみにdiscardというのは、初期化直後の乱数値とシードの値との関連性をある程度断ち切るためにある、スキップ用メンバ関数です。
一応乱数によって読み飛ばすべき量は変わってきますが、面倒なので今回は横着して適当な数を使います。

もちろんC++標準のrandomを用いても問題ありません。
その場合、こちらが参考になるでしょう。

C++0xの新しい乱数ライブラリ、random
https://cpplover.blogspot.jp/2009/11/c0xrandom.html

今回は自前で組み立ててみました。

三角形のデータ型

一応後々のコードのために三角形のデータ型を用意します。

    struct Triangle {
        Triangle() {}
        Triangle(const Vec3 &v0, const Vec3 &v1, const Vec3 &v2) :v{ v0, v1, v2 } {

        }
        Vec3 &operator[](std::size_t i) {
            return v[i];
        }
        const Vec3 &operator[](std::size_t i) const {
            return v[i];
        }
        std::array<Vec3, 3> v;
    };

単一三角形内に一様に分布させる

ググるとすぐに出てきたため、ひとまずこちらを使うことにします。

0 <= \varepsilon_1 < 1\\
0 <= \varepsilon_2 < 1

の一様乱数を使って、点ABC上の点Pは、

P = A(1 - \sqrt{\varepsilon_1}) + B(\sqrt{\varepsilon_1}(1 - \varepsilon_2)) + C\sqrt{\varepsilon_1}\varepsilon_2 

で得ることができます。
なんでこれでいいのかは・・・後々調べることにします。

実装は以下のようにしました。

    template <class Generator>
    inline Vec3 uniform_on_triangle(RandomEngine<Generator> &e, const Triangle &tri) {
        double eps1 = e.continuous();
        double eps2 = e.continuous();

        double sqrt_r1 = glm::sqrt(eps1);
        lc::Vec3 p =
            (1.0 - sqrt_r1) * tri[0]
            + sqrt_r1 * (1.0 - eps2) * tri[1]
            + sqrt_r1 * eps2 * tri[2];
        return p;
    }

複数三角形に対応する

複数三角形に対応するにはどうしたらいいでしょうか？

単純に複数の三角形の中から一様に一つ選ぶ

とするとひとまずはばらまくことができます。
ですが、すぐに気づきますが、面積が小さな三角形の点密度が高く、逆に面積の大きいポリゴンは点密度が低くなってしまいます。これは微妙ですね。
もちろんポリゴンの面積がほとんど均一のものを扱う場合で、分布の精度がさほど必要とされない場面では、これでよいでしょう。妥当な実装であるといえます。

面積に比例する確率分布で三角形の中から一つ選ぶ

一歩進んで、面積に対して確率を比例させれば、ポリゴンの大きさがいくら偏っていても、均一にばらまくことができそうです。

三角形の面積を求める

面積に比例させるためにはまずは面積を求める必要があります。
これには、
二つのベクトルの外積のノルムは、そのベクトルのなす平行四辺形の面積に等しい
という性質を使います。
なので、三角形の場合は半分にすればよいだけです。

Area = \frac{\|\vec{a}\times\vec{b}\|}{2}

実装は以下のようにしました。

    inline double triangle_area(const Vec3 &p0, const Vec3 &p1, const Vec3 &p2) {
        auto va = p0 - p1;
        auto vb = p2 - p1;
        return glm::length(glm::cross(va, vb)) * 0.5;
    }

確率分布を制御する

面積に対して確率を比例させるには、それぞれの面積を累積的な数値に変換して、テーブルに格納し、
それを参照して偏らせればよさそうです。
より具体的に例をあげます。

三角形Aの面積: 10
三角形Bの面積: 20
三角形Cの面積: 30

のとき、累積面積テーブルは、

Aの累積面積: 10 (10)
Bの累積面積: 30 (10 + 20)
Cの累積面積: 60 (10 + 20 + 30)

となります。
そして、
乱数を0 <= x < 60で発生させます。

そのとき、xの値と、それにより選ばれる三角形との関係性は、

0 <= x < 10 のときA
10 <= x < 30 のときB
30 <= x < 60 のときC

このように決めてあげればよさそうです。
ただ、ABCどれを選ぶか選定するときに、２分探索したほうが、ポリゴンの個数が多い場合は有利になるでしょう。C++には、ある数より大きい数値をもつ最初の要素を探索するためのstd::upper_bound関数があるため、実装を省略できます。
余談ですが、連続的な値と離散的な値を連携させるときに、なんだかんだとupper_boundというのは便利に使えますね。

さて、それを考慮にいれて実装するとたとえば以下のようになるでしょう。

    class UniformOnTriangle {
    public:
        void build() {
            if (_triangles.empty()) {
                _cumulative_areas.clear();
                _area = 0.0;
                return;
            }

            _cumulative_areas.resize(_triangles.size());
            double area = 0.0;
            for (int i = 0; i < _triangles.size(); ++i) {
                const Triangle &tri = _triangles[i];
                area += triangle_area(tri[0], tri[1], tri[2]);
                _cumulative_areas[i] = area;
            }
            _area = area;
        }

        template <class Generator>
        Vec3 uniform(RandomEngine<Generator> &e) const {
            double p = e.continuous(0.0, _area);
            auto it = std::upper_bound(_cumulative_areas.begin(), _cumulative_areas.end(), p);
            std::size_t index = std::distance(_cumulative_areas.begin(), it);
            index = std::min(index, _cumulative_areas.size() - 1);
            return uniform_on_triangle(e, _triangles[index]);
        }

        void set_triangle(const std::vector<Triangle> &triangles) {
            _triangles = triangles;
        }
        double get_area() const {
            return _area;
        }
        std::vector<Triangle> _triangles;
        std::vector<double> _cumulative_areas;
        double _area = 0.0;
    };

結果

以上をcinderやopenframeworks上でOpenGLでビジュアライズテストします。
やっぱりこういうときは生のGLを扱うと環境構築だなんだとまったく面倒なので、助かりますね。win, mac選ばないのもいいところです。
ちょっと余談でした。結果は以下のようになります。均一に分布していることがわかります。

一応Cinderでのテスト用実装コードを載せておきます。

#include "cinder/app/App.h"
#include "cinder/app/RendererGl.h"
#include "cinder/gl/gl.h"
#include "cinder/CameraUi.h"

#include "random_engine.hpp"
#include "triangle_area.hpp"
#include "uniform_on_triangle.hpp"

#include <boost/format.hpp>

#define TINYOBJLOADER_IMPLEMENTATION
#include <tiny_obj_loader.h>

using namespace ci;
using namespace ci::app;
using namespace std;

class TriangleRandomApp : public App {
public:
    void setup() override;
    void keyDown(KeyEvent event) {
        _index++;
    }
    void mouseDown(MouseEvent event) override;
    void update() override;
    void draw() override;

    CameraPersp _camera;
    CameraUi _cameraUi;
    gl::BatchRef _plane;

    int _index = 0;

    lc::UniformOnTriangle _uniformTri;
};

void TriangleRandomApp::setup()
{
    _camera.lookAt(vec3(0, 0.0f, 4.0f), vec3(0.0f));
    _camera.setPerspective(40.0f, getWindowAspectRatio(), 0.01f, 100.0f);
    _cameraUi = CameraUi(&_camera, getWindow());

    auto colorShader = gl::getStockShader(gl::ShaderDef().color());
    _plane = gl::Batch::create(geom::WirePlane().size(vec2(10.0f)).subdivisions(ivec2(10)), colorShader);


    std::vector<tinyobj::shape_t> shapes;
    std::vector<tinyobj::material_t> materials;
    std::string err;
    std::string path = (getAssetPath("") / "poly.obj").string();
    bool ret = tinyobj::LoadObj(shapes, materials, err, path.c_str());

    const tinyobj::shape_t &shape = shapes[0];
    std::vector<lc::Triangle> triangles;
    for (size_t i = 0; i < shape.mesh.indices.size(); i += 3) {
        lc::Triangle tri;
        for (int j = 0; j < 3; ++j) {
            int idx = shape.mesh.indices[i + j];
            for (int k = 0; k < 3; ++k) {
                tri.v[j][k] = shape.mesh.positions[idx * 3 + k];
            }
        }
        triangles.push_back(tri);
    }

    _uniformTri.set_triangle(triangles);
    _uniformTri.build();

    // 面積検証コード
    /*
    lc::RandomEngine<lc::Xor> e;
    e.discard(100);
    for (int i = 0; i < 100; ++i) {
        lc::Vec3 p0(0.0, 0.0, 0.0);
        lc::Vec3 p1(10.0, 0.0, 0.0);
        lc::Vec3 p2(10.0, 10.0, 0.0);
        double base_area = 50.0;

        float scale = e.continuous(0.1f, 10.0f);

        // 回転、移動が行われても、正しく面積を求めることができるはずである
        lc::Vec3 translation(e.continuous(), e.continuous(), e.continuous());
        lc::Quat rotation = glm::angleAxis(e.continuous(0.0, glm::pi<double>()), e.on_sphere());
        p0 = rotation * (p0 + translation);
        p1 = rotation * (p1 + translation);
        p2 = rotation * (p2 + translation);

        p0 *= scale;
        p1 *= scale;
        p2 *= scale;

        double area = lc::triangle_area(p0, p1, p2);

        // スケールがx倍されるとき、面積はx^2倍になるはずである
        double d = area - base_area * scale * scale;
        if (0.000001 < d) {
            abort();
        }
    }
    */
}

void TriangleRandomApp::mouseDown(MouseEvent event)
{
}

void TriangleRandomApp::update()
{
}

void TriangleRandomApp::draw()
{
    gl::clear(Color(0, 0, 0));

    // Set up the camera.
    gl::ScopedMatrices push;
    gl::setMatrices(_camera);

    {
        gl::ScopedColor color(Color::gray(0.2f));
        _plane->draw();
    }

    gl::VertBatch vb(GL_POINTS);

    static lc::RandomEngine<lc::Xor128> e_triin;
    for (int i = 0; i < 10000; ++i) {
        lc::Vec3 p = _uniformTri.uniform(e_triin);
        vb.vertex(p);
    }

    vb.draw();

    for (int i = 0; i < _uniformTri._triangles.size(); ++i) {
        auto tri = _uniformTri._triangles[i];
        for (int j = 0; j < 3; ++j) {
            gl::drawLine(tri[j], tri[(j + 1) % 3]);
        }
    }
}

CINDER_APP(TriangleRandomApp, RendererGl, [](App::Settings *settings) {
    settings->setConsoleWindowEnabled(true);
})

単一三角形内に一様に分布させる（さらに効率の良い方法）

gilbert_yumu さんより、もっとシンプルで高速な方法があると教えていただきまして、そちらも紹介いたします。式は以下のようなシンプルなものです。

\alpha \leftarrow min\left( \varepsilon _{ 1 },\varepsilon _{ 2 } \right) \\ \beta \leftarrow max\left( \varepsilon _{ 1 },\varepsilon _{ 2 } \right) \\ P=A\alpha +B\left( 1-\beta  \right) +C\left( \beta -\alpha  \right) 

これでも等価な処理だというのですから驚きです。

[出典]
gilbert_yumu さん 三角形内の一様乱数
William E. Stein, Matthew F. Keblisb, "A new method to simulate the triangular distribution"

手元の雑ベンチマークですが、10,000,000個ほど点を両方で生成してみると、

	時間
sqrtを使った方法	160ms
min, maxを使った方法	120ms

※乱数生成 32bit xorshiftのコストも入ってしまっています

とやはりSqrtを削ることのできるメリットが見て取れました。
これは積極的に乗り換えていきたいところです。

gilbert_yumu さんありがとうございました。