質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.50%
Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

Q&A

解決済

1回答

10165閲覧

python スペクトル分布

homel15k

総合スコア19

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

0グッド

0クリップ

投稿2018/11/12 02:04

やりたいこと

pythonでパワースペクトル画像の特徴を分析したいと思っています.
具体的には画像上の角度方向のスペクトル分布を求めてピークの位置を求めたいです.
イメージ説明

やったこと

二次元フーリエ変換でパワースペクトルまでは求めました。

python3

1def main(): 2 # 入力画像をグレースケールで読み込み 3 img_orig = cv2.imread('ceiling2.orig.png') 4 img_recon = cv2.imread('ceiling2.recon150.png') 5 # グレースケール変換 6 gray = cv2.cvtColor(img_orig, cv2.COLOR_BGR2GRAY) 7 grayr = cv2.cvtColor(img_recon, cv2.COLOR_BGR2GRAY) 8 h, w= grayr.shape 9 # 高速フーリエ変換(2次元) 10 fimg = np.fft.fft2(gray) 11 fimgr = np.fft.fft2(grayr) 12 # 第1象限と第3象限, 第2象限と第4象限を入れ替え 13 fimg = np.fft.fftshift(fimg) 14 fimgr = np.fft.fftshift(fimgr) 15 #振幅スペクトルを計算 16 amp = np.abs(fimg) 17 ampr = np.abs(fimgr) 18 # 位相 19 phase = np.angle(fimgr) 20 # パワースペクトルの計算 21 mag = 20*np.log(np.abs(fimg)) 22 magr = 20*np.log(np.abs(fimgr)) 23 # 入力画像とスペクトル画像をグラフ描画 24 plt.subplot(121) 25 plt.imshow(mag, cmap = 'gray') 26 plt.subplot(122) 27 plt.imshow(magr, cmap = 'gray') 28 plt.show() 29 # 表示 30 plt.plot(phase, ampr) 31 plt.grid() 32 plt.show() 33 plt.plot(amp) 34 plt.grid()

パワースペクトルの角度方向はどのように抽出すればよいでしょうか?
アドバイスよろしくお願いします.

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

guest

回答1

0

ベストアンサー

ご質問の図からの想像ですが、多分以下のレポートを参照しておられるのではないでしょうか?

http://www.edu.kobe-u.ac.jp/eng-arch-sled/dat/research/adachi-kigami-1/1.pdf

この中に

q(θ)は水平軸から角度θの微小な線形領域内のパワースペクトルの和を表す。

とありますので、DFTの結果の中からq(θ)の範囲を示す扇型の内側の点のみを選択的に集計すればよいと思います。

以下、煩雑さを避けるために解析対象の画像が正方形であり、垂直方向の符号は数学でおなじみの座標系(複素平面)に合わせる前提とします。

画像の辺の画素数をWとすると、dftshift後のスペクトル配列spectrumの要素spectrum[W // 2 + re, W // 2 - im]は、複素数 c = re + 1j * im を用いて以下のような波の成分(スペクトル)となります。

・方向: cの偏角
・波長: W/|c| (単位は画素)

numpyには複素数の偏角を求めるangle、複素数の絶対値(複素平面原点からの距離)を求めるabsなどの関数があります。これらの機能を利用すると次のように集計できます。(本件では同じ波長範囲を集計すればよく、コードでは波長ではなく周波数で範囲をチェックしています)

python

1import numpy as np 2 3 4# 複素平面を構成する値の配列を生成(原点は中心) 5def complex_plane(width): 6 half_width = width // 2 7 re = np.array(range(width)) - half_width 8 im = - re 9 re, im = np.meshgrid(re, im) 10 return re + im * 1j 11 12 13# スペクトルを与えられた角度の範囲内だけ集計 14def aggregate_in_angle(agg_fun, spectrum, min_angle, max_angle): 15 width = spectrum.shape[0] 16 assert spectrum.shape[1] == width 17 min_radius, max_radius = 0, width // 2 18 assert 0 <= min_angle < max_angle <= 180 # 単位はdegree 19 20 cp = complex_plane(width) 21 22 cp_mag = np.abs(cp) 23 in_radius = np.logical_and(min_radius < cp_mag, cp_mag < max_radius) 24 25 cp_angle = np.angle(cp, deg=True) 26 in_angle = np.logical_and(min_angle <= cp_angle, cp_angle < max_angle) 27 28 return agg_fun(spectrum[np.logical_and(in_radius, in_angle)]) 29 30 31def test(spectrum, angle_gap): 32 # 指定角度ごとの平均を求める 33 means = [aggregate_in_angle(np.mean, spectrum, angle, angle + angle_gap) 34 for angle in range(0, 180, angle_gap)] 35 return means

ちなみにレポートでは「パワースペクトルの和」と書いてありますが・・・

  • 集計方法

DFTの結果を集計するなら和ではなく平均の方がよい気がします。連続関数の積分ならまだしも離散的スペクトルを集計するのですから扇型の範囲に含まれるデータ数が角度によって増減することを考慮した方がよいのではないでしょうか。

  • パワースペクトル

質問者さんのコードにあるデシベル(?)換算は「ダイナミックレンジが大きなデータの傾向を目視で確認しやすくするための処理」であって「集計などのデータ処理をする前に施すべき変換ではない」気がします。

自然画像では結果を評価しにくいので単純なスペクトルを持つサンプル画像を生成して上記を試してみました。

イメージ説明

ちなみに、いろいろな波の成分で試してみたところ、fft結果のスペクトルの絶対値に対して平均や最大値をとるとそこそこの結果が出るように見えました。ただ、目的によっては集計方法を変えたほうがよいかも知れません。そのあたりの知識が不十分なので自分にははっきりしたことが言えませんが。

投稿2018/11/20 04:01

編集2018/11/21 12:16
KSwordOfHaste

総合スコア18392

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

homel15k

2018/11/21 07:45

ありがとうございます。とてもわかり易いです。
guest

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問