en_US English ja 日本語 zh_CN 中文

pythonを使ってjpeg画像をまとめて圧縮する方法

久々にwordpressの投稿関連の記事。とは言ってもpython記事だけれども。

さてwordpressユーザーの方はご存知の通り、wordpressのブログにアップできる画像の容量は1Mバイトまで。というわけでデジカメの画像はもちろん、携帯のカメラで撮影した写真でも多少サイズを小さくする必要があることが多い。

数枚ならどんな方法で圧縮してもいいのだけど、これが大量の画像になってくるとそれなりに面倒臭くなっていくもの。
いくつか良いフリーソフトウェアもあるようだけど、こんな時こそpythonを使って見るのも悪く無い。書き出しのファイルサイズを多少はコントロールできたり利点もある。

そういうわけで、この投稿ではpythonを使ってjpeg画像をまとめて1M以下にサイズ圧縮する方法について書いていきたい。

Pillowを使ってjpeg画像圧縮

画像の読み込みはpillowで。画像の読み込みができるpythonパッケージはいくつかあると思うのだけど、簡単な処理をするだけなら大して変わらないんじゃ無いかと信じている。
そんなことないよ。これのが全然早いよってのがあったら教えてください。

必要パッケージのインポート

まずは必要なパッケージをインポート。ファイル名の読み込みようにtkinter各種。それから画像処理にPIL、ファイル名の処理にosをインポート。
ioは画像サイズの処理の関係で利用。詳細は後ほど。

import os
from tkinter import Tk, ttk
import tkinter.filedialog as tkfd
import tkinter.messagebox as tkms
from PIL import Image
from io import BytesIO

tkinterで読み込みファイル名の取得

取り立てて変わったこともないけれど、一応jpgでもjepgでも読み込めるようにしてはある。画像ファイルは./ORIGINALっていうディレクトリに入れておいて一気に読み込めるような設定。

tk = Tk()
tk.withdraw()
print ('select a data files')
filenames = tkfd.askopenfilenames(filetypes= [('jpg','*.jpg'),('jpeg','*.jpeg')], initialdir='./ORIGINAL')
tkms.showinfo('file paths are',filenames)
tk.destroy()

画像サイズを1M以下になるまでjpeg圧縮

というわけでメインパート。画像を縮小するのではなく圧縮する場合だと、jpeg保存のqualityを変えていけばサイズが変わる。この場合qualityを5ずつ落としていっているけれど、そこらへんは適当に調整していただきたい。
for文でqualityを下げながら一回一回ファイルを書き出し、osのstat.st_sizeでファイルサイズのチェックをしている。
ちなみに書き出しは./REDUCEDというディレクトリの中に。

for i in range(len(filenames)):
  img = Image.open(filenames[i])
  basename = os.path.basename(filenames[i])
  for j in range(10):
    img.save('./REDUCED/'+basename+'.jpeg',quality=90-j*5,optimize=True)
    if os.stat('./REDUCED/'+basename+'.jpeg').st_size < 1000000:
    break
  img.close()

しかしこのやり方だと毎回ファイルを書き出して、読み込んでっていうのが少し気持ち悪い。

というわけで以下は実際にはファイルを書き出さない例。ioのByteIOを利用することでそんなことが可能になる。初めて使ったのでこのパッケージの詳細についてはよく把握していない。

qual = [0]*len(filenames)
  for i in range(len(filenames)):
  img = Image.open(filenames[i])
  for j in range(10):
    img_file = BytesIO()
    img.save(img_file,'jpeg',quality=90-j*5,optimize=True)
    #print (img_file.tell())
    if img_file.tell() < 1000000:
      qual[i] = 90-j*5
      break
  basename = os.path.basename(filenames[i])
  img.save('./REDUCED/'+basename+".jpeg",quality=qual[i],optimize=True)
  img.close()

実際にファイルを書きださないで処理しようと思った理由は速度があげられるかなっていうところもあった。

というわけで二つのスクリプトを1.3M、3.5M、4.5Mのjpeg画像を使って速度チェック。
ちなみに書き出されたファイルのサイズは674、989、962[KB]。元のサイズ差のせいなのか、qualityを5ずつ落としているせいなのか目が荒い。

ファイルを書き出す最初のプログラム:

4.87; 4.89; 4.91; 4.87; 4.84

ファイルをメモリに保存しとく2個目のプログラム:

5.34; 5.35; 5.61; 5.22; 5.23

というわけで意外なことに、一々ファイルを書き出すプログラムの方が早かった。メモリを食う方が速度が落ちるのかな。使っているパソコンのシステム次第で変わりそうな気もするけどどうなんだろう。
大した速度の違いがあるわけではないので、個人的には2個目の方が良いかなって思ってるんだけど、どんなもんでしょうか?

変換前後

最後に圧縮前後の画像の一例を載せようと思ったら、当たり前だけど圧縮前の画像がサイズの都合でアップできなかった。なんのためにこの投稿を書いていたのやら。画像を並べたスクリーンショットを参考までに。

関連記事

1. pythonのまとめ

2. pythonでデジカメのraw画像を一括読み込みして別形式で保存する

D

en_US English ja 日本語 zh_CN 中文

pythonでデジカメのraw画像を一括読み込みして別形式で保存する

以前Raw画像を取り込んで処理するっていう投稿を書いた。

内部リンク: pythonでraw形式の画像を読み込んで処理する方法

この投稿はX線データなどの画像処理を書いたものだった。しかし後々考えるとRaw画像の取り込みっていうタイトルだと、写真画像が取り込みたくて訪れてくれてる人もいるのかなとちょっと気になっていた。

それが理由というわけでもないのだけど、この投稿では写真の画像を取り込んで処理するプログラムについて書いていきたい。
こんなプログラムを書いた理由は、容量がきつくなってきた古いパソコンのお掃除がしたかったこと。嫁が大量の写真データを保存していたので、Raw形式で保存しておく必要のないものを一括で読み込んで、jpegかなんかで小さめのファイルで書き出そうという試み。

使用した画像はNikonのデジタルカメラのRaw画像。拡張子はNEF。

rawpyを使ってRaw画像の読み込みと保存

というわけで画像の読み込みに使ったのはLibRawという様々なデジカメ画像を読み込んでくれるライブラリ。pythonの場合だとrawpyという素敵なwrapperがあるのでとても便利。installはpipで。

外部リンク1: LibRawのウェブページ

外部リンク2: rawpyのウェブページ

pip install rawpy

必要なライブラリなどのインポート

まずは必要なパッケージをインポートする。この投稿で使っているのはtkinter, rawpy, imageio, PILのImage,それとos。

from PIL import Image
from tkinter import Tk, ttk
import tkinter.filedialog as tkfd
import tkinter.messagebox as tkms
import rawpy
import imageio
import os

tkinterを使ってファイルの一括読み込み

まずはいつも通りtkinterを使ってファイルを一括読み込みする。

def main():
  tk = Tk()
  tk.withdraw()
  print ('select data files')
  filenames = tkfd.askopenfilenames(filetypes= [('raw','*.NEF')], initialdir='./')
  tkms.showinfo('file paths are',filenames)
  tk.destroy()

これでファイル名の読み込みは完了。読み込みの詳細は別投稿で。

内部リンク: pythonでファイルを読み込む時にファイル名を取得する方法のまとめ

rawpyを使っての読み込みとフォーマットを変えて保存

ファイル名の取得ができたらメインパート。

#forを使って読み込んだファイルを1個1個処理。
for i in range(len(filenames)):
  rawData = rawpy.imread(filenames[i])
  basename = os.path.basename(filenames[i])
#postprocessを使うとnumpyのarrayにデータを詰め込める。最初のオプションはホワイトバランスをカメラの設定で。2個目は明るさを勝手に調整しないっていう設定をオフにしている。
  rgb = rawData.postprocess(use_camera_wb=True,no_auto_bright=False)
#imageioを使ってjpegに書き出すパターン。
  imageio.imsave('test.jpeg',rgb)
#同じくPILを使って。
  img = Image.fromarray(rgb)
  img.save(basename+".jpeg",quality=85,optimize=True)
#保存フォーマットは適当に目的に合わせて。
  img.save(basename+".tiff")

  rawData.close()

というわけでこんな感じでrawpyを使うと、とてもあっさりと写真画像の読み込みができるということがわかった。
上のrgbはnumpyのarrayに詰まっているので、色々と処理をしたかったらここですることができる。
書き出しはまあ何を使ってもすぐ簡単にできるので問題はないでしょう。

rawpyの難点としては少し読み込みが遅いこと。
初期設定を使うと1ファイルの読み込みで4秒ほど、2ファイルをforで回して7秒ほどだった。特にforで遅くなるわけでもなさそう。
demosaic_algorithmを変えたらスピードが結構変わった。どこかウェブの片隅で早いらしいと読んだLINEARを使ったら約2倍の速度になった。何か速さの代償があるのかは知らない。

それでも1ファイル2秒ほどかかるので結構遅い。しかし一回動かしたら放っておけば良いタイプのプログラムなので、別に多少遅くても問題はないでしょう。

左:macのプレビューで書き出し 右:rawpy経由で書き出し

もう一つ問題点としては、ホワイトバランスなど画像読み込みの調整。
これは別にrawpyの問題というわけではなく、Rawデータを開くとどんなソフトだろうと起こる問題だそうな。
Rawデータから元々のデジカメがjpegなどに変換している詳細は企業秘密。カメラのホワイトバランスの値は使うことができるけど、同じ画像にしたかったら詳細をかなり頑張って調整しないといけないのだとか。
以下はそんな感じで色々と設定をいじりながら試していた時の一つ。一番よく映ったオプションというわけではない。
読み込みのオプションの詳細についてはrawpyのホームページに乗っている。

外部リンク: rawpyのオプション

rgb = rawData.postprocess(output_color=rawpy.ColorSpace.sRGB, demosaic_algorithm=rawpy.DemosaicAlgorithm.LINEAR,output_bps=8, use_camera_wb=True,no_auto_bright=True,bright = 1.8,gamma=(1.8,4.0), no_auto_scale=False)

カメラのホワイトバランスを使うuse_camera_wbか、rawpyにホワイトバランスを調整してもらうuse_auto_wbが一番劇的に結果が変わる。
基本的にはこの二つのどちらか好きな方を使ってあとは微調整をする感じだろうか。
それっぽく画像が変換できるオプションを知っている方がいたら教えていただきたいところ。

というわけでこの投稿はひとまずここでおしまい。時間ができたらもう少し他の読み込み方法など試してみたいところ。

関連記事

1. pythonのまとめ

2. pythonを使ってjpeg画像をまとめて圧縮する方法

D

en_US English ja 日本語 zh_CN 中文

pythonでファイルを読み込む時にファイル名を取得する方法のまとめ

久々のpython投稿。

さてだいぶ前にpythonでのテキストファイルの入出力について書いた。

内部リンク: pythonでのテキストデータの入出力方法いろいろ

上記の投稿はどちらかというと開いたファイルの中のテキストデータをどう処理するかということに焦点を絞って書いた。

今回の投稿はちょっとそれとかぶってもしまうところもあるのだけど、最初にファイルを開くところ、ファイル名(もしくはファイルパス)を取得するということに焦点を絞って書いてみたい。

tkinterでファイル名取得

最近私がもっぱらこれを使っているので最初はtkinterから。
自分で使うのにはいちいち窓が立ち上がるし、クリックしなきゃいけないこともあるので鬱陶しさもある。しかし他の人にプログラムを渡すときには、どうやってファイルを読み込むかをあまり説明しなくていいので楽。

また大量のファイルを一度に読み込みたい場合だと、窓でファイルを選択するときにshiftやcontrol/commandを使って大量選択が簡単にできるのも利点。

というわけで以下はtkinterを使ったファイルパス取得の例。

from tkinter import Tk, ttk
from tkinter import filedialog as tkfd

tk = Tk()
tk.withdraw()
print ('select data files')
filenames = tkfd.askopenfilenames(filetypes= [('text','*.xy'),('EXCEL','*.xls')], initialdir='./')
tkms.showinfo('file paths are',filenames)

print ('select a background file')
bkgname = tkfd.askopenfilename(filetypes=[('text','*.xy')], title = 'bkg file', initialdir='./')
tkms.showinfo('bkg file path is', bkgname)

tk.destroy()

大量のデータファイルを読み込むときはaskopenfilenamesを使って、一つのファイルの場合はaskopenfileを使う。
initialdirで窓を開いた時の最初の位置を指定できるので、適当に./DATAみたいなフォルダを作っておいてデータをそこに入れておけば、command+Aで簡単に全部データが開けるかな。

ただし本当に大量のデータを読み込むって場合は、別の方法を使った方が楽かもしれないけど。

tkinterのその他詳しい使い方については別投稿を参照していただきたい。

内部リンク: pythonのtkinterを使ってユーザー入力を取得する方法

argvsでファイル名取得

C言語の時は一番お世話になっていたargvs。

自分で2、3個のデータファイルを読み込もうとする時には一番楽な方法。
Macのターミナルでtabでの名前補完と組み合わせれば、ファイルの読み込みがサクサクとできる。
ただ大量に読み込ませたいファイルがある場合は大変。せいぜい5個くらいまでの読み込みに使う感じだろうか。

import sys
argvs =sys.argv
argc = len(argvs)
if (argc != 3):
  print ("input filename and output name")
  quit()
file = open(argvs[1], 'rb')

pythonではsysをimportしておく必要がある。

if文は読み込みたいインプットの数が正しい数であるかを確認。間違ってたら正しいインプットを要求して落とす。

一番目に読み込んだ名前はargvs[1]になるので注意が必要。argvs[0]はpythonの.pyファイルになる。argvsは他にも単純に適当な文字列を追加するのにも使える。上の例を書いた時には、argvs[2]に出力ファイルネームを入れていた。

そうしてファイル名を取得したらopenやnp.loadtxtなどで開いてあげると良い。

インプット用のファイルを作っておいてファイル名はプログラム内に

プログラムに直接ファイル名を書き込んでしまっておく。必ず何かしらパラメータの入力を受け付けるタイプのプログラムでは便利。input.txtなどのファイルを作っておいて、プログラムを走らせる前にその中のパラメータだけ書き換える感じ。

プログラムの方を書き換えないで良いのが利点と言える。

ファイル名はパスでも指定できるし、同じ階層にいるなら名前を書いておけば良い。

f = open('./input/input.txt')
f = open('input.txt')

など。

ファイルを出力するための名前を取得する方法

ファイルを書き出すのに名前を設定しておきたい場合ってのも結構ある。数少ないなら上のargvsの例のようにユーザーに名前を決めてもらうのも良い。
もちろんプログラム内に直に出力名を決めておいても良い。ただしこの場合プログラムを回すたびに上書きされるのでそれだけ注意が必要。

たくさんファイルをまとめて読み込んだ時に便利なのは、入力ファイル名に何かしら名前を足して出力ファイルの名前に変えてしまう方法。

以下はargus[1]で入力ファイル名を取得して、入力ファイルにsubtという文字を足して出力するという例で。Subtってのは何かしらデータファイル。

name = 'subt'+argvs[1]
np.savetxt(name, Subt)

これを例えばtkinterで大量に読み込んだファイル名に対して同じことをしてあげれば、入力ファイル名に次々と自動的にファイル名が割り振られていくことになる。

とここまで書いて、思い出した。tkinterの場合だとパス取得になってるから、このままだとおかしなことになるんだ。

basename = os.path.basename(filenames[i]+str(i)+'.txt')

こんな感じにos.path.basenameを使ってファイル名に変換してあげる必要がある。

今の所思い浮かんだ名前の取得方法はこれくらい。何か他のものを使う機会があれば後ほど追記していくことにする。

関連記事

pythonのまとめ

D

en_US English ja 日本語 zh_CN 中文

pythonを使ったバックグラウンド強度の評価法いろいろ

色々と科学データを扱うとき、バックグラウンド強度の見積もりをすることは多い。たいていの場合は付属のソフトがついていて、バックグラウンドの計算をして減算してくれる機能がついている。
しかしそんな場合でも勝手に自動でやってくれるならいいのだけど、意外と1個1個マニュアルでやらなければならないこともある。10個くらいのデータならマニュアルでやっても良いのだけど、100個から1000個とかになってくるとプログラムを走らせてまとめて見積もりをさせたくなる。

まとめてファイルを読み込む方法はいろいろあるけれど、下記の投稿などを参照していただきたい。

内部リンク:pythonのtkinterを使ってユーザー入力を取得する方法

ということでこの投稿ではまとめて読み込んだデータを、まとめて処理するのに便利なバックグラウンド評価や減算方法についていくつか書いて行きたい。私は主にX線散乱に使っているけれど、ラマンなど他のスペクトルデータでもなんでも良い。

数式を使ったバックグラウンド計算

以下データはnumpyのndarrayに入った散布データということで。0列にX軸の値、そして1列に対応する強度データなどが入っているとする。下の例ではpwdという名前になっているが、なんでこんな名前をつけたかは忘れた。ndarrayへのデータの収納などに興味があれば下記投稿などで。

内部リンク:pythonでのテキストデータの入出力方法いろいろ

さてどうバックグラウンドを処理するかということだけど、一番都合が良いのはバックグラウンドが理由があるタイプのもので数式に従うはずのケース。理論式をそのまま使っても良いし、フィッティングをしても良い。
さてそれをどうpythonで引くかってことだけど、普通に式を作ってデータポイントごとに引くだけで良い。
まずは意外とよく使うコンスタントな最小値をバックグラウンドとして引く場合。これは最小値をnumpyのnanminで見つけて引くだけで良い。範囲指定(この例の場合XからY)の最小値や平均値を拾ってもいいし、どこかの定点の値を恣意的に引っ張ってしまっても大した問題はないだろう。

#XからYの最小値を探して各データポイントから引く
const = np.nanmin(pwd[X:Y,1])
pwd[:,1] = pwd[:,1]-const

続いて以下はデータをPower Lawでフィッティングをしてから引く場合の例。この例ではpower lawだけど、gaussianでもpolynomialでも好きな関数を使ってフィットすれば良い。フィッティングにはlmfitを使っている[1]。

内部リンク:pythonのlmfitを使ったカーブフィッティング

from lmfit.models import PowerLawModel

amp, exp = fit_PL(pwd)

def fit_PL(pwd):
#データファイルの44から70と574から584のエリアを使ってフィットしたかったのでnumpyのconcatenateで結合。astypeで変換してたのはデータファイル読み込みの都合だったと思う。たいていの場合なくて良いはず。
y = np.concatenate((pwd[44:70,1],pwd[574:584,1]),axis=0)
y = y.astype(float)
x = np.concatenate((pwd[44:70,0],pwd[574:584,0]),axis=0)
x = x.astype(float)

#lmfitのデフォルトのPowerLawModelを使ってフィット
PL_mod = PowerLawModel(prefix ='PL_')
pars = PL_mod.make_params()
pars['PL_amplitude'].set(10000)
pars['PL_exponent'].set(-1.2, min=-4, max = -1)
res = PL_mod.fit(y, pars, x=x)
print (res.fit_report())

#Power Lawのパラメータを取得
amp = res.best_values.get('PL_amplitude')
exp = res.best_values.get('PL_exponent')

#プロットするならmatplotlibなどで
plt.plot(pwd[:,0],pwd[:,1])
plt.plot(pwd[:,0],amp*(pwd[:,0]**exp))

#Power Lawパラメータをリターン
return amp, exp

この例の場合だと、データの最初と最後を使ってpower lawにフィットしたかったので、numpyのconcatenateを使ってそのようなデータセットを作ってからフィットした。
そしてamplitudeとexponentの二つのパラメータを関数から返す。

pwd[:,1] = pwd[:,1] - amp*(pwd[:,0]**ex)

この二つのパラメータを使って各データポイントの値を計算し、引き算をすればPower Lawの引き算がおしまい。こんな感じで式を使う場合は、簡単に引くことができる。

バックグラウンドがシンプルな場合のベースライン補正的なバックグラウンド演算

上記のような計算できるバックグラウンドならば良いのだけど、計算式によらないいろいろ混ざったバックグラウンドが引きたいってことも結構多い。以下はそんな時に使えるバックグラウンドの見積もり方について。こういった手法って結構いろいろ転がっていて、pythonだと割と簡単に導入できることが多い。

まずはバックグラウンドがどんどん下がって行くや、どんどんと上がっていくのようなシンプルなケースに使えるバックグラウンドの引き方。アイディアは少し古い1997年のクレイのX線研究の論文から[2]。たぶんオープンアクセス。
私がやったのは論文の手法と全く同じではないのかな。少しわかりにくかった場所があったのでそこらは適当に書いたので。結果はほぼ一緒になっていると思うけど。

大雑把にいうと、ある地点から直線を引く。直線の傾きを少しずつ変えていってデータのカーブと接したところでバックグラウンドラインとする。でその接点からまた直線でのバックグラウンド評価をやり直すということを繰り返す。基本的にはピーク以外のところがフラットになる操作。

装置的要因や、アモルファス的なバックグラウンドなどを考えて引くものではない。どんな実験したんだっていうような安定感のないデータ群を、とりあえずフラットなベースラインに揃えたい場合に使う感じ。
アモルファスなどの非ピーク情報を得たい場合にはあまりよろしくない手法。根元をがっつり削っちゃうので。
あとはスペクトルデータのベースライン補正などにも使えるかな。

#まずはスタート地点(init_ang)を指定する。conv_lineは値を行番号に変換する定義関数。意外と便利というかよく使う。
init_ang = 9.5
init_ang = conv_line(pwd,init_ang)
lipbkg = fit_lip(pwd,init_ang)

#バックグラウンド演算の関数fit_lip()。lipbkgにバックグラウンドのデータを詰めていく。
def fit_lip(pwd,init_ang):
  pwd = pwd.astype(float)
  plt.plot(pwd[:,0],pwd[:,1])
  step = 1
  init_loop = init_ang-1
  lipbkg = pwd[:,1]

#バックグラウンドを決める直線を引き直すためのループ。回数はデータ次第だけどある程度大きければ良いはず、適当。
  for j in range(100):
    init_loop = init_loop+step+1

#直線の開始地点がデータポイントの行数を超えたらループ終了。この場合582だったけどもちろんデータ次第。
    if init_loop > 582:
      print ('break')
      break
#バックグラウンドの直線とデータのカーブが接しているかチェックするループ。直線は単純に一次関数で定義して、傾きを10ずつ変えている。データ次第で刻みは適宜調整する必要がある。
    for i in range(40):
      slope = -(i*10)
      sec = pwd[init_loop,1] - slope*pwd[init_loop,0]
      linear = slope*pwd[init_loop:,0]+sec
      #plt.plot(pwd[init_loop:,0],linear,alpha=0.2)
      diff = pwd[init_loop:,1]-linear
      step = 1
#生データと直線の差分の最小値がプラスになったところで直線とカーブが接していると判断。argpartitionで何個目のデータポイントで接しているのかIndexを抜き出して、バックグラウンドの強度を得るのと次のループでのスタート地点を決めるのに利用。
      if np.nanmin(diff[1:]) > 0:
        step = np.argpartition(diff[1:],0)[0]
        print('i,j',i,j)
        lipbkg[init_loop:init_loop+step+2] = linear[0:step+2]
        break

  plt.plot(pwd[:,0],lipbkg[:])
  plt.show()

  return lipbkg

#データの値と行番号を変換する関数。差分の一番小さくなるIndexをargpartitionで抜き出すだけ。よく使うので比較的速度の速いargpartitionで。
def conv_line(plot,num):
  num = float(num)
  dif = np.absolute(plot-num)
  mini = np.argpartition(dif,0,axis=0)
  return mini[0][0]

というわけで、ある点から引いた直線がカーブと接するたびにその地点からバックグラウンド評価をやり直すことで、直線的なベースラインを得る手法でした。
評価範囲を限定しながら、直線の角度をもっとフレクシブルにすればデコボコしたバックグラウンドにも使えるだろう。しかしこのプログラムを書いた時は一様に減少するバックグラウンドだったのでこれだけで大丈夫だった。forループたくさんだし早くもないけれど、困るほど遅くもなかったと思う。

ピークvs滑らかなアモルファスカーブのスムージングを利用した分離法

こちらはシャープなピークとそうじゃないものをざっくばらんに分けたいという場合に有用。元々はX線のバックグラウンドで開発された手法だけれども[3]、スペクトロスコピーのベースラインの平滑化とかにも使える。
必ずしも正確にバックグラウンドを表現できる手法とはいえないかもしないが、いくつか利点がある。
比較的早い。プログラムの途中に組み込みやすい。評価手法次第だがある程度の定量性がある。演算の回数やスムージングの条件次第でブロードピークの除去から、スペクトルのベースラインを平滑化したい時まで利用できる。比較的複雑なバックグラウンドにも対応できる。などなど。
アイディアとしてはピークを含んだデータに対してスムージングを行うと、ピークは潰れるけどピークじゃないところは大体そのままになる。そしてスムージング後のデータとオリジナルデータを各データポイントごとに比べて値の小さな方を残す。これをどんどんと繰り返すことで、ピークが速やかに潰れていき最終的にバックグラウンドが残ることになる。弱点は大体どんな評価法でもそうだけどピークのオーバーラップ。

というわけで以下はスクリプト。オリジナルの論文とは違う論文のものを参考に、savgol_filterを利用したもの[4]。

#下に書いたbkg_estimate関数を使って最初にバックグラウンドを潰す。
bkg_init = bkg_estimate(Imin, Iave, pwd)
#bkg_iteration関数を一度回して1回目のバックグラウンドを得る。
BackG = bkg_iteration(bkg_init,pwd,lin1,lin2,win)
#後は好きなだけスムージングサイクルを回せば良い。
last = 599
for j in range(1,last+1):
  BackG = bkg_iteration(BackG,pwd,lin1,lin2,win)

#バックグラウンドの初期値を決める関数。のちのバックグラウンド減算の回数を減らせる効果があるのだけど、今のコンピュータだと誤差範囲。オリジナルの論文が書かれた2000年には効果があったのだろう。
def bkg_estimate(arr1, arr2, arr3):
  Imin = arr1
  Iave = arr2
  Int = copy.deepcopy(arr3[:,1])
  plt.plot(arr3[:,0],Int)
  for i in range(len(arr3[:,1])):
    if arr3[i,1] < (Iave +2*(Iave-Imin)):
      Int[i] = arr3[i,1]
    else:
      Int[i] = Iave +2*(Iave-Imin)

  plt.plot(arr3[:,0],Int)
  return Int

#ピークをスムージングでつぶしながら、バックグラウンドを計算する関数。
#引数はarr2は諸事情でプロット用に入れていたので、普通はなくても良い。それからlin1,lin2もデータ範囲の指定に使ってるだけなのでなくても良い。winは結構この引き方に致命的に効くパラメータなので、メインの関数から書き直せるようにしているので引数になっている。もちろん直接関数の中に書けば省略できる。
def bkg_iteration(arr1,arr2,lin1,lin2,win):
  bkg = arr1
  area = 0
#scipyのsavgol_filterを使ってピークをスムージングで潰す。
  Nbkg = savgol_filter(bkg,win+1,1,mode='nearest')
#あるデータ範囲でスムージング前後の小さい方の値をバックグランドとして残すことでピークがどんどん潰れて行く。
  for i in range(lin1,lin2):
    if (bkg[i]-Nbkg[i])<0:
      bkg[i] = bkg[i]
    else:
      bkg[i] = Nbkg[i]
  plt.plot(arr2[:,0], bkg)
  return bkg

状況によっては千回とか回すので、forループの中でいちいち保存するよりも自分自身(BackG)を書き換えた方が良いのかなとプログラムを書いたのだった。
しかしこういう書き方は絵面があまりよろしくない。どういう風に書くと美しいのだろうか?なんかPython的な書き方もありそうなものだけど。

いずれにしろこの手法の一番の問題点はいつスムージングを終わらせるかというところ。
特にバックグラウンドを使うわけでなければ目で見て好きなところで止めれば良いけれど、バックグラウンドを使いたい場合はある程度定量性のある止め方を考えないといけない。しかしこの投稿の本題ではないので深くは突っ込まない。

参照文献
[1] lmfit, https://dx.doi.org/10.5281/zenodo.11813
[2] Clays and Clay Minerals, 45, 132-146, 1997.
[3] Journal of applied chemistry. 33, 977-979, 2000.
[4] Carbohydrate polymers. 78, 543-548, 2009

関連記事

pythonのまとめ

D

en_US English ja 日本語 zh_CN 中文

pythonを使って座標変換計算 〜回転行列を使って座標回転

座標計算はプログラムを書いておくと便利。 一度書いてしまえば使い回しもしやすいし。
無料のソフトウエアもあるけど融通が効かない場合もあるので。

ということでこの投稿ではpythonでの座標計算、主に回転行列を使って任意軸周りに座標を回転変換する方法について書いていきたい。

pythonでgenfromtxtを使った座標の読み込み

読み込みは色々とあるでしょうけれどご自由に。
xyz座標をfloatなりdoubleなりで取得しておく。 私は以下のような文字と数値の混ざった結晶pdbファイルからの読み込み。

ATOM   1  C4  N   1    -0.094  -0.316   0.079 1.000  2.63
ATOM   2  C3  N   1    -0.818   0.983  -0.299 1.000  2.63
ATOM   3  O3  N   1    -1.800   1.360   0.658 1.000  2.63

残りの文字列などを残したまま編集したかったのでgenfromtxtを使用してnumpyのarrayに投入。 ヘッダーが四行ついていたのでskip_header=4で読み飛ばす。
書き出す時ににヘッダーを元に戻したいならば、処理した座標を書き出す前に別にヘッダーを読み込んでファイルに書き込んでおけば良い。

orig = np.genfromtxt(argvs[1],dtype=['U8','int','U8','U8','int','float','float','float','float','float'],skip_header=4)

pythonでの座標の移動

とりあえず座標の移動は足し算でも引き算でもしとけば良い。
のだけれどgenfromtxtで違うデータ型を読み込んだ場合だと、少しややこしいことがある。
今上の例でorigをプリントすると(3,)というshapeを持つ、一見タプルを突っ込んだ三行のリストのように出力される。 でも実はこれはタプルではなくstructured arrayというもの。

[('ATOM', 1, 'C4','N', 1, -0.094, -0.316,  0.079, 1., 2.63)
 ('ATOM', 2, 'C3','N', 1, -0.818,  0.983, -0.299, 1., 2.63)
 ('ATOM', 3, 'O3','N', 1, -1.8  ,  1.36 ,  0.658, 1., 2.63)]

さてこれ一見タプルだけどタプルじゃない。 そこでそのまま値の変更が可能なのが良い点。
しかし同じデータ型で読み込んだ時とは異なり、普通のarrayではない。 そのためスライシングなどを使って、一気に値を変更する方法が見つけられなかった。

私の例の場合はデータの5列目から7列目の値を変更したかったのだ。 とりあえず下の例ではAdd[i][5]からAdd[i][7]に対してfor文を回して書き換えをしている。

Adcord = addcoords(orig, 6.8,1.8,42.0)
def addcoords(orig,num1,num2,num3):
 Add = copy.deepcopy(orig)
 for i in range(len(orig)):
 Add[i][5]=num1+Add[i][5]
 Add[i][6]=num2+Add[i][6]
 Add[i][7]=num3+Add[i][7]
 return Add

しかしこの投稿を書くのにいろいろ調べていたのだけど、こういう風にindexをつけてデータを抜き出してる例ってほとんどない。
というかgenfromtxtでデータ型の混ざったテキストを読み込んでそのまま処理する例自体があまり見当たらず。

きっともう少しスマートなやり方はあるのだろう。
だけどこのくらいならfor文1万回回したとしてもそんなに遅くはないのでとりあえず良しとしておく。

それからここでdeepcopyを使ってるのは、関数内で引数データに対して処理をして返すのに必要だったのだと思う。
これもベターな書き方があるんだろうけど、調べるのが面倒臭くてスキップしてしまっているものの一つ。

pythonで回転行列を使っての座標の回転

さてようやく本題の座標の回転へ。
私がやりたかったのはある2点を通る線を軸にしてある座標を回転させるということ。

まずは任意軸周りの回転行列の定義。 ロドリゲスの公式というのだったかな。 引数のAngが角度でnが任意軸の単位ベクトル。

def rotmat(Ang,n):
 R = np.array([[np.cos(Ang)+n[0]*n[0]*(1-np.cos(Ang)), n[0]*n[1]*(1-np.cos(Ang))-n[2]*np.sin(Ang), n[0]*n[2]*(1-np.cos(Ang))+n[1]*np.sin(Ang)],
[n[1]*n[0]*(1-np.cos(Ang))+n[2]*np.sin(Ang), np.cos(Ang)+n[1]*n[1]*(1-np.cos(Ang)), n[1]*n[2]*(1-np.cos(Ang))-n[0]*np.sin(Ang)],
[n[2]*n[0]*(1-np.cos(Ang))-n[1]*np.sin(Ang), n[2]*n[1]*(1-np.cos(Ang))+n[0]*np.sin(Ang),np.cos(Ang)+n[2]*n[2]*(1-np.cos(Ang))]])
 return R

上記の行列を使って回転させる関数の定義。 引数は座標と角度。 角度はどっかでラジアンに直しておく必要あり。
ちょっと大きめの座標だったのだけど、七十七行目と0行目の座標を軸(n1)に残りの座標を全部回転させている。 linalg.normでベクトルの大きさを出力し、それで割って単位行列にしてある。

for文の範囲が半分なのは座標の半分だけ回転させたかったから。 ここら辺は任意で。

一度全座標を軸ベクトルの始点で引いてから回転させている。 がちょっと前に作ったスクリプトなので、どうしてこうしていたのかは忘れてしまった。 こうしないとうまく動かなかったんだっけ?

回転自体はnp.dotを使って各座標に対して行列積を計算してあげるだけ。

def rotresi1(orig,thmo):
 new = copy.deepcopy(orig)
 Ang = (thmo*np.pi)/180
 n1 = np.array([new[77][5],new[77][6],new[77][7]])-np.array([new[0][5],new[0][6],new[0][7]])
 n1 = n1/np.linalg.norm(n1)
 ROT = rotmat(Ang,n1)
 for i in range(0,int((len(orig)/2))):
  Sc = np.array([orig[0][5],orig[0][6],orig[0][7]])
  Coord = np.array([orig[i][5],orig[i][6],orig[i][7]])
  Coord = Coord -Sc 
  new[i][5]=np.dot(ROT,Coord)[0]+orig[0][5]
  new[i][6]=np.dot(ROT,Coord)[1]+orig[0][6]
  new[i][7]=np.dot(ROT,Coord)[2]+orig[0][7]
return new

んでたくさん違う角度で回した座標とかを一気に吐き出したかったら、下のようにfor文を回してあげれば良い。 for文は避けても良いけれど、このくらいの計算量ならまいっかなってくらいの速度。
データ数が莫大ならばちょっと工夫した方が良いかもしれない。

最後は回転した座標をheaderを一行だけつけてフォーマットを揃えてsavetxtで出力。 ここら辺も任意で。

for i in range(36):
 rot = 10*(i+1)
 New = rotresi1(orig,rot)
 np.savetxt('adj'+str(rot)+'.pdb',New,fmt='%4s %6d %3s  %-4s  %2d %11.3f %7.3f %7.3f %4.3f %5.2f',header = 'CRYST1    5.130   19.280   10.660  90.00  90.00  97.00',comments='')

というわけでpythonを使った座標計算の一例でした。
最近こういう計算から離れてしまっているので、もう少し取り組んでいきたいところ。

関連記事

pythonのまとめ

D

 

 

en_US English ja 日本語 zh_CN 中文

pythonを使って画像の外れ値探しと修正処理

取得した画像に望まれないそして理由のはっきりしている例外が含まれていることはよくある。 そんな時は見た目的にも解析的にも修正できるなら修正しておきたいところ。
理由がわからない場合はそれが事象かもしれないので、修正してしまうのは問題だけど。

この投稿では二次元画像にはっきりと区別できる望まれない何かしらが写り込んでしまった場合の修正の一例について書いてみたい。

画像の取得

raw画像の取得について、以前の投稿で書いたので詳細に興味がある方は参照していただきたい。

リンク: pythonでraw形式の画像を読み込んで処理する方法

tiffやjpegなどだったらもっとシンプルに読み込めるのじゃないかな。
いずれにしろPILなどを用いて画像を読み込み、それを処理のためにnumpyのndarrayに突っ込んでおくのが基本方針だろう。

画像の中から例外値の取得

X線画像

ここではX線の画像に、X線を止める物質が混入したという例で。 まあビームストップなんだけど。
画像は360x1500ピクセルの画像で、上記したようにnumpyのndarrayにピクセルごとに入れてある(上の写真はトリミングしてるので縦800ピクセルになっているけれど)。
ndarrayのshapeがそのまま(1500, 360)になるのでどう処理していけば良いかは直感的にわかりやすい。

この例の場合はビームストップの周りはピクセル強度があからさまに落ちるので(画像の縦に走る黒い線)、強度の低いデータ場所を抜き出すというのが最初の方針になる。
単純に弱いのを抜き出すのでよければ、適当に決めた閾値よりも小さい値を一気に抜き出せば良いので簡単。

ただこの例の場合だと画像の下に行くほど比例的に強度が落ちて行くため、その方法だと画像下部が大量にマスクされてしまう。 それはあまりよろしくない。 いや実はこの画像の例ではよろしかったんだけど、よろしくないことにしておく。

さてじゃあどうするか。 とりあえず閾値を横方向の1データセットごとにセットすることにした(画像を横方向にスライスして最小値を拾う)。
savgol_filterでスムージングフィルターをかけているのは、横方向スライスした時に極端な外れ値があった場合に取り除けるように。 だけど私の画像は割と平坦なのでほとんど働いていなかった。
ちなみにビームストップの強度を拾わないように画像左半分[:,0:180]の範囲から最小値を拾っている。

#最初に閾値を1行あたりに決めるためのarrayを作成
threshold = np.linspace(0,1499,1500)

#0列から180列まででデータ処理をする。 polar_dataが画像を格納したndarray。
polar_data_smo= polar_data[:,0:180]

#1500行について1行ずつスムージングフィルターをかける
for i in range(1500):
  polar_data_smo[i] = savgol_filter(polar_data[i,0:180],11,4)

#閾値に1行ごとの最小値を格納していく
  threshold = np.nanmin(polar_data_smo,axis=1)

というわけでこれで閾値が取得できた。
実用的にはthresholdから多少値を引いてあげて余裕をもたせても良いかもしれない。 拾いたい値がどの程度通常の値と離れているか次第。

例外値をどうするか

均一な画像の場合

私の例はそうではないのだけど、均一なデータの場合はNANに値を変えてしまうのが楽。 均一データは平均値で物が扱えるので、平均を取るときにnanmeanなどを使えば良い。

arrayの中である値より小さな値を探す処理はnumpyのwhereが便利だと思う。

polar_data = msk_small(img2,origin,polar_data,threshold)

#閾値以下の値をNoneに変える関数
def msk_small(img2,origin,polar_data,threshold):

#1行ずつ閾値より小さな値を探してNoneを代入していく
 for i in range(1500):

#whereを使って閾値より小さな値を探す
  polar_data[i,np.where(polar_data[i,0:360]<threshold[i])] = None

#データをPILのImageで画像に出力してうまく処理できているか確認
  img = Image.fromarray(polar_data)
  img.save("test.tif")
  return polar_data

不均一だけど対称的な画像の場合

私の例の画像はこの例。
対称的だから例外値のところを使わなければ良いのだけど、ちょっと例外値の場所の周りも使いたい事情がある場合。
とりあえず最小値をthresholdに取得してあるので、それを流用して最小値で埋めてみる。 見た目的には多少マシになる。

 polar_data[i,np.where(polar_data[i,0:360]<threshold[i])] = threshold[i]

というわけでできあがった図は下のような感じ。 左側が修正後の画像。
ぱっと見はそんなに悪くはないような気もするけれど。 並べてみればまだまだはっきりとわかる。
まあ絵としては別にいいけれどデータ解析には使えない。

 

ついでにこれだとビームストップの下に特徴(ピーク)があるような場所は全然ダメ。
やっぱり適当に対称の場所をフィットしてそれを強度を合わせるように持ってきた方が良いのか。 しかしあんまりごちゃごちゃ処理するのもこういう画像処理としてはやりすぎかな。

ビームストップの位置を変えて2回データを取るなど、外れ値を作らない努力をする方が幸せかもしれない。

まあこういう処理でうまくいくような画像もきっとあるのだろう。 きっと。

不均一で何の規則もないような画像

というのもちょっと考えてみたけれど試していない。 そもそもこういう場合だと、例外値の取得がまずもって難しそうだし。 ちょうど良い画像の例もみつからなかったので諦めた。

というわけで何やら失敗談の投稿になってしまった。 せっかく書いたので投稿はしておくが後でもう少し取り組むチャンスがあれば修正・加筆するかもしれない。

何かこういう画像処理が簡単にできるパッケージ・ライブラリなどをご存知の人がいたら教えていただけるとありがたい。

関連記事

pythonのまとめ

D

en_US English ja 日本語 zh_CN 中文

pythonを使って網羅的計算 〜高速計算を求めて

というなかなか自己矛盾なタイトルでの投稿。 pythonはスクリプトの書く時間を短縮するのが利点であって、pythonに計算をやらせるのはそんなに早くないことが多いので。
だからと言ってふと思い立ったときにCを書くかと言うとちょっと手が渋るのも事実。 この投稿ではそんな大量計算をpythonで頑張ってやってみた時の経験について書きたい。

私的結論から書いておくと道は二つ。 一つはやはり面倒臭がらずにC系統で書く。 pythonを使うことにこだわるならpython的な高速化は求めずに、愚直にスクリプトを書いてからcythonに書き直したりnumba@jitで高速化を目指す。 このどちらかの方が結果として早くなるんじゃないかなと思った。

しかしこの投稿で書くのはそう言う方法ではなく、なるべくpython的に頑張ろうとした時の試行錯誤をいくつか書いていきたい。 もちろんpythonianな人たちはもっとエクセレントな方法を知っているだろう。 気づいたことがあったら教えていただけるととても助かる。

pythonでどのような網羅的計算をしたかったか

ちょっとここの説明が長いので、この章を飛ばしてもらってもなんとなくわかると思う。
4つのパラメータa,b,c,dがある式。 この式は3つの整数変数h,k,lの総組み合わせに対してそれぞれ解を吐き出すような式。 例えば

s = h*a+k*b+c*l+d

といった感じで。 実際はもう少し複雑な計算をさせているが。 いずれにしろ10個ずつ整数変数があるとしたら一つのa,b,c,dに対して1000個のs値が計算されることになる。

一方でこの時の実験値が8個。 この実験値と上記で計算された値sとが良く合うパラメータa,b,c,dのセットを探したいというのがやりたかったこと(どのsと実験値を合わせても良い)。
ある程度情報があって最小二乗法などで解けるなら何も問題がない。 そうでない場合はちょっと考えないといけない。

いくつかある方法の中で一番愚直な方法は、とりうるパラメータa,b,c,dを細かい刻みで動かして全ての値のを網羅的に計算する。 その中で一番実験値と合う値を引き算で算出。 その値を叩き出したパラメータa,b,c,dを求めるという方法だ。

プログラムの構成

a,b,c,dのパラメータ変数を適当な間隔で定義する。 各a,b,c,dから式の値を計算して、実験値との差分を計算する関数を回す。 差分が例えば0.001以下のパラメータセットを出力する。 そんな単純プログラム。 しかしa,b,c,dを増やせば増やすほどガンガンと計算量が増えるプログラムだった。

for文を避けてnumpyで頑張る

上記のようなプログラムを書きたい時はCだったらfor文をガスガス回すと言うのは正しい戦略だと思う。 わかりやすいし早いし。
ただpythonはfor文が割と遅いという特徴がある。 for文を如何に避けるかはpythonで重たい計算をするときの一つの方針にはなる。
と言うわけでまず4つのパラメータに対して式を計算するのにfor文を回すのではなく、numpyのarrayに収納してみることをトライ。 どの程度速度が変わるかをまずは単純計算で比較。 4つのパラメータを10個ずつ取ったとき。 まずは単なる網羅的足し算での比較を。

for文:

import time
import numpy as np

def main():
 start = time.time()
 for a in range(10):
 for b in range(10):
 for c in range(10):
 for d in range(10):
 a+b+c+d
 elapsed_time=time.time()-start
 print ('elapsed_time',elapsed_time)
 return
結果
elapsed_time 0.00111 [s]

array文:

import time
import numpy as np

def main():
 start = time.time()
 a = np.linspace(0,9,10)
 b = np.linspace(0,9,10)
 c = np.linspace(0,9,10)
 d = np.linspace(0,9,10)

 ar1 = np.array(np.meshgrid(a,b,c,d)).T.reshape(-1,4)
 np.sum(ar1,1)

 elapsed_time=time.time()-start
 print ('elapsed_time',elapsed_time)
 return
結果
elapsed_time 0.00107 [s]

あれ? ちょっとだけ早いけど10000回の計算だと大して変わらない。 しかもこのnumpyのarrayの組み方は弱点がある。 網目を大きくするとarrayがとっても大きくなる。 その分だけメモリを食うという弱点。
というわけで網目を30x30x30x30にした場合。

for: elapsed_time 0.0681 [s]
array: elapsed_time 0.1256 [s]

網目を50x50x50x50にした場合(6250000回計算)。

for: elapsed_time 0.51124 [s]
array: elapsed_time 1.07741 [s]

というわけでこれだけだと実はforの方が良かった。

さてこれをちょっとだけ複雑な計算にしてあげるとどうなるか。 各30パラメータの例で、変更は二乗の和をとるだけ。 変更と結果は下記の通り。
for:
a**2+b**2+c**2+d**2
elapsed_time 0.90076 [s]

array:
ar1 = ar1*ar1
ar2 = np.sum(ar1,1)
elapsed_time 0.13798 [s]

forは二乗を取るだけでかなり遅くなる。 一方でarrayは大して変わらない。 というわけでここでだいぶarrayの方が早くなった。

とはいえこういう単純計算だったらもっと早く書く方法があるのかな。 pythonianな方に教えて頂きたいところ。

ループで使用する関数内での変数宣言を避ける

これは回数の少ないループだとむしろ中に書いておいた方がわかりやすいかもしれない。 しかし100万回計算するループの中に変数宣言が入っていたりするとpythonはとても遅くなる。
普段はあまり使いたくないgrobal変数だけど、この場合は使用を検討してみると早くなるかもしれない。 grobal変数は、関数の外でgrobalに続けて宣言すればよい。

grobal h
h = np.linspace(0,5,6,dtype=int)

など。

argminよりargpartitionを使う

本題に戻って私がしたかったことは、実験値と計算値を比べること。 これは単純に差分を取ってあげて、そこから最小値を抜き取ればよい。 ここでパッと考えるとargminを使いたくなるところ。 しかしargpartitionという関数があり、こちらの方が早い。

dif2 = np.absolute(slist - dobs)
min =np.argpartition(dif2,0,axis=0)

こんな感じで差分を取ったarrayの最小のインデックスが抜き出せる。

@numba.jit

時折とても強力に高速化してくれる。 けれどむしろ遅くなる場合もある。 全体にかけて効くときもあるし、関数1個だけにかけた方がよく効くときもある。 よくわからないけれど、効いてくれたら嬉しい機能。 試す価値はあり。

数学関数はmathかnumpyか

物によってはmathのが早いし、numpyの方が早かったこともある。 というかあまり変わらない。 最後の最後にどっちのが早いか見てみてもよいかも。 100万回も計算すれば違いが出るかもしれない。

anacondaのpythonを使う

私がanacondaのpython3に乗り換えたのはこれが原因。 なぜかちょっとだけ計算が速くなる。 もっと早いのがあるなら教えて欲しい。

計算量を減らす努力

小さな努力が大きな結果に跳ね返る場合もある。 計算に対称性などがあって計算が減らせる場合は減らしておく。 計算量が多い場合だとパラメータが減ったり式が単純化されるだけで結構な違いに跳ね返る。 なんだかんだでマニュアル努力が大切。

終わりに

というわけでいくつか私がpythonの高速化を計った時に使った方法を書いてみた。 私は結局目的に使い物になるくらいのところ(1回の試行が帰宅前に回し始めて朝には終わっている)で諦めてしまったのだけど。
最終盤のプログラムでもCで適当に書いた似たようなプログラムに数十倍の速度差をつけられていた。 本職pythonの人がどのくらい早くできるのかは気になるところではある。

cythonへの変換も試した。 しかしいろいろnumpyで書いていたせいか、変数の宣言などcythonを単純に使っても速度を上げることはできなかった。
というわけである程度計算量が必要なプログラムを書く時は、最初からcythonを意識した書き方をしておく。 もしくはcythonフォーマットで書き始めるというのが正しい道なのかな。

関連記事

pythonのまとめ

 

D

en_US English ja 日本語 zh_CN 中文

エクセルにスクリプト言語としてpythonが導入される?

エクセルにpythonが導入される?

最近エクセルにpythonが導入されるかもしれないという記事をいくつか読んだ。 実際に導入されるかはまだ決定はされていない模様。

python as an Excel scripting language

ことの発端は上記のマイクロソフトのユーザーフォーラム。 pythonの導入がだいぶ前にリクエストされていたそうなのだ。 もともとのポストは2015年の11月4日だから2年前。 それに対してマイクロソフトチームが昨年末の2017年11月14日に突如コメントをつける。 どういう需要があるのかを調査したいのでアンケートに協力してねっ、てことなのでポジティブなコメントと言えるだろう。 しかしアンケートしてるぐらいだしまだまだ初期も初期段階。

マイクロソフトエクセルチームからのアンケート

せっかくなのでアンケートに答えてきた。

  1. どこ住んでるの?
  2. お仕事の種類はお金目的?
  3. 会社の種類はどんな感じ?
  4. 研究開発的な専門職についてるの?
  5. どんだけの間python使ってる?
  6. pythonは何に使ってるの?
  7. エクセルとpythonが統合するとしたらどんなことに一番便利かな?
  8. 7と同じ質問だけど、こっちは一番じゃないので何個も答えて良い
  9. エクセルにpython導入されたら君の仕事のどんな役にたつの? 特に7、8の質問と関連づけてくれるとわかりやすいな。
  10. pythonスクリプトは自分で使ってるの? 他の人に提供してるの?
  11. いままでに何かエクセルとpythonを繋げるようなソフト使ってた?
  12. 他にどんなプログラミングする?
  13. 何か他にも意見ある?
  14. もうちょっと協力してくれるならemailアドレス教えて

重要な質問は7−9と11くらいかな。

エクセルとpython統合に関する私見

さてエクセルもpythonもデータ解析にちょろっとずつ使っている私の意見。
自分用にだけならば今のところこの二つを統合する利点は感じていない。 大量計算を簡単にしたいだけならpythonを使う。 もっと大量に早く計算したいとなるとむしろCとかにいく必要がある。
元データがエクセルデータなら簡単な計算はエクセルでそのまま計算。 そうじゃないならばpandasでpythonにインポートして計算させる。

ただし他の人のためのスクリプトを書くとなるとちょっと話しは別かな。 サイエンス界隈の人たちでもエクセルを使ってる人は多い。 裾野を広げたらもっと多くなるだろう。 何かしらスクリプトを書くとしたらデータがエクセルフォーマットでやってくることは多い。 しかも中身が統一されていない様々なフォーマットだったり。 そんな時にいちいちスクリプトを書き直すのは面倒臭い。 エクセルにネイティブにpythonが組み込まれていれば、そこらへんの融通は効きやすくなるんじゃないかなと期待はしている。

といわけでエクセルとpythonの統合で何かが便利になりそうという期待はあまりしていない。 しかし統合されたならば使ってみたい、と思うくらいの興味はあるというのが現状。 上のアンケートなどでとても便利でエクセレントなアイディアを出す人たちがいると良いのだけど。
またエクセルのような大御所がpythonを導入するならば、pythonコミュニティーがまた広がるんじゃないかなと思う。 pythonってわかりやすいし書きやすい良い言語だと思うのだよね。 遅いけど。

関連記事

1. pythonでのテキストデータの入出力方法いろいろ

2. pythonのまとめ

D

en_US English ja 日本語 zh_CN 中文

pythonを使って積分計算とシグマ計算

pythonを使って積分計算をする機会があった。 備忘にここにまとめておく。

方位角プロット

データは上の図の青ラインのようなガビガビの散布データ。 numpyの360行のarrayに入っている。 横軸は360度の角度で縦軸は強度。 やりたかったことは1度刻みで強度を持つデータに、その角度ごとに三角関数をかけてからその総和を取るということ。
過去の文献をのぞいて見るとシグマで計算している人とインテグラルで計算している人がいる。 もともとの式の定義はインテグラル。 なのでインテグラルを選択するのが順当だろうけど、せっかくなのでpythonを使ってできる計算を色々と書いてみようと思い立った。
まずできるのはデータポイントごとに、三角関数の計算をしてそれからnumpy.sumをとってのシグマ計算。 それからやっぱり積分計算の方が良いかと思い、同じプロセスからnumpy.trapzを使っての積分計算。 とはいってもポイントごとの値のばらつきが大きいデータなのでやっぱりピークフィットしてから、その関数を使って積分計算をかけた方が良いかなとも思い、scipy.integral.quadを使っての積分計算。
というわけでこれらの方法での結果と実行速度を比べてみることにした。

生データでのシグマ計算

まずは横軸を0度から359度まで作っておく。 ついでに三角関数用にradianも準備しておく。 例のごとくimportは省略。 numpyとscipy、lmfit、それからTimeくらいかな今回は。

pix = np.linspace(0,359,360)
AzTheta = np.radians(pix)

コンスタントバックグラウンドの準備。 ガビガビのデータなので外れ値を引いてもしょうがない。 numpy.nanminを使って適当な範囲で適当な最小値を選択。

const = np.ones(90)*np.nanmin(Int[140:220])

バックグラウンドを引いたピークの裾からピークトップまでのデータ。 それから各データポイントにそれぞれ0から90度までのsinとcosをかけた式。

eqsum = (Inttheta[88:179]-const)
eqsin = (Inttheta[88:179]-const)*(np.sin(AzTheta[0:91]))
eqsin = (Inttheta[88:179]-const)*(np.cos(AzTheta[0:91]))

そうしたらsumを使ってarrayの総和をとる。

sum1 = np.sum(eqsum)
sum2 = np.sum(eqsin)
sum3 = np.sum(eqcos)

printで出力。

sum1 66885.0328331
sum2 61156.5549622
sum3 15573.1626916
Calculation time: 0.0002231597900390625[s]

生データからのインテグラル計算

上と同じものをsumではなくnumpy.trapzで計算。 trapzはtrapezoidal ruleでの積分計算。 関数じゃなくても生データから積分計算ができる。

sum1 = np.trapz(eqsum)
sum2 = np.trapz(eqsin)
sum3 = np.trapz(eqcos)

1度刻みのデータなので少しだけ総和よりも小さくなるのは想定通り。 しかしtrapzの方が速度が早かったのは少し意外。

sum1 63719.308712
sum2 58081.6139909
sum3 15482.3795419
Calculation time: 0.00015211105346679688[s]

フィッティングからのインテグラル計算

こちらはscipyのintegrate.quadを用いて。  フィッティングは以前書いたので省略。 私はいろいろ便利なのでlmfitを使っている。 投稿下にカーブフィッティングへの関連記事リンクを貼っておいた。 ここではconstantバックグラウンドとareaGaussianでフィットした(上の図の赤点線)。
それから積分用の関数を書く。 シンプルなX依存の関数ならそのままintegrate.quadの中に式を書けば良い。
他のパラメータが入る式を使いたい時は先に式を定義しておく。 それからintegrate.quadのargsでパラメータの値を指定することができる。 今回のように最初に関数のフィットをして、その結果の値を積分に突っ込みたいときにはこの方法が便利だと思う。

def eq1(x, amp, cen, hf):
return amp*(2/hf)*np.sqrt((np.log(2)/np.pi))*np.exp((-4*np.log(2))*((cen-x)/hf)**2)

def eq2(x, amp, cen, hf):
return (amp*(2/hf)*np.sqrt((np.log(2)/np.pi))*np.exp((-4*np.log(2))*((cen-x)/hf)**2))*(np.sin(np.radians(x)))

def eq3(x, amp, cen, hf):
return (amp*(2/hf)*np.sqrt((np.log(2)/np.pi))*np.exp((-4*np.log(2))*((cen-x)/hf)**2))*(np.cos(np.radians(x)))

式のパラメータにカーブフィッティングのベストの値をresult.best_values.getを使って代入。 それから0から90まで積分計算を行う。

s1 = sp.integrate.quad(eq1, 0, 90, args = (result.best_values.get('g1_amp'),90,result.best_values.get('g1_hf')))
s2 = sp.integrate.quad(eq2, 0, 90, args = (result.best_values.get('g1_amp'),90,result.best_values.get('g1_hf')))
s3 = sp.integrate.quad(eq3, 0, 90, args = (result.best_values.get('g1_amp'),90,result.best_values.get('g1_hf')))

生データはピークの裾が持ち上がっているので、フィットしてからの積分値はだいぶ値が異なっている。 それをどう解釈するかはデータ次第。
フィットをもうちょっと真面目にするとしたらVoigtなどを使ってフィットしても良いし。 このデータの場合だと生データのバックグラウンドの取り方を変えても値は近くなるかな。

sum1 48043.715584253885
sum2 47717.13305405471
sum3 4456.932291343271

時間はフィッティングも合わせるとだいぶ遅い。 積分だけでも生データより10倍遅くなっている。 とはいえよっぽど大量のデータ処理じゃなければ問題はない速度。

total time: 0.011394023895263672[s]
integration time: 0.002457857131958008[s]

というわけで結構値が変わるので、どう計算するかを結構慎重に考える必要がある。 大体の場合先行研究をいくつか見れば答えが書いてあるのだけれど。

関連記事

1. pythonのlmfitを使ったカーブフィッティング

2. pythonのまとめ

D

en_US English ja 日本語 zh_CN 中文

pythonのtkinterを使ってユーザー入力を取得する方法

tkinter、wxpython、pyqtあたりが有名どころのようである。 私はシンプルなことがしたいだけで、見た目も割とどうでも良いということでtkinterを使うことにした。 見た目にこだわったり複雑な処理をしたい場合には、後者二つの方が良いらしい。
使う機会があればこれらについても後ほど書くかもしれないがそんな複雑なスクリプトを書く機会には恵まれないんじゃないかなとも思う。 今の所はtkinterでハッピー。

tkinterをimportするときのエラー

最初にtkinterを使うにあたって一つ最初に躓いたことある。 たぶんインストール設定次第だと思うのだけど、matplotlibと一緒にインポートするとクラッシュする。 matplotlibでTkAggを使うように指定することによって回避できた。

import matplotib
matplotlib.use('TkAgg')
from matplitlib import pyplot as plt
from tkinter import *

といった感じで。

tkinterを使ってユーザー入力補助のスクリプト

ここでは適当なスペクトルデータを表示するプログラムを例に。

tkinterでファイルの読み込み

さてtkinterを使ってまずしたかったのは入力データファイルの読み込み。 自分で使うのにはlsなりdirなりで確認してargvsで読み込めば良い。 もしくはインプットファイル名を決めておいて、勝手に読み込ませるとか。
しかし今回プログラムをphdの学生に渡すのにちょっと手間取った。 特に彼女はあまりpythonというかスクリプトを使わないので。 またwindowsユーザーなのでコマンドライン操作にも慣れがない。
そうなるとプログラムを走らせた時にウインドウを開き、ファイルをクリックして開けるようにした方が間違いがないかなと思ったのだ。
前置きが長くなったが、まずはtkinterの必要なものをインポート。 他のインポートは省略。 brewer2mpl以外はありふれたものしかインポートしてないはず。

from tkinter import Tk, ttk
from tkinter import filedialog as tkfd

それからダミーのルートファイルを開いて隠しておく。 これをしておかないと意味のない小窓が開くので煩わしい。

tk = Tk()
tk.withdraw()

続いてaskopenfilenamesを使ってファイルのパスを取得しfilenamesに格納。 ファイルタイプは適当に。 初期ディレクトリは私はファイルと同じ場所にプログラムを持ってくことが多いので./でその場を指定。 もちろんデスクトップなど好きなところで良い。

filenames = tkfd.askopenfilenames(filetypes= [('txt','*.txt')], initialdir='./')

この操作でこの場合だとテキストファイルを複数開くことができる。 このfilenamesはタプルになるので、必要がある場合はリストに変換する。 取得するのはファイルパス。 なのでファイル名などが欲しい場合はosを使ってコンバージョンする。 他の方法もあるかも。
さてファイルパスが取得できたら一度メインウインドウを閉じておく。 これをしないで他の関数で再びメインウインドウを開くと止まってしまうことがある。

tk.destroy()

filedialogは色々あるのでちょうど良いものを使うと良い。

tkinterでユーザー入力の数値の取得

さてファイルを取得したら、そのファイルをどうしたいかなどをユーザーに決めてもらいたい場合がある。 ここではプロットする範囲を取得するという例で。

num1,num2 = get_range('plot')

get_rangeは以下の自己定義関数。 インデントは省略。

def get_range(word):
#インプット窓の動作を設定するクラス。
class Getr:
#クラス変数の初期化。
input1 =0
input2 =0
#コンストラクタでルートウインドウの修飾をする。
def __init__(self,root):
#インプット窓のラベルの設置。 
#場所の指定はgridが便利だと思うのだけど他はどうなんだろう。
label0 = ttk.Label(root,text = 'input')
label0.grid(row=0,column=1)
#インプット窓のラベルと入力欄の設置。
word1 = 'low ' + word + ' limit'
label = ttk.Label(root,text = word1)
label.grid(row=1,column=0)
self.ent = ttk.Entry(root)
self.ent.grid(row=1,column=1)
self.ent.focus_set()
#同じようにもう一つラベルと入力欄の設置。
word2 = 'high ' + word + ' limit'
label2 = ttk.Label(root,text = word2)
label2.grid(row=2,column=0)
self.ent2 = ttk.Entry(root)
self.ent2.grid(row=2,column=1)
self.ent2.focus_set()
#ボタンの設置と、ボタンを押した時の動作の設定。
#lambdaの理由は忘れたけど、これないと動かなかったはず。
but = ttk.Button(root, text='OK',width=10,command = lambda: self.get_quit())
but.grid(row=3,column=0)

#ボタン動作用のメソッドをいくつか。
def gettext(self):
Getr.input1 = self.ent.get()
Getr.input2 = self.ent2.get()
def quit(self):
root.destroy()
def get_quit(self):
self.gettext()
self.quit()
#クラスはここでおしまい。
#インプットルートウインドウ。 ここら辺はテンプレ。
root = Tk()
#メインウインドウでクラスGetrを。
range=Getr(root)
root.mainloop()
#インプットの取得。
return Getr.input1, Getr.input2

さてこうしてユーザー入力で得られる数値とプログラムで使いたいデータの列番号は異なっていたりするので、列番号に直しておく。 ここでargpartitionを使っている理由は早いからだけど、この場合はあまり関係ない。

def conv_line(plot,num):
num = float(num)
dif = np.absolute(plot-num)
mini = np.argpartition(dif,0,axis=0)
return mini[0][0]

こうして得られたライン1からライン2までのプロットを打つ。 ちなみに想定してるのはデータはx軸とy軸の2列データ。

def plot_raw(filenames,num1,num2):
num = len(filenames)
print (num," files read")
if num >2 and num <11:
bmap = brewer2mpl.get_map('spectral','diverging', num)
bmap2 = brewer2mpl.get_map('OrRd','sequential', num)
colors =bmap.mpl_colors
colors2 =bmap2.mpl_colors
for i in range(num):
plot = np.loadtxt(filenames[i])
lin1 = conv_line(plot,num1)
lin2 = conv_line(plot,num2)
plt.plot(plot[lin1:lin2,0],plot[lin1:lin2,1], '-',color=colors[i],mec=[1,0,0],ms=8,lw=0.5, alpha=1)
return

といった感じのスクリプトでした。 あえてpythonでプロットを打つ必要はないのだけど、打てると便利な場合もある。 私はバックグラウンド補正とか強度のキャリブレーションを読み込んだ複数ファイルにまとめてかけて、それからプロットを出力するのに使っている。 というわけでtkinterの一例でした。

関連記事

pythonのまとめ

D