Showing posts with label research. Show all posts
Showing posts with label research. Show all posts

2013-04-27

Rhinoceros で、複数のファイル間の layer の色を統一する

高速度ビデオで撮影した動画を1フレームごとに切り出し、1フレームを1ファイルとして background bitmap を Rhinoceros (ver 4) で読ませ、その上で物体の輪郭や特徴的な部位を複数の 2D 曲線としてトレース、という作業を行なっている。

曲線の区別のため、1本を 1 layer として色分けをしている。この色がフレーム(ファイル)によってまちまちなのを統一したいが、layer 数 20・ファイル数 200 以上にもなるため手作業ですべての layer の色を変えるのは大変すぎる。多少楽にする方法がないか考えてみた。

まず、データファイル (Fdata) 以外に、ファイルを2種類用意する。Fcolor, Ftemp と呼ぶことにする。
  • Fdata: 入力済みデータ(曲線群)がある。layer の色がバラバラになっている
  • Fcolor: 色テンプレート用ファイル - 必要な layer 数分の適当なダミーの線を引き、色分けを行っておく。各 layer の名前は Fdata のものに合わせておく
  • Ftemp: 空っぽのファイル。一時避難所として使う
以上3つのウィンドウを常時開いた状態で作業する。手順は、
  1. Fdata の曲線を全て選択し、cut & paste で Ftemp に退避する。非表示 layer や locked layer がないか注意。
  2. Fdata では layer の残骸が残っているので、これを全て消す。Default layer は消せない。
  3. Fcolor のダミー曲線群を Fdata へ copy & paste する。この時点で layer の色が更新される。用済みなのでこれらはすぐに消す。つまり Ctrl-V の直後に Del でよい。
  4. Ftemp に退避してあった曲線を Fdata へ戻して完了。
問題点として、sublayer があった場合にはその構造は破壊されてしまう(上位レイヤが消滅して、全部フラットな関係になる)。解決策は不明。
⇒ 解決ではないが、上位レイヤ名のダミーレイヤを Fcolor に用意しておくことはできる。もちろん結局は手作業で曲線レイヤをサブレイヤ化しないといけないので面倒だが…

UPDATE 2013-12-17
Rhino 5 SR6 を使っているが、上記のサブレイヤ構造が破壊される問題は直っているようだ。
ただ、いま別の問題が…色ではなくレイヤの名前を一気に変えたいのだが、これはもう手作業しか無いのか…

2013-03-12

DippMotionで低解像度の特徴点を扱うときは画像の解像度を増やせ

DippMotion でプロットするときの最低単位は「1画素 (pixel)」です。画像が高解像度だったりして、プロット対象が沢山のpixelから成っている(たとえば特徴点の直径が 30 pixel もある、とか)なら、特に問題はありません。しかし、画像が低解像度だったりして、プロット対象が数 pixel しかないときはこれが問題になりえます。

例:
特徴点の直径がちょうど 4 pixel だった。真ん中のどっちの画素に点を打つべきだろうか…? どちらかを選ばないといけないが、いずれを選んでも「真の中心」からのずれが 0.5 pixel あり、これは直径の約 12% に達する。見た目でも明らかに「ずれてるなぁ」(というか、「偏ってるなぁ」)とわかる。

対策:
irfanview などで全てのフレームの画像を resize または resample し、解像度を 4, 8, 16 倍などにする。たとえば解像度を4倍(縦横2倍ずつ)にすれば、特徴点の直径は 8 pixel となり、真ん中のどちらを選んでも「真の中心」からのずれは直径の約 6% と半減する。8倍(縦横4倍ずつ)にすればさらに半減。当然ながら、calibration の画像も同様に拡大する必要がある。

Drawback:
画像を大きくするほど、画像も DippMotion も重くなるし、あまりにも大きいとそもそも読み込めなくなる。自分の場合は、たしか 4000x4000 程度が限界だった(ビデオカードによるかも)。

どうしても解像度が足りない!という場合は、画像の一部のみを切り出して resize/resample → Dipp し、再構築の前に切り取った分の座標値を(手動で)足してやる、とかは可能と思われる。

2013-01-23

Why "evolution of flight" and only flight?

最近 evolution of flight に興味があるんだけど、そうなるとなんで flying だけ特別扱いするの?そんなことしていいの?とか考えたりしてた。してる。

今のところ key かなと思ってるのは、メカニクスからの要求というか制約が、他のさまざまな淘汰圧に比べて非常に大きいのが飛翔体の特徴だろう、というあたり。Swimming になると、鉛直方向の力という制約が緩くなるので、メカニクス以外のファクタが入り込む余地が増える、はず。

飛ぶ奴はまずみんな当然翼がないといけない。一番わかり易い convergent evolution(といっていいのか?)。あと、軽くないといけない。鳥類は哺乳類に比べて形態の広がりが著しく小さい(みんな似ている)。飛翔に適した基本形/基本スペックから大きく外れると…ダチョウ・キーウィ・etc. のように飛べなくなってしまう。重すぎて飛べなくなったが、鉛直方向の流体力を浮力に切り替えたお陰で泳げるようになったのがペンギン。まぁフリッパーの流体力についての論文まだ読んだこと無いから実は結構フリッパーでも上下に力だしてるかもしんないけど。

昆虫が鳥やコウモリに比べるとバリエーションに富んでいるのは、low Reynolds number なので流線形が重要でないことが効いている。でも軽さの縛りはかなり厳しい。

遊泳の場合、浮力があるので、重さの制約がなくて鉛直上向きの流体力をがんばって作らなくていい。シロナガスクジラ重いよね。ただ、魚とイルカ・クジラは Re 高いので、流れ方向の抗力係数を減らすための流線形ということで形が似てる。この制約はかなり強いんだな。その意味ではメカニクスかなり効いてる。でも一方で、エビとかタコとかイカとかクラゲとか色んなやつがいる。揚力作らなくても浮けるし生きていけるので、バリエーションの余地がある。

というあたりで、flight を特筆することは(evolution の観点からも)正当化できそうなんだけど、逆にこういうのを横断的に見渡してなんかうまいこと論文を書けないかなぁ。本とか review article にはなりそうだけど普通の article にはしづらいかな。まぁ当たり前すぎるかな。

2012-11-19

RhinocerosでBackground Bitmapを使う場合、画像のピクセル数は縦横とも2の累乗でないといけない(横長・縦長はOK)

http://wiki.mcneel.com/rhino/backgroundbitmap
http://wiki.mcneel.com/ja/rhino/backgroundbitmap
に書いてあった(日本語訳は「2の2乗」となっていておかしいが正しくはもちろん累乗)。

ヘルプ大事だね…気づくまでかなり時間を無駄にした。

つまり、1024x1024 や 1024x512 は OK だが、800x600はダメ。
これを守らないとピクセル間が勝手に補間されて、にじむ。

なお、
If this criteria is not met, then the image will not display at all.
ともあるがこれは嘘。そうであればもっと早く気づけたよ…。

2012-10-07

Rhino で線を分割するときの向き

Rhinoで線を分割して点列を作るとき、どっちからどっちの向きに点列が出力されるかは、調べて変えることができる。

面から法線ベクトルが出ているようなアイコンの Analyze Direction/Flip Direction を使う。実際には左クリックの Analyze Direction だけでよい。アイコンを左クリックすると object を選択せよと言われるので、曲線を選んで Enter. すると白い矢印で現状がどちらの向きかが表示される。ここで曲線を左クリックするたびに向きが逆へ変わる。どちらかに決めたら Enter で確定。

実際に向きを変えてから曲線を分割し、点列を Export selected してみれば、向きが変わっていることがわかる。ただし、なぜかはよくわからないが、出力されたテキストファイルでは、一列目が矢印の終点側になっている…つまり、直感とは逆向きに並ぶようだ。慣れれば問題ないかとは思うが。

2012-09-05

複数の画像に対して一度に余白を追加する

現在の作業の都合で、画像に余白を追加する必要がある。具体的には、「800x600 の下に 800x200 を追加して正方形にしたい」とか、「512x512 の周囲に余白を追加して 1024x1024 にしたい」のふたつ。コマンドラインでやるしかないかなぁ、と思っていたら、GUIでできそうなのがいくつかあった。

[2012-10-23 追記]
ふつうに IrfanView でできたわ…。メニューの File > Batch conversion/Rename... から、Output format の右下の Advanced ボタンでダイアログを開き、右下の MISCELLANEOUS の一番上 Canvas size にチェックし、その右の Settings ボタンをクリック。上下左右の余白と色を選択。こっちの方が早いので Irfan 使いならオススメ。
[/追記]

新しそうな Ralpha を最初に試してみたらできたので、これでいいや。
http://nilposoft.info/ralpha/

入力は drag & drop で直感的に。

設定は一瞬戸惑ったけれど、ヘルプ ( http://nilposoft.info/ralpha/help/rotation/ ) 見たらなんとなくわかった。「左、上、右、下のサイズを指定(px)」にして、条件はよくわからないので「ALL」に。その下のところに余白の数字を記入する。
  • 「800x600 の下に 800x200 を追加して正方形にしたい」→「0,0,0,200」
  • 「512x512 の周囲に余白を追加して 1024x1024 にしたい」→「 256, 256, 256, 256」
でいけた。

出力ファイル(群)の保存のされかたは、デフォルトでは「そのファイルがあるディレクトリに resize というディレクトリが作られてその中に同名で保存」になっていた。

試してないけど他にググって出てきたのは
  • XnView - これはたぶんできそう
  • PhotoScape
  • PhotoShifter
いずれも Windows 向けだと思う。

2012-09-04

GIMPでDBPを使って複数の(連番)ファイルにunsharp maskをかける

準備

GIMPをダウンロード
http://www.gimp.org/downloads/

GIMP 用プラグインである DBP (David's Batch Processor) をインストール
http://members.ozemail.com.au/~hodsond/dbp.html
GIMP 2.6 では問題なかったが、最新版の 2.8 でも動くかは知らない。このプラグインの Sharpen は実際には unsharp mask を使っている。ソフトによっては sharpen と unsharp mask は別ものなので、紛らわしい。

DBP の使い方

GIMP でメニューの「フィルタ > Batch process...」から DBP を起動する。

Input タブで Add Files する。このとき、左上のアイコンでアドレスバーを表示させて、エクスプローラのアドレスバーのパスをコピペすると移動が楽(これは FIELDVIEW など他のアプリでも使える tip)。

Sharpen タブで unsharp mask を Enable する。Radius, Amount, Threshold は元画像で試して先に決めておく。いちおう Test ボタンでテストできるが、いちいち Show Images を押さないとやり直せないので若干面倒。プレビューしながらフィルタ適用ということもできない。

Rename タブで出力先ディレクトリと名前を決める。まず Select Dir ボタンから出力先のディレクトリを決める。ソースと同じ場所でよければ Source Dir ボタンでいい。出力するファイル名を変えたければ、Add Prefix か Add Postfix に何かを書いておく。どういう名前になるかはすぐ下の Example のところでシミュレートされている。他にもチェックボックスでグレースケールへの変換とか選べるぽい(Colour タブにもグレースケール変換チェックボックスがあるが…)。

Output タブでは出力ファイルのファイル形式(拡張子)しか選べない。これは Rename タブと統合すべきだと思う……。

準備OKとなったら下の方にある Start ボタンでバッチ処理開始。このとき Show Images トグルボタンを押していると、処理中(処理後?)のファイルがいちいち表示される。押さない方がよさそう。出力先ディレクトリに同名のファイルが既にある場合は、エラーメッセージが出て処理が止まる。

関連エントリ

作業メモ:Phantom -> DippMotionで飛翔体の特徴点から位置情報取得して3次元再構築
Phantom で連番ファイルを保存する方法(ファイルに「+桁数」を付ける)はこちらに。

Rhinocerosで画像からwing gridを作る

2012-06-11

作業メモ:Phantom -> DippMotionで飛翔体の特徴点から位置情報取得して3次元再構築

Phantom でのキャプチャ設定

画像処理

下の方にある Image processing ボタンから(PCCでは上の方にあるパレットのボタン。Image Toolsとある)。gamma が最初1.41になっていることに注意。Set default を押すと 1.0 になる(PCCでは2.2)。

フィルタは、元画像をそのまま見るなら Sharpness を、輪郭を抽出するなら Edge Hipass 3x3 か Edge Laplacian 3x3 がよさそう(いま見てる動画だと、どっちでも同じに見える)。ただ、結局、「人の勘」に頼る場面が出てきたため、シャープネスだけに落ち着いた。もっと鮮明な画像が得られている場合は輪郭抽出だけでいいのかも。

保存

右のほうの Convert から(PCCでは、右のペインでPlayタブ > 下のSave Cineボタン。動画でもこのボタン)。必要な範囲よりやや長めの範囲を切り出す。まずは短く切り出しただけのcine形式で一度保存しておくと、後の作業が楽。

いろいろ試したけど、背景をちゃんと白くしていて照明も安定している(蛍光灯などのチラツキない)自分の場合、brightness と gamma は 0 と 1.0 のままにして、gain だけを飛翔体の輪郭がギリギリ見えるまで上げる、というのがよかった。brightness = 0 & gamma = 1.0 の状態でプリセットを保存しておくと便利。

動画の切り出し

avi にする場合はフォーマット選択ダイアログが出てくるので適当に選ぶ。常識的に考えると、最高画質ということで非圧縮だろうか。プレゼン等で見せるためにはファイルサイズ大きすぎ & 重すぎるので、Cinepak などで適当に圧縮したバージョンも作っておくといいかも。非圧縮AVIを適当なコンバータで変換してもいいし。

静止画の切り出し

ビットマップの静止画を切り出して、GIMPとかIrfanとかでバッチ処理したい場合。連番ファイルにするために、ファイル名に +3 を付ける。たとえば N1_cam1_+3.bmp という名前にすると、 N1_cam1_001.bmp, N1_cam1_002.bmp, ... という連番ファイルが生成される。最初はこっちでやって、ガンマとかアンシャープマスクとかかけてたけど、結局 Phantom 内でのシャープネスだけにしたから今はやってない。

[2012-09-04追記]
と思ったがやっぱり使いそうなので、GIMP でのバッチ処理について別エントリを書いた:
GIMPでDBPを使って複数の(連番)ファイルにunsharp maskをかける

DippMotionPRO での作業メモ

フォルダ関係

いろんなファイルを保存したり読み込むダイアログで、デフォルトが

  • 前回の作業フォルダ
  • C:\Program Files\Ditect\DippMotionProKP
  • C:\Documents and Settings\ユーザ名\Application Data\Ditect\DippMotionProKP

のどれか(あるいはその下のフォルダ)なっていることがある。このせいでリンク切れになって「なんか動画動かないなぁ?」とかよくなる。というかこのへんいじったらさっき画像が全部消滅したよ・・・

そこで、DippMotion を起動したらまずメニューの「オプション > 一時ファイル格納場所指定」から、自分の目の届くところにしておく。

それから、エイリアス(ショートカット)から起動する時の「作業フォルダ」に、いま作業中のフォルダの近くを指定していても、開くダイアログ等で「C:\Program Files\Ditect\DippMotionProKP」に飛ばされることがよくあるので、この場所に作業中フォルダへのショートカットを置いておくとはかどる。同様に、マイドキュメントに飛ばされることもある(いま2D座標入力でそうなった…)、ここにもショートカットを置いてもいいかも。

レンズキャリブレーションについて

全部の撮影で同じ組み合わせのはずなので、一度作ったファイル(カメラに対応した3つの .dxc ファイル)を使い回せばいい。理想的には、キャリブレーションフレーム・被写体を撮影したままの条件で、つまり撮影ごとに毎回ドットパタンを撮影するのがよさそう。

既存のレンズキャリブレーションファイルを使いまわす場合は、Dippを開いた初期画面の右にあるアウトラインペインで、それぞれのカメラを展開してから、「キャリブレーションデータ」を右クリックして「ロード」。なぜかうまくロードできない場合があるが、アプリを再起動してみたり、上記の作業フォルダを変えたりしたら直った。

3Dキャリブレーションについて

キャリブーレションパラメタの計算を行ったら、条件が同じ複数回の撮影に使いまわすために、保存しておく。これはカメラごとなので注意。計算方法について、説明書には漸近法の方が良さそうに書いてあるが、少なくともこれまでに行ったケースでは線形法の方が明らかに良かった。計算は線形法が圧倒的に速い。両方試して、よりよさそうな方を使う。

UPDATE 2014-09-06
Win 7 64 bit だからなのか、いつのまにか DLT の係数がプロジェクト保存しても .dxc ファイルにかかれなくなっている (Dipp Motion Pro 2.24d)。どれか、バイナリに保存してるのだろうか?ともあれ、解決策は2つ。ひとつは、3Dキャリブレーションを実行した後で、「オプション > カメラパラメータ確認」の画面を見ること。でも選択できないので手で書き写すことになる。面倒。もう一つは、作業用フォルダを見ること。3Dキャリブレーション実行直後なら、AnaDataCh1 などの中に TmpDXCamParam.dxc というのができており、ここには従来通り [MATRIX] T1=... などがあるので、これを保全すればよい。
と思ったら普通に保存で表示されるようになった。相変わらず挙動がヘンだ…

2Dデータ入力について

基本操作

左手:
  • 矢印の左右キーでコマ送り/戻し
  • 矢印の上下キーで部位移動

右手:
  • マウス左クリックで点を打つ
  • マウス右クリックでキャンセル・終了
  • マウスホイールで画像の上下移動(拡大中)

カメラを開いてからまずやること

-画面下の「画像フィット」チェックボックスを切り、メニューの「表示 > 拡大」から 4x4 あたりを選択
-メニューの「属性 > 座標入力方向の指定」で、デフォルトが「カラム単位」(縦に進む)になっているのを「部位単位」(横に進む)に変更
残念ながら、こういうインタフェース部分の設定は保存されないので、カメラを開くたびに毎回やる必要がある。なぜかカーソルの設定だけは保存されるようだが…。
-「編集 > 部位数の設定」から、部位数を多めに確保する(下記参照)

重要:保存について

2Dデータ入力画面から保存する場合、デフォルトのフォルダは「一時ファイル格納場所指定」の下のフォルダになっている(名前が同じなので騙される)。自分のプロジェクト下のフォルダに保存されるのは、一度2Dデータ入力画面を閉じてホーム画面に戻り、プロジェクトの保存を行った瞬間。もしプロジェクトの保存を行わずにDippを終了すると一時ファイルの方にしか保存されていない状態になる。このままもし他のプロジェクトを開くと、データは全て吹っ飛ぶ。「ホーム画面に戻ってプロジェクトを保存」は癖にしたほうがよさそう(というか、一時ファイルでなく直接保存すりゃいいのか?→と思ったけど、複数のプロジェクトをい行ったり来たりしてる場合は不便なので、やはり「プロジェクトで保存」が現実的)。

重要:部位の数と命名について

メニューの「属性 > 最大表示セル数」から、「部位方向(ライン)」の数を必ず現在の部位数よりも多くしてから、左下の部位名の欄をダブルクリックして命名する。もし、「部位方向(ライン)」の数が現在の部位数よりも少ない場合は、正しく命名できないだけでなく、5つ上の部位名が上書きされてひどいことになる。謎仕様というかどうみてもバグ。とにかく部位数は早めに増やしておくこと。ただ、あまりに多いとアプリの動作が重くなる。カメラ1, 2, 3, ... で共通の部位名を当然使いたいので、早めに使う点を全部決めて、いずれかのカメラの2D座標入力画面において、「ファイル > 部位名ファイル...」から「保存」で部位名ファイル(拡張子 .bui)を保存する。他の2D入力画面で「ファイル > 部位名ファイル...」から「開く」。このとき、既存の部位名は問答無用で上書きされた気がする。

重要:左下の「追尾」チェックボックスについて

修正追尾での挙動はまだハッキリわからないが、ここにチェックが付いているものを優先的に(?)追尾していくようで、間違うと非常に危険。これから打っていく点のみにチェックが付いていることを常に確認すること。ミスるとセーブした時点からやりおなしになる(ので、少なくとも部位一つ終わるごとにセーブ必須)。

チョウの翅の模様追尾について。サインペンなどで打った点を追尾するのとはわけが違う。まずは何度か再生してみて、どういう方向に飛ぶかをつかむ。慣れると自動追尾ができるかどうかだいたいわかる。全部手動でやると、却ってガタガタになると思われる。面倒なようでも自動追尾を併用したセミオート的なのが結局はよさそう。 ←これはなんともいえない。相関係数下げた場合には自動でも揺らいだり、deviation 載ったり(とりたい点よりも「ずっと左に1ドットずれたまま追尾」とか)することはある。

常時修正追尾で行う。右手はマウス、左手は矢印の左右キー。通常は自動でやって、追尾できない部分はどんどんとばして、とにかく最後までいく。それから必要に応じて、追尾方向を逆(巻き戻し再生)にしたり、手動で打ったりする。あと、テクニックとして、「1ピクセルだけ間違った点に打ったかも」というときは、マーカ自体のせいで見えないので困るが、いちいち削除するのは時間がかかるので、一旦手動に切替えて画面上の適当なところをクリック、で再度自動に変えて正しい点をとり直すと若干早い。

自動での、相関の閾値は取る点の目的と要求される正確さに応じて変えるが、基本的には 0.9 で初めて、進まなくなったら 0.85 とか 0.8 まで下げる。探索範囲については、XY方向へのベクトル的なものを指定する「指定方向」を推奨。ただし、追尾点の移動方向に合わせてときどき微調整が必要。追尾範囲(赤の四角)は小さめにしたほうが高い相関値でもうまく追尾できることが多いようだ。

拡大表示で打ちたいところに打つためには、カーソルの中心を打ちたいピクセルに完全に重ねるのではなく、0.5ピクセル程度(?)左上にずらすように打つ必要がある、ようだ(キャリブレーション時は逆に右下にずらす必要がある…なんなのこれ…)。

ポイント表示モードについて。最初は「全ポイント表示」で打っていく。混んできて見づらくなったら「ポイント表示」も併用。「軌跡表示」も、ガタつきを見つけるのに便利。ある程度打ったら切り替えてチェックするといい。

削除

メニューの「編集 > ポイントの部分削除」から。このソフトウェアにUndo 機能は存在しない(これ…商用ソフトなんだぜ…)ので、削除は慎重に。間違っても上の方の「座標データの全削除」や「部位の削除」をやってはいけない。特に後者は「ポイントの部分削除」と近いので要注意。たぶん誰もが一度はやると思う。けっこう泣ける。メニューのUI作った人に軽く殺意を覚える。というかメニューだけでなく全体に(ry

3D再構築したとき

「すべて」で表示したときに、なぜか2Dとの同期再生ができなかった。2D解析結果を一つ(カメラ1でよさそう)開いて、めんどうだけどメニューから画像と2D座標ファイルを読み込む。それから閉じてまた「すべて」を開くと、なぜか他のカメラも含めてうまく同期再生できてた。謎だ・・・。

カメラ2台以上で追えていない(入力していない)時刻があったら、その時刻は点が消える。変なところに飛んでいたら、キャリブレーションがおかしいかもしれない。あと、そもそも同期撮影した画像でなく、カメラごとに別のケースを使っていたという場合も…(「すべて」表示にすれば変だとすぐわかるはずだが、「3Dのみ」にしてると気づきづらい)。

2011-12-18

ストーリ思いついた

FCT 関係で Biol. Lett. にいけそうなストーリ思いついた。2つの不安定性があってそのどちらかが卓越すると不安定になる、という感じ。検証はこれからするというのに、なんか直感的にこの方向でいけるという気がしてならない。

問題はコードの実装に手こずってるってこと。開発のモチベーションは高まったんだけども。

2011-12-16

プレゼンで気をつけたいこと

Tipというか基本なんだけど。

学会発表だけでなくゼミやミーティングでのプレゼンというかスライドによる発表も想定している。

  • スライドにはビジュアルを詰め込む。文字は、結局はすべて口で言えるので極力減らす。
    • 特に学会発表では限界まで減らす。
    • ゼミでは、持ち時間にもよるが多くてもいい/多いほうがいい。
  • グラフを見せたら必ず縦軸と横軸が何かを言う。
    • 最も基本的なことなのにできてない人が非常に多い。特にゼミで多い。
  • 数値計算をした場合、主要な計算条件を示す。
    • 何が主要かというのは audience による。たとえば毎回同じ顔ぶれのゼミで、先週の発表とほとんど同じ計算をした、というのであれば前回と違う点についてのみ触れればよい。前回とかなり違うことをしたり、初めての計算を発表するのであれば、xyz軸がどっちか、とかの基本から丁寧に説明すべき。自分で計算したなら当然条件は言えるだろうし、誰かに手伝ってもらった場合でも、自分の研究の総責任者は当然自分なので、条件は把握しておかないといけない。
    • 学会発表は時間制限があるので話がかなり変わってくる。CFDの場合であれば無次元数(Reとか)とそれに使った代表長さ・代表速度、境界条件あたりがだいたい必須で、計算領域のサイズ、メッシュ数、初期条件なども必要に応じて入れる感じになるが、持ち時間と自分のストーリによってはかなり大胆に端折ることも重要。しかたなく削る、のではなく、あえて削る。生物屋さんに「Navier-Stokes equations の離散化が有限差分法でも有限要素法でもなく、有限体積法で、しかも strong conservation form と pseudo-compressibility を使っているんですよ!」なんて言っても相手寝ちゃう。流体計算について生物屋が知りたいのは、まず「2Dなのか3Dなのか」「blade element method のような簡易的なモデルなのか、ガチのCFDなのか」というところなので、たとえば「3次元(非圧縮の)Navier-Stokes equations を解いている」くらいでいい。解いている対象が明らかなら、「非圧縮の」は入れなくてもいいくらい。スライドの文字は削って口で言う、のが理想かな。
まだ増えそうだけどとりあえずここまで。

2011-12-15

修士院生に大切なこと

ってタイトルにしたけど他の立場の研究者にもある程度は当てはまるかも。でもいま自分がはまってるのが後輩の修士の面倒の見方をちょっとまずったなぁということなので修士を念頭に。

具体的な goal の設定と、それを常に意識することが必要。具体的なってのを具体的にいうと、

「こういう背景・目的があって、こういう実験・計算をしました。こういう結果が得られたので、こういうことが言えます」

という短いストーリが書けること。結果まだ出てないって?じゃあそこは(仮)でいいから。

「こういう背景・目的があって、こういう実験・計算をします。こういう結果が得られたとして、こういうことが言えます」

これなら書ける?書けない?じゃあ最終的に修論書けないよね。

特に大事なのは、「こういうことが言えます」の部分かも。それ、目的に対応してる?自分も含め、そもそも「こういう結果がでました」で終わりで、「こういうことが言えます」まですぐ出てこない人も多いんじゃないかな。でもそれって端的に言うと「なんでそんな研究したの?」ってことになるんだよね。「いやー、なんか面白い結果が出るかと思って」としか返せない感じ。それは卒論までだよね、と言われるかも。いや、実際に面白い結果が出ればそれでいいっちゃいいんだけど、そこは後付でいいから「実はこういう目的に対応してます(こういう嬉しいことがあります)!」ってのを自信を持って言えなきゃ(書けなきゃ)いけない。

時期的には、M1の4月に設定するのは無理だろうから、その時点ではテーマの大枠でもいいと思う。なので、たとえばM1の夏休みとか、遅くとも中間発表までにはこれを書く。書いたらすぐ先生にどうですかってメールして、必要なら練り直す。

で、この goal をずっと意識することも重要。「どういうストーリで修論書こうかなぁ」ってのを、単純作業以外の時間は意識しておいて、「あれ?このストーリは微妙かな」ってなったら修正。「Aという実験をしました」はちょっとできそうにないので、代わりにBという実験をしよう。あ、でもそうするとこれ結果出ても当初の目的満たしてなくね?じゃあ目的修正するか?みたいな。

----
面倒を見る(固く言うと「指導する」)立場からすると、どこまで口をだすか、どんなことに口をだすか、が難しい。さらに作業内容の分担も難しい。

(ミカン(おいしい))

2011-02-11

FCT (flapping counter torque) について

今のところ押さえるべき主要な papers は3つ。

Cheng et al. (2009) IEEE conference paper
Hedrick et al. (2009) Science
Cheng et al. (2010) JEB


ここで、Cheng et al. と Ty Hedrick の大きな違いは、
  • Cheng = Constant angular velocity simulation
  • Hedrick = 実際の飛翔体の計測値と、数式による予測(←ここは理解イマイチ)
ということ。だと思う。でないと Cheng の方あんなに綺麗にトルクでるはずがない。
正直 constant でやるのに意味あんのかよくわからんけど…。

なお関連する文献として、
Fry et al. 2003 Science
Ramamurti & Sandberg 2007 JEB
がある。



追記 (2011-12-13):

上の const. ang. vel. というのは多分勘違い。Body ang. vel. を入れた wing kinematics でやってる。問題は Ramamurti & Sandberg の結果を active torque w/o body rotation としている(ようにしか読めない)ということ。実際は入ってる。

Hedrick (2011) JEB が出ている。Commentary という要は review paper. 先輩に教えてもらった。

なおイギリスのO大方面では FCT にはわりと懐疑的らしい (personal communication at a pub : p)

ところでこんなのここに日本語で書いても読む人いないよね(それこそ先輩くらいじゃん)

2010-12-07

Rhinocerosで画像からwing gridを作る

準備

元画像(たとえば翅を広げた蝶を背中から見た平面図、など)を用意する。フォーマットはビットマップ (.bmp)。

前処理

たとえば、Irfan Viewで必要な範囲だけ選択 → crop (Ctrl+Y) し、Ctrl+Eでフィルタ選択画面へ。Find Edgesをlevel 5とかでかけてエッジを強調する。Sを押して別名で保存 → ビットマップ (.bmp) で保存。

あるいは、GIMPでDBPを使って複数の(連番)ファイルにunsharp maskをかける など。

Rhinoceros での作業

色の設定

白黒画像に対して黒い線は見づらいので、まず色を変える。メニューの File > Properties... から、

Rhino Options > Appearance > Colors > Selected objects
の右の色をクリックして変更。

Rhino Options > Appearance > Advanced Settings > Wireframe > Objects > Control polygon color
を Use fixed color にし、すぐ下の色をクリックして変更。

Rhino Options > Appearance > Advanced Settings > Wireframe > Objects > Curves > Curve color usage
を Use single color for all curves にし、すぐ下の色をクリックして変更。

グリッドの設定

後で書く。

基本操作

Rhinocerosを起動。Top とある左上のペイン(上面図)をクリック。メニューから「View > Maximize」で最大化。マウスホイールで視点を拡大縮小できる。キャンバスが全部見えるくらいまで縮小しておく。ウィンドウ下の Snap をクリックしておく(グリッドにスナップする)。右ドラッグで画面移動できるのを覚えておく。また、Spaceキーで直前の動作を再実行できる。これを知らないと死ぬ。

画像の配置

メニューの「View > Background Bitmap > Place」で画像を選ぶ。すると上の方のコマンドバー(?)に「画像を置く四角の始点を指定しろ」と出る。始点をキャンバス
左上中央、原点に指定し、終点は
画像が最大になるようにとるキャンバスの右端になるように、右下に向かって図を拡大する。こうすると1ピクセルが1mmに対応して、DippMotionなど他のソフトとの連携がしやすくなる、はず。ビットマップが(グレースケールで)置かれる。必要なら「View > Background Bitmap > Grayscale」を外せばカラーになる。

プロット

メニューの Curve を多用する。

外形のみ抽出する場合

TBD.

芯線と半径を取得する場合 (UPDATE 2014-11-28)

レイヤーを適宜利用すると楽。
  1. 外形を Control Point Curve で適当にトレースする。点数はいくら多くてもいいので、それなりにリアルに引いておくとよい
  2. 基準線を引く。下記の垂線が取れるように。ただし、CFDするのだから前縁後縁をどう配置するかも考えながら決める(この取り方次第ではうまくいかなくて、やり直しもあり得る)
  3. 基準線を適当な数に分割。60とか。これは目安
  4. 基準線から、垂線を適当に一本引く
  5. (オプショナル)背景を隠す
  6. 垂線を、両側に十分長く(どこに平行移動しても外形線を超えるように)延長する
  7. 垂線を、offset (平行移動コピー)しまくる。コマンドバーの"Distance="の部分をクリックすると距離を変えられるので使う(コマンドでdをタイプしてもいい)。Spaceキーあるいは右クリックによる (redo) も当然利用。このとき、基準線上の点を目安に使う。厳密に従う必要はない
  8. 基準線と、基準線上の点を別レイヤーにして不可視化しておく
  9. 垂線の、外形線からはみ出た部分を、外形線を利用して trim する
  10. (オプショナル)外形の曲率が大きな箇所(直線的な変化でない箇所)に垂線をたくさんひいたりして手直し
  11. 外形線をレイヤーで不可視化
  12. 残った直径にあたる線を全選択し、 Curve > Point Object >Divide Curve by > Number of segments で MarkEnds=Yes にして 2 で分割。これで一気に端点と中点が得られる
  13. 点を全て選択して「Export selected > .txt」すれば座標データが得られる。あるいは芯線データと外径点を分けて export してもよい

2010-11-23

definitions of wing areas for aeroplanes and insects

特に意味はないんだけどさっき共同研究者に訊かれたので僕らはこういう定義使ってるよっていう。

飛行機の翼面積(主翼面積) S_w はこういう感じ。
上から見た投影面積で、胴体部分も含める。
C_L = L / ( 0.5 * rho * U_ref^2 * S_w ) みたいに無次元化に使う。


昆虫は胴体とrigidに繋がってないから、こう。
こうというか要は一枚の翅の片面の面積。キャンバとかギザギザは考慮してない。
C_f = F / ( 0.5 * rho * U_ref^2 * S_w1 ) みたいに使ってる。
F は空気力ね。

2010-09-06

ParaViewで可視化

2011-12-15追記:

なぜか「ParaView 可視化」でググったらこのエントリが3番目に表示される謎。中身ほとんどないのでググって来た人はそっとタブを閉じて次以降の検索結果見た方が幸せになれるよ。

そういえば ParaView の日本語の本も買ったんだけどまぁ予想通り自分には意味ない感じだった。素直に英語の wiki とかマニュアル読んだ方がいい。
----


データ形式について
Plot3D は読めるが、xyz と q のみで、function に対応してない(追記:そうでもないかも?)ので使い勝手が悪い。Tecplot は読めないいつからか読めるようになっていた (2014-05-27)。標準の VTK classic (Legacy VTK?) 形式で作るのが楽そう(書き方は→ http://visitusers.org/index.php?title=ASCII_VTK_Files ←ParaViewのライバルであるVisItのwikiだったw)。

やってみた。



流れが左から当たっていて、渦度の画面に垂直な軸成分を表示している(って言われなくてもわかる人しかこの記事読んでないだろうけど)。似たような条件について FIELDVIEW で作ったのは→ http://twitpic.com/2jc680

アニメの作り方
単に時系列に連番のファイルを用意するだけ。
  • data000.vtk
  • data001.vtk
  • ...
と用意してやれば、ParaView で開こうとしたときに「data..vtk」という連番ファイルとして自動認識してくれる。あとは適当に見た目を調節して(後で書く予定)、File > Save Animation でOK.

VTKファイル出力用Fortran 90ソース

  !! ファイルの出力
  write(filename,"('data',1i3.3,'.vtk')") n_file !! 出力ファイル番号をファイル名に書き込む
  open(200,file=filename,status="unknown",form="formatted",position="rewind")
  
    write(200,"('# vtk DataFile Version 3.0')")
    write(200,"('2D flow')")
    write(200,"('ASCII ')")
    
    write(200,"('DATASET STRUCTURED_GRID')")
    write(200,"('DIMENSIONS ',3(1x,i3))") im, jm, 1
    
    write(200,"('POINTS ',i9,' float')") im*jm
    do j=1,jm
    do i=1,im
      write(200,"(3(f9.4,1x))") x_plt(i), y_plt(j), 0.0d0
    enddo
    enddo
    
    write(200,"('POINT_DATA ',i9)") im*jm
    
    !! velocity vector
    write(200,"('VECTORS velocity float')")
    do j=1,jm
    do i=1,im
      write(200,"(3(f9.4,1x))") u_plt(i,j), v_plt(i,j), 0.0d0
    enddo
    enddo
    
    !! pressure
    write(200,"('SCALARS pressure float')")
    write(200,"('LOOKUP_TABLE default')")
    do j=1,jm
    do i=1,im
      write(200,"(f9.4)") p(i,j)
    enddo
    enddo
    
    !! z_vorticity
    write(200,"('SCALARS z_vorticity float')")
    write(200,"('LOOKUP_TABLE default')")
    do j=1,jm
    do i=1,im
      write(200,"(f9.4)") z_vorticity(i,j)
    enddo
    enddo
    
  close(200)


いろいろ適当だけど気にスンナ。do ループが外側から i, j でなく j, i の順になっているのは、i と j の数が違うときに初めて気づいたがこうやらないとグリッドがぐちゃぐちゃになる=vtkのグリッド認識がこの順なので。

2010-01-23

3x3 eigenvalue problem: solving cubic function for obtaining principle axes

結局使わなかったけどメモ


背景

http://jikosoft.com/cae/engineering/strmat05.html
みたいな感じ。
ただ胴体が曲がってる場合(とか)を想定して、
応力じゃなくて慣性モーメントの主軸を求めようかなと。

やることは全く同じのようだ。



結論メモ

A をもとの行列とする。
R^{-1} A R = A' とすると、
対角行列 A' の対角項 = 固有値 = 主軸座標系での値
R が固有ベクトル(固有列行列)をまとめた行列 = 座標変換行列(でいいはず)、なので R^{-1} は R^T でOK. まぁ確かめは AR = RA' になってるかでいいので R^{-1} は要らないけど。

んでグリッドを座標変換するときは、グリッド位置を p とすると、
p' = R p
でいいはず。その前に止めたからやってないけど。


というのは、そもそも非対角成分の値が二桁も小さかったので、
A' が A とほとんど変わらなかった。

無視して全然OKっていうか無視すべき。
A っていうか慣性モーメントなので J という記号使ってるけど。



手法

固有値問題なんだけど、一般の n x n のメイトリックスでなくて 3 x 3 なので、
「直説法の方が誤差小さいんじゃね(妄想)」という方針だった。
でも本とかウェブには n x n に対応できるような Jacobi method とか QR method とかの反復法しか出てない。
なので手で解こうとした。

要は det( A -λI) = 0 を解きゃいいはず。
3x3 だから三次方程式。

http://en.wikipedia.org/wiki/Cubic_equation
この解の公式そのままでもできるかもだけど、なんかうまく行かなかった(どうも公式の問題じゃなくて三次方程式の展開をミスってたっぽいorz)

ともあれ、同じこと考えてやってる人がいた(Cで)
http://d.hatena.ne.jp/malibu-bulldog/20081009/1223516459

  http://ja.wikipedia.org/wiki/%E4%B8%89%E6%AC%A1%E6%96%B9%E7%A8%8B%E5%BC%8F
  http://www004.upp.so-net.ne.jp/s_honma/urawaza/eigenvector.htm

これを Fortran に書き直したらうまくいったっぽいけど、
要素 (1,2) と (2,1) が完全に 0 だとおかしくなる(ちゃんと読んでないけどゼロ割っぽい?
1.0d-5 とか微小な非ゼロいれとけばいける。
あと倍精度より4倍精度の方がよかった。


related articles

http://www22.atwiki.jp/linearalgebra/
あと『プログラミングのための線形代数』が完全にこの問題のためにあるんじゃね的なシンクロ率だった。いまさらだけど。しかも使わないけど。

感想
「n が大きい場合には大変なので~」ってよくあるけど「n > 2 の場合には」じゃね?

2009-12-07

A先輩との対話メモ

  • GCL: geometric conservation law だいじょぶ?(変形時。歪ませてるだけならいいが)
  • turbulence をどこまで考慮するか
    1. まずは NACA0012 とかで how much high Re までいけるかを
    2. 実は LE 丸い方が小さな vortices がいっぱい出て turbulence の影響がでかい?
    3. 逆に insect wings/bird feathers のようなテキトーな LE のがいいのやも
  • wing swap は B&B の方がいいかも。飛行状態は、trimmed hover? どうしよう?→ ふってやれば?(あ!
  • htop
  • journals を教えていただいた
  • 先行研究の引用時、批判する必要はない。未知の領域が広大に広がっており、そのうち彼らは一部を、私もまた別の一部をやりましたよ、という感じで
  • Re の u_ref をどこにとるか。Fwd Flt のとき。St = 1 が境目?
  • stiffness と流体の密度との比を保ってない FSI はアウツ
  • hb の SPA, BA は・・・ビデオが正しい?

2009-10-21

Xeon 3500系 (Nahalem-WS) と 5500系 (Nahalem-EP) の違いがよくわからん

Intel Xeonについてあれこれ
http://blog.motdotnet.com/?eid=782763
んじゃ3500系と5500系では何が違うのと言うと、3500系はi7のbloomfieldベースでXeon向けにNehalem-WSとして開発され、5500系はXeonとしてNehalem-EPとして開発されています。
モデルの番号にWだのXだのEだのLだのとアルファベットが付加されていますが、W=ウルトラハイエンド X=ハイエンド E=メインストリーム L=低電力モデルと言う意味らしいのです。
で、3500系はと言うと全てWです。何がWかと言うと多分消費電力なんだと思います。5500系のWと消費電力が同じです。
多分i7のコアベースで電圧をかけてクロックをブンブンと振り回しているのが3500系なのかな?と
なるほど。じゃあ安い 3500 系にしよう。
電圧かけても大丈夫ってことで売ってるんだろうから。

----
ていうか、基本的に 3500系は dual processor 非対応っぽい。
Intel の公式サイト参照