Groups

Chapter 13

シーングラフを作る。

f:id:mtXTJocj:20210703202735p:plain

Tree 構造

シーングラフといっても処理を複雑にしないよう、実際には Graph ではなく Tree 構造を用いる。今回の目的に合致する Tree は

  • Node が持つ情報は Shape によって異なる
  • root を除く各 Node は親 Node へのリンクを持つ
  • Tree を構成する Node は任意の数の子 Node を持つ
  • Node に対して子 Node(sub tree) を加えることが可能

となる。

これを Rust で実現することを考える。

Rust で Tree を作るとなると、Node に持たせる情報は Generics を使ってユーザが決められるようにすることが多い。全 Node が同一の情報を持つのであれば問題ないが、今回のように、この NodeSphere, あの NodeCylinder といったように Node ごとに持たせる情報を変えたい場合には使えない。もちろん、enum を使えば可能ではあるが、その場合全ての Node のサイズがとり得る最大のものになってしまい、大規模な Tree には向かない。

なので、今回は Node のメンバに trait object である Shape を持たせることにした。trait object はコンパイル時にサイズがわからないので boxing する。

    shape: Box<dyn Shape>,

他、Shape の種類によらず共通して持つ必要のある情報は Node に持たせる。具体的には Transform と親 Node へのリンクがそれにあたる。

結果として Node のメンバは以下のようになった。

pub struct Node {
    /// 親 Node
    parent: Option<NonNull<Node>>,
    /// 親 Node の座標系への変換
    transform: Transform,
    /// 本体
    shape: Box<dyn Shape>,
}

ここで親 Node へのリンクである parent はポインタを使用する。

親 Node へのリンク

Rust での Tree 構造の実現方法を調べると簡単なのは親 Node へのリンクを持たないものになる。常に root からトラバースするのであれば、親 Node へのリンクも不要。そうでない場合に親 Node へのリンクを持つようにすると、ownership まわりが複雑になる。一番自然なのは親 Node が子 Node を所有する。となると、子 Node から 親 Node へのリンクは循環参照を避けるためにも ownership を持たない形にしなくてはならない。

ownership を持たない参照を実現するために最初は Weak を使ってみたけれど、そのためには Rc::downgrade()Weak を作れるよう、NodeRc で包まなくてはならない。また Rc だと owner が複数ある可能性からコンパイル時に borrow check が行えず、そのままでは read only になってしまう。それを避けようとするとさらに RefCell を使うことになり、Rc<RefCell<Node>> とややこしい。

これ自体は Rust のイディオムのようなものだから我慢するという選択肢もあるけれど、親 Node へのリンクが欲しいだけの割にどうにも大仰な感じがする。root 以外の Node の owner は 親 Node のみという点は変わらないわけだし、子 Node が親 Node より長生きしてリンクが無効になることもない。なので、代わりにポインタを使った方がすっきり書けそうな感じがした。C や C++ に毒されすぎと言われそうな気もするが。

Node へのリンクにポインタを使うのであれば、親 Node のアドレスが変わってしまっては困る。なので Node を直接 move できないようにする。具体的には常に Box<Node> の形で使うようにする。これなら Node 自体は heap に確保され Box を move しても Node のアドレスは変わらない。というわけで、new() の時点で boxing してしまう。

    pub fn new(shape: Box<dyn Shape>) -> Box<Self> {
        Box::new(Node {
            parent: None,
            transform: Transform::identity(),
            shape,
        })
    }

Node のメンバは private だから外で直接 Node を作ることはできず、必ず new() で作ることになり、boxing されている。

問題は Node を作ってしまうと、 shape が元々何であったか(Sphere, Cube, etc)がわからなくなってしまい、固有の操作が行えなくなってしまうことだが、今回は一部例外を除き Node を作る前に設定を終えておくことで対処する。

子 Node の追加

先に挙げた Node は子 Node を保持するメンバを持っていない。これは Sphere 等は常に Leaf Node で子 Node を持たないことによる。そこで子 Node を持つ Shape として別途 Group を用意する。

pub struct Group {
    /// 子 Node
    children: Vec<Box<Node>>,
}

何だか Vec で heap を使っているのに、更に Box まで使っていると無駄っぽいが、子 Node に加える前後で Node のアドレスを変えないようにするためなので我慢する。

Node を加えるには add_child() を使用するが、これは Group のメソッドではなく Shape に用意する。これは Node から add_child() できるようにするため。そうしないと子 Node に親 Node のポインタをセットできない。

    pub fn add_child(&mut self, mut child: Box<Node>) {
        child.parent = NonNull::new(&mut *self);
        self.shape.add_child(child);
    }

なので、Shapeadd_child() を直接呼ぶことはない。とはいえ、現在 Group 以外の Nodeadd_child() を呼ぶと panic してしまうため、Error を返すくらいはしたい。

その他

NodeTransform 移した際、transform_mut() で新しい Transform をセットするのではなく、 set_transform() で設定するようにした。前者は &mutTransform を得るのに対し、後者はセットする Transform を引数にとる。
こうした理由は Transform をセットしたタイミングに合わせて他の処理をする場合に備えるため。

例えば、現在は法線ベクトルを求める normal_at() の中で行う World -> local 変換は各 Node の変換を root に向かってたどることで実現しているが、毎回これでは効率が悪い。対応策として各 Node には World 変換をキャッシュとして持たせることが考えられるが、Transform セットのタイミングでこのキャッシュは invalid なものになってしまう。となるとキャッシュの更新なり、キャッシュが無効である旨を示すフラグのセットなりを行う必要がある。transform_mut() ではそれができないため set 用のメソッドに変更した。

Rust で hoge というメンバに対する getter, setter は hoge()hoge_mut() とすることが多いようだが、上記要求には対応できない。こうしたときに Rust ならこうするというパターンがあるのだろうか。
hogeset_hoge()hogerahogera_mut() を使うなんて面倒で仕方ない。

雑感

Node の実装方針が本当にこれで良かったのか今でも自信がない。Rust による Tree の実装を色々と見てみたが、parent へのリンクを持たないとか、Node は全て同一の構造とか、Node は root からのみ追加可能で中間 Node に対して部分木を追加しないとか、子 Node の数は固定とか、どれも限定的な使用法を想定しているものばかりで、これぞ決定版と思えるものにはついぞめぐり会えなかった。そのため、今回の用途に適用しようとすると、そのままではどれも使えなかった。

他にも Box<Node> とそのメンバである Box<Shape> はライフサイクルが一致するため、別々に heap の allocation をするのではなく、一回の allocation で両方合わせたサイズの heap を確保して NodeShape に分配したいと思ってしまうが、Rust ではどのようにすれば良いのか分からない。

賢い人ならこの程度のことはすぐに解決できるのかもしれないが、ちょっと自分には無理だった。Rust は頭の良い人しか使ってはいけない言語なのかも。Haskell みたいに。

印象として実行時に自由な編集を許すタイプのプログラムを Rust で作ろうとすると、多大な困難を伴う感じがする。そういった意味では今回の raytracer なんてかわいいものだが、3D CG のモデリングツールを Rust で作るといったときにはどんな感じになるのだろう。 他、自由に構造をいじれるグラフ構造の実装なんて考えたくもない。Typed Arena だと全 Node の lifetime が同じになってしまうので、グラフの一部を削除とかはしづらいし。

f:id:mtXTJocj:20210627115251p:plain

Cylinders

Chapter 13

Cylinders

intersection

ここも Cube のときと同様にあまり書くことがない。ただ、Cube では割とていねいに intersection の解説がされていたのに、こちらでは天下り式にやり方が書かれていて写経するだけになっていたため、その部分をここに残しておく。

Cylinder

円筒部分と Ray との交点を求める。Ray  \mathbf{r}


\mathbf{r}  = \mathbf{d}t + \mathbf{o}

で、y 軸方向にのびる半径 1 の円筒は任意の y において


x^2 + z^2 = 1

となる。よって


\eqalign{
  {r_x}^2 + {r_z}^2 \cr
  &= {({d_x}t + {o_x})}^2 + {({d_z}t + {o_z})}^2 \cr
  &= {d_x}{t}^2 + 2{d_x}{o_x}t + {o_x}^2 + {d_z}{t}^2 + 2{d_z}{o_z}t + {o_z}^2 \cr
  &= {({d_x}^2 + {d_z}^2)}t^2 + 2({d_x}{o_x} + {d_z}{o_z})t + {o_x}^2 + {o_z}^2 \cr
  &= 1
}

これを  t二次方程式  a{t}^2 + bt + c = 0 と見れば、


\eqalign{
  a &= d_x^2 + d_z^2 \cr
  b &= 2(d_x o_x + d_z o_z) \cr
  c &= o_x^2 + o_z^2 - 1
}

としてあとは解の公式に当てはめるだけ。 a = 0 の時には円筒と Ray が平行なので、交点なしとする

y 軸方向に無限に長い円筒ならこれで終わりだけれど、今回実装した Cylinder は minimum と maximum を設定でき、その場合には有限の長さを持つ。
とはいえ、単に y の有効範囲があるというだけ。つまり、上記の二次方程式を解いて得た  t_0, t_1 を Ray の式に代入して  r_y を求める。


r_y  = d_yt + o_y

この  r_y (mininum, maximum) の範囲内にあれば Ray と円筒が交差する。

最後に cap 部分だけど、これもそんなに難しくはない。Ray の  y minimum, maximum になる場所での  x, z が半径 1 の円の中に収まっているかを見るだけ。

 y minimum, maximum になる場所を知るには Ray の  y 成分のみに注目すれば、


\eqalign{
  {d_y}t + o_y &= maximum \cr
  {d_y}t &= maximum - o_y \cr
  t &= \frac{maximum - o_y}{d_y}
}

で求まる。ここでは  maximum のみだけれど、  minimum についても同様。こうして求めた  t を Ray の式にあてはめればその時点での  x, z が求まるので、


  x^2 + z^2 <= 1

で判定できる。

Cone

Cone もほとんど Cylinder と変わらない計算で求まるから、ここでは側面との交点計算についてのみ書いておく。Cylinder の場合、側面は円筒で  y に関係なく常に半径は 1 だった。Cone になるとそうではなくなり、半径は  y になる。そのため、


\eqalign{
  {r_x}^2 + {r_z}^2 &= {r_y}^2\cr
  {(d_xt + o_x)}^2 + {(d_zt + o_z)}^2 &= {(d_yt + o_y)}^2\cr
  {({d_x}^2)}t^2 + 2d_xo_xt + {o_x}^2 + {d_z}^2{t}^2 + 2 d_zo_zt + {o_z}^2 &= {d_y}^2{t}^2 + 2 d_yo_yt + {o_y}^2\cr
  {({d_x}^2 - {d_y}^2 + {d_z}^2)}{t}^2 + 2(d_xo_x - d_yo_y + d_zo_z)t + {o_x}^2 - {o_y}^2 + {o_z}^2 &= 0
}

になり、


\eqalign{
  a &= {d_x}^2 - {d_y}^2 + {d_z}^2 \cr
  b &= 2({d_x}{o_x} - {d_y}{o_y} + {d_z}{o_z}) \cr
  c &= {o_x}^2 - {o_y}^2 + {o_z}^2
}

として二次方程式を解く。 a = 0 のときは  bt + c = 0 の一次方程式として


  t = -c / b

で ok。

結果

f:id:mtXTJocj:20210522145110p:plain

f:id:mtXTJocj:20210529140109p:plain

Cubes

Chapter 12

結果

本の通りにやっただけで、書くことがない。なので今回は結果だけ。

f:id:mtXTJocj:20210503140151p:plain

本の中にある画像を出力するための設定とかわかるとうれしいんですけどね。

f:id:mtXTJocj:20210503001650p:plain

Reflection and Refraction

Chapter 11

反射、屈折。

Reflection

今まではある Shape 上の 1 点の色は ambient, diffuse, および specular で決まっていた。これに reflection (反射)の分が加わる。

ambient を除けば、diffuse は Light の方向、specular は Light の反射方向(と視線ベクトル)で決まる。これに対して、新たに加わる reflection は視線(Ray)の反射方向を用いる。

    fn reflected_color(
        &self,
        is: &IntersectionState,
        remaining: usize,
    ) -> Color {
        if is.object.material().reflective == 0.0 {
            return Color::BLACK;
        }
        if remaining <= 0 {
            return Color::BLACK;
        }

        let reflect_ray = Ray::new(is.over_point.clone(), is.reflectv.clone());
        let color = self.color_at(&reflect_ray, remaining - 1);

        &color * is.object.material().reflective
    }

ここで、反射先でまた反射となる場合に再帰が無限に深くなってしまう可能性があるため、一定の深さで再帰を打ち切っている。

Transparency and Refraction

透明な Shape に対して更に refection (屈折)の分を加える。

refraction の場合も基本的な考え方は reflection と変わらない。ただ、reflection は視線の反射方向が Shape によらず視線と法線ベクトルから簡単に求まるのに対して、refraction は Shape ごとに異なる屈折率(refractive index)を用いて計算しなくてはならない。

しかも、単一の Shape ではなく、視線が出ていく Shape と次に視線が入る Shape の 2 つを考慮するため、reflection よりも大分複雑になる。

        let n_ratio = is.n1 / is.n2;
        let cos_i = is.eyev.dot(&is.normalv);
        let sin2_t = n_ratio * n_ratio * (1.0 - cos_i * cos_i);
        if sin2_t > 1.0 {
            return Color::BLACK;
        }

        let cos_t = (1.0 - sin2_t).sqrt();
        let direction =
            &(&is.normalv * (n_ratio * cos_i - cos_t)) - &(&is.eyev * n_ratio);

それ以外は reflection と変わらないため、方向ベクトルを求める部分を除けば reflected_color() とほぼ同じになる。

Fresnel Effect

反射と屈折の割合を近似的に計算する。本のままなので書くことがない。

結果

f:id:mtXTJocj:20210429181957p:plain

本の中で図示されている画像と同じ条件でレンダリングしたかったのだけれど、それが書かれないなかったので、てきとーにやった。 なので、これで正しくレンダリングできているかがいまいちはっきりしない。個別の単体テストは通っているので、大丈夫だと思うことにする。

f:id:mtXTJocj:20210429153207p:plain

Patterns

Chapter 10

まあ Procedual Texture のことだよね。

Pattern trait

基本的には Shape trait の時と変わらない。

pub trait Pattern: Debug {
    fn transform(&self) -> &Transform;
    fn transform_mut(&mut self) -> &mut Transform;

    fn pattern_at(&self, p: &Point3D) -> Color;
    fn pattern_at_shape(&self, shape: &dyn Shape, p: &Point3D) -> Color {
        let local_p = shape.transform().inv() * p;
        let pattern_p = self.transform().inv() * &local_p;
        self.pattern_at(&pattern_p)
    }
}

この Pattern trait を実装した。

pattern_at() が Pattern の計算を行う場所で、Pattern 座標系内の点 p における Color を返すように実装する。それ以外は Shape のローカル座標系と Pattern 座標系間の変換がらみ。

Pattern を所有するのは Material。複数の Material で 1 つの Pattern を共有しても問題ないから Rc を使うことも考えたけれど、とりあえずは Shape のときと同じ Box にしておいた。なので、Box を使っていることに根拠はなく、あとで Rc なりに変更するかもしれない。
というか、複数の Shape で Material を共有してしまっても良いよネ。

pub struct Material {
    pub color: Color,
    pub ambient: FLOAT,
    pub diffuse: FLOAT,
    pub specular: FLOAT,
    pub shininess: FLOAT,
    pub reflective: FLOAT,
    pub transparency: FLOAT,
    pub refractive_index: FLOAT,
    /// パターン。None の場合は使用しない。
    pattern: Option<Box<dyn Pattern>>,

Pattern は使用しない場合もあるため、Option にしている。

あとは各パターンに応じた色の計算を実装してやるだけでほぼ完了。

f:id:mtXTJocj:20201228175736p:plain

結果

f:id:mtXTJocj:20201229173046p:plain

Planes

Chapter 9

Shape trait

Plane の導入にあたって、今まで Sphere を使用していた場所で Plane や他の形状も使えるようにする。そのため trait Shape を用意して各形状で実装する。

pub trait Shape: Debug {
}

以降、Sphere を扱っていた場所を Shape trait object を使うように変更する。なお、trait object なので、コンパイル時に object size が決まらない。というわけで World のように Shape を所有する個所では直接所有する代わりに Box<dyn Shape> を所有する。

pub struct World {
    /// ライト
    lights: Vec<Light>,
    /// オブジェクト
    shapes: Vec<Box<dyn Shape>>,
}

また、Shape やその参照を所有する struct が Debug print できるように Shape は Debug trait を継承するようにした。そのため、Shape trait を実装する struct は必ず Debug trait も実装しなくてはならなくなったが、他に良いやり方がわからなかったのでとりあえずはこのようにしておく。

Shape trait の interface

基本的には本で書かれていた interface をそのまま定義するだけで問題ない。

pub trait Shape: Debug {
    fn transform(&self) -> &Transform;
    fn transform_mut(&mut self) -> &mut Transform;

    fn material(&self) -> &Material;
    fn material_mut(&mut self) -> &mut Material;

    fn intersect(&self, r: &Ray) -> Vec<Intersection> {
        let local_ray = self.transform().inv() * r;
        self.local_intersect(&local_ray)
    }
    fn local_intersect(&self, r: &Ray) -> Vec<Intersection>;

    fn normal_at(&self, p: &Point3D) -> Vector3D {
        let p_in_local = self.transform().inv() * p;
        let n = self.local_normal_at(&p_in_local);

        self.transform().apply_to_normal(&n)
    }
    fn local_normal_at(&self, p: &Point3D) -> Vector3D;
}

ここで intersect() と normal_at() はそれぞれ World-Local 変換をしてから local_intersect() と local_normal_at() を呼ぶという処理で、全ての Shape に共通していることから Shape で default の実装をする。各 Shape 固有の実装は local_*() の側で行う。

なので直接呼ばない local_*() は本来公開する必要はない。というか可能なら公開したくないのだけれど、trait の 一部メソッドのみを非公開にすることが可能なのか良くわからず現状は全て公開になってしまっている。

test の変更

Sphere から Shape への変更に伴い、一部テストを変更する必要がでてきた。

Intersection は hit した Shape の参照をメンバに持つが、それを確認するテストでは Shape 間の同一性を調べる必要が出でくる。以前はPartialEq を実装した Sphere のみであったため、単純な == による比較が使えたが、trait object となると同じ手法は使えない。

代わりにアドレスによる比較を行う。要は

        assert!(std::ptr::eq(&s as &dyn Shape, i.object));

のようにすれば万事 ok ... のはずだったのだけれど、うまくいかないケースがあった。アドレスは一致しているにもかかわらず、std::ptr::eq() が false を返すことがある。

どうも std::ptr::eq() は単純なアドレス比較以上のことを行っているみたいで、false を返す原因は使用する codegen Unit が異なる場合には同一アドレスの trait object でも vtable が異なる可能性があるためらしい。
std::ptr::eq() は trait object の fat pointer に対しては vtable のアドレス比較まで行うらしく、たとえ同一オブジェクトであったとしても、異なる module 内で trait object にした場合には std::ptr::eq() で比較すると false になり得る。

今回は vtable のアドレスの違いを気にする必要はないため fat pointer を regular pointer にキャストしてから比較することで、この問題を回避した。

        assert!(std::ptr::eq(
            xs[0].object as *const _ as *const (),
            &s as &dyn Shape as *const _ as *const (),
        ));

ところで、object のアドレスは一致するけれど vtable のアドレスが違うことを識別したいケースというのはどのくらいあるものなのだろう。ちょっと想像がつかない。

Plane

こちらが本題になるはずなのだけれど、特に書くことがない。

Plane は y 軸正方向を向いた xz 平面なので、local_intersect() では Ray の y 成分のみを考慮すれば良い。direction が平面と平行(direction の y 成分がほぼ 0)な場合は hit しない。

f:id:mtXTJocj:20201227125024p:plain

結果

f:id:mtXTJocj:20201227130930p:plain

前の見た目はほとんど変わっていないけれど、背景の壁が Sphere から Plane になっている。

Shadows

Chapter 8

この章で扱う内容は今までの道具立てを利用して実現できることから、そんなに書くことがない。

FLOAT 型の導入

今までは浮動小数点に f32 を用いてきた。だがこの章で self intersection を回避するために intersection 位置から法線ベクトル方向に少しだけ offset して shadow ray を生成するところで問題が発生した。
この本で acne と呼んでいる現象が offset を追加しても消えない。

原因は f32 の精度が足りないことで、これを f64 にするだけで問題は解決する。とはいえ、壁として用いている Sphere はかなり極端な scaling をしているためにこのようになってしまっただけで、そうでなければ f32 でも充分使える。

というわけで FLOAT という新しい型を導入して f32 と f64 を切り替えられるようにする。

pub type FLOAT = f64;

本来ならビルド時にオプションで切り替えられるようにすべきではあるが、とりあえずはこのようにしておく。

結果

f:id:mtXTJocj:20201222201702p:plain

f:id:mtXTJocj:20201222121134p:plain