pythonでfor文を使わずに任意の値でconvolutionを行いたい．

Question

pythonにおいてfor文を使わずにconvolutionをしたいと考えています．

画像サイズを64*64としています．

行いたい処理として，3*3の注目画素範囲に対して，
その中心画素（convolution[4]）と
周辺８近傍（convolution[4]-convolution[0~8]: 中心画素同士の減算はしない）
との差を２乗したものを，
その中心画素からの画素の距離に基づいた重みづけをして加算したものを新たな画素値になる処理です．
（中心画素に対して上下左右なら重みは1,斜め方向なら√2となります．）
scipyなど調べましたが，何を使用したらよいかわかりませんでした．
ご教授のほどよろしくお願いいたします．

import numpy as np

y_size=64
x_size=64
x=np.empty(3)
y=np.empty(3)

original_image=np.arange(64*64)
convolution=np.empty(9)
last_image=np.empty(64*64,dtype="double") 

for i in range(y_size):
      
        for j in range(x_size):
            
            x[0]=(j+x_size-1) % x_size
            x[1]=j
            x[2]=(j+1) % x_size
                        
            y[0]=(i+y_size-1) % y_size
            y[1]=i
            y[2]=(i+1) % y_size
                
            for m in range(9):
                #配列convolutionに3*3の画素値を代入
                A=x[int(np.floor(m%3))]
                B=y[int(np.floor((m/3)))]*x_size
                convolution[m]=original_image[int(B)+int(A)]
            #中心画素とその8近傍の画素値の差を２乗して，画素間の距離の重みづけ    #改行
            value=np.power((convolution[4]-convolution[0]),2.0)/np.sqrt(2.0) \
               +  np.power((convolution[4]-convolution[1]),2.0)              \
               +  np.power((convolution[4]-convolution[2]),2.0)/np.sqrt(2.0) \
               +  np.power((convolution[4]-convolution[3]),2.0)              \
               +  np.power((convolution[4]-convolution[5]),2.0)              \
               +  np.power((convolution[4]-convolution[6]),2.0)/np.sqrt(2.0) \
               +  np.power((convolution[4]-convolution[7]),2.0)              \
               +  np.power((convolution[4]-convolution[8]),2.0)/np.sqrt(2.0) \
               
            last_image[i*x_size+j]=value   

print(last_image)
#後で64*64にreshapeをする

Answer 1

以下でどうでしょう

original_image = np.arange(64*64).reshape(64, 64)

pad_image = np.pad(original_image, (1, 1), 'wrap')

f = np.array([[1/2**0.5, 1, 1/2**0.5],
              [1, 0, 1],
              [1/2**0.5, 1, 1/2**0.5]])
arr_for_conv = np.lib.stride_tricks.sliding_window_view(pad_image, (3, 3))
last_image = np.sum((arr_for_conv - original_image[:, :, None, None])**2*f, (-2, -1))

---

手順

配列の端の値を畳み込み演算するため、まずnp.pad()で外周を埋めます。質問文のコードを見る限り'wrap'モードでのパディングと思われます。

original_image = np.arange(64*64).reshape(64, 64)  # オリジナル画像、形は64*64
# array([[   0,    1,    2, ...,   61,   62,   63],
#        [  64,   65,   66, ...,  125,  126,  127],
#        [ 128,  129,  130, ...,  189,  190,  191],
#        ...,
#        [3904, 3905, 3906, ..., 3965, 3966, 3967],
#        [3968, 3969, 3970, ..., 4029, 4030, 4031],
#        [4032, 4033, 4034, ..., 4093, 4094, 4095]])

pad_image = np.pad(original_image, (1, 1), 'wrap')  # 形は66*66
# array([[4095, 4032, 4033, ..., 4094, 4095, 4032],
#        [  63,    0,    1, ...,   62,   63,    0],
#        [ 127,   64,   65, ...,  126,  127,   64],
#        ...,
#        [4031, 3968, 3969, ..., 4030, 4031, 3968],
#        [4095, 4032, 4033, ..., 4094, 4095, 4032],
#        [  63,    0,    1, ...,   62,   63,    0]])

np.lib.stride_tricks.sliding_window_view()を使って3*3のブロックに変換します。

arr_for_conv = np.lib.stride_tricks.sliding_window_view(pad_image, (3, 3))  # 形は64*64*3*3
# array([[[[4095, 4032, 4033],
#          [  63,    0,    1],
#          [ 127,   64,   65]],  # この3*3それぞれが質問文のコードのconvolutionに相当
# 
#         [[4032, 4033, 4034],
#          [   0,    1,    2],
#          [  64,   65,   66]],
# 
#         [[4033, 4034, 4035],
#          [   1,    2,    3],
#          [  65,   66,   67]],
# 
#         ...,

ブロードキャストを用いて、中心との差を計算し、

arr_for_conv - original_image[:, :, None, None]

それを２乗、

(arr_for_conv - original_image[:, :, None, None]) ** 2

さらに重み付け

# 重み付け用のフィルター
f = np.array([[1/2**0.5, 1, 1/2**0.5],
              [1, 0, 1],
              [1/2**0.5, 1, 1/2**0.5]])
(arr_for_conv - original_image[:, :, None, None]) ** 2 * f  # 形は64*64*3*3

その3*3の部分をnp.sum()で合計

last_image = np.sum((arr_for_conv - original_image[:, :, None, None])**2*f, (-2, -1))  # 形は64*64

# array([[3.96381285e+07, 3.92578213e+07, 3.92578213e+07, ...,
#         3.92578213e+07, 3.92578213e+07, 3.88966733e+07],
#        [2.93616653e+04, 1.97820659e+04, 1.97820659e+04, ...,
#         1.97820659e+04, 1.97820659e+04, 2.93616653e+04],
#        [2.93616653e+04, 1.97820659e+04, 1.97820659e+04, ...,
#         1.97820659e+04, 1.97820659e+04, 2.93616653e+04],
#        ...,
#        [2.93616653e+04, 1.97820659e+04, 1.97820659e+04, ...,
#         1.97820659e+04, 1.97820659e+04, 2.93616653e+04],
#        [2.93616653e+04, 1.97820659e+04, 1.97820659e+04, ...,
#         1.97820659e+04, 1.97820659e+04, 2.93616653e+04],
#        [3.88966733e+07, 3.92578213e+07, 3.92578213e+07, ...,
#         3.92578213e+07, 3.92578213e+07, 3.96381285e+07]])

Comments

[masaosa]

非常にわかりやすい説明ありがとうございます！
ちなみにsliding_window_viewはnumpy1.20の最新のバージョンからの関数でしょうか？

上記のコードで行い，値としてはほぼ間違いなかったのですが，私のコードのnp.sqrt(2.0）で除算する場合と，kirara0048様の1/2**0.5で乗算するところで誤差（非常に小さい値なのですが，，，）が出る箇所がありました．これはやはり，小数点の乗算と除算での誤差によるものということでしょうか？

[kirara0048]

>sliding_window_viewはnumpy1.20の最新のバージョンからの関数でしょうか？
そうです。なお、scikit-imageをインストールしていればskimage.util.view_as_windows()という関数で同じことができます。

>小数点の乗算と除算での誤差によるものということでしょうか？
まあそうですね。
`np.allclose(配列1, 配列2)`で小数点の誤差を気にせず２つの配列が同じかどうか比較することができます。

[masaosa]

ありがとうございました！非常に処理が速くなりました！

Answer 2

ネットで調べると、

scipy.signal.correlateが使えると書いてありました。
5年前の記事絵古いかもしれませんが、ご参考まで。

Pythonのconvolve2dについて

Comments

[masaosa]

ご回答ありがとうございます！
任意のフィルター係数のコンボリューションのやり方は上記の記事で理解しました．

フィルターの係数を決める際，中心画素とその８近傍の画素値の差を係数にしたい場合，for文なしで書くことはできますでしょうか？