勾配方向画像のHSV表示スクリプト

画像の勾配方向をHSVに変換するPythonスクリプトです．
OpenCV2.3.1で動作します．2.2系では動かないと思います． 

使い方は，
-i オプションで勾配方向画像を得たい画像ファイルを与え
-o オプションで変換後の画像ファイルを指定します．
-s オプションをつけると変換後に表示してくれます．

 使用例）
 python gradient_direction_hsv.py -i source.jpg -o result.jpg -s

 変換例）

ダウンロード

NIPALSとかPLSとか

昨今の認識アルゴリズムは，マシンパワーにものを言わせて判別に使う情報量をどんどんと増やしていく傾向にあると思うのです． 例えばHOGとかがいい例ではないでしょうか．私の5年前のノートPCでは遅すぎて使いもになりません．

そもそも自分が学習機ならとても覚えられないような情報量を無理やり食わせて，性能をだすことができるんでしょうか． もっとシンプルにやる方法が必ずあるに違いないと心の奥底で私は信じているのですが，今は具体案を提案できている状態ではないので私の主張には説得力はありません．

実際に人間の脳みそは，体重比率で大きなウエイトを占めているわけで，生物の知能が発達すると共に脳みそは大きくなってきたのを鑑みると，高度な情報を処理するためには複雑な仕組みが必要なのだと，暗に示唆されているようにも思えます．高度な知性を実現するには脳みそのような複雑なハードが必要なのだという証だといわれれば反論するのが難しいです．

だけれどもほんの一部でもいいから，シンプルな手法によって人間の知性を計算機の上で実行できないものだろうかと思うのです． 映像から人がどこを歩いているのかだけでも人並みに判別できれば，悲惨な交通事故をもっと減らすことのできる役立つ機械を作ることが出来るはずです．
もしシンプルに計算量の少ないアルゴリズムでそれが実現できれば，安いハードウエアで製品を構成することができるようになり，普及が期待されます．

そこで手っ取り早い方法として性能は保証されているけれど膨大な情報を収集する必要のあるアルゴリズムを使いつつ，その中からあまり重要ではない情報をふるいにかけて情報量を減らして，計算量を減らすアプローチを考えます．このようなことをやる代表的な手法として主成分分析(PCA)があります．PCA-SIFTなどが応用例としていい例でしょう．

でも実際にPCAを分類器に入力するデータの次元圧縮手段として使ってみると，期待したほど性能は良くないんじゃないかという印象を持っています． もちろんそれはケースバイケースなんだと思いますが，画像認識に使うようなデータでは判別結果に使うパラメータの線形独立性を期待できないことがほとんではないでしょうか．

そのようないい加減なデータでもパワフルに働いてくれる次元圧縮方法がないものかと探していたら，先日PLSという手法(Partial Least Square)を見つけました．化学分析分野でよく用いられる方法のようで，例えば，サンプルの各波長における吸光度から，どのような物質がどのような割合で含まれているかなどを分析するために開発された手法のようです．

日本語で書かれたフリーで読める文献がありました．こちらです「PLS 回帰におけるモデル選択」．でも日本語で読んでもなんだかよく分かんなかったです．ちょっと苦しいけれど英語の文献の方が丁寧に書かれていて，理解しやすかったりすることも多いです．これもいいかも知れません「A Beginner’s Guide to Partial Least Squares Analysis」
実は自分が参考にした最もわかりやすかった文献のリンクが今探せないです．どこいっちゃったものだか．あとで見つけたら追記しておきます．

このPLSをお手軽にためすには，「R」のplsパッケージを使うのが良いみたいで，PLSの紹介論文などでは，Rに付属のgasolineという吸光度とオクタン価の関係を測定したデータが良く紹介さrています．Rはあんまり好きじゃないので，早速，ガソリンデータを使ってPLSアルゴを試すpythonコードを書いてみました．

def NIPALS_internal(X, Y, n, M, N, u, epsilon):
    u0 = u
    while True:
        w = X.T * u / (u.T * u)
        w = w / norm(w)
        t = X * w
        c_ = Y.T * t / (t.T * t)
        c = c_ / norm(c_)
        u = Y * c
        if norm(u0 - u) < epsilon: break
        u0 = u
    p = X.T * t / (t.T * t)
    q = Y.T * u / (u.T * u)
    Xnew = X - t * p.T
    Ynew = Y - t * c_.T
    return Xnew, Ynew, t, u, p, q, w

def NIPALS(X, Y, k = 1, epsilon = 1e-12):
    n, N = X.shape
    M = Y.shape[1]
    X_mean = X.mean(axis=0)
    X_std = X.std(axis=0)
    Y_mean = Y.mean(axis=0)
    Y_std = Y.std(axis=0)
    X_ = (X - X_mean) / X_std
    Y_ = (Y - Y_mean) / Y_std
    X__, Y__ = X_, Y_
    W = matrix(zeros( (N, k), float32))
    T = matrix(zeros( (n, k), float32))
    U = matrix(zeros( (n, k), float32))
    P = matrix(zeros( (N, k), float32))
    Q = matrix(zeros( (M, k), float32))
    for i in range(k):
        u = Y_[:,0]
        Xnew, Ynew, t, u, p, q, w = NIPALS_internal(X__, Y__, n, M, N, u, epsilon)
        X__, Y__ = Xnew, Ynew
        W[:,i], T[:,i], U[:,i], P[:,i],Q[:,i] = w,t,u,p,q
    B = X_.T * U * (T.T * X_ * X_.T * U).I * T.T * Y_
    return (B, X_mean, X_std, Y_mean, Y_std)

Rのplsに付属のガソリンデータの読み込みと，既述のNIPALS(), NIPALS_predict()の利用シーンは次のような感じです.

def load_gasoline():
    f = open("gasoline.txt")
    reader = csv.reader(f, delimiter=' ', quoting=csv.QUOTE_NONE)
    octane, NIR = [], []
    for i, row in enumerate(reader):
        if i > 0:
            octane.append(int(row[0].strip('\"')))
            NIR.append([float(j) for j in row[1:]])
    return matrix(octane).T, matrix(NIR)

if __name__ = '__main__':
    Y, X = load_gasoline()
    tt = []
    t_y0, t_y1 = [],[]
    for h in range(2,26):
        s = 0
        for i in range(60):
            p = NIPALS(X, Y, h)
            Y_ = NIPALS_predict(p, X[i,:])
            t_y0.append(Y[i,0])
            t_y1.append(Y_[0,0])
            s += abs(Y[i, 0] - Y_[0,0])
        tt.append(s/60)

    gp = Gnuplot.Gnuplot(debug = 1)
    gp.xlabel('X')
    gp.ylabel('Y')
    gp('set grid')
    s = Gnuplot.Data(range(len(tt)), tt, title='error',with_='points 3 3')
    s1 = Gnuplot.Data(range(len(t_y0)), t_y0, title='raw',with_='line')
    s2 = Gnuplot.Data(range(len(t_y1)), t_y1, title='predict',with_='line')
    gp.plot(s, s1, s2)
    raw_input()

これを足がかりにして認識アルゴリズムに応用して性能を試すようなことを正月休みにやってみたいと思っていたけど，子供たちの喧騒の中でブログに書くのが精一杯でした．

休みも残すところ今日の半日と明日1日になってしまいました. そろそろ家出ゴロゴロしているのがだるくなってきたので，これからどこかへ出かけてみたいと思っていますが，人ごみは嫌いであり，はて　どこへ行こうか，子供たちは雪遊びできるところへ行けば喜ぶに違いありませんが，心の底から雪遊びを楽しむ心を忘れてしまった私は，寒いばかりなので他にいいところは無いものかと考えるのですが，人ごみも苦手だし，あてもなく困ったものです．

漠然とした不安とmilkモジュール

皆様おはようございます．
今朝の滝沢村は，雨だれの音が音がこだましていて
どんよりとした雲が低く垂れこんでおります．

地震から日が経つに連れ余震の頻度は落ちてきた
はずなのに，正直なところ気持ちは落ち込む一方です．
なにか悪いことが起こるのではないかという，
なんともいえない気持ち悪さが心から離れません．
そんなことは無いさと思ってみても，気を緩める頃には
決まって地震に大地が震え，私は再び悪い予感に
とらわれるという循環になっています．

自分の心に素直に向き合い，本能の声に耳を傾けるならば
「これだけでは済まされない，もう一つ何かが起こるよ」
というささやきが聞こえてくるようです．
その声を無視して日常をやり過ごせばよいのか，それとも
理屈では説明のつかない悪い予感に従えばいいのか，
まだ割り切れずにいながら，時間が過ぎていきます．

ところで，ここまでの話と全く脈略がありませんが，
ここからガラリと話題を変えまして．．．

先日，少し知的な処理をコンピュータにやらせたいと思いまして，
pythonの機械学習モジュールについて調べてみる機会がありました．
そのとき見つけた「milk」というモジュールについてメモを披露したいと思います．

「milk」は名前の印象とは関係なくって複数の機械学習
アルゴリズムを集めたpythonパッケージのようです．
SVMやboostingなどをお手軽に使えるところが気に入りました．
複数のアルゴリズムを組み合わせて使うことも考慮されています．
残念なことに，milkのドキュメントは２つくらいしか見つかりません．
pypiとluispedroです．
インストールは「easy_install milk」でいけました．

使い方ですが，サンプルを引用します．

import numpy as np
import milk
features = np.random.rand(100,10)
labels = np.zeros(100)
features[50:] += .5
labels[50:] = 1
learner = milk.defaultclassifier()
model = learner.train(features, labels)
example = np.random.rand(10)
print model.apply(example)
example2 = np.random.rand(10)
example2 += .5
print model.apply(example2)

milk.defaultclassifier()は，お手軽に使えるように
あらかじめデフォルト設定された学習機となっていて
SVM=Support Vector Machineを使った学習となっています．
教師信号として入力する特徴量(features)と対応するクラスタラベル（label）を
入力してやります．
特徴量の正規化や線形従属な要素の除去なども全自動でやってくれます．
defaultclassifier()では複数のクラスタに分類するためにone-versus-rest手法が
用いられるようになっています．
もちろんこれらの組み合わせを自分で自由に設定することも可能です．
defaultclassifier()の中身は，次のようになっているようです．

defaultclassifier = ctransforms(
   chkfinite(),
   interval_normalise(),
   featureselector(linear_independent_features),
   sda_filter(),
   gridsearch(one_against_one(svm.svm_to_binary(svm.svm_raw())),
   params={
     'C': 2.**np.arange(-9,5),
     'kernel': [svm.rbf_kernel(2.**i) for i in np.arange(-7,4)],
   }
))

ctransforms()というのが複数のフィルタ（学習器）を束ねてくれるmilkのAPIです．
そのなかで１つづつ与えられている引数の要素がそれぞれ別個の学習機や
フィルタとなっています．
・chkfiite()は,特徴量やラベルに無限やNaNなどが含まれていないかチェックするフィルタ．
・interval_normalise()は，特徴量を正規化するフィルタ．
・sda_filter()は，判別分析を行うフィルタ．
・gridsearch()は，多値問題をグリッド検索に基づいて与えられた学習器で学習するフィルタ．
・svm.rbf_kernel()は，Support Vector Machineのradial basis function kernel．

milkのフィルタ（学習器）は，皆 learn_func(features, label)のような引数を持つように
設計されていて，ctransforms()で束ねて与えると，パイプライン処理により１つの
学習機として動作するようになっているようです．

冒頭のサンプルコードに戻りますが，学習器のインスタンスを生成した後，
「train」で特徴量とラベルを与えて学習させます．

learner = milk.defaultclassifier()
model = learner.train(features, labels)

処理には数秒を要します．学習が終わると学習モデルが得られますので，
未知の特徴量を与えれば，学習結果に基づいて分類してくれます．

print model.apply(example)

学習モデルは，pickle化できるのでディスクから読み込むことも可能です．

OpenCVと組み合わせたり，あるいはネットから集めたデータを
自動分類したい時など，お手軽に使えるのではないでしょうか．

近況報告と64bit windows環境でのCython

皆様こんにちは．
どうやら２ヶ月ぶりの更新だったようです．
２月から３月までは，ただブログ更新をサボっていただけなのですが，
３月から４月にかけての１ヶ月間は　リアルな生活で精一杯でした．
ようやく余震は，頻度も規模も小さくなってきて3.11以前の日常に戻りつつあります．
大震災で亡くなられた皆様のご無念を思うと切なくてなりません．
亡くなった皆様のご冥福をお祈り致しております．
また被災者の皆様には一日も早く落ち着いた生活を取り戻せますように祈っております．

日本政府の遅々とした対応には憤りを感じるばかりです．
奮闘している現場と対照的に政府にはリーダーシップのかけらも感じられません．
今回の大震災は被災地が広範囲に渡っているとはいえ散々な対応ぶりに見えます．
震災以前は，日本政府の能力は高く国民主体であると漠然と思って疑いませんでしたが，
これが大きな勘違いであったと今回の出来事で思いました．
普段あまり意見の合わないカミさんとさえも，この点においては意気投合しておりましたので，
世の中の皆さんにもこのように感じた方が少なくないのではないかと
勝手ながら想像しているところです．
どうか日本に一刻も早く良きリーダーが与えられて人々が希望を胸に
前に進むことができるようにと強く願い，毎日祈る日々です．

　私の近況ですが．．地震の後　３日ほど停電，
その後の余震でも１日程度停電があった事などで
仕事にも多少の影響はでておりましたが，
事前の計画を遅れなく実現すべく開発業務に勤しんでいます．
　この頃は，もっぱら物理現象を観察しその現象をある程度再現する
なるべく簡略なモデル式とパラメータを探すといった地味な仕事をしております．
地味な仕事とはいえ　このような最適化の作業は，物理の教科書をほっくりかえしつつ，
過去にサボったツケを実感しながら，
試行錯誤の繰り返しと計算の待ち時間のやり過ごし方など
ストレスの多い作業となっています．
　数時間かかった最適化計算が収束せず，良く調べてみたら計算プログラムのバグに
行き着いたときに感じる無気力感を　気持ちの隅に追いやりながら
作業を継続せねばなりません．
　以前は最適化計算をさせるためにC++やCを使っていましたが，
この頃は　Python+SciPy+Numpyの組み合わせが気にいております．
特に「scipy.optimizeパッケージ」に感謝です．
　私は普段，Intel T7600(Core2duo 2.3GHz)＋32bitVistaのノートPCを
持ち歩いて仕事しているのですが，最適化の計算だけは時間がかかってしまうので，
その時だけちょっと速そうなデスクトップパソコンを起動して計算させています．
例えばAMD PhenomII3GHz＋Vista64bitなどにPythonスクリプトを走らせると，
T7600ノートPCでは及びもつかないスピードで計算してくれます．
当初は，私もこれだけで満足していたのですが，どんどんフィッティング対象の式が
複雑怪奇なものに成長するにしたがってPhenomマシンでも待ち時間が長くなり
どうしたものかと悩むようになりました．
　JITで早いとウワサのPyPyを試してみようと思いましたが，残念なことに
NumpyもScipyもPyPyでは動きませんので箸にも棒にもかかりません．
　そこで，Cythonを試してみることにしました．
Cythonは，Pythonもどきのスクリプトを書いてあげるとCに変換してくれて
pythonからimport可能なDLLを生成してくれるものですが，
現時点におけるpython高速化の確実なソリューションといえると思います．
　しかし普段やらないことにとりかかるのは億劫なものです．
ドキュメントも読まねばなりません．
こういうパッケージを追加していざ取り組んでみると落とし穴にハマって
抜け出せないことが多々あります．Cythonも例外ではありませんでした．

　結論からいうと，Vista32bit環境でCythonを使うことは全く問題ありませんでした．
Vista32bitではVisual Studio2005をインストールしているけれども，
mingwでコンパイルする設定を追加しておりますが問題もなく動きました．
きっとVCでもサクっとコンパイル出来るに違いありません．
私のVista32bit環境python27にあるLIB/distutils/distutils.cfgの中身を書いておきます．

[build]
compiler = mingw32

Cythonスクリプトをコンパイルしたときのsetup.pyは次のとおりです．

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
ext_modules = [Extension("hoge_func", ["hoge_func.pyx"], include_dirs=["C:/cygwin/usr/local/Python27/Lib/site-packages/numpy/core/include"])]
setup(
  name = 'hoge function',
  cmdclass = {'build_ext': build_ext},
  ext_modules = ext_modules
)

これでhoge_func.pyxにCythonスクリプトを書いておいて
"python setup.py build_ext --inplace"を実行すればサクっとコンパイルされて
pythonから"import hoge_func"できて恩恵に預かることができました．
ただしVista32bit環境だけです．
　Vista64bit環境のCython設定でハマりました．
正しく言えばdistutilの設定ということでしょうか．
Vista32環境＋Cythonで十分早くって64bitで動かすことをしなくてもいいぐらい
高速化されたんですが，欲が出まして Vista64でも scipy+numpy+cythonの
最適化計算環境がどうしても欲しくなりました．
そこで 64bitネイティブなmingw64でもないかと探して"x86_64-w64-mingw32-gcc"なるものを
見つけインストールしsetup.pyを走らせましたがエラーを吐きます．
エラーの内容はもはや思い出せませんが，多数のエラーにブチ当たりました．
行き着いた先が，「Compiling 64-bit extension modules on Windows」というドキュメント．
64bitでcythonを楽しむためには次の条件を揃えねばならないようです．
1.Python2.6か2.7もしくは3.1
2.Microsoft Windows SDK for Windows 7 and .NET Framework 3.5 SP1が必須
　つまり，mingw64は使えないということです．
先程のドキュメントにはSDKバージョンは"Microsoft Windows SDK for Windows 7 and .NET Framework 3.5 SP1"以降となっていますが，最新版をインストールしてもコンパイルできませんでした．
SDKのダウンロードも結構かかった上に２種類もインストールしたり，
インストール中にフリーズしたりとハマリにハマってやっとコンパイルできました．
コンパイル時の呪文も書いておきます．
これをバッチファイルやMakefileなどに書いておくといいかもしれません

C:\Program Files\Microsoft SDKs\Windows\v7.0>set DISTUTILS_USE_SDK=1
C:\Program Files\Microsoft SDKs\Windows\v7.0>setenv /x64 /release
python setup.py build_ext --inplace

　苦労した割に，振り返ってみると大した事ではなかったのに
いったい何でハマってしまうのかと毎回思うのですが，
知らないもんは知らないのですからしょうがないと諦めます．
　ちなみにCythonの高速化は，ハンパではありませんでした．
正確にベンチマークとったら報告したいと思いますが，
体感スピードでは100倍以上早くなったといっても大げさではないと思っています．

FT2232HをPythonで使う

ストロベリーリナックスからFT2232Hを載せたUSB高速シリアル変換基板が発売されております．
FT2232H(2ch)高速USBシリアル変換モジュールキット
これは単純に言えば，パソコンから自作マイコンやFPGAボードなどに480Mbpsでデーターを簡単に送れてしまえるスグレモノなのです．
パソコン側のソフトを開発するためのライブラリなどはFTDI様が配っておられます．
FTDI CHip D2XX Driver Download site
筆者はこの頃，年をとったせいか開発中にコンパイラを起動するのがストレスになり，
単にアウトプットだけが問題なときやテストメインの時は，Pythonでお手軽に処理するのが最近の趣向なのです．
そこで上記のボードもPythonから使えるようにFTDI様のDLL用にctypesラッパーモジュールを書きました．
これを書いているほうが断然キーボードを叩く量が多い気が途中でしましたが，一度作れば480Mbpsの通信を
Pythonで制御できるので嬉しい気分となります．
ftd2xx-py.zip
スクリプトの中身は整理もなんにもしてなくて，コメントもありませんがご容赦下さい．

使い方は，例えばこんな感じです．３２メガバイトを送信する例となっています．

from ftd2xx import *
handle = c_ulong()
ret=FT_Open(0, byref(handle))
ret=FT_SetBitMode(handle, c_ubyte(0xff), c_ubyte(0x0))
print "FT_SetBitmode=", ret
ret=FT_SetBaudRate(handle, 6000000)
print "FT_SetBaudRate=", ret
data = range(1024)*1024*32
ret=easy_write(handle, data)

Pythonで非線形最適化計算する

私自身がPythonを便利に思って重点的に利用させてもらっている技術分野は，科学技術計算とUSBやシリアルなどとの通信に関するもの．

特にお世話になっている最適化計算についてちょっとご紹介．

最適化で使っているパッケージは, scipy.
他に Numericが必要.

excel等で近似多項式を求める事が出来るけど，奇数次数のみの多項式に
フィットさせるなどは不可能．

そういうときにPythonの scipyは便利．
次のスクリプトは９次項までの奇数次数のみの多項式へのフィッティングを行うもの．
私は魚眼カメラの射影式を扱うとき重宝している．

import scipy
import scipy.optimize
import Numeric as n

def func(p, x):
    k1,k2,k3,k4,k5 = p
    p0=[k5,0.,k4,0.,k3,0.,k2,0.,k1,0.]
    return scipy.polyval(p0, x)

def residue(p,y,x,sigma):
    #print "residual_func ", y, x
    err = y - func(p, x)
    #print "error ", err
    return err

def fitting(x, y):
    sigma = 0
    r = scipy.optimize.leastsq(
        residue,
        [1.0,0,0,0,0],args=(y, x, sigma), full_output=1, ftol=1e-32, xtol=1e-32, maxfev=10000)
    return r[0].tolist()

if __name__== '__main__':
    sample = n.array([
    0.16581,0.128181366,
    0.33161,0.258930716,
    0.49742,0.393306338,
    0.66323,0.530482793,
    0.82903,0.666571512,
    0.99484,0.793397256,
    1.16064,0.899001855,
    1.32645,0.970166103,
    1.49226,0.99910374,
    1.65806,0.991575672])

    xy = n.reshape(sample, (10,2))
    x = xy[:,0]
    y = xy[:,1]
    optimized_parameter = fitting(x, y)

    print optimized_parameter

scipy.optimize.leastsqは minpackを移植したもののようである．
Levenberg-Marquardt法を用いた最適化をお手軽に行うにはもってこいだと思う．