レイトレーシング(8): 輝度推定による画像生成
前回、画像生成プログラムの大枠を作成した。今回はフォトンマップから各画素の輝度を 計算して画像生成までやってみよう。
各画素の処理
前回の記事でtraceScreenとした処理の実際の中身を考える。フォトンマッピング法は
光線追跡法をベースにして、（大雑把に言うと）拡散反射光を求めるのにフォトンマップを
使うやりかただ。光源から届く光の量を、光源への光線追跡により求めるかフォトンマップ
から推定するかの違いだ（と思う）。よって、処理の概略は次のとおり。
- 視線と物体との交点を求める。
- 交点における輝度をフォトンマップから推定する。
1の交点を求めるのは、以前にフォトンマップを作成するときに作ったのでそれをそのまま 流用すれば良い。問題は2の輝度の推定だ。
放射輝度推定
1) 交点に近い方からn個のフォトンを抽出する
kdtライブラリはkNearestという関数でn個の近傍点を抽出することができる。
まさにこれを使えば良い。抽出には中心点(交点)を他のフォトン情報と同じ
型で渡さないといけないので交点からダミーのフォトンを作って渡す。pは交点。
ps = kNearest pmap nPhoton $ photonDummy p photonDummy :: Position3 -> PhotonInfo photonDummy p = PhotonInfo Red p ex3
2) 推定に不要なフォトンを除く(フォトンの入射方向と交点の法線が反対側とか)
このpsから不要なものを除くが、ここでは交点の平面に対し裏側から入射したフォトンを
除外することにしよう。要するに、交点の法線とフォトンの入射方向の内積が負なら省くと
いうことだ。nは交点の法線ベクトル。
ps' = filter (isValidPhoton n) ps isValidPhoton :: Direction3 -> PhotonInfo -> Bool isValidPhoton n pi = n <.> (photonDir pi) > 0 photonDir :: PhotonInfo -> Direction3 photonDir (PhotonInfo _ _ d) = d
ここで、「不要なフォトンを除外したあとでn個のフォトンを使って推定すべきでは？」 という指摘があるだろう。その通りだと思うが、ちょっと面倒くさいので今回は無視だ。
3,4) 交点と各フォトンの距離を求める+一番遠いフォトンの距離を選ぶ
これは簡単だ。フォトンから位置を取り出し、交点との差を取って大きさを求めればいい。
rs = map (\x -> norm ((photonPos x) - p)) ps' rmax = maximum rs
haskellだと、リスト処理の関数を使えるので簡単だな。
5,6) 各フォトンによる照度を計算して集計する+交点の輝度を計算する
このうち、このステップで使うのはの部分だ。
は、今の所は完全拡散反射しか扱わないことにしたのでとても簡単だ。
これは、完全拡散反射が微小平面で全ての方向へ均等に反射するため半球全体で積分したもの
で割って「輝度」にするためだ。
各フォトンから得られた全照度にこれをかければいいので
の外に出せる。
あとは
の集計だ。各フォトンに一個分のエネルギーを掛けて
合計してやればよい。
(Qiitaの方でshocker-0x15さんにご指摘いただきました。ありがとうございました。修正しました。)
rad = foldl (+) radiance0 $ map (photonInfoToRadiance pw) ps' photonInfoToRadiance :: Double -> PhotonInfo -> Radiance photonInfoToRadiance pw (PhotonInfo Red _ _) = Radiance pw 0 0 photonInfoToRadiance pw (PhotonInfo Green _ _) = Radiance 0 pw 0 photonInfoToRadiance pw (PhotonInfo Blue _ _) = Radiance 0 0 pw
各フォトンは周波数情報しか持たないので、ここでは赤緑青それぞれに一個のエネルギーを
かけた「輝度」を返しそれを足し合わせている。radiance0は初期値で「黒」だ。
得られた集計値にを掛け、最後に推定に使った全フォトンが含まれる円の面積で
割って完了だ。円の面積は先の求めた最大半径(
rmax)を使うのだ。なお
を掛けるところは後々の拡張も考え、
brdfという関数に切り出す。
radiande = (1.0 / (pi * rmax * rmax)) *> (brdf m rad) brdf :: Material -> Radiance -> Radiance brdf m rad = (1.0 / pi2) *> ((diffSpec m) <**> rad) diffSpec :: Material -> Color diffSpec (Material d) = d
やっとできた。
サンプルシーンの描画
プログラムができたので、早速実行してみよう。以前作った二種類のフォトンマップで 画像生成した結果は次のとおり。比較のため、使ったフォトンマップ（を可視化したもの)と 古典的光線追跡法で描画したものも並べてみる。
・・・とても微妙な出来栄えだ・・・。今作っているのは「バージョン1」で、点光源だけ、 拡散反射のみ、フォトン追跡は反射を無視（相互拡散反射は計算しない）という条件がついて いるので、ある程度予想していたが、それにしても。。。ちなみに、この条件では古典的 光線追跡法は（たぶん）「完璧な画像」になるはずだ。いわゆる環境光成分も入れていないし。 だから、理想は一番右の画像にすることだ。
例2の生成画像から、ダメダメな点をいくつかピックアップしてみよう。
- 光源から放射される光にバラツキがあるので綺麗な同心円にならない(A)
- 床面の輝度が光が少ない壁面に滲んでいる(B)
- 球の影がくっきり出ていない、ほとんど影になっていない(C)
- 壁が均一の色になっていない(D)
A)について
有限個のフォトンしか使えないのである程度まばらになるのは仕方がないが、もう少し 均一にしたい。そういえば光源から放射されるフォトンは本当に均一なのか疑問だった。 また、層化サンプリングも有効かもしれない。現状は球状にばらまいているだけだが、 球面をいくつかの部分に分けてその部分ごとにランダムに放射すれば均一性が高まるだろう。 次回検証しよう。フォトンの均一性は後述のD)にも影響してくるはずだ。
B)について
交点に近いフォトンを集めて輝度を推定するため致し方ない結果ではある。これを軽減する 手法について例の本に記載があるが。。。今はバージョン1のため特に奥の壁にはフォトンが 到達しないようになっている。機能拡張をしてフォトン追跡で反射も扱うようにすれば、 奥の壁にもフォトンが届き、自ずとこの問題は目立たなくなるだろう。そもそも床面の フォトンを推定に使うことは「問題」ということはなく、それだけのフォトンが近くに 到達しているのだから、計算に入れるのは「妥当」と考えるべきだろう。 この懸案はバージョン2まで保留する。
C)について
一番ショックなのがこれ。影がほとんどない。フォトンマップはそのアルゴリズム上、 フォトンがないと推定できない。このように「まったくフォトンが届かない場所」は 正しく描画できなくて当たり前かもしれない。相互拡散反射に対応したら緩和されるかも しれないが、一方、例の本にある「円錐フィルタ」を適用してその効果を見てみようと おもう。これも次回以降だ。
D)について
これはフォトンマッピング法における「誤差」成分なので、フォトン数を増やしたり、 輝度推定に使うフォトン数を増やせば低減されるものだ。ただフォトン数を増やすにも 限界があるし輝度推定時にフォトンを増やすと推定に使うフォトンの分布範囲が広がって しまうため余計なフォトンを推定に使ってしまい具合が悪い（Cの問題もその一つと言える)。 上の画像は推定に200個のフォトンを使ったが、下記に50,100,300も示そう。
個数を増やすとにじみはかなり消えてくる一方、影はなんだかわからなくなっている。 逆にn=50にすると影は良くなってきているが壁の滲みはかなりひどい。 推定に使うフォトン数の最適値を考えるのは難しいのでとりあえず保留しよう。 まずは他の部分の改良だ。
まとめ
今回、曲がりなりにも画像生成までこぎつけ、古典的光線追跡法と比較することができた。 思ったほど綺麗な画像にならなかったが、それも改善のしがいがあると捉えよう。 次回は光源からのフォトン生成と推定時のフィルタ導入だ。
（今回のソースはこちらから）
レイトレーシング(7): 光線追跡処理の大枠
忙しさにかまけて更新をサボってしまった。。。気を取り直して再開しよう。 前回まででフォトンマップが出来上がったので、今回からはそれを使った光線追跡の プログラムに取り組む。
処理の大枠
フォトンマップを使う以外は古典的な光線追跡のアルゴリズムである。 一番大枠の処理は次の通り。
- 仮想カメラを用意（定義）する（カメラ位置、スクリーン、解像度など）。
- フォトンマップを読み込む。
- スクリーンの各画素について光線追跡を行い輝度を得る。
- 輝度値を画像ファイルに書き出す。
3が本プログラムのメインとなるが、今回は外堀を埋めよう。1、2、4を定義する。 3も画素をループするところまで。
仮想カメラの用意
1について必要なことを書き出してみる。流儀などがあるだろうが、自分が 作るときには以下の要素を用いている。
- 視線
- 視点
- 視線方向: スクリーンの中心方向
- 垂直方向: スクリーンの「上」方向を決めるための情報
- スクリーン
- フォーカス: スクリーンまでの距離としている
- 解像度(x, y)
なお、スクリーンの中心点は視点から視線方向にフォーカス分の距離を進んだ地点とし、 そこから縦横±1.0の領域と固定している。この正方形を解像度の設定値で細かく 分割してそれぞれに光線追跡するわけだ。なお、スクリーンは視線方向に垂直な平面である。
これらの情報は、ちゃんとしたプログラムにするときには設定ファイルなどから読み込んで いろいろ変更できるようにするのだろうが、とにかく「フォトンマッピング」法の効果を 見たいので決め打ち＝ソース中に直書きしてしまう。具体的には前回示したフォトン マップの視覚化と同じ設定を使う(あとで画像を比較できるし)。今回使う設定は以下だ。
eyepos = initPos 0 2 0 -- 視点 eyedir = ez3 -- 視線方向 (0, 0, 1)=Z軸 upper = ey3 -- 上方向を表すベクトル (0, 1, 0)=Y軸 focus = 1.0 :: Double -- フォーカス(スクリーンまでの距離) xres = 256 :: Int -- 解像度(横方向) yres = 256 :: Int -- 解像度(縦方向)
ただ、これだけでは光線追跡処理には不十分なので、いくつかの定義を追加する。 最初の画素の位置や次の画素の場所（縦横方向）などを決めてやらないといけないから あらかじめ計算しておく。
stepx = 2.0 / fromIntegral xres :: Double -- 画素の大きさ(横) stepy = 2.0 / fromIntegral yres :: Double -- 画素の大きさ(縦) eex = ex3 -- スクリーンの横方向を表す方向ベクトル(大きさ1) eey = negate ey3 -- スクリーンの縦方向を表す方向ベクトル(大きさ1) origin = focus *> eyedir + ((-1.0 + 0.5 * stepx) *> eex) - (( 1.0 - 0.5 * stepy) *> eey) -- スクリーンの左上の点(最初の画素)
これらを使えば、位置(X, Y)の画素pointは
point = origin + X * stepx * eex + Y * stepy * eey
で表せるわけだ。
フォトンマップの読み込み
さてほとんど忘れかけていたが、前回作ったフォトンマップのフォーマットを思い出してみよう。
次のような形だった。フォトンキャッシュ情報はPhotonCache型である。
フォトン数(改行) フォトン1個のエネルギー(改行) 1番目のフォトンキャッシュ情報(改行) 2番目のフォトンキャッシュ情報(改行) : n番目のフォトンキャッシュ情報(改行)
よって読み込みはこの逆をしてやればよい。最終的なデータ構造にはk-d treeを使う。 これは例の本(フォトンマッピング)にk-d treeが良いと書いてあったのでそのまま採用 しただけだ。幸い、haskellにはk-d treeのライブラリがいくつかあるので自作する必要は ない。実はここでちょっとつまづいた。ライブラリによって、処理性能（速度）が段違い なのだ。ここでは最終的に採用した kdtを元に進める。最初使ったものは 遅すぎてとても実用にならなかったので（多分、処理速度が数十倍違う）。インストールは cabalから。
cabal install kdt
kdtで使うため、読み込んだデータを変換してPhotonInfo型のListにしよう。
PhotonInfo型の定義、変換用関数は以下のとおり（PhotonCache型を使えば
いいのだが、フォトンの入射方向を逆転させる必要があるのでこれはあとで取り組む）。
また、k-d treeのデータ構造へ変換するときに必要になるので、フォトンの位置を取り出す
関数も定義しておく。
data PhotonInfo = PhotonInfo Wavelength Position3 Direction3 deriving (Show, Eq) convertToInfo :: PhotonCache -> PhotonInfo convertToInfo (wl, (rp, rd)) = PhotonInfo wl rp (negate rd) infoToPointList :: PhotonInfo -> [Double] infoToPointList (PhotonInfo _ (Vector3 x y z) _) = [x, y, z]
あとは全フォトン情報を読み出してk-d treeに変換すれば良い。なお、フォトンマップの データは標準入力で与えることにする。
各画素の輝度値を求める
具体的な処理は次回以降に回すとして、今回は大枠なので適当な輝度値をでっち上げ、
次の処理でファイルに書き出せるようにしよう。入力は各画素の位置、出力は輝度値とする。
輝度値は、各周波数の値を持った型を定義することにしよう。Radianceとする。
data Radiance = Radiance Double Double Double deriving (Read, Show)
輝度値も光線追跡の最中に頻繁に演算する。なのでVectorとほとんど同じように
加減などの演算関数を定義することになる。が、まずは型だけ。各画素の位置はxとyの
組みで表そう。そうすると、レイトレーシングのメイン処理は大体次のような形に
なるだろう。全画素の組はあらかじめ作っておき、mapで各組を輝度値に
変換する。
scrmap = [(y, x) | y <- [0..(yres - 1)], x <- [0..(xres - 1)]] : image = map (traceScreen ??) srcmap traceScreen :: <<必要な引数>> -> (Int, Int) -> Radiance traceScreen ?? scrmap = ... :
ひとまず、RGB各周波数の輝度値が 0.05 のRadianceのリストを返すようにしておこう。 あとでわかるが、0.05にしておくと灰色になるはずだ。
輝度値の書き出し
輝度値を画像ファイルへ書き出すため、各周波数成分を256段階の整数に変換したい。 次にような関数を用意すれば良いだろう。
clip :: Double clip = 0.1 gamma = 1.0 / 2.2 rgbmax = 255.0 radianceToRgb :: Double -> Int radianceToRgb d = floor (r * rgbmax) where d' = d / clip r = (if d' > 1.0 then 1.0 else d') ** gamma
なお、clipは輝度値と最大値"255"(=rgbmax)の対応付けのための閾値。
0.1という設定は、0.1[W]以上のエネルギーを持っていたら最大値とするということ。
この値は適当なので、出来上がる画像の明暗を見て適当に変更するつもりだ。gammaは
ガンマ補正のための設定値で、一般的(?)な を用いる。
出力する画像データだが、単純なPPM形式を用いよう。フォーマットはググればすぐ わかる。これを標準出力へ吐き出すようにする。以上を踏まえると、書き出し処理の 関数は次のようなものになるだろう。
outputImage :: [Radiance] -> IO () outputImage rs = do putStrLn "P3" putStrLn "## test" putStrLn (show xres ++ " " ++ show yres) putStrLn "255" forM_ rs $ \i -> do putStrLn convertOneCell i
各画素の輝度値のリストを引数にして書きだすだけだ。なおconvertOneCellは
一画素の結果を文字列に直す関数。
まとめ
上記をまとめ、ソースにしてみた。
(RT.hs) module Main where import System.IO import Control.Monad import Data.KdTree.Static import NumericPrelude import Ray.Algebra import Ray.Geometry import Ray.Physics import Ray.Optics -- PARAMETERS -- -- for camera eyepos = initPos 0 2 0 eyedir = ez3 upper = ey3 focus = 1.0 :: Double xres = 256 :: Int yres = 256 :: Int stepx = 2.0 / fromIntegral xres :: Double stepy = 2.0 / fromIntegral yres :: Double eex = ex3 eey = negate ey3 origin = focus *> eyedir + ((-1.0 + 0.5 * stepx) *> eex) - (( 1.0 - 0.5 * stepy) *> eey) scrmap = [(y, x) | y <- [0..(yres - 1)], x <- [0..(xres - 1)]] -- FUNCTIONS -- main :: IO () main = do (power, photonmap) <- readMap let image = map traceScreen scrmap outputImage image -- readMap :: IO (Double, KdTree Double PhotonInfo) readMap = do np' <- getLine pw' <- getLine let np = read np' :: Int let pw = read pw' :: Double pcs <- forM ([1..np]) $ \i -> do l <- getLine return $ (read l :: PhotonCache) let pmap = build infoToPointList (map convertToInfo pcs) return (pw, pmap) -- traceScreen :: (Int, Int) -> Radiance traceScreen (y, x) = Radiance 0.05 0.05 0.05 -- outputImage :: [Radiance] -> IO () outputImage rs = do putStrLn "P3" putStrLn "## test" putStrLn (show xres ++ " " ++ show yres) putStrLn "255" forM_ rs $ \i -> do putStrLn $ convertOneCell i convertOneCell :: Radiance -> [Char] convertOneCell (Radiance r g b) = (show $ radianceToRgb r) ++ " " ++ (show $ radianceToRgb g) ++ " " ++ (show $ radianceToRgb b)
(Ray/Physics.hs) data Wavelength = Rad | Green | Blue deriving (Show, Read, Enum, Eq)
(Ray/Optics.hs) data Radiance = Radiance Double Double Double deriving (Read, Show) clip = 0.1 :: Double gamma = 1.0 / 2.2 rgbmax = 255.0 radianceToRgb :: Double -> Int radianceToRgb d = floor (r * rgbmax) where d' = d / clip r = (if d' > 1.0 then 1.0 else d') ** gamma type PhotonCache = (Wavelength, Ray) data PhotonInfo = PhotonInfo Wavelength Position3 Direction3 deriving (Show, Eq) infoToPointList :: PhotonInfo -> [Double] infoToPointList (PhotonInfo _ (Vector3 x y z) _) = [x, y, z] convertToInfo :: PhotonCache -> PhotonInfo convertToInfo (wl, (rp, rd)) = PhotonInfo wl rp (negate rd)
コンパイルして実行してみる。
$ ghc -o rt0 RT.hs $ rt0 < photonmap > out1.ppm $ convert out1.ppm out1.png <= ImageMagickを使う
結果は以下のとおり。当たり前だが面白くもなんともない・・・。
ということで、次回は光線追跡の処理を考えていこう。
レイトレーシング(6): やっとフォトン追跡
シーン作成
前回までで、物体を定義するための準備ができた。それを用いて簡単なサンプルシーンを 作ってみる。閉じた空間でないと、せっかく放射したフォトンが無限遠に行ってしまって無駄に なるので、箱型の空間を用意しよう。これで平面が6つ。それだけではつまらないので、 1つだけ球面を置こう。その定義を以下に示す。
m_ball = Material 0.8 0.3 0.3 m_wall = Material 0.8 0.8 0.8 m_ceil = Material 0.4 0.2 0.02 m_flor = Material 0.8 0.6 0.4 wall_bt = initObject (Plain ey3 0) m_flor -- bottom wall_tp = initObject (Plain (negate ey3) 4) m_ceil wall_rt = initObject (Plain (negate ex3) 2) m_wall wall_lt = initObject (Plain ex3 2) m_wall wall_bk = initObject (Plain ez3 1) m_wall wall_ft = initObject (Plain (negate ez3) 5) m_wall ball1 = initObject (Sphere (initPos 1 0.8 3) 0.8) m_ball objs = [wall_bt, wall_tp, wall_rt, wall_lt, wall_bk, wall_ft, ball1]
横幅と高さが4[m]、奥行き6[m]の部屋に直径1.6[m]の球が置いてあるという想定だ。 壁は白、床はオレンジ、球は少し薄めの赤とした。まだ描画していないので色は よくわからないけど。
最終的にobjsが物体のリストである。これを使ってフォトンを追跡しよう。
フォトンの追跡
前回
示したtracePhotons関数を定義しよう。といっても、物体のリストとフォトンのリストを
引数にするだけだ。
tracePhotons :: [Object] -> [Photon] -> IO [PhotonCache] tracePhotons objs photons = do pcs <- mapM (tracePhoton objs) photons return $ concat pcs
うーん、何のことはない。フォトンの数だけtracePhoton(複数形でないことに注意)を
呼び出して、結果をフラットなリストにして返しているだけだ。tracePhotonが具体的な
追跡処理なのだが、実行することは次のようになるだろう。
これをコードにしたら次のようになった。
tracePhoton :: [Object] -> Photon -> IO [PhotonCache] tracePhoton os (wl, r) = do let iss = filter (\x -> fst x > nearly0) (concat $ map (calcDistance r) os) (t, s) = head $ sortBy (comparing fst) iss return [(wl, initRay (target t r) (getDir r))] calcDistance :: Ray -> Object -> [(Double, Object)] calcDistance r o@(Object s m) = zip ts (replicate (length ts) o) where ts = distance r s
しつこいようだが"バージョン1"ではフォトンの反射・屈折は扱わないので、上記コードは
若干無駄である。例えばtracePhotonの返値はこの場合1個しかないので
IO [PhotonCache]とリストにする必要はない。しかし今後、反射・屈折に対応したら
1個のフォトンの追跡で複数のフォトンキャッシュ情報が得られるようになる（ハズ）だから
今のうちにその対応をしておこうということだ。また内側で使われているcalcDistance
についても、今必要なのは交点の位置ベクトルだけなのでtarget t rを呼んで距離を
位置ベクトルに変換するだけでよいハズだ。しかし、反射・屈折を扱うには交点での材質の
情報が必要だ。また反射・屈折方向を求めるためには衝突した物体の「形」も要る。
ということで、あえて衝突した位置までの距離と衝突した物体をセットにして返している
わけだ。
フォトンキャッシュ情報の書き出し
さあ、一連の処理の最後だ。各フォトンを追跡して得られたフォトンキャッシュ情報を 標準出力に書き出そう。これをファイルに書き込めばフォトンマップのデータファイルとなる。 ファイル形式は、
フォトン数(改行) フォトン1個のエネルギー(改行) 1番目のフォトンキャッシュ情報(改行) 2番目のフォトンキャッシュ情報(改行) : n番目のフォトンキャッシュ情報(改行)
とする。以前示したメインルーチンより少しシェイプアップして次のようになった。
putStrLn $ show nphoton putStrLn $ show power forM_ (photoncaches) $ \i -> do putStrLn $ show i
やっとできた!これまでのプログラムをまとめ、コンパイルして実行だ。 結果はかなり大きなファイルになる。10万個のフォトンで約130MBぐらい。そうか、 反射・屈折なしでこれぐらいということは、それを実装すると簡単に1GBに到達する。 まあ、ファイル形式の再考はあとにしよう。
実行結果
ここまでくると、フォトンマップがどんな結果になったのかどうしても見てみたい。 もちろん、ちゃんと動いているのか気になる所であるし。フォトンキャッシュ情報にある 交点座標を、適当な二次元平面に投影させてフォトンの色も踏まえて描き出してみた結果が 以下だ。
おお、なんとなくそれっぽいかも。ちゃんと部屋になって、右下に球も見える。 黒いところは球の影でフォトンが届かないところだ。
物は試しに、複数の光源に変えて実行してみる。天井の中央ではなく、左右奥に白と橙色の 2つの光源だ。ただ、全体の光強度は同じ、フォトン数も同じだ。
なかなかよい。とりあえず、フォトンマップは完成した、ということにしよう。
メモリ使用量の改善
・・・いや、終わらなかった・・・。
フォトン数を100万個にして処理時間を計測中、たまたま眺めていたtopコマンドの 表示に驚いた。処理中の最大メモリ使用量が約1.5GBだったのだ! 改めて作成したメインルーチンを見てみよう。
main = do (power, photons) <- generatePhotons nphoton lgts photoncaches <- tracePhotons objs photons putStrLn $ show nphoton putStrLn $ show power forM_ (photoncaches) $ \i -> do putStrLn $ show i generatePhotons :: Int -> [Light] -> IO (Double, [Photon]) generatePhotons nphoton lights = do let power = (sum $ map flux lights) / (fromIntegral nphoton) ns = map (calcN power) lights photons <- mapM (mapM generatePhoton) (zipWith replicate ns lights) return $ (power, concat photons)
フォトン数に応じていくつも巨大なリストを作っている。フォトンを生成するのに まず、n個の光源のリストを生成し、そのリスト全体に処理をかけている。 さらにそのフォトンのリストを処理してフォトンキャッシュのリストを生成、 最後にリスト全体を表示している。とにかく大きなリストをいくつも作っている。 全フォトンを並列に処理しているのでフォトン数が増えるほどメモリを食うのも仕方なしか。
ここで、フォトンを1個だけ追跡する処理を改めて考えてみよう。
最後にフォトンキャッシュを出力してしまえば、何も情報を保存しておく必要が ないので、並列ではなく直列に1個ずつ処理を終わらせていけばいいということか。
ような「二重ループ」にすればよいだろう。コードにしてみたのが以下。
main :: IO () main = do putStrLn $ show nphoton let power = (sum $ map flux lgts) / (fromIntegral nphoton) ns = map (calcN power) lgts putStrLn $ show power zipWithM_ outputPhotonCaches ns lgts calcN :: Double -> Light -> Int calcN power light = round (flux light / power) outputPhotonCaches :: Int -> Light -> IO () outputPhotonCaches n lgt = mapM_ (outputPhotonCache lgt) [1..n] outputPhotonCache :: Light -> Int -> IO () outputPhotonCache lgt _ = generatePhoton lgt >>= tracePhoton objs >>= mapM_ (putStrLn.show)
mainの中で求めているnsが、各光源から放出されるフォトン数n(i)のリスト。
最終行のzipWithM_でそれぞれの光源について処理している(外側のループ)。
outputPhotonCahcesではn(i)回フォトンキャッシュを計算して出力するように
している(内側のループ)。ちなみに2つ目の引数は使わないので捨てている。
このように「同じ処理をn回実行する」というのはどう書くのがスマート
なんだろうね。
`outputPhotnCache'(単数形)が実際に一個のフォトンを生成から追跡、最後の 書き出しまでやっている関数だ。
これをコンパイルして実行してみると、出力結果は同じだが実行時のメモリ使用量が 「数十kB」まで激減した!アルゴリズムは大事だなぁ。（というか最初考えた アルゴリズムがあまりにもイケていなかっただけか・・・）
ようやく完成した。
レイトレーシング(5): 物体の定義
前回までで、光源から放射されるフォトンが生成できたので、次はそれを 追跡してフォトンマップを作ることになるが、そのためには描かれる「物体」を 準備しないといけない。
「物体」定義
次に考える処理は、メインルーチン中の次の部分だ。
photoncaches <- tracePhotons objs photons
objsはシーン中の物体のリスト、photonsは前回生成したフォトンのリストである。
だから、tracePhotonsの処理を考える前にobjsがどういうものかを定義する必要がある。
本プログラムでは物体をObject型で定義しよう。物体に必要な情報は何だろうと考えると
「形」と「材質」だろう、と考えてみた。今後他にも出てくるかもしれないが、今は
その二つにしておく。
data Object = Object Shape Material
材質は木とかガラスとかで、情報としては物体の色や反射率、透明度などなどだ。 考え出すととても複雑な構造になるが、まず作ろうとしている"バージョン1"では 物体は「表面が拡散反射のみ」で「フォトンの追跡は反射を無視」としたのだった。 だからフォトンの追跡は「物体にぶつかったらその位置を記録して終わり」という ことになる。実はこれだとフォトンマップを作るためには材質として何の情報も要らない のだが、あとあとレイトレーシングで画像を作る段になったらさすがに色の情報が 必要なのでその分だけ定義しておこう。「拡散反射率」だ。
data Material = Material Double Double Double
3つのDoubleは赤緑青それぞれの波長での拡散反射率を表す。0から1の間の数値を
設定する。色や輝度をどのような型で表すかまだ決めていないが、その検討次第では
型コンストラクタの引数は型が変わるだろう。
さて、今大事なのは「形」の方だ。レイトレーシングでは、無限平面や球面、二次曲面、 ポリゴンなど様々な「形」が使われる。光源と同じだ。光源の定義では失敗したので、 形の定義では最初から多相性を意識して定義しよう。ただしひとまず無限平面と球面のみ 扱うことにする。（ちょっと先のことも考え、大きさのない"点"も入れておくが)
data Shape = Point Position3 | Plain Direction3 Double | Sphere Position3 Double
ここで、無限平面と球面は次の方程式を満たす三次元空間中の点
(位置ベクトル)の集合である。
ここで、は平面の法線ベクトル、
は平面の位置に
関係するパラメータ、
は球の中心座標、
は球の
半径だ。このあたりの詳しいところはその筋の文献などを参照のこと。たとえば、
- 作者: Fletcher Dunn,Ian Parberry,松田晃一
- 出版社/メーカー: オライリージャパン
- 発売日: 2008/10/04
- メディア: 大型本
- 購入: 21人 クリック: 141回
- この商品を含むブログ (40件) を見る
などに記載がある。上記のPlainとSphereの型コンストラクタの引数は、それぞれ
と
である。
「形」に必要な関数
フォトンの追跡は、まず物体と衝突する場所(交点)を求めることから始まる。
バージョン1では反射は考えないので交点計算がやることのすべてと言っていい。
Shape型は上記の通り方程式で記述できるので、交点を求めるのは容易だ。光線(Ray)
との連立方程式を解けばよい。交点を
とすると、
連立させて
を求めれば、位置
も解るわけだ。
よって、
を計算する関数
distanceを用意しよう。
distance :: Ray -> Shape -> [Double] -- Point distance r (Point p) = [] -- Plain distance (pos, dir) (Plain n d) | cos == 0 = [] | otherwise = [(d + n <.> pos) / (-cos)] where cos = n <.> dir -- Sphere distance (pos, dir) (Sphere c r) | t1 <= 0.0 = [] | t2 == 0.0 = [t0] | t1 > 0.0 = [t0 - t2, t0 + t2] where o = c - pos t0 = o <.> dir t1 = r * r - (square o - (t0 * t0)) t2 = sqrt t1
交点がない場合、複数の場合があるので、結果はのリストとする。
今後反射や屈折を考えたり、輝度計算をするときには交点での法線ベクトルを求める 必要が出てくる。今はいらないが簡単なので定義しておこう。
getNormal :: Position3 -> Shape -> Maybe Direction3 -- Point getNormal p (Point p') = Nothing -- Plain getNormal p (Plain n d) = Just n -- Sphere getNormal p (Sphere c r) = normalize (p - c)
「点」の場合法線ベクトルがないので、関数getNormalの結果をMaybe型に
している。
まとめ
レイトレーシング(4): フォトンの生成
光源の定義をやりなおす
前回の記事で、光源について型クラスLightを定義した上で
各種光源の型(例えば点光源PointLight)をそのインスタンスと
するようにした。その後、シーン情報内で光源のリストを作ろうと
した時に問題に気づいた。Haskellでは"同じ型"しかリストに
できない。新たに面光源を用意したとして、「光源リスト」に
点光源と面光源を混ぜることはできない。
オブジェクト指向でのスーパークラスの概念とごちゃまぜに
してしまったのだ。。。
改訂した定義は以下の通り。
data Light = PointLight Color Flux Position3 flux :: Light -> Flux flux (PointLight _ f _) = f generatePhoton :: Light -> IO Photon generatePhoton (PointLight c _ p) = do : (あとは同じ)
引数でPointLightをパターンマッチさせるところがHaskellらしい。
今後、平面光源や球光源を追加するときは、dataの定義や各関数の
パターンを増やしていけばよい。
フォトンの生成
前回の話に立ち返り、メイン処理の前半を詳細化していこう。メイン処理を次に 再掲する。
main = do photons <- generatePhotons nphoton lgts photoncaches <- tracePhotons objs photons a <- forM (photoncaches) $ \i -> do putStrLn $ show i return ()
このphotons <- generatePhotons nphoton lgtsの部分だ。
nphotonとlgtsはそれぞれ、追跡したいフォトン数と光源のリスト。
次のように定義しよう。
nphoton = 100000 :: Int pl1 = PointLight (initColor 1 1 1) 2.0 (initPos 0 3.95 2.5) -- 2W lgts = [pl1]
フォトン数は生成される画像の品質を左右するパラメータだ。とりあえず10万個としよう。
点光源1個で試してみる。この点光源の定義は、赤緑青が同じ比率＝白色光
(initColor 1 1 1のところ)であり、発光の強さ(光束、ざっくり言えばエネルギー）が
2[W]、光源の位置は x,y,z=0, 3.95, 2.5 [m]ということ。ちなみにSI単位系とし、
長さは[m]とする。強さが2[W]だとかなり暗いのではと思われがちだが・・・。
光の強さというか明るさについては、白熱電球だと100[W]とか60[W]とかがメジャーだ。
ちなみに白熱電球の場合、エネルギー変換効率は2.3〜2.6%ぐらいだそうで、ほとんどが
熱などになってしまうということだ。地球に優しくないとして製造中止になるわけだ。
一方、最近LED電球ではlumen(lm)を明るさを示すために使っているが
白熱電球の100[W]クラスは大体1520[lm]ぐらいらしい。それで、lumenをエネルギーという
意味でワットに換算するとき、波長555[nm]の光では、
1[W] = 683[lm]
だそうなので、波長の違いをとりあえず無視して、1520 ÷ 683 ≒ 2[W] としたのだ。 この換算が正しいかどうかよくわからないが。 ただし、波長が違うとこの数字が大きく異なるので、本当は無視するとダメかもしれない。
ではあらためてgeneratePhotonsの定義を考えよう。要は、複数の光源からn個の
フォトンが放出される、という状況を作り出したい。Lightの定義で、ある光源から1個の
フォトンを放出する関数は定義した(generatePhoton :: Light -> IO Photon)。
これを使って全光源から放出される数をnphotonにしたいのだ。そうすると次のことを
考えていくことになろう。
ということで、コードにしてみた。
generatePhotons :: Int -> [Light] -> IO (Double, [Photon]) generatePhotons nphoton lights = do let power = (sum $ map flux lights) / (fromIntegral nphoton) ns = map (calcN power) lights photons <- mapM (mapM generatePhoton) (zipWith replicate ns lights) return $ (power, concat photons) calcN :: Double -> Light -> Int calcN power light = round (flux light / power)
ここで、powerがフォトン一個のエネルギー(W)でありnsは各光源が放出すべき
フォトン数のリストだ。
ここまでくれば後段はその数だけフォトンの生成を繰り返せば良い。具体的には各光源を
放出するフォトンの数だけ並べたリストを作り、そのリストの要素毎にgeneratePhotonを
呼び出しているだけだ。これで欲しい数のフォトンが得られた。
生成されたフォトンの品質
さて、これで10万個のフォトンを生成したが、果たしてちゃんとランダムに生成できている のだろうか?放出される偏りはないのだろうか? 点光源なので、無作為に球状のどの方向にも偏りなく放出されて欲しい。もし完璧に偏り なければ、フォトン数を無限に近づけることで全フォトンの放出方向のベクトル（大きさは 1に正規化されている）を足し合わせれば「0ベクトル」になるだろう。既出のメインルーチンの 処理を「途中下車」して、生成されたフォトン群を出力させた。さらにそれを下記の プログラムで読み込み、全ベクトルを足し合わせたものを出力してみた。
main :: IO () main = do np <- getLine pw <- getLine let nphoton = read np :: Int power = read pw :: Double dat <- getContents pcs <- forM (lines dat) $ \i -> do return $ (read i :: Photon) let tv = sumVectors pcs putStrLn $ show tv sumVectors :: [Photon] -> Position3 sumVectors [] = o3 sumVectors (pc:pcs) = getDirVec pc + sumVectors pcs getDirVec :: Photon -> Position3 getDirVec (_, (_, d)) = d
5回実行し、その結果のベクトルの長さを計算したところ下表の通りとなった。 もちろん理想は長さゼロだ。
| #photon | min | max | avg | err/photon | time |
|---|---|---|---|---|---|
| 10000 | 80.5 | 129.7 | 101.3 | 1.01% | 0.62 s |
| 100000 | 47.2 | 373.4 | 216.1 | 0.22% | 5.70 s |
| 1000000 | 531.8 | 1232.8 | 792.3 | 0.079% | 57.18 s |
これが統計的にみて「十分に無作為」と言えるのかどうか正直微妙だ。 しかし今の所、何か改善できるのか、どうしたらできるのかよくわからないので とりあえず放置しよう。
一方、生成されるフォトンの波長毎の数は無作為か? 測ってみたところ下記のようになった。（フォトン数10万個で5回試行）
| light | red | green | blue |
|---|---|---|---|
| 白色(2W) | 33.3% | 33.3% | 33.3% |
| 白色(2W)+橙色(1W) | 44.4% | 33.3% | 22.2% |
橙色の光源は波長の割合を赤2:緑1:青0とした。 こちらは綺麗に無作為な生成ができていると言える。まあ、波長については 一様乱数で0-1の値を生成してそれを各色の割合で直接選択しているのだから うまくいって当たり前か。
まとめ
とりあえず、フォトンを生成することはできた。次回は各フォトンを追跡して フォトンマップを作ることにしよう。
※ これまでのソースはここ。 ただし、記事より先に進んでいる場合があるので注意。
レイトレーシング(3): フォトンマップ生成の大枠を考える（ついでに光源も）
以前に書いたように、フォトンマッピング法は(1)フォトンマップの生成、 (2)レイトレーシングの二段階で画像を生成する手法だ。だから、まずは フォトンマップを作らないと始まらない。今回からそのプログラムを作っていこう。
仕様とメインルーチン
まずは仕様から。「シーン」の情報が与えられているとして、次のような ステップで生成したらいいだろう。「シーン」とは描こうとしている三次元世界の ことで、情報とは「ここにこういう球がある」とか「ここにスポットライトを配置 して」とかいう設定のこととする。
- 光源のリストから、追跡したいn個のフォトン（光子）を作り出す。
- その各フォトンについて、光源からの経路を追跡して物体との衝突記録（フォトンキャッシュと呼ぶことにする）を作成する。
- フォトンキャッシュを標準出力へ書き出す。
大枠はこんなに単純だ（本当かどうか知らない、これは実験だから）。 さっそくコードにしてみよう。
main = do photons <- generatePhotons nphoton lgts photoncaches <- tracePhotons objs photons a <- forM (photoncaches) $ \i -> do putStrLn $ show i return ()
ほとんど上の仕様そのまんま（本当はここに至るまでに、特に最後のforMのところで
いろいろ試行錯誤したのだが、出来上がってみると単純だった）。
ここで、lgtsは光源のリスト、objsは物体のリストだ。あとは各関数を詳細化
していけばよい。Haskellだとこういうトップダウンでの開発がやりやすいように
感じる。もちろん、CでもJavaでもそういうアプローチをしているのだが、何をしたいか
の次に「どういう手続きにしたらいいか」を考えないといけないのが「素直じゃない」
と感じるところかもしれない。あくまで私見だが。Haskellだとこの次に
「generatePhotonsは何をするべきか」を記述していけばいいはず。
光源とフォトン
ちょっと横道に逸れて、シーンの情報をどうするか考えておく。上のlgtsやobjsの
実際の定義だ。いまのところ"登場人物"は光源と物体だけだ。レイトレーシングの段階に
なると視点や投影面など他にも必要だが、それは後で考えよう。ではまず光源だ。
以前の記事で、まずは点光源のみ扱うとした。それはそうなのだが、後々面光源なども
扱いたいので、型クラスLightを定義して、それを継承する形でPointLight型を作ろう。
フォトンマップを作るために、Lightには何が必要か。
一つは、その光源の光の強さを提供すること。フォトンが放出される量を知るためだ。
もう一つは放出されるフォトンの生成。フォトンの性質（色とか方向とか）は光源の状態に
依存するからだ。
前者の関数をflux、後者をgeneratePhotonとしよう。型クラスは次のように
なるだろう。
class Light a where flux :: a -> Flux generatePhoton :: a -> IO Photon
なにやらよく分からない型が登場している。いろいろ試した後だからサクッと
書くが、もちろん試行錯誤を繰り返した結果だ。Flux型は発光強度(放射束
(radiant flux)、単位は[W])を表すが、実際にはDoubleの別名、Photon型は
光源から放出される一つのフォトンを表す。フォトンをどう定義するか大変悩んだが、
最終的には次のような仕様にした。
特に「特定の波長の」が重要だ。一般にレイトレーシング法では光は白色光など 「色付き」で処理するはず。様々な波長が混ざった結果を「光」として扱っている。 しかし、「複数の波長が混ざった光が物体に衝突すると、果たして 反射光や屈折光はどう計算したらいいのか？」という疑問がわく。光が物体に衝突すると、 反射、屈折、吸収のいずれかが起こるが、どれぐらい反射するか、どれぐらい吸収するか などは反射率などで表される。問題はその率が「波長ごとに違う」ことだ。
例えば赤い物体が赤く見えるのは赤以外の光子が吸収されて反射しないから。 では、白色光が赤い物体に衝突したら、次は反射したとして追跡を続けるのか 吸収として停止するのか？ もちろん、各波長の反射率を平均して反射かどうか決めることはできるし、 白色光に各波長の反射率を掛け算して反射光としてもよい。しかし、各フォトンの強さが 一定のほうがレイトレーシング時の品質に有利と本に書いてあるので掛け算はNGだ。 一方で反射率等を平均したら「プリズムで分光された光は表現できるのか？」 「赤い球からの相互拡散反射の効果を表現できるのか？」などが疑問だった。
フォトンを特定の波長にするということは、全体的な「色」は複数のフォトンで 表現されるということ。白色なら赤、緑、青の少なくとも3つのフォトンがほぼ同数 必要だ。白色光ならフォトン1個ですむところ、3個必要なのだ。ならば、追跡する フォトンの数を増やさないと「色がまばらな」画像になってしまうかも知らない。 数が増えると大変だが、今回は実験なのだから計算量や時間にはこだわらず、素直な実装で行こう。 どうせ、現実世界は特定波長の光子の集まりで照らされているのだから。
あと、返り値がIO Photonなのは、generatePhoton関数の中で乱数を
使うから。乱数を使うということは「純粋」じゃなくなるのでIOをつけないと
ダメらしい。筆者はこれまでHaskellでどのように乱数を生成したらいいのか、
IO型を返す関数を実際のプログラムでどのように扱えばいいのかわかってなかった
ので今回のプログラムに手を出せなかったのだ。ここにきて、randomRIOの存在と
IOの使い所がほんの少し分かったので今回手を出してみた。
フォトンの定義
話が長くなった。結局Photon型は、波長＋放射の起点＋放射方向の情報があれば
よいだろう。波長はとりあえず赤、緑、青のいずれかとしてWavelength型を定義
しよう。
data Wavelength = Red | Green | Blue deriving (Show, Enum)
Enumクラスに属しておけば、色と番号を相互に変換でき、あとあと使えそうだ。
放射の起点と方向だが、これは要するに「直線」のベクトル表現である。
ある点を通り方向が
である
直線
は、
を任意の実数とすれば
と表される。これはレイトレーシング法で重要な「光線」そのものだ。
これをRay型として定義しておこう。Haskellではタプルを使えばよいだろう。
type Ray = (Position3, Direction3) initRay :: Position3 -> Direction3 -> Ray initRay p d = (p, d) target :: Double -> Ray -> Position3 target t (p, d) = p + t *> d getPos :: Ray -> Position3 getPos = fst getDir :: Ray -> Direction3 getDir = snd
ということでPhoton型はこれらのタプルで定義できる。合わせてフォトンの
衝突記録（フォトンキャッシュ）も定義しておこう。衝突記録には、波長、衝突場所、
光子の来た方向があればよい。光子の来た方向は、あとあと画像生成時に放射輝度を
計算するのに必要そうなので入れておく。結局、フォトンキャッシュもPhotonと
同じ構造だ。Rayの意味するところが違うから別の型PhotonCacheとしておこう。
type Photon = (Wavelength, Ray) initPhoton :: Wavelength -> Ray -> Photon initPhoton l r = (l, r) type PhotonCache = Photon
色
まだ、型の定義が続く・・・。光源の発する光の「色」を指定したい。
柔らかいオレンジがかったルームライトもあれば、緑のスポットライトもあるだろう。
これをColor型とする。各波長の比率として表せばよさそうだ。
data Color = Color Double Double Double initColor :: Double -> Double -> Double -> Color initColor r g b | mag == 0 = Color (1/3) (1/3) (1/3) | otherwise = Color (r'/mag) (g'/mag) (b'/mag) where r' = clipColor r g' = clipColor g b' = clipColor b mag = r' + g' + b' clipColor :: Double -> Double clipColor a | a < 0 = 0 | a > 1 = 1 | otherwise = a decideWavelength :: Color -> Double -> Wavelength decideWavelength (Color r g b) p | p < r = Red | p < r + g = Green | otherwise = Blue
比率なので、赤緑青の要素を足して1.0になるようにする。黒色は、色の
比率はなんでもよくて、光強度がゼロと考える。だから赤緑青が全部ゼロという
ColorはNGなのだが、万が一指定されたら「白」にする。
あと、decideWavelengthは0から1の間の実数を与えたら対応する波長を返す。
フォトンの波長をランダムに決めるときに使うものだ。だから、「比率」で表しておくのだ。
（エラー処理が面倒なので、0未満なら0、1以上なら1として処理する）
再び光源
長かった。下ごしらえができたので、PointLightの実装に移ろう。点光源なので、
持つべき情報は光色、発光強度、位置でいいだろう。
data PointLight = PointLight Color Flux Position3 instance Light PoingLight where flux (PointLight _ f _) = f generatePhoton (PointLight c _ p) = do theta <- randomRIO (0, pi) phi <- randomRIO (0, pi * 2) wl <- randomRIO (0, 1.0) let d = initDirFromAngle theta phi r = initRay p (fromJust d) w = decideWavelength c wl return (w, r)
ここで、initDirFromAngleが初出だが説明しておこう。方向ベクトルを
生成する関数の極座標版(?)だ。x, y, zを指定する代わりに2つの角度と
を与える。両方ともradianだ。
はY軸との角度、
はX軸からZ軸方向への角度とする。点光源は、その位置を中心にあらゆる
方向へ均一にフォトンが放出される。ランダムに放出方向を決めるには角度を乱数で
指定したほうが楽そうだったのでこのようにした。x, y, zの3つの乱数を使い、条件に
合わない場合は棄却する、という方法でランダムな方向ベクトルを得る方法もあるらしい。
ただ、よくわからないので無視。
まとめ
今回はフォトンマップ生成プログラムのメインルーチンを考えてみた。 たった数行だ。これを今後肉付けしていこう。そのための下ごしらえとして幾つかの型を 定義したので、次はこれらを使ってまずはフォトンをn個作る処理を考えようと思う。
レイトレーシング(2): `Algebra`モジュールをいじる
今回はテストについて書くつもりだったが、テストをいろいろ「テスト」していて
なかなか確認事が多そうなので、後回しにする。そこで、前回作った
Algebraモジュールをちょっといじろうと思う。
位置ベクトルと方向ベクトル
三次元ベクトルVector3を定義したが、実用的にはもう少し分類しておきたい。
具体的には、位置ベクトル(positional vector)と方向ベクトル(directional vector)に
分けて扱えるようにしたい。そこでVector3に別名をつけよう。
type Position3 = Vector3 type Direction3 = Vector3
これだけだと名前だけの話だ。今回のプログラムでは方向ベクトルは必ず
正規化されているもの、という制約をつけよう（ただ、邪魔くさそうなので
厳密にはやらないが）。そのため、生成時に必ず正規化されるようにinitDirを用意する。
形だけだが、位置ベクトルにもinitPosを用意しよう。
また、後々のため、ゼロベクトル、単位ベクトルも合わせて定義しておこう。
initDir :: Double -> Double -> Double -> Maybe Direction3 initDir x y z = normalize $ Vector3 x y z initPos :: Double -> Double -> Double -> Position3 initPos x y z = Vector3 x y z o3 = initPos 0 0 0 -- ゼロベクトル ex3 = fromJust $ initDir 1 0 0 -- 単位ベクトル(x) ey3 = fromJust $ initDir 0 1 0 -- 単位ベクトル(x) ez3 = fromJust $ initDir 0 0 1 -- 単位ベクトル(x)
initDirの結果にMaybeを使っているのは、引数が全部ゼロの場合、すなわち
ゼロベクトルがあり得るから。ゼロベクトルは方向ベクトルにできないので「値なし」
ということでNothingを返す。
関数名を演算子にしたい
ベクトルの加減算のため、maddやmsubなどを定義した。使い方は次のように
なるだろう。例として反射ベクトルの計算式
を示す。
r = msub d (mscale (2 * (dot n d)) n) -- 前置 もしくは r = d `msub` ((2 * (n `dot` d)) `mscale` n) -- 中置
なんと見難くて醜いことか。今回のネタではベクトル演算を多用するため、このままでは 見た目もタイプ量もデバッグにもよろしくない。なんとかできないものか。
幸いHaskellでは演算子も関数として定義できるらしい。代わりに+や-を定義してみよう。
class (Show a, Eq a) => BasicMatrix a where (+) :: a -> a -> a (-) :: a -> a -> a (中略) instance BasicMatrix Vector3 where (Vector3 ax ay az) + (Vector3 bx by bz) = Vector3 (ax + bx) (ay + by) (az + bz) (Vector3 ax ay az) - (Vector3 bx by bz) = Vector3 (ax - bx) (ay - by) (az - bz)
でもこれだと山ほどコンパイルエラーが出る。
$ ghc -o Algebra Algebra.hs
[1 of 1] Compiling Ray.Algebra ( Algebra.hs, Algebra.o )
Algebra.hs:50:57:
Ambiguous occurrence ‘+’
It could refer to either ‘Ray.Algebra.+’,
defined at Algebra.hs:13:3
or ‘Prelude.+’,
imported from ‘Prelude’ at Algebra.hs:5:8-18
(and originally defined in ‘GHC.Num’)
(以下同様)
加算の定義中で使われている'+'と、定義した'+'がぶつかっているようだ。修飾子をつけないと
ダメらしい。定義中のものにPrelude.を追加してみる。
instance BasicMatrix Vector3 where (Vector3 ax ay az) + (Vector3 bx by bz) = Vector3 (ax Prelude.+ bx) (ay Prelude.+ by) (az Prelude.+ bz) (Vector3 ax ay az) - (Vector3 bx by bz) = Vector3 (ax Prelude.- bx) (ay Prelude.- by) (az Prelude.- bz)
これでエラーが出なくなった。Vector3の定義はあまり綺麗ではないが、他のところで
スッキリ書けるならまあいいだろう。と思って、以下のようなサンプルを書いてみた。
module Main where import Ray.Algebra main :: IO () main = do let a = Vector3 1 2 3 let b = Vector3 3 4 5 let c = a + b putStrLn $ show c
そしたら、エラーが出た。
$ ghc -o ray Main.hs
[2 of 2] Compiling Main ( Main.hs, Main.o )
Main.hs:13:13:
Ambiguous occurrence ‘+’
It could refer to either ‘Prelude.+’,
imported from ‘Prelude’ at Main.hs:5:8-11
(and originally defined in ‘GHC.Num’)
or ‘Ray.Algebra.+’,
imported from ‘Ray.Algebra’ at Main.hs:7:1-18
だめだ、Algebraモジュール外でもエラーになる！引数がVector3型なのだから
どちらを使うかは自明と思うが? Haskellは型推論が優れていると自慢しているのに、
なぜこのぐらいの判断ができない?こんな修飾子を毎回書くのなら、せっかく式を簡略化
しようとしたのに本末転倒だ。ちょこちょこ調べたところnobsunさんの
コメント
を発見した。結局のところ、'+'とか'-'とかがNumクラスで定義されているせいだと。
しかしベクトル型をNumクラスのインスタンスにするのは無理がある。乗算とか。
だいたい、なぜ'+'を数値型前提で定義するのだろう。文字列でもなんでも数値以外にも
使いようがいっぱいあるのに。数学者が寄って仕様を作ったのかと思ってた。。。
と愚痴っても仕方ないので調べたところ、 NumericPrelude というのがあるそうな。これを組み入れてみよう。cabalでインストールする。
$ cabal install numeric-prelude
これを使うために、Algebraのソースを少々（いやかなり）いじる。
Additiveは加減算、Moduleはスカラー倍を定義しているクラス。
{-# LANGUAGE NoImplicitPrelude #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FlexibleContexts #-} : module Ray.Algebra where : import NumericPrelude import qualified Algebra.Additive as Additive import qualified Algebra.Module as Module : class (Show a, Eq a, Additive.C a, Module.C Double a) => BasicMatrix a where : instance Additive.C Vector3 where zero = Vector3 0 0 0 (Vector3 ax ay az) + (Vector3 bx by bz) = Vector3 (ax + bx) (ay + by) (az + bz) (Vector3 ax ay az) - (Vector3 bx by bz) = Vector3 (ax - bx) (ay - by) (az - bz)
これに倣い、他の関数名も変更してみた。
| 変更前 | 変更後 |
|---|---|
| madd | + |
| msub | - |
| mscale | *> |
| mdiv | /> |
| nearlyEqual | .=. |
| dot | <.> |
| cross | <*> |
mainの方も少し追加が必要だ。
{-# LANGUAGE NoImplicitPrelude #-} module Main where import NumericPrelude import Ray.Algebra main :: IO () main = do let a = Vector3 1 2 3 let b = Vector3 1 (-1) 1 putStrLn $ show (a + b) putStrLn $ show $ norm d putStrLn $ show (a <.> b) putStrLn $ show ((1.2::Double) *> a) putStrLn $ show (b /> 3) let w = fromJust $ normalize b let r = w - (2 * (w <.> ey3)) *> ey3 putStrLn $ show r
(詳しいことはソースを)
これだとVector3を使う方のソースは簡潔だし、正しく計算も出来ている！
$ ghc -o ray Main.hs [2 of 2] Compiling Main ( Main.hs, Main.o ) Linking ray ... $ ./ray [2.0,1.0,4.0] : [0.5773502691896258,0.5773502691896258,0.5773502691896258]
記述方法を比較してみる。
r = msub d (mscale (2 * (dot n d)) n) -- 前置 r = d `msub` ((2 * (n `dot` d)) `mscale` n) -- 中置 r = d - (2 * (d <.> n)) *> n -- 改善後
ああ、見やすくなった。（と自分では思う）
公開するもの、しないもの
中には他のモジュールには見せたくない関数なども含まれている（としよう）。
Javaでいうprivateメソッドとかのようなものだ。
このように外部に晒したくないものがある場合は、見せてよいものだけを
列挙したらいいらしい。モジュールの先頭で列挙するだけなので簡単だった。
特にAlgebraでは位置ベクトルと方向ベクトルを定義したので、Vector3で
直接生成したり要素を取り出したりできないようにしておく。また後で必要に
迫られたら考えたらいい。
module Ray.Algebra ( nearly0 , o3 , ex3 : , (+) , (-) : , , initPos , initDir ) where
これで、ここに書いてある以外の定義、関数などは他から使えなくなる。
関数や定数(o3とか)も並べるだけでよい。
今回はここまで。
