【論文】p = k ( ρ - ρ0 )でのρ0の意味

Mullerらの"Particle-Based Fluid Simulation for Interactive Applications"では、圧力pを密度ρから求める際に、定常密度ρ0を引いて

p = k(ρ - ρ0)

としている（kは定数）。ρ0を引く理由としてMullerらは、

... Since pressure forces depend on the gradient of the pressure field, the offset mathematically has not effect on pressure forces. However, the offset does influence the gradient of a field smoothed by SPH and makes the simulation numerically more stable. ...

と述べている。

なんのこっちゃ？

オフセット（ρ0）が数学的に圧力に影響を与えないといいつつ、それによってシミュレーションがより安定化するといっている。しかし、数学的に影響を与えないのであれば、安定化するもなにもないではないか。

そこで、引用元である、Desbrunの"Smoothed Particels: A new paradigm for animating highly deformable bodies"の中で、なぜρ0が導入されているのかを見てみた。

Desbrunらによると、

... In astrophysics applications, pressure forces were often combined with gravitational forces balancing the expansion phenomenon.
  In contrast, we would like to animate materials with constant density at rest. Consequently, the material sould exhibit some internal cohesion, resulting in attraction-repulsino forces as in the Lennard-Jones model. ...

とある。

なるほど。

もともとSPHは天体のシミュレーションで用いられていた手法なので、真空中での粒子の運動を対象としている。一方で、Desbrunらが行おうとしている大変形体のシミュレーションは、大気中での運動が対象となる。

大気中での運動なので、本来であれば、大変形体の周りに空気があって、その気圧と大変形体の圧力が釣り合っていないといけない。

それを論文でどう扱っているかというと、周りの空気も粒子としてシミュレーションに含めるのではなく、その代わりに、ρからρ0を引いて気圧の影響を補正した圧力を用いるという方法をとっている。つまり、定常状態での圧力を仮想的に0気圧となるようにオフセットをとって計算している。

つまり、定常状態での圧力を仮想的に0気圧とすることで、気圧と大変形体の圧力を釣り合わせているということになる。

実際に、ρ0を変えてシミュレーションをしてみた。

【ρ0=600】


【ρ0=300】


【ρ0=0】


確かに、ρ0が小さくなるにつれて、粒子が空気中に広がっていっている様子がわかる。ρ0=0に至っては、まさに真空中に拡散しているように見える結果が得られた。

結論
圧力pを密度ρから求める際に定常密度ρ0を引いてp = k (ρ - ρ0)とするのは、気圧の影響を加味して粒子が拡散してしまわないようにするためだった。SPHを宇宙から地球上に持ってくるわけですな。

■関連する記事
粒子法のプログラム第１回（概要）
粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）
【粒子法】粒子を流体としてレンダリング
３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…

ClickTaleの無料利用分を使い切ってしまった

おととい入れてみたClickTale。

２日で1週間分のレコーディング数を使い切ってしまった…。これじゃ使えないなぁ。

１人１人のユーザの動きを追える素晴らしいツールなんだけど、月$99は払えないや。

自分で作ろうか？？

ClickTaleの再生がAnalyticsに計測されてしまわないか？

このブログにClickTaleを入れてみた。

ClickTaleはユーザのマウスの動きをレコーディングしてくれるアクセス解析ツールで、レコーディングされたマウスの動きを再生できる。

再生機能で１つ気になったのが、再生したときのページ読み込みが解析対象ページに埋まっているGoogle Analyticsで計測されてしまわないかという点。

Analyticsで確認してみると、ClickTaleからのアクセスは計測されないみたいでよかった。

OneApp

MicrosoftがOneAppという携帯電話向けのアプリケーションフレームワークを発表した。

OneAppは、途上国の人々がTwitterやFacebookのようなサービスに携帯電話からアクセスできるようにするサービスで、プロセッサが非力でメモリが少なくてもJavaが動きさえすれば使うことができる。

大型Webアプリケーションを動かすために必要な処理の一部をMicrosoftのサーバがクラウドとして代行する仕組みで、デベロッパはJavaScriptとXMLによる比較的簡単な作業でアプリをOneAppに対応させることができる。

Microsoftが途上国向けにサービスを提供する理由には、人道的理由もあるが、途上国が潜在的に巨大な市場であるというビジネス上のものもあるだろう。

iTunesでムービーのサムネイルを変えるには

iTunesのムービーのサムネイルは好きなタイムラインのものに変えることができる。

▼変え方

1. サムネイルにしたいタイミングでムービーを停止する
2. ムービーを右クリック
3. 「ポスターフレームを設定」を選ぶ

これで好きなカットをサムネイルにできる♪

粒子法のシーンを２倍のサイズにしてみた

「粒子法のシーンを２倍のサイズにしてみたが…」で失敗したものをもう一度計算した。

【２倍のサイズ】


【元のサイズ】


初期配置や境界の形はまったく同じにして、xyz各方向に２倍している。一見同じように見えるけど、床のチェック模様の細かさが違うところに注目して欲しい。

２倍サイズにすると、より高いところからさらに多くの流体が落下するため、盛大なしぶきがあがった。上のレンダリング結果だと黒くなってしまってわかりづらいが、陰関数曲面をそのままレンダリングするとその様子がよくわかる。

粒子法(SPH)のプログラム

粒子法(SPH)のプログラムを解説したシリーズです。ソースコードも公開しています。

粒子法(SPH)のプログラム一覧

動画

シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。

動画の一覧

--

【論文】Particle-Based Fluid Simulation for Interactive Applications

Matthias Mullerらの「Particle-Based Fluid Simulation for Interactive Applications」を読んだ。

この論文では、SPHをベースに、自由表面を持つ流体をインタラクティブにシミュレーションしレンダリングする方法を提案している。

この論文におけるMullerらの貢献は、

• Desbrunらの大変形体向けの手法を流体に拡張
• インタラクティビティの実現

の２点。

Desbrunの手法を流体に拡張するにあたり、Mullerらは流体が受ける力をナビエ=ストークス方程式から直接導出するようにし、さらに、表面張力をモデル化した項を加えることを提案している。

また、インタラクティブに流体を描画できるように、高速に計算できるカーネルを設計すると共に、キャッシュヒット率が向上するようデータ構造を工夫している。

レンダリングにはPoint SplattingとMarching Cubesの２通りの手法を用いており、2003年当時のPC(Pentium4+GeForce4)で2000個程度の粒子を計算した例では、Point Splattingで20fps、Marching Cubesで5fpsでの描画を実現している。

ハンニバル戦記

塩野七生の「ローマ人の物語」を読んでいる。

今日は、２巻の「ハンニバル戦記」を読み始めた。ハンニバル戦記は、そのタイトルの通りハンニバルのアルプス越えで有名なポエニ戦役について書かれている。

ポエニ戦役は紀元前264年〜同133年までローマとカルタゴの間で行われた戦争で、今日読んだのはそのうちポエニ戦役のきっかけとなったメッシーナでの戦い。

紀元前265年、シチリア半島の東端に位置するメッシーナは、南に位置するシラクサから攻め込まれていた。

メッシーナは、海峡を挟んで対岸にあるローマへ救援を求める。頼られたローマは迷ったあげく、クラウディウス率いる軍を派遣した。

ローマはシラクサ、カルタゴを撃破。シラクサとは同盟関係を結ぶことでシチリア島に影響力を伸ばすことに成功する。

しかし、これによってカルタゴの勢力圏と直接接することになり、その後のポエニ戦役へと繋がっていく…

次はいよいよポエニ戦役の始まり！

由比ヶ浜

夕方に鎌倉の由比ケ浜までふらっと行ってきた。

鎌倉には何度かいったことがあるけど、由比ヶ浜にいったのは初めて。

行ってみると、広いし、水もきれいだし、設備もあるし、駅から近い。すごく良い！

海というと子供のころの印象が強くて、車で何時間もかけて頑張っていかないといけないというイメージがあったけど、由比ケ浜は全然そんなことがない。

すごくいいところを見つけた♪

粒子法のシーンを２倍のサイズにしてみたが…

３次元の粒子法シミュレーションを、さらにシーンを２倍のサイズにして計算してみた。



あ…カメラの位置を調節していなかった…
再計算します。

追記：再計算しました。こちら


■関連する記事
粒子法のプログラム第１回（概要）
粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）
【粒子法】粒子を流体としてレンダリング
３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…
粒子法のシーンを２倍のサイズにしてみた

３次元の粒子法シミュレーション

粒子法のプログラム第１回（概要）や【粒子法】粒子を流体としてレンダリングでは粒子を２次元で計算していたので、今回はそれを３次元に拡張して計算してみた。平面だったものを、そのままz方向にも粒子を配置している。


【粒子表現】


【流体表現】



■流体シミュレーションに関するエントリ
粒子法のプログラム第１回（概要）
粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）
【粒子法】粒子を流体としてレンダリング
３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…
粒子法のシーンを２倍のサイズにしてみた
SPHによる巻き波のシミュレーション２
Haskell、OCamlでSPH法
このあとやりたいこと
カメラ位置を変えて流体をレンダリング
固液連成シミュレーション

■剛体シミュレーションに関するエントリ
粒子ベース剛体シミュレーション（プレビュー）
粒子ベース多体衝突シミュレーション
引き続き、粒子ベース剛体シミュレーション

粒子法(SPH)のプログラム

粒子法(SPH)のプログラムを解説したシリーズです。ソースコードも公開しています。

粒子法(SPH)のプログラム一覧

動画

シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。

動画の一覧

--

Yahoo! BOSS

米Yahoo!のサービスに、BOSSというのがある。BOSSはウェブ検索サービスを構築するためのオープンプラットフォームで、APIを通じてサードパーティーが自由にYahoo! Searchのインデックスやランキングアルゴリズムなどの技術を使うことができる。

TechCrunchによると、今年の５月現在でBOSSは1日あたり3000万の検索クエリを扱っており、さらに利用量が伸びている。この規模は、Bingと互角かやや抜いているレベルの規模だ。

BOSSを採用している企業には、セマンティック検索エンジンを手掛けるhaikaやソーシャルサービスのリアルタイム検索エンジンOneRiot、より理解しやすい検索結果の提供を目指すCluuzといったものがある。TechCrunchのpower searchにも使われているようだ。

BOSSを使えば大きな投資資金がなくても検索エンジンを作ることができる。その計算資源を使って、何かビジネスができないだろうか？

Haskellでレイトレーシング（第５回）〜タプルを返せる

Haskellでは、関数の戻り値としてタプルを返すことができます。

今回は、タプルを返す場合についてCommon LispとHaskellでの書き方の違いを書きます。

Common Lispでタプルを返す場合は、以下のようにvalues関数で値を返して、multiple-value-bindでその値を受け取ります。

(defun sendray (pt xr yr zr)
  (multiple-value-bind (s int) (first-hit pt xr yr zr)
    (if s
        (* (lambert s int xr yr zr) (surface-color s))
        0)))

(defun first-hit (pt xr yr zr)
  (let ...略...
    (values surface hit)))

first-hitがレイと交差した面とその交点の組を返す関数です。それをsendrayから呼び出して、multiple-value-bindで受け取っています。

一方、Haskellではこのように書きます。

sendray :: [Surface] -> Light -> Point ->
           Double -> Double -> Double -> Color
sendray world light pt xr yr zr =
  case first_hit world pt xr yr zr of
    Just (s@(Sphere (Color col) _ _ _), int)
      -> Color ((lambert s int light) * col)
    Nothing
      -> Color 0.0

first_hit :: [Surface] -> Point ->
             Double -> Double -> Double ->
             Maybe (Surface, Point)
first_hit world pt xr yr zr = nearest (map hit world)
  where nearest :: [Maybe (Double, Surface, Point)]
                   -> Maybe (Surface, Point)
        nearest hits = case foldr cmp Nothing hits of
                         Just (d, s, h) -> Just (s, h)
                         otherwise      -> Nothing
     
        ...略...

first_hitでMaybe (Surface, Point)型を返し、それをsendrayでパターンマッチで受け取っています。Haskellでは、タプルを返すのに関数をはさむ必要はなく、そのまま書くことができます。

次回は、「5.衝突点の計算」について書きます。


■他の記事
Haskellでレイトレーシング（第１回）〜導入
Haskellでレイトレーシング（第２回）〜ループは使わない
Haskellでレイトレーシング（第３回）〜型による分岐にはパターンマッチを用いる
Haskellでレイトレーシング（第４回）〜NilはMaybeモナドで表現する
Haskellでレイトレーシング（第５回）〜タプルを返せる

Haskellでレイトレーシング（第４回）〜NilはMaybeモナドで表現する

今回は、関数の戻り値に正規の値とそうでない値がある場合に、Haskellではどのように表現するのかを書きます。

関数が正規の結果とそうでない結果を返す場合の例として、視点から投げられたレイと空間中の物体の交点を返す関数を考えます。

この関数ではレイと物体の交点があるときにはその交点座標を返します。では、レイと物体が交差しない場合には何を返すのでしょうか？

Common Lispのコードでは、レイと物体が交差しない場合には、（暗黙的に）Nilを返すようになっています。

(defun sphere-intersect (s pt xr yr zr)
  (let* ...略...
    (if n
        (make-point :x (+ (x pt) (* n xr))
                    :y (+ (y pt) (* n yr))
                    :z (+ (z pt) (* n zr))))))

レイと物体の交点を求めるための方程式の解nがある場合はmake-pointで作った交点を返し、解nがない場合には暗黙的にNilを返しています。

実用上はこのようにNilを返す形でも問題ないのですが、HaskellではこれをMaybeモナドを使ってより適切に表現することができます。

intersect :: Surface -> Point -> Double -> Double -> Double -> Maybe Point
intersect (Sphere _ r (Point cx cy cz) _) (Point px py pz) xr yr zr =
    let ...略...
    in case minroot a b c of
         Just n  -> Just (Point (px+xr*n) (py+yr*n) (pz+zr*n))
         Nothing -> Nothing

minrootで求めた方程式の解nがある場合にはJust Pointを返し、解nがない場合にはNothingを返しています。

つまり、関数の戻り値をMaybe Point(交点があるかもしれないよ)型と定義し、交点がある場合にはJust Point(交点はコレ！)を返し、交点がない場合にはNothing(交点がなかった…)を返すのです。

いかがでしょう？やりたいことをすごく自然に表現できているのではないでしょうか？

次回は、「4.タプルを返せる」について書きます。

■他の記事
Haskellでレイトレーシング（第１回）〜導入
Haskellでレイトレーシング（第２回）〜ループは使わない
Haskellでレイトレーシング（第３回）〜型による分岐にはパターンマッチを用いる
Haskellでレイトレーシング（第４回）〜NilはMaybeモナドで表現する
Haskellでレイトレーシング（第５回）〜タプルを返せる

Haskellでレイトレーシング（第３回）〜型による分岐にはパターンマッチを用いる

今回は、パターンマッチを使った分岐について書きます。

Common Lispのコードでは、レイとSurfaceの交点を計算するときに、typecase関数を使ってSufaceの型によって処理を分岐させています。分岐用の関数を用意しなければなりません。

(defun intersect (s pr xr yr zr)
 (funcall (typecase s (sphere #'sphere-intersect))
          s pt xr yr zr))

(defun sphere-intersect (s pr xr yr zr)
 ...)

一方、Haskellではパターンマッチを使うことで、型による分岐を自然に表現することができます。

intersect :: Surface -> Point -> Double -> Double -> Double -> Maybe Point
intersect (Sphere _ r (Point cx cy cz)) (Point px py pz) xr yr zr =
   let a = (sq xr) + (sq yr) + (sq zr)
       b = 2 * ((px-cx)*xr + (py-cy)*yr + (pz-cz)*zr)
       c = (sq (px-cx)) + (sq (py-cy)) + (sq (pz-cz)) - (sq r)
       n = minroot a b c
   in case n of
        Just n  -> Just (Point (px+xr*n) (py+yr*n) (pz+zr*n))
        Nothing -> Nothing

今回はSphereのみですが、その他の型を追加した場合にはマッチさせるパターンを追記するだけです。

このように、Haskellでは引数のパターンによって分岐することができるのです。

次回は、「3.NilはMaybeモナドで表現する」について書きます。

■他の記事
Haskellでレイトレーシング（第１回）〜導入
Haskellでレイトレーシング（第２回）〜ループは使わない
Haskellでレイトレーシング（第３回）〜型による分岐にはパターンマッチを用いる
Haskellでレイトレーシング（第４回）〜NilはMaybeモナドで表現する
Haskellでレイトレーシング（第５回）〜タプルを返せる

Haskellでレイトレーシング（第２回）〜ループは使わない

いま、Haskellでレイトレーサを書いています。

第１回のブログで、Common LispのコードをHaskellで書くときにポイントとなったところをあげました。

今回は、そのポイントのうち、「1.ループを使わない」について書きます。

ループを使わない

各ピクセルごとに色を計算するために、Common Lispのコードでは以下のように手続き的にループを回しています。

(defun tracer (pathname &optional (res 1))
 (with-open-file (p pathname :direction :output)
   (format p "P2 ~A ~A 255" (* res 100) (* res 100))
     (let ((inc (/ res)))
       (do ((y -50 (+ y inc)))
           ((< (- 50 y) inc))
         (do ((x -50 (+ x inc)))
             ((< (- 50 x) inc))
           (print (color-at x y) p)))

一方、Haskellでは手続き的なループは使わないので、このような書き方はしません。

それではどのように書くのかというと、以下のようにリストで表現するのです。

tracer :: FilePath -> [Surface] -> Int -> IO()
tracer pathname world res =
   let
       cols = [color_at world (pos i) (pos j)  | j <- [0..(res*100)-1],                                                   
                                                 i <- [0..(res*100)-1]]
   in
       do putPGM pathname (res*100) (res*100) cols
   where
         putPGM :: FilePath -> Int -> Int -> [Int] -> IO()
         putPGM pathname w h cols =
             let header = "P2 " ++ show w ++ " " ++ show h ++ " 255\n"
                 body   = foldr (++) "" [(show c) ++ "\n" | c <- cols]
             in do writeFile pathname (header ++ body)
         
         pos :: Int -> Double
         pos i = (fromIntegral i) / (fromIntegral res) - 50.0

putPGMがpgm画像を出力するアクションで、それにピクセルの色のリストcolsを渡すのです。

このような書き方をすると、「すべてのピクセルをリストにしてから処理するなんて、メモリを無駄に使いすぎるのではないか」と思ってしまいます。しかし、Haskellは遅延評価なので、必要な要素から順々に計算されていくため、リスト全体をメモリに持っておくということは起きません。

このように、リストと遅延評価を使うことで、手続き型のループと同様の処理をHaskellで表現することができるのです。

次回は、「2.型による分岐にはパターンマッチを用いる」について書きます。

■他の記事
Haskellでレイトレーシング（第１回）〜導入
Haskellでレイトレーシング（第２回）〜ループは使わない
Haskellでレイトレーシング（第３回）〜型による分岐にはパターンマッチを用いる
Haskellでレイトレーシング（第４回）〜NilはMaybeモナドで表現する
Haskellでレイトレーシング（第５回）〜タプルを返せる