2015-05-24

レイトレーシング(6): やっとフォトン追跡

graphics

シーン作成

前回までで、物体を定義するための準備ができた。それを用いて簡単なサンプルシーンを作ってみる。閉じた空間でないと、せっかく放射したフォトンが無限遠に行ってしまって無駄になるので、箱型の空間を用意しよう。これで平面が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が具体的な追跡処理なのだが、実行することは次のようになるだろう。

フォトンをRayとして各物体と衝突する位置=交点までの距離を求める。
フォトンの放射された位置から見て前にある交点群について、最も近いものを選ぶ。
フォトンキャッシュ情報として返す。

これをコードにしたら次のようになった。

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に到達する。まあ、ファイル形式の再考はあとにしよう。

実行結果

ここまでくると、フォトンマップがどんな結果になったのかどうしても見てみたい。もちろん、ちゃんと動いているのか気になる所であるし。フォトンキャッシュ情報にある交点座標を、適当な二次元平面に投影させてフォトンの色も踏まえて描き出してみた結果が以下だ。

f:id:eijian:20150523204302p:plain

おお、なんとなくそれっぽいかも。ちゃんと部屋になって、右下に球も見える。黒いところは球の影でフォトンが届かないところだ。

物は試しに、複数の光源に変えて実行してみる。天井の中央ではなく、左右奥に白と橙色の 2つの光源だ。ただ、全体の光強度は同じ、フォトン数も同じだ。

f:id:eijian:20150523204303p:plain

なかなかよい。とりあえず、フォトンマップは完成した、ということにしよう。

メモリ使用量の改善

・・・いや、終わらなかった・・・。

フォトン数を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個のフォトンを生成する
フォトンを追跡してフォトンキャッシュを求める
フォトンキャッシュを出力する

最後にフォトンキャッシュを出力してしまえば、何も情報を保存しておく必要がないので、並列ではなく直列に1個ずつ処理を終わらせていけばいいということか。

まずi番目の光源から放出するフォトン数n(i)を計算し、
各光源iについて一個のフォトン追跡をn(i)回処理する

ような「二重ループ」にすればよいだろう。コードにしてみたのが以下。

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」まで激減した!アルゴリズムは大事だなぁ。（というか最初考えたアルゴリズムがあまりにもイケていなかっただけか・・・）

ようやく完成した。

2015-05-23

レイトレーシング(5): 物体の定義

graphics

前回までで、光源から放射されるフォトンが生成できたので、次はそれを追跡してフォトンマップを作ることになるが、そのためには描かれる「物体」を準備しないといけない。

「物体」定義

次に考える処理は、メインルーチン中の次の部分だ。

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

ここで、無限平面と球面は次の方程式を満たす三次元空間中の点 $\boldsymbol x$ (位置ベクトル)の集合である。

${ 無限平面: \boldsymbol n \cdot \boldsymbol x = d }$

${ 球面: ||\boldsymbol c - \boldsymbol x|| = r }$

ここで、 $\boldsymbol n$ は平面の法線ベクトル、 $d$ は平面の位置に関係するパラメータ、 $\boldsymbol c$ は球の中心座標、 $r$ は球の半径だ。このあたりの詳しいところはその筋の文献などを参照のこと。たとえば、

実例で学ぶゲーム3D数学

作者: Fletcher Dunn,Ian Parberry,松田晃一
出版社/メーカー: オライリージャパン
発売日: 2008/10/04
メディア: 大型本
購入: 21人クリック: 141回
この商品を含むブログ (40件) を見る

などに記載がある。上記のPlainとSphereの型コンストラクタの引数は、それぞれ $\boldsymbol n, d$ と $\boldsymbol c, r$ である。

「形」に必要な関数

フォトンの追跡は、まず物体と衝突する場所(交点)を求めることから始まる。バージョン1では反射は考えないので交点計算がやることのすべてと言っていい。 Shape型は上記の通り方程式で記述できるので、交点を求めるのは容易だ。光線(Ray) との連立方程式を解けばよい。交点を $\boldsymbol x = \boldsymbol p + t \cdot \boldsymbol d$ とすると、連立させて $t$ を求めれば、位置 $\boldsymbol x$ も解るわけだ。よって、 $t$ を計算する関数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

交点がない場合、複数の場合があるので、結果は $t$ のリストとする。

今後反射や屈折を考えたり、輝度計算をするときには交点での法線ベクトルを求める必要が出てくる。今はいらないが簡単なので定義しておこう。

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型にしている。

まとめ

今回はフォトン追跡の下準備として、「物体」を定義した。次回はフォトンを追跡することにしよう。

2015-05-11

レイトレーシング(4): フォトンの生成

graphics

光源の定義をやりなおす

前回の記事で、光源について型クラス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の値を生成してそれを各色の割合で直接選択しているのだからうまくいって当たり前か。

まとめ

とりあえず、フォトンを生成することはできた。次回は各フォトンを追跡してフォトンマップを作ることにしよう。

※ これまでのソースはここ。ただし、記事より先に進んでいる場合があるので注意。

2015-05-01

レイトレーシング(3): フォトンマップ生成の大枠を考える（ついでに光源も）

graphics

以前に書いたように、フォトンマッピング法は(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クラスに属しておけば、色と番号を相互に変換でき、あとあと使えそうだ。

放射の起点と方向だが、これは要するに「直線」のベクトル表現である。ある点 $\boldsymbol{p}$ を通り方向が $\boldsymbol{d}$ である直線 $\boldsymbol{r}$ は、 $t$ を任意の実数とすれば

${ \boldsymbol{r} = \boldsymbol{p} + t \cdot \boldsymbol{d} }$

と表される。これはレイトレーシング法で重要な「光線」そのものだ。これを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つの角度 $\theta$ と $\phi$ を与える。両方ともradianだ。 $\theta$ はY軸との角度、 $\phi$ はX軸からZ軸方向への角度とする。点光源は、その位置を中心にあらゆる方向へ均一にフォトンが放出される。ランダムに放出方向を決めるには角度を乱数で指定したほうが楽そうだったのでこのようにした。x, y, zの3つの乱数を使い、条件に合わない場合は棄却する、という方法でランダムな方向ベクトルを得る方法もあるらしい。ただ、よくわからないので無視。

まとめ

今回はフォトンマップ生成プログラムのメインルーチンを考えてみた。たった数行だ。これを今後肉付けしていこう。そのための下ごしらえとして幾つかの型を定義したので、次はこれらを使ってまずはフォトンをn個作る処理を考えようと思う。

2015-04-22

レイトレーシング(2): `Algebra`モジュールをいじる

graphics

今回はテストについて書くつもりだったが、テストをいろいろ「テスト」していてなかなか確認事が多そうなので、後回しにする。そこで、前回作った 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などを定義した。使い方は次のようになるだろう。例として反射ベクトルの計算式 $\boldsymbol{r} = \boldsymbol{d}-2(\boldsymbol{n} \cdot \boldsymbol{d}) \boldsymbol{n}$ を示す。

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とか)も並べるだけでよい。

今回はここまで。

2015-04-12

レイトレーシング(1): バージョン1の定義、ベクトル演算

graphics

Haskellは数学と関連があるというような話をちょくちょく見ることがある。圏論がどうとか数学的な概念が…といったところは筆者にはわからないが、ソース（見た目）はかなり数学っぽいと思う。実際これが一番Haskellにはまっている理由かもしれない。この簡潔さは素晴らしい。数学とくればレイトレーシング？ということで作ってみる。

手元にこんな本がある。

フォトンマッピング―実写に迫るコンピュータグラフィックス

作者: Henrik Wann Jensen,苗村健
出版社/メーカー: オーム社
発売日: 2002/07
メディア: 単行本
購入: 4人クリック: 58回
この商品を含むブログ (9件) を見る

古典的なレイトレーシングソフトは作ったことがあるので、今回はフォトンマッピングに手を出してみよう。そういうことだから、完成できる保証はない。また、途中でときどき別のネタに脱線すると思う。なお、ここでは理論の詳細には触れない。レイトレーシングのアルゴリズムや実装については上記の本やここを参照するとよいかもしれない。

"バージョン1"の定義

フォトンマッピング法に詳しいわけではないので、大ウソの連発かもしれないことはあらかじめ言い訳しておく。

さてこの手法は第一フェーズでフォトンマップを作成し、第二フェーズでレイトレーシングする、二段階アルゴリズムである。ただし第二フェーズでは普通のレイトレだと光源が見えるかどうか調べるが、この手法は光が到達する量をフォトンマップから「推定」する。これがみそ。

最初から超リアルな画像を生成できるものは無理なので、簡単なものを作って少しずつ肉付けしていけばいいだろう。バージョン1ではだいぶ単純化した仕様にして、まずは動くものを作る。以下が最初の仕様だ。

光源は点光源だけ
物体表面は拡散反射のみ（鏡面反射・屈折は無視）
フォトンの追跡は反射を無視（=相互拡散反射による効果はお預け）
物体は球と無限平面のみ
材質（というか色）は単色

この仕様で画像が生成できるのか今の時点ではよくわからないが、とりあえずは進めてみよう。

ベクトル演算

レイトレーシングの処理は、ほとんどが3次元ベクトルか光量(輝度）の演算で占められている。まずはベクトル演算のモジュールを作ろう。代数に関するモジュールなので名前をAlgebraにしよう。開発用ディレクトリの基本的な構造は以前に書いた通り (GitHub)。今回はsrcディレクトリの下にRayというディレクトリを作ってその中にソースファイルを作ることにする。トップディレクトリから見ると src/Ray/Algebra.hsだ。

そうそう、三次元座標系は筆者の好みで「左手系」かつy軸が上(x軸は右、 z軸は奥)を正方向とする。

主要なベクトル演算は型クラスで定義しておくと2次元ベクトルや行列など似たような型を定義するのにも使えそうである。どちらかというとベクトルは行列の特殊なものと考えれば、型クラスは

 BasicMatrix --> Matrix, Vector

という親子関係にしたほうがよさそう。BasicMatrixに行列やベクトルに共通な基本的な演算(関数)を定義し、特有の演算はそれぞれMatrixと Vectorクラスに定義するようにしよう。まずBasicMatrixで基本的な演算を定義する。加減算、スカラー倍、スカラー除算、ノルムにしよう。他に必要なものがあれば出てきてから追加する。なお、こっそりShowクラスと Eqクラスの子にしておく。 (実はこの歳になって初めてノルムにも色々な種類があると知った。ただここでは一般的(?)な、ベクトルで長さを意味するノルムとしよう。）

class (Show a, Eq a) => BasicMatrix a where
  madd :: a -> a -> a                 -- 加算
  msub :: a -> a -> a                 -- 減算
  mscale :: Double -> a -> a          -- スカラー倍
  mdiv   :: a -> Double -> Maybe a    -- スカラー除算
  mdiv a s                                                                      
    | s == 0    = Nothing                                                       
    | otherwise = Just ((1 / s) `mscale` a)
  norm :: a -> Double                 -- ノルム
  nearlyEqual :: a -> a -> Bool       -- ≒

class (BasicMatrix a) => Vector a where
  dot :: a -> a -> Double             -- 内積
  normalize :: a -> Maybe a           -- 正規化
  normalize a = a `mdiv` (norm a)
  square :: a -> Double               -- 二乗
  square a = a `dot` a

mdivは逆数を掛けるのと等しいことをクラス定義で記述しておく。除算なのでsが0の場合はエラーだ。ここでは解をMaybe型にし、エラーなら Nothingを返すようにしている。 Vectorクラスのnormalize(正規化)とsquare(二乗)についても同じくクラス定義で実装してしまおう。normalEqualはベクトル同士の比較用関数である。ご存知の通りコンピュータで実数を扱うと誤差が生じるので、理論上同一になる筈の結果がそうならないことがある。「誤差の範囲なら同じとみなす」ような比較用だ。使うかどうかわからないが。またこれに付随して(?)、nearly0::Doubleも定義しておく。

のちのちMaybe型に関連する関数を使うためにはData.Maybeモジュールをimportしないといけない。今のうちに入れておく。

import Data.Maybe

なお、Matrix(行列)は将来的には使うが、とりあえず今は無視する。

次に三次元ベクトル型を定義しよう。Vector3だ。

data Vector3 = Vector3 Double Double Double

(中略)

instance Matrix Vector3 where
  madd (Vector3 ax ay az) (Vector3 bx by bz) = Vector3 (ax + bx) (ay + by) (az + bz)
  msub (Vector3 ax ay az) (Vector3 bx by bz) = Vector3 (ax - bx) (ay - by) (az - bz)
  
(中略)

cross :: Vector3 -> Vector3 -> Vector3                                          
cross (Vector3 ax ay az) (Vector3 bx by bz) = Vector3 (ay * bz - by * az) (az * bx - bz * ax) (ax * by - ay * bx)       


(以下続く…)

外積だけは三次元ベクトル特有の演算なので（本当かどうか知らない）、クラス定義には含められず独立した関数crossとして定義した。動作確認のため、対話環境(ghci)で試す。以下はsrcディレクトリ内で実行した場合である。

$ ghci
Prelude> :l Ray.Algebra
[1 of 1] Compiling Ray.Algebra      ( Ray/Algebra.hs, interpreted )
Ok, modules loaded: Ray.Algebra.
*Ray.Algebra> let a = Vector3 1 2 3
*Ray.Algebra> a
[1.0,2.0,3.0]
*Ray.Algebra> let b = Vector3 4 5 6
*Ray.Algebra> putStrLn $ show $ madd a b
[5.0,7.0,9.0]
*Ray.Algebra> putStrLn $ show $ msub a b
[-3.0,-3.0,-3.0]
*Ray.Algebra> 
*Ray.Algebra> putStrLn $ show $ mscale 5 a
[5.0,10.0,15.0]
*Ray.Algebra> putStrLn $ show $ mdiv a 5
Just [0.2,0.4,0.6000000000000001]
*Ray.Algebra> putStrLn $ show $ norm a 
3.7416573867739413
*Ray.Algebra> putStrLn $ show $ dot a b
32.0
*Ray.Algebra> let c = normalize a
*Ray.Algebra> c
Just [0.2672612419124244,0.5345224838248488,0.8017837257372732]
*Ray.Algebra> let d = fromJust c
*Ray.Algebra> d
[0.2672612419124244,0.5345224838248488,0.8017837257372732]
*Ray.Algebra> putStrLn $ show $ norm d
1.0
*Ray.Algebra> let x = Vector3 1 0 0
*Ray.Algebra> let y = Vector3 0 1 0
*Ray.Algebra> let z = cross x y
*Ray.Algebra> z
[0.0,0.0,1.0]

それなりにうまくいっているようだ。が、いくつか適当な値で試しても、正直な所ちゃんとテストできているかどうかわからない。ということで次回はユニットテストを考える。

2015-03-22

同一画像検索(6): 改良して完成

graphics

このネタの最後に幾つか確認と改良をして完成させよう。

総当たり処理を少し改善

roundRobinは再帰で定義しているが、場合分けが格好悪い。isSameの戻り値がBoolなところが問題か。

roundRobin x (y:ys)
  | isSame x y == False = roundRobin x ys
  | otherwise           = (snd y):(roundRobin x ys)

同一ならyのFilePathを、そうでなければ空リストを返せばよさそう。ついでに、前回isSame内で同一かどうかを判定するのにfindを使ったのをanyに変えておこう。「同一」かどうかBoolで返してくれればよいのでfindである必要はない。（単にanyを知らなかっただけ）

roundRobin x (y:ys) = isSame x y ++ roundRobin x ys                       

isSame :: Image -> Image -> [FilePath]                                 
isSame x y = if any differ (zip (fst x) (fst y)) then [] else [snd y]         
  where                                                                         
    differ :: (Word8, Word8) -> Bool                                            
    differ (a, b) = (if a > b then a - b else b - a) > threshold

これでだいぶすっきりした。

重複した出力を取り除く

前回の出力結果を再掲する。

probably same: work/IMG_0309-2.jpg, work/IMG_0309-3.jpg, work/IMG_0309-4.jpg, work/IMG_0309.jpg                                                                
probably same: work/IMG_0309-3.jpg, work/IMG_0309-4.jpg, work/IMG_0309.jpg      
probably same: work/IMG_0309-4.jpg, work/IMG_0309.jpg                           
probably same: work/sample1.jpg, work/sample7.jpg                               
probably same: work/sample2.jpg, work/sample5.jpg

2行目、3行目は1行目の部分集合であることがわかる。部分集合かどうかを調べるのはHaskellなら簡単にできそう。 Webを検索したらやはり、Data.List の中にそのものズバリisInfixOfがあった！matchImageから返ってくるリストについて、ある要素が他のすべての要素と比較してどれの部分集合でもなければ自身を返すようにすればよいだろう。 findSameの中でmatchImageの結果を渡すようにする。

findSame fs = do                                                            
  fps <- mapM (getFingerPrint 4) fs                                             
  let ps = matchImage $ zip fps fs                             
  return $ deduplicate ps ps                                                    

deduplicate :: [[FilePath]] -> [[FilePath]] -> [[FilePath]]                     
deduplicate _ [] = []                                                           
deduplicate xs (y:ys)                                                           
  | any (y `isInfixOf`) xs = deduplicate xs ys                                   
  | otherwise              = y:deduplicate xs ys

動かしてみると。。。だめだ、一件も同じ画像とみなされなくなった!? ここで30分ほどハマった。確認したいのは自分が「他の要素に」含まれているかどうかだ。しかし上記のxsには「自分自身も」含まれている! これでは自分自身にマッチしてしまって全要素が消える。自分以外でという条件をつけよう。

deduplicate xs (y:ys)                                                           
  | any isProperSubset xs = deduplicate xs ys                                   
  | otherwise             = y:deduplicate xs ys                                 
  where                                                                         
    isProperSubset :: [FilePath] -> Bool                                        
    isProperSubset x = x /= y && y `isInfixOf` x

簡潔に書く方法がよく分からないのでちょっと面倒くさくなってしまったが、とりあえず無駄な要素は除去できた!

probably same: work/IMG_0309-2.jpg, work/IMG_0309-3.jpg, work/IMG_0309-4.jpg, work/IMG_0309.jpg                                                                
probably same: work/sample1.jpg, work/sample7.jpg                               
probably same: work/sample2.jpg, work/sample5.jpg

パラメータを引数で与える

ここまでのところ、ソース中に2つの定数が埋め込まれている。fingerprintの解像度と画像比較時の差の閾値だ。何度か試して適した値を埋め込んでおくのも良いが、試すためにもコマンド実行時にいろいろ変えて与えたい。そこで、"-p"オプションをサポートしよう。ただし、細かいエラーチェックは面倒なので割愛する。第一引数が"-p"で始まっていたらオプションが指定されたとし、そうでなければ第一引数も処理対象のディレクトリとみなす。オプションは"-pR,T"で、Rが解像度、Tが閾値。両パラメータとも正整数である前提だ。よって、変な引数を与えたときの動作は保証されない。

まず、オプションをちゃんと理解できたと仮定して、その後の処理ができるように改造しよう。両パラメータともfindSameに渡す必要があるので関数定義を変更する。

findSame :: Int -> Int -> [FilePath] -> IO [[FilePath]]                         
findSame r t fs = do                                                            
  fps <- mapM (getFingerPrint r) fs                                             
  let ps = matchImage (fromIntegral t) $ zip fps fs

第一引数が解像度、第二が閾値だ。解像度をgetFingerPrintの引数にそのまま渡せば良い。閾値は初お目見えなのでmatchImageに渡して最終的にはisSameで条件判定に使われるようにしておく。 (ソースはこちら)

なおfromIntegral tとしているのはisSame内ではWord8として比較しているからIntのまま渡せないため。

下準備ができたところで引数処理に移ろう。引数全部を渡してオプションの有無、両パラメータを処理するparseOptを定義する。戻り値は解像度、閾値、処理対象ディレクトリのリスト、の3つ。

parseOpt :: [String] -> (Int, Int, [FilePath])                                  
parseOpt (d:ds)                                                                 
  | "-p" `isPrefixOf` d = (r, t, ds)                                            
  | otherwise           = (8, 8, d:ds)                                          
  where                                                                         
    [r, t] = map (read :: String -> Int) (splitOn "," (drop 2 d))

当初、どうしたら引数の有無やパラメータを取り出せるかだいぶ悩んだが、まずは第一引数が"-p"で始まっていなければ、与えられたリスト(=d:dsだ)をそのまま、パラメータはデフォルト値(両方とも8とした)を返せば良いとした(otherwiseの行)。"-p"で始まるかどうかは、先に "部分集合"の判定を考えていた時にisPrefixOfもチェックしていたのでそれが使えると判断。あとはそう難しくない。最初の二文字(="-p")を除き、","(カンマ)で分割、それぞれの文字列を「正整数と仮定」して Intに変換すればよい。最後の変換のところ、ちょっと立ち止まったが、最終的には上記の通りread関数でなんとかなった。（エラー処理を無視すれば）

$ ghc -o picf Main.hs
$ ./picf -p16,4 ~/work
    :

ちゃんと動く！パラメータを変えると条件がきつくなって同一と表示されなくなる。スバラシイ。

cabalでコンパイル

この件の最初の回でcabalを使う準備をしていながら最後まで何も使わないのはもったいないので、本プログラムをcabalでコンパイルしてみる。 cabalを使った一連の流れはここを参考にした。cabal buildとすればよいらしい。cabalファイルはプロジェクトディレクトリのトップにあるのでそこで実行する・・・と、エラーが出た。

$ cabal build
./picfinder.cabal has been changed. Re-configuring with most recently used
options. If this fails, please run configure manually.
Warning: The package list for 'hackage.haskell.org' is 103 days old.
Run 'cabal update' to get the latest list of available packages.
Resolving dependencies...
Configuring picfinder-0.1.0.0...
Building picfinder-0.1.0.0...
Preprocessing executable 'picfinder' for picfinder-0.1.0.0...
cabal: can't find source for Main in .

ソースが見つからないだと。指定していないから当たり前だ。 hs-source-dirsで指定するらしい。気を取り直して。

$ cabal build
./picfinder.cabal has been changed. Re-configuring with most recently used

(中略)

Preprocessing executable 'picfinder' for picfinder-0.1.0.0...

src/Finder.hs:7:8:
    Could not find module ‘Data.ByteString’
    It is a member of the hidden package ‘bytestring-0.10.4.0’.
    Perhaps you need to add ‘bytestring’ to the build-depends in your .cabal file.
    Use -v to see a list of the files searched for.

(以下略)

よくわからないがbuild-dependsにXXXXを足せとある。今度は成功した！

$ cabal build
./picfinder.cabal has been changed. Re-configuring with most recently used

(中略)

Resolving dependencies...
Configuring picfinder-0.1.0.0...
Building picfinder-0.1.0.0...
Preprocessing executable 'picfinder' for picfinder-0.1.0.0...

実行ファイルはdist/build/picfinder/picfinderとしてできているらしい。実行してみる。

$ dist/build/picfinder/picfinder -p4,4 ~/work
probably same: work/IMG_0309-2.jpg, work/IMG_0309-3.jpg, work/IMG_0309-4.jpg, work/IMG_0309.jpg
  :

やっとここまでたどり着いた。最終回ということで詰め込みすぎた感はあるがよしとしよう。

次は何をしようか。

Haskellでいってみよう

日曜プログラマにも満たないレベルでもHaskellで何かソフトウェアを作りたい！

レイトレーシング(6): やっとフォトン追跡

シーン作成

フォトンの追跡

フォトンキャッシュ情報の書き出し

実行結果

メモリ使用量の改善

レイトレーシング(5): 物体の定義

「物体」定義

「形」に必要な関数

まとめ

レイトレーシング(4): フォトンの生成

光源の定義をやりなおす

フォトンの生成

生成されたフォトンの品質

まとめ

レイトレーシング(3): フォトンマップ生成の大枠を考える（ついでに光源も）

仕様とメインルーチン

光源とフォトン

フォトンの定義

色

再び光源

まとめ

レイトレーシング(2): `Algebra`モジュールをいじる

位置ベクトルと方向ベクトル

関数名を演算子にしたい

公開するもの、しないもの

レイトレーシング(1): バージョン1の定義、ベクトル演算

"バージョン1"の定義

ベクトル演算

同一画像検索(6): 改良して完成

総当たり処理を少し改善

重複した出力を取り除く

パラメータを引数で与える

cabalでコンパイル