2009年7月24日金曜日

Processingを使ってWebカメラを監視カメラにする

このブログ記事をはてなブックマークに追加

Processingvideoライブラリを利用して、手持ちのWebカメラを防犯・監視カメラにしてみることにした。監視カメラとして、以下の2点を満たすこととする。

  • 動くものを捕らえたときだけ写真を残す。
  • 写真は指定のメールアドレスに直ちに送信される。

今回、コードを書くに当たって、Webカメラ関連については建築発明工作ゼミ2008、メール送信についてはE-mail Processingを参考にさせてもらった。特に建築発明工作ゼミ2008については他にもProcessing関連の有益な情報があって重宝する。

因みに、今回用いたWebカメラはずいぶん前に購入したLogitechのQuickCam for Notebooks Proだ。多分、大抵のWebカメラで問題なく動作すると思う。

それでは、早速作成してみる。

下準備

まず、ProcessingでWebカメラを利用するには、QuickTimeが必要になるので予め用意しておく。また、WindowsについてはQuickTimeでWebカメラを利用するためのVDIG (QuickTime-compatible video digitizer)をインストールする必要がある。このとき、最新バージョン(1.04/1.05)では動作時にエラーが起こるようなので、バージョンは1.01を用意すること。

次に、Processingでメール送信を行うために、JavaMailJavaBeans Activation Framework (JAF)をダウンロードする。JavaMailからmailapi.jarとsmtp.jar、JAFからactivation.jarを取り出す。そして、Processingの作業フォルダ(コーディングしているファイルの存在するフォルダ)にcodeという名前のフォルダを作成し、その中にその3つのjarファイルをコピーする。

これで下準備は完了だ。

コーディング

ライブラリについては、ビデオ関連でprocessing.video、電子メール関連でjavax.mailをインポートする。setup関数でvideoオブジェクトを作成し、draw関数で映像を表示する。今回は320×240の解像度に設定した。次にカメラの映像から、以前の画像と現在の画像との差をピクセルごとで比較し、その差が許容値(torelance)以上であれば動体とし、それらすべてのピクセルを数え上げる。そして、その動体として認識されたピクセルの数が一定数(今回は100に設定)以上であれば写真を撮る。フレームごとに写真を撮ってしまうと膨大なファイル数になってしまうので、シャッター間隔は最小で5秒とした。また、写真を保存するフォルダはPATH変数に設定しておく。

次に、撮った写真をメールで送信するために、sendmail関数を実装し、ファイル添付についてはMIMEを利用した。今回のコードではGmailのSMTPサーバで送信するので、認証のためのアカウントとパスワードを予め設定しておく必要がある。もし、自前でSMTPサーバを用意できるのであれば、認証の必要はない。その場合、

props.put("mail.smtp.port", PORT); props.put("mail.smtp.auth", "true"); props.put("mail.smtp.starttls.enable", "true"); Session session = Session.getDefaultInstance(props, new Auth(FROM, PASS));

を、

Session session = Session.getDefaultInstance(props,null);

に変更する。また、Authクラスも、PORT/PASS変数も必要なくなる。

今回のコーディングは結構面白かった。メールの送信で若干もたつくけど、まあ実用にはなると思う。それにしても、こんなに簡単に、短いコードでWebカメラを監視カメラに仕立て上げることができるとは、Processingの可能性を再認識させられた。



watching_camera.pde

import processing.video.*; import javax.mail.*; String SMTP = "smtp.gmail.com"; String FROM = "your_account@gmail.com"; String TO = "your@mail.address"; // Gmail Authentication String PORT = "587"; // or "465" String PASS = "your_password"; // 画像フォルダ. String PATH = "C:/your_folder/Processing/watching_camera/"; // 映像オブジェクト. Capture video; int W = 320; int H = 240; color[] ex_color = new color[W * H]; int tolerance = 50; // 許容値. float timer = 0; void setup() { size(W, H); video = new Capture(this, W, H, 30); // 30fps noStroke(); } void draw() { if (video.available()) { background(0); video.read(); set(0, 0, video); loadPixels(); int num_pixels = 0; for (int i = 0;i < W * H; i++) { // 前回と現在のピクセルの差. float diff_R = abs(red(ex_color[i]) - red(video.pixels[i])); float diff_G = abs(green(ex_color[i]) - green(video.pixels[i])); float diff_B = abs(blue(ex_color[i]) - blue(video.pixels[i])); // 許容値以上の差があれば動体. if(diff_R > tolerance && diff_G > tolerance && diff_B > tolerance) num_pixels++; // 現在のピクセルを保存. ex_color[i] = video.pixels[i]; } // 撮影間隔は5秒以上. 差のあるピクセルの数が100以上あれば撮影. if (millis() - timer > 5000 && num_pixels > 100) { DecimalFormat formatter = new DecimalFormat("00"); String date = str(year()); date += formatter.format(month()); date += formatter.format(day()); date += formatter.format(hour()); date += formatter.format(minute()); date += formatter.format(second()); String filename = "img" + date + ".jpg"; save(PATH + filename); timer = millis(); println("shot: " + date); // 画像ファイル添付のメールを送信. sendmail(date, filename, PATH); } } } void sendmail(String date, String filename, String path) { Properties props = System.getProperties(); props.put("mail.transport.protocol", "smtp"); props.put("mail.smtp.host", SMTP); // Gmail Authentication props.put("mail.smtp.port", PORT); props.put("mail.smtp.auth", "true"); props.put("mail.smtp.starttls.enable", "true"); Session session = Session.getDefaultInstance(props, new Auth(FROM, PASS)); MimeMessage mimeMessage = new MimeMessage(session); try { mimeMessage.setFrom(new InternetAddress(FROM)); mimeMessage.setRecipients(Message.RecipientType.TO, TO); mimeMessage.setSubject("SHOT: " + date, "iso-2022-jp"); Multipart mp = new MimeMultipart(); MimeBodyPart mbp; // メール本文. mbp = new MimeBodyPart(); mbp.setText("撮影日時: " + date, "iso-2022-jp"); mbp.setHeader("Content-Type", "text/plain"); mp.addBodyPart(mbp); // ファイル添付. mbp = new MimeBodyPart(); FileDataSource fds = new FileDataSource(path + filename); DataHandler imgdh = new DataHandler(fds); mbp.setDataHandler(imgdh); mbp.setFileName(MimeUtility.encodeWord(filename)); mp.addBodyPart(mbp); mimeMessage.setContent(mp); mimeMessage.setSentDate(new Date()); // メール送信. Transport.send(mimeMessage); } catch (Exception e) { e.printStackTrace(); } }

Auth.pde

import javax.mail.Authenticator; import javax.mail.PasswordAuthentication; public class Auth extends Authenticator { private String username, password; public Auth(String username_, String password_) { super(); username = username_; password = password_; } public PasswordAuthentication getPasswordAuthentication() { return new PasswordAuthentication(username, password); } }

続きを読む...

2009年7月21日火曜日

Google App Engine: 簡単にグラフ・チャートを作成する

このブログ記事をはてなブックマークに追加

ブログなどでちょっとしたグラフを表示させたい場合、Google Chart APIを使っていた。しかし、手作業で入力するのはかなり面倒だ。ウェブ上にはグラフやチャートを作成するサービスなどもあるが、登録が必要だったり、手順が面倒だったりと、個人的にはあまり手軽だと思えない。数値をコピー&ペーストしてクリック一つでグラフを作成したいのだ。そこで、使いたいサービスは自分で作ってしまえということで、Google App Engine (GAE)とGoogle Chart APIを使って簡単にグラフ・チャートを作成するウェブアプリを作ってみた。



即席で作成したウェブアプリなので不備な点もあるが、取り敢えず自分で使う分にはこの程度で十分なのでGAEに登録しておいた。気が向いたら機能を拡張していくかもしれない。

以下に使い方を示す。

使い方

まず、以下のデータをテキストエリアに入れてみる。

1.2 2.2 3.5 5.6 2.8

次に、「作成」ボタンを押す。チャートの種類が「折れ線グラフX」になっていれば、入力データをもとに以下のようなグラフが描かれる。


グラフのX軸、Y軸の下限と上限、およびその目盛り分割は自動で行われる。もちろん、手動での設定もできるし、自動で設定された数値を修正することも可能だ。また、出力画像ファイルの大きさも、グラフの色(16進数6桁[RGB]で表示)も変更できる。タイトルを表示したければタイトル欄に入力すればよい。因みにタイトルを改行させるには | を入れる。さらに、下記のようなGoogle Chart APIによるURLも表示するようにした。このURLを画像URLとしてブログなどで利用すれば、上のグラフを表示させるのも簡単だ。

http://chart.apis.google.com/chart?chs=300x200&chd=t:1.2,2.2,3.5,5.6,2.8&chds=1,6&cht=lc&chco=4d89f9,c6d9fd&chxt=y&chxl=0:|1|2|3|4|5|6

現時点で利用できるグラフの種類は、「折れ線グラフX」、「折れ線グラフXY」、「スパークライン」、「積上横棒グラフ」、「積上縦棒グラフ」、「集合横棒グラフ」、「集合縦棒グラフ」、「円グラフ」、「3D円グラフ」、「散布図」の10種類だ。2カラム以上のデータが必要なグラフは、「折れ線グラフXY」と「散布図」で、ほかは1カラム以上のデータとなる。また、「転置」チェックボックスにチェックを入れると行と列を入れ替えることができ、「ラベル」チェックボックスにチェックを入れると、入力データの1カラム目がラベルとなる。ラベル内に空白を入れたい場合は、空白を+に置き換えること。

データA 0.113 データB 0.243 データC 0.543 データD 0.443

上記のデータを入力し、「集合横棒グラフ」で作成すると、以下の棒グラフが描画される。因みに横棒グラフについては、X軸とY軸が入れ替わるので注意して欲しい。


また、「円グラフ」、「3D円グラフ」で作成すると以下の出力になる。



2カラムの数値データとして以下のデータを入力する。

0.113 0.223 0.343 0.285 0.443 0.386 0.543 0.192 0.723 0.256

これを、「折れ線グラフXY」で作成すると、以下の折れ線グラフを出力する。


さらに、同じデータまま、「積上横棒グラフ」を作成。


そして、「散布図」は次のようになる。


以下にソースコードを示す。

# 簡易グラフ・チャート作成 class Chart(webapp.RequestHandler): def __init__(self): self.chart_url = "" self.chart_html = "" self.chart_form = "" self.width = 300 self.height = 200 self.has_label = False self.is_trans = False self.title = "" self.color = "4d89f9,c6d9fd" self.is_clear_x = False self.is_clear_y = False self.min_x = None self.max_x = None self.step_x = None self.min_y = None self.max_y = None self.step_y = None self.html = u""" <html> <head> <link type="text/css" rel="stylesheet" href="/stylesheets/chart.css" /> <title>簡易グラフ・チャート作成</title> </head> <body> <div class="NoBoxedItem"> <form action="/chart" method="post"> <div><select name="chart"> <option value="lc">折れ線グラフ X</option> <option value="lxy">折れ線グラフ XY</option> <option value="ls">スパークライン</option> <option value="bhs">積上横棒グラフ</option> <option value="bvs">積上縦棒グラフ</option> <option value="bhg">集合横棒グラフ</option> <option value="bvg">集合縦棒グラフ</option> <option value="p">円グラフ</option> <option value="p3">3D円グラフ</option> <option value="s">散布図</option> </select> 幅:<input type="text" name="width" size="4" value="%d"></text> 高さ:<input type="text" name="height" size="4" value="%d"></text></div> <div>タイトル: <input type="text" name="title" size="40" value="%s"></text></div> <div><textarea name="content" rows="10" cols="60">%s</textarea></div> <div>X (クリア<input type="checkbox" name="clear_x" value="clear_x" />): 下限:<input type="text" name="min_x" size="3" value="%s"></text> 上限:<input type="text" name="max_x" size="3" value="%s"> 分割:<input type="text" name="step_x" size="2" value="%s"></div> <div>Y (クリア<input type="checkbox" name="clear_y" value="clear_y" />): 下限:<input type="text" name="min_y" size="3" value="%s"></text> 上限:<input type="text" name="max_y" size="3" value="%s"> 分割:<input type="text" name="step_y" size="2" value="%s"></div> <div>色: <input type="text" name="color" size="40" value="%s"></text></div> <div><input type="submit" value="作成" /> <input type="checkbox" name="label" value="label" />ラベル <input type="checkbox" name="trans" value="trans" />転置 <span><a href="./chart" onclick="window.open('./chart'); return false;">新規</a> <a href="http://handasse.blogspot.com/2009/07/google-app-engine.html" onclick="window.open('http://handasse.blogspot.com/2009/07/google-app-engine.html'); return false;">ヘルプ</a></span></div> </form> </div> <div class="NoBoxedItem"> %s %s </div> <div class="Contents"><span style="font-size:x-small;color:#666666;">Copyright &copy; 2009 <a href="http://handasse.blogspot.com/">nox</a>. All rights reserved.</span></div> </body> </html>""" def get(self): self.response.out.write(self.html % (self.width, self.height, "", "", "", "", "", "", "", "", self.color, "", "")) def get_scale(self, min_data, max_data, axis): stepsize_log = math.log10(max_data - min_data) stepsize = 10**int(stepsize_log) if stepsize_log < 0: stepsize = Decimal(str(stepsize)) * Decimal("0.1") scalesize = float(Decimal("100") / stepsize) stepsize = float(stepsize) else: scalesize = 1 min_scale = min_data // stepsize * stepsize redun = 1 if Decimal(str(max_data)) % Decimal(str(stepsize)) == Decimal("0"): redun = 0 max_scale = float((Decimal("%.10f" % max_data) // Decimal("%.10f" % stepsize) + redun) * Decimal("%.10f" % stepsize)) if stepsize_log >= 0 and stepsize_log == int(stepsize_log): stepsize *= 0.1 if stepsize < 1: scalesize = 100.0 / stepsize if axis == "x": if self.is_clear_x or self.min_x == None or self.max_x == None or self.step_x == None: self.min_x = min_scale self.max_x = max_scale self.step_x = int(Decimal(str(max_scale - min_scale)) / Decimal(str(stepsize))) else: min_scale = self.min_x max_scale = self.max_x stepsize = float(Decimal(str(max_scale - min_scale)) / Decimal(str(self.step_x))) if stepsize < 1: scalesize = 100.0 / stepsize elif axis == "y": if self.is_clear_y or self.min_y == None or self.max_y == None or self.step_y == None: self.min_y = min_scale self.max_y = max_scale self.step_y = int(Decimal(str(max_scale - min_scale)) / Decimal(str(stepsize))) else: min_scale = self.min_y max_scale = self.max_y stepsize = float(Decimal(str(max_scale - min_scale)) / Decimal(str(self.step_y))) if stepsize < 1: scalesize = 100.0 / stepsize return min_scale, max_scale, stepsize, scalesize def post(self): try: content = cgi.escape(self.request.get('content')) cht = cgi.escape(self.request.get('chart')) self.width = int(cgi.escape(self.request.get('width'))) self.height = int(cgi.escape(self.request.get('height'))) self.title = cgi.escape(self.request.get('title')).strip() if self.title: chtt = "&chtt=%s" % self.title.replace(" ", "+") else: chtt = "" min_x = cgi.escape(self.request.get('min_x')) if min_x: self.min_x = float(min_x) max_x = cgi.escape(self.request.get('max_x')) if max_x: self.max_x = float(max_x) step_x = cgi.escape(self.request.get('step_x')) if step_x: self.step_x = int(step_x) min_y = cgi.escape(self.request.get('min_y')) if min_y: self.min_y = float(min_y) max_y = cgi.escape(self.request.get('max_y')) if max_y: self.max_y = float(max_y) step_y = cgi.escape(self.request.get('step_y')) if step_y: self.step_y = int(step_y) self.color = cgi.escape(self.request.get('color')).replace(" ", "") if self.request.get('clear_x') == 'clear_x': self.is_clear_x = True if self.request.get('clear_y') == 'clear_y': self.is_clear_y = True if self.request.get('label') == 'label': self.has_label = True else: self.has_label = False if self.request.get('trans') == 'trans': self.is_trans = True else: self.is_trans = False if content: lines = [x.strip().replace(",", " ").replace("\t", " ") for x in content.rstrip().split("\n")] sdata = [[d for d in l.split()] for l in lines] if not self.is_trans: sdata = map(list, zip(*sdata)) if self.has_label: labels = sdata[0] sdata = sdata[1:] data = [[float(d) for d in sd] for sd in sdata] if cht == "lxy" or cht == "s": min_data_x, max_data_x = data[0][0], data[0][0] d = data[0] min_data_x, max_data_x = min(min_data_x, min(d)), max(max_data_x, max(d)) min_scale_x, max_scale_x, stepsize_x, scalesize_x = self.get_scale(min_data_x, max_data_x, "x") min_data_y, max_data_y = data[1][0], data[1][0] d = data[1] min_data_y, max_data_y = min(min_data_y, min(d)), max(max_data_y, max(d)) min_scale_y, max_scale_y, stepsize_y, scalesize_y = self.get_scale(min_data_y, max_data_y, "y") min_data_x, max_data_x = "%f" % min_scale_x, "%f" % max_scale_x min_data_y, max_data_y = "%f" % min_scale_y, "%f" % max_scale_y elif cht == "bhs" or cht == "bvs": tdata = map(sum, map(list, zip(*data))) min_data_y, max_data_y = tdata[0], tdata[0] for d in tdata: min_data_y, max_data_y = min(min_data_y, d), max(max_data_y, d) min_data_y = 0 min_scale_y, max_scale_y, stepsize_y, scalesize_y = self.get_scale(min_data_y, max_data_y, "y") min_data_y, max_data_y = "%f" % min_scale_y, "%f" % max_scale_y else: min_data_y, max_data_y = data[0][0], data[0][0] for d in data: min_data_y, max_data_y = min(min_data_y, min(d)), max(max_data_y, max(d)) if cht == "bhg" or cht == "bvg" or cht == "p" or cht == "p3": min_data_y = 0 min_scale_y, max_scale_y, stepsize_y, scalesize_y = self.get_scale(min_data_y, max_data_y, "y") min_data_y, max_data_y = "%f" % min_scale_y, "%f" % max_scale_y sdata = [",".join(d) for d in sdata] chd = "t:" + "|".join(sdata) if cht == "lxy" or cht == "s": chds = "%s,%s,%s,%s" % (min_data_x.rstrip("0").rstrip("."), max_data_x.rstrip("0").rstrip("."), min_data_y.rstrip("0").rstrip("."), max_data_y.rstrip("0").rstrip(".")) else: chds = "%s,%s" % (min_data_y.rstrip("0").rstrip("."), max_data_y.rstrip("0").rstrip(".")) chco = self.color if cht == "bhg" or cht == "bhs": if self.has_label: chxt = "x,y" chxl = "0:|" + "|".join([str(n / scalesize_y) for n in range(int(min_scale_y * scalesize_y), int(max_scale_y * scalesize_y) + 1, int(stepsize_y * scalesize_y))]) chxl += "|1:|" + "|".join(labels) else: chxt = "x" chxl = "0:|" + "|".join([str(n / scalesize_y) for n in range(int(min_scale_y * scalesize_y), int(max_scale_y * scalesize_y) + 1, int(stepsize_y * scalesize_y))]) elif cht == "lxy" or cht == "s": chxt = "x,y" chxl = "0:|" + "|".join([str(n / scalesize_x) for n in range(int(min_scale_x * scalesize_x), int(max_scale_x * scalesize_x) + 1, int(stepsize_x * scalesize_x))]) chxl += "|1:|" + "|".join([str(n / scalesize_y) for n in range(int(min_scale_y * scalesize_y), int(max_scale_y * scalesize_y) + 1, int(stepsize_y * scalesize_y))]) elif cht == "p" or cht == "p3": if self.has_label: chxt = "x" chxl = "0:|" + "|".join(labels) else: chxt = "" chxl = "" else: if self.has_label: chxt = "y,x" chxl = "0:|" + "|".join([str(n / scalesize_y) for n in range(int(min_scale_y * scalesize_y), int(max_scale_y * scalesize_y) + 1, int(stepsize_y * scalesize_y))]) chxl += "|1:|" + "|".join(labels) else: chxt = "y" chxl = "0:|" + "|".join([str(n / scalesize_y) for n in range(int(min_scale_y * scalesize_y), int(max_scale_y * scalesize_y) + 1, int(stepsize_y * scalesize_y))]) self.chart_url = "http://chart.apis.google.com/chart?chs=%dx%d&chd=%s&chds=%s&cht=%s&chco=%s&chxt=%s&chxl=%s%s" % (self.width, self.height, chd, chds, cht, chco, chxt, chxl, chtt) self.chart_html = '<div><img src="%s" /></div>' % self.chart_url self.chart_form = u""" <form name="form" action="/chart" method="post"> <div><textarea name="content" rows="5" cols="60">%s</textarea></div> <div><input type="button" onclick="document.form.content.select()" value="選択" /></div> </form>""" % self.chart_url if self.min_x == None or self.max_x == None or self.step_x == None: min_x = max_x = step_x = "" else: min_x, max_x, step_x = (str(self.min_x), str(self.max_x), str(self.step_x)) if self.min_y == None or self.max_y == None or self.step_y == None: min_y = max_y = step_y = "" else: min_y, max_y, step_y = (str(self.min_y), str(self.max_y), str(self.step_y)) if self.has_label: self.html = self.html.replace(u"label\"", u"label\" checked") if self.is_trans: self.html = self.html.replace(u"trans\"", u"trans\" checked") self.response.out.write(self.html.replace("\"%s\"" % cht, "\"%s\" selected" % cht) % (self.width, self.height, self.title, content, min_x, max_x, step_x, min_y, max_y, step_y, self.color, self.chart_html, self.chart_form)) except: self.response.out.write(self.html.replace(u"\"%s\"" % cht, u"\"%s\" selected" % cht) % (self.width, self.height, self.title, content, min_x, max_x, step_x, min_y, max_y, step_y, self.color, u"エラーが発生しました。データおよび設定を確認してください。", ""))

続きを読む...

2009年7月9日木曜日

Python: 階層的クラスタリングのデンドログラム描画と閾値による区分け

このブログ記事をはてなブックマークに追加

前回、「C++: マルチコアCPUを利用した並列化による高速な階層的クラスタリング」でクラスタリングを行ったのだが、ここではその出力データをPython+PILにより、デンドログラムを描画したり、指定した閾値で区分けを行ってみる。

まず、クラスタリングの出力データは前回のtest.outを使用する。内容は以下の通りだ。

0.0833487 4 14 0.11483 3 15 0.123895 0 5 0.126783 7 11 0.144271 16 -4 0.14854 8 9 0.253505 -5 -6 0.264889 -1 -3 0.301108 2 10 0.366858 6 -2 0.382649 13 -7 0.439469 17 19 0.588505 18 -8 0.648837 -10 -11 0.80762 -9 -13 1.03717 1 12 1.16488 -14 -15 1.46078 -12 -16 2.92199 -17 -18

これを下記のPythonのコードで処理すると、以下のようなデンドログラムが作成される。


使い方は、

draw_clusters.py test.out test.jpg

または、

draw_clusters.py test.out test.jpg labels.txt

とする。test.outは入力するクラスタデータ、test.jpgは出力されるデンドログラムの画像ファイルだ。また、二つ目に示した使い方では、labels.txtにデータ数分の文字列を入れておくことで、デンドログラムの数字を任意のラベルにすることができる。因みにこのデンドログラム表示プログラムについても、「集合知プログラミング」を参考にさせてもらった。

draw_clusters.py

#!/usr/bin/env python """ Draw dendrogram for hierarchical clustering by nox, 2009.07.02 """ import sys, os from PIL import Image, ImageDraw HEAD_LENGTH = 50 TAIL_LENGTH = 150 WIDTH = 1200 HEIGHT = 20 class bicluster: def __init__(self, left=None, right=None, distance=0.0, id=None): self.left = left self.right = right self.id = id self.distance = distance def find_cluster(clust, id): for i, c in enumerate(clust): if c.id == id: return i return None def read_clusters(fname, clust, labels): has_labels = False if labels: has_labels = True lines = file(fname).readlines() num = len(lines) + 1 for i in range(num): clust.append(bicluster(id=i)) if not has_labels: labels.append(i) for i, l in enumerate(lines): i = -(i + 1) data = l.strip().split() left_id = find_cluster(clust, int(data[1])) right_id = find_cluster(clust, int(data[2])) clust.append(bicluster(id=i, left=clust[left_id], right=clust[right_id], distance=float(data[0]))) del clust[right_id] del clust[left_id] def get_height(clust): if clust.left == None and clust.right == None: return 1 return get_height(clust.left) + get_height(clust.right) def get_depth(clust): if clust.left == None and clust.right == None: return 0 return max(get_depth(clust.left), get_depth(clust.right)) + clust.distance def draw_dendrogram(clust, labels, image_file): h = get_height(clust) * HEIGHT depth = clust.distance scaling = float(WIDTH - (HEAD_LENGTH + TAIL_LENGTH)) / depth img = Image.new("RGB", (WIDTH, h), (255, 255, 255)) draw = ImageDraw.Draw(img) draw.line((0, h / 2, HEAD_LENGTH, h / 2), fill=(255, 0, 0)) draw_node(draw, clust, HEAD_LENGTH, (h / 2), scaling, labels) img.save(image_file) def draw_node(draw, clust, x, y, scaling, labels): if clust.id < 0: h1 = get_height(clust.left) * HEIGHT h2 = get_height(clust.right) * HEIGHT top = y - (h1 + h2) / 2 bottom = y + (h1 + h2) / 2 ll = clust.distance * scaling lll = ll - clust.left.distance * scaling llr = ll - clust.right.distance * scaling draw.line((x, top + h1 / 2, x, bottom - h2 / 2), fill=(255, 0, 0)) draw.line((x, top + h1 / 2, x + lll, top + h1 / 2), fill=(255, 0, 0)) draw.line((x, bottom - h2 / 2, x + llr, bottom - h2 / 2), fill=(255, 0, 0)) draw_node(draw, clust.left, x + lll, top + h1 / 2, scaling, labels) draw_node(draw, clust.right, x + llr, bottom - h2 / 2, scaling, labels) else: draw.line((x, y, WIDTH - TAIL_LENGTH, y), fill=(255, 0, 0)) draw.text((WIDTH - TAIL_LENGTH + 5, y - 7), str(labels[clust.id]), (0, 0, 0)) def main(args): if len(args) < 2: print >>sys.stderr, "Draw dendrogram for hierarchical clustering." print >>sys.stderr print >>sys.stderr, " Usage: %s cluster_data_file [image_file=clusters.jpg label_data_file=[]]" % os.path.basename(args[0]) sys.exit(1) labels = [] if len(args) > 2: image_file = args[2] else: image_file = "clusters.jpg" if len(args) > 3: for l in file(args[3]): labels.append(l.strip()) clust = [] read_clusters(args[1], clust, labels) draw_dendrogram(clust[0], labels, image_file) if __name__ == "__main__": main(sys.argv)

コード内の変数、WIDTHを変更することで画像ファイルの幅を、HEIGHTを変更することで1データの高さを指定できる。初期値では幅1,200ピクセル、1データの高さを20ピクセルとしている。

次に、任意の閾値でクラスタデータを分けてみる。使い方は、

group_clusters.py test.out test.grp 1.0

のようにする。test.outを入力とし、出力ファイルにtest.grp、閾値を1.0としている。実行後の画面には以下のように表示され、

Input cluster file: test.out Output group file : test.grp Threshold : 1.000000 Number of elements: 20 Number of clusters: 5

出力されるtest.grpは、以下の内容となる。

001: 1 002: 12 003: 6 3 15 13 16 7 11 8 9 004: 2 10 18 4 14 0 5 005: 17 19

ソースコードを以下に示しておく

group_clusters.py

#!/usr/bin/env python """ Grouping clustering data by nox, 2009.07.02 """ import sys, os, math, re def expand_group(clust, n): group = [] for x in clust[n][1:]: if x >= 0: group.append(x) else: group.extend(expand_group(clust, -(x + 1))) return group def arrange_cluster(clust, threshold): tt = [[i, data[1:]] for i, data in enumerate(clust) if data[0] > threshold] t = [] for d in tt: for x in d[1]: if -x < tt[0][0] + 1: t.append(x) groups = [] for x in t: g = [] if x >= 0: g = [x] else: g = expand_group(clust, -(x + 1)) groups.append(g) if not groups: groups = [range(1, len(clust) + 2)] return groups def read_clust(fname, clust_data): for l in file(fname): data = l.strip().split() clust_data.append((float(data[0]), int(data[1]), int(data[2]))) def output_group(fname, groups): f = file(fname, "w") count = 1 for c in groups: f.write("%03d:" % count) for d in c: f.write(" %d" % d) f.write("\n") count += 1 f.close() def main(args): if len(args) < 4: print >>sys.stderr, "Grouping clustering data." print >>sys.stderr print >>sys.stderr, " Usage: %s cluster_data_file output_group_file threshold" % os.path.basename(args[0]) sys.exit(1) clust_file = args[1] group_file = args[2] threshold = float(args[3]) print "Input cluster file: %s" % clust_file print "Output group file : %s" % group_file print "Threshold : %f" % threshold clust_data = [] read_clust(clust_file, clust_data) groups = arrange_cluster(clust_data, threshold) print "Number of elements: %d" % sum([len(x) for x in groups]) print "Number of clusters: %d" % len(groups) output_group(group_file, groups) if __name__ == "__main__": main(sys.argv)

続きを読む...

2009年7月8日水曜日

C++: マルチコアCPUを利用した並列化による高速な階層的クラスタリング

このブログ記事をはてなブックマークに追加

比較的データ数の多い階層的クラスタリングを行う必要があり、手元にあるRで処理を始めたのだが思ったよりも遅かった。そこで、マルチコアCPUを利用した並列化で高速化することにした。RでもGPGPUを使って高速化したプログラムがあるようなのだが、すぐに使える高性能GPUを用意できなかったし、それに、TBBライブラリを使った並列化は手間も時間も掛からないので作ってしまった方が良いと判断した。

尚、このプログラムで作成したクラスタデータのデンドログラム描画や閾値による区分けについては一つの記事で書くには大きすぎるので、記事を分けて、「Python: 階層的クラスタリングのデンドログラム描画と閾値による区分け」に回すことにする。

まずは、「集合知プログラミング」にPythonによる階層的クラスタリングのソースコードが載っていたのでそれをC++で書き直した。アルゴリズム自体はそれほど複雑なものではないのでコーディングは簡単だった。ただ、クラスタ同士の距離を最遠隣法(complete link)にしたかったので、それについては変更を加えた。また、TBBライブラリを利用するために、分割して計算を行う部分については関数オブジェクトで実装した。

余談だが、この「集合知プログラミング」はかなり良い書籍だと思う。Pythonを使っているのでコードも読みやすい。興味のある方は是非読んでみて欲しい。

完成したソースコードは最後に示しておいた。赤で示した部分がユークリッド距離を計算する関数であり、この部分を変更することで、べき乗距離にすることも、マンハッタン距離にすることも、別のどんな指標にすることも、簡単にできる。因みに自分はRMSDの2乗で計算を行った。また、青で示した部分が最遠隣法の計算部分であり、ここを書き換えることで群平均法や重心法などの他の結合方法に変更することができる。

入出力データについてだが、任意の次元数の1データをスペース区切りの一行とし、複数データを格納したファイルを入力ファイルとしている。そして、結合する距離、データの番号のペアを一行としたデータが出力ファイルとして書き出される。番号は入力データを正の整数(0~)として、データの集合体であるクラスタを負の整数(-1~)として扱っている。

それでは実際に計算してみよう。まず、以下のような入力ファイルを用意する。この例では次元数3でデータ数は20となる。

test.dat

3.790 7.255 10.259 3.848 7.352 9.208 3.869 7.320 9.771 3.937 7.418 10.621 3.848 7.050 10.359 3.840 7.232 10.370 3.944 7.509 10.918 3.833 7.168 10.598 3.908 7.260 10.822 3.904 7.112 10.834 3.786 7.336 10.060 3.856 7.155 10.722 3.723 6.424 9.654 3.796 6.936 10.652 3.843 7.069 10.440 3.829 7.417 10.582 3.863 7.289 10.669 3.905 6.078 8.719 3.830 6.668 10.246 3.811 6.202 8.308

これを、以下のように実行する。ここでの実行ファイルは clusters_tbb としている。

clusters_tbb test.dat test.out

出力ファイル test.out の内容は以下の通り。

test.out

0.0833487 4 14 0.11483 3 15 0.123895 0 5 0.126783 7 11 0.144271 16 -4 0.14854 8 9 0.253505 -5 -6 0.264889 -1 -3 0.301108 2 10 0.366858 6 -2 0.382649 13 -7 0.439469 17 19 0.588505 18 -8 0.648837 -10 -11 0.80762 -9 -13 1.03717 1 12 1.16488 -14 -15 1.46078 -12 -16 2.92199 -17 -18

一番最初の行は、距離0.11483で3番(先頭は0番)と15番のデータが結合してクラスタを作っていることを示し、最後の行では、距離2.92199で17番(先頭は1番、-1と表示される)と18番のクラスタが結合していることを示している。

最後にRとの速度比較を行った。入力として次元数45の2000データを使用した。また、Rについては、Rscriptを利用して以下のように実行している。

Rscript test.R

test.R

data = read.table("test.dat") d = dist(data) hc <- hclust(d, method="complete") write.table(hc["merge"], "merge.out", row.names=F, col.names=F) write.table(hc["height"], "height.out", row.names=F, col.names=F)

結果は以下の通り。

Fedora 8 / Intel Core2 Duo E8400@3.00GHz x 1 CPU / 2 GB R 2.6.1 Rscript / gcc 4.1.2, tbb 2.1 : g++ clusters_tbb.cpp -O3 -ltbb R 17.46秒 1 core 13.02秒 2 cores 6.91秒

Red Hat Enterprise Linux 5 / Quad-Core AMD Opteron 8356 x 8 CPU / 196 GB R 2.9.0 Rscript / icpc 11.0, tbb 2.1 : icpc clusters_tbb.cpp -O3 -ltbb R 36.11秒 1 core 27.25秒 2 cores 14.58秒 4 cores 7.45秒 8 cores 4.50秒

Red Hat Enterprise Linux 5 / Intel Xeon E5450@3.00GHz x 2 CPU / 8 GB icpc 11.0, tbb 2.1 : icpc clusters_tbb.cpp -O3 -ltbb 1 core 11.16秒 2 cores 5.45秒 4 cores 3.25秒 8 cores 2.42秒

Microsoft Windows XP SP3 / Intel Core2 6600@2.40GHz x 1 CPU / 2 GB R 2.9.1 Rscript / cl 15.00 (Visual Studio 2008), tbb 2.1 : cl clusters_tbb.cpp /EHsc /Ox /MD R 23.65秒 1 core 43.53秒 2 cores 23.14秒

Linuxでは1コアでの計算でも今回作成したプログラムの方がRの階層的クラスタリングよりも若干速い結果を出している。また、並列化に関しては、2コアでは1コアのほぼ倍の性能を出しており、SMP環境での実行速度を見ると、4コアで3.6倍以上、8コアでも6倍以上で、十分な性能が出ていると思う。ただ、OpteronのSMP環境ではRでも自作プログラムでも遅いが何故だろう。あと、Windows環境では1コアでRの方が倍ぐらい速く、2コアでトントンだった。メモリ関連の実装で効率が違うのだろうか。

Python: 階層的クラスタリングのデンドログラム描画と閾値による区分け」に続く。

clusters_tbb.cpp

// hierarchical clustering using the TBB library // by nox, 2009.07.02 // // Linux - Intel C++ // icpc clusters_tbb.cpp -O3 -ltbb -o clusters_tbb // // Linux - GCC // g++ clusters_tbb.cpp -O3 -ltbb -o clusters_tbb // // Windows - Microsoft Visual Studio 2008 // cl clusters_tbb.cpp /EHsc /Ox /MD #include <iostream> #include <fstream> #include <sstream> #include <vector> #include <algorithm> #include <utility> #include <cmath> #include "tbb/task_scheduler_init.h" #include "tbb/blocked_range.h" #include "tbb/parallel_for.h" #include "tbb/parallel_reduce.h" #include "tbb/parallel_sort.h" using namespace std; using namespace tbb; double calc_distance(const vector<double>& vec1, const vector<double>& vec2) { vector<double>::const_iterator p, q; double s = 0.0; for (p = vec1.begin(), q = vec2.begin(); p != vec1.end(); ++p, ++q) s += (*p - *q) * (*p - *q); return sqrt(s); } static vector<int> table_i; struct bicluster { bicluster* left; bicluster* right; int id; vector<int> vec_ids; double distance; bicluster(bicluster* left_, bicluster* right_, int id_, vector<int> vec_ids_, double distance_) : left(left_), right(right_), id(id_), vec_ids(vec_ids_), distance(distance_) { } }; typedef vector<bicluster*>::iterator vec_t; typedef pair<vec_t, vec_t> pair_t; typedef pair<int, int> distance_pair_t; typedef vector<double> distances_t; class CompleteLink { public: double dist; const vector<int>& clust; const distances_t& distances; CompleteLink(double d_, const vector<int>& clust_, const distances_t& distances_) : dist(d_), clust(clust_), distances(distances_) { } CompleteLink(const CompleteLink& cl, split) : dist(cl.dist), clust(cl.clust), distances(cl.distances) { } void operator()(const blocked_range<vector<int>::const_iterator>& r) { for (vector<int>::const_iterator p = r.begin(); p != r.end(); ++p) { for (vector<int>::const_iterator q = clust.begin(); q != clust.end(); ++q) { int i = min(*p, *q); int j = max(*p, *q); dist = max(distances[table_i[i] + j], dist); } } } void join(const CompleteLink& cl) { dist = max(dist, cl.dist); } }; double calc_cluster_distance(const vector<int>& clust_i, const vector<int>& clust_j, const distances_t& distances) { CompleteLink cl(0.0, clust_j, distances); parallel_reduce(blocked_range<vector<int>::const_iterator>(clust_i.begin(), clust_i.end(), 100), cl); return cl.dist; } class ReduceHCluster { public: const vector<vector<double> >& data; const vector<bicluster*>& clusters; distances_t& distances; pair_t lowest_pair; double closest; ReduceHCluster(const vector<vector<double> >& data_, const vector<bicluster*>& clusters_, distances_t& distances_, pair_t lowest_pair_, double closest_) : data(data_), clusters(clusters_), distances(distances_), lowest_pair(lowest_pair_), closest(closest_) { } ReduceHCluster(const ReduceHCluster& hc, split) : data(hc.data), clusters(hc.clusters), distances(hc.distances), lowest_pair(hc.lowest_pair), closest(hc.closest) { } void operator()(const blocked_range<vec_t>& r) { for (vec_t p = r.begin(); p != r.end(); ++p) { for (vec_t q = p + 1; q != clusters.end(); ++q) { int i = min((*p)->id, (*q)->id); int j = max((*p)->id, (*q)->id); if (distances[table_i[i] + j] < 0) distances[table_i[i] + j] = calc_cluster_distance((*p)->vec_ids, (*q)->vec_ids, distances); double d = distances[table_i[i] + j]; if (d < closest) { closest = d; lowest_pair = make_pair(p, q); } } } } void join(const ReduceHCluster& hc) { if (closest > hc.closest) { closest = hc.closest; lowest_pair = hc.lowest_pair; } } }; class HClusterDistances { private: const vector<vector<double> >& data; int n; distances_t& distances; public: HClusterDistances(const vector<vector<double> >& data_, int n_, distances_t& distances_) : data(data_), n(n_), distances(distances_) { } void operator()(const blocked_range<int>& r) const { for (int i = r.begin(); i != r.end(); i++) for (int j = i + 1; j < n; j++) distances[table_i[i] + j] = calc_distance(data[i], data[j]); } }; inline int index_i(int n, int max_n) { int s = 0; for (int i = 0; i < n; i++) s += max_n - i; return s - n - 1; } void hcluster(const vector<vector<double> >& data, vector<bicluster*>& clusters) { int n = clusters.size() * 2; for (int i = 0; i < n; i++) table_i.push_back(index_i(i, n - 1)); // distance between i and j: distances[table_i(i) + j] distances_t distances(n * (n - 1) / 2, -1.0); int current_clust_id = clusters.size(); parallel_for(blocked_range<int>(0, clusters.size() - 1, 100), HClusterDistances(data, clusters.size(), distances)); while (clusters.size() > 1) { cerr << clusters.size() << endl; pair_t lowest_pair = make_pair(clusters.begin(), clusters.begin() + 1); double closest = calc_cluster_distance((*lowest_pair.first)->vec_ids, (*lowest_pair.second)->vec_ids, distances); ReduceHCluster hc(data, clusters, distances, lowest_pair, closest); parallel_reduce(blocked_range<vec_t>(clusters.begin(), clusters.end(), 50), hc); closest = hc.closest; lowest_pair = hc.lowest_pair; vector<int> vec_ids((*lowest_pair.first)->vec_ids); vec_ids.insert(vec_ids.end(), (*lowest_pair.second)->vec_ids.begin(), (*lowest_pair.second)->vec_ids.end()); vector<int>().swap((*lowest_pair.second)->vec_ids); vector<int>().swap((*lowest_pair.first)->vec_ids); bicluster* new_cluster = new bicluster(*lowest_pair.first, *lowest_pair.second, current_clust_id, vec_ids, closest); current_clust_id++; clusters.erase(lowest_pair.second); clusters.erase(lowest_pair.first); clusters.push_back(new_cluster); } } void read_data(const string fname, vector<vector<double> >& data) { fstream fs(fname.c_str(), ios_base::in); string line; stringstream ss; vector<double> v; double value; while (getline(fs, line)) { for (ss.str(line), v.clear(); ss >> value; v.push_back(value)); ss.clear(); data.push_back(v); } } void store_clusters(const bicluster* cluster, vector<pair<double, distance_pair_t> >& cluster_data, int max_n) { static int n = 0; if (cluster->left && cluster->right) cluster_data.push_back(make_pair(cluster->distance, make_pair( cluster->left->id < max_n ? cluster->left->id : -(cluster->left->id - max_n + 1), cluster->right->id < max_n ? cluster->right->id : -(cluster->right->id - max_n + 1)))); if (cluster->left) store_clusters(cluster->left, cluster_data, max_n); if (cluster->right) store_clusters(cluster->right, cluster_data, max_n); } void write_clusters(const string fname, const vector<pair<double, distance_pair_t> >& cluster_data) { fstream fs(fname.c_str(), ios_base::out | ios_base::trunc); for (vector<pair<double, distance_pair_t> >::const_iterator p = cluster_data.begin(); p != cluster_data.end(); ++p) { fs << p->first << " " << p->second.first << " " << p->second.second << endl; } } int main(int argc, char** argv) { if (argc < 3) { cerr << "Hierarchical clustering using the TBB library." << endl; cerr << endl; cerr << " Usage: " << argv[0] << " input_data_file output_cluster_file" << endl; exit(1); } task_scheduler_init init; vector<vector<double> > data; vector<bicluster*> clusters; read_data(argv[1], data); for (int i = 0; i < data.size(); i++) clusters.push_back(new bicluster(NULL, NULL, i, vector<int>(1, i), 0.0)); hcluster(data, clusters); vector<pair<double, distance_pair_t> > cluster_data; store_clusters(clusters[0], cluster_data, data.size()); parallel_sort(cluster_data.begin(), cluster_data.end()); write_clusters(argv[2], cluster_data); return 0; }

続きを読む...