maya mental rayでメモリに制限かかってレンダリングできなかったら

走らなかったら次の3ステップで走るかも．
① 左下のMEL欄で
select mentalrayGlobals;
と叩く

② 追加のアトリビュートのタブをクリックし，開いたタブの
自動メモリのチェック
を外す

③ 自動メモリのチェックのすぐ下の
メモリゾーン
を高くする

ちなみにメモリ制限でレンダリングできないときの警告は，
(Mayatinr.script): Computed memory limit of -XXXMB is below threshold, setting to XXX MB insread

XXXはたぶんものによって違う．知らんけど．

レンダリングメモ：指定した物体のフラグを一括で変えたい

maya 2011でフラグを一括で変えたいときは，アトリビュートスプレッドシートを使うと楽．
アトリビュートスプレッドシートの表示は，メニュー→ウィンドウ→一般エディタ→アトリビュートスプレッドシート．
変えたいフラグを指定して，1を入力すればON，0ならOFFになる．

知らんがmayaのバージョン違いでもできるかも．

レンダリングメモ：サブサーフェイススキャッタリングが変

maya mental ray 等でレンダリングするとき，光源を天空にして，SSSを設定すると，スクリーンで合成（だったかな？名前忘れた)フラグが立っていると，想像した質感と違うことがある．フラグ外すと想像通りかもしれない．

レンダリングメモ：透過物体越しに物体が見えない

maya　 mental ray等のレンダリングで，ある物体Aとカメラの間に透過物体を挟むと，Aが見えないときがある．
Aに屈折に現れる，というようなフラグが立っていない可能性あり．
立てると見えるようになる．
もちろん屈折後の絵が出る．

反射も同様．

レンダリングメモ：白っぽいとき

maya mental ray等でのレンダリングメモ．
レンダリング結果が自分の予想よりも白っぽいとき，反射のハイライトのみ合成すると予想通りになるかもしれない．

シコシコ活動中

ぶへぇ（吐血）大量のモデルのレンダリングが遅すぎる．
いちいち修正に時間がかかるため，ほぼ徹夜で絶賛活動中．
やっぱり，劣化がないLoDもしくは高速で大量のデータを扱う技術や所望の質感を一発設定できる技術が必須だと感じる．慣れろとか言ってる間はずっとしんどい作業を強いられるわ．

つくば旅

さて，くだらない文句言ってないで，某つくばの旅でもまとめてみます．

二日間いましたが，とにかく寒かったです．
パーカーとジーンズ，マフラーで過ごしましたが，震えるくらい寒かったです．

つくばの道路はこんな感じでした．

車がすごく多い印象．
また，運転が荒らそうに見えました．
めちゃくちゃ寒いのにオープンカーを運転しているおっさんとかを見ました．

この次の写真あたりから筑波大学の構内．

ヘリポート？

という感じです．
大学内には学生が多くいました．

今回は某研究会に参加してきました．
CGとCVの合同の研究会でした．
合同だったおかげで，ふたつの研究会の文化の違いを感じることができました．
CG関係で見たことある先生の多くはスーツを着ていませんでした．発表中に発表者にいきなり質問したりしませんでした．

一方，CV関係の先生の多くがスーツを着ていました．発表中にいきなり質問する人もいました．

それ以外は，さほど違いは感じませんでした．

発表に関してですが，多くの発表者が時間内に発表を終えていないようでした．

時間は守ろうぜ，と念じながらの参加になりました．

懇親会では，いつも通りおとなしくしておりましたが，数名の方とお話ししました．
はじめてお話しする人もいました．
顔見知りが増えてよかったと思います．

このような懇親会は，国内トップ層と直に，研究内容ではなく，人間を見れる珍しいチャンスなので面白いです．

研究会には，私が所属している研究所の人も数名参加しており，少し話をしました．

ドライな冗談を聞きましたが，それくらい仲が良く，風通しがいいということでしょう．

こんなもんでしょうか．

とにかく発表は面白いので，周りの学生も活発に発表するようになればいいと思います．

それとmasafumiさんはやっぱり有名人で，いろいろ聞かれたので，聞いてきたやつらに嘘8億吹き込んでおきました．これで超絶いい人だと勘違いした人が増えると思います．

驚懼疑惑が起こるなんて当たり前だから常に正常な判断できるようにしておけという理解

驚懼疑惑を克服するということは，驚いたり恐怖したり疑心を抱いたり惑ったりということが心に起こると，正常な判断ができなくなるから，心を鍛えて恐懼疑惑を起こらないようにしておく，ということなのだけれど，恐懼疑惑しない人なんてあまりいない．
どういうわけか，この世では想像もつかないことが簡単に起きる．
いくら準備しておいても，恐懼疑惑するのだから，しょうがない．
だから，恐懼疑惑を避けるよりも，恐懼疑惑したときでも正常な判断できるように，正常な判断の基準を自分の中で持っておいて，いかなるときにも，その基準を満たすように活動しておくことが大切だと理解しています．
今のところ．

32bitグレースケールはfloatで[0,1]の値

グレースケール形式にハマりました．

一般的なグレースケール画像は8bitグレースケールですが，より分解能が高い形式もあり，16, 32, 64bitグレースケールがあります．（たぶん今後はもっと分解能上がる）

8bitグレースケールは，画素あたり[0,255]整数で表現されるので，一般的に整数で表現します．
プログラミングでの書き出しもこの範囲の整数で書きだします．
しかし，同じノリで16,32,64も整数で記述するものと思っていたのですが，違ったようです．
32,64は整数で扱うソフトウェアがかなり少なく，浮動小数点形式が一般的のようです．
このような情報をどこかに明記してくれていればありがたいのですが，見つからなかったので，photoshopでグレースケール画像を生成して，tiff形式で保存して，tiff画像のタグを見ることで，知ることができました．

はじめに，16bitと32bitグレースケールをtiff形式で作ります．
次に，TIFFの分析ソフトを使ってタグを見ます．私はTiffAnalyzerを使いました．
http://www.vector.co.jp/soft/mac/art/se136821.html

タグを見た結果，16bitグレースケールのときは，タグなしだったのに，32bitグレースケールでは，sampleformatが3になっていました．
sampleformatが3の意味は，浮動小数点形式で画素を表現するということです．
下のページは，TIFFのタグ「sampleformat」について書かれたページです．
http://www.awaresystems.be/imaging/tiff/tifftags/sampleformat.html


32bitで浮動小数点といえば，floatなので，保存はデータをfloat形に変換して出力することがわかります．
また，実験的に生成したグラデーションのデータを観察すると[0,1]表現だったので，32bitグレースケールはfloatの[0,1]で表現することがわかります．

この条件を満たすように32bitグレースケールを生成して，いろんなソフトで開けるかを実験します．

作った32bitグレースケール画像は，ウィンドウズ7の標準のビューワーやphotoshopなどなど複数のソフトウェアで開くことができました．

一方，実験の中で作った，整数型のグレースケールは32bit指定しても，多くのソフトウェアで開けませんでした．

ということで，既存のソフトウェアで32bitグレースケール画像を開けるような画像を生成したときは，floatに変換して書きだす必要があります．

桁違い

例にもよって，このブログでは研究活動の内容を書かないようにしているのですが，ちょろっとメモ．

今，M先生と論文原稿テニスをしているのですが，その中で桁違いの実力の差を感じました．
先生はものすごく編集が速いのです．
今回のテニスでは，編集にMS Wordを使っています．
いつからか知りませんがWordでは，リボンのファイルタブの情報という項目で編集時間が見れるようになっています．
その機能を使って，意図的ではないですが，先生の編集時間を知ってしまいました．
ものすごく短いのです．
6から7ページの原稿を概ね1時間で編集しています．

私は一文の書き方が決まらずに1時間過ぎてしまうことがあるというのに，先生はそうこうしている間に一編を適切に編集してしまうわけです．
自分の場合は一時間かけても不適切なことがあるのに，と思うと，もう桁違いだと諦めるほかありません．
やはり先生と呼ばれる先生は，あーでなくっちゃ駄目だと思いました．

GI含めたレンダリングは時間がかかる

毎度毎度レンダリングするたびに思いますが，レンダリングは時間がかかる．
600万ポリゴン程度でもグローバルイルミネーションつけるだけで1時間かかる．
待ってる間はPCで遊んでろという話かもしれませんが，CPUの使用率100%の中ではさすがにPCを使う気にはなりません．

あー，毎度思う．

maya 2011 plugin作成メモ

VS2010で開発してますが，VS2010ではmayaのプラグイン開発用ウィザードが走ってくれないため，手動でプロジェクトのプロパティを編集するほかありません．
ネット上でどこもまともなこと書いてないので，実験しながら適当にプロパティ変更しました．
そして走ったので，変更点のメモです．

プロジェクトをwin32コンソールアプリケーション（DLL）として作成します（空のプロジェクトの必要なし）．
プロジェクトを作成後，cppファイルがなければ作成します．
その後，プロパティページの「構成プロパティ」内の以下4点を変更．
1. 「全般」→「ターゲットの拡張子」をdllからmllに変更．
2. 「C/C++」→「全般」→「インクルードディレクトリ」と「リンカー」→「全般」→「追加のライブラリディレクトリ」はmayaのincludeとlibフォルダを指す．
3. 「リンカー」→「コマンドライン」→「追加オプション」欄に，
/export:initializePlugin /export:uninitializePlugin
を記述．
4. 「リンカー」→「入力」→「追加の依存ファイル」欄に
Foundation.lib
OpenMaya.lib
OpenMayaUI.lib
OpenMayaAnim.lib
OpenMayaFX.lib
OpenMayaRender.lib
Image.lib
を記述

以上4点の変更でmaya2011に読んでも走るpluginができました．
変更点以外はデフォルトのままです．
VS2008のときはもっと変更点があるような文書もありましたがそんなの知りません．
コードはapiを呼んで，initializePluginとuninitializePluginを書いておけば問題ありませんでした．

グローバルイルミネーションのチェック

maya使っててmental rayでレンダリングしようとしているのに，絵が出ないときがあって，調べてもわからなかったので，レンダリングの設定をさわっていたら出るようになりました．
「ウィンドウ」→「レンダリングエディタ」→「レンダー設定」で使用するレンダラをmental rayに設定して，「機能」タブ→「二次エフェクト」→グローバルイルミネーションをチェックです，絵が出ました．理由がよくわからない素人ですが，結果的に出たということでメモです．

内部の例外を見えなくする件

しばらくCGALを触っているとVCのデバックにて， CGAL::Uncertain_conversion_exception という例外のメッセージを目にします．このメッセージはCGALのFAQによると無害のメッセージとのことです．ということは見えなくても問題ではなく，ただただ鬱陶しいだけです．同様に思っている人が表示させない方法を紹介してくれてるのでメモです．

http://www.helixoft.com/blog/archives/24

計算誤差のせいで現れるセルの排除

3次元ドロネー図を計算させるとき，計算誤差のせいでありもしないドロネーセルが呼ばれることがありますが，それは体積を計算して，極々小さいときに排除することで解決できます．
原因がわからなかっただけに二日も費やしてしまいました．

某舘山寺温泉に行ってきました

いってまいりました，某ワークショップに，某舘山寺温泉に．
受講の都合で，後半からの参加になりましたが想像以上に楽しかったです．
他大学の学生との交流もほどほどにできました．
それに温泉にもつかってきました．その場では，家で入ってる風呂と大差ない気がしましたが，疲れも無くなりましたし，実際は効果あったような気がします．
夕飯のコースもおいしかったです．

発表は最後のセッションの最後でして，皆疲れている様子でしたが，お構いなしでよくしゃべれました．おしゃべりの内容はテクスチャです．詳細は教えません．会場からの質問も見た中では最も多かったですし，興味を持っていただいたようでした．質問の内容は，発表した手法をより一般的に使えるようにするための議論（高さのみならず色や材質，メッシュに依存しないようにするためには等々の議論）や入力情報を足して操作性を上げるための議論がメインでした．
実はこれらの質問は，事前にM研究室で行った発表練習のときに修士2年の学生から出た質問と似たような内容でした．つまり国内のCG界のトップ層と同等の質問ができるくらいに，ここの修士が育っていることがわかりました．

はい，これで発表が終わり，会場を後にし，某M先生と一緒に浜松に行って鰻を食べてきました．行った店は，某清水屋です．某早稲田の博士後期の学生(自称M先生ファン)が先生に教えてくれたそうです．
で，しばらく待っていたのですが，出てきたうな重は最高においしかったです．これまで食べてきた鰻とは別格においしかったような気がしました．最近良いものを口にしていないので，特にそのように思ったのかもしれませんが，とにかくおいしかったです．サブのくせにビビビッと来たのが肝吸いでして，一口すすっただけで，ニヤリです．ちなみにおごってもらいました．また頑張ります．
あんなにおいしいものを食べて思ったのが，あんなのばっかり食ってたら人間ダメになる．

そういえば，食べ終わったあと，先生が手のにおいを嗅いで「鰻くさい」と言ってましたがその時は理解できませんでした．電車乗って自分が鰻くさいことに気付きました．

そんなこんなで，時間もないので，さっさと帰ってきました．

cmake→slnファイル→ビルド

インストールメモ．
CMAKEしたあと，作られたslnソリューションファイルを開いて，ビルドしなければ，ライブラリが生成されません．ビルド時は，ソリューション構成が4つ選択できるので，それぞれのビルドします．これでライブラリが生成され，参照できるようになります．

Polyhedron_incremental_builder_3でのインプットエラー

Polyhedron_incremental_builder_3でのインプットに関して，エラーが出たのでメモ．
facetを指定するとき，すべて同じ方向でインデックスを用意する必要があるようです．
たとえば下の例では，7654をコメントにして，//でコメントされた4567でfacetを指定すると，3267と反対になるため，エラーがでます．

 B.begin_facet();
        B.add_vertex_to_facet(3);
        B.add_vertex_to_facet(2);
        B.add_vertex_to_facet(6);
        B.add_vertex_to_facet(7);
 B.end_facet();
 B.begin_facet();
        B.add_vertex_to_facet(7);
        B.add_vertex_to_facet(6);
        B.add_vertex_to_facet(5);
        B.add_vertex_to_facet(4);
        //B.add_vertex_to_facet(4);
        //B.add_vertex_to_facet(5);
        //B.add_vertex_to_facet(6);
        //B.add_vertex_to_facet(7);
 B.end_facet();

ちなみにエラーはpolyhedron_incremental_builder_3.hの440行目でm_errorがtrueになることで発生します．
同じ方向の面を検出するチェッカーがこの行の前に入っていて，m_errorがtrueになることで次からの入力を拒絶しています．

Polyhedron_incremental_builder_3を使ってみる

これまで私も任意の多角形メッシュをCGALで扱ったことなかったので，知らなかったのですが，CGALでは Polyhedron_incremental_builder_3 を使って，閉じた多角形メッシュを作ることができます．
たとえば，examplesのソースをいじって六面体を作ると以下のようになります．

#include < CGAL/Simple_cartesian.h >
#include < CGAL/Polyhedron_incremental_builder_3.h >
#include < CGAL/Polyhedron_3.h >

template < class HDS >
class Build_mesh : public CGAL::Modifier_base< HDS > {
public:
    void operator()( HDS& hds) {
        // Postcondition: `hds' is a valid polyhedral surface.
        CGAL::Polyhedron_incremental_builder_3< HDS > B( hds, true);
        B.begin_surface( 8, 6 );
        typedef typename HDS::Vertex   Vertex;
        typedef typename Vertex::Point Point;
        B.add_vertex( Point( 0, 0, 0));
        B.add_vertex( Point( 0, 1, 0));
        B.add_vertex( Point( 1, 1, 0));
        B.add_vertex( Point( 1, 0, 0));
        B.add_vertex( Point( 0, 0, 1));
        B.add_vertex( Point( 0, 1, 1));
        B.add_vertex( Point( 1, 1, 1));
        B.add_vertex( Point( 1, 0, 1));
  B.begin_facet();
        B.add_vertex_to_facet(0);
        B.add_vertex_to_facet(1);
        B.add_vertex_to_facet(2);
        B.add_vertex_to_facet(3);
  B.end_facet();
  B.begin_facet();
        B.add_vertex_to_facet(3);
        B.add_vertex_to_facet(2);
        B.add_vertex_to_facet(6);
        B.add_vertex_to_facet(7);
  B.end_facet();
  B.begin_facet();
        B.add_vertex_to_facet(7);
        B.add_vertex_to_facet(6);
        B.add_vertex_to_facet(5);
        B.add_vertex_to_facet(4);
  B.end_facet();
  B.begin_facet();
        B.add_vertex_to_facet(4);
        B.add_vertex_to_facet(5);
        B.add_vertex_to_facet(1);
        B.add_vertex_to_facet(0);
  B.end_facet();
  B.begin_facet();
        B.add_vertex_to_facet(0);
        B.add_vertex_to_facet(3);
        B.add_vertex_to_facet(7);
        B.add_vertex_to_facet(4);
  B.end_facet();
  B.begin_facet();
        B.add_vertex_to_facet(1);
        B.add_vertex_to_facet(5);
        B.add_vertex_to_facet(6);
        B.add_vertex_to_facet(2);
  B.end_facet();
  B.end_surface();
    }

};

add_vertexで頂点を追加して，begin_facetとend_facetの間に面を定義します．
面の定義はadd_vertex_to_facetでインデックスを追加します．
これらをbegin_surfaceとend_surfaceではさめば，定義完了です．

呼び出して使うときは，まず以下を記述します．

typedef CGAL::Simple_cartesian < double >  Kernel;
typedef CGAL::Polyhedron_3 < Kernel >      Polyhedron;
typedef Polyhedron::HalfedgeDS          HalfedgeDS;

呼び出し側は，以下のようにPolyhedronを用意して，delegateするようにします．

Polyhedron P;
Build_mesh< HalfedgeDS > mesh;
P.delegate( mesh );

この例から任意の形状用に作りかえれば，細分割に任意の形状を突っ込むことができるようになります．
exampleがこのような構成になっているので，きっとこれでいいのだと思います．

subdivisionの例

CGALには頻繁に使われる細分割手法が実装されていて，かつそこそこ高速に計算してくれるので非常に使いやすいです．ただマニフォールドのみという制約があった(と記憶している)ので，そこは使用者が実装しなければならない(はず)です．
使い方はPolyhedronと分割数を関数に渡すだけです．

たとえば以下のコードで3回のDoo-Sabinです．
まず，関数を使えるようにヘッダーをインクルードします．

#include < CGAL/Subdivision_method_3.h >

呼び出し側で以下のコードです．

Polyhedron P;
Build_mesh< HalfedgeDS > mesh;
P.delegate( mesh);

CGAL::Subdivision_method_3::DooSabin_subdivision(P, 3);

ここでは，Build_meshで六面体のPolyhedronを生成して，delegate関数でPにデータを渡しています．
0-3回の例が以下です．






以上．
簡単でしょ．

VSでCGALを読み込んでコンパイルが通るまで

CGALを使ったプログラムのビルドが通らない原因はCGALのインストールがうまくいっていないことだと思います．私も結構苦労しました．2時間くらい調べたりしてようやくわかった感じです．ですので，ここでは，CGALが本当にインストールできたかをテストする方法を記述します．
下で示す手続きでビルド（コンパイル）を実行して通れば，CGALのインストールができていると言えます．逆にビルドが通らなかったら，CGALのインストールができていないと言えます．
ここではVisual studio 2010で実験しますので，ほかの開発環境ではどうなのか知りません．でも，だいたい似たようなものだと予想します．

さて，テスト用の手続きですが，ざっくりいうと次の5ステップです．
1 プロジェクトの作成
2 ソースのコピぺ
3 インクルードディレクトリーの設定
4 ライブラリディレクトリーの設定
5 ビルド
以上です．

1 プロジェクトの作成
VSではまず，プロジェクトを作成しなければビルドも何もできません．
ということで，作ります．
ここではテストなので，win32コンソールアプリケーションを選びます．
で，何も考えずに，完了です．





2 ソースのコピぺ
次にソースのコピペですが，こんなものは何でもいいので，CGALのサンプルからコピーして，貼り付けます．
ここでは，私がこれから着手するサブディビジョンのテスト用に，3D Surface Subdivision Methodsのサンプルをとってきます．



赤い波下線があり，コンパイルエラーが起こることがばっちりわかります．

3 インクルードディレクトリーの設定
次にプロジェクトのインクルードディレクトリーを設定します．
手順は，プロジェクトのプロパティを開いて，「C/C++」→「全般」→「追加のインクルードディレクトリ」にCGALのインクルードディレクトリとboostのインクルードディレクトリを追加します．



ここでのディレクトリは環境に依存するので，適宜調整してください．

4 ライブラリディレクトリーの設定
インクルードディレクトリー同様，リンカーの設定として，ライブラリを設定します．
手順は，プロジェクトのプロパティを開いて，「リンカー」→「全般」→「追加のライブラリディレクトリ」にCGALのライブラリとboostのライブラリを追加します．



これもインクルードディレクトリの設定と同様に，適宜調整してください．
この例では，デバック版を設定していますが，リリース版でビルドする場合は，リリース版を設定します．

5 ビルド
以上で準備は整いました．
ビルドします．



赤い波下線も消え，ログからビルドが正常終了したことがわかります．
これが確認できれば，CGALのインストールは成功しています．

もし，以上の手続きで，エラーが出るようでしたら，失敗しています．
経験上，失敗の例として，CMAKEを実行していない，CMAKE後にコンパイルしていない，CMAKEがミスってる，インストールしたboostのバージョンが悪い，そもそも必要な他のライブラリをインストールしていない等々です．


ユーザー名がKSakuraiってなっていますが，そのあたりは気にせず，無意味に特定しようとか思わないでください．

CGAL3.7リリースされてた

CGAL3.7がリリースされています．VS2010でも，インストールがしやすくなっています．

cgalをw7, vs2010にインストール

windows 7, visual studio 2010でのインストールの確認ができました．
正直知らなきゃ入らないし，もうしたくない気持ちです．
おおむね公式サイトのFAGに書いてある通りすれば，OKですが，ひとつ罠がありました．
http://www.cgal.org/FAQ.html#vc10

上記のサイトの通りすると，C:\Program Files (x86)\CGAL-3.6.1\auxiliary\gmp\libのフォルダ内に，libgmp-10.lib,libgmp-10.dll,libmpfr-4.lib, libmpfr-4.dllがある状態になります．
このlibgmp-10.libとlibmpfr-4.libの名前をgmp-vc100-mt.libとmpfr-vc100-mt.libに変更しなければ，コンパイルが通りませんでした．

私はコンピューターなんか全然詳しくないので，何が何でこういう変更が必要なのか，何が何でCGALはこの名前を指定するようにしているのか知りませんが，結構げんなりする罠でした．

あー

cgalがvisual studio 2010でcmakeできないんですけど

できるっぽいです．
boostがひっかかってるだけか．

cgalでのオブジェクトへのアクセス順

cgal全般に言えるのかもしれないけれど，オブジェクトをinsertした順にイタレータでアクセスできるわけではないみたいですね．
順にアクセスしたいときは，locate関数で位置を指定するのが定石なのか．

あかんわ

某blenderの使い方で細かいところがわからなすぎるので，プロ仕様のソフトを使うことにしました．
でも，ライセンス管理している人が研究室に来てないので，手を動かせないので，つまらんのです．
ちょいイライラします．投稿まで時間が短いからでしょうね．
でも，某M先生みたいに気長にならなければ，良い仕事できないんでしょうね．勉強になります．はい．

Faceのinfo, index

近隣のフェイスを参照するとき，これらのインデックスが着いていたほうが扱いやすいです．しかし，CGALでは標準でidを持ったfaceの形式がないため，別の情報で代替する必要があります．
ここで，適当なのはinfoで，どんな形式でも追加できる情報です．
今回は，Delaunay図でのインデックスを得る実験をします．
infoを持つfaceを使うために，以下のようにインクルードを指定します．

#include < CGAL/Exact_predicates_inexact_constructions_kernel.h >
#include < CGAL/Delaunay_triangulation_2.h >
#include < CGAL/Triangulation_vertex_base_with_id_2.h >
#include < CGAL/Triangulation_face_base_with_info_2.h >

最後のTriangulation_face_base_with_info_2.hがfaceにinfoを付加した形式です．
3番目はid付の頂点情報です．faceのインデックスの場合はこれである必要はありません．
typedefは以下のように設定します．

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Triangulation_vertex_base_with_id_2< K > Tvb;
typedef CGAL::Triangulation_face_base_with_info_2< int, K > Tfb;
typedef CGAL::Triangulation_data_structure_2< Tvb,Tfb > Tds;
typedef CGAL::Delaunay_triangulation_2< K, Tds > Triangulation;

Triangulation_face_base_with_info_2< int, K > Tfbの最初の引数でinfoの形式を指定します．この場合は，intを指定します．

 int cnt;
 cnt = 0;
 for(Finite_vertices_iterator fvi = t.finite_vertices_begin (); 
fvi!=t.finite_vertices_end(); fvi++)
 {
  fvi->point().x();
  fvi->point().y();
  fvi->id() = cnt;
  cnt++;
 }

 cnt = 0;
 for(Finite_faces_iterator ffi = t.finite_faces_begin (); 
ffi!=t.finite_faces_end(); ffi++)
 {
  int f[3];
  f[0] = ffi->vertex(0)->id();//vertexから引いてくる
  f[1] = ffi->vertex(1)->id();
  f[2] = ffi->vertex(2)->id();
  ffi->info() = cnt;
  cnt++;
 }

前半は，vertexのidで，後半がこの話題の中心であるfaceのインデックスです．idの指定同様
ffi->info()に直に通し番号を入れていきます．
これで，faceのインデックスを参照できます．

VertexにID, index

CGALでデータを扱うときは，イタレーターと通すので，頂点(以下vertex)のインデックスなんて気にしません．ただ，そのデータをファイルに出力するときは，概ねvertexのインデックスを扱う必要が出てきます．たとえば，wavefrontのOBJ形式はフェイスをインデックスで指定します．
このようにvertexのインデックスを得るためには，CGALでは，idを用いればよいっぽいです．まだ，ほかの方法があるのかもしれませんが，よくわからないので，私はこの方法でインデックスを得ることにします．
今回は，Delaunay図でのインデックスを得る実験をします．
vertexにidを持った形式を使うために，インクルードは以下のように指定します．

#include < CGAL/Exact_predicates_inexact_constructions_kernel.h >
#include < CGAL/Delaunay_triangulation_2.h >
#include < CGAL/Triangulation_vertex_base_with_id_2.h >

この最後のインクルードであるTriangulation_vertex_base_with_id_2.hがid()関数を持ったvertexの形式を扱えるようにするヘッダーファイルです．
typedefは以下のように設定します．

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Triangulation_vertex_base_with_id_2< K > Tvb;
typedef CGAL::Triangulation_face_base_2< int, K > Tfb;
typedef CGAL::Triangulation_data_structure_2< Tvb,Tfb > Tds;
typedef CGAL::Delaunay_triangulation_2< K, Tds > Triangulation;

Triangulation_data_structure_2でTriangulation_vertex_base_with_id_2形式を指定してやります．
以上で，vertexにてid()が使えるようになります．
さて，次にこのid()を活用してインデックスを指定します．

イタレーターでvertexにアクセスして，そのid()に整数を入れます．0からループごとにインクリメントして，id()にいれます．
フェイスでは，vertexのid()を参照して，格納します．
以下が簡易的なサンプルです．

 int cnt;
 cnt = 0;
 for( Finite_vertices_iterator fvi = t.finite_vertices_begin (); 
fvi!=t.finite_vertices_end(); fvi++)
 {
  fvi->point().x();
  fvi->point().y();
  fvi->id() = cnt;
  cnt++;
 }

 cnt = 0;
 for( Finite_faces_iterator ffi = t.finite_faces_begin ();
 ffi!=t.finite_faces_end(); ffi++)
 {
  int f[3];
  f[0] = ffi->vertex(0)->id();//vertexから引いてくる
  f[1] = ffi->vertex(1)->id();
  f[2] = ffi->vertex(2)->id();
 }

以上で，vertexの通し番号のインデックスを得ることができます．

2次元三角形のトラバース

CGALで，2次元三角形による領域分割が必要になったので，2次元ドロネー図の三角形をトラバースできるか実験してみました．
まず，インクルードファイルとtypedefの宣言します．

#include < stdlib.h >
#include < math.h >

#include < CGAL/Exact_predicates_inexact_constructions_kernel.h >
#include < CGAL/Delaunay_triangulation_2.h >

#include < fstream >

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;

typedef CGAL::Delaunay_triangulation_2< K > Triangulation;
typedef Triangulation::Finite_faces_iterator Finite_faces_iterator;
typedef Triangulation::Edge_iterator Edge_iterator;
typedef Triangulation::Point Point;

以上で二次元ドロネー図を構成する準備ができました．
次にコードですが，簡単にいうと，ドロネー図に点を追加していって，できたグラフをトラバースするという手続きです．
点を追加する時点で，図を構成していますので，構成の関数を呼ぶことはないです．

  Triangulation t;
  for(int i=0;i < N;i++){
   t.insert(Point((double)x[i][0], (double)x[i][1]));
  }

  int faceN = t.number_of_faces();
  double *face = new double[faceN*3*2];//x 3 vertices x 2D

  int i = 0;
  for(Finite_faces_iterator ffi = t.finite_faces_begin (); ffi!=t.finite_faces_end(); ffi++)
  {
   face[i*3*2+0] = ffi->vertex(0)->point().x();
   face[i*3*2+1] = ffi->vertex(0)->point().y();
   face[i*3*2+2] = ffi->vertex(1)->point().x();
   face[i*3*2+3] = ffi->vertex(1)->point().y();
   face[i*3*2+4] = ffi->vertex(2)->point().x();
   face[i*3*2+5] = ffi->vertex(2)->point().y();
   i++;
  }

xはランダムに配置した点です．
トラバースは三角形ごとに見ているので，faceのイタレーターで参照しています．

CGALでのメッシュ生成

以下のサイトのGenerating Meshes with CGALに，関数からメッシュを作ってメッシュにアクセスする方法が書いてある．
http://schmid.dk/wiki/index.php/CGAL

この方法で，メッシュ生成すると，ちゃんとメッシュが作れた．
ユーザーの気持ちを知っているすばらしいサイトだと思う．

Facet_iteratorでcomplex上のfaceにアクセスし，アクセスしたfaceのfirst(四面体)からsecond以外の頂点にアクセスする．そのsecond以外の頂点は，complex上のfaceなので，ちゃんと頂点が取れてくる．ということです．
知らないとできないです．こんなの．

C2t3でメッシュ生成後，以下のコード(上のサイトからの転載)で，メッシュにアクセスできます．

    // for each facet (triangle), print the vertices
    for(C2t3::Facet_iterator fi = c2t3.facets_begin();
            fi != c2t3.facets_end(); fi++) {

       C2t3::Cell_handle cell    = fi->first;
       int opposite_vertex_index = fi->second;
       for(int i = 0; i < 4; i++)
           if(i != opposite_vertex_index) {
               std::cout << "(" << cell->vertex(i)->point().x()
                         << "," << cell->vertex(i)->point().y()
                         << "," << cell->vertex(i)->point().z() << ")"
                         << std::endl;
           }
    }

ありがとうございました．

CGAL理解できず

さっぱりわからなくなってきた．
handleとiteratorが理解できればほしいデータにアクセスできるのはわかるのだけれど，それらをどうやって使えばいいのかさっぱりわからない．そのせいで，ほしい三角形の頂点にアクセスできない．リファレンスにいろいろ載っているのだけれど，意味が理解できない．無駄に時間を過ごしているような気がしてきた．

陰関数でのメッシュ生成は多少時間がかかる

陰関数からメッシュ生成することが素晴らしいと思って開発していたのだけれど，計算時間がかかりすぎることを忘れていた．
ひとつの形状を作るのには，優れているかもしれないが，めちゃくちゃ多くのものを作るには，面白くない．
それに，今回の目的に対して陰関数である必要性が，弱い．
他の方法で，簡単にメッシュ作れるならばそっちを適用した方がかっこいいので，そっちに乗り換える．

Oracleによって

CGALのリファレンスにパラメータはOracleによって設定されるというような記述があったので面白い表現だということで，メモ

ドロネー図の点の参照

CGALは，creationの時点で，グラフを作ります．
ですので，宣言して，insertで点を入れていくだけでグラフが形成されます．
問題はどうやってそのグラフを参照するかだったんですが，cellかvertexのiterator経由でハンドルつかんで，点を順に参照していけばよいです．
ドロネー図をFinite_cells_iterator経由で，四面体の4頂点を自作四面体クラス(Tetrahedron)に格納していく例を挙げます．
自分用メモみたいな感じですが，以下コードです．

 Delaunay dt;
 for(int i=0; i < _siteN; i++){
  dt.insert(Point(_site[i][0], _site[i][1], _site[i][2]));
 }
 //std::cout << dt.number_of_finite_cells() << "\n";
 tetrahedraN = dt.number_of_finite_cells();
 tetrahedra = new Tetrahedron*[tetrahedraN];
 for(int i=0; i < tetrahedraN; i++){
  tetrahedra[i] = new Tetrahedron();
 }

 int cnt = 0;
 for(Delaunay::Finite_cells_iterator fci = dt.finite_cells_begin(); fci != dt.finite_cells_end(); fci++){
  for(int j=0;j < 4;j++){
   Vertex_handle v = fci->vertex(j);
   tetrahedra[cnt]->set(j, 0, v->point().x());
   tetrahedra[cnt]->set(j, 1, v->point().y());
   tetrahedra[cnt]->set(j, 2, v->point().z());
  }
  cnt++;
 }

siteは母点です．tetrahedraのset関数は，インデックスと次元，位置を格納します．

IntelliSenceでは見えないフィールドがたくさんある

ようやくCGALのリファレンスの意味がわかってきました．
CGALのライブラリはIntelliSenseがあると，見えないフィールドがあるので，使いたい関数が使えないと思っていました．でも実際は，使えるわけで，それを使えば，リファレンスに書かれた機能がサクッと使えます．
例はまた後で書くとして，わかったということだけメモっておきます．明日くらいにコード挙げようと思います．

CGALビルド時の設定 intelliSenseのため

CGALをVC用にCMAKEするとslnファイルが出力されますが，このままビルドすると読み込み先でintellSenseが働かないライブラリが出力されます．
プロパティで\cln(共通言語ランタイムのコンパイル)に設定しないとだめなようです．
なぜかはよくわかりません．

それ以外にCGALをシステムパスに登録しておかなければならないのは，前回の通りです．

CGALのインクルードディレクトリとライブラリディレクトリ

CGALが3.6にバージョンアップしてたので，インストールしてみました．

若干構造が変わった気がしますが，誤差の範囲です．
インストール方法は，以前と変わりません．
管理者として実行でインストールする方が楽ってことを気付いたくらいです．

走らせるときのパスは，以下です．
追加のインクルードディレクトリは
C:Program File\CGAL-3.6\include
C:Program File\CGAL-3.6\auxiliary\gmp\include
C:Program File\CGAL-3.6\auxiliary\tausc\include
C:Program File\boost\boost_1_14

追加のライブラリディレクトリは
C:Program File\CGAL-3.6\lib
C:Program File\CGAL-3.6\auxiliary\gmp\lib
C:Program File\CGAL-3.6\auxiliary\tausc\lib
C:Program File\boost\boost_1_14\lib

これで同様に走ったので，問題ないでしょう．

Vistaは管理者として実行か

右クリックから管理者として実行すれば問題なくインストールできた．
CGALもCMAKEも全部管理者として実行です．
とくに，CMAKEの実行は絶対に管理者として実行しなければなりません．
動けば何でもいいです．

思い通りのレンダリングが出ない1例

テクスチャが設定した通りにレンダリングされない．
と1日中悩んでいたが，原因はmipmapで，これが適用される設定になっていた．
そのため，場合によってぬるいテクスチャが適用される．

はじまりました

さっそく研究用のソフトウェアを開発し始めました．
さくっといい結果が出てくれるとうれしいです．

CGALのVoronoi図計算実験

CGALで2次元のボロノイ図を計算させてみました．母点をランダムに配置した後，ボロノイ図を計算し，ボロノイエッジを描画しました．上の図のように描画ができますが，ソースコードはいまいちな気がします．

ボロノイ図の計算にあたって，リファレンスの2D Voronoi Diagram Adaptorを読んで使い方を学ぼうと思ったのですが，思考がついていかず，細部の設計が理解できませんでした．ですので，付属のサンプルのソースコードを追っかけてなんとなく使う作戦にしました．
ボロノイ図を計算させるexamples/Voronoi_diagram_2/vd_2_point_location.cppのコンパイルと実行が確認できたので，このソースを元に簡単に変更していきました．
まずは，インクルード周りですが，元にしたソースのままです．

#include < CGAL/basic.h >

// includes for defining the Voronoi diagram adaptor
#include < CGAL/Exact_predicates_inexact_constructions_kernel.h >
#include < CGAL/Delaunay_triangulation_2.h >
#include < CGAL/Voronoi_diagram_2.h >
#include < CGAL/Delaunay_triangulation_adaptation_traits_2.h >
#include < CGAL/Delaunay_triangulation_adaptation_policies_2.h >

// typedefs for defining the adaptor
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Delaunay_triangulation_2 < K > DT;
typedef CGAL::Delaunay_triangulation_adaptation_traits_2 < DT > AT;
typedef CGAL::Delaunay_triangulation_caching_degeneracy_removal_policy_2 < DT > AP;
typedef CGAL::Voronoi_diagram_2 < DT,AT,AP > VD;

// typedef for the result type of the point location
typedef AT::Site_2 Site_2;
typedef AT::Point_2 Point_2;

typedef VD::Locate_result Locate_result;
typedef VD::Vertex_handle Vertex_handle;
typedef VD::Face_handle Face_handle;
typedef VD::Halfedge_handle Halfedge_handle;
typedef VD::Ccb_halfedge_circulator Ccb_halfedge_circulator;

次にボロノイ図の計算周りですが，ランダムな点を母点とし，母点の位置にある要素をとってきて描画しました．

 double w = 2.5;
 double h = 2.0;

 int n = 100;
 double** sitePos = new double*[n];
 for(int i = 0; i < n; i++){
  sitePos[i] = new double[2];
  sitePos[i][0] = math::myRandom(-0.50)*2.0*w;
  sitePos[i][1] = math::myRandom(-0.50)*2.0*h;
 }

 //2D Voronoi adaptor
 VD vd;
 for(int i = 0; i < n; i++){
  //insertion with insert() func
  vd.insert(Site_2(sitePos[i][0], sitePos[i][1]));
 }

 for(int i = 0; i < n; i++){
  //母点の指定
  Point_2 p(sitePos[i][0], sitePos[i][1]);

  //指定母点上のボロノイ図の要素
  Locate_result lr = vd.locate(p);

  //Face_handleを受け取る
  Face_handle* f = boost::get< Face_handle >(&lr);

  //Face_handleからhalfedgeを受け取る
  Ccb_halfedge_circulator ec_start = (*f)->outer_ccb();

  //Voronoi edgeのスタート地点をecに渡す
  Ccb_halfedge_circulator ec = ec_start;

  do {
   //Voronoi edgeの有無
   if(ec->has_source()){
    //voronoi edgeの端点その1の位置
    double vxs = ec->source()->point().x();
    double vys = ec->source()->point().y();
    if(ec->has_target()){
     //Voronoi edgeの端点その2の位置
     double vxt = ec->target()->point().x();
     double vyt = ec->target()->point().y();
     //描画
     glVertex2d(vxs,vys);
     glVertex2d(vxt,vyt);
    }
   }
   else{
    //無限遠
    std::cout << "inf" << std::endl;
   }
  } while ( ++ec != ec_start );//スタート地点に戻ってくるまでのdo while文
 }

 for(int i = 0; i < n; i++)
  delete[] sitePos[i];
 delete[] sitePos;

ボロノイエッジの描画だけならもっとさらっとしたコードになりそうですが，今はよくわからないことが多いので，これが今の実力では最善です．

インテリセンスとシステムパス

ライブラリを用いた時に，インテリセンス(intellisense)が働かないことが起こる理由のひとつに，そのライブラリがシステムパスに登録されていない可能性がある．

http://msdn.microsoft.com/ja-jp/library/ms235517.aspx

知らなかったから，インテリセンス使えなかった．

CGALインストールメモ

SIGGRAPH ASIA 2009のコースで初めてCGALの存在を知り，使ってみたいと思いつつも，使わないままこの時期になってしまいました．
遅くなったが，とにかく始めてみよう．
CGAL
http://www.cgal.org/

まずはインストールしなければなりません．
しかし，インストールするだけが，結構大変です．

いろいろとハマったので，次からはハマらないようにメモしておきます．


今回の最新バージョンは3.5.1なので，これをダウンロードしました．環境はwindows + VCなので，インストーラーがあります．
ダウンロードしたインストーラーを実行して，Userフォルダに展開します．
私の環境はWindows Vistaですので，Program Fileフォルダ内では，勝手な書き換えは許していません．ですので，のちのcmake時にアクセスの問題がでます．それを避けるためにUserフォルダ内に展開しました．

展開すると，インストールマニュアルがあるので，順に進めていくことにします．

CGALを使うために，必要なライブラリがあったので，以下にメモ．

・cmakeとboostが必要
cmakeは以下のサイトからダウンロード
http://www.cmake.org/

boostはBoostProでwindowsのインストーラーが手に入る(登録制)
http://www.boostpro.com/download
boostは，boost_threadを含めてインストールする．含めないと，cmake時に問題になる．

また，バージョン1.42以上だと定義エラーが出るので，1.41をインストールする．
http://n4.nabble.com/compiling-CGAL-with-cmake-and-Visual-Studio-td1559113.html

・計算精度のためにGMPとMPFR，LEDAを使用
GMPとMPFRはVCだと適切でないことがあるから，CGALチームがプリコンパイルしたのを使う．
LEDAは下記のサイトで見つかる(ダウンロードは登録制)
http://www.algorithmic-solutions.com/leda/index.htm

・デモの多くを実行するにはその他もろもろのライブラリが必要

以上で必要なものはそろいます．


次に，cmakeでCGALを構築します．
CGALのルートフォルダで，cmake-gui .を打ち込んでConfigureを押して，1回目の処理が終わったらもう一回押してやります．
アクセス拒否やboostの問題を避けていれば，問題なく構築できるはずです．
構築できたら，今まで消えていたGenerateのボタンが押せるようになります．押しましょう．

生成されたslnファイルでソリューションを開きます．
開いたソリューションでDebugやReleaseでビルドするとlibファイルが作られます．
バージョン違いがなければ，エラーなくビルドできます．

これで，CGALのインストールは完了です．


プロジェクトを作って，インクルードファイルとリンクファイルを設定すれば，関数を使えるようになります．
インクルードファイルの設定で，CGAL直下のincludeフォルダとCGAL下のauxiliary\include，boost_1_41直下のincludeフォルダを追加すればコンパイルは通ります．
リンカの設定で，CGAL直下のlibフォルダとCGAL下のauxiliary\lib，boost_1_41直下のlibフォルダを追加すればリンカは通ります．
インストールは以上です．

ISVD2010気になる

ボロノイ図に関する国際会議ISVD2010の論文募集がcg-mlで流れていた．
ボロノイ図がらみというだけでマニア心をくすぐる．
ボロノイ図がらみならば，募集の領域をCGや計算幾何に限定していない，というのも気になる．

委員を見てたら某JAIST某浅野先生がいた．気になる．

シグラフあじあ終了

長い四日間が終わった．
今日の午前は，CGALのコースを半分聞いて，アトムの話を聞いた．午後は，Joe Rohdeの話を聞いて，モデリングの話を聞いた．
今日は最高だった．
まずCGALを使おうと思った．これはメッシュのブーリアンとかメッシュ構築とかができるライブラリで，汎用性やデータ構造の説明を聞いただけで，かなり魅力を感じた．使おう．
Joe Rohdeの話は物語の作り方についてだけど，内容は秘密なのだけれど，彼の話し方・プレゼンの仕方が最高だった．これまで聞いたプレゼンで一番素晴らしいと思った．ゆっくりと話し，強調したいことは強く，繰り返して言うし，言葉の韻を踏むし，聞きやすかった．こういうプレゼンを目指そうと思った．デザインされたプレゼンだった．
といった具合に最高でした．
このあと，仲間たちと飲みに行った．テーブルの座席時計回り順でいうと，某宮田先生と某高橋氏，某柿原氏，ざくさん，某せい氏とテーブルを囲んだ．
皆相変わらずのクオリティで，素晴らしかったです．


ザクさんのクオリティの高さを示した写真．

ということで，最高でした．

丸一日行動②

暦本先生が話すセッションでの質問は，話した内容とあまり関係ないような質問が多かったように思う．
聞き手によって，聞こえた内容が違うということなのかな，と思った．
このセッションで，昔のVRの展示や実験のビデオがいくつか流されていたけど，どれも反応が早かった．10年くらい前でも，バーコード読んで，オブジェクトを描画するだけなら，ほぼリアルタイムだということに驚いた．

論文で名前をみたことがあるけれど，それがどんな顔なのか知らないことが多い．
Shape Analysisセッションで，アメンタが質問するために名乗って，はじめてその人がアメンタだと知った．

とにかく疲れた．

丸一日行動

朝から夜までしぐらふあじあでした．
朝から昼は，論文セッションを聴講して，夜はComputer Animation Festicalを見てました．どれも素敵なセッションです．

昼間は，ポスター発表でしたので，作ったポスターの横に立っていました．
色んな人と話ができたのでよかったです．

夕方に来年の研究概要を宮田先生に見てもらって指摘をもらいました．時間もないので，すぐに直したいと思います．

色んなことしたので，とにかく疲れました．

HoppeとLoopと並んだ

SIGGRAPH ASIA 2009に参加して良かった．
ついにHoppeとLoopと並ぶことができました．
実力じゃなくて，物理的に並んだのですが．

感激しました．
握手とかしてもらったりして．
緊張して声とか手とか振るえてしまいました．

Technical Papers Fast Forwardで，某柿原氏と某せい氏が前のほうにいるからおいで，というメールを送りつけてきたので，私も前のほうに行って，彼らの後ろくらいに座ったの始まりです．
このとき，左には，某松藤氏がいて，ちょっと右のほうには某宮田先生がいました．すぐ右には，二人で向かい合いながら盛り上がっている白人男性がいました．そして，座った時には気付きませんでしたが，その二人がHoppeとLoopだったのです．
そのことに気付かずにいると，宮田先生から電話がかかってきて，君の右手に，Hoppeいるよ，と受話器から聞こえてきました．それを聞いて，右の男性の顔を見ると，確かにwebサイトで見てた顔と同じやつが座っているわけです．
驚きました．まじかよ．というような感じで．
そして，その話している相手はLoopでした．
どういうことよ，と思いました．

このチャンスを逃すのも惜しいので，ビビりながらも話しかけました．
大したことは話していないですが．
でも，こんな時のために，もっとまともな英語は常に準備しておかなきゃなとは思いました．

とにかくこのとき，ビッグネーム二人と肩を並べたわけです．
本当に感激しました．
並んでいることに気付いて電話してくる宮田先生にもさすがだな，と感銘を受けました．私もあんなふうに心に余裕がある大人になりたいです．

以上のように刺激的な時間を体感して，会社で素人に毛の生えたような人たちと話していてもいかんなと，このままじゃ腐ってしまうな，と思いました．
別の求めるなら，求めるで自分に厳しくいなければならないし，会社にいるなら，いるで強く生きなきゃならないな，と．

この出会いは，最高の思い出になりました．
ありがとう．



並んだ写真．あまりにも緊張したため，微妙な笑顔になってますね．

Predictice-Corrective Incompressible SPHはすばらしいと思う

Siggraph 2009で発表されたSolenthalerさんとPajarolaさんの流体シミュレーション手法の論文Predictice-Corrective Incompressible SPHはすばらしいと思った．
これまでのラグランジェパーティクルを用いた流体シミュレーションでは，濃度が指定した濃度より濃くなることがあり，濃度を一定に保てなかった．
この論文では，濃度が指定した値以上に高くなると，プレッシャーを積算し，次のタイムステップでの濃度が一定になるようにする．これにより潰れることがない流体シミュレーションができる．
この論文の素晴らしいところは，この計算方法を実装するために，これまで書いたSPHのコードをほとんど変更しないでよいことだと思う．
プレッシャーの一部を書き換えるだけで，そのまま動く．はず．すばらしい．
が，ひとつ問題を発見した．

LibTIFFのTIFFOpenの設定によるファイルアクセスの挙動の違い

LibTiffというライブラリを使ってプログラムすることがあるが，TIFFOpen関数のモードの設定の"r"と"r+"の違いではまったことがあったので，メモ．
環境はwindows xp pro, visual studio 2008 vc++．

TIFF* TIFFOpen(const char *filename, const char *mode)
この関数はfilenameであるファイルのアドレス，modeであるファイルアクセスのモードを設定することで，tiffファイルを開く機能が実装されている．
それで，modeの違いでちょっとした挙動の違いを示す．
"r"は読み込み，"w"は書き込み，"a"は追記とリファレンス(http://www.remotesensing.org/libtiff/man/index.html)のTIFFOpenの項目には記述してあるのだが，それ以外にも，"r+"などが設定できる．
今回は"r"と"r+"の若干の違いについて気付いたので，メモっとく．
TIFF* tiff = TIFFOpen(filename, "r");
TIFFClose(tiff);
とするとtiffはアクセスができなくなっているが，他のファイルからfilenameのファイルにアクセスしようとしてもアクセスエラーがでる．
たとえば，
TIFF* tiff = TIFFOpen("filename.tiff", "r");
TIFFClose(tiff);
system("del filename.tiff");
としても，filename.tiffは消えない．

つまりは，複数の実行ファイルをひとつの一括したようなプログラムでtiffを読む場合，この読み込みでは，一番最初に読んだやつだけがアクセス可能となり，それ以降ではCloseしているのにも関わらず，アクセス不可能となる．

一方，
TIFF* tiff = TIFFOpen("filename.tiff", "r+");
TIFFClose(tiff);
とした場合は，このfilename.tiffはどれからもアクセス可能になっている．つまり，
TIFF* tiff = TIFFOpen("filename.tiff", "r");
TIFFClose(tiff);
system("del filename.tiff");
とすると，filename.tiffは消える．

つまり複数の実行ファイルをまとめたようなプログラムでは，どいつもこいつもfilename.tiffにはアクセス可能となっている．Closeしているので，これが正しい形に思える．
セキュリティの違いで組みかたを変えればよいと思うが，この違いを知らないと，ずっとはまることになる．

ちなみに私は3時間くらいはまっていた．

openGLの光源の位置設定でありがちなミス

openGLでの光源の位置は，4つの要素配列を使って設定するべし．

void glLightfv(GLenum light, Glenum pname, const GLfloat *params)
この関数を使って光源を設定する．
lightにはGL_LIGHT0からGL_LIGHT7までの8個のライトの名前を入れる．
pnameにはGL_POSITION,GL_SPOT_DIRECTION,GL_AMBIENT,GL_DIFFUSE,GL_SPECULARなどから設定したいパラメータの名前を入れる．
paramsには設定するパラメータを入れる．

この設定でGL_POSITIONは4要素の配列で設定しないと思いもしないレンダリング結果が出てくることが多い．
3要素のx,y,zだけを設定した気になっていると，4要素目が0と見なされて，平行光源の光源の方向設定していることになる．
スポットライトを設定しているはずなのに，平行光となり，予想していた絵と全然違うということになる．
回避するためには，4要素[x,y,z,w]の配列を用意して，wは0以外の値（私は1を入れている）を入れてglLightfvを呼んでやる．

覚えている間は，ミスることはないけど，たまにしかライティングのコードを書かない私のような間抜けは，同じミスを繰り返す．
ちなみに私は，こんな設定の話をすっかり忘れて1日くらいはまったことがある．昨日くらいの話だけど．

Deformationテスト

デフォメーションのプログラムを組んでみた．
詳しいことは後でまとめるとして，まずは絵だけ．


初期状態の四面体です．ここからてっぺんの頂点を右にずらします．
以下に有限要素法を適用していない結果と適用した結果を並べます．


有限要素法を適用していない結果．てっぺんの頂点が移動しただけです．


有限要素法を適用した結果．てっぺんの頂点が移動するにつれて，他の頂点にも力が加わり移動しています．

今回の古典力学のシミュレーションはオイラー法で計算しているので，精度はすげー悪いはずです．論文に出ているようなタイムステップを実装しているわけではないので，そのうち誤差がたまって，四面体が爆発してしまいます．次はまともなタイムステップの計算をしてみます．

追記(2009/9/6):このデフォメーションはオイラー法でしか，古典力学の計算をできないかもしれない．ルンゲクッタ法では導関数を求めなければならない．しかしデフォメーションでは次のタイムステップの加速度は次の計算を行わない限り求まらない．つまり，ルンゲクッタ法に用いるための導関数は次の計算まで求まらないが，次の計算をするためには，次の導関数が必要となり，計算ができないということになる．だから，オイラー法で逐次計算していくしかない．ということは，みんな計算誤差なんて無視しているのか．

メモ
ちなみに，2007年のシグラフで公開されたBargteilさん達の論文に書いてあるFEMをそのまま実装しても全然うまくいきません．理由は，若干の記述を間違えている箇所があるからです．まじめに読んで実装しても，頂点は変な回転をしながら移動します．この論文が参照している2003年に公開されたTeranさん達のFinite Volume Methods for the Simulation of Skeletal Muscleか，2004年に公開されたIrvingさん達のInvertible Finite Elements For Robust Simulation of Large Deformationに記述されているFEMを実装するべきです．

F=XBはXが先に決まる

先日のF=XBで，deformation gradient Fと変形後の座標を意味するXのどちらが先に与えられるのかという疑問があったが，変形後の座標Xだということに気がついた．
Xは，力が加わり，加速度，速度が頂点を移動させた結果であるため，タイムステップがひとつ前の状態で決まる．変形後，Fが求まることで，さらに四面体の各頂点に力が伝搬するような仕組みになっている．
有限要素法の意味を考えてみるとすぐにわかるはずだった．

A Finite Element Method for Animating Large Viscoplastic Flowを再度読んでみる

Irvingさんたちの論文Invertible Finite Elements For Robust Simulation of Large Deformationを読んで，Bargteilさんたちの論文A Finite Element Method for Animating Large Viscoplastic Flowで記述されている有限要素法の触りの部分がわかった．

まず四面体を3x3行列で表現する．x₀とx₁，x₂，x₃を頂点とすると行列X_ij = x _ij - x_0jのように定義する．iは1,2,3でjは次元．
X₀を変形前の形状として，B=X₀^-1を基底行列とする．
変形後をXとして，deformation gradient FはF=XBとする．
SVDでF=UF^V^Tのように分解し，応力P^をP^=2u(F^-I)+λ Tr(F^-I)Iで計算する．uとλは素材の係数とのこと．
i番目の頂点に加わる力は，g_i=UP^V^T∑_jA_jN_jで計算する．

ここでトーシロー（素人）の私がわからないのが，変形後のXがあっての話なのか，変形させるためのFがあっての話なのか．このあと読み進めてみる．

Deformation Gradientに関するメモ

Deformation Gradientという言葉に遭遇したが，調べずに放置していたので，このあたりでメモ程度残しておく．

とりあえずググって出てきたのが，wikipedia．
http://en.wikipedia.org/wiki/Finite_strain_theory

これによると，デフォメーション前後の位置を見たとき，この近隣のサンプル点との相対的な位置の変化量を示した行列のことを指しているらしい．
具体的にはどうすればいいのかな．

剛体のアニメーションの計算では，物体内部のある任意の点は，近隣の点との相対的な距離は変化しない．しかし，デフォメーションでは，相対的な距離が変化する．準備として近隣ノードとの隣接となる点を見つけておきたい．

詳しくは，リンク先に書いてある．


Deformation Gradientが1のとき，体積は変わらないというようす．

BargteilさんらのViscoplasticモデル

BargteilさんたちのA Finite Element method for Animating Large Viscoplastic Flowの3章をさっと読んでみた．
前回読んで意味がわからず途中で断念したIrvingさんたちの論文を参考に，粘塑性の物体の変形と動きをモデル化している．今のところ，この論文に関しても，ほとんど理解できていない．
Irvingらの考えも理解できていないので，この分野の知見が足りないのだろう．

まぁ，よくわからなかったが，私の知見と照らし合わせて，何が書いてあったかまとめてみる．
3.1 Finite Element Methodでは，オブジェクトを四面体にして，四面体を変形するということが記述されていた．
四面体を行列X（ラージエックス）で表記している．この行列Xは3x3行列で，3次元上のあるノードから他の3つのノードの相対的な位置を表している．このとき，行をノードのインデックスで，列を次元としている．
シミュレーションをはじめるために，basis matrix（基底行列，と訳すのかな）を計算する．この基底行列βは初期のX0の逆行列とする．
このとき，与えられた変形状態がXのとき，Deformation gradient(変形勾配，と訳すかな)Fは，F=Xβで計算できる．これは，初期の四面体X0にDeformation gradient Fという変形行列をかけることで，変形後の座標Xが求まるという見慣れた計算．論文ではF=∂x/∂u=Xβと記述されている．uの意味が不明だが，きっと変形ベクトルを意味するんだろう．
あとの計算で使うため，FはF=UF^Vのように対角化する．
ここでは，Deformation gradient Fを計算するために，1st Piola-Kirchhoff stress(何これ?と思ってググったらwikipediaに書いてあった)を計算する．
http://en.wikipedia.org/wiki/Stress_(physics)
どうやら，素材と応力を関係させるとのことらしい．その応力P^は，P^ = 2μ(F^-I)+λTr(F^-I)I，で表されている．ここでは，λとμは材質のパラメーター．Trは行列のトレース．
それで力の計算は，g=UP^V^T∑ANで計算できるとのこと．ANは面法線の重み付き面積とする．
ここまでで，材質が関係する応力を計算し，座標変換の行列Fを求めたいということがわかった．たぶん1st Piola-Kirchhoff stressを知っていれば，ここまで読む必要なかったように思えた．

次に，3.2 Plasticity Modelについて読んだが，これまたよくわからない．
かかっている力eを可塑性e_pと弾性e_eに分けられると仮定する．体積の変化によるひずみによって加速が起こるとしている．（で，何よって感じがする）
可塑性のF_Pと弾性F_eの掛け算モデルは適切であり，F=F_P・F_eで表す．このとき，det(F_P)=1とすることで，体積の変化が起こらないことを保障するとのこと．
F^* =(det(F^))^-1/3 F^を計算するために対角化されたDeformation gradientのF^を使う．このとき，det(F^*)=1とする．そして，F^_P = (F^*)^γとする．γはねばりの係数．
つまり，(det(F^))^-1/3 F^に適当な数をかけた行列を座標変換の行列としてbasis matrix βにかけてやれば，よいということだろうか．

Deformation gradient Fの行列式が1になるように，かつ，弾性と可塑性に分けて計算すれば，それなりに見えるのではないかと，思ってきた．

Plasticモデル

今回読む論文は，Bargteilさん達の
A Finite Element Method for Animating Large Viscoplastic Flow
です．
PlasticをYahoo!辞書で引いてみますと，「プラスチックの」「ビニールの」等の素材を意味する形容詞と「可塑性の」という変形の性質を意味する形容詞が出てきます．まだ論文はちゃんと読んでいないので，はっきり言えませんが，変形について論じられているため，ここでのPlasticは「可塑性の」なのだと思います．
この論文を読み切った後は，可塑性の変形を実装しようと思います．

この論文は読めない

Volume conserving finite element simulations of deformable modelsを読んでいましたが，全然内容を理解できません．ですので，次の論文読みに移ろうと思います．
別の機会に，運悪くこの論文を読む羽目にならないように，ここでメモって残しておこうと思います．

この論文に書かれている有限要素法は，2章 Time Discretizationと3章 Spatial Discretization，4章 Collisions and Contactを理解しなければ実装はできない気がしているのですが，正直2章の時点でつまずいています．この章には，1タイムステップの処理内容が記述されています．1から9ステップで構成されていて，位置と速度を計算します．次のタイムステップの位置は現在の位置と速度（なんか色々ややこしい処理したあとの速度）から計算しますが，次のタイムステップの速度は，式を見る限り(ステップ8の式に7の式を代入することで)，加速度にタイムステップをかけた値になるはずです．この速度では，現在の速度は失われてしまいますので，加えられた力の分しか計算しません．しかし，その加速度も明確に計算式が記述されておらず，同著者の論文をちょっと漁ってみても，加速度の計算の記述はありません．何を計算すればよいのかさっぱりわかりません．
体積を一定にするために，位置を決定する速度を適当な値に調整し，次の速度を適当な値に調整するのは，わかりましたが，これらの計算方法がわからないので，どうしたもんか，といったところです．論文の売りは，Divergence freeとone-ringの体積の保持でしょうから，動きの計算は適当でいいのかもしれませんが，この論文に書いてあることは，実装できないので，もう読むのをやめます．次の論文に移ろうかと思います．

次はElasticモデル

次は，Elasticモデルについて，IrvingさんらのVolume conserving finite element simulations of deformable modelsを読もうと思います．

球の衝突判定実験

Bridsonさんらの論文に書いてある衝突判定法を用いてFEM法で作った球で衝突判定の実験をし，この衝突判定法で判定ができることを確認した．
実験のために以下の準備をした．
・球のメッシュの頂点を移動可能にし，速度と加速度を設定できるようにした．高い精度が必要ないので，オイラー法でシミュレーション．
・球の下に面を配置し，衝突判定をする．ここでのリアクションは，衝突していると判定された頂点は動かないようにする．
・球に下向きに加速度をかけ続ける．なんとなく重力加速度を意識した．
この実験の結果が上の数枚の画像．左側はメッシュのエッジの描画で，右側はそれに面をつけてレンダリングした結果．

ここでは，Bridsonさんらの論文の説明が不親切な部分があるので，少しメモを残しておく．

この衝突判定法は，三角形x₁x₂x₃と点x₄の衝突を検知する手法．というか，結構一般的な手法なので，Bridsonさんらの論文を読まなくても，CGの数学関係の本には書いてあること．

まず，点x₄と三角形x₁x₂x₃を含む平面の距離を計算し，接触の距離h以下の点x₄を判定する．準備として2個のベクトルを計算する．1個は，頂点x₃（実はx₁でもx₂でも十分)から点x₄までのベクトルを計算したベクトル．もう1個は，三角形の法線ベクトルnを計算したベクトル．これら2ベクトルの内積は平面からの距離であり，条件(x₄-x₃)・n < hが満たされるとき，x₄と三角形x₁x₂x₃の平面は接触しているとする．

次に，点x₄を平面上に射影した点x'₄が，三角形x₁x₂x₃の内側か判定する．内側ならば，三角形と点は接触している．次の連立方程式を解き，その解を条件式にはめることで判定できる．
(x₁ - x₃)・(x₁ - x₃)w₁ + (x₁ - x₃) ・ (x₂ - x₃)w₂ = (x₁ - x₃)・(x₄ - x₃)
(x₁ - x₃)・(x₂ - x₃)w₁ + (x₂ - x₃) ・ (x₂ - x₃)w₂ = (x₂ - x₃)・(x₄ - x₃)
1 = w₁ + w₂ + w₃
点x'₄は，x'₄ = w₁ x₁ + w₂ x₂ + w₃ x₃である．また，w1とw2がw1>=0, w2>=0, 1-w1-w2>=0を満たしていれば，点x'₄は三角形の内側である．
以上のように，接触の判定を行う．

ちなみに，今回のリアクションは，適当に設定したため，論文のモデルに則っていない．