Gnuplotによるカメラ画像のPSF推定

カメラ画像のボケぐあいであるPSF(Point Spread Function、点拡がり関数)を求める１つの手法をまとめる。 Gnuplotによる区分関数の最小二乗法 を用いるとプログラミングの必要なく非常に簡単にPSFを求めることができる。 ここで仮定するPSFは、次の式で表現されるピルボックス関数（Pillbox function）である。

PSFとしてピルボックス関数を仮定する

\[ R(x, y, r)= \left\{ \begin{array}{cl} \frac{1}{\pi r^2} & x^2+y^2 \leq r^2\\ 0 & {\rm otherwise} \end{array} \right. \]

階段関数のテストパターン

テストパターンとして次の図のような白黒2色のパターンを撮影する。 これは水平方向に見れば階段関数（step function, staircase function）と見なすことができる。

\[ S(x,y)= \left\{ \begin{array}{cl} 1 & (x \leq 0)\\ 0 & (x < 0) \end{array} \right. \]

観測パターンとそのモデル

上記のパターンがピンぼけて下図（原画の4倍スケール）のようなパターンが観測される。

ピルボックス関数とテストパターンの畳み込み積分を考えると、 ピルボックス関数の一部がテストパターンの黒い領域に重なり、その残りの部分が白い領域に重なる瞬間がある。 この２つ重なり領域のピルボックス関数の体積比が画像の輝度として観測されていると考えることができる。 こう考えると観測画像の水平方向の輝度は下式で表されることが知られている。 $-1 \leq t \leq 1$の範囲がピルボックス関数とテストパターンの黒白境界が重なる範囲である。

\[ f(x)= \left\{ \begin{array}{ll} I_{\rm min} & t<-1 \\ g(t)(I_{\rm max} - I_{\rm min}) + I_{\rm min} & -1 \leq t \leq 1 \\ I_{\rm max} & t>1 \end{array} \right. \] \[ g(t) = \frac{1}{2} + \frac{1}{\pi}(t \sqrt{1-t^2}+\sin^{-1}t) \] \[ t = \frac{x - c}{r} \]

PSFのパラメータ推定

上記関数のパラメータは何か？ パラメータは白い領域の輝度$I_{\rm max}$、黒い領域の輝度$I_{\rm min}$、中心位置$c$、ボケ幅$r$だけど、 PSFとして意味のあるパラメータは$r$だけである。 4つのパラメータを求めてそのうちボケ幅$r$だけがピルボックス関数の形を決めるパラメータとなる。

一方データとしては、画像の水平方向に走査したときの輝度値を用いる。 垂直方向に輝度値を平均すれば画像のノイズの影響を軽減することができる。 このデータに対して先ほどの関数をGnuplotを用いてフィッティングする。 フィッティングされた関数とデータの例を下に示す。

reset

# 関数の定義
t(x,c,r)=(x-c)/r
g(x,c,r)=0.5+(t(x,c,r)*sqrt(1-t(x,c,r)**2)+asin(t(x,c,r)))/pi
f(x,c,r,imax,imin)=(t(x,c,r)<-1) ? imin :((t(x,c,r)>1) ? imax : g(x,c,r)*(imax-imin)+imin)
 
# 初期値の設定
c=50
r=25
imax=95
imin=5
 
# フィッティング
fit f(x,c,r,imax,imin) 'sample.txt' via c,r,imax,imin
 
# データと関数描画
set xr[-5:105]
plot 'sample.txt' with points pt 3
replot f(x,c,r,imax,imin) linewidth 3

Gnuplotの出力は以下の通り。 以下の結果からボケ幅は$r=21.298$画素であることが分かる。

 Iteration 0
 WSSR        : 587.799           delta(WSSR)/WSSR   : 0
 delta(WSSR) : 0                 limit for stopping : 1e-005
 lambda   : 0.859572
 
initial set of free parameter values
 
c               = 50
r               = 25
imax            = 95
imin            = 5
/
 
 Iteration 1
 WSSR        : 68.9751           delta(WSSR)/WSSR   : -7.5219
 delta(WSSR) : -518.824          limit for stopping : 1e-005
 lambda   : 0.0859572
 
resultant parameter values
 
c               = 49.3057
r               = 21.7139
imax            = 94.3184
imin            = 6.75013
/

　　　：
　　（省略）
　　　：
 
 Iteration 5
 WSSR        : 63.5581           delta(WSSR)/WSSR   : -8.02681e-007
 delta(WSSR) : -5.10168e-005     limit for stopping : 1e-005
 lambda   : 8.59572e-006
 
resultant parameter values
 
c               = 49.278
r               = 21.298
imax            = 94.2626
imin            = 6.80794
 
After 5 iterations the fit converged.
final sum of squares of residuals : 63.5581
rel. change during last iteration : -8.02681e-007
 
degrees of freedom    (FIT_NDF)                        : 96
rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 0.813673
variance of residuals (reduced chisquare) = WSSR/ndf   : 0.662063
 
Final set of parameters            Asymptotic Standard Error
=======================            ==========================
 
c               = 49.278           +/- 0.07501      (0.1522%)
r               = 21.298           +/- 0.1477       (0.6935%)
imax            = 94.2626          +/- 0.1459       (0.1548%)
imin            = 6.80794          +/- 0.1517       (2.228%)
 
correlation matrix of the fit parameters:
 
               c      r      imax   imin   
c               1.000 
r              -0.013  1.000 
imax            0.432  0.325  1.000 
imin            0.449 -0.338 -0.013  1.000 

参考

• Gnuplotによる区分関数の最小二乗法
  http://buaiso.blogspot.jp/2013/11/gnuplot.html
• B. Okumura, M. Kanbara, and N. Yokoya: "Augmented reality based on estimation of defocusing and motion blurring from captured images", Proc. IEEE and ACM Int. Sympo. on Mixed Augmented Reality (ISMAR 06), pp. 219-225, 2006.

左翼先生の栄えた時代

自分だけの稀な体験なのか、それとも全国的に似たような傾向があるのか全く分からないけど、 自分が義務教育時代に経験した、というより違和感を感じた、出来事をここで紹介したい。

はっきり言って高校に入るまで心から尊敬できる先生は居なかった。 その最大の原因は、先生があまり勉強していないことが、生徒である自分にもなんとなく分かったから。 専門の教科ですら、質問しに行けばろくに答えられずにごまかす先生も多かった。 生徒は子供だから大丈夫だろうというのは全くの間違いで、その時はそういうものだと思っても必ず覚えていて、 大人になってやっぱりあの先生があの簡単な質問に答えられなかったのは、おかしいじゃないかと気づくもの。 自分の場合も、そういえば、今考えてみたら、当時の先生は変な人が多かったなという印象をもつ。 今思えば、今の自分の知識からすれば、当時自分が感じていた以上にその先生方は変な先生だった。 どう変かと言えば極端に左に偏った思想をもち、ろくに勉強もせずに自分の思い込みを、 めちゃくちゃな論理で押し付けてくる。そんな先生があちこちに居た時代だった（今は知らない）。

めちゃくちゃな論理

小学校時代の遠足や運動会の時だけお弁当を持って行くことになっていた。 ただ、なぜかお弁当の中身に関して変なルールがあった。

お弁当はおにぎり２個と漬物以外持ってきてはならない。

この理由は今の人には全く分からないだろう。 要するに貧乏な人もいるからその格差が目立ってはいけないので、 全員貧乏な飯にすべしということだった。 いかにも平等が大好きな左翼先生が考えそうなことだ。 当時の自分は、お弁当も持ってこれないような貧乏がいるかよと思った。 確かに芸人のダウンタウンの子供時代の話を聞けば、塩むすびとみょうがしか持ってこない生徒も居たとか。 それに、共働きでお弁当を作ってくれない親もいるかもしれない。 そういう配慮かもしれないけど、全員を最下層に合わせるという発想が貧困というか、だったらパンでもいいじゃないかと思ってしまう。 子供の頃は人との違いを非常に気にして、ちょっと変わった人がいたらいじめてしまうものである。 でも、それは何も貧乏かそうでないかだけじゃなくて、体型だとか、顔だとか、名前ですら。 ありとあらゆるものに関してである。 決して、体型が分からないようにダボダボの服を着ましょうとか、名前はやめて優劣の無い記号にしましょう、とは言わない。 にも関わらず、弁当だけが平等の象徴かのように取り上げてくだらないルールを作るのは、 平等が大好きなあまり何が大事なのかということを考えなくなっている証拠じゃないかな。 この弁当ルールに関して、「育ち盛りの子供に米＋漬物のみ食事をさせるなんてありえない」と言い、 大きなおにぎりの中に大量の具材を入れてくれた、母親の方がよっぽどまっとうな考え方のように思える。

中学校時代にも無意味な校則が存在した。 今でも納得行かないものは、冬場の手袋とウィンドブレーカーに関する校則である。 なぜか、手袋を着用できるようになる日が決まっていた。 しかも、かなり寒くなってから。 自分たちは手がかじかむのを必死で耐えて自転車で登校していた。 校則にはなっていないが、ポケットに手を入れて自転車に乗ってはいけないと、熱心に生活指導する先生もいる。 なぜ手袋だめなの？

ウィンドブレーカーも同様に無意味なルールで縛られていた。 黒か紺のみ、という。 いやいや、安全第一って何時も言ってるやん！とおもった。 絶対派手な色の方が自転車乗るにも、道路脇を歩くにも安全でしょ。 未だに納得いかない。

さらには、こうした校則の改正は、生徒全員が校則を守ってからじゃないとだめだ！と言い出す、お馬鹿な先生も出てきた。 どこからそんな考えが生まれてくるのだろう。 問題のある法律が存在する時、国民が全員守るまで改正すべきでない、という考え方でもあるのだろうか？ とにかく生徒を平らに平らに、個性をなくし、順序付けをなくし、平等にするためなら、本末転倒もやむなし、というのが左翼先生の基本のようだった。

ソ連と社会党が大好き

授業中、露骨に自分の政治心情を生従に押しつける先生もいた。 とにかくその先生はソ連（今のロシア）と社会党（現在の社民党の母体）が大好きであった。 例えばこうである。

「ソ連は、全員が平等に扱われている。私達が見習わないといけませんね。」
「社会党も平等が理想で本当に良い政党だと、先生は思います。」
「皆さんはどう思ぃますか?」

とにかくこういう事を事あるごとに言う。 自分は、親から先生の言う事が必ずしも正しく無いので鵜呑みにしてはいけないと、釘を差されていたため全て信用した訳では無かった。 しかし、こういう聞き方をされれば小学生くらいであれば「私もそう思います」っ言えば、先生に誉められると思いそう言ってしまうものだ。 本当に信じるかどうかは別の話だが。

他にもこういう事を覚えている。 その先生が授業中にテレビをつけ、湾岸戦争のニュースを見せた。 当時は、自衛隊がPKO参加すべきか否かで議論された結果、結局参加出来ないまま、戦争が始まってしまった直後ことである。 公平に言って、生の二コースを見せて自分たちが勉強している内容との対応を示すのは非常に良い事だと思う。 しかし、先生はその後も露骨に偏った事を言う。

「自民党は、この戦争に参画しようとしてたのだから本当に悪い政党ですね。」
「PKOに反対する社会党は、戦争を阻止しようとしていますね。」
「皆さんはどう思ぃますか?」

遠足はヤマギシズム教育

遠足の行き先はヤマギシだった。正確には、「ヤマギシズム社会実顕地」（通称、ヤマギシの村）というらしい。 記憶が定かでないので、正確にこれかどうかは分からないが、バスに乗って連れて行かれ、 当地の担当者が「村」という呼び方をしていたのは覚えている。 ヤマギシをご存じない方のために簡単に説明する。ヤマギシは要するに、自由な国日本国の中にある共産主義組織である。 共産主義なので全ての人は平等に暮らすというのが理想だが、入会と同時に財産没収、 その後は「ヤマギシズム」という考え方を教育され、給料なし休日なしで農業に従事させられる。 政治団体というよりどちらかとカルト団体と見られることが多く、Amazonで「ヤマギシ会」と検索しても 「洗脳」、「カルト」などというネガティブな言葉が並ぶ。

公平に言うために、ヤマギシの良い印象も伝えておかねばなるまい。 ヤマギシが頻繁に農産物や加工品をトラックで売りに来る時期があった。 風の谷のナウシカの音楽をかけながらくるのはそのイメージが近いからだと思った。 ヤマギシが売るパンやシュークリームは当時のスーパーなどでは売ってなかったようなクオリティでめちゃくちゃ美味かったことを記憶している。今では結構あたりまとなったがシュークリームはでかくクリームが濃厚で美味かった。 天然酵母パンも当時のパン屋には売ってなかったような、手作り感あふれる味で、明らかに群を抜いて美味かった。

ヤマギシの村での体験は、前半は楽しかった。 ヤマギシでは、鶏も通常の養鶏所のようにただ餌を与えられて卵を生むだけではなく、 囲いの中を放し飼いにされていて、オスもメスも一緒くたに育てられていた。 卵はあちこちで生むので探しながらとらないといけない。 より自然な形で育てられているから良いのだという。

しかし、後半見せられたものには、今でも鮮明に記憶するほどの違和感を感じた。 １つの部屋に生徒全員が集められ、ヤマギシの村がどういうところかという説明を受けた。 まるでNHK教育テレビで出てくる体操のお兄さんのような若いお兄さんが登場して挨拶した。

「ハレハレー」

顔の両側に両手を広げて手のひらをこちら側に向け、カサゴのようにヒレをひらひらさせながら「ハレハレー」と声を放つのである。

「皆さんも一緒に！せぇ～のハレハレー」
「みなさ～～ん、ヤマギシでは喧嘩をする人は一人も居ません。誰もが笑って暮らしています。みんないつも笑顔です。」

はっきり言ってドラマ「トリック」に出てくるカルト団体そっくりである。 それ以外も色々説明してたけど、殆ど覚えていない。 この喧嘩をしないという点に子供ながら、嘘だと思ったことをはっきり覚えている。 根拠はないけど喧嘩しないやつがいる訳がない。 喧嘩の絶えない家族に生まれた自分が特殊なのかもしれないけど、とにかく率直にそう思った。

しかし、まぁ今思えば、自分の子供がこのような遠足に連れて行かれると知ったら反対するだろうな。 オオム真理教に遠足に行きますと言って怒らない親がいるだろうか？ 自分の思想は自由だけど子供に偏ったことだけを見せる事自体が洗脳と同じじゃないか。 しかし、当時はそれがまかり通ってた訳だから一人の特殊な左翼先生が居ただけではないということが、今になって分かる。

反戦教育のためなら暴力もやむなし

中学時代の修学旅行は当然、左翼先生が大好きの広島になった。 反戦教育ができるからである。 確か初日は広島の手前の町で観光し、翌日大きなホールで被爆者が講演する体験談を聞いた。

中学生にとっては全く面白くないけど、それは元気盛りの中学生、友達同士が同じ宿に泊まるというだけで楽しくて仕方ない。 当然、どこの学校でも同じ、前日の夜は興奮して全く眠れない。 予め枕投げは禁止と釘を刺されたけど、たった１つのゲームを奪われただけで、中学生の自由は奪えない。 先生が居ない時は寝ないでふざけまわり、先生が来ると寝たふりをする遊びで十分楽しめた。

ほぼ寝れないまま被爆者の体験談を聞く時間になった。 薄暗いホールで座らされ、特にしゃべり上手でもない被爆者の話を聞かされれば、どんな人でも眠たくなるのと思う。 自分は真ん中からやや後ろに座っていたのだけど、前の方の男の子も女の子も多くの首が傾いているのが見えた。

自分も、うとうととし始めたその時、周りの生徒がざわめき初めて目が覚めた。 生活指導などで生徒に暴力を振るうことで有名なあの先生が寝ている生徒めがけて突進してきた（ように見えたが実際は歩いてきた）。 体も手も丸々太った１９０ｃｍくらいの先生が生徒の頬を順番に バッチーーン、バッチーーン、バッチーーンと会場に音が鳴り響くくらいに連続ビンタしているのが見えた。 この先生、同じクラスのやんちゃな生徒の顔をめちゃくちゃに殴って本当に顔の形が変わるくらいボコボコにしたことがある。 実際の殴られた生徒の顔を見たが両目が腫れて戦った後のボクサー見たいになっていた。 担任の先生が大丈夫かと心配した程である。 そんな先生に連続ビンタされるのは嫌だから、その後も必死で寝ないことだけを考えて、長い講演を乗り切った。

中学生ながら後から考えてどう考えてもおかしいとおもった事を覚えている。 反戦教育は被爆者の方が可愛そうだと同情するための教育ではなくて、 暴力自体がダメなことなんだとと教えたかったに違いない。 にも関わらず、その教育をやってる目の前で教師が暴力を振るってどうするんだと。 しかも、それは警告射撃なしの即射殺のような連続ビンタである。 ヘタすれば死んでしまう生徒だっているかもしれない。 一体何を考えているのか、この先生はと思った。

こんな先生方が今義務教育組織の長になっている

これまで紹介した先生は一人ではない。 その中には、自分の得意分野を見つけてくれた、自分の人生にとって大事な大事な恩師も含まれる。 それでもやはり、偏った教育がなされていたように思う。 なぜ、義務教育の先生が特に左に偏っていたのか？ 日教組の影響が強いかもしれない。 平等が正義だと考える左翼先生にとっては、とにかく他の都合の悪い事実は無視して、 というより勉強不足で知りもせず、平等であればすべてが良いと信じ込んでいた。 しかし、少し勉強すれば、事実をより多く知れば、物事はそれほど単純ではないということが分かる。 心配なのは、ここで紹介した先生方は、当時30代後半から40代前半、今は丁度教育組織の長になっている。 こうした先生が世代交代により総入れ替えするのは、50年近くかかる。 教育システムを真剣に考えるということは、本当に大事だと痛感する。

参考

• 幸福会ヤマギシ会（公式ホームページ）
  http://www.koufukukai.com/
• Wikipedia「幸福会ヤマギシ会」
  http://ja.wikipedia.org/wiki/%E5%B9%B8%E7%A6%8F%E4%BC%9A%E3%83%A4%E3%83%9E%E3%82%AE%E3%82%B7%E4%BC%9A

Tracking関連論文のエレガントな変数定義例

PWP3D

Prisacariu, Victor and Reid, Ian: "PWP3D: Real-Time Segmentation and Tracking of 3D Objects," International Journal of Computer Vision, vol. 98, no. 3, pp. 1-20, 2012.

Let the image be denoted by $I$ and the image domain by $\Omega\subset\mathbb{R}^2$. An image pixel $x = [x,y]$ has a corresponding image value $I(x) = y$ (in our experiments an RGB value), a 3D point $X = [X,Y,Z]^T = RX_0 +T \in \mathbb{R}^3$ in the camera coordinate frame and a point $X_0 = [X_0,Y_0,Z_0]^T \in \mathbb{R}^3$ in the object coordinate frame. $R$ and $T$ are the rotation matrix and translation vector respectively (representing the unknown pose) and are parameterised by 7 parameters (4 for rotation and 3 for translation), denoted by $\lambda_i$. We use quaternions to represent rotation, so we use 4 parameters to capture the 3 degrees of freedom of the rotation.
We assume the camera calibration parameters (for each camera) to be known. Let $(f_u,f _v)$ be the focal distance expressed in horizontal and vertical pixels and $(u_o, v_o)$ the principal point of the camera. In the case of multiple cameras, the extrinsics were obtained using the Matlab Calibration Toolbox.

DTAM

Newcombe, Richard A. and Lovegrove, Steven J. and Davison, Andrew J.: "DTAM: Dense tracking and mapping in real-time", Proc. Int. Conf. on Computer Vision (ICCV2011), pp. 2320-2327, 2011.

We refer to the pose of a camera $c$ with respect to the world frame of reference $w$ as \[ T_{wc} = \left( \begin{array}{cc} R_{wc} & c_{w}\\ 0^T & 1 \end{array} \right) \] where $T_{wc} \in \mathbb{SE(3)}$ is the matrix describing point transfer between the camera's frame of reference and that of the world, such that $x_w = T_{wc}x_c$. $R_{wc} \in \mathbb{SO(3)}$ is the rotation matrix describing directional transfer, and $c_w$ is the location of the optic center of camera $c$ in the frame of reference $w$. Our camera has fixed and pre-calibrated intrinsic matrix $K$ and all images are pre-warped to remove radial distortion. We describe perspective projection of a 3D point $x_c = (x, y, z)^T$ including dehomogenisation by $\pi (x_c) = (x/z, y/z)^T$. Our dense model is composed of overlapping keyframes. Illustrated in Figure 1, a keyframe $r$ with world-camera frame transform $T_{rw}$, contains an inverse depth map $x_{r}: \Omega \rightarrow \mathbb{R}$ and RGB reference image $I_r : \Omega \rightarrow \mathbb{R}^3$ where $\Omega \subset \mathbb{R}^2$ is the image domain. For a pixel $u := (u, v)^T \in \Omega$, we can back-project an inverse depth value $d = \xi (u)$ to a 3D point $x = \pi^{-1}(u, d)$ where $\pi^{-1}(u, d) = \frac{1}{d} K^{-1}\dot u$. The dot notation is used to define the homogeneous vector $\dot u := (u, v, 1)^T$.

KinectFusion

Newcombe, Richard A. and Davison, Andrew J. and Izadi, S. and Kohli, P. and Hilliges, Otmar and Shotton, J. and Molyneaux, David and Hodges, Steve and Kim, David and Fitzgibbon, A.: "KinectFusion: Real-time dense surface mapping and tracking," Proc. 10th IEEE Int. Symp. on Mixed and Augmented Reality (ISMAR2011), pp. 127-136, 2011.

We represent the live 6DOF camera pose estimated for a frame at time $k$ by a rigid body transformation matrix: \[ T_{g,k} = \left[ \begin{array} R_{g,k} & t_{g,k}\\ 0 & 1 \end{array} \right] \in \mathbb{SE}^3 \] where the Euclidean group $\mathbb{SE}^3 := \{R, t | R \in \mathbb{SO}^3, t \in \mathbb{R}^3g\}$. This maps the camera coordinate frame at time $k$ into the global frame $g$, such that a point $p_k \in \mathbb{R}^3$ in the camera frame is transferred into the global co-ordinate frame via $p_g = T_{g,k} p_k$. We will also use a single constant camera calibration matrix $K$ that transforms points on the sensor plane into image pixels. The function $q = \pi(p)$ performs perspective projection of $p \in \mathbb{R}^3 = (x,y,z)^T$ including dehomogenisation to obtain $q \in \mathbb{R}^2 = (x/z,y/z)$. We will also use a dot notation to denote homogeneous vectors $\dot u := (u^T|1)^T$

Gnuplotによる区分関数の最小二乗法

不等号で定義域を分けて異なる関数を組み合わせた作った関数を区分関数（Piecewise Function）という。 例えば絶対値を用いた関数$f(x)=|x|$も$x=0$を境に$x>=0$における$f(x)=x$と、$x<0$における$f(x)=-x$の組み合わせと考えられる。 誤差を含むデータ点から区分関数のパラメータを求める最適化はGnuplotのfitコマンドで意外と簡単にできる。 面白いのは区分関数の区切り位置（2つの関数の境目）自体もパラメータとすることでその位置も求められること。

V字関数のフィッティング

データの作成が容易なので、次の簡単な区分関数のフィッティング方法を紹介する。 正確な名前は分からないがここではV字関数と仮名を付けておく。 この関数どういう形かというと、$a<0,\ b>0$のとき$x$軸上の$x=c$で$y=0$となり、折り返すまさにV字の関数となる。

\[ f(x)= \left\{ \begin{array}{ll} a(x-c) & (x < c)\\ b(x-c) & (x \geq c)\\ \end{array} \right. \]

上記関数に対して誤差を含むデータ点が与えられたとき、 Gnuplotのfitコマンドでパラメータ$a,b,c$を同時に非線形最適化して求めてみる。 使用したGnuplotのコマンドとデータは以下のとおり。

Gnuplotコマンド

reset

#関数の定義
f(x,a,b,c)= (x<c) ? a*(x-c) : b*(x-c)
 
#初期値の設定
a=0
b=0
c=0

#関数のフィッティング
fit f(x,a,b,c) 'sample.txt' via a,b,c
 
#表示設定
set xr[-10:10]
set samples 100
set terminal wxt

#関数とデータの描画
plot f(x,a,b,c) 
replot 'sample.txt' lc 3

データ "sample.txt"

-4 11
-3 7.5
-2 6
-1 3
-0 2
1 0
2 5
3 12
4 16
5 20
6 27

このデータは、$a=-2, b=5, c=1$の時の関数の値に対してわずかな誤差を人手により与えている。 最適化の初期値は意図的に正解値から遠ざけて$a=0, b=0, c=0$としたが、おおよそ上の正解値に近づいた結果$a=-2.0544, b=5.21874, c=0.936731$が得られている。

データとフィッティングされた関数

フィッティング結果

gnuplot> load 'fitting_sample.txt'
 
 Iteration 0
 WSSR        : 1780.25           delta(WSSR)/WSSR   : 0
 delta(WSSR) : 0                 limit for stopping : 1e-005
 lambda   : 1.91485
 
initial set of free parameter values
 
a               = 1e-030
b               = 1e-030
c               = 1e-030
/
　　　：
　　（省略）
　　　：
 
 Iteration 6
 WSSR        : 5.75979           delta(WSSR)/WSSR   : -9.82704e-009
 delta(WSSR) : -5.66017e-008     limit for stopping : 1e-005
 lambda   : 1.91485e-005
 
resultant parameter values
 
a               = -2.0544
b               = 5.21874
c               = 0.936731
 
After 6 iterations the fit converged.
final sum of squares of residuals : 5.75979
rel. change during last iteration : -9.82704e-009
 
degrees of freedom    (FIT_NDF)                        : 8
rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 0.848513
variance of residuals (reduced chisquare) = WSSR/ndf   : 0.719974
 
Final set of parameters            Asymptotic Standard Error
=======================            ==========================
 
a               = -2.0544          +/- 0.1336       (6.501%)
b               = 5.21874          +/- 0.1976       (3.787%)
c               = 0.936731         +/- 0.1152       (12.3%)
 
correlation matrix of the fit parameters:
 
               a      b      c      
a               1.000 
b               0.403  1.000 
c               0.490  0.822  1.000・

参考

• Walfram MathWorld "Piecewise Function"
  http://mathworld.wolfram.com/PiecewiseFunction.html
• 米澤進吾先生のページ
  http://www.ss.scphys.kyoto-u.ac.jp/person/yonezawa/contents/program/gnuplot/wxt.html

cv::gpu::flipの入出力に同じ変数は使えない

OpenCVのGpuMatは非常に便利だけど、GPGPUの仕組みを全く理解せずにBlack Boxとして使うと時々問題になる。 例えばcv::gpu::flip()関数。 CPUのみを使う関数cv::flipでは入出力に同じ変数を用いても問題にならなかったが、 cv::gpu::flip()の場合は正常にflipされない（OpenCV2.6.4で確認）。 次の有名な標準画像を入力として試してみた。

#include <opencv2/opencv.hpp>
#include <opencv2/gpu/gpu.hpp>
 
void main(int argc, char *argv[])
{
    cv::Mat image = cv::imread("lena_std.png");
    cv::gpu::GpuMat image_gpu;

    //画像をCPUからGPUに送る
    image_gpu.upload(image);
 
    //CPU上で画像の上下反転、保存
    cv::flip(image, image, 0);
    cv::imwrite("fliped_lena_cv.png", image);
 
    //GPU上で画像の上下反転（srcとdstを同じにしている）
    cv::gpu::flip(image_gpu, image_gpu, 0);
 
    //画像をGPUからCPU側に戻しファイルに保存
    image_gpu.download(image);
    cv::imwrite("image_gpu.png", image);
}

上記コードで出力される画像は次のとおり。 左はcv::flip()、右はcv::gpu::flip()の結果。 GPUの方は、正常に画像が生成されていないことが分かる。

次のコードは上記コードを少し修正して、GpuMatを2つ用意して、 cv::gpu::flip()の入出力の変数を変えてみた。

#include <opencv2/opencv.hpp>
#include <opencv2/gpu/gpu.hpp>

void main(int argc, char *argv[])
{
    cv::Mat image = cv::imread("lena_std.png");
    cv::gpu::GpuMat image_gpu;
    cv::gpu::GpuMat image_gpu2;

    //画像をCPUからGPUに送る
    image_gpu.upload(image);
 
    //CPU上で画像の上下反転、保存
    cv::flip(image, image, 0);
    cv::imwrite("fliped_lena_cv.png", image);
 
    //GPU上で画像の上下反転（srcとdstを変えた）
    cv::gpu::flip(image_gpu, image_gpu2, 0);
 
    //画像をGPUからCPU側に戻しファイルに保存
    image_gpu2.download(image);
    cv::imwrite("image_gpu22.png", image);
}

左はcv::flip()、右はcv::gpu::flip()の結果。 今回はCPU、GPU共に同じ結果が得られた。

OpenGLで超簡単なオフスクリーンレンダリング

やりたい事

やりたい事はタイトルの通り。 GLUTを使うとOpenGLを簡単に扱える反面、glutMainLoop()に全部制御を持って行かれてしまう。 FBO(frame buffer object)のような難しいものは触りたくない。 ウィンドウを表示せずにOpenGLで描いた画像を取得する方法を探しても意外と見つからなかったのでここにメモしておく。 要は、glutCreateWindow()関数を呼んでウィンドウを作る準備をしておけば、 バッファは利用可能な状態になっているようなので、 実際にウィンドウが表示されなくても描画したり、描画したものを画像として取得したりできるということ。 この方法だと複数のウィンドウも可能だけど、１つのウィンドウを表示してもう一つのウィンドウは表示しないということは出来ないかも。

サンプルコード

#include <gl/glut.h>
#include <opencv2/opencv.hpp>
 
void main(int argc, char *argv[])
{
    int w = 256;
    int h = 128;
    cv::Mat_<uchar> image(h, w);
    glutInit(&argc, argv);
    glutInitWindowSize(w, h);
    glutCreateWindow(argv[0]);
 
    //描画
    glClearColor(0.1, 0.1, 0.1, 1.0);
    glClear(GL_COLOR_BUFFER_BIT);
    glutSolidCube(0.5);
    glFlush();
 
    glReadPixels(0, 0, w, h, GL_RED, GL_UNSIGNED_BYTE, image.data);
 
    cv::imwrite("test.png", image);
}

出力画像