レイトレーシング(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
とか)も並べるだけでよい。
今回はここまで。
レイトレーシング(1): バージョン1の定義、ベクトル演算
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]
それなりにうまくいっているようだ。が、いくつか適当な値で試しても、 正直な所ちゃんとテストできているかどうかわからない。ということで 次回はユニットテストを考える。
同一画像検索(6): 改良して完成
このネタの最後に幾つか確認と改良をして完成させよう。
総当たり処理を少し改善
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 :
やっとここまでたどり着いた。最終回ということで詰め込みすぎた感は あるがよしとしよう。
次は何をしようか。
同一画像検索(5): 同一判定方法の再考
前回簡易版を作ってみたが、縮小した画像が「同一」と見做されなかった。 今回はこれにどう対処するか検討してみる。 この問題は、各画像を4x4解像度に変換した結果(=fingerprint)を取った時に、 元画像と縮小画像とで微妙に異なるのが原因だった。 簡易版は「同一画像ならfingerprintが全く同じになる」ことを前提に作って あったからだ。
そこで、同一画像と見做せるようにする対策案を考えてみた。
- 解像度をもっと減らす
- 色階調を減らす
- fingerprintの微妙な違いを許容する
a. 解像度をもっと減らす
前回のサンプル画像で試してみる。解像度を2x2に減らしたところ、2つの画像の 差はなくなった。これは行けるかもと思ったが、別の画像で試したらやはり差が生じた。 また、そもそも2x2程度では、異なる画像がたまたま同一のfingerprintになる 可能性も高まる。この方法は却下だ。
b. 色階調を減らす
各色256階調、3原色で1600万色であるから、画像の縮尺が変われば同じ画像でも 多少の差が出てきてしかるべき。であれば、階調を減らせば多少の差は吸収される のではないか。階調を減らすには256段階(=1byte)の下数桁をマスクして 0にしてやれば良いだろう。たとえば、
<色値> AND 0xFC (0xFCは下2桁を除外するマスク→64段階にする)
といった処理をfingerprintの全バイトへ適用するということ。
しかしこの方法は根本解決にならない。階調を減らしたところで近い色値が「同一」 になる保証はない。極端な例だと、2つの画像のある点で色値がそれぞれ127と128と なるような場合、きわめて近い色なのにどのようなマスクを設定しても「同一」と 見做されない。ゆえにこの方法も却下だ。
c. fingerprintの微妙な違いを許容する
微妙な違いを許容して「同一」と考えるために、2画像のfingerprintの各点が閾値 以下の差なら「同一」とみなそう。(当初のロジックで、2段階目に解像度16x16で やろうとしたことだ)
元々、fingerprintが一致する前提ならMap
を使って計算量を抑えられそうと
見込んだ。Map
なら計算量は (?)。しかし、差を調べるには総当たり
するしかないので計算量は だ。画像数が数千枚レベルになると
大変そうなので嫌なのだが・・・仕方ない。
総当たりなので基本的なロジックは、、、
- リスト中の最初の画像とそれ以降の画像のfingerprintを比較して同じと みなせるものがあればそれらをリストで返す => (n-1回の比較)
- 二番目の画像とそれ以降の・・・(以下略)=> (n-2回の比較)
- 三番目の・・・(以下略)
- ・・・n-1番目の・・・(以下略)=> (1回の比較)
- => しめて 回の比較
だろう。返ってくるのは「同一」と判断されたリストの集合(リスト)。
ちなみに、比較がやり易いようにFingerPrintの型をByteString
から
[Word8]
に変更し、画像を表すあらたな型Image
も定義してみた。
type FingerPrint = [Word8] type Image = (FingerPrint, FilePath)
先のロジックをコードにすると次のような感じ。
matchImage :: [Image] -> [[FilePath]] matchImage [] = [] matchImage (x:[]) = [] matchImage (x:xs) | ps == [] = matchImage xs | otherwise = (snd x:ps):(matchImage xs) where ps = roundRobin x xs roundRobin :: Image -> [Image] -> [FilePath] roundRobin x [] = [] roundRobin x (y:ys) | isSame x y == False = roundRobin x ys | otherwise = (snd y):(roundRobin x ys)
matchImage
でn-1回の繰り返し、roundRobin
である画像とその他の
画像を比較して同じとみなせるものがあればリストにして返す、とやっている。
なおisSame
が同一かどうかの判別関数であり、同一の場合はその画像を
リストに加えている。
画像(のfingerprint)を同じとみなす方法もいくつかある(fingerprintをn次ベクトル
とみて距離を取るとか、同色の点が何割以上とか)が、今回は「各点の色の差が
すべて閾値以下」を条件とした。上記の通りfingerprintをWord8
(1バイト符号なし整数)のリストと定義しなおしたので、二つのリストを先頭から順に
比べていけばよいだろう。
閾値を超えたところで「違う」のだから比較をやめればよい。
isSame :: Image -> Image -> Bool isSame x y = d == Nothing where d = find differ (zip (fst x) (fst y)) differ :: (Word8, Word8) -> Bool differ (a, b) = d' > threshold where d' = if a > b then a - b else b - a threshold = 4 :: Word8
何かもっとスマートな書き方がありそうだが…。これは今後の課題にしよう。
threshold
は各点の各色の差の許容範囲(閾値)。ここでは4にしたが、
適当に決めた。色階調が256段階だから1.5%ぐらい。どれぐらいが適当か不明だが
5%ぐらいまでは許容範囲かなと思う。
あと、getFingerPrint
にバグがあった。修正したものは以下の通り。
getFingerPrint :: Int -> FilePath -> IO FingerPrint getFingerPrint r f = do (sin, sout, serr, ph) <- runInteractiveCommand command waitForProcess ph fp <- BS.hGet sout size return $ B.unpack fp where geo = (show r) ++ "x" ++ (show r) size = r * r * 3 command = "convert -filter Cubic -resize " ++ geo ++ "! " ++ f ++ " PPM:- | tail -c " ++ (show size)
まず、コマンドの出力を取り込むところ、BS.hGetLine
をBS.hGet
に訂正。
出力にたまたま改行コードと同じものがあるとそこで切れてしまうバグだった。
次にImageMagickのコマンド部分、引数から-define
を抜いた。これをつけると
縮小処理が大幅に高速化するのだが、逆に画像の品質が悪くなった。今回の場合
fingerprintの誤差が本来よりかなり大きくなってしまった。これはいただけない。
(気付くのにだいぶ時間がかかった…)
実行してみた結果がこれ。
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
前回試した画像に対し、50%、25%縮小画像でもちゃんと同一とみなしてくれている。 他にも2種類の画像(sample?.jpg)を試したが、いずれも認識できた。やれやれ。 ただ、一行目と二行目を比べるとわかるように、3つ以上の画像A,B,Cがあって、 それらが互いに「同一」と判断された場合は複数行に結果が出てくる。しかも 後の方は前のリストに包含されているので無駄だ。これは何とかしたい。
ということで、次回はこの件の最終回。上記の無駄を省く処理と若干の ソースの改善をしたいと思う。
同一画像検索(4): 簡易版を作ってみた
これまでの確認などを踏まえ、今回は簡易版を作ろう。仕様では4x4解像度の画像比較をして 一致したら次に16x16解像度でもう少し詳細に比較する、としているところ、まずは4x4解像度での 比較だけで作ってみようと思う。
最初に、「外枠」で示したfindSame
関数の型シグネチャに問題があったので訂正する。
これまで示したのは下記の形だった。
findSame :: [String] -> [[String]]
画像比較のためImageMagickを呼び出す仕様なので、どうしても副作用が発生する。
haskellではそれが含まれる場合はどうやっても純粋な(副作用のない)型シグネチャには
できなさそうだ。(言われてみれば、その関数の処理のどこか一部に副作用があるなら
全体として関数の処理結果が引数だけで一意に決まらないのは当然か。)
下記が訂正版とその呼び出し元main
の変更部分。
findSame :: [String] -> IO [[String]] : main = do ds <- getArgs fs <- mapM getFileLists ds ss <- findSame $ concat fs putGroups ss
さて、肝心のfindSame
をどうするかだが、処理の流れは次の通りかと。
- 各画像のfingerprint(前々回説明)を取得(ここでImageMagickを使う)
- 同じfingerprintを持つ画像群をListにまとめる
- 同じものが見つかった(=Listの長さが2以上)ら、そういった画像群のListを返す
これを関数にしてみる。プログラムが上の"処理の流れ"そのまんまなのはさすが haskellといったところか。なお、以後のコード片ではファイル名をStringではなく 別名のFilePath(Preludeで定義済み)、fingerprintもByteStringではなく FingerPrint型とした。
findSame :: [FilePath] -> IO [[FilePath]] findSame fs = do fps <- mapM getFingerPrint4 fs let es = Map.elems $ foldl insertItem Map.empty (zip fps fs) return $ filter (\x -> length x > 1) es
getFingerPrint4
は4x4解像度のfingerprintを取得する関数だが、実際は
解像度を引数で与えるgetFingerPrint
に第一引数を4としたもの。この辺も
関数型言語らしい。
getFingerPrint :: Int -> FilePath -> IO FingerPrint getFingerPrint r f = do (sin, sout, serr, ph) <- runInteractiveCommand command waitForProcess ph BS.hGetLine sout where geo = (show r) ++ "x" ++ (show r) size = r * r * 3 command = "convert -define jpeg:size=" ++ geo ++ " -filter Cubic -resize " ++ geo ++ "! " ++ f ++ " PPM:- | tail -c " ++ (show size) getFingerPrint4 = getFingerPrint 4
insertItem
は前回書いたままである。万が一、全ソースに興味がある方は
GitHubをどうぞ。
早速実行して試してみよう。以下がサンプルで使った画像(の25%縮小版)。
オリジナル(IMG_0309.jpg)、その単純なコピー(IMG_0309-2.jpg)、
50%縮小版(IMG_0309-3.jpg)、25%縮小版(IMG_0309-4.jpg)を用意して試した。
結果は惨敗。
$ ghc -o picf Main.hs [1 of 2] Compiling Finder ( Finder.hs, Finder.o ) [2 of 2] Compiling Main ( Main.hs, Main.o ) Linking picf ... $ ./picf ~/work probably same: work/IMG_0309.jpg, work/IMG_0309-2.jpg
同じサイズの画像しか同一とみなしていない。オリジナルと50%版のfingerprintを 比べてみる。
$ convert -define jpeg:size=4x4 -filter Cubic -resize 4x4! ~/work/IMG_0309.jpg PPM:- | tail -c 48 | od -x 0000000 4f47 5d52 6561 9996 bd9b c6c1 3d3e 4d36 0000020 3b43 6171 7e59 7476 555d 5a40 3b4c 596f 0000040 7f47 5970 707a 7b4f 4f6e 7485 8c54 5a7c 0000060 $ convert -define jpeg:size=4x4 -filter Cubic -resize 4x4! ~/work/IMG_0309-3.jpg PPM:- | tail -c 48 | od -x 0000000 4f47 5d52 6561 9996 bd9b c5c1 3d3e 4d37 0000020 3b43 6171 7e59 7476 555d 5a40 3b4c 596f 0000040 7f47 5970 707a 7b4e 4f6e 7485 8c54 5a7c 0000060
違う箇所が3つあるので3つの画素でRGBのうち一色だけ1/256の差があると。。。 本当に微妙な差だが、fingerprintが違うのだから「別の画像」ですね・・・。
さてどうしたものか。ちょっと考えてみて、結果は次回で報告。