忍者ブログ
物理学者(ポスドク)による日々の研究生活のメモ書きです ( python/emacs/html/Japascript/シェルスクリプト/TeX/Mac/C言語/Linux/git/tmux/R/ポケモンGO)
import statsmodels.api as sm
a = sm.tsa.acf(ts, nlags=n
nlagsは書かないと、すごく少ない数しか結果が出てこないので必須
nlagsの値はデータのサンプル数で置き換えてOK


gwpyでプロットしたいときは、こんな感じ
import statsmodels.api as sm
from gwpy.timeseries import TimeSeries

a = sm.tsa.acf(ts, nlags=int(T*fs))
a = TimeSeries(a, sample_rate=fs)
a.plot()

自己相関関数のプロットは、コレグラムとか言うらしい
へぇ〜

https://qiita.com/hadacchi/items/ea61d6d714c5a7481461

の計算速度の比較を見ると、自作するよりもそれなりに早いらしい






PR
img_new = cv2.addWeighted(src1=img1, alpha=0.5, src2=img2, beta=0.5, gamma=0)

img1とimg2を透明度0.5ずつにして足し合わせる

もしimg1の明るくて、img2が暗い場合はimg1ばかりが見えてしまうので比率を0.5じゃなくて3と0.1とかにしないといけない場合もあるので要調整

gammaはオフセットみたいなもの
とりあえず0でOK





cv2.resize(img, Size,fx=0,fy=0,interpolation=INTER_LINEAR)

でresizeを行う
Sizeは(width, height)の順番

fx, fyは拡大率とかだけど、比率を同じにしたままresizeするならなくて良い





いつか使うかもしれないのでメモ

これににたので、アフィン変換というのもある
これは、平行移動とある点を中心にした回転を同時に行える変換らしい

■ 参考 : サイズ変更・回転・切り抜き

■ 参考 : Python, OpenCVで幾何変換(アフィン変換・射影変換など)

2枚の画像を比較する際に、片方はdtype=np.uint8、もう片方がdtype=np.uint16だと比較できないと怒られる
ダイナミックレンジが異なる

こういうときは、小さい方のダイナミックレンジの画像を大きい方に合わせる

img_oldがnp.uint8だとすると、
new_fig = np.uint16(img_old)

で OK




2枚の画像を比較する際に、解像度が違うことがある。
アスペクト比が同じであれば、単純にresizeして縮小・拡大すればいいが、アスペクト比が違う場合は余白を追加してそれを揃える必要がある
import cv2

img_orig = cv2.imread(fname, cv2.IMREAD_UNCHANGED)
img_gige = img_orig.copy()
height_gige, width_gige = img_orig.shape
fnameというファイル名の画像は読み込む
ここで扱う画像は、白黒だと考えている。カラーの場合は最後のshapeはもう1つ要素が増えているはず

これで縦横の大きさを取得する


img_gigeという画像に上下にmargin_height、左右にmargin_widthを追加する。
cv2.copyMakeBorder(src, top, bottom, left, right, borderType, value)
というのが使い方

■ 参考 : OpenCVでzero-paddingを1行でする方法


valueというのは追加する余白の色、ここでは黒なので[0, 0, 0]にしているが白とか赤にすることもできる
# 余白を上下に追加する

img_gige2=cv2.copyMakeBorder(img_gige, margin_height, margin_height, margin_width, margin_width, cv2.BORDER_CONSTANT, value=[0, 0, 0])







condaよりも圧倒的に早い・・・

まずはcondaのアンインストール
anaconda-cleanというパッケージをインストールしてあれこれやっているけど、関係してそうなディレクトリをそのまま削除するのでも良いと思う・・・

■ 参考 : Anaconda3をきれいにアンインストールする

■ 参考 : Mambaを使った高速condaパッケージ管理 (python)

■ 参考 : condaの代わりに高速なmambaを使う


condaでできることは、(自分が知る限りでは)mambaでもできる
そのまま
conda install hoge

mamba install hoge
と置き換えるだけでOK








png画像の一部を楕円でフィットする要請があったので色々と調べた時のメモ

結局opencvで、画像の上に楕円を書いて、一番あっているものを採用することにした
(要するに手でフィットした)

楕円と一言で言っても回転している可能性がたかいので、それも踏まえて解析しないといけない・・・
関数形どうなるんだよ・・・・・




回転していない楕円のパラメーターは、中心(x0, y0)と長軸、短軸のパラメーター2つ(a, b)と定数c
の合計5つある

なので、楕円の円周上の5点をサンプルできれば、それを解けば上記の5つのパラメーターがわかるのかなと思った・・・・

pythonで、数式を解いてくれるパッケージがsympyというのがある
使い方などは以下を参照

■ 参考 : Python, SymPyの使い方(因数分解、方程式、微分積分など)

■ 参考 : Python (SymPy) で方程式・連立方程式を解く、数列を求める

■ 参考 : 非線形方程式をpythonで解く




開店した楕円の式とか、関数形とか

■ 参考 : 複素数平面を利用して回転した楕円の方程式を得る

■ 参考 : 2次曲線の回転移動、標準化、判別式




5点から楕円のパラメーターを求めるページ

■ 参考 : 5点を通る楕円の方程式を求める方法

■ 参考 : 5点を通る楕円の方程式を求め標準形を計算し概形より描画




opencvで楕円を描く方法

■ 参考 : 【Python】OpenCVでの図形や文字列の描画まとめ(四角形・線分・矢印・円・楕円・マーク・多角形)


■ 参考 : OpenCVとPILで図形を描写する


■ 参考 : 【Python】OpenCVでの図形や文字列の描画まとめ(四角形・線分・矢印・円・楕円・マーク・多角形)



hist, bins = np.histogram(input_data, nbin)
delta = (bins[1]-bins[0])/2
xr = np.linspace(np.min(bins)+delta,np.max(bins)-delta,len(bins)-1)
x_max = xr[np.argmax(hist)]
input_data というのがヒストグラムを書くデータ
x_maxが最頻値
ただし、ヒストグラムなので、階段状になっている。ここのx_maxはその階段の真ん中の値になっていると思う

■ 参考 : numpyでヒストグラムの最大地点を取得する関数




自分のスクリプトが現状どの程度完了しているかのために、プログレスバーを表示できる

これをcgi-binと組み合わせて使いたかったけど、それはできないっぽいので諦めた
ひとまず、どこかで使うかもしれないのでメモしておく


■ 参考 : Python の progress bar いろいろ

■ 参考 : Pythonで進捗表示したい!

■ 参考 : Pythonで進捗可視化のためにtqdmを使う


これが一番楽そう



自作で、プログレスバーを表示させたいときは

■ 参考 : 改行しない文字列出力
# print関数の引数endを指定する
print('文字列',end='')

# 先頭にキャリッジリターンをつけると、上書きできる
print('\r文字列',end='')


へぇ〜〜〜
これは便利そう








numpyのload("hoge.npz")とかでデータを読み込んだ時に、小数点以下の桁に0じゃない値が出てきて困った・・・

その値を他のものと比較する必要があったので、
こちらで丸めて、対応した

import numpy as np
b = np.around(a, 2)

でaを小数点以下2桁まで丸める
小数点以下じゃなくて、整数部分を何桁かで丸めたいときは、-2とかを入れればいい

他にもいろんな種類の丸めがnumpyには実装されてる↓

■ 参考 : numpy – 浮動小数点に関係する関数について

■ 参考 : NumPy配列ndarrayの表示形式(桁数や指数表記、0埋めなど)を指定






opencvでできる

まずは、うまくいった例をメモ


■ 参考 : OpenCVをつかった特徴点マッチングについて少しだけ掘り下げる

基本的に↑の記事通りに動かした

■ 参考 : OpenCV3とPython3で特徴点を抽出する(AgastFeature, FAST, GFTT, MSER, AKAZE, BRISK, KAZE, ORB, SimpleBlob, SIFT)
import cv2
from IPython.display import Image
from IPython.display import display
from matplotlib import pyplot as plt

#
# ここで2枚の画像読み込み
#
# (お好きな画像を読み込んでください)

akaze = cv2.AKAZE_create()
kp1, des1 = akaze.detectAndCompute(img1, None)
kp2, des2 = akaze.detectAndCompute(img2, None)

#
# 内容のチェック
#
print('##### 特徴点の数 #####')
print(len(kp1))

print('##### 特徴量記述子 #####')
print(des1)

print('##### 特徴ベクトル #####')
print(des1.shape)

#
# マッチング
#
bf = cv2.BFMatcher(cv2.NORM_HAMMING, crossCheck=True)
matches = bf.match(des1, des2)

#
# 特徴点間のハミング距離でソート
#
matches = sorted(matches, key=lambda x: x.distance)

#
# 2画像間のマッチング結果画像を作成
#
img1_2 = cv2.drawMatches(img1, kp1, img2, kp2, matches[:40], None, flags=2)
plt.figure(figsize = (width*2/100, height/100), dpi=100)
plt.imshow(img1_2)
plt.xticks([]), plt.yticks([]) # to hide tick values on X and Y axis
plt.savefig("diff_features.png")
その後で、2枚の画像のうち、どことどこが一致しているかを調べた

特徴量のうち、一致度上位n個のみをプロットする
距離が遠いものはNoneとして無視する(これであってるのかわからんが・・・)
import matplotlib as mpl
mpl.rcParams['xtick.labelsize'] = 16
mpl.rcParams['ytick.labelsize'] = 16
mpl.rcParams["axes.labelsize"] = 20
mpl.rcParams["axes.titlesize"] = 20
mpl.rcParams["legend.fontsize"] = 20
mpl.rcParams['axes.linewidth'] = 2
mpl.rcParams['figure.facecolor'] = "white"
mpl.rcParams['mathtext.fontset'] = 'stix'
mpl.rcParams['font.family'] = 'STIXGeneral'

import numpy as np

n=100
x_diff = np.zeros(n)
y_diff = np.zeros(n)
for i in range(n):
x = kp1[matches[i].queryIdx].pt[0] - kp2[matches[i].trainIdx].pt[0]
if abs(x) > 20:
x_diff[i] = None
else:
x_diff[i] = x

x = kp1[matches[i].queryIdx].pt[1] - kp2[matches[i].trainIdx].pt[1]
if abs(x) > 20:
y_diff[i] = None
else:
y_diff[i] = x

fig = plt.figure(figsize=[12, 12])

x_diff = x_diff[~np.isnan(x_diff)]
y_diff = y_diff[~np.isnan(y_diff)]

plt.subplot(2, 1, 1)
plt.title("Comparison %s\nand %s" % (fname1, fname2))
plt.hist(x_diff, bins=20, label="horizontal difference\nmean=%.1f pixel\nstd=%.1f pixel" % (np.mean(x_diff), (np.std(x_diff))))
plt.legend(loc="upper right")
plt.xlabel("horizontal difference on features (x-value)")
plt.grid(linestyle='dotted', linewidth=1)

plt.subplot(2, 1, 2)
plt.hist(y_diff, bins=20, label="vertical difference\nmean=%.1f pixel\nstd=%.1f pixel" % (np.mean(y_diff), (np.std(y_diff))))
plt.legend(loc="upper right")
plt.xlabel("vertical difference on features (y-value)")
plt.grid(linestyle='dotted', linewidth=1)

fig.savefig("hist.png")



検索ワード 「python opencv 位置合わせ qiita」とか

■ 参考 : [OpenCV] いまさら局所特徴量で物体検出!?

■ 参考 : OpenCVで画像マッチングをする

■ 参考 : python+opencvで画像処理の勉強8 パターン・図形・特徴の検出とマッチング

■ 参考 : 2つの画像を比較して違いを見つける


■ 参考 : 特徴点のマッチング





(2023/10/06 追記)

■ 参考 : OpenCV-Pythonチュートリアル

opencvで画像解析する人へのチュートリアル記事





(2023/12/08 追記)

■ 参考 : 特徴点のマッチングと対応点の座標のCSV出力

■ 参考 : OpenCVで特徴量の座標を取得する

■ 参考 : 特徴量マッチングによるテンプレートマッチング

■ 参考 : Python+OpenCVを利用したマッチング処理

■ 参考 : 特徴マッチングによる物体検知






適当にグラフを書きたい時はsubplotsよりもsubplotのほうが簡単に使えそうな気がする・・・・

■ 参考 : subplot? add_subplot? subplots?

■ 参考 : matplotlibの描画の基本 - figやらaxesやらがよくわからなくなった人向け


from matplotlib import pyplot as plt

plt.subplot(2, 1, 1)
plt.hist(x, bins=20, label="horizontal")
plt.legend(loc="upper right")
plt.xlabel("horizontal")
plt.grid(linestyle='dotted', linewidth=1)

plt.subplot(2, 1, 2)
plt.hist(y, bins=20, label="horizontal")
plt.legend(loc="upper right")
plt.xlabel("vertical")
plt.grid(linestyle='dotted', linewidth=1)

fig.savefig("hist.png")
グリッドは
plt.grid(linestyle='dotted', c='purple', linewidth=3)
のような感じ





np.nanじゃなくて、pythonのNoneに対しても使えたのでメモ


■ 参考 : NumPy配列ndarrayの欠損値np.nanを含む行や列を削除



xというarrayに対して、Noneだけを取り除きたい時
x = x[~np.isnan(x)]







それをsvgで保存して、legendをクリックしたら、そのデータが画像から消えるようなパッケージを発見したので利用していた

まず、matplotlibでプロットを作成する
そして、すでに作成した fig = plt.figure(figzsize=[12, 8])の中からデータの部分を抜き出して使っているらしい

保存する時に
f = BytesIO()
plt.savefig(f, format="svg")
なることが行われているが、これがよくわからない・・・
このときに、元々のfigのサイズが消えてしまうため、pngでsavefigしたときとsvgでこのパッケージを使ってsavefigしたときで大きさが違ってしまう


なので、svgで保存する時にfigsizeを手で設定する
plt.gcf().set_size_inches(18, 8)
plt.tight_layout()
最初は[12, 8]なのに、後のところは[18, 8]にすると大体同じ大きさになった
空白の大きさ? よくわからん








そのまま、不等号とかで大小の比較ができる
from datetime import datetime
date_beg = datetime(2023, 8, 20)
date_end = datetime(2023, 8, 28)
if date_beg < date_end:
 print("True")
みたいな感じ





AttributeError: module 'numpy' has no attribute 'int'.
`np.int` was a deprecated alias for the builtin `int`. To avoid this error in existing code, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use
eck the release note link for additional information.
The aliases was originally deprecated in NumPy 1.20; for more details and guidance see the original release note at:
https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations


np.intの代わりにintを使えばいいっぽい

Numpyの1.20では、numpy.intは非推奨になったとのこと






import cv2

photo = cv2.imread('hoge.png')
grays = cv2.cvtColor(photo, cv2.COLOR_RGB2GRAY)
rgb = cv2.cvtColor(grays, cv2.COLOR_GRAY2RGB)
みたいな感じ。

グレースケールの明るさがそのままRGBのそれぞれに割り当てられるらしい
一度グレースケールにしたら、元のRGBにはそのままでは戻せないので注意

自分がやりたかったのは、グレースケールの画像を使って(0, 0, 0)と(255, 255, 255)の2色のフィルターみたいなものを作りたかったけど、
グレースケールの画像だと 0〜255の1つの値しかなくて、RGBの3色の画像と重ねたりはできない
なので、グレースケールの画像をRGBの3つの値に変換したかった



■ 真っ黒な画像を用意する方法
import cv2
import numpy as np

size=(512,512)
black_img=np.zeros(size,np.uint8)

みたいな感じ
これで全部のピクセルに0が代入されている。


これと↑のグレースケールをRGBに変換を組み合わせると、RGBで真っ黒な画像が作れる





plt.grid(which='major',color='black',linestyle='-')
plt.grid(which='minor',color='black',linestyle='-')
■ 参考 : 初心者メモ:PythonのList/Tuple/Set/Dictのチートシート

毎回、あれどうやるんだっけ・・・って調べてたものが全部一覧されとる・・・・

(これのnumpy array版もほしい)







■ 参考 : 知らなかったPython3の言語仕様まとめ


■ 参考 : Pythonでflatten(多次元リストを一次元に平坦化)

sumでもできるらしい
importしなくても使えるのが便利
l_2d = [[0, 1], [2, 3]]

print(sum(l_2d, []))
# [0, 1, 2, 3]



evalとexecというのが見つかった

execはそのまま式が実行される

evalは式を評価する、その返り値は取り出せる
func="1+1"
a = eval(func)

func="x=5"
exec(func)

みたいな感じ
aには1+1の結果の2が代入されてる
execで実行されたので、xには5が代入されてる






関数の外で定義した変数を、関数の中でいじったら関数外の変数の値も変わると思ってたがそうではないらしい・・・・

関数の外で定義した変数はグローバル変数と呼ばれる
それを関数内で呼び出すときは、
global hoge
みたいな感じで頭にglobalをつければOK

関数の内外で値を共有する方法として一番簡単な方法は、
def func(hoge):
 # なんかする
 return hoge
で返り値にしてしまう方法

今回はこれではちょっと都合が悪かったのでglobalを使った





(2023/06/09 追記)

mainのコードとmodule内で共通の変数を使いたい、かつ関数の引数として渡さない方法もあるっぽい

■ 参考 : Pythonで全モジュール共通のグローバル変数を扱う方法

モジュール内で空のmoduleを読み込む
mainのコード内でも同じ空のmoduleを読み込む
モジュールの関数内で、空moduleの変数として定義すればmain関数でもそれを使えるようになる



他にも
■ 参考 : Global変数をmodule間やmainとmoduleの間で共有する方法

というのがあるっぽい








jsonファイルを作るときに、
辞書型のkeyをシングルコーテーションで括っているとこのエラーが表示される
(一応jsonファイル自体はできていたが・・・・)

■ 参考 : Python JSONデータ読み込み時にエラーが発生する (Expecting property name enclosed in double quotes)








リポジトリ名とパッケージ名が同じ場合は・・・・
pip install git+https://github.com/acc_name/pkg_name.git




comments="%"でOK
import numpy as np
a = np.loadtxt(fname, delimiter=",", unpack=True, comments="%")





プロフィール
HN:coffee
職業:物理屋(自称)
趣味:映画鑑賞、登山
出身:大阪府の南の田舎
自己紹介:
import MyProfile
import coffee_pote from TWITTER
import amazonのほしい物リスト from WISH_LIST

print "先月子供が産まれました!"

# 最終更新 2022/10/25
カウンター
カウンター カウンター
ブログ内検索
ツイートするボタン
リンク
相互リンク募集中です (Twitterにてお知らせください)

Demo scripts for gnuplot version 5
(gnuplotのさまざまなデモ画像と作り方がまとめられている、眺めているだけでできるようになった気分になれる)

gnuplotスクリプトの解説
(米澤進吾さんの個人ページ、gnuplotと言えばこのかた)

gnuplot のページ
(Takeno Lab、うちのブログがリンクされていたのでリンク返し)

Twitterから映画の評価が分かる & 映画の鑑賞記録が残せる coco
(映画の感想をまとめられるサイト、いつもお世話になっています)

Astronomy Picture of the Day Archive
(天文や宇宙関連の最新の話題について画像とともにNASAが説明しているページ)

今日のほしぞら
(任意の時刻の空で見える星を表示してくれる、国立天文台が管理している)

GNUPLOTとアニメーション
(応用の項目の「見せてあげよう!ラピュタの雷を!!」あたりからすごすぎる)

読書メーター
(読んだ本をリストできる便利なサイト)

flickr難民の写真置き場
(20XX年、flickrは有料化の炎に包まれた。あらゆるflickr無料ユーザーは絶滅したかに見えた。 しかし、tumblr移住民は死に絶えてはいなかった。)

教授でもできるMac OS X へのLaTeX, X11, gccのインストレーションと環境設定
(阪大の山中卓さんのwebページ、タイトルにセンスが溢れている、内容は超充実してる、特にTeX関連、学振DCとかPDの申請書類作成時にはお世話になっております)

英語論文執筆用の例文検索サービス
(とんでもないものを見つけてしまった・・・・ arXivに収録されている 811,761報の 英語論文から,例文を検索するための検索エンジン)


Template "simple02" by Emile*Emilie
忍者ブログ [PR]