あー、もう泣きそう。胸が締めつけられる。
死とひきかえに誠人に愛される女性になろうとする、静流が健気すぎる。
2009年10月5日月曜日
2009年9月30日水曜日
SPHによる巻き波のシミュレーション第３回
前回のエントリで、巻き波のシミュレーションができた。今回は、それを綺麗にレンダリングしてみる。まず、粒子でレンダリングしたものがこれ。
これを、マーチングキューブで陰関数曲面を抽出してレンダリングするとこうなる。
さらに、屈折や反射を考慮してレンダリングしたものがこれ。
いい感じ！次は何をしようかな？
■関連するエントリ
SPHによる巻き波のシミュレーション第１回
SPHによる巻き波のシミュレーション第２回
-
2009年9月27日日曜日
SPHによる巻き波のシミュレーション第２回
前回のエントリでは、可動壁によって波を起こすことで巻き波のシミュレーションをしようとしたが、うまく巻き波が起きなかった。そこで、巻き波を起こすにはどうしたらいいのかいろいろと調べてみたところ、RadovitzkyとOritzの"Lagrangian finite element analytics of newtonian fluid flows."という論文を見つけた。
この論文は、有限要素法によってラグランジュ的に自由表面をシミュレーションする手法を提示している。その適用例の１つとして砕波を扱っており、そこに孤立波の理論解が示されていた。これを適用すれば波の岸に打ち寄せる波を再現することができそう。
孤立波の初期条件は、水面の高さをη、x方向の初速をu、y方向の初速をvとすると、
のようになる（x方向のオフセットを省略）。ここで、gは重力加速度で9.8[m/s^2]、Hとdはそれぞれ基準となる水面の高さと波の高さである。
また、tanhとsechは双曲線関数で、その定義は
である。
この初期条件を適用して、またシミュレーションしてみると、
お！巻き波が起きた！
シミュレーション結果を見ると、左から進んできた孤立波が、斜面にさしかかったところでせり上がって、巻き波になってくずれる様子が見てとれる。また、巻き波が打ち寄せたあとにさーっと引いていく様子も、実際の砂浜とよく似ている。
いい感じに巻き波が起きるようになったので、次は、このエントリのように綺麗にレンダリングしてみます。
■関連するエントリ
SPHによる巻き波のシミュレーション第１回
SPHによる巻き波のシミュレーション第３回
-
2009年9月16日水曜日
SPHによる巻き波のシミュレーション第１回
波の一種に巻き波というものがある。浅海域に進入した波が対称性を失って、前面が切り立った壁のようになり、頂部の水が前に飛び出す形で巻くように砕けるものをいう。
粒子法では、しぶきがあがるような激しい変形を表現することができる。そこで、巻き波をシミュレーションしたいと思った。
左側の境界条件を周期的に動かすことで波を起こしている。また、真ん中から右は斜面になっている。
ところが、結果をみてみるとうまくいかない…
波が減衰して斜面まで届かないのと、斜面の底部分で不自然に左向きの流れが起きているのが気になる。
波が減衰してしまう点については、容器の形や可動壁の動きを調整しようと思う。
斜面の底部分で不自然な流れが起きている点については、理由がよくわからない…境界条件にはペナルティ法を使っているが、そのせいで斜面方向に変に力がかかっているんだろうか？
もうすこし調整が必要そうだ。
■関連するエントリ
SPHによる巻き波のシミュレーション第２回
SPHによる巻き波のシミュレーション第３回
-
2009年8月29日土曜日
【論文】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を宇宙から地球上に持ってくるわけですな。
■関連する記事
粒子法のプログラム第１回（概要）
粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）
【粒子法】粒子を流体としてレンダリング
３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…
2009年8月28日金曜日
ClickTaleの無料利用分を使い切ってしまった
おととい入れてみたClickTale。
２日で1週間分のレコーディング数を使い切ってしまった…。これじゃ使えないなぁ。
１人１人のユーザの動きを追える素晴らしいツールなんだけど、月$99は払えないや。
自分で作ろうか？？
2009年8月27日木曜日
ClickTaleの再生がAnalyticsに計測されてしまわないか？
このブログにClickTaleを入れてみた。
ClickTaleはユーザのマウスの動きをレコーディングしてくれるアクセス解析ツールで、レコーディングされたマウスの動きを再生できる。
再生機能で１つ気になったのが、再生したときのページ読み込みが解析対象ページに埋まっているGoogle Analyticsで計測されてしまわないかという点。
Analyticsで確認してみると、ClickTaleからのアクセスは計測されないみたいでよかった。
2009年8月25日火曜日
OneApp
MicrosoftがOneAppという携帯電話向けのアプリケーションフレームワークを発表した。
OneAppは、途上国の人々がTwitterやFacebookのようなサービスに携帯電話からアクセスできるようにするサービスで、プロセッサが非力でメモリが少なくてもJavaが動きさえすれば使うことができる。
大型Webアプリケーションを動かすために必要な処理の一部をMicrosoftのサーバがクラウドとして代行する仕組みで、デベロッパはJavaScriptとXMLによる比較的簡単な作業でアプリをOneAppに対応させることができる。
Microsoftが途上国向けにサービスを提供する理由には、人道的理由もあるが、途上国が潜在的に巨大な市場であるというビジネス上のものもあるだろう。
iTunesでムービーのサムネイルを変えるには
iTunesのムービーのサムネイルは好きなタイムラインのものに変えることができる。
▼変え方
- サムネイルにしたいタイミングでムービーを停止する
- ムービーを右クリック
- 「ポスターフレームを設定」を選ぶ
これで好きなカットをサムネイルにできる♪
2009年8月24日月曜日
粒子法のシーンを２倍のサイズにしてみた
「粒子法のシーンを２倍のサイズにしてみたが…」で失敗したものをもう一度計算した。
【２倍のサイズ】
【元のサイズ】
初期配置や境界の形はまったく同じにして、xyz各方向に２倍している。一見同じように見えるけど、床のチェック模様の細かさが違うところに注目して欲しい。
２倍サイズにすると、より高いところからさらに多くの流体が落下するため、盛大なしぶきがあがった。上のレンダリング結果だと黒くなってしまってわかりづらいが、陰関数曲面をそのままレンダリングするとその様子がよくわかる。
粒子法(SPH)のプログラム
粒子法(SPH)のプログラムを解説したシリーズです。ソースコードも公開しています。動画
シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。--
2009年8月23日日曜日
【論文】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)のプログラム
動画
--
2009年8月20日木曜日
Yahoo! BOSS
米Yahoo!のサービスに、BOSSというのがある。BOSSはウェブ検索サービスを構築するためのオープンプラットフォームで、APIを通じてサードパーティーが自由にYahoo! Searchのインデックスやランキングアルゴリズムなどの技術を使うことができる。
TechCrunchによると、今年の５月現在でBOSSは1日あたり3000万の検索クエリを扱っており、さらに利用量が伸びている。この規模は、Bingと互角かやや抜いているレベルの規模だ。
BOSSを採用している企業には、セマンティック検索エンジンを手掛けるhaikaやソーシャルサービスのリアルタイム検索エンジンOneRiot、より理解しやすい検索結果の提供を目指すCluuzといったものがある。TechCrunchのpower searchにも使われているようだ。
BOSSを使えば大きな投資資金がなくても検索エンジンを作ることができる。その計算資源を使って、何かビジネスができないだろうか？
2009年8月16日日曜日
Haskellでレイトレーシング（第５回）〜タプルを返せる
Haskellでは、関数の戻り値としてタプルを返すことができます。今回は、タプルを返す場合についてCommon LispとHaskellでの書き方の違いを書きます。
Common Lispでタプルを返す場合は、以下のようにvalues関数で値を返して、multiple-value-bindでその値を受け取ります。
(defun sendray (pt xr yr zr)first-hitがレイと交差した面とその交点の組を返す関数です。それをsendrayから呼び出して、multiple-value-bindで受け取っています。
(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)))
一方、Haskellではこのように書きます。
sendray :: [Surface] -> Light -> Point ->first_hitでMaybe (Surface, Point)型を返し、それをsendrayでパターンマッチで受け取っています。Haskellでは、タプルを返すのに関数をはさむ必要はなく、そのまま書くことができます。
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
...略...
次回は、「5.衝突点の計算」について書きます。
■他の記事
Haskellでレイトレーシング（第１回）〜導入
Haskellでレイトレーシング（第２回）〜ループは使わない
Haskellでレイトレーシング（第３回）〜型による分岐にはパターンマッチを用いる
Haskellでレイトレーシング（第４回）〜NilはMaybeモナドで表現する
Haskellでレイトレーシング（第５回）〜タプルを返せる
2009年8月11日火曜日
Haskellでレイトレーシング（第４回）〜NilはMaybeモナドで表現する
今回は、関数の戻り値に正規の値とそうでない値がある場合に、Haskellではどのように表現するのかを書きます。関数が正規の結果とそうでない結果を返す場合の例として、視点から投げられたレイと空間中の物体の交点を返す関数を考えます。
この関数ではレイと物体の交点があるときにはその交点座標を返します。では、レイと物体が交差しない場合には何を返すのでしょうか？
Common Lispのコードでは、レイと物体が交差しない場合には、（暗黙的に）Nilを返すようになっています。
(defun sphere-intersect (s pt xr yr zr)レイと物体の交点を求めるための方程式の解nがある場合はmake-pointで作った交点を返し、解nがない場合には暗黙的にNilを返しています。
(let* ...略...
(if n
(make-point :x (+ (x pt) (* n xr))
:y (+ (y pt) (* n yr))
:z (+ (z pt) (* n zr))))))
実用上はこのようにNilを返す形でも問題ないのですが、HaskellではこれをMaybeモナドを使ってより適切に表現することができます。
intersect :: Surface -> Point -> Double -> Double -> Double -> Maybe Pointminrootで求めた方程式の解nがある場合にはJust Pointを返し、解nがない場合にはNothingを返しています。
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
つまり、関数の戻り値をMaybe Point(交点があるかもしれないよ)型と定義し、交点がある場合にはJust Point(交点はコレ！)を返し、交点がない場合にはNothing(交点がなかった…)を返すのです。
いかがでしょう？やりたいことをすごく自然に表現できているのではないでしょうか？
次回は、「4.タプルを返せる」について書きます。
■他の記事
Haskellでレイトレーシング（第１回）〜導入
Haskellでレイトレーシング（第２回）〜ループは使わない
Haskellでレイトレーシング（第３回）〜型による分岐にはパターンマッチを用いる
Haskellでレイトレーシング（第４回）〜NilはMaybeモナドで表現する
Haskellでレイトレーシング（第５回）〜タプルを返せる
2009年8月2日日曜日
Haskellでレイトレーシング（第３回）〜型による分岐にはパターンマッチを用いる
今回は、パターンマッチを使った分岐について書きます。
Common Lispのコードでは、レイとSurfaceの交点を計算するときに、typecase関数を使ってSufaceの型によって処理を分岐させています。分岐用の関数を用意しなければなりません。
(defun intersect (s pr xr yr zr)一方、Haskellではパターンマッチを使うことで、型による分岐を自然に表現することができます。
(funcall (typecase s (sphere #'sphere-intersect))
s pt xr yr zr))
(defun sphere-intersect (s pr xr yr zr)
...)
intersect :: Surface -> Point -> Double -> Double -> Double -> Maybe Point今回はSphereのみですが、その他の型を追加した場合にはマッチさせるパターンを追記するだけです。
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
このように、Haskellでは引数のパターンによって分岐することができるのです。
次回は、「3.NilはMaybeモナドで表現する」について書きます。
■他の記事
Haskellでレイトレーシング（第１回）〜導入
Haskellでレイトレーシング（第２回）〜ループは使わない
Haskellでレイトレーシング（第３回）〜型による分岐にはパターンマッチを用いる
Haskellでレイトレーシング（第４回）〜NilはMaybeモナドで表現する
Haskellでレイトレーシング（第５回）〜タプルを返せる
2009年8月1日土曜日
Haskellでレイトレーシング（第２回）〜ループは使わない
いま、Haskellでレイトレーサを書いています。
第１回のブログで、Common LispのコードをHaskellで書くときにポイントとなったところをあげました。
今回は、そのポイントのうち、「1.ループを使わない」について書きます。
ループを使わない
各ピクセルごとに色を計算するために、Common Lispのコードでは以下のように手続き的にループを回しています。
(defun tracer (pathname &optional (res 1))一方、Haskellでは手続き的なループは使わないので、このような書き方はしません。
(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)))
それではどのように書くのかというと、以下のようにリストで表現するのです。
tracer :: FilePath -> [Surface] -> Int -> IO()putPGMがpgm画像を出力するアクションで、それにピクセルの色のリストcolsを渡すのです。
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
このような書き方をすると、「すべてのピクセルをリストにしてから処理するなんて、メモリを無駄に使いすぎるのではないか」と思ってしまいます。しかし、Haskellは遅延評価なので、必要な要素から順々に計算されていくため、リスト全体をメモリに持っておくということは起きません。
このように、リストと遅延評価を使うことで、手続き型のループと同様の処理をHaskellで表現することができるのです。
次回は、「2.型による分岐にはパターンマッチを用いる」について書きます。
■他の記事
Haskellでレイトレーシング（第１回）〜導入
Haskellでレイトレーシング（第２回）〜ループは使わない
Haskellでレイトレーシング（第３回）〜型による分岐にはパターンマッチを用いる
Haskellでレイトレーシング（第４回）〜NilはMaybeモナドで表現する
Haskellでレイトレーシング（第５回）〜タプルを返せる
登録:
投稿 (Atom)