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

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

水玉コラ作成ソフトの詳細情報 : Vector ソフトを探す！

一応こんな感じに動作する．

水玉コラ作成ソフトを作ってるのだがサンプル動画用に使えるフリー画像が見つからない．
サンプルだからできればあんまりエロく無くて金髪ツインテールの水着姿スレンダー系美少女が良いのだけれど．絵心ないと辛い(´Д｀)ﾊｧ… pic.twitter.com/XQ00t3h7Au

— もりとにー (@TonyMooori) 2016, 2月 13

これを使うことで慣れれば2,3分で水玉コラを作ることができるようになる．やったぜ．

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

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

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

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

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

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

これ( https://t.co/uMm0inKXeV )にインスパイアされて、円が１０個のバージョンを作ってみたらキモくなった pic.twitter.com/lUkBNNldy9

— どね (@donnay1224) 2016, 2月 5

k_1(x)=のところに好きな関数（数列）を入れて遊べるフーリエ級数視覚化マシーンを作りましたhttps://t.co/GmQo5NoZbz pic.twitter.com/vHrQ32FdWw

— 鯵坂もっちょ (@motcho_tw) 2016, 2月 6

という経緯です．

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

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

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

・ソースコード

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

	// フーリエ級数可視化クラス
	class FourierPlotter {
	int n_data; // データの個数
	int n_approx; // 近似次数
	float[][] radius; // 半径．radius[0][k]にsin成分のk次の成分が入る
	float[] data; // データ
	float centerX; // 中心X座標
	float centerY; // 中心Y座標
	float phase; // 位相のズレ
	public float x; // 先端部分のX座標
	public float y; // 先端部分のY座標
	// コンストラクタ
	FourierPlotter(float[] data, int n_approx,
	float centerX, float centerY, float phase) {
	// 引数をメンバ変数に代入した後，フーリへ級数展開
	this.data = data;
	this.n_data = data.length;
	this.n_approx = n_approx;
	this.phase = phase;
	this.centerX = this.x = centerX;
	this.centerY = this.y = centerY;
	this.calcRadius();
	}
	// フーリエ級数展開を計算する関数
	void calcRadius() {
	// フーリエ級数の係数は半径に当たる
	this.radius = new float[2][this.n_approx];
	float dt = 2.0*PI / this.n_data; // 微小長さdt
	// バイアス項(0次の項)は無視している
	for ( int i = 1; i < this.n_approx; i++ ) {
	for (int j = 0; j < this.n_data; j++ ) {
	float t = 2.0*PI*j/this.n_data;
	// 普通に矩形積分する
	this.radius[0][i] += sin(i*t)*this.data[j]*dt/PI; // sin成分の計算
	this.radius[1][i] += cos(i*t)*this.data[j]*dt/PI; // cos成分の計算
	}
	}
	}
	// 描画を行う関数
	void visualize(float theta) {
	// x,y座標をとりあえず中心に持ってくる
	this.x = this.centerX;
	this.y = this.centerY;
	float tempX, tempY;
	for (int i = 0; i < this.n_approx; i++ ) {
	for (int j= 0; j<2; j++) {
	// 円を描く
	ellipse(x, y, 2.0*this.radius[j][i], 2.0*this.radius[j][i]);
	tempX = x;
	tempY = y;
	// 先端部分の位置座標を計算(j=1の時，余弦なので位相がPI/2ずれる)
	x += this.radius[j][i]*cos(i*theta+PI*j/2 + this.phase);
	y += this.radius[j][i]*sin(i*theta+PI*j/2 + this.phase);
	line(x, y, tempX, tempY); // 各円の中心をを直線で結ぶ
	}
	}
	}
	}
	float scale = 1.0; // 拡大率
	float dx=0.0; // x方向の移動
	float dy=-0.0; // y方向の移動
	int n_data = 100; // データの個数
	FourierPlotter fp; // 描画クラス
	float[] points; // 過去に通った座標
	int count; // ループカウンタ
	int n_approx = 20; // 近似次数
	void setup() {
	size(640, 480);
	// 関数funcを離散的なデータにしてFourierPlotterに渡す
	float[] data = new float[n_data];
	for (int i = 0; i < n_data; i++ )
	data[i] = func( (float)(i)/n_data );
	fp = new FourierPlotter(data, n_approx, -100, 0, 0);
	// 諸々を初期化
	count = 0;
	points = new float[1000];
	}
	void draw() {
	background(0);
	translate(width/2-dx, height/2-dy); // 原点を中心に移動
	scale(scale, scale);
	fill(0, 0, 0, 0);
	strokeWeight(3);
	stroke( 0, 255, 0);
	float theta = count * 2*PI / n_data;
	fp.visualize(theta); // 描画
	points[count%points.length] = fp.y; // y座標を保存
	// 過去のデータの描画
	stroke(255,0,0);
	line(fp.x,fp.y,100.0,fp.y);
	plot(points,100,0,count%points.length);
	count +=1; // ループカウンタを増やす
	// フレームを保存する場合はコメントインする
	//saveFrame("./frames/#####.png");
	}
	// [0,1)区間で定義される区分的に滑らかな関数
	// 近似する関数の対象なのでここを変えればいい
	float func(float x) {
	// 矩形波
	return x < 0.5 ? -50:50;
	// サイン波
	//return 100*sin(2*x*PI);
	}
	// 配列arrをプロットする関数
	void plot(float[] arr,float x,float y,int startIndex){
	int N = arr.length;
	float y2 = arr[startIndex%N];
	for(int i = 1 ; i < N ; i++ ){
	line(x+N-i+1,y2,x+N-i,y+arr[(startIndex+i)%N]);
	y2 = arr[(startIndex+i)%N] + y;
	}
	}
	void keyPressed() {
	// キー操作による拡大・移動
	if (key == CODED) {
	if (keyCode == UP) {
	dy -= 10;
	} else if (keyCode == DOWN) {
	dy += 10;
	} else if (keyCode == LEFT) {
	dx -= 10;
	} else if (keyCode == RIGHT) {
	dx += 10;
	}
	} else if (key == '+') {
	scale *= 1.05;
	} else if (key=='-') {
	scale /= 1.05;
	}
	}

view raw VisualizingFourierSeries.pde hosted with ❤ by GitHub

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

1.概要

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

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

今回サンプルとして画像はこいつを使います．勝手に使ってください．

Processingしかできない人や，Pythonの実行環境ないひとは一筆書きの経路を探索した結果の座標ファイルを置いておきます．こいつをProcessingのソースコードと同一フォルダに入れておけば一応はできます．ちなみにこんな感じ．

ブログの記事( https://t.co/1ygaJ4I2WO )用
これでプログラムは全部公開したかな pic.twitter.com/ytK70SiiTM

— もりとにー (@TonyMooori) 2016, 2月 16

プログラムが長過ぎるので注意してください．もっとスマートにかければ良いのだけれど，無駄にクラスを作成して再利用しやすい形にした……(´・3・`)ﾄﾞｳｾﾂｶﾜﾅｲﾉﾆ

2.プログラム

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

	#coding:utf-8
	import cv2
	import numpy as np
	import time
	import matplotlib.pyplot as plt
	import pandas as pd
	from scipy.spatial import distance as dis
	"""
	参考URL
	[1] 蟻コロニー最適化 - Wikipedia https://ja.wikipedia.org/wiki/蟻コロニー最適化
	[2] 任意の確率密度分布に従う乱数の生成（von Neumannの棄却法） | Pacocat's Life http://pacocat.com/?p=596
	[3] Algorithms with Python / 巡回セールスマン問題 http://www.geocities.jp/m_hiroi/light/pyalgo64.html
	"""
	class TSP:
	def __init__(self,path=None,alpha = 1.0,beta = 1.0,Q = 1.0,vanish_ratio = 0.95):
	""" 初期化を行う関数 """
	self.alpha = alpha # フェロモンの優先度
	self.beta = beta # ヒューリスティック情報(距離)の優先度
	self.Q = Q # フェロモン変化量の係数
	self.vanish_ratio = vanish_ratio # 蒸発率
	if path is not None:
	self.set_loc(np.array(pd.read_csv(path)))
	def set_loc(self,locations):
	""" 位置座標を設定する関数 """
	self.loc = locations # x,y座標
	self.n_data = len(self.loc) # データ数
	self.dist = dis.squareform(dis.pdist(self.loc)) # 距離の表を作成
	self.weight = np.random.random_sample((self.n_data,self.n_data))+1.0 # フェロモンの量
	self.result = np.arange(self.n_data) # もっともよかった順序を保存する
	def cost(self,order):
	""" 指定された順序のコスト計算関数 """
	order2 = np.r_[order[1:],order[0]]
	return np.sum(tsp.dist[order,order2])
	def plot(self,order=None):
	""" 指定された順序でプロットする関数 """
	if order is None:
	plt.plot(self.loc[:,0],self.loc[:,1])
	else:
	plt.plot(self.loc[order,0],self.loc[order,1])
	plt.show()
	def solve(self,n_agent=1000):
	""" 巡回セールスマン問題を蟻コロニー最適化で解く """
	order = np.zeros(self.n_data,np.int) # 巡回経路
	delta = np.zeros((self.n_data,self.n_data)) #フェロモン変化量
	for k in range(n_agent):
	city = np.arange(self.n_data)
	now_city = np.random.randint(self.n_data) # 現在居る都市番号
	city = city[ city != now_city ]
	order[0] = now_city
	for j in range(1,self.n_data):
	upper = np.power(self.weight[now_city,city],self.alpha)*np.power(self.dist[now_city,city],-self.beta)
	evaluation = upper / np.sum(upper) # 評価関数
	percentage = evaluation / np.sum(evaluation) # 移動確率
	index = self.random_index2(percentage) # 移動先の要素番号取得
	# 状態の更新
	now_city = city[ index ]
	city = city[ city != now_city ]
	order[j] = now_city
	L = self.cost(order) # 経路のコストを計算
	# フェロモンの変化量を計算
	delta[:,:] = 0.0
	c = self.Q / L
	for j in range(self.n_data-1):
	delta[order[j],order[j+1]] = c
	delta[order[j+1],order[j]] = c
	# フェロモン更新
	self.weight *= self.vanish_ratio
	self.weight += delta
	# 今までで最も良ければ結果を更新
	if self.cost(self.result) > L:
	self.result = order.copy()
	# デバッグ用
	print("Agent ... %d,\t Now Cost %lf,\t Best Cost ... %lf" % (k,L,self.cost(self.result)))
	return self.result
	def save(self,out_path):
	""" 最もコストが低かった順序で保存する関数 """
	points = self.loc[ self.result ]
	f = open(out_path,"w")
	f.write("x,y\n")
	for i in range(len(points)):
	f.write(str(points[i,0]) + "," + str(points[i,1])+"\n")
	f.close()
	def random_index(self,percentage):
	""" 任意の確率分布に従って乱数を生成する関数 """
	n_percentage = len(percentage)
	arg = np.argsort(percentage)
	while True:
	index = np.random.randint(n_percentage)
	y = np.random.random()
	if y < percentage[index]:
	return index
	def random_index2(self,percentage):
	""" 精度低めで任意の確率分布に従って乱数を生成する関数 """
	n_percentage = len(percentage)
	arg = np.argsort(percentage)[::-1]
	n_arg = min(n_percentage,10)
	percentage = percentage / np.sum( percentage[arg] )
	while True:
	index = np.random.randint(n_arg)
	y = np.random.random()
	if y < percentage[arg[index]]:
	return arg[index]
	def save_edge_points(img_path,out_path):
	# 画像を読み込んでエッジ処理
	img = cv2.imread(img_path)
	edge = cv2.Canny(img,100,200)
	# エッジになっているx,y座標を取り出す
	h,w = edge.shape
	x = np.arange(w)
	y = np.arange(h)
	X,Y = np.meshgrid(x,y)
	# 255になっている部分がエッジ部分
	X_true = X[ edge > 128 ]
	Y_true = Y[ edge > 128 ]
	# エッジの点になっている座標が入っている
	index = np.array([X_true,Y_true]).T
	# 保存
	f = open(out_path,"w")
	f.write("x,y\n")
	for i in range(len(index)):
	f.write(str(index[i,0]) + "," + str(index[i,1])+"\n")
	f.close()
	if __name__=="__main__":
	# エッジを検出し保存．img.pngが対象の画像
	save_edge_points("img.png","edge_points.csv")
	# TSPで巡回セールスマン問題として一筆書きの手順を計算・保存
	tsp = TSP(path="edge_points.csv",alpha=1.0,beta=16.0,Q=1.0e3,vanish_ratio = 0.8)
	tsp.solve(100)
	tsp.save("best_order.csv")
	tsp.plot(tsp.result)

view raw ImageEdgeTSP.py hosted with ❤ by GitHub

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

	class FourierPlotter {
	int n_data; // データの個数
	int n_approx; // 近似次数
	float[][] radius; // 半径
	float[] data; // データ
	float centerX; // 中心X座標
	float centerY; // 中心Y座標
	float phase; // 位相のズレ
	public float x; // 先端部分のX座標
	public float y; // 先端部分のY座標
	// コンストラクタ
	FourierPlotter(float[] data, int n_approx,
	float centerX, float centerY, float phase) {
	// 引数をメンバ変数に代入した後，フーリへ級数展開
	this.data = data;
	this.n_data = data.length;
	this.n_approx = n_approx;
	this.phase = phase;
	this.centerX = this.x = centerX;
	this.centerY = this.y = centerY;
	this.calcRadius();
	}
	// フーリエ級数展開を計算する関数
	void calcRadius() {
	// フーリエ級数の係数は半径に当たる
	this.radius = new float[2][this.n_approx];
	float dt = 2.0*PI / this.n_data; // 微小長さdt
	// バイアス項(0次の項)は無視している
	for ( int i = 1; i < this.n_approx; i++ ) {
	for (int j = 0; j < this.n_data; j++ ) {
	float t = 2.0*PI*j/this.n_data;
	// 普通に矩形積分する
	this.radius[0][i] += sin(i*t)*this.data[j]*dt/PI;
	this.radius[1][i] += cos(i*t)*this.data[j]*dt/PI;
	}
	}
	}
	// 描画を行う関数
	void visualize(float theta) {
	// x,y座標をとりあえず中心に持ってくる
	this.x = this.centerX;
	this.y = this.centerY;
	float tempX, tempY;
	for (int i = 0; i < this.n_approx; i++ ) {
	for (int j= 0; j<2; j++) {
	// 円を描く
	ellipse(x, y, 2.0*this.radius[j][i], 2.0*this.radius[j][i]);
	tempX = x;
	tempY = y;
	// 先端部分の位置座標を計算(j=1の時，余弦なので位相がPI/2ずれる)
	x += this.radius[j][i]*cos(i*theta+PI*j/2 + this.phase);
	y += this.radius[j][i]*sin(i*theta+PI*j/2 + this.phase);
	line(x, y, tempX, tempY); // 各円の中心をを直線で結ぶ
	}
	}
	}
	}
	class Fourier2d {
	float[][] points; // 過去に通った座標
	int n_data; // データ数
	FourierPlotter fpX; // FourierPlotterクラスX方向
	FourierPlotter fpY; // FourierPlotterクラスY方向
	int count; // ループカウンタ
	int centerX;
	int centerY;
	// コンストラクタ
	Fourier2d(String filename,int centerX,int centerY) {
	float[][] data = load_data(filename);
	this.n_data = data[0].length;
	this.points = new float[2][n_data];
	this.count = 0;
	this.centerX = centerX;
	this.centerY = centerY;
	this.fpX = new FourierPlotter(data[0], n_data, centerX+150, centerY+150, -PI/2);
	this.fpY = new FourierPlotter(data[1], n_data, centerX-150, centerY-150, 0);
	}
	// x,yのデータを読み込む関数
	float[][] load_data(String path) {
	// 文字列を読み込む
	String lines[] = loadStrings(path);
	FloatList x = new FloatList();
	FloatList y = new FloatList();
	// 先頭行(ヘッダー)は読み飛ばすので1始まり
	for (int i=1; i < lines.length; i++) {
	// コンマで区切って2つ以上の要素があったらリストに加える
	String [] element =split(lines[i], ',');
	if (element.length >= 2 ) {
	x.append(parseFloat(element[0]));
	y.append(parseFloat(element[1]));
	}
	}
	float[][] data = new float[2][x.size()];
	data[0] = x.array();
	data[1] = y.array();
	return data;
	}
	void visualize( ) {
	float theta = this.count * 2*PI / this.n_data;
	fill(0, 0, 0, 0);
	strokeWeight(3);
	// 円を描画
	stroke( 0, 255, 0);
	this.fpX.visualize(theta);
	stroke( 255, 0, 0);
	this.fpY.visualize(theta);
	// 色を変える
	stroke(0, 255, 255);
	fill(0, 255, 255);
	strokeWeight(1);
	// 円の先端からの直線を描画
	line(fpX.x, fpX.y, fpX.x, fpY.y);
	line(fpY.x, fpY.y, fpX.x, fpY.y);
	rect(fpX.x, fpY.y, 1, 1);
	// 軌跡を描画
	stroke(255, 255, 255);
	fill(255, 255, 255);
	strokeWeight(3);
	this.points[0][this.count%n_data] = fpX.x;
	this.points[1][this.count%n_data] = fpY.y;
	for (int i = 0; i < min(this.n_data, this.count) - 1; i++ )
	line(points[0][i], points[1][i], points[0][i+1], points[1][i+1]);
	this.count+=1; // stepをすすめる
	}
	}
	float scale = 0.75; // 拡大率
	float dx=0.0; // x方向の移動
	float dy=-0.0; // y方向の移動
	Fourier2d f1; // 描画クラス
	void setup() {
	size(640, 480);
	// ここでデータを読み込むので同一フォルダにcsvファイルが必要
	f1 = new Fourier2d("best_order.csv", -100, 100);
	}
	void draw() {
	background(0);
	translate(width/2-dx, height/2-dy); // 原点を中心に移動
	scale(scale, scale);
	// 描画
	f1.visualize();
	//saveFrame("./frames/#####.png");
	}
	void keyPressed() {
	// キー操作による拡大・移動
	if (key == CODED) {
	if (keyCode == UP) {
	dy -= 10;
	} else if (keyCode == DOWN) {
	dy += 10;
	} else if (keyCode == LEFT) {
	dx -= 10;
	} else if (keyCode == RIGHT) {
	dx += 10;
	}
	} else if (key == '+') {
	scale *= 1.05;
	} else if (key=='-') {
	scale /= 1.05;
	}
	}

view raw DrawEdge.pde hosted with ❤ by GitHub

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

1.動画紹介

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



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

2.手法とか

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


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


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

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


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

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

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

巡回セールスマン問題(Traveling Salesman Problem,TSP)は有名だけれどなにげに解いたことがなかったのでちょいと実装してみました．

手法としてはいろいろあると思うのですが，今回は遺伝的アルゴリズム(GA)と蟻コロニー最適化(ACO)の2つを実装・比較しました．

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

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

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

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

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

• 巡回セールスマン問題 - Wikipedia
• Algorithms with Python / 巡回セールスマン問題
• 巡回セールスマン問題と遺伝的アルゴリズム(GA)

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

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

蟻コロニー最適化に関する数式は今回全てWikipediaに載っていたものを利用しました（蟻コロニー最適化 - Wikipedia）．

一応『心はプログラムできるか』という本を読んだ時に出てきてまして非常にわかりやすい解説がなされていましたので読んでみると良いかと思います．

以下ソースコード

	#coding:utf-8
	import numpy as np
	import matplotlib.pyplot as plt
	import pandas as pd
	from scipy.spatial import distance as dis
	"""
	参考URL
	[1] 蟻コロニー最適化 - Wikipedia https://ja.wikipedia.org/wiki/蟻コロニー最適化
	[2] 任意の確率密度分布に従う乱数の生成（von Neumannの棄却法） | Pacocat's Life http://pacocat.com/?p=596
	"""
	class TSP:
	def __init__(self,path=None,alpha = 1.0,beta = 1.0,Q = 1.0,vanish_ratio = 0.95):
	""" 初期化を行う関数 """
	self.alpha = alpha # フェロモンの優先度
	self.beta = beta # ヒューリスティック情報(距離)の優先度
	self.Q = Q # フェロモン変化量の係数
	self.vanish_ratio = vanish_ratio # 蒸発率
	if path is not None:
	self.set_loc(np.array(pd.read_csv(path)))
	def set_loc(self,locations):
	""" 位置座標を設定する関数 """
	self.loc = locations # x,y座標
	self.n_data = len(self.loc) # データ数
	self.dist = dis.squareform(dis.pdist(self.loc)) # 距離の表を作成
	self.weight = np.random.random_sample((self.n_data,self.n_data)) # フェロモンの量
	self.result = np.arange(self.n_data) # もっともよかった順序を保存する
	def cost(self,order):
	""" 指定された順序のコスト計算関数 """
	n_order = len(order)
	return np.sum( [ self.dist[order[i],order[(i+1)%n_order]] for i in np.arange(n_order) ] )
	def plot(self,order=None):
	""" 指定された順序でプロットする関数 """
	if order is None:
	plt.plot(self.loc[:,0],self.loc[:,1])
	else:
	plt.plot(self.loc[order,0],self.loc[order,1])
	plt.show()
	def solve(self,n_agent=1000):
	""" 巡回セールスマン問題を蟻コロニー最適化で解く """
	order = np.zeros(self.n_data,np.int) # 巡回経路
	delta = np.zeros((self.n_data,self.n_data)) #フェロモン変化量
	for k in range(n_agent):
	city = np.arange(self.n_data)
	now_city = np.random.randint(self.n_data) # 現在居る都市番号
	city = city[ city != now_city ]
	order[0] = now_city
	for j in range(1,self.n_data):
	upper = np.power(self.weight[now_city,city],self.alpha)*np.power(self.dist[now_city,city],-self.beta)
	evaluation = upper / np.sum(upper) # 評価関数
	percentage = evaluation / np.sum(evaluation) # 移動確率
	index = self.random_index(percentage) # 移動先の要素番号取得
	# 状態の更新
	now_city = city[ index ]
	city = city[ city != now_city ]
	order[j] = now_city
	L = self.cost(order) # 経路のコストを計算
	# フェロモンの変化量を計算
	delta[:,:] = 0.0
	c = self.Q / L
	for j in range(self.n_data-1):
	delta[order[j],order[j+1]] = c
	delta[order[j+1],order[j]] = c
	# フェロモン更新
	self.weight *= self.vanish_ratio
	self.weight += delta
	# 今までで最も良ければ結果を更新
	if self.cost(self.result) > L:
	self.result = order.copy()
	# デバッグ用
	print("Agent ... %d , Cost ... %lf" % (k,self.cost(self.result)))
	return self.result
	def random_index(self,percentage):
	""" 任意の確率分布に従って乱数を生成する関数 """
	n_percentage = len(percentage)
	while True:
	index = np.random.randint(n_percentage)
	y = np.random.random()
	if y < percentage[index]:
	return index
	if __name__=="__main__":
	#tsp = TSP(path="input.csv")
	tsp = TSP()
	tsp.set_loc(np.random.random_sample((30,2)))
	tsp.plot() # 計算前
	tsp.solve(n_agent=1000) # 1000匹の蟻を歩かせる
	tsp.plot(tsp.result) # 計算後

view raw TSP_ACO.py hosted with ❤ by GitHub

とりあえず評価のちゃんとした結果は遺伝的アルゴリズムの実装の後に書きます．

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

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

ソースコード

	#coding:utf-8
	import numpy as np
	import matplotlib.pyplot as plt
	import pandas as pd
	from scipy.spatial import distance as dis
	class TSP:
	def __init__(self,path=None,n_gene = 256,n_parent = 10,change_ratio = 0.1):
	""" 初期化を行う関数 """
	self.n_gene = n_gene # 一世代の遺伝子の個数
	self.n_parent = 10 # 親として残す個体数
	self.change_ratio = change_ratio # 突然変異で変化させる場所の数
	if path is not None:
	self.set_loc(np.array(pd.read_csv(path)))
	def init_genes(self,):
	""" 遺伝子をランダムに初期化 """
	self.genes = np.zeros((self.n_gene,self.n_data),np.int)
	order = np.arange(self.n_data)
	for i in range(self.n_gene):
	np.random.shuffle(order)
	self.genes[i] = order.copy()
	self.sort()
	def set_loc(self,locations):
	""" 位置座標を設定する関数 """
	self.loc = locations # x,y座標
	self.n_data = len(self.loc) # データ数
	self.dist = dis.squareform(dis.pdist(self.loc)) # 距離の表を作成
	self.init_genes() # 遺伝子を初期化
	def cost(self,order):
	""" 指定された順序のコスト計算関数 """
	return np.sum( [ self.dist[order[i],order[(i+1)%self.n_data]] for i in np.arange(self.n_data) ] )
	def plot(self,order=None):
	""" 指定された順序でプロットする関数 """
	if order is None:
	plt.plot(self.loc[:,0],self.loc[:,1])
	else:
	plt.plot(self.loc[order,0],self.loc[order,1])
	plt.show()
	def solve(self,n_step=1000):
	""" 遺伝的アルゴリズムで解く関数 """
	for i in range(n_step):
	print("Generation ... %d, Cost ... %lf" % (i,self.cost(self.genes[0])))
	self.step()
	self.result = self.genes[0]
	return self.result
	def step(self):
	""" 遺伝子を一世代分進化させる関数 """
	# 突然変異
	for i in range(self.n_parent,self.n_gene):
	self.genes[i] = self.mutation( np.random.randint(self.n_parent) )
	self.sort() # ソートする
	def sort(self):
	""" 遺伝子を昇順にする関数 """
	# コストを計算し，ソート
	gene_cost = np.array([self.cost(i) for i in self.genes])
	self.genes = self.genes[ np.argsort(gene_cost) ]
	def mutation(self,index):
	""" 突然変異を起こした遺伝子を返す関数 """
	n_change = int(self.change_ratio * self.n_data)
	gene = self.genes[index].copy()
	for i in range(n_change):
	# n_changeの個数だけ値を入れ替える
	left = np.random.randint(self.n_data)
	right = np.random.randint(self.n_data)
	temp = gene[left]
	gene[left] = gene[right]
	gene[right] = temp
	return gene
	if __name__=="__main__":
	#tsp = TSP(path="input.csv")
	tsp = TSP()
	tsp.set_loc(np.random.random_sample((30,2)))
	tsp.plot()
	tsp.solve()
	tsp.plot(tsp.result)

view raw TSP_GA.py hosted with ❤ by GitHub

4.評価

一応同じデータを用いてやった結果です．30個の[0,1)のランダムな二次元の点について行いました．
パラメータはデフォルトです．まずは遺伝的アルゴリズム．

直線が交差してますね．コストは6.056765でした．続いて蟻コロニー最適化．

直線が交差すること無く，コストも5.701502です．しかも計算時間が遺伝的アルゴリズムよりも遥かに実行時間が早いです．

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

離散コサイン変換(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もコメントアウトしてありますが一応実装してあります．
ベンチマークとして乱数で作った配列データを試しに変換・逆変換してみます．

	#coding: utf-8
	import numpy as np
	import matplotlib.pyplot as plt
	"""
	参考:
	[1]『画像処理とパターン認識入門』酒井幸市 著
	[2] scipy.fftpack.dct http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.fftpack.dct.html
	"""
	class DCT:
	def __init__(self,N):
	self.N = N # データ数．
	# 1次元離散コサイン変換の変換行列を予め作っておく
	self.phi_1d = np.array([ self.phi(i) for i in range(self.N) ])
	def dct(self,data):
	""" 1次元離散コサイン変換を行う """
	return self.phi_1d.dot(data)
	def idct(self,c):
	""" 1次元離散コサイン逆変換を行う """
	return np.sum( self.phi_1d.T * c ,axis=1)
	def phi(self,k):
	""" 離散コサイン変換(DCT)の基底関数 """
	# DCT-II
	if k == 0:
	return np.ones(self.N)/np.sqrt(self.N)
	else:
	return np.sqrt(2.0/self.N)*np.cos((k*np.pi/(2*self.N))*(np.arange(self.N)*2+1))
	# DCT-IV(試しに実装してみた)
	#return np.sqrt(2.0/N)*np.cos((np.pi*(k+0.5)/self.N)*(np.arange(self.N)+0.5))
	if __name__=="__main__":
	N = 100 # データ数を100とします
	dct = DCT(N) # 離散コサイン変換を行うクラスを作成
	x = np.random.random_sample(N) # N個の乱数データを作成
	c = dct.dct(x) # 離散コサイン変換を実行
	y = dct.idct(c) # 離散コサイン逆変換を実行
	# 元のデータ(x)と復元したデータ(y)をグラフにしてみる
	plt.plot(x,label="original")
	plt.plot(y,label="restored")
	plt.legend()
	plt.title("original data and restored data")
	plt.show()

view raw dct_1dim.py hosted with ❤ by GitHub

実行結果

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

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次元のフーリエ変換は実装されてますが．

サンプルとして適当な画像を変換・逆変換してみます．

	#coding: utf-8
	import numpy as np
	import matplotlib.pyplot as plt
	"""
	参考:
	[1]『画像処理とパターン認識入門』酒井幸市 著
	[2] scipy.fftpack.dct http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.fftpack.dct.html
	"""
	class DCT:
	def __init__(self,N):
	self.N = N # データ数．
	# 1次元，2次元離散コサイン変換の基底ベクトルをあらかじめ作っておく
	self.phi_1d = np.array([ self.phi(i) for i in range(self.N) ])
	# Nが大きいとメモリリークを起こすので注意
	# MNISTの28x28程度なら問題ない
	self.phi_2d = np.zeros((N,N,N,N))
	for i in range(N):
	for j in range(N):
	phi_i,phi_j = np.meshgrid(self.phi_1d[i],self.phi_1d[j])
	self.phi_2d[i,j] = phi_i*phi_j
	def dct(self,data):
	""" 1次元離散コサイン変換を行う """
	return self.phi_1d.dot(data)
	def idct(self,c):
	""" 1次元離散コサイン逆変換を行う """
	return np.sum( self.phi_1d.T * c ,axis=1)
	def dct2(self,data):
	""" 2次元離散コサイン変換を行う """
	return np.sum(self.phi_2d.reshape(N*N,N*N)*data.reshape(N*N),axis=1).reshape(N,N)
	def idct2(self,c):
	""" 2次元離散コサイン逆変換を行う """
	return np.sum((c.reshape(N,N,1)*self.phi_2d.reshape(N,N,N*N)).reshape(N*N,N*N),axis=0).reshape(N,N)
	def phi(self,k):
	""" 離散コサイン変換(DCT)の基底関数 """
	# DCT-II
	if k == 0:
	return np.ones(self.N)/np.sqrt(self.N)
	else:
	return np.sqrt(2.0/self.N)*np.cos((k*np.pi/(2*self.N))*(np.arange(self.N)*2+1))
	# DCT-IV(試しに実装してみた)
	#return np.sqrt(2.0/N)*np.cos((np.pi*(k+0.5)/self.N)*(np.arange(self.N)+0.5))
	if __name__=="__main__":
	N = 10 # データの次元は10x10とする
	dct = DCT(10) # 離散コサイン変換を行うクラスを作成
	# サンプル画像を作る
	img = np.array([
	[0,0,0,0,0,0,0,0,0,0],
	[0,0,1,1,1,1,1,1,0,0],
	[0,1,0,0,0,0,0,0,1,0],
	[0,1,0,0,0,0,0,0,1,0],
	[0,0,0,0,0,0,1,1,0,0],
	[0,0,0,0,1,1,0,0,0,0],
	[0,0,0,1,0,0,0,0,0,0],
	[0,0,1,0,0,0,0,0,0,0],
	[0,1,1,1,1,1,1,1,1,0],
	[0,0,0,0,0,0,0,0,0,0]
	])
	c = dct.dct2(img) # 2次元離散コサイン変換
	y = dct.idct2(c) # 2次元離散コサイン逆変換
	# 元の画像と復元したものを表示
	plt.subplot(1,2,1)
	plt.imshow(img,cmap="Greys")
	plt.title("original")
	plt.xticks([])
	plt.subplot(1,2,2)
	plt.imshow(y,cmap="Greys")
	plt.title("restored")
	plt.xticks([])
	plt.show()

view raw dct_2dim.py hosted with ❤ by GitHub

実行結果です．

画像がちゃんと元の状態に戻っているのがわかります．ちなみにDCT係数（周波数成分）はこんな感じでした．

概ね右下に行くほど色が薄くなっています．ここから低周波成分(左上）が多くを占めていることがわかります．
高周波成分を消していったものがこちらです．

DCT係数cについてc[i:,i:] = 0とした時に復元された画像です．i=3でだいたい認識できそうですね．結構データを圧縮できてるのではないでしょうか．まぁまぁうまく行って良かったです．

相関係数を用いて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.ソースコード

	#coding: utf-8
	import numpy as np
	import matplotlib.pyplot as plt
	from sklearn.datasets import fetch_mldata
	from sklearn.cross_validation import train_test_split
	from sklearn.metrics import confusion_matrix, classification_report
	def corrcoef(x,y):
	""" ベクトルx,yの相関係数を求める関数 """
	# ノルムが0になる場合はないものとします
	# 要素の平均を引く
	x = x - np.average(x)
	y = y - np.average(y)
	# ベクトルの内積はベクトルを転置して行列の積となる
	return x.T.dot(y)/(np.linalg.norm(x)*np.linalg.norm(y))
	# mnistの手書き数字データをロード
	# カレントディレクトリ（.）にない場合は、
	# Webから自動的にダウンロードされる（時間がかかるので注意！）
	# 70000サンプル，28x28ピクセル
	mnist = fetch_mldata('MNIST original', data_home=".")
	"""
	# 100個適当に表示してみる
	for i in range(10):
	for j in range(10):
	img = mnist.data[ mnist.target == i ][j]
	plt.subplot(10, 10, 10 * i + j + 1)
	plt.axis('off')
	plt.imshow( img.reshape(28, 28), cmap=plt.cm.gray_r, interpolation='nearest')
	plt.show()
	"""
	# 訓練データを作成
	X = mnist.data
	y = mnist.target
	# 訓練データとテストデータに分解
	X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=1./7.)
	# 代表値を作成する
	avg_img = np.zeros((10,X_train.shape[1]))
	for i in range(10):
	# ラベルがiの画像を取り出して，平均をとる
	i_imgs = X_train[ y_train == i ]
	avg_img[i] = np.average(i_imgs,axis=0)
	"""
	# 代表値(平均画像)を表示してみる
	for i in range(10):
	img = avg_img[i]
	plt.subplot(1, 10, i + 1)
	plt.axis('off')
	plt.imshow( img.reshape(28, 28), cmap=plt.cm.gray_r, interpolation='nearest')
	plt.show()
	"""
	# 予測
	predictions = np.zeros(y_test.shape)
	for i in range(len(y_test)):
	test = X_test[i].astype(np.float)
	# 内積が最も大きいものを予想ラベルとする
	predictions[i] = np.argmax( [ corrcoef(test,img) for img in avg_img ] )
	# 評価
	print (confusion_matrix(y_test, predictions))
	print (classification_report(y_test, predictions))

view raw mnist_corrcoef.py hosted with ❤ by GitHub

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] 多層パーセプトロンで手書き数字認識 - 人工知能に関する断創録

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

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

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

1.振り子

こんなの

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

2.検出プログラム

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

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

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

	# -*- coding: utf-8 -*-
	import numpy as np
	import cv2
	import time
	cap = cv2.VideoCapture(0)
	f = open("output.csv","w")
	f.write("t,x,y\n")
	start = time.time()
	while(True):
	# フレームをキャプチャする
	ret, frame = cap.read()
	# hsv色空間に変換
	hsv = cv2.cvtColor(frame,cv2.COLOR_BGR2HSV)
	# 0 <= h <= 179 (色相)　OpenCVではmax=179なのでR:0(180),G:60,B:120となる
	# 0 <= s <= 255 (彩度)　黒や白の値が抽出されるときはこの閾値を大きくする
	# 0 <= v <= 255 (明度)　これが大きいと明るく，小さいと暗い
	# ここでは緑色を抽出するので60±20を閾値とした
	low_color = np.array([40, 75, 75])
	upper_color = np.array([80, 255, 255])
	# 色を抽出する
	ex_img = cv2.inRange(hsv,low_color,upper_color)
	# 縦，横方向の和の中で最も大きいインデックスを抽出
	x = np.argmax(np.sum(ex_img,axis=0))
	y = np.argmax(np.sum(ex_img,axis=1))
	# 時間と座標の保存
	t = time.time() - start
	f.write(str(t)+","+str(x)+","+str(y)+"\n")
	print(str(t)+","+str(x)+","+str(y)+"\n")
	# 抽出した座標に丸を描く
	cv2.circle(frame,(x,y),10,(0,255,0))
	# 画面に表示する
	cv2.imshow('frame',frame)
	if t > 120:
	break
	key = cv2.waitKey(1) & 0xFF
	# qが押された場合は終了する
	if key == ord('q'):
	break
	# 後始末
	cap.release()
	cv2.destroyAllWindows()
	f.close()

view raw DetectGreen.py hosted with ❤ by GitHub

3.プロット

データがoutput.csvに入ったら次はmatplotlibで可視化してみます．可視化プログラムはこんな感じです．

	#coding:utf-8
	import numpy as np
	import matplotlib.pyplot as plt
	import pandas as pd
	# 最初のデータはWebカメラの調子が悪いため無視
	data = pd.read_csv("output.csv")
	data = data[50:]
	# 計算しやすいように取り出す
	t = data['t']
	x = data['x']
	y = data['y']
	# x,yは平均(振動の中心)を0に持ってこさせる
	x = x - np.average(x)
	y = y - np.average(y)
	# 時間はt[0]が0になるようにしとく
	# (pandasを使うとt[0]という表現ができないのでこうしました)
	t = t - np.min(t)
	# 凡例とともに表示
	plt.plot(t,x,label="x")
	plt.plot(t,y,label="y")
	plt.legend()
	plt.show()

view raw PendulumPlot.py hosted with ❤ by GitHub

グラフはこんなんが出てくる．これ横からではなく，下から撮影したのでこうなってます．

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


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

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