カテゴリー: 科学

久々にpython投稿。

ちょっと論文を書くのにy軸二つ並べて、モデルから計算した値と実験値を並べたかったので、pythonで描写。

シンプルに公式ホームページの通りに書いていけば問題はないと思う。

外部リンク：　Plots with different scales

importはmatplotlibの設定くらい。

import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt

まずは適当に好みの図の描写設定。

plt.rcParams['font.family']='Arial'
plt.rcParams['font.size']= 14
plt.rcParams['axes.linewidth']=2.5
plt.rcParams['xtick.major.width']=2.5
plt.rcParams['xtick.labelsize']=14
plt.rcParams['ytick.major.width']=2.5
plt.rcParams['ytick.labelsize']=14
plt.rcParams['figure.figsize']=(3,4)

それから適当にプロットするモデルのデータ。

X1 = np.array([-0.5])
SLD1 = np.array([15.8])
SLD2 = np.array([6.25])
SLD3 = np.array([11.53])
SLD4 = np.array([2.9])

いつもの。

fig, ax1 = plt.subplots()

ガスガス１軸目をプロット。それから目盛りのサイズやらラベルやらの調整。

ax1.plot(X1,SLD1,marker='v',ms=11,color='orange',fillstyle='none')
ax1.plot(X1,SLD2,marker='D',ms=9,color='g',fillstyle='none')
ax1.plot(X1,SLD3,marker='x',ms=9,color='m')
ax1.plot(X1,SLD4,marker='+',ms=9,color='b')
ax1.tick_params(axis ='x', which ='both', top='off',bottom='off', pad=10)
ax1.set_xticklabels([])
ax1.set_ylabel(r"$\Delta$"+'SLD$^{2}$'+'x'+'10$^{20}$'+ ' [cm$^{-2}$]',fontname='Arial')
ax1.set_ylim(0,20)

２軸目追加。

ax2 = ax1.twinx()

２軸目用の実験データ。

X2 = np.array([0.5])
I1 = np.array([388])
I2 = np.array([187])
I3 = np.array([313])
I4 = np.array([35])

２軸目データのプロットと軸のメモリとラベルの設定。

ax2.plot(X2,I1,marker='v',ms=11,color='orange',fillstyle='none')
ax2.plot(X2,I2,marker='D',ms=9,color='g',fillstyle='none')
ax2.plot(X2,I3,marker='x',ms=9,color='m')
ax2.plot(X2,I4,marker='+',ms=9,color='b')
ax2.tick_params(axis ='x', which ='both', top='off',bottom='off', pad=10)
ax2.set_xticklabels([])
plt.xlim(-1,1)
ax2.set_ylabel('INT [cm$^{-1}$]')
ax2.set_ylim(0,400)

最後に絵を出力。figを使っても使わなくても出力は変わんない？

fig.tight_layout()
fig.savefig('FIGURE/test1.png', bbox_inches='tight')
plt.savefig('FIGURE/test2.png', bbox_inches='tight')
plt.show()

絵はこんな感じ。

というわけで、そんなにスケールがあってなかったよ、という絵。

関連記事

1. pythonのまとめ

D

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

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

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

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

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()

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

というわけでメインパート。画像を縮小するのではなく圧縮する場合だと、jpeg保存のqualityを変えていけばサイズが変わる。この場合qualityを５ずつ落としていっているけれど、そこらへんは適当に調整していただきたい。
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を５ずつ落としているせいなのか目が荒い。

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

4.87; 4.89; 4.91; 4.87; 4.84

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

5.34; 5.35; 5.61; 5.22; 5.23

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

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

関連記事

1. pythonのまとめ

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

D

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

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

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

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

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

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

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

外部リンク１：　LibRawのウェブページ

外部リンク２：　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を使って読み込んだファイルを１個１個処理。
for i in range(len(filenames)):
  rawData = rawpy.imread(filenames[i])
  basename = os.path.basename(filenames[i])
#postprocessを使うとnumpyのarrayにデータを詰め込める。最初のオプションはホワイトバランスをカメラの設定で。２個目は明るさを勝手に調整しないっていう設定をオフにしている。
  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の難点としては少し読み込みが遅いこと。
初期設定を使うと１ファイルの読み込みで４秒ほど、２ファイルをforで回して７秒ほどだった。特にforで遅くなるわけでもなさそう。
demosaic_algorithmを変えたらスピードが結構変わった。どこかウェブの片隅で早いらしいと読んだLINEARを使ったら約２倍の速度になった。何か速さの代償があるのかは知らない。

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

左：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

久々の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。

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

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

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

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

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

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

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

以下データはnumpyのndarrayに入った散布データということで。０列にX軸の値、そして１列に対応する強度データなどが入っているとする。下の例では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):
#データファイルの４４から７０と５７４から５８４のエリアを使ってフィットしたかったので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だと割と簡単に導入できることが多い。

まずはバックグラウンドがどんどん下がって行くや、どんどんと上がっていくのようなシンプルなケースに使えるバックグラウンドの引き方。アイディアは少し古い１９９７年のクレイの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

#直線の開始地点がデータポイントの行数を超えたらループ終了。この場合５８２だったけどもちろんデータ次第。
    if init_loop > 582:
      print ('break')
      break
#バックグラウンドの直線とデータのカーブが接しているかチェックするループ。直線は単純に一次関数で定義して、傾きを１０ずつ変えている。データ次第で刻みは適宜調整する必要がある。
    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関数を一度回して１回目のバックグラウンドを得る。
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)

#バックグラウンドの初期値を決める関数。のちのバックグラウンド減算の回数を減らせる効果があるのだけど、今のコンピュータだと誤差範囲。オリジナルの論文が書かれた２０００年には効果があったのだろう。
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的な書き方もありそうなものだけど。

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

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

関連記事

pythonのまとめ

D

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

ということでこの投稿では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ではない。　そのためスライシングなどを使って、一気に値を変更する方法が見つけられなかった。

私の例の場合はデータの５列目から７列目の値を変更したかったのだ。　とりあえず下の例では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文１万回回したとしてもそんなに遅くはないのでとりあえず良しとしておく。

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

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

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

まずは任意軸周りの回転行列の定義。　ロドリゲスの公式というのだったかな。　引数の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

上記の行列を使って回転させる関数の定義。　引数は座標と角度。　角度はどっかでラジアンに直しておく必要あり。
ちょっと大きめの座標だったのだけど、七十七行目と０行目の座標を軸（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

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

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

画像の取得

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

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

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

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

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

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

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

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

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

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

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

＃閾値に１行ごとの最小値を格納していく
  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):

＃１行ずつ閾値より小さな値を探して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]

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

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

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

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

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

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

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

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

関連記事

pythonのまとめ

D