2016年2月24日水曜日

水玉コラ作成ソフトを公開した話

Bloggerで過去の記事に追記をすると最新の記事の扱いになる現象が勃発している今日この頃です.先日「水玉コラ作成ソフト」をVectorさんで公開しました


一応こんな感じに動作する.
これを使うことで慣れれば2,3分で水玉コラを作ることができるようになる.やったぜ.

一度フリーソフトを公開してみたかったので夢というか目標が達成できてよかった.Vectorさんもフリーソフトの審査を人の目が見てくれているようでとても良い対応だったのでよかった.

ただ現状ダウンロードして実行すると次のような画面が出る……(´・ω・`)
しばらくの辛抱らしい.

Windows SmartScreenは認識されないアプリの起動を停止しました……と.一応弁明させていただきますと僕がこの画面を作ってパスワードを盗み取ろうとしているわけではないです.

この画面,沢山ダウンロードされると消えるらしいのだけど俺だったらこんなソフト使わないな(笑).こういう詐欺もあるだろうし,気をつけねば.

フーリエ級数の可視化プログラムの実装[Processing]

もはや2週間前のネタを真似ただけで新規性に乏しいのですが……(´・ω・`)
一応実装にあたって元ネタを書いておかないと申し訳ないのでツイートを貼らせていただきます……



という経緯です.

今回の実装はProcessingによるものなので上記ツイートのような数式で描くことができませんが,プログラミング言語というだけあって自由度は高いものとなっています.

円の数とか任意に決められますし.関数も自由に決められますし.

追記:鰺坂もっちょ先生がトゥギャってくれました→世はまさに大フーリエ時代 - Togetterまとめ

・ソースコード


例として矩形波を近似しています.関数funcの中身を[0,1)区間で変えてやればその形になりますが,無限大に発散させたりするとダメです.
あと数字が小さすぎる(1.0とか)だと小さすぎて見えませんので気をつけてください.


2016年2月17日水曜日

画像を一筆書きで描くプログラムを実装した[Python,Processing]

1.概要

この動画の最後のやつの,ほぼ完璧なトレースができます.
 

とりあえず順序としては
1.画像を読み込む
2.エッジ部分の座標を取り出す
3.巡回セールスマン問題として一筆書き経路を探索
4.フーリエ級数展開
5.画像を描画
という手順です.

今回サンプルとして画像はこいつを使います.勝手に使ってください.
Processingしかできない人や,Pythonの実行環境ないひとは一筆書きの経路を探索した結果の座標ファイルを置いておきます.こいつをProcessingのソースコードと同一フォルダに入れておけば一応はできます.ちなみにこんな感じ.


プログラムが長過ぎるので注意してください.もっとスマートにかければ良いのだけれど,無駄にクラスを作成して再利用しやすい形にした……(´・3・`)ドウセツカワナイノニ

2.プログラム

まず上記の手順の1-3を行うソースコードです.
蟻コロニー最適化のプログラムなのでこの前書いたやつの書き直しみたいな感じです.
実際はopt2を使って更に最適化してますが,それのコードはネットのやつを丸パクリしたやつなので流石に上げられません(参考URL先に借りたソースコードがあります).


次にProcessingで上記の手順の4,5の描画を行うコード

初音ミクを三角関数で描いてみた

1.動画紹介

フーリエ級数の視覚化を見て動画を作った.



まぁ二番煎じなんですけど.その辺の詳細はこちらで.

2.手法とか

まず序盤のやつはこいつ.Processingで作ったフーリエ級数の視覚化プログラム.こいつでいろんな関数を近似している.ここの先で先駆者たちのツイートが見られます.
フーリエ級数の可視化プログラムの実装[Processing]


次に出てくる二次元平面での三角関数的な奴ってのはこれ.実はフーリエ変換ではなく離散コサイン変換の二次元の基底関数なので嘘といえば嘘だけど,まぁそんなに違いはないし説明が面倒なので省略した.
離散コサイン変換(DCT)をPythonで実装した[Python]


その次の巡回セールスマン問題についてはこいつで解いた.基本的に蟻コロニー最適化を用いているけど,2000個の点で蟻コロニー最適化をすると遅すぎるので,パラメータを変えて貪欲法みたいな方法で経路探索して,その下のサイトで紹介されているopt2法のプログラムを用いてる,

・巡回セールスマン問題を蟻コロニー最適化と遺伝的アルゴリズムで解いてみる[Python]
Algorithms with Python / 巡回セールスマン問題


そんで最後のやつはフーリエ級数の可視化プログラムを使って,巡回セールスマン問題でもとめた一筆書きのルートを描いてる.
・画像を一筆書きで描くプログラムを実装した[Python,Processing]

結構大掛かりなプログラムだった割に動画を見てみると正直あまり面白くない……(´・ω・`)
アイデアおもついて実装したら面白くないのよくあるから仕方ないか…….

2016年2月14日日曜日

巡回セールスマン問題を蟻コロニー最適化と遺伝的アルゴリズムで解いてみる[Python]

巡回セールスマン問題(Traveling Salesman Problem,TSP)は有名だけれどなにげに解いたことがなかったのでちょいと実装してみました.
手法としてはいろいろあると思うのですが,今回は遺伝的アルゴリズム(GA)と蟻コロニー最適化(ACO)の2つを実装・比較しました.

1.巡回セールスマン問題

よく出される例としてはゴミ収集車の経路でしょうか.ゴミ収集所から出発して,100軒の家を最短で全て回ってゴミ収集所にまた戻ってくるまでの経路を求める話です.

まぁ僕がここで云々説明せずともネット上には沢山関連するサイトが有りますからそちらを参照してください.とりあえずWikipediaと検索したらトップに出てきたサイトのURLを貼っときます.

真ん中のサイトは手法について非常に詳しく解説されている上,Pythonで実装されたソースコードがありますのでそちらもごらんになると良いかと思います.

一番下のサイトは巡回セールスマンについての説明が詳しく,フリーソフトの紹介もされています.


この問題はまだ多項式時間で最適解を解くことができるアルゴリズムが無く,うちの大学でも研究してる研究室がたしかあった気がします.

2.蟻コロニー最適化による実装

蟻コロニー最適化に関する数式は今回全てWikipediaに載っていたものを利用しました(蟻コロニー最適化 - Wikipedia).
一応『心はプログラムできるか』という本を読んだ時に出てきてまして非常にわかりやすい解説がなされていましたので読んでみると良いかと思います.

以下ソースコード
とりあえず評価のちゃんとした結果は遺伝的アルゴリズムの実装の後に書きます.

3.遺伝的アルゴリズムによる実装

ぶっちゃけ蟻コロニー最適化よりも精度低いです.30個点でさえも蟻コロニー最適化に負けています.1000個ぐらいの点になるともう最適解求めるなんて夢のまた夢でした.
原因としては恐らく突然変異がランダムすぎて,ほとんど改悪になるばかりという結果でしょうか.遺伝子は巡回経路の順番を表した変数(5つの点なら[2,4,3,1,0]という配列)です.突然変異で中身を入れ替えてます.
※追記:先日図書館にてほかの人の実装例を見ましたがこれとは異なる手法を用いていました.それを踏まえるとこのプログラムは一般的な実装例とは異なると思われますので注意してください.

ソースコード

4.評価

一応同じデータを用いてやった結果です.30個の[0,1)のランダムな二次元の点について行いました.
パラメータはデフォルトです.まずは遺伝的アルゴリズム.
直線が交差してますね.コストは6.056765でした.続いて蟻コロニー最適化.
直線が交差すること無く,コストも5.701502です.しかも計算時間が遺伝的アルゴリズムよりも遥かに実行時間が早いです.

とりあえず今回の実装では蟻コロニー最適化のほうが早かったですが,遺伝的アルゴリズムも突然変異をもう少し工夫してやればもっと良い結果が出るような気がします.いつか焼きなまし法とかやってみたいです.

2016年2月10日水曜日

離散コサイン変換(DCT)をPythonで実装した[Python]

ダブルミーニングでテストが終わりました.春休みに入ったので今年の目標の「数学を勉強する」を早速実行してみます.将来勉強するはめになりそうなものから学ぼうと思うので,まずは主軸変換から勉強します.

フーリエ変換はnumpyやscipyで実装されているので,されていないであろう離散コサイン変換(DCT)を今日は実装してみました……と思っていたらscipy.fftpackにdctの実装ありました(公式サイト).こちらはWikipediaで紹介されている基底関数DCT-I,DCT-II,DCT-IIIが実装されているようです.

1.離散コサイン変換(DCT)の概要

離散コサイン変換(Discrete Cosine Transform)はWikipediaさんを見る限りではいろんな基底関数があるようですが,今手元にある本『画像処理とパターン認識入門』では離散コサイン変換の基底関数を次のように定義しています.Wikipediaで言うところのDCT-IIですね.
\[ \phi_k[i] = \begin{cases} \frac{1}{\sqrt{N}} \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,(k=0) \\ \sqrt{\frac{2}{N}} \cos {\frac{(2i+1)k\pi}{2N}} (k=1,2,...,N-1) \end{cases} \]

この基底関数ですがこんな感じのグラフになります.N=64,横軸がiに相当し上からk=0,1,2,...,9の時の基底関数のグラフです.この定義だとほぼ普通の三角関数ですね.
離散化された信号をf_iとすると,離散コサイン変換,離散コサイン逆変換の式は以下のとおりです.
\[ C[k] = \sum_{i=0}^{N-1}f_{i}\phi_k[i] \,\,(k=0,1,...,N-1)\\ f_{i}=\sum_{k=0}^{N-1}C[k]\phi_{k}[i] \,\,(i=0,1,...,N-1) \]
とりあえずこいつを実装します.

2.離散コサイン変換(1次元)の実装例

まぁこんな感じです.関数phiは基底関数です.本で紹介されていた数式以外にWikipediaにあったDCT-IVもコメントアウトしてありますが一応実装してあります.
ベンチマークとして乱数で作った配列データを試しに変換・逆変換してみます.

実行結果

まぁ問題なくデータが復元されてますね.乱数なので高周波数成分があるので圧縮は無理ですが.

3.離散コサイン変換(2次元)の概要

さて手元の本のタイトルには「画像処理」の文字が書かれています.すなわち二次元データに対する解説も書かれています(というかむしろそれが本題).なのでpythonで二次元データに対しても実装してみます.

本によれば二次元のDCT係数F_k,lは次のように表されるとのこと.また画像の濃度値f_i,jはDCT係数を用いてこのように逆変換できるとのこと.
\[ F_{k,l}=\sum_{j=0}^{N-1}\sum_{i=0}^{N-1}f_{i,j}\phi_k[i]\phi_l[j] \\ f_{i,j}=\sum_{l=0}^{N-1}\sum_{k=0}^{N-1}F_{k,l}\phi_k[i]\phi_{l}[j] \]
ここでこの基底ベクトル(基底画像)がどんな形をしているか見てみます.
この画像,離散コサイン変換(DCT)で検索すると出てくるやつですね.これらの画像は10x10ピクセルの画像を100個並べたものですが,二次元離散コサイン変換ではこの100個の画像をうまいこと足し合わせることで様々な画像を表現するというのがテンションが上ります.(つまりこの画像の中にあんな画像,こんな画像まで,ありとあらゆる画像が含まれていると言っても過言ではない……?)

4.離散コサイン変換(2次元)の実装例

では実装です.scipyにも2次元の離散コサイン変換は実装されてないのでちょっとテンションが上ります.まぁ2次元のフーリエ変換は実装されてますが.
サンプルとして適当な画像を変換・逆変換してみます.
実行結果です.
画像がちゃんと元の状態に戻っているのがわかります.ちなみにDCT係数(周波数成分)はこんな感じでした.
概ね右下に行くほど色が薄くなっています.ここから低周波成分(左上)が多くを占めていることがわかります.
高周波成分を消していったものがこちらです.
DCT係数cについてc[i:,i:] = 0とした時に復元された画像です.i=3でだいたい認識できそうですね.結構データを圧縮できてるのではないでしょうか.まぁまぁうまく行って良かったです.

2016年2月5日金曜日

相関係数を用いてmnistのラベルを認識してみる[Python]

線形代数はわかるけど,ディープラーニングとかわかる気がしねぇ(゜∀。)ワヒャヒャヒャヒャヒャヒャ
という人向けです.
各ラベルごとの平均画像との相関関係からmnistのラベルを予想します.

1.使用する画像

MNIST は 0 から 9 の手書き数字が書かれた,70,000 点の白黒画像です.こんな感じ.
画像は28×28の784ピクセル.今回は学習データを評価用のデータに含めちゃダメなので,この画像データ70,000点のうち,60,000点を学習用のデータとし,10,000点を評価用のデータとします。
とは言うものの学習という学習はしないのであまり関係がないけれど.

2.相関係数の求め方

まぁ今回は一応相関係数を求める部分も実装したのでTeXの練習がてら書いときます.普通はnumpy.corrcoef関数とか使えばいいと思います.
次のようなベクトルx,yがあったとします.
このとき,x,yの相関係数は次のようになります.

ただし,xバー,yバーは各要素の平均です.
同じことですが,これはx,yの平均ベクトルを次のように定義してやるともう少し見通しが良くなります.
とすると,相関係数rは
と表すことができ,x,yベクトルから要素の平均値を引いたベクトルの方向余弦(cosθ)になるわけです.(割とこれが言いたいだけの記事だったりする.線形代数の先生がドヤ顔して話してた)

今回作ったプログラムではxは調べたい画像のベクトル,yは0,1,2,3...9の平均画像に当たります.最も似ている(相関がある)ものを探していきます.

3.ソースコード


4.結果

こんな感じ.
[[ 878    0    3    6    4   40   36    0   14    2]
 [   0 1095   16    4    1    6    2    3   33    1]
 [  26   37  732   44   23    1   43   24   45    4]
 [   8   14   35  809    2   45    6   17   73   25]
 [   4   14    5    0  818    1   24    4   17  117]
 [  29   41   13  127   21  554   24    9   36   23]
 [  21   23   22    1    9   25  862    0    4    0]
 [  17   27   15    0   17    1    0  823   18   68]
 [   7   44   14   85    8   31   10    8  728   38]
 [  19   19    9   12   66    4    1   53   33  820]]
             precision    recall  f1-score   support

        0.0       0.87      0.89      0.88       983
        1.0       0.83      0.94      0.88      1161
        2.0       0.85      0.75      0.79       979
        3.0       0.74      0.78      0.76      1034
        4.0       0.84      0.81      0.83      1004
        5.0       0.78      0.63      0.70       877
        6.0       0.86      0.89      0.87       967
        7.0       0.87      0.83      0.85       986
        8.0       0.73      0.75      0.74       973
        9.0       0.75      0.79      0.77      1036

avg / total       0.81      0.81      0.81     10000
統計はよく知らない民なのでprecision,recallでググるとwikipediaの画像が非常にわかりやすいので見ると良いかと思います.まぁ普段生活する中で言う「的中率」というのはrecallだと思う….

幸いprecision,recall,f1-scoreも81%なので,とりあえず結果は81%ってことですね.mnistの公式サイトを見ると圧倒的最下位ですね.現状1位はConvolutional neural networkの99.77%ッスか.もうこの域まで来ると人よりも精度がいいらしいのでもう関係ないような気がしますがね.

5.感想とか

こうしてみるとやっぱりニューラルネットすげーってなるわけですが,もう少しニューラルネットを使わずに精度を上げられないかと思われるわけですが,ミスした画像を見ればどうすればいいかはある程度わかります.
一番左は平均画像.周りが判定ミスをした画像たちです.数字が書かれている位置が違っていたり,ある程度傾いていたりすると平均画像との相関係数だけだと限界があるわけですね.
もし更に精度を上げるには傾きや位置によらない特徴量を抽出してやる必要がありそうです.その特徴量ベクトルに対して相関係数を用いるのもよし,ユークリッド距離を用いるのもよし,別の手法でやってみるのも良しです.

ちなみにKaggleにはDigit Recognizerというのがあって,Forumとかを見れば結構わかりやすいコードだとかがありますし,自分で実装してmake a submissionで結果を送信するのも楽しいです.
でも大体みんなForumにある頭のいい人が書いたソースコードを実行して送信してるだけみたいな人が多いのでオリジナルのコードで上位に食い込むにはやっぱりディープラーニングとかを実装するしかない気がします.

6.参考

[1] MNIST handwritten digit database, Yann LeCun, CorinnaCortes and Chris Burges
[2] 酒井幸市 著,『画像処理とパターン認識入門』,森北出版株式会社
[3] 多層パーセプトロンで手書き数字認識 - 人工知能に関する断創録

2016年2月2日火曜日

振り子の動きを画像認識で検出してみる[Python]

前やった特定の色の物体の座標を検出の続きみたいなもんです.
振り子を作ってWebカメラを使って振り子の位置を検出して,グラフにしてみます.

まぁ今日のテストがイミフで,水曜日の第二外国語の勉強したくないだけですけど.

1.振り子

こんなの

文句はなしで行きましょう(笑).
暇つぶしクオリティなので重りはキャップを2つくっつけただけです.
糸も細くて変形しにくい針金などではなく,裁縫道具に入ってたやつです.

2.検出プログラム

こんな感じ.output.csvに時間,x座標,y座標が保存されます.

検出アルゴリズムとしては,HSV色空間にWebカメラの画像を変換し,緑色の成分のみをOpenCVのinRange関数を使って取り出しています.
なので以下のプログラムは「緑色の物体」の座標を認識しているので緑色のものがないところでやるようにしてください.

座標は緑色が一番たくさんあった行と列の要素番号をx,yの座標としてます.
まぁそのへんはプログラムを参照してください.

3.プロット

データがoutput.csvに入ったら次はmatplotlibで可視化してみます.可視化プログラムはこんな感じです.
グラフはこんなんが出てくる.これ横からではなく,下から撮影したのでこうなってます.

x,y座標の振り幅の極大極小はそれなりに取れてますね.振り幅が小さくなるとノイズが目立ちますが.(まぁあとカメラの傾きとか重力の方向とかどうなってんのみたいな話は置いといて……)


今日は以上です.
振り子の長さは大体37cmだったので,T=2*pi*sqrt(l/g)で,周波数はおおよそ0.8[Hz]になっているはずで,本当はここでフーリエ変換で周波数を求めようと思ったのですが……これ一定間隔でサンプリングしてないのです.

numpyやscipyのフーリエ変換の関数はたぶんサンプリング間隔が一定という条件でやってる気がする……時間のデータを与えられないし.
複素数をPythonで扱ったりとかはできず,現状実装できる能力がないので近いうちにフーリエ変換を実装しなければ…….