画像を一筆書きで描くプログラムを実装した[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]

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