いんちょき通信

動画ファイル作成をスクリプトにしてみた

以前に書いた「OpenGLのアニメーションを動画ファイルにしてみた」でのPPMファイル→JPEG→AVI→MPEGという流れが覚えていられないので、シェルスクリプトにしてみた。

1. #! /bin/bash
3. MOGRIFY=mogrify ←パスが通っていなければフルパスで記載
4. FFMPEG=ffmpeg ←パスが通っていなければフルパスで記載
6. IMGDIR=.
8. FPS=20
10. while getopts "d:f:" OPT; do
11.   case $OPT in
12.     \?) ERROR=11; break;;
13.     d) IMGDIR=$OPTARG;;
14.     f) FPS=$OPTARG;;
15.   esac
16. done
17. shift $(( $OPTIND - 1 ))
19. FILENAME=$1
21. $MOGRIFY -format jpeg $IMGDIR/*.ppm
22. $FFMPEG -r $FPS -i ${IMGDIR}/%d.jpeg -vcodec mjpeg -sameq　${IMGDIR}/${FILENAME}.avi
23. $FFMPEG -i ${IMGDIR}/${FILENAME}.avi -b 900 -ab 96 ${IMGDIR}/${FILENAME}.mpeg
25. exit 0

-dオプションで画像が置いてあるディレクトリ、-fオプションでフレーム/secが指定できる。出力ファイル名を最後につけて実行する。実行方法は以下のような感じ

$ ./mkmov -d ./img -f 20 outmovie

mkmovをダウンロード

微少時間刻みの衝突を考えてみる

過去3回、質点の運動(等速直線運動、等加速度運動、放物運動)を微少時間刻みで考えたが、今回は衝突を考えてみようと思う。衝突のモデルとしてKelvin–Voigtモデルというものを導入する。

 出典：Wikipedia

今回はひとまず壁との衝突を考えるが、要素と壁が衝突したときに、上図のモデルが接触箇所に存在すると仮定して計算する。接触の判断は要素中心と壁との距離が、要素半径より小さければ接触していることになる。つまり、

これを満たすとき、要素と壁は衝突している。ここで、は要素の座標、は壁の座標、は要素の半径。

とし、をモデルバネの縮みと見なす。このとき、ばね定数をとすると、バネによる力は

と表せる。また、ダッシュポッドの粘性係数をとすると、ダッシュポッドによる力は速度に比例し、速度と逆方向に働くので、

と表せる。ここで、は現在の速度である。

とをあわせたものが、要素に働く力になるので、

との向きは考える必要がある。

現在の力がわかったので、要素の質量をとすると、要素の現在の加速度は

となる。

ここまでできれば、後は等加速度運動の時と同じ処理になる。この部分をコードにすると以下のような具合。

1.     if(wp1 - xc < r){ // 右壁
2.       double dx = r - (wp1 - xc); // 壁と要素の重なり
3.       fc = -dx * k - c * vc;
4.     }
5.     if(xc - wp2 < r){ // 左壁
6.       double dx = r - (xc - wp2); // 壁と要素の重なり
7.       fc = dx * k - c * vc;
8.     }
9.     ac = (fc + fp) / (2 * m);

ばね定数k=10000、粘性係数d=20としたときの結果は次のようになる。

 

動画ファイルをダウンロード

放物運動を微少時間刻みで考えてみる

前回、前々回と1次元の運動をシミュレートしたが、今回は2次元に拡張して放物運動をシミュレートした。要は、水平方向を等速直線運動、鉛直方向を等加速度運動として計算すればいいということになる。等速直線運動は等加速度運動の加速度が"0"のときともいえるので、コードとしては水平、鉛直ともに等加速度運動としている。

1. #include <fstream>
2. #include <cstdlib>
4. using namespace std;
6. int main(int argc, char* argv[])
7. {
8.   double dt = 0.0001; // 1ステップの微少時間
9.   double po[2] = {0.0, 400.0}; // 初期座標
10.   double pc[2]; // 現ステップの座標
11.   double pp[2]; // 1ステップ前の座標
12.   double vc[2]; // 現ステップの速度
13.   double vp[2]; // 1ステップ前の速度
14.   double a[2] = {0.0, -9.81}; // 加速度
15.   int totalStep = 100000; // 総ステップ数
16.   double p[2][200]; // 保存用座標
18.   ofstream outf("./output.txt"); //結果出力ファイル
19.   int count = 0;
20.   p[0][0] = pp[0] = po[0];
21.   p[1][0] = pp[1] = po[1];
22.   vp[0] = 30.0; // 水平方向初速
23.   vp[1] = 0.0; // 鉛直方向初速
24.   for(int i = 1; i < totalStep; i++){
25.     vc[0] = vp[0] + a[0] * dt;
26.     vc[1] = vp[1] + a[1] * dt;
27.     pc[0] = pp[0] + (vc[0] + vp[0]) * dt / 2;
28.     pc[1] = pp[1] + (vc[1] + vp[1]) * dt / 2;
29.     count++;
30.     if(count == 500){
31.       count = 0;
32.       p[0][i/500] = pc[0];
33.       p[1][i/500] = pc[1];
34.     }
35.     pp[0] = pc[0];
36.     pp[1] = pc[1];
37.     vp[0] = vc[0];
38.     vp[1] = vc[1];
39.   }
40.   for(int i = 0; i < 200; i++){
41.     outf << p[0][i] << endl;
42.     outf << p[1][i] << endl;
43.   }
45.   return EXIT_SUCCESS;
46. } 

各時間刻みで水平方向についてと鉛直方向について同じ計算を行っているが、水平方向の加速度は"0"なので、放物運動となる。この結果は以下のようになる。

 

やっぱり、まだまだ地味だ。

動画ファイルダウンロード

等加速度運動を微少時間刻みで考えてみる

前回の等速直線運動に続き、等加速度運動を考えてみる。要は物体の落下とかに当たるもの。これは次のように表される。

前回と同じく、微少時間刻みで考えると、

となるので、これを計算してみた。コードにすると次のようになる。

1. #include <fstream>
2. #include <cstdlib>
4. using namespace std;
6. int main(int argc, char* argv[])
7. {
8.   double dt = 0.0001; // 1ステップの微少時間
9.   double xo = 0.0; // 初期座標
10.   double xc; // 現ステップの座標
11.   double xp; // 1ステップ前の座標
12.   double vc; // 現ステップの速度
13.   double vp; // 1ステップ前の速度
14.   double a = 0.02; // 加速度
15.   int totalStep = 100000; // 総ステップ数
16.   double x[200]; // 座標保存用配列
18.   ofstream outf("./output.txt"); // 出力ファイル
19.   int count = 0;
21.   x[0] = xp = xo;
22.   vp = 0.1; // 初速
24.   for(int i = 1; i < totalStep; i++){
25.     vc = vp + a * dt;
26.     xc = xp + (vc + vp) * dt / 2;
27.     count++;
28.     if(count == 500){
29.       count = 0;
30.       x[i/500] = xc;
31.     }
32.     xp = xc;
33.     vp = vc;
34.   }
35.   for(int i = 0; i < 200; i++){
36.     outf << x[i] << endl;
37.   }
39.   return EXIT_SUCCESS;
40. }

この結果を先ほどと同じく動画にすると、次のようになる。

  

前回と今回を組み合わせると、放物運動も表現できるはず。

動画ファイルダウンロード

等速直線運動を微少時間刻みで考えてみる

力学の最初と言えば、等速直線運動。力が加わっていない質点は静止するか、等速直線運動をするというやつの等速直線運動。ある時間の質点の位置は初期位置に速度×時間を加えることで求められる。

この式から任意の時間における質点の位置を求めることができる。

今回はあえて、微少時間刻みで質点の移動を計算してみた。つまり、を初期位置とし、を速度、を微少時間すると、開始からdt時間後の位置はとなる。これから、任意のに対して、と表現できる。これを使って順番に計算していこうということ。

コードにするとこんな感じになる。

1. #include <fstream>
2. #include <cstdlib>
4. using namespace std;
6. int main(int argc, char* argv[])
7. {
8.   double dt = 0.0001; // 1ステップの微少時間
9.   double xo = 0.0; // 初期座標
10.   double xc; // 現ステップの座標
11.   double xp; // 1ステップ前の座標
12.   double v = 0.01; // 速度
13.   int totalStep = 100000; // 総ステップ数
14.   double x[200]; // 座標保存用配列
16.   ofstream outf("./output.txt"); //結果出力ファイル
17.   int count = 0;
19.   x[0] = xp = xo;
20.   for(int i = 1; i < totalStep; i++){
21.     xc = xp + v * dt;
22.     count++;
23.     if(count == 500){
24.       count = 0;
25.       x[i/500] = xc;
26.     }
27.     xp = xc;
28.   }
30.   for(int i = 0; i < 200; i++){
31.     outf << x[i] << endl;
32.   }
34.   return EXIT_SUCCESS;
35. }

単位系をMKSとすると、原点からスタートして、速度1cm/sで等速直線運動する質点の位置を、微少時間0.0001秒ごとに10秒間ぶん計算していることになる。そして出力は500回に1回なので、0.05秒ごとの位置をファイルに記録している。

ちなみに、この結果をアニメーションにしたものが「OpenGlのアニメーションを動画ファイルにしてみた」の記事になる。ちょっと円の半径は大きくしすぎたと思う。

OpenGLのアニメーションを動画ファイルにしてみた

OpenGLで作成したアニメーションを動画ファイルとして保存できないと、素材として使いにくい。とりあえずGoogle先生に聞いてみると、「ffmpegを使ってOpenGLのアニメーションから動画を作る方法(Mac)」というページが出てきたので、まねしてみた。

1. バッファをファイルに書き出す。
2. 書き出したファイルをjpegとかに変換する。(ImageMagick)
3. 2.のファイルをつなげて動画にする。(ffmpeg)

という流れになるようだ。先ほどのページではMac Portsを使っているが、あれもこれもと入れるのは面倒なので、Gentoo Prefixで何とかした。以前の記事でGentoo PrefixのMaskを外そうとしていたのは、このためだったというわけ。

はじめと言うことで、円を等速直線運動させるOpenGLアニメーションを作成し、描画関数(display関数)内でglutSwapBuffersする直前に、バッファをPPMとしてファイルに落とした。バッファをファイルに書き出す実装は「ffmpegを使ってOpenGLのアニメーションから動画を作る方法(Mac)」に詳しく書かれている。

PPM→jpeg変換をImagemagickで行い、ffmpegで連結しaviファイルとした。この部分も上記のページ通りに実行している。

ただ、aviだとサイズが大きくなりがちなので、さらにmpegに変換した。これもffmpegを使っている。コマンドはこんな感じ

$ ffmpeg -i out.avi -b 900 -ab 96 sample.mpeg

別にQuickTimeでmovとして書き出してもいいと思う。

できあがったものはこんな感じ。

 

これだけだとなんてことないけど、まあできたってことで。

ブログに数式を埋め込んでみる

ブログに数式を表示したいと思い、調べてみると奥村先生のブログに行き当たった。Google Chart APIというものが紹介されていたので、早速試してみる。

<img src="http://chart.apis.google.com/chart?cht=tx&amp;chl=x=\frac{-b\pm\sqrt{b^2-4ac}}{2a}" alt="" />

ようは"chl="の後ろにTeXコマンドで数式を書けばいいようだ。

<img src="http://chart.apis.google.com/chart?cht=tx&amp;chl=\b{\nu}_c=\b{\nu}_p%2B\b{a}dt" alt="" />

URL中に使えない文字は%16進2桁で指定する必要がある。上記の例だと"+"が"%2B"になっている。詳細はMathematical (TeX) Formulasに載っている。

ちなみにこのGoogle Chart APIは名前の通り、本来はチャートをWEB上に表示するためのものらしい。

Gentoo Prefixにautounmaskとかを導入してみる

Twitterで「#gentoo #prefix では、 #ImageMagic もマスクされてるじゃないか！ #MacPortsも導入しようかな。」とつぶやいたら、"autounmask"というものを紹介されたので、ついでにちょっといろいろ調べてみた。

Google先生に問い合わせて、最初にたどり着いたところはここ。せっかくなので、色々試してみた。

まずは、/etc/make.confのCFLAGSの設定から。どうも最近のgccには"-march=native"というものが存在しているらしい。参考サイトでは分散コンパイルを考慮して、展開したものを使っているが、うちでは別にdistccとかは使わないので、そのまま設定してみることにした。

その後は、MAKEOPTSやFEATURESの説明が続くが、そんなに頻繁にコンパイルしないし、ということで放っておくことにした。

USEフラグも色々書いてもごちゃごちゃしてよくわからなくなるだけなので、必要があるときに必要なものをの精神で特に追加はしなかった。

make.confの現状はこんな感じ。

1. USE="unicode nls cxx examples bzip2 zlib gzip png jpeg pdf truetype
2. "CFLAGS="-march=native -O3 -fomit-frame-pointer -pipe"
3. CXXFLAGS="${CFLAGS}"

たいした書いていません。

パッケージ探しが楽になりそうなので、portage-utilsは導入することにした。

# emerge portage-utils

でなんてことなく、インストール完了...と思ったら、"-march=native"なんて知らないよとエラーで終了してしまった。Gentoo Wikiをみると、GCC4.2.3以上で使えると書いてあるが、いまインストールされているgccは4.2.1だった。"make.conf"の"-march=native"を"-march=nocona"に変更して再度トライしたら今後はうまくいった。

本題からずれてしまっているが、色々導入するにしてもGentoo PrefixではMaskされているパッケージが多いので、さっさと本題のautounmaskを導入する。

# emerge autounmask

とすればOKみたい。依存パッケージがいくつかあるが、すんなり入ってくれた。

さっそく当初の目的のImageMagickをautounmaskでMaskを外してみることにした。

# autounmask -n media-gfx/imagemagick-6.6.5.6

"-n"オプションをつけると、すべてのバージョンのMaskを外してくれるようだ。実行結果では、このファイルに追加したよというのも表示される。インストールも無事できた。

バックグランドでファイル出力をするクラスを作ってみた

大量のデータをファイル出力するときに、バックグランドでファイル書き出しをしてほしいなということで、そんなクラスを作ってみた。

仕様は次の通り。

• ファイル名を渡してインスタンスを生成する。
• ライト関数に文字列を渡すと、キューに登録する。
• インスタンスは別スレッドで動作し、キューに文字列が入ったらファイルに書き出す。
• TERMINATEコマンドを送ると、キューに残っているものをすべて書き出し、終了する。
• SUSPENDコマンドを送ると、書き出し処理を中断する。
• RESUMEコマンドを送ると、書き出し処理を再開する。

クラスの定義は次の通り。

1. class DBGWriter{
2. private:
3.   boost::thread th_; //ライタ関数が動作するスレッド
4.   std::queue<std::string> strQueue_; //書き出し文字列のキュー
5.   boost::mutex mutex_;
6.   boost::condition_variable_any condition_;
8.   std::string outputFileName_; //出力ファイル名
9.   std::ofstream ofStream_; //ファイル出力ストリーム
11.   bool isWritable_; //書き込み可のフラグ
12.   bool isSuspend_; //サスペンドのフラグ
14.   DBGWriter(){};
15.   void writer(); //ライタ関数
16. public:
17.   DBGWriter(std::string outputFileName);
19.   void write(std::string str); //書き出し文字列登録関数
20.   void sendCmd(DBGW_CMD cmd); //コマンド送信関数
21. };

ライタ関数は次のようになっている。

1. void DBGWriter::writer()
2. {
3.   boost::mutex::scoped_lock lock(this->mutex_);
5.   while(this->isWritable_){
6.     while(this->isSuspend_){
7.       this->condition_.wait(lock);
8.     }
9.     while(this->strQueue_.empty()){
10.       this->condition_.wait(lock);
11.     }
12.     this->ofStream_ << this->strQueue_.front() << std::flush;
13.     this->strQueue_.pop();
14.   }
15.   while(!this->strQueue_.empty()){
16.     this->ofStream_ << this->strQueue_.front() << std::flush;
17.     this->strQueue_.pop();
18.   }
19. }

5〜14行目が通常の書き込み可能状態の処理となっている。6〜8行目、9〜11行目でそれぞれサスペンド時とキューが空の時の待ち状態に入る。キューに文字列が入ったら、12行目でファイルに出力し、13行目でキューから取り除いている。15〜19行目はTERMINATEコマンド発行時の処理となっている。キューが空になるまで、キューの中身をファイルに書き出している。

RESUMEコマンドが発行されるか、キューに文字列が入れられると、boost::condition_variable_any::notify_all()が発行され、待ち状態を解除するようになっている。

他はたいしたことないので、割愛する。

全体のソースはこちら→ bgw.tar.gzをダウンロード

アクティブなスレッドをプールするクラスを作ってみた

スレッドを事前に生成しプールしておくと、実際に使用するときにスレッド生成のオーバーヘッドを避けることができる。このようなパターンをWoker threadパターンと言うらしいので、実際に実装してみた。

仕様はこんな感じ。

• スレッドプールを作成するときにプールサイズを指定できる。
• スレッドプールを作成するときにプールサイズを指定しない場合は、デフォルト値のスレッドが用意される。
• リクエストメソッドでスレッドにタスクを渡す
• クローズメソッドで残っているタスクが完了次第、すべてのスレッドを削除する。
• スレッド追加メソッドで引数個のスレッドをプールに追加する。
• スレッド削除メソッドで引数個のスレッドをプールから削除する。引数が現在のスレッド数を超えている場合はfalseを返す。
• スレッドを削除するときはアイドル状態のスレッドから削除する。

実装はboost C++ Librariesのthreadを使用した。

クラスの定義は次の通り。

1.   class DThreadPool : boost::noncopyable  {
2.   private:
3.     std::queue<boost::function0<void> > queue_; //タスクを入れるキュー
4.     std::map<boost::thread::id, boost::thread *> threadMap_; //スレッドIDとスレッドのマップ
5.     std::map<boost::thread::id, bool> idleMap_; //スレッドIDをスレッド状態のマップ
6.     boost::thread_group thGrp_;
8.     boost::mutex mutex_;
9.     boost::condition_variable_any condition_;
11.     bool isAcceptable_; //タスク受付可能かどうかのフラグ
13.     void run(); //スレッドで実行するタスク待ち関数
14.     void setIdleFlag(boost::thread::id threadID, bool idleFlag);
15.   public:
16.     DThreadPool();
17.     DThreadPool(int size);
18.     ~DThreadPool();
20.     void request(boost::function0<void> func); //タスクをキューに入れる関数
21.     void close(); //クローズ関数
22.     int size();
24.     void addThread(int n); //スレッド追加関数
25.     bool removeThread(int n); //スレッド削除関数
26.   };

スレッドの生成は次のようになっている。

1. DThreadPool::DThreadPool(int size)
2. {
3.   this->isAcceptable_ = true;
4.   for (int i = 0; i < size; i++){
5.     boost::thread* th = thGrp_.create_thread(boost::bind(&DThreadPool::run,   this));
6.     this->threadMap_[th->get_id()] = th;
7.   }
8. }

まず、3行目でタスクを受付可能とする。4〜6行目でsize個のスレッドにrun関数を渡して生成する。それぞれのスレッドをスレッドIDをキー、スレッドへのポインタを値としてマップに入れる。

run関数は次のようになっている。

1. void DThreadPool::run()
2. {
3.   while (isAcceptable_){
4.     boost::mutex::scoped_lock lock(this->mutex_);
6.     setIdleFlag(boost::this_thread::get_id(), true);
7.     while (this->queue_.empty()){
8.       this->condition_.wait(lock);
9.     }
10.     setIdleFlag(boost::this_thread::get_id(), false);
12.     boost::function0<void> func = this->queue_.front();
13.     this->queue_.pop();
14.     func();
15.   }
16.   while (!this->queue_.empty()){
17.     boost::function0<void> func = this->queue_.front();
18.     this->queue_.pop();
19.     func();
20.   }
21. }

3〜15行目は受付可能な状態での処理。まず、ロックをかけキューにタスクが入るの待機する状態になる。このとき自分のスレッドにアイドルフラグをセットする。キューに何かが入ってくると、アイドルフラグを外し、キューから取り出して実行する。16〜20行目は受付不可の状態での処理。キューが空になるまで、キューからタスクを取り出し実行する。キューが空になった時点で、関数が終了する。

タスクをキューに入れる処理は次のようになっている。

1. void DThreadPool::request(boost::function0<void> func)
2. {
3.   boost::mutex::scoped_lock lock(this->mutex_);
4.   this->queue_.push(func);
5.   this->condition_.notify_all();
6. }

ただロックかけて、キューにプッシュして、ロックを外しているだけ。

クローズ処理は次のようになっている。

1. void DThreadPool::close()
2. {
3.   this->isAcceptable_ = false;
4.   for(std::map<boost::thread::id, boost::thread*>::iterator iter =  this->threadMap_.begin();
5.         iter != this->threadMap_.end();
6.         ++iter){
7.     (iter->second)->interrupt();
8.     (iter->second)->join();
9.     this->thGrp_.remove_thread(iter->second);
10.   }
11. }

まず、3行目で受付不可状態にする。4〜10行目ですべてのスレッドに割り込みをかけて、終了を待つ。終了したスレッドはプールから削除する。デストラクタ内でこの関数を呼ぶことで、安全に終了することができる。

スレッドの追加は次の通り。

1. void DThreadPool::addThread(int n)
2. {
3.   if(!this->isAcceptable_){
4.     this->isAcceptable_ = true;
5.   }
6.   for (int i = 0; i < n; i++){
7.     boost::thread* th = thGrp_.create_thread(boost::bind(&DThreadPool::run,   this));
8.     this->threadMap_[th->get_id()] = th;
9.   }
10. }

要はコンストラクタ内でやっていることと同じことをしている。違う点は受付状態を確認してから受付可能にするところ。

スレッドの削除は次の通り。

1. bool DThreadPool::removeThread(int n)
2. {
3.   if(n > this->thGrp_.size()){
4.     return false;
5.   }
6.   int count = 0;
7.   while(count < n){
8.     for(std::map<boost::thread::id, boost::thread*>::iterator iter =  this->threadMap_.begin();
9.           iter != this->threadMap_.end();
10.           ++iter){
11.       if(count < n){
12.         if(this->idleMap_[iter->first]){
13.           boost::thread* th = this->threadMap_[iter->first];
14.           th->interrupt();
15.           th->join();
16.           this->thGrp_.remove_thread(th);
17.           this->idleMap_.erase(iter->first);
18.           count++;
19.         }
20.       }
21.     }
22.   }
23.   if(this->thGrp_.size() == 0){
24.     this->isAcceptable_ = false;
25.   }
26.   return true;
27. }

まず、3〜5行目で削除するスレッド数が現在のスレッド数より多いかをチェックする。多い場合はfalseを返す。次に7〜20行目で削除する個数分だけ、アイドルフラグが立っているスレッドを探して、割り込み処理→終了待ち→削除としていく。このときスレッド状態マップから該当するスレッドIDを削除する。23〜25行目でサイズが0になった場合は受付不可にする。

ソースはこちら→ threadpool.tar.gzをダウンロード

参考サイト

• Let's Boost
  
• metl日記 Boost によるデザインパターン - Worker Thread