Python

Python③Jupyterで木星の画像解析

Pythonの画像解析ライブラリpillowとフィッティングライブラリlmfitを使ってJupyterでJupiterの画像を解析します. 第一近似では所詮点光源といってもガウス関数でのフィッティングは上手くいかず, 必要に応じて補正項を取り入れたローレンツ関数に適切なパラメータ初期値と範囲を与えることでなんとかそれっぽい結果を得ることが出来ます.

使用環境

WSL2 Ubuntu-20.04
Anaconda3 4.10.3
Python3.8.11

pillowによる画像の下処理

これが木星とその衛星の画像です. 私の師匠が望遠鏡とスマートフォンで撮影したものです. 4つある衛星のうち3つが肉眼でも確認できます.

この画像をJupiter_and_sattelites.jpgという名前で適当なディレクトリに保存し, 同じディレクトリでノートブックを作成し以下のセルを実行していきます.

from PIL import Image
import numpy as np
#画像の読み込み
im = Image.open("./Jupiter_and_sattelites.jpg")
#RGBに変換
rgb_im = im.convert('RGB')
#画像サイズを取得
size = rgb_im.size

#im.show()
#print(im.format, im.size, im.mode)
array = np.asarray(im)
array
R    = array[...,0]    # Ellipse notation!
G    = array[...,1]
B    = array[...,2]
S    = (R + G + B) 
im_R = Image.fromarray(R)
im_G = Image.fromarray(G)
im_B = Image.fromarray(B)
im_S = Image.fromarray(S)
#im_R.show()
#im_G.show()
#im_B.show()    # Halo can be seen!
item = G
Image.fromarray(item).show()

とりあえず人間の目が感知しやすい緑色成分に関して解析していきます. 各ピクセルのGreen成分の強度を格納した3次元配列Gをitem変数に代入し画像を確認します.

Green成分のみ抽出しグレースケールで改めて出力した画像.

少し寄り道です. 以下のコードでこの画像のピクセルごとの値をz軸に持つ3次元プロットを出力出来ます.

import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d.axes3d as axes3d
fig  = plt.figure(dpi=200)
ax   = fig.add_subplot(111, projection='3d')
ax.set_facecolor('whitesmoke')
plt.rcParams['font.family'] ='sans-serif'# 使用するフォント
plt.rcParams['font.size'] = 10 # フォントの大きさ

x  = np.array(range(item.shape[1]))   # not [0]
y  = np.array(range(item.shape[0]))   # not [1]
xx,yy = np.meshgrid(x,y)
print(item.shape)
ax.plot_surface(xx,yy,item, cmap='plasma') # cmap: spring,summer,autumn,winter,gray,hot,cool,etc...
ax.set_xlabel('y pixel #', size = 14)
ax.set_ylabel('x pixel #', size = 14)
ax.set_zlabel("pixel intensity", size = 14)
ax.set_zlim(0,255)
ax.set_title('3D plot of the pixel intensity')
ax.view_init(elev=30, azim=210)
Green成分のピクセル値をz軸に表示した3次元プロット.

以下のコードでは, x方向に関してはy方向のピクセル平均値を, y方向に関してはx方向のピクセル平均値を計算し散布図として表示しています. pillowで取り込んだ今回の画像ファイルはどういうわけかx方向とy方向がそれぞれ縦方向と横方向に対応しているようなので, x,yが逆の様に扱われています...

今回はこのように各方向に平均化したデータに対してフィッティングを行うことで, 1つのピークを持つ関数で画像の全体的な特徴を掴もうという趣旨です.

fig  = plt.figure(dpi=100)
ax   = fig.add_subplot(111)
ax.set_facecolor('whitesmoke')
plt.rcParams['font.family'] ='sans-serif'# 使用するフォント
plt.rcParams['font.size'] = 10 # フォントの大きさ
# x軸に補助目盛線を設定
ax.grid(which="major",axis="x",color="grey",linestyle = "--",linewidth=1)
# y軸に目盛線を設定
ax.grid(which="major",axis="y",color="grey",linestyle="--",linewidth=1)
y_average = np.average(item, axis=0)
x_average = np.average(item, axis=1)
plt.scatter(x,y_average, label='y_averaged data',s=1)
plt.scatter(y,x_average, label='x_averaged data',s=1)

plt.title('Averaged intensity')
plt.xlabel('Pixel #'); plt.ylabel('Intensity')
plt.xticks(np.arange(0,1000,50), rotation=45)
#plt.xlim(0,1000); plt.ylim(0,100)
plt.legend(loc='upper right', bbox_to_anchor=(1, 1))

lmfitによるフィッティング

以下のコードではlmfitに与えるモデル関数の候補としてガウス関数とローレンツ関数の2つを用意します.
ただし通常のガウス関数, ローレンツ関数ではなく, saturation係数の掛っている補正項を加えてモデルの表現力を高めています. この項は1978年のNewton, Andrews等による水素原子2S-2P Lamb shift分光におけるモデリングを参考にしています.

def Gaussian(x, height, cen, width, saturation, offset):# width equals sigma in the case of Gaussian
    L1   = (x-cen)/width
    Q1   = np.exp(-(((x-cen)/width)))**2
    val1 = height * (Q1 + saturation*(Q1*(Q1-1/2)*(Q1-1)))
    return val1 + offset
def Lorentzian(x, height, cen, width, saturation, offset):
    L1   = (x-cen)/width
    Q1   = 1/(1+L1**2)
    val1 = height * (Q1 + saturation*(Q1*(Q1-1/2)*(Q1-1)))
    return val1 + offset 

fitting_function = Lorentzian # choose fitting model

以下のコードでlmfitによるフィッティングを実装します.

from lmfit import Model
fit_model = Model(fitting_function, nan_policy='raise')  # raise error on NaNs
# configuration of fitting parameters initial values or so
par_init   = [30,400,50,0,0]                
par_limits = [[0,0,0,0,0],[255,1000,100,10,30]]
par_names  = ['height','cen','width','saturation','offset']
par_val    = dict(zip(par_names,par_init))
par_min    = dict(zip(par_names,par_limits[0]))
par_max    = dict(zip(par_names,par_limits[1]))
par_vary = {'height':True,'cen':True,'width':True,'saturation':True,'offset':True}

# 定義したモデル用のParameters クラスオブジェクトの宣言
params = fit_model.make_params()
for name in par_names: 
    params[name].set(
        value=par_val[name], # 初期値
        min=par_min[name], # 下限値
        max=par_max[name], # 上限値
        vary=par_vary[name] # パラメータを動かすかどうか
)
# data to fit 
x_data = y_average
y_data = x_average

result_x = fit_model.fit(x=x, data=x_data, params=params, method='leastsq')
print(result_x.fit_report()) #fitting result summary
#fig, gridspec = result.plot(data_kws={'markersize': 3},
#                           numpoints=1000,show_init=False,
#                            xlabel="x pixel #",ylabel="x pixel intensity") #fitとresidualsをax.subplotで出力,詳細の設定
result_y = fit_model.fit(x=y, data=y_data, params=params, method='leastsq')
print(result_y.fit_report()) #fitting result summary
#fig, gridspec = result.plot(data_kws={'markersize': 3},
#                           numpoints=1000,show_init=False,
#                            xlabel="y pixel #",ylabel="y pixel intensity") #fitとresidualsをax.subplotで出力,詳細の設定

以下のコードでフィッテイング結果と元のデータを同じグラフに出力します.

# Graphing
fig  = plt.figure(dpi=100)
ax   = fig.add_subplot(111)
ax.set_facecolor('whitesmoke')
plt.rcParams['font.family'] ='sans-serif'# 使用するフォント
plt.rcParams['font.size'] = 10 # フォントの大きさ
# x軸に補助目盛線を設定
ax.grid(which="major",axis="x",color="grey",linestyle = "--",linewidth=1)
# y軸に目盛線を設定
ax.grid(which="major",axis="y",color="grey",linestyle="--",linewidth=1)

# output fitted array by eval() method 
plt.plot(x,result_x.eval(), label='fitted x averaged values')
plt.plot(y,result_y.eval(), label='fitted y averaged values')
plt.scatter(x,y_average, label='x averaged data',s=1)
plt.scatter(y,x_average, label='y averaged data',s=1)

plt.title('Fitted averaged intensity')
plt.xlabel('Pixel #'); plt.ylabel('Intensity')
plt.xticks(np.arange(0,1000,50), rotation=45)
#plt.xlim(0,1000); plt.ylim(25,60)
plt.legend(loc='upper right', bbox_to_anchor=(1, 1))

良い感じです.
ここではGreen成分に関するフィッテイングのみ紹介しましたが, item変数をR,B,Sなどに変更することでRed, Blue及び各成分のSumに対するフィッティングも同様に行えます. これらのおおよそ正しそうなフィッティング結果を与えるparametersのinitial valuesとlimitsを参考までに貼っておきます. 色々いじって遊んでみてください!

# best parameters initial values and limits when fitting R with Lorentzian
par_init   = [30,400,50,0,0]                
par_limits = [[0,0,0,0,0],[255,1000,100,10,30]]

# best parameters initial values and limits when fitting G with Lorentzian
par_init   = [30,400,50,0,0]                
par_limits = [[0,0,0,0,0],[255,1000,100,10,30]] # saturation correction is not necessarily needed...

# best parameters initial values and limits when fitting B with Lorentzian
par_init   = [30,400,50,0,30] 
par_limits = [[0,0,0,0,20],[255,1000,100,10,40]] 

# best parameters initial values and limits when fitting S with Lorentzian
par_init   = [40,400,130,0,0]                
par_limits = [[0,0,0,0,0],[255,1000,1000,3,30]]

-Python

Copyright© ExpPhys_hacks , 2021 All Rights Reserved Powered by AFFINGER5.