2014年11月25日火曜日

DICOM



今回はかなり専門的な話です。

DICOMというものをご存知でしょうか。

たぶん、ほとんどの人は見る機会さえないので知らないと思いますが、医療用に用いられる画像形式(のようなもの)です。


例えば、レントゲンとかCT、MRIとかです。

昔はフィルムで撮影したものも技術の発展とともにデジタル化され、

今では撮影・撮像機械と病院内のコンピュータシステムが連携しています。

最近は都市部の病院は小さな開業医でも電子カルテを使っていたりして、急速に電子化が進んでいます。


DICOM形式はアメリカで開発されたもので、今では世界的な標準仕様になっています。

なぜ、普通のJPEGやPNGではだめなのかというと、医療用の画像には画像ピクセルそのものの情報だけではなく、撮影した状況、部位、患者情報、日付などが記されています。


仕事の関係でDICOMのヘッダ部分を分析する必要があったため
備忘録として簡単な導入部分をまとめてみました。
pythonはライブラリが豊富なので、DICOMライブラリも読み込むだけで使えます。


その前に、まずPIPでPydicomを インストールします。



$ pip install pydicom

例: 引数で与えるファイル名のDICOMファ イルを読み込み、患者IDなどの情報を表示するスクリプト


ファイル名:showDicomInfo.py


#!/usr/bin/env python

import sys,dicom

d = dicom.read_file( sys.argv[1] )
print d


実行例:


$ ./showDicomInfo.py example.dcm


患者IDのみを表示する場合、

#!/usr/bin/env python

import sys,dicom

d = dicom.read_file( sys.argv[1] )
print d.PatientID



ほかの値についてはdirメ ソッドで見てください。
つまり、上の例でいうと

print dir(d)

とすると
メソッド・メンバ変数の一覧が表示されます。


2014年11月10日月曜日

JSON文法チェック方法

データの送受信または記述のためよく使われるJSON(ジェイソン)という形式があります。
テキストなので、windowsなら標準で入っているメモ帳のようなテキストエディタで編集できます。

例えば、
{
   "name" : ["shimotomai", "takayuki"],
   "degree":"PhD"
}
こんな感じです。

入れ子構造にもでき、XML同様柔軟性が高い記述方式として最近よく使われています。

ただ、コンマ忘れやすいので、よく文法エラーになります。
もしくは最後の要素の後にコンマを入れたり、
コロンがなかったり、ダブルクォーテーションがシングルになっていたり。

ロボットの音声認識システムで、採用していたのでJSON形式よく利用していましたが、

記述が長いと、目で見てエラーをチェックするとか拷問に近いです。
あり得ません。

ということで、

自動的にチェックしてくれる方法をいくつか利用していました。
そのうちPythonを使って行う方法をご紹介します。
簡単です。

$ python -mjson.tool dummy.json



これだけです。
今回はプログラムさえ書きません。

エラーの場合

Expecting property name: line 12 column 2...
という風に行番号が表示されます。つまりこの場合は、12行目にエラーがありますよということです。


demjsonがインストールされていれば
jsonlint -v hoge.json
でもチェック可能です。

demjsonのインストールはubuntuの場合、
sudo aptitude install python-demjson
でできます。


pythonを使ったjson形式の取り扱いについてはそのうちまた詳しく書きます。




2014年11月8日土曜日

ラン ダム文字列生成

ご無沙汰しております。
最近仕事で、任意の長さのランダム文字列を生成するスクリプトを作りました。

CGIなどでセッション情報などをpythonで管理する場合、重複する可能性の低いランダムな文字列が必要になることがあります。

それで、英数字大小文字混合の26+26+10=62種類をまぜて任意の長さの文字列を作ります。
下のソースが作成したスクリプトです。

import string
import numpy
m=100
s = ''.join(numpy.random.choice( list( string.digits + string.letters ), m ))
print s


pythonはたくさんのライブラリがあって便利な反面、簡単なことをするにしても、ライブラリを複数読み込む必要があります。
今回もnumpyとstringを読み込んでいます。
numpyは以前も説明した通り、python標準ではありませんので、別にインストールする必要があります。

ポイントとしては
listで文字列をリストに変換しているところと、
numpy.random.choiceで配列要素から100個サンプリングしているところです。
string.lettersは'abcd.....'というようなすべての英数字を含んだ文字列が欲しかったので使いました。

もちろん、他の方法でもできますが、これが一番理解しやすくて簡潔だと思います。
是非使ってみてください。





2014年7月22日火曜日

mysql

今回はMySQLとの連携です。
PythonからMySQLに接続するためのライブラリはいくつかあります。

ここではMySQL標準の
mysql-connector-python
を使います。

MySQLのConnectorのページに行けばインストーラがあるはずです。
http://dev.mysql.com/downloads/connector/
ここからインストーラをダウンロードしてインストールします。



呼び出すコードは


import mysql.connector

con = mysql.connector.connect( user='x', password='xyz', host='127.0.0.1', database='test' )
cur = con.cursor()
sql = 'CREATE TABLE hoge( id int, value double, comment varchar(32) )'
cur.execute( sql )

con.commit()
cur.close()
con.close()


という感じです。

connectの時の引数は適宜入れ替えて使って下さい。
特にuserとpasswordです。ここでは適当な値を入れていますが、それぞれ、設定した通り入力してください。
sqlの中身に関しても同様です。今回はテーブルを作成していますが、SQLクエリを作成してexecuteに入れます。
ここはsqliteとだいたい一緒です。

データを引き出してくるときは

cur.execute( 'SELECT * from hoge;' )
for r in cur:
   print r

という感じで取ってこれます。
これもsqliteの時とだいたい一緒です。
sqliteとは違うのは、DBのサーバ設定が必要なところです。
そっちはMySQLのサイトなどを参照してください。


補足ですが、
MySQLに保存したデータの実体は
C:\ProgramData\MySQL\MySQL Server 5.6\data
というところにあります。
この下にSchema名のフォルダがあります。そのフォルダ下にtable名のファイルがあります。
commitすると更新されるようです。
逆に言うと、commitを実行しないと、更新されません。



2014年6月23日月曜日

OpenCV 4

Contourとは等高線を指しますが、輪郭という意味もあります。
OpenCVではcontourを輪郭情報という意味で使用しています。

cv2の中にはfindContoursというcontourを計算する関数があります。
つまり、画像から輪郭を抽出することができます。
たとえば、
ret, thresh = cv2.threshold( imgray, 127, 255, 0 )
c = cv2.findContours( thresh, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE )
とすることによりimgrayの画像を閾値処理して輪郭抽出します。

ただし、imgrayというのは2値化された画像で、通常のカラー画像を入れるとエラーになります。
カラー画像imを2値画像gimに変換する場合は、
gim = cv2.cvtColor( im, cv2.cv.CV_BGR2GRAY )
とします。


contourを描画する場合、
cv2.drawContours( im, c, -1, (0,255,0), 3 )
とすると緑色で輪郭情報cを画像imに描画します。

輪郭情報を最小の長方形領域に納める場合
x,y,w,h = cv2.boundingRect(c)
とします。
この長方形領域は
cv2.rectangle(im,(x,y),(x+w,y+h),(0,255,0),2)

とすることによって描画できます。



以上を踏まえて、
WEBカメラから取得した画像から輪郭を抽出して描画するプログラムを作ってみました。
ENTERキーを押すと終了します。


import cv2
import numpy as nmp

cv2.namedWindow('CAMERA')
v = cv2.VideoCapture(0)
v.grab()
im = v.read()

### to binary
gim = cv2.cvtColor( im[1], cv2.cv.CV_BGR2GRAY )
### thresholded
ret, thresh = cv2.threshold( gim, 127, 255, 0 )

c = cv2.findContours( thresh, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE )
for i, cc in enumerate(c[0]):
 cv2.drawContours( im[1], cc, -1, (50*i/255,255,i%255), 3 )
 print '%d,255,%d' % ( 20*i/255, i%255)


cv2.imshow( 'CAMERA', im[1] )
cv2.waitKey(0)


ラボの画像から輪郭抽出してみました。

たくさんの小さい輪郭ができてしまいました。
実際に使う場合は、いかにして目的とするものの輪郭を計算するかという工夫が重要だと思います。

2014年6月9日月曜日

フラクタル 第4回

フラクタルシリーズの4回めです。

今回も、マンデルブロのアレンジを作ってみました。
ただ、毎回関数を別に作るのは無駄なので、P乗で一般化してみました。
つまり、

z(n+1) = z(n)**p + c
z(0)=0

とした時のcと最終のzの値を出力します。


import numpy as np

def makeMandelbrotP( itrn, p=2, n=100, xrange=None, yrange=None ):
if xrange == None:
xrange = (-2.5, 2.5)
if yrange == None:
yrange = (-2.5, 2.5)
r = np.random.uniform( xrange[0], xrange[1], n) + 1j*np.random.uniform( yrange[0], yrange[1],n)

c = r.copy()
z = 0
for i in range(itrn):
z = z**p + c
return c, z

今回はp=4として4乗にした時、このzの値が発散しているかどうかを複素平面上に表示しました。


予想通り、120度の回転対称です。
前回の3乗で180度対称だったので、p乗で306/(p-1) 度の回転対称になるものと想像できます。
この図は16回の繰り返しです。


繰り返し回数については、8回ぐらいでだいたいフラクタルっぽくなります。


繰り返し2回

繰り返し3回

繰り返し4回

繰り返し5回

繰り返し6回

繰り返し7回
繰り返し8回
という感じです。
一枚あたりのデータ量が多く、リスト内表記を使うとメモリアロケーションのエラーになるので、
仕方なくforループを使いました。
通常書類作成用に使っているPC(Core2Quad 2.4GHz)で、上の画像作成時間は一枚あたり6分位かかりました。






2014年6月6日金曜日

pythonからsqlite(DB)の使用

今日はsqliteの話です。
簡単なことしか書いていませんが、しばらく使わないと手順を忘れてしまうこともありますので自分自身の備忘録です。


シミュレーションや分析で、パラメータが多いとデータの管理に手間がかかります。
一度のシミュレーションデータが膨大な場合もどのように保存するかで後々の分析の効率が悪くなったり問題が起こったりします。

一度のシミュレーションにかかる時間が長い(1日とか)と次にどのパラメータでシミュレーション(または分析)をするべきか悩むことがあります。
また、観点を変えて分析をしてみたいと思ったときに、データがそろっているとは限りません。判断をするために足りないデータが何なのかを知る必要があります。


そこで、データの管理のため、データベースを使うと便利なことがあります。


python2.6では標準モジュールとしてsqliteが入っています。
sqliteはフリーのデータベースです。

使い方は


データベースに接続
c = sqlite.connect( 'hogehoge.db' )

切断
c.close()

SQLクエリ送信
sql = 'create table ccc ( iteration integer, realpart float, imagpart float );'
c.execute( sql )

結果受信
sql = 'select * from ccc'
cc = c.cursor()
for r in cc.execute( sql ):
   print r

という感じで使います。

大事なのが、コミット
c.commit()
です。

とりあえずは「これで保存します」の意味だと思って下さって構わないと思います。
これがないと保存されません。

SQLに関してはネット上を検索すればわかりやすいサイトがたくさんあると思いますので、そちらに任せます。




2014年5月29日木曜日

フラクタル 第3回


今度は前回のマンデルブロ集合の変則的な形として、試してみたものです。


今回は

z(n+1) = z(n)^3 + c

という漸化式による収束性について複素平面の図を作ってみました。

できたのは下のような左右対称な図です。


各c(複素数)の点で、zの発散しなかった領域が黒い部分です。緑部分は発散領域です。
この場合もフラクタルのような図形になりました。
図の中央が原点で、左右対称な図形となりました。
範囲は実数軸、虚数軸ともに-1.5から1.5までです。

データ点は約14M 個生成して作成しました。
漸化式の繰り返し回数は10回です。

今回はpythonからデータを生成して、MySQLを使ってデータを管理しましたが、30分程度かかりました。
今回の計算結果はWEBを検索しても簡単に見つからなかったので、ちょっとした課題としていいかもしれません。



作った関数はこんなかんじです。

def makeMandelbrot3( itrn, n=100, xrange=None, yrange=None ):
if xrange == None:
xrange = (-2.5, 2.5)
if yrange == None:
yrange = (-2.5, 2.5)
c = np.random.uniform( xrange[0], xrange[1], n) + 1j*np.random.uniform( yrange[0], yrange[1],n)

z = 0
for i in range(itrn):
z = z**3 + c
return c, z


itrnは繰り返し回数で、データ数はnです。
戻り値はcとzなので、
この関数を使用するときは

c, z = makeMandelbrot3( 10, 1000 )

のようにします。
これで繰り返し10回のデータcとzが1000点分生成されます。
メモリ確保の問題さえなければ、nの値を大きくすればその分データが生成されます。



2014年5月21日水曜日

bool

ブール型の演算について

どの言語もそうですが、整数型や浮動小数点型の演算はだいたい同じですが、ブール型の演算は少しずつやり方が異なっているような気がします。

pythonでの扱い方について記録しておきます。

a = True
b = False
a and b

とすると
False

となります。
a+b
とすると
 1
となり、
a+a
とすると
 2
になります。

つまり、ブール型はTrueが1、Falseが0の値になっていると考えて構わないと思います。

それで、今度はnumpy.arrayでの扱いです。
aa = numpy.array( [ a, a, b] )
とすると
array([ True,  True, False], dtype=bool)

となり、bool型として扱われていることがわかります。
で、今度は
aa and aa
とするとエラーになります。

実行した結果のエラーメッセージは下のようになっていました。
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: The truth value of an array with more than one element is ambiguous.
 Use a.any() or a.all()

numpy.arrayは他の型ではそれぞれの要素に関して個別に演算するので、ブール型もそうかと思うとそうでもないようです。
で、今度は

aa = numpy.array( [True, True, False], dtype=bool)
bb = numpy.array( [False, True, False], dtype=bool)
aa + bb
aa * bb
とすると、それぞれ下のようになりました。

>>> aa+bb
array([ True,  True, False], dtype=bool)
>>> aa*bb
array([False,  True, False], dtype=bool)

つまり、要素個別に演算する場合は+と*でandとorの代用する必要があるということです。
andとor,notは単独のTrue, Falseといったブール値にしか使えません。


notで実現していた、反転をさせる場合はどうするかというと、
~aa
とします。
結果は
array([False, False,  True], dtype=bool)
となります。

ブールのarrayは結構有用です。
たとえば、

x = numpy.random.randint( 0, 10, 12 )
c = x>5
print x[c]

というように、条件にあったものだけ抽出する際に使います。
cは
array([ True, False,  True, False,  True, False,  True,  True,  True,
       False, False, False], dtype=bool)
といった風にbool型のarrayになっています。







2014年5月14日水曜日

フラクタル 第2回

以前、マンデルブロ集合の図を作成しました。
今回pythonはあまり関係ないです。


マンデルブロ集合の画像を線分データのベクター画像(svg)として変換しました。
PDF


そのベクターデータを元に彫刻してみました。

工作関連は別なブログに載せてありますので、詳しくはそちらを見て下さい。


2014年5月12日月曜日

フラクタル

pythonで複素数の練習を兼ねて
フラクタル図形を作ってみました。
マンデルブロ集合というものです。


フラクタルが何かについては専門書または他のサイトなどで調べてみてください。
マンデルブロ集合の定義は簡単な複素数漸化式で記述されます。
つまりマンデルブロ図形というのはある複素数cについて、

z(n+1) = z(n)*z(n) + c
z(0) = 0

とした場合、n→∞で発散するか否かを色分けして複素平面に表示したものです。

複素数平面の点をランダムに選び、n回漸化式を解いてz(n)を作ります。
その時の、
log |z(n)|
の値を画素の輝度値としました。
forループを2重に使うと時間がかかるので、numpy.arrayで一気に計算します。

m = 10000
c = numpy.random.uniform( -2.5,2.5,m ) + 1j * numpy.random.uniform( -2.5,2.5,m )
z = numpy.zeros( m, dtype=numpy.complex )
for i in range(16):
     z = z**2 + c

と言う感じです。

opencvを使って画像化します。
しかし、案外時間が掛かるので、zの漸化式計算と、画像化を分けて処理しました。

Mandelbrot set

30Mポイントのサンプリングをして、そこから画像化しましたが、埋め切れていません。画像を拡大してみるとポツポツと黒い点が見えます。
実は細くつながっています。が、データ点数(サンプリング)で埋め切れていないのと画像化の解像度が荒いため、つながっているようには見えません。






2014年5月9日金曜日

複素数

python numpyでの複素数の扱い方についてです。


複素数は型としてはnumpy.complexと定義されています。
これを利用してみます。

つまり、

numpy.array( [3, 2, 1], dtype=numpy.complex )

とすると

 array([ 3.+0.j,  2.+0.j,  1.+0.j])

というふうに複素数配列となります。


複素数の演算もできます
a = numpy.array( [3, 2, 1], dtype=numpy.complex )
a ** 2

とすると2乗します。
実行結果は

 array([ 9.+0.j,  4.+0.j,  1.+0.j])

となります。

指数計算はそのまま掛けて計算されます。
つまり、python風な数式で表すと
z=x+jy
とすると
z**2 = (x+jy) * (x+jy) = x**2 - y**2 + j 2*x*y
となります。


複素数には色々と実数にはなかった計算方式が出てくるので注意が必要です。
実数部は
numpy.real(z)
虚数部は
numpy.imag(z)
とします。
これらの関数はzがarrayでもOKです。その場合、出力もarrayになります。

複素共役は
numpy.conj(z)
とします。
複素共役が何かわからない場合は検索してみて下さい。
誰かがわかりやすく説明してくれていると思います。

絶対値はnumpy.abs(z)
で計算できます。


複素数の乱数発生は簡単にはできませんでした。
一様分布の場合
z = numpy.random.uniform( -5, 5, 10 ) + 1j * numpy.random.uniform( -5, 5, 10 )
としました。
結果は

array([-4.48658518-4.19169279j,  2.39542899+4.87030019j,
        4.03412469-1.00727465j,  1.72313031+0.99345025j,
        0.68053779-0.385337j  ,  2.08333535-0.32240792j,
       -0.65638447-0.27364433j,  1.36415791-1.89435328j,
        0.19786534-4.51657661j,  2.44853396+3.04850488j])

となりました。
もっとスマートな方法があれば、是非ご教示下さい。





2014年5月7日水曜日

svgwriteのエラー

pythonでSVGファイルを生成する必要があったので、
svgwriteというモジュールをインストールしてみました。

インストールは簡単にできます。
pip install svgwrite
です。



サンプルを実行してみようとすると
エラー…

import下だけなのにエラーが出ました。
エラーメッセージを貼り付けると下のようにxml関連らしいです。

>>> import svgwrite
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "c:\Python26\Lib\site-packages\svgwrite\__init__.py", line 51, in <module
>
    from svgwrite.drawing import Drawing
  File "c:\Python26\Lib\site-packages\svgwrite\drawing.py", line 29, in <module>

    from svgwrite.container import SVG, Defs
  File "c:\Python26\Lib\site-packages\svgwrite\container.py", line 27, in <modul
e>
    from svgwrite.base import BaseElement
  File "c:\Python26\Lib\site-packages\svgwrite\base.py", line 12, in <module>
    from svgwrite.etree import etree
  File "c:\Python26\Lib\site-packages\svgwrite\etree.py", line 33, in <module>
    original_serialize_xml = etree._serialize_xml
AttributeError: 'module' object has no attribute '_serialize_xml'


下のサイトを参考に少し修正しました。
https://bitbucket.org/mozman/svgwrite/pull-request/3/fix-etreepy/diff


pythonインストールフォルダに
Lib/site-packages
というフォルダがあります。
この下にいわゆるライブラリがインストールされます。

インストールされていれば、
この下にsvgwriteというフォルダが有ると思います。
その中にetree.py
というファイルがあるので、

33行目
original_serialize_xml = etree._serialize_xml

を削除して、最終行に移動します。
現状のエラーでは循環呼び出しをしてしまったらしいです。

これで、ひとまずimportは出来ました。

使い方などはまた次の機会にレポートします。



2014年4月25日金曜日

arrayデータの16進パターン

numpyのarrayにはいくつかのメソッドがあります。

覚えておくとたまに便利になる機能がいくつかあります。
まず、tostringです。
動作確認しつつまとめてみます。

tostringメソッドは文字通り、arrayのデータをstringデータとして出力します。
数値を文字列に変換するわけではありません。
どういうときに使用するかというと、場合それぞれですが、
たとえば、シリアル通信やTCP/IP通信などで文字列として受け取ったデータを、数値として扱いたい場合などに使えます。
最近では、画像データの変換として使いました。


例を挙げます。
import numpy
x = numpy.arange(-3,3)
とするとxの中身は
array([-3, -2, -1,  0,  1,  2])
となります。
そこで、x
s = x.tostring()
とすると、sの中身は
'\xfd\xff\xff\xff\xfe\xff\xff\xff\xff\xff\xff\xff\x00\x00\x00\x00\x01\x00\x00\x00\x02\x00\x00\x00'
となり、16進のバイナリパターンが文字列として出力されています。
つまり、arangeで作成した段階で、4[byte](32bit)のint型として扱われ、リトルエンディアン方式で保管されているようです。

これを文字列として保存して、arrayとして戻すためには、
numpy.frombuffer( s, dtype=int )
とします。結果は、
array([-3, -2, -1,  0,  1,  2])
となります。
浮動小数点の場合はdtypeをfloatに指定します。
その他、32bitや64bit、または整数型か浮動小数点型か、符号有無などの細かい設定が必要な場合は、numpyの下に定義されている型を利用します。
たとえば

符号無整数
numpy.uint8
numpy.uint16
numpy.uint32
numpy.uint64

符号付整数
numpy.int8
numpy.int16
numpy.int32
numpy.int64

浮動小数点
numpy.float16
numpy.float32
...

複素数
numpy.complex
numpy.complex64
...

などです。



ちなみに、数値データをそれぞれ文字列に変換する場合は、いくつか方法がありますが、
a.astype( str )
とすればできます。
listとして出力するなら
map( str, a )
としてもよいかもしれません。
今度、それぞれの方法について処理時間などの比較をしてみます。



例で使用したnumpyは1.8.1 (python 2.6)を使用しています。




2014年4月24日木曜日

数理統計学・データサイエンスとpython


pythonを使って統計アルゴリズムを記述することができます。
現在、統計アルゴリズムをコーディングする上でいくつかの比較的簡便な方法が開発されてきました。統計言語Sを基にしたR言語もその一つです。

ここでは最近急速に発展してきたpythonを使った方法について考察してみます。
最近の流れとしては大規模データ、いわゆるビッグデータの分析に注目が集まっているようです。
pythonはこれまでにも書いてきたようにスクリプト言語で基本インタプリタですので、書き方によっては計算時間がネックになってしまうことがあります。
さらには、もともと統計分析を目的として開発されたSまたはR言語のような機能の豊富さには負けます。
しかし、pythonはそれを上回るほどの開発効率やソースコードの簡潔さが特徴で、ライブラリも多岐にわたるため応用範囲が広いという利点があります。

科学技術計算で使うためのライブラリとしては、有名なものだけでも

  • Numpy
  • Scipy
  • Sklearn
  • Pandas
  • Jubatus

などがあります。

Numpyについてはこれまでも何回か説明してきましたが、統計に限らず数値演算の基本的な機能を提供します。
Numpyのもっとも重要な機能はArrayだと思います。
Arrayというのはつまり配列ですが、python標準のlistより、ベクトル的な計算を意識したもので、これを効果的に使うことにより、効率的な計算が可能です。
NumpyにはMatrixという型もあります。
Matrixを使うと行列演算が簡単にかけるので便利ではありますが、
Arrayの方が使い勝手が良いと思います。


ScipyはNumpyではカバーしきれなかった科学技術寄りの機能を実現しています。
たとえば、統計アルゴリズムの例でいうと、乱数発生はNumpyの機能が使えます。つまり、
numpy.random.randn(100)
などとすると100個の正規乱数がnumpy.array形式で生成されます。
しかし、正規分布の確率密度関数のグラフが書きたいとなれば、
scipy.statsの機能が使えます。

import numpy
import scipy.stats
p = scipy.stats.norm.pdf( numpy.arange(-3,3,0.1) )

とすれば、pに-3から3までの正規密度関数の値が入ります。
正規分布のほかにもcauchy分布やgammaやvon misesなどの分布もあります。
統計の論文に出てきそうな標準的確率分布はあらかたあります。


残りのライブラリの説明はまた今度します。
興味があれば、それぞれググってみてください。







2014年3月17日月曜日

C言語との連携 ctypes


pythonは高機能で大変便利ですが、処理が遅いので困る場合があります。
そういう場合は、いろいろと工夫して高速化する方法があります。

ctypesによるDLLの利用はその高速化手法の一つです。
つまり、C言語(など)でDLLを作成して、pythonから読み込みます。
そのDLLを呼びだすためのライブラリがctypesです。

import ctypes

d = ctypes.cdll.LoadLibrary("libx.dll")


C言語と言いましたが、実際C++で作成しても構いません。
但し、呼びだす関数部分はC言語の方式で書かなければなりません。

つまり裏側ではC++のclass定義やインスタンスを生成しても構いませんが、
pythonから呼び出す部分は

extern "C"

などとしてC言語方式で、関数にしなければならないということです。




2014年3月13日木曜日

グラフ

グラフはプレゼンするうえでは必須で、見やすいグラフをどうやって作るかは重要な課題です。
Excelは最近のバージョンで少し見やすくなりましたが、昔はかなりひどいものでした。

これまでMATLABを使ったりRを使ったりしていましたが、
現在は
pythonのグラフ描画ライブラリであるmatplotlibを使っています。
結構グラフがきれいに描画できるのでお勧めです。

import matplotlib.pyplot as ppl
import numpy

t = numpy.arange(100)
ppl.plot( t, numpy.sin( t ) )

ppl.show()




こんな感じで簡単にグラフが描画できます。
行列の値に応じて色分けして表示する場合は、
a = numpy.random.randn(10,15)
ppl.pcolor( a )
とすると表示できます。




2014年3月12日水曜日

並列処理

pythonによるプロセスベースの並列処理に関するメモです。
ライブラリのmultiprocessingを使います。



import multiprocessing as mul

def getLL( d ):
    return -d*d/2

pool = mul.Pool( processes=2 )

data = [ 0.5, 0.2, 0.1 ]
res = pool.map( getLL, data )



getLLは自分で定義した関数です。
dataはリストで、この中身が一つずつgetLLに渡されます。
ここでこの関数は
プロセス毎に実行されます。
ただし、一度に流れるプロセスは最初に決めた値(processes=2)です。

これで簡単に並列計算することができます。
関数への引数は一つしか受け付けないようです。


2014年3月10日月曜日

OpenCV3: WEBカメラ画像の取得

pythonでOpenCVライブラリを利用し、画像処理をします。
今回は単純にWEBカメラ画像を取得して表示するプログラムを書いてみます。


import cv2

cv2.namedWindow('CAMERA')
v = cv2.VideoCapture(0)
v.grab()
im = v.read()
cv2.imshow( 'CAMERA', im[1] )
key = cv2.waitKey(0)


以上で終了です。
C言語などに比べると驚くほど短いです。

これを実行すると、WEBカメラから取得した画像をウィンドウに表示します。
ウィンドウ上でエンターキー入力するとウィンドウが消えてプログラムが終了します。


継続的にカメラ画像を取得して表示したい場合はgrabからwaitKeyまでをループさせてください。
waitKey命令を入れないと画面が更新されません。
waitKeyの引数はキーの入力待ち時間[msec]です。
ループさせると、カメラ画像を継続的に取得表示します。
つまりカメラの映像を、そのまま表示する事ができます。

ソースコードは

import cv2

cv2.namedWindow('CAMERA')
v = cv2.VideoCapture(0)
key = 0
while key != 32:
  v.grab()
  im = v.read()
  cv2.imshow( 'CAMERA', im[1] )
  key = cv2.waitKey(1000)


こんな感じです。
今度はスペースキーで終了します。

スペースキーに32という番号が振られているので、key=32になったらループから抜け出すようになっています。


実行すると
指定したCAMERAと言う名前のウィンドウが出てきます。







2014年3月5日水曜日

コマンドライン1行でpythonのコマンドを実行する方法

これもたまにしか使わないので、書いておかないとよく忘れてしまいます。

linuxなどの端末(terminal)で、スクリプトファイルを作るまでもない些細なことで、
スクリプトなら簡単にできる処理がしたいとします。
且つ、linuxコマンドだけでは複雑すぎて面倒な処理などに、1行スクリプトが効果を発揮します。

例をお見せします。

python -c 'print "\n".join([ "mv %d %d.jpg" % ( i, i ) for i in range(5) ])'

以前紹介したリスト内包表記を使っています。
これを実行すると下のように出力表示されます。

mv 0 0.jpg
mv 1 1.jpg
mv 2 2.jpg
mv 3 3.jpg
mv 4 4.jpg

こうしてコマンド文字列を生成することができます。
文字列がうまくできていることを確認できたら今度は

python -c 'print "\n".join([ "mv %d %d.jpg" % ( i, i ) for i in range(5) ])' | bash

などとして、bashの標準入力に流し、実際に実行できます。
こうしてpythonの機能を使って、効率よく処理できます。
挙げた例でいうと、順番に番号が付けられたファイルをいっぺんに名前を変更するときなどに便利です。


python を-cというオプションで実行するとそのあとの文字列が実行されます。
細かいことはググればヘルプなども出てきます。
実はこれはpearlやrubyなど他のスクリプト言語でも似たようなことが出来ますが、
pythonは文法が少し違うので、慣れが必要です。




2014年3月4日火曜日

OpenCV2 画像の変換


pythonはライブラリが豊富ですが、ライブラリ間のフォーマットは同じであるとは限りません。
画像処理ライブラリに関しても主なものでOpenCVとPILというライブラリがあります。
その間の変換方法について解説します。
この機能、実際使う人がどのくらいいるのかは不明ですが、自分の備忘録も兼ねています。
もし、お役にたつならどうぞ。


CVからPIL
pim = Image.fromarray( cv2.cvtColor( cvim, cv2.cv.CV_BGR2RGB ) )


CVにおけるHSVとRGB,BGRの変換
xx = cv2.cvtColor( x, cv2.cv.RGB2HSV )

実際のところ、画像はnumpy.arrayとして表現されているので、
numpyのスライス機能を使えば同様の処理が可能なことがあります。
たとえば、RGBとBGRの変換では1層目と3層目を入れ替えるので、


tmp = x[:,:,0]
x[:,:,0] = x[:,:,2]
x[:,:,2] = tmp


とします。

反転もnumpy.arrayなら簡単にできます。
y = x[:,range(x.shape[1]-1,0,-1),:]
とすると左右反転します。
スライス機能を使うとなぜかうまくいかなかったのでrange関数使いました。

あと、ロボット用ソフトウェアのためにOpenRTミドルウェアというものを使っています。
pythonでロボットのモジュールが記述できます。
そのなかで使用する画像形式への変換方法についても載せておきます。


RTC.CameraImageからCVへの変更
## ci : cameraImage
### x : numpy.array


ci = RTC.CameraImage( tm, .... )
x = numpy.fromstring( ci.pixels, dtype=numpy.uint8 ).reshape( ci.height, ci.width, -1 )

CVからRTC.CameraImageへの変換
ci.pixels = x.tostring()
ci.width = x.shape[1]
ci.height = x.shape[0]

RTC.CameraImageからPILへの変換cam = RTC.CameraImage(...)
cvim = nmp.fromstring( cam.pixels, dtype=nmp.uint8 ).reshape( cam.height, cam.width, 3 )
pim = Image.fromarray( cv2.cvtColor( cvim, cv2.cv.CV_BGR2RGB ) )

2014年3月3日月曜日

OpenCV 1



pythonにはPILとOpenCVという画像関連ライブラリがあります。

GUIとしてはWxやTkがあり、グラフなどに特化したmatplotlibがあります。

ここでは、OpenCVを中心に説明します。




OpenCV

OpenCVはオープンソースで開発される画像処理ライブラリです。

従来はC++のみであったが、最近C#やPythonにも対応しています。

画像データの形式はnumpyのarray型を使用しています。

numpyは数値計算用のライブラリでいくつも便利な関数が定義されています。

なので、実際C++版よりも使いやすいです。




まだ、windowsインストーラのみでは完全にインストールはされないので、

インストール後手動でcv2.pyd(ダイナミックリンクライブラリ)をコピーする必要があります。

python2.6なら

c:\Python26\lib\site-packages\にコピー

python2.7なら

c:\Python27\lib\site-packages

という具合です。

2014年2月28日金曜日

ctypesによるarrayの受け渡し

計算処理時間がpythonを使っていると問題になることがあります。
pythonはスクリプト言語なので仕方ないといえばそれまでですが、一部分だけでもC言語で書いて実行形式で計算したいと思うことがあります。
ただ、そう思うのが、開発を始めた頃ではなく、開発がほとんど終盤に差し掛かった頃で、移行が簡単ではないことが多いです。

自作class やscipyやnumpyを使い、引数をやりとりしていると、簡単には置き換えしにくくなります。
今回はnumpyのctypesを使ってnumpy.arrayを引数として渡す方法についてまとめてみます。

手順は

  1. c言語でダイナミックリンクライブラリを作ります
  2. pythonプログラムでライブラリを読み込む設定を書き込みます
    1. load_libraryで読み込み
    2. 引数の型設定
    3. 戻り値の型設定
    4. テスト

という感じです。



C言語によるダイナミックリンクライブラリ作成
arrayで渡したい部分は配列にします。
たとえば
double[] hoge( double x[], int n ){
...
}

と言う感じです。
ファイル名はここでは、libhoge.soと言う名前でつくります。
ウィンドウズの場合は拡張子dllにする必要があります。
gccの場合、
 gcc -shared hoge.c -o libhoge.so
とします。


2.1 python側
ダイナミックリンク読み込み
dll = numpy.ctypeslib.load_library( 'libhoge.so', '.' )
通常ctypesというライブラリが標準でついているのそちらを使っていますが、numpyのctypeslibはarrayを使うときに便利です。
numpyの場合、実行ディレクトリを指定する必要があります。


2.2. 引数型の設定
ndtype = numpy.ctypeslib.ndpointer( dtype=numpy.float64, ndim=1, flags='C_CONTIGUOUS' )
dll.hoge.argtypes = [ ndtype, numpy.ctypeslib.ct.c_int,  ]


2.3.戻り値型指定
dll.hoge.restype = numpy.ctypeslib.ndpointer( dtype=numpy.float64, shape=(100,), flags='C_CONTIGUOUS' )


2.4 実行テスト
x = dll.hoge( numpy.array([ 0.1, 0.3, 0.5],dtype=float), 3 )
として実行してみます。
これでxがnumpy.arrayになっていればOKです。


2.2や2.3の型の設定がややこしいところです。



2014年2月24日月曜日

確率統計:モンテカルロ

モンテカルロ計算を考えてみます。

統計を含め、科学技術計算には、numpyが使えます。
行列操作ができ、乱数や確率に関するライブラリが豊富です。
numpyは数値計算用のライブラリであり、
中身の逆行列計算などのアルゴリズムはCなどで書かれています。
そのため、うまく書けばpythonとは言えども高速な計算が可能です。

モンテカルロと言うのはモナコにある町の名前です。
計算機科学の世界でいういわゆるモンテカルロ法というのは確率的なサンプリングを行うことによって、計算を行うことを指し、
フォンノイマンによって命名されたと言われています。

今回はnumpyで簡単なモンテカルロ計算を行います。

numpyをまずimportします。
import numpy

中心0分散1の正規乱数を発生させるには、

x = numpy.random.randn(10000)

とします。

同一正規分布に従う二つの独立な変数x1,x2があったとします。
常にx1<x2であるという制約があった場合、この二つはどのような分布となるでしょうか。

x1 = numpy.random.randn(10000)
x2 = numpy.random.randn(10000)
index = x1 < x2

x1とx2をプロットすると


この場合、切断正規分布という分布になります。
こうするとx1とx2は独立ではなくなり、p(x1,x2)と書くと、x1の確率は

p(x1) = ∫p(x1,x2) dX2

となります。
x1とx2のヒストグラムをとるとこんな感じです。

このようにサンプリングを通して、面積や分布などを調べることをモンテカルロ法といいます。
よく円の面積などの例が教科書などに載っていると思います。
モンテカルロは概念が簡単で容易に理解できますが、奥が深いです。


プロットはmatplotlibを使いました。
matplotlibについてはまたそのうちレポートします。



2014年2月21日金曜日

擬似unixコマンド1

特に仕事でウィンドウズ使っているとunixコマンドを使いたくなることがあります。

そんなの時のために擬似的に似たような動作をする簡易的なpythonプログラムを作っています。
例えばgrep
第1引数の文字列が含まれている行のみ表示するコマンドです。

grep x file

こんな風に使います。
つまり、この場合fileの中身でxが入っている行のみを表示します。
といってもgrepにはもっといろんな機能が付いていますが、簡易版なので、この機能のみを実装してみます。
pythonで書くとかなり簡単に書けます。


import sys

for line in sys.stdin:
    if sys.argv[1] in line:
        print line.rstrip()


処理が遅いと思うのならばリスト内包表記で

import sys
print '\n'.join( [ line.rstrip()   for line in sys.stdin if sys.argv[1] in line ] )

とか書くと少し処理が速くなります。
これをpgrep.pyとかいう名前で保存し、

pgrep.py import pgrep.py
などとコマンドラインで実行すると
文字列importのある行のみ抜き出して表示されます。
いまは作成したプログラムpgrep.pyに対して処理してみましたが、3番目のファイルはテキストならなんでも構いません。
書式は下のとおりです。

pgrep.py [検索文字列] [対象ファイル名]


はじめに

数年前からpythonを使って仕事しています。
科学技術計算のためのソフトウェアとして最近注目されてきているように思います。
ロボット分野だとそれなりにユーザがいて、
ライブラリも充実しています。
ロボット用だけではなくWEB関連もあり、画像処理や統計計算など種類も豊富です。

しかし、新しいことをやろうとすると日本語ドキュメントがないことも多いです。
なので、備忘録もかねて、pythonの使い方についてレポートします。
内容は専門的なことを書きます。

文法構造等の基本的なことは、もっとわかりやすく書ける人がいるので、他のサイトを参考にしてください。



pythonは
使っているライブラリの対応の関係から、version 2.6を使っています。