はじめに

3DCGのモデリングにあたっては、多角形の面をたくさん集めたポリゴン形式のデータが良く使われています。多くの3Dモデリングソフトも、ポリゴンメッシュに最適化したデータ構造や表示方法を確立していると思います。
一方でポリゴンのような点と面による表現方法は、曲面を少ない情報で表すのが苦手です。下の画像はそれぞれ面数800, 200, 50で表現した球を表示しています。スムーズシェーディングで表示すると面数800のモデルはほとんど球と変わりないですが、面数200のモデルは違和感を感じ、面数50だとだいぶ厳しいです。

しかし、球は関数を使って簡潔に表すことが出来ます。原点を中心とする半径1の球であれば $ x^2 + y^2 + z^2 = 1 $ でしたね。このように、球の情報が欲しければ中心座標と半径（とそれが球であること）さえ分かれば十分なのです。それだけあればモデルのデータサイズもずっと抑えることができますし、むしろ離散的に点を取っているよりも正確です。実際OpenGLも球を描画したいときは、glTranslateで中心を決めてglutSolidSphereで半径と分割数を指定しますし、一部のファイル形式は球の情報を保持するものがあるはずです（FBXとか？自信ないです）。

また、たとえ球としてモデルの情報を保持していたとしても、多くのソフトでは関数をポリゴンに分割してから表示しています（先程のglutSolidSphereの分割数はまさにそれです）。ポリゴンの表示を高速に行う手法が確立しているのでそれはそれでいいのですが、一方で関数の情報をそのまま表示に使う方法も存在します。

ということで、今回は関数で表現されたモデルをそのまま表示する、ということをやってみようと思います。

陰関数曲面（Implicit Surface）

関数といっても、陽関数（例：$y = ax + b$ ）やパラメータ関数（例：$x = cos\theta , y = sin\theta $）、陰関数（例：$x+y = 0$）などいろいろあります。その中でも陰関数は曲面の内外を明確にできる、勾配が求めやすいなどの利点があるので（そのありがたみは後述します）陰関数を使います。今回は扱う世界が3次元なので、 $f(x, y, z) = 0$ を満たす点 $(x, y, z)$ の集合が面となるような陰関数曲面を取り扱います。

例えば、半径 $r$ の球面であれば

    f(x, y, z) = x^2 + y^2 + z^2 - r^2 = 0

トーラス面（円をある軸のまわりに1回転させてできるドーナツみたいなもの、以下の式は $z$ 軸中心に $R$ 離れた $xz$ 平面上の半径 $r$ の円を回転したもの）であれば

    f(x, y, z) = (x^2 + y^2 + z^2 + R^2 - r^2)^2 - 4R^2 (x^2 + y^2) = 0

という陰関数で表せます。

モデリングから表示までの流れ

前置きが長くなってしまい、この後の話の展開が見えづらくなってきたのでここでどういう流れで実装したかを書いておきます。

環境

Ubuntu 20.04
C++14　CMake 3.16.3
OpenCV 4.2.0 (レンダー画像作成用)

流れ

各モデルの設計
- 今回はいわゆるプリミティブ型である球・直方体・円柱・トーラスと、二次曲面を用意しました。
- Sphere TracingとMatcap shadingのために、各モデルに対して符号付き距離関数と法線ベクトルを設定します。
カメラの位置・光線の方向を決める
符号付き距離関数を計算し逐次的に光線とモデルとの衝突位置を判定（Sphere Tracing）
衝突位置でのモデルの法線を計算し、光線の方向と加味して色付け（Matcap shading）

それでは、実装していきましょう。

陰関数モデリング

陰関数はそのまま用いるわけではなく、上でも説明した通り符号付き距離関数と法線ベクトルがレンダリングのために使われます。

符号付き距離関数

英語ではSigned Distance Functionと表記し、SDFという略記もよく用いられます。これはある点 $\boldsymbol v$ から曲面までの最短距離を返す関数であり、符号はモデルの外部の点なら正、内部の点なら負が返されます。

    \rm{sdf} (\boldsymbol v) = \min_{\boldsymbol x \in S} d(\boldsymbol v, \boldsymbol x)

ところが一般的に $ f(\boldsymbol v) $ と $ \rm{sdf} (\boldsymbol v) $ は一致せず、かつ $ f(\boldsymbol v) $ から $ \rm{sdf} (\boldsymbol v) $ を解析的に出すのも難しいです。かといって数値計算を用いるのは計算量が増えてしまいますので、今回はこちらを参考に、SDFの近似解を解析的に出すことにしました。

    \rm{sdf} (\boldsymbol v) := \delta \approx \frac{f(v)}{\| \nabla f(v) \|} := \delta_1

これは $\boldsymbol v$ から曲面 $S$ 上の最も近い点 $\boldsymbol x$ 周りのテイラー展開を１次で打ち切ったものであり、従って $\boldsymbol v$ が $S$ に近いほど正確になります。また、論文では $\delta_1$ が $\delta$ を超えないことについて記されていました。後述のSphere Tracingでは、距離関数を多めに少なめに見積もる方が多めに見積もってしまうよりも弊害が小さいので、今回こちらを採用しています。

法線ベクトル

こちらはもっと簡単です。$f(\boldsymbol v)$ に対して、 $ \boldsymbol n(\boldsymbol v) = \nabla f(\boldsymbol v) = [\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}, \frac{\partial f}{\partial z}]$ を求めて正規化するだけですね。

モデル設計

というわけで、作っていきましょう。全てのモデルはSDFと法線ベクトルを定義してほしいので、最初に抽象クラスを作っておくと見通しがいいですね。SDFはある頂点からモデルまでの最短距離を返し、法線ベクトルはその位置のモデルのベクトルを返すのでこんな感じです。

model.h

class Model {
public:
    Model() {}
    virtual double sdf(Vector3 v) {}
    virtual Vector3 normal(Vector3 v) {}
}

Vector3 クラスについての説明を忘れていました。配列や std::array よりも扱いやすいよう、いくつかの演算子をオーバーロードし、メンバ関数（ノルム・内積・外積・正規化など）を登録しているクラスです。

球

一番分かりやすいです。中心の位置ベクトル $\boldsymbol c$ と半径 $r$ がパラメータになります。SDFはその点から中心までの距離、から半径を引いたものです（球面上なら0になり内部は負、外部は正となる）。法線ベクトルは球の方程式の勾配を取ります。

model.h

class Sphere : public Model {
public:
    Sphere(double x, double y, double z, double r);
    Sphere(Vector3 c, double r);
    double sdf(Vector3 v) override;
    Vector3 normal(Vector3 v) override;
private:
    Vector3 c;
    double r;
};

model.cpp

Sphere::Sphere(double x, double y, double z, double r) : 
    c(Vector3(x, y, z)), r(r) {}
Sphere::Sphere(Vector3 c, double r) : c(c), r(r) {}
double Sphere::sdf(Vector3 v) {
    return v.dist(c) - r;
}
Vector3 Sphere::normal(Vector3 v) {
    return Vector3(2 * (v.x - c.x), 2 * (v.y - c.y), 2 * (v.z - c.z)).normalize();
}

直方体

Axis alignedに限定しました。この場合は座標が最小となる点と、各辺の長さがあればモデリングできます。Axis alignedでない場合は、隣り合う2辺のなす角が直角という制約があるので、最低限のパラメータがいくつで十分なのかを考える必要がありそうです。ちなみに4つのVector3を用意すれば、平行六面体がモデリングできますね。

SDFは、$x, y, z$ でそれぞれ場合分けして計算して求めます。

model.cpp

double Rectangular::sdf(Vector3 v) {
    double d_x = std::max(v.x - a.x, o.x - v.x);
    double d_y = std::max(v.y - b.y, o.y - v.y);
    double d_z = std::max(v.z - c.z, o.z - v.z);
    if(d_x < 0.0 && d_y < 0.0 && d_z < 0.0){  // inside
        return std::max(std::max(d_x, d_y), d_z);
    }
    else{  // outside
        d_x = std::max(d_x, 0.0);
        d_y = std::max(d_y, 0.0);
        d_z = std::max(d_z, 0.0);
        return Vector3(d_x, d_y, d_z).norm();
    }
}

円柱

軸によらずに形状を作れるようにしました。2つの底面の中心と底面となる円の半径がパラメータです。長くなってしまうのでSDFだけ載せます。コメントアウトにあるとおり、円柱の軸を法線ベクトルとして点 $\boldsymbol v$ を通る平面がどこで交わるかを判定しています。それが円柱の内か外かで場合分けしています。

model.cpp

double Cylinder::sdf(Vector3 v) {
    Vector3 dir = c2 - c1;
    double height = dir.norm();
    /* compute the intersection of line and plane
        line: x = c1 + t * d
        plane: d・(v - x) = 0 */
    double t = dir.dot(v - c1) / (dir.dot(dir));
    double side_dist = v.dist(c1 + dir * t) - r;
    
    if(t >= 0 && t <= 1){
        return side_dist >= 0 ? side_dist : std::min(side_dist, std::min(height * t, height * (1-t)));
    }
    else{
        if(t > 1){
            return side_dist >= 0 ? std::sqrt(side_dist * side_dist + height * (t-1) * height * (t-1)) : height * (t-1);
        }
        else{
            return side_dist >= 0 ? std::sqrt(side_dist * side_dist + height * (0-t) * height * (0-t)) : height * (0-t);
        }
    }
}

トーラス

関数は複雑ですが、場合分けが要らないので記述量は直方体や円柱より少なく書けました。上に書いてある方程式を微分するだけです。

2次曲面

ax^2 + by ^2 + cz^2 + dxy + eyz + fzx + gx + hy + iz + j = 0

で表される関数です。10のパラメータだけで様々な曲面を表現できます。しかし閉曲面とならない場合は多分うまく表示できないです。楕円面はうまくできると思います。

ブーリアン演算

陰関数モデリングでは、形状同士を足したり引いたり共通部分を取ったり、という操作をSDFの最小または最大をとるだけで実現できます。

    // 和集合
    double Union::sdf(Vector3 v){
        return std::min(m1->sdf(v), m2->sdf(v));
    }
    // 共通部分
    double Intersection::sdf(Vector3 v){
        return std::max(m1->sdf(v), m2->sdf(v));
    }
    // 差集合
    double Difference::sdf(Vector3 v){
        return std::max(m1->sdf(v), -m2->sdf(v));
    }

これらも Model を継承して実装しました。ツリーが深くなるとSDFやnormalで無駄な計算が増えてしまう実装になっているので、工夫したいです。

boolean.h

class Union : public Model {
public:
    Union(Model* m1, Model* m2): m1(m1), m2(m2) {}
    double sdf(Vector3 v){
        return std::min(m1->sdf(v), m2->sdf(v));
    }
    Vector3 normal(Vector3 v){
        if(m1->sdf(v) > m2->sdf(v)) return m2->normal(v);
        else return m1->normal(v);
    }
private:
    Model *m1, *m2;
};

レンダリング

モデリングだけでは完了しないので、それを表示させるためにレンダリングの部分を実装していきます。今回は画像処理のライブラリとしてOpenCVを用いています。Model 型のmodel とmatcap用の画像を用いて、image の画素を埋めていくという流れです。

render.cpp

void render(cv::Mat &image, Model &model, cv::Mat &matcap_img){
    int height = image.rows;
    int width = image.cols;

    std::vector<std::vector<Vector3>> ray_d(height, std::vector<Vector3>(width, Vector3(0,0,0)));
    decide_ray_direction(ray_d, width, height);

    std::vector<std::vector<Vector3>> nrms(height, std::vector<Vector3>(width, Vector3(0, 0, 0)));
    std::vector<std::vector<bool>> collided(height, std::vector<bool>(width, false));
    sphere_tracing(model, ray_d, nrms, collided, width, height);

    matcap_texture(image, matcap_img, ray_d, nrms, collided, width, height);

光線の方向を決定

Ray tracingと同様に、各画素に入ってくる光線の動きを逆向きに考えたいので、各画素の光線の方向を定めます。今回は透視投影、つまりある1点を源としてそこからの方向を求めます。中心画素は光源から原点に向かうところを映し、そこから画素を広げていくイメージです（変数名のセンスの無さは見逃してほしいです）。

render.cpp

Vector3 camera_o = Vector3(1, 1, 1);  // 適当
Vector3 camera_t = Vector3(0, 0, 0);

void decide_ray_direction(std::vector<std::vector<Vector3>> &ray_d, int width, int height){
    Vector3 v_view = camera_t - camera_o;
    Vector3 v_right = (v_view.x != 0 || v_view.y != 0) ? Vector3(-v_view.y, v_view.x, 0).normalize() : Vector3(1, 0, 0);
    Vector3 v_up = v_right.cross(v_view).normalize();

    double real_h = 5;  // 適当
    double real_w = real_h * width / height;
    for(int i = 0; i < height; i++){
        Vector3 h_translate = v_up * (real_h * (i-height/2) / height);
        for(int j = 0; j < width; j++){
            Vector3 w_translate = v_right * (real_w * (j-width/2) / width);
            ray_d[i][j] = ((camera_t + h_translate + w_translate) - camera_o).normalize();
        }
    }
}

Sphere Tracing

これは、現在地でのモデルのSDFを求め、その値だけ方向ベクトルに沿って光線を進めていくRay marchingの1種です。私はこの方法を初めてみた時、陰関数と光の直線との方程式を求めて衝突点を求める方法でいいのではと思っていたのですが、その解を求める数値計算のステップよりこの方法のステップの方が収束が早くなるケースが多いのかなと解釈しています。光線とモデルが垂直に交わらない場合はSDFがずーっと微小のまま0にならないので、ある閾値より小さくなったら収束したと判断します。

（引用元：[Keinert et al. 2014]）
以下のコードでは最小値のほかに、（前ステップのSDFとの差の）最大値とステップ数も、収束条件として設定しています。これら2つは慎重に選ぶべきで、特に最大値は何ステップ連続で増えている（モデルから遠ざかっているので衝突の見込みがない）、みたいな条件にすることもできると思います。

render.cpp

const int MAX_STEP = 100;
const double FINISH_MINIMUM = 0.001;
const double FINISH_MAXIMUM = 100;

void sphere_tracing(Model &model, std::vector<std::vector<Vector3>> &ray_d, std::vector<std::vector<Vector3>> &nrms, 
    std::vector<std::vector<bool>> &collided, int width, int height){

    std::vector<std::vector<Vector3>> coords(height, std::vector<Vector3>(width, camera_o));

    for(int i = 0; i < height; i++){
        for(int j = 0; j < width; j++){
            // John C. Hart 1995
            double sdf = model.sdf(coords[i][j]);
            double old_sdf = 1e18;
            int step = 0;
            while(step < MAX_STEP && sdf > FINISH_MINIMUM && (sdf - old_sdf) < FINISH_MAXIMUM){
                coords[i][j] += ray_d[i][j] * sdf;
                old_sdf = sdf;
                sdf = model.sdf(coords[i][j]);
                step += 1;
            }
            coords[i][j] += ray_d[i][j] * sdf;
            if(std::abs(sdf) <= FINISH_MINIMUM){
                collided[i][j] = true;
                nrms[i][j] = model.normal(coords[i][j]);
            }
        }
    }
}

なお元の論文がこちらにあるので、あわせてご参照ください。

Matcap texture

実装はしましたが、既存のコードのほぼ模倣でまだあまり理解できてないので、省略させてください。すごくざっくり説明すると、光線の方向ベクトルと衝突点におけるモデルの法線ベクトル（どちらも単位ベクトル）のなす角を計算し、テクスチャ画像のUV座標から色を取ってきています。

結果

わりかし綺麗にできた円筒と今回オリジナルで作った2次曲線のモデルの結果を載せます。

他にもGitHubに見た目がよかったものをアップロードしています（フルのコードもあります）。カメラの角度によっては、球の見た目が変になることがありました。モデリングは間違っていないはずなので、色付けのコードに変なところがありそうです。もしくは光源を用意して光量をちゃんと計算したり、投影面をしっかり作ったりということが必要なのかもしれません。他にもエイリアシングが発生しているなどの見た目の悪さもあるので、レンダリングをもっとちゃんと勉強しようと思います。

参考