風景から歩行者を消す手軽な方法


固定したカメラで撮った動画で、画素ごとに時間について平均を取れば、(適当な速度で)動くものを消せます。Mathematicaだとこんな感じです。(参照:フリーソフトウェアを使う方法

Export["result.jpg",
 Image[Mean[Map[ImageData,
    Import["movie.mov", "ImageList"]]]]]

おまけ:フレームの平均を計算していく過程(最初の5秒を30秒で)
詳細:風景から歩行者が消えていく様子(リアルタイム版)

追記:画質的には平均ではなく中央値や最頻値を使った方がいいかもしれませんが、「手軽」ではなくなります。「平均でもできるんだ」という「手軽」さの実例だと理解していただければと思います。

中央値:MeanMedianに置き換えるだけで試せますが、計算時間・消費メモリともに増大します。平均なら約90秒で終わるこの動画(1280x720x372フレーム)の処理に約450秒かかります。消費メモリは2倍くらいになるようです(Core i7 4700MQ、メモリ16GB、Windows 7 64bit、Mathematica 9.0.1)。

最頻値:中央値と同様、計算負荷が高くなります。数値の揺らぎも心配です。Mathematicaの最頻値(Commonest[])は戻り値がリストなので、コードの書き換えも面倒です。

カメラからの入力を「手軽」にリアルタイム処理する場合にも、やはり平均を使うのがいいでしょう。

縦と横にしか動けない世界で(1994年東京大学入学試験理系数学第6問)


1994年の東京大学の入学試験、理系数学第6問は次のようなものでした。

平面上の2点P, Qに対し、PとQをx軸またはy軸に平行な線分からなる折れ線で結ぶときの経路の長さの最小値をd(P, Q)で表す。

(1) 原点O(0, 0)と点A(1, 1)に対し、
d(O, P)=d(P, A)を満たす点P(x, y)の範囲をxy平面上に図示せよ。

(2) 実数a>=0に対し、点Q(a, a^2+1)を考える。
次の条件(*)を満足する点P(x, y)の範囲をxy平面上に図示せよ。
(*) 原点O(0, 0)に対し、d(O, P)=d(P, Q)となるようなa>=0が存在する。

この問題は、d(O, P)=Abs[x]+Abs[y]であることがわかれば解けます。d(O, P)=Sqrt[x^2+y^2]ではありません。

Mathematicaで試します。

(*1*)
d[p_, q_] := Total[Abs[p - q]]
expr = d[{0, 0}, {x, y}] == d[{x, y}, {1, 1}];
cond = Reduce[expr, {x, y}, Reals];
reg = ImplicitRegion[cond, {x, y}];
RegionPlot[reg, PlotRange -> {{-2, 2}, {-2, 2}}]

(*2*)
expr = Exists[a, a >= 0,
   d[{0, 0}, {x, y}] == d[{x, y}, {a, a^2 + 1}]];
cond = Reduce[expr, {x, y}, Reals];
reg = ImplicitRegion[cond, {x, y}];
RegionPlot[reg, PlotRange -> {{-2, 2}, {-2, 2}}]

描画領域の境界にも実線が引かれている、という問題がありますが、とりあえずはこれでいいでしょう。

Mathematica 10.3.1と11.2, 11.3には,直線部分が描かれないというバグがあります(紛らわしいことに10.4, 11.0.1は大丈夫)。11.2以降なら,RegionImageを使うといいかもしれません。

RegionImage[reg, PlotRange -> {{-2, 2}, {-2, 2}}]

関連:数学まちがい大全集

数学ガールのミルカさんが「こうなる」とひと言発する間にやっていること(ウラムの螺旋)


4797374152結城浩『数学ガールの秘密ノート 整数で遊ぼう』(SBクリエイティブ, 2013)(索引あり・参考文献リストあり)にウラムの螺旋が出てきます。

100くらいまで描いたところで、

「あ、あたしっ、この先にも興味があります。ぐるぐるぐるぐるぐるぐるぐる……とずっと続けたらどんな図形が表れるんでしょう」

と言うテトラちゃんに応えて、

「こうなる」

と言ってミルカさんがかなり大きなウラムの螺旋を見せています。(p.70)

「ウラムの螺旋」で検索すればぱっと出せそうな気もしますが(インタラクティブなものなんかも)、本書で描かれているのは1から始まるよくある螺旋ではなく、0から始まる螺旋なので、そうでもないかもしれません。そもそもミルカさんは、「こうなる」といいながらウェブ検索をするキャラではありません。

ですから彼女は、「こうなる」とひと言発するわずかの時間でウラムの螺旋を描くスクリプトを書いているんだと思うのです。(彼女はノートPCを持ち歩いている設定でしたっけ?)

そんなことができるのかと思って最初に思いつくのはこんなスクリプトです(Mathematica)。

spiral[size_, start_] := With[{center = Ceiling[(2 size - 1)/2]},
  Module[{
    m = Table[start, {2 size - 1}, {2 size - 1}],
    i = start},
   Do[
    Do[m[[r, center + k]] = ++i, {r, center + k - 1, center - k, -1}];
    Do[m[[center - k, c]] = ++i, {c, center + k - 1, center - k, -1}];
    Do[m[[r, center - k]] = ++i, {r, center - k + 1, center + k, 1}];
    Do[m[[center + k, c]] = ++i, {c, center - k + 1, center + k, 1}],
    {k, 1, size - 1}];
   ArrayPlot[PrimeQ[m] /. {True -> 1, False -> 0}]]]

これで「spiral[4, 0]」などとすればウラムの螺旋を描けますが、会話のテンポはかなり遅くなるでしょう。

会話のテンポを自然なものに保つためには、もっとコンパクトなスクリプトを書けなければなりませんが、すぐには思いつきません。(かつて書いたスクリプトがノートPCに保存してあった?)

というわけで、彼女は優秀な高校生だという物語の設定は、そのとおりだと思いました。

会話形式で書かれている数学ガールには、会話のペースに合わせてできるものだろうかと考えながら読む楽しみがあります。

せっかくなので、サイズを大きくしたときの変化を見てみましょう。

Animate[spiral[size, 1], {size, 1, 100, 1}]

サイズを大きくしたときの変化

中心の数を変化させたときの様子も見てみましょう。

Animate[spiral[50, start], {start, 0, 100, 1}]

中心の数を0から100まで変化させた様子

Wolfram CDF Playerがインストールされていれば、インタラクティブに調べられます。


一様乱数から正規乱数を作る方法?


ソフトウェア関係の雑誌を読んでいたら、次のようなことが書かれていました。

「0以上1未満の一様乱数を12回足して6を引くと正規乱数になる」というテクニックを覚えておけば、一様乱数しか提供されていない開発環境で大いに役立つはずだ

これは中心極限定理のいい例ですね。実際に乱数を大量(下の絵では100万個)に発生させてそのヒストグラムを描いてみると、正規分布になっているように見えます(太線は正規分布)。これで十分な場合も多いでしょう。

しかし、正規性の検定をすると、正規分布であるという仮説を棄却したくなるようなp値が得られます。Mathematicaで試します(参考:DistributionFitTest)。

In[1]:= randomGaussian[] := Sum[RandomReal[], {12}] - 6

In[2]:= data = Table[randomGaussian[], {1000000}];

In[3]:= DistributionFitTest[data, NormalDistribution[0, 1]]

                  -8
Out[3]= 8.33717 10

この方法で生成される乱数の分布は、足す数(上の例では12)を増やせば正規分布に収束しますが、正規分布ではありません。ですから、「正規乱数になる」という表現はいけません。(例えば「randomGaussian[] := (Sum[RandomReal[], {24}] - 12)/Sqrt[2]」ならもっと正規分布に近づきますが、依然として正規分布ではありません。)

一様乱数しか提供されていない環境では、ボックス=ミュラー法などを使うのがいいと思います。Excelでは「=NORM.INV(RAND(),0,1)」が簡単です。

4874085601Pressほか『ニューメリカルレシピ・イン・シー』(技術評論社, 1993)にこんな記述がありました。

悪いrand()のため結果が疑わしいすべての科学論文を図書館の棚から消せば、各棚にこぶしが入るほどの隙間ができるであろう。(p.203)

「悪いrand()」ではなく「悪い乱数」を対象にすると、隙間はどのくらいになるのでしょう。

Mathematica 10への期待(とその後)


昨年ついに無料になったMathematicaですが、思ったほど盛り上がってはいないようです。盛る上がらない原因の一つに、性能の悪さがあるのは間違いないでしょう。特にGUIはかなり重く、Manipulateなどでインタラクティブに遊ぶのは無理です。(無料のMathematicaが動くのは、700MHzのCPU、最大でも512MBしか主記憶のないRaspberry Pi上。)

実際にベンチマークを取ってみると、Windows上のMathematica 9.01(CPUはCore i7 4700MQ、主記憶は16GB)で10秒で終わるMathematicaMark9が、Raspberry Pi上では約3300秒かかりました(800MHzにオーバークロックして約3000秒)。Core i7で動かすエミュレータ上でもおそらく同じ程度でしょう。

しかし不思議なことに、Raspberry Pi上のMathematicaのほうが早く解ける問題があります。

たとえば、A^2 + B^2 == 1かつP^2 + Q^2 == 1という条件で、以下のsの最小値を単純にMinValueで「解析的に」求めようとすると、上述のWindowsマシンではなんと約20000秒かかるのに対して、Raspberry Piでは約2200秒で終わります。追記:Raspberry Pi 2では約800秒です。

s = 1/24 (3 Sqrt[3] Abs[A P] + Sqrt[3] Abs[A P + 2 Sqrt[2] Q] + 
     Abs[Sqrt[3] A P - 3 Sqrt[2] B P - Sqrt[6] Q] + 
     Abs[Sqrt[3] A P + 3 Sqrt[2] B P - Sqrt[6] Q]);

MinValue[{s, And[A^2 + B^2 == 1, P^2 + Q^2 == 1]}, {A, B, P, Q}]

(これは比較のために作った問題ではなく、実際に私が遭遇した問題です。実は、単純にMinValueを使うというのが間違っているのですが、本来Mathematicaはそういうことを気にせず使うソフトウェアのはずです。)

Raspberry Piの方が速い原因は、Windows版Mathematicaのバージョンが9.01なのに対して、Raspberry Pi版Mathematicaのバージョンが10であり、バージョン9.01と10の間で、このあたりのアルゴリズムが改良されたためだと思われます。

非力なRaspberry Piでの話なので、比較的速いCPUがあれば、エミュレートされたRaspberry Piでも同程度の時間でできるでしょうし(Windowsの場合Ubuntuの場合)、ユーザーモードエミュレーションならこの半分くらいの時間でできるでしょう(ライセンス違反なので想像しているだけです)。

というわけで、Windows版Mathematica 10のリリースを待っています。

追記:Windows版Mathematica 10がリリースされたので試したところ54秒でした(CPUはCore i7 4700MQ、主記憶は16GB)。Raspberry Piは約2200秒なので、遅すぎです