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

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

ただいまの
回答率

88.10%

python スペクトル分布

解決済

回答 1

投稿

  • 評価
  • クリップ 0
  • VIEW 3,964

score 19

 やりたいこと

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

 やったこと

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

def main():    
    # 入力画像をグレースケールで読み込み
    img_orig = cv2.imread('ceiling2.orig.png')
    img_recon = cv2.imread('ceiling2.recon150.png')
    # グレースケール変換
    gray = cv2.cvtColor(img_orig, cv2.COLOR_BGR2GRAY)
    grayr = cv2.cvtColor(img_recon, cv2.COLOR_BGR2GRAY)
    h, w= grayr.shape
    # 高速フーリエ変換(2次元)
    fimg = np.fft.fft2(gray)
    fimgr = np.fft.fft2(grayr)
    # 第1象限と第3象限, 第2象限と第4象限を入れ替え
    fimg =  np.fft.fftshift(fimg)
    fimgr =  np.fft.fftshift(fimgr)
    #振幅スペクトルを計算
    amp = np.abs(fimg)
    ampr = np.abs(fimgr)
    # 位相
    phase = np.angle(fimgr)
    # パワースペクトルの計算
    mag = 20*np.log(np.abs(fimg))
    magr = 20*np.log(np.abs(fimgr))    
    # 入力画像とスペクトル画像をグラフ描画
    plt.subplot(121)
    plt.imshow(mag, cmap = 'gray')
    plt.subplot(122)
    plt.imshow(magr, cmap = 'gray')
    plt.show()
    # 表示
    plt.plot(phase, ampr)
    plt.grid()
    plt.show()
    plt.plot(amp)
    plt.grid()


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

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

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

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

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 過去に投稿した質問と同じ内容の質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

回答 1

checkベストアンサー

+2

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

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などの関数があります。これらの機能を利用すると次のように集計できます。(本件では同じ波長範囲を集計すればよく、コードでは波長ではなく周波数で範囲をチェックしています)

import numpy as np


# 複素平面を構成する値の配列を生成(原点は中心)
def complex_plane(width):
    half_width = width // 2
    re = np.array(range(width)) - half_width
    im = - re
    re, im = np.meshgrid(re, im)
    return re + im * 1j


# スペクトルを与えられた角度の範囲内だけ集計
def aggregate_in_angle(agg_fun, spectrum, min_angle, max_angle):
    width = spectrum.shape[0]
    assert spectrum.shape[1] == width
    min_radius, max_radius = 0, width // 2
    assert 0 <= min_angle < max_angle <= 180  # 単位はdegree

    cp = complex_plane(width)

    cp_mag = np.abs(cp)
    in_radius = np.logical_and(min_radius < cp_mag, cp_mag < max_radius)

    cp_angle = np.angle(cp, deg=True)
    in_angle = np.logical_and(min_angle <= cp_angle, cp_angle < max_angle)

    return agg_fun(spectrum[np.logical_and(in_radius, in_angle)])


def test(spectrum, angle_gap):
    # 指定角度ごとの平均を求める
    means = [aggregate_in_angle(np.mean, spectrum, angle, angle + angle_gap)
             for angle in range(0, 180, angle_gap)]
    return means

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

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

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

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

イメージ説明

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

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2018/11/21 16:45

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

    キャンセル

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

  • ただいまの回答率 88.10%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

関連した質問

同じタグがついた質問を見る