やりたいこと

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

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

行動規範の内容に同意します

回答1件

ベストアンサー

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

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