Home

画像のウェーブレット解析の試み

OpenCVとndarray型配列

OpenCVのインストール（Windows、Linux）

　ウェーブレット解析が金谷（2003；36刷2023）に分かり易く説明されていた。これは１変数の関数が扱われているが、２変数の関数（画像データ）の場合を、Nason (2008, sec. 2.13)を参考にして、以下のように考えた。

　２変量の関数において、領域の外では０であるものについて考える。さらに、各矩形領域の中では定数であるとする。すなわち、

である。画像はピクセル単位で値が与えられるので、各矩形内で定数であるという条件は本質的には制約ではない。

　いま、スケーリング関数とウェーブレット母関数を次のように置く。

　このとき、さらに次のように置く。

これらの関数、、、は互いに直交している。したがって、ノルムを１に調整した以下の関数は正規直交系を構成する。

　いま、

と置けば、

である。

　さらに、

と置けば、

　

となる。は、レベルの空間の直和になっている。

　いま、が、空間において、次式のように表されているとする。

　

関数は、０あるいは１を値とする関数であるので、はの関数値である。関数を、空間、、、の成分に分解するために、それらの正規直交基底との内積を求めると、以下の漸化式を得る。詳しくは、PDFファイルWaveAnal.pdfを参照されたい。

上記の漸化式（下降サンプリング）により、式（１）に基づいて、式（２）における係数から空間、、、における係数、、、を求めることができる。

ウェーブレット解析のプログラムをリスト１のように作成した。ファイルは、圧縮ファイルwa_files.zipにまとめたので、自由にダウンロード展開して利用されたい。

リスト１のプログラムはモノクロ画像対応であり、カラー画像を読み込んだ場合はモノクロ化されて処理が行われる。

カラー画像対応のプログラムも作成してみた。入力画像ファイルもサイズは自由とした。読み込み後、512×512に調整している。カラー画像対応のプログラムは、このウェブサイトの後半で説明している。

リスト１は、Windowsにおいては、ライブラリopencv-pythonとnumpyおよびmatplotlibがインストールされておれば実行できる。すなわち、

pip install opencv-python numpy matplotlib

とインストールしておけばよい。

リスト１　画像のウェーブレット解析

import cv2 # pip install opencv-python

import numpy as np

import matplotlib.pyplot as plt

img1 = cv2.imread('Sample_4.png') # 画像512×512ファイルの読み込み

plt.imshow(img1)

plt.show()

f_0 = np.sum(img1, axis=2)/3 # モノクロ化

#

#　　　ウェーブレット解析結果を辞書型として保存

#

WaveletAnal = {}

WaveletAnal['0'] = {'f0': f_0}

n = 9

N = 2**n

#

# レベルｊの解析結果の図示

#

def disp_map(j):

Cj = WaveletAnal[f'{j}']['f0'] # 係数C_h,k

DHj = WaveletAnal[f'{j}']['DHj']

DVj = WaveletAnal[f'{j}']['DVj']

DDj = WaveletAnal[f'{j}']['DDj']

#

# 画像を配列imgsに設定

# 大きさ512×512を４分割して使用

#

imgs = np.full((512, 512, 3), 0.5)

for h in range((2**(n-j))):

for k in range((2**(n-j))):

for ii in range((2**j)//2): # スケールに合わせて引き延ばし

for jj in range((2**j)//2):

for kk in range(3): # 赤、緑、青の３チャンネル

#

# C_hkの値

#

imgs[512//2 + ((2**j)//2)*h + ii]\

[512//2 + ((2**j)//2)*k + jj][kk] \

= Cj[h][k] / 255

#

# DH, DV, DDは、値０を黒に設定

#

imgs[512//2 + ((2**j)//2)*h + ii][((2**j)//2)*k + jj][kk] \

= 0.0

imgs[((2**j)//2)*h + ii][(512//2) + \

((2**j)//2)*k + jj][kk] \

= 0.0

imgs[((2**j)//2)*h + ii][((2**j)//2)*k + jj][kk] \

= 0.0

if DHj[h][k] > 0: # 正ならば緑に設定

imgs[((2**j)//2)*h + ii][(512//2) + \

((2**j)//2)*k + jj][1] \

= DHj[h][k] / 255

elif DHj[h][k] < 0: # 負ならば赤に設定

imgs[((2**j)//2)*h + ii][(512//2) +\

((2**j)//2)*k + jj][0] \

= -DHj[h][k] /255

if DVj[h][k] > 0:

imgs[512//2 + ((2**j)//2)*h + ii][((2**j)//2)*k + jj][1] \

= DVj[h][k] / 255

elif DVj[h][k] < 0:

imgs[512//2 + ((2**j)//2)*h + ii][((2**j)//2)*k + jj][0] \

= -DVj[h][k] /255

if DDj[h][k] > 0:

imgs[((2**j)//2)*h + ii][((2**j)//2)*k + jj][1] \

= DDj[h][k] / 255

elif DDj[h][k] < 0:

imgs[((2**j)//2)*h + ii][((2**j)//2)*k + jj][0] \

= -DDj[h][k] /255

plt.imshow(imgs) # 画像配列imgsの表示

plt.xticks([])

plt.yticks([])

plt.title(f'j = {j}')

plt.show()

for j in range(1,n):

# レベルｊの解析

print(f'Level {j} started.')

Cj_1 = WaveletAnal[f'{j-1}']['f0']

Cj = np.empty((2**(n-j), 2**(n-j)))

DHj = np.empty((2**(n-j), 2**(n-j)))

DVj = np.empty((2**(n-j), 2**(n-j)))

DDj = np.empty((2**(n-j), 2**(n-j)))

#

#　　　　下降サンプリング

#

for h in range(2**(n-j)):

for k in range(2**(n-j)):

Cj[h][k] = (Cj_1[2*h][2*k] + Cj_1[2*h+1][2*k] +

Cj_1[2*h][2*k+1] + Cj_1[2*h+1][2*k+1]) / 4

DHj[h][k] = (Cj_1[2*h][2*k] - Cj_1[2*h+1][2*k] +

Cj_1[2*h][2*k+1] - Cj_1[2*h+1][2*k+1]) / 4

DVj[h][k] = (Cj_1[2*h][2*k] + Cj_1[2*h+1][2*k] -

Cj_1[2*h][2*k+1] - Cj_1[2*h+1][2*k+1]) / 4

DDj[h][k] = (Cj_1[2*h][2*k] - Cj_1[2*h+1][2*k] -

Cj_1[2*h][2*k+1] + Cj_1[2*h+1][2*k+1]) / 4

WaveletAnal[f'{j}'] = {'f0': Cj}

WaveletAnal[f'{j}']['DHj'] = DHj

WaveletAnal[f'{j}']['DVj'] = DVj

WaveletAnal[f'{j}']['DDj'] = DDj

disp_map(j) # 解析結果の表示

while True:

print(f'0 < j < {n}')

j = int(input('j = '))

if j <=0 or j >= n:

break

disp_map(j)

リスト１を実行すると、画像ファイルSample_4.pngが読み込まれて表示される（図１）。画像サイズは512ピクセル×512ピクセルであり、である。画像ファイルを512ピクセル×512ピクセルのサイズに変換するサンプルスクリプトも本サイトの後半に用意している。

図１

図１のフォームを閉じると、レベル１のウェーブレット解析が始まる。レベル１の解析が終わると、解析結果が図示される（図２）。

図２

表示画面の右下にのマップが表示されている。右上は、ウェーブレットの垂直成分、左下はウェーブレットの水平成分、左上が対角成分である。プログラム中では、画像はnumpyのndarray型で表されているが、表示はOpneCVで行っている。OpenCVでの縦座標はndarray型の行に、横座標はndarray型の列に対応しているが、このことは本稿の最後の附記で確認している。

図２のフォームを閉じると、次のレベル、すなわちレベル２の解析が始まり、解析が終了すると、そお結果が表示される（図３）。

図３

図３のフォームを閉じると、次のレベル３の解析へと進む。この手順を繰り返すと、レベル８まで進む。画像の大きさが512ピクセル×512ピクセルのときは、であるので、レベル８の解析が最終レベルとなる（図４）。

図４

図４のフォームを閉じると、端末画面は以下のようになっている。

= RESTART: D:/WaveletAnal/wa_files/WaveAnal.py

Level 1 started.

Level 2 started.

Level 3 started.

Level 4 started.

Level 5 started.

Level 6 started.

Level 7 started.

Level 8 started.

0 < j < 9

j =

ｊに解析結果のレベル値ｊを設定すると、そのレベルの結果が表示される。例えば、５を設定すると、下図（図５）のようになる。

図５

図５のフォームを閉じて、ｊに範囲外の値、例えば、０を設定すると、プログラムの実行終了である。

0 < j < 9

j = 5

0 < j < 9

j = 0

附記　OpneCVの座標軸とnumpyのndarray型の行・列との関係

画像を表すNumpyのndarray型配列の行は、OpenCVの画像において縦軸に対応し、列はOpenCVの画像の横軸に対応している。このことを確認するスクリプトをリスト２のように用意した。

リスト２は、ライブラリopencv-pythonとnumpyおよびmatplotlibがインストールされておれば実行できる。

Windowsの場合のインストールは、例えば

pip install opencv-python numpy matplotlib

とインストールしておけばよい。

Importは

import cv2

である。

Linux（Ubuntu/CentOS Stream）の場合のOpenCV 2のインストールは、AnacondaのPython3.10の仮想環境において、

conda install opencv

でインストールできた。

ImportはWindowsの場合と同じく

import cv2

である。

リスト２　Numpyのndarray型配列とOpneCVの画像との対応

import cv2 # pip install opencv-python

import numpy as np

import matplotlib.pyplot as plt

img = np.full((300, 500, 3), 128, dtype = np.uint8)

#

# Numpyのndarray型配列の要素に値を設定

#

for i in range(130, 170):

for j in range(500):

img[i][j][1] = 255

#

# OpenCVの関数で描画

#

cv2.rectangle(img, (225, 0), (275, 300), (0, 0, 255), thickness = -1)

#

# 画像を表すndarray型配列の表示

#

plt.imshow(img)

plt.show()

リスト２を実行すると図Aのフォームが表示される。

図A

横の薄緑色の帯が、ndarray型配列imgの要素に直接値を設定して描いたものである。配列imgの行がOpenCVの画像の縦座標に、配列imgの列がOpenCVの画像の横軸に対応していることが分かる。青色の縦の帯は、OpenCVの関数rectangleで描いたものである。第１座標が横座標に、第２座標が縦座標に対応していることが分かる。

画像サイズの変換

画僧のサイズを512ピクセル×512ピクセルのサイズに変換するスクリプトをリスト３のように作成した。

リスト３　画像サイズの変換スクリプト

import cv2

import matplotlib.pyplot as plt

img = cv2.imread(input('File name = ')) # 画像の読み込み

plt.imshow(img)

plt.title('Read-in Image')

plt.show()

img1 = cv2.resize(img, dsize = (512, 512)) # 画像のリサイズ

plt.imshow(img1)

plt.title('Resized Image')

plt.show()

fn_img = input('Output file name = ')

cv2.imwrite(fn_img, img1) # 画像の保存

print(fn_img, 'was saved.')

リスト３のスクリプトを実行すると、入力画像ファイル名の設定が求められる。

File name =

サイズを変換したい画像ファイルの名前を設定する。

File name = Sample_p0.png

名前を設定すると読み込まれて、読み込まれた画像が表示される（図S.1）。

図S.1

フォーム右上角のX印をクリックして閉じると、サイズ変換後の画像が表示される。サイズの変換は、図S.1の座標軸の値と座右S.2の座標軸の値が異なることに反映されている。

図S.2

図S.2のフォームを閉じると、変換された画像を保存するファイル名の設定が求められる。

Output file name = Sample_rs.png

ファイル名を設定すると、画像が保存されてプログラムの終了である。

端末画面には、次のメッセージが表示される。

Sample_rs.png was saved.

保存された画像ファイルを、リスト１のスクリプトの入力画像の設定のコードを

img1 = cv2.imread('Sample_rs.png')

に変えて解析してみると以下のようになる。

・

・

・

引用文献

金谷健一（2003）これなら分かる応用数学教室．共立出版

Nason, G. P. (2008) Wavelet methods in statistics with R. Springer

なお、Pythonによるデータ分析の解説を＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で行っている。

カラー画像のウェーブレット解析

カラー画像のウェーブレット解析のスクリプト例をリストC1に示す。2次元ウェーブレット解析の基本的考え方はWaveAnal.pdfに示した。カラー画像のウェーブレット解析を、RGBの各チャンネルに対応して行った。入力画像の大きさは任意としたが、読み込み後、512ピクセル×512ピクセルに調整した。リストC1のスクリプトファイルとサンプル画像ファイルは、まとめて圧縮ファイルwa_color_files.zipとした。ダウンロード展開して自由に利用できる。

リストC1のスクリプトを次のコマンドで実行する。

PS ******\wa_color_files> python wa_color.py

入力画像ファイル名を聞いてくるので、設定する。

File name = sample_pctr.png

上の例では、sample_pctr.pngを設定している。入力画像ファイルは、スクリプトファイルと同じフォルダに置いておくと指定が簡単である。

入力画像ファイルを設定すると、読み込まれた画像が表示される（図C1）。

図C1

図C1のフォームを閉じると、512ピクセル×512ピクセルに調整された画像が表示される（図C2）。

図C2

画像の座標軸の値から、512ピクセル×512ピクセルに調整されていることが分かる。

図C2のフォームを閉じると、レベル１でのウェーブレット解析が始まり、終了すると解析結果が表示される（図C3）。

図C3

図C3のフォームを閉じると、レベル２の解析が行われ、解析結果が表示される（図C4）。

図C4

図C4のフォームを閉じると、次のレベルの解析が行われ、結果が表示される。このフォームを閉じて、次の解析の結果が表示されるという手順が繰り返されて、図C5のレベル８の結果の表示まで進む。

図C5

図C5のフォームを閉じると、すべてのレベルの解析が終了となり、どのレベルでの解析結果を表示するか聞いてくる。

0 < j < 9

j =

レベル５の結果を表示するときは、次のように５と設定すればよい。

0 < j < 9

j = 5

レベルとして５を設定すると、レベル５での解析結果が表示される（図６）。

図C6

レベルとして範囲外の値を設定すると、プログラムの実行終了である。次の例では、０を設定している。

0 < j < 9

j = 0

PS ******\wa_color_files>

リストC1のスクリプトは、サンプル画像ファイルとともに圧縮ファイルwa_color_files.zipにまとめた。ダウンロード展開すれば自由に利用できる。

リストC1　カラー画像対応のウェブレット解析（wa_color.py）

"""

Wavelet Analysis of Color Image

Yasuharu Okamoto, 2023.11

"""

import cv2 # pip install opencv-python

import numpy as np

import matplotlib.pyplot as plt

img1 = cv2.imread(input('File name = ')) # 画像ファイルの読み込み

img1 = cv2.cvtColor(img1, cv2.COLOR_BGR2RGB) # RGBに変換

plt.imshow(img1)

plt.title('Input image')

plt.show()

img1 = cv2.resize(img1, dsize = (512, 512)) # 512×512に調整

plt.imshow(img1)

plt.title('Resized image')

plt.show()

img1 = np.array(img1, dtype='float') / 255 # int8からfloatへの変更

#

#　　　ウェーブレット解析結果を辞書型のリストWaveletAnalに保存

#

d1 = {'0': {'f0': img1[:,:,0]}} # R-channel

d2 = {'0': {'f0': img1[:,:,1]}} # G-channel

d3 = {'0': {'f0': img1[:,:,2]}} # B-channel

WaveletAnal = [d1, d2, d3]

n = 9

N = 2**n

#

# レベルｊの解析結果の図示

#

def disp_map(j):

Cj = [0]*3

DHj = [0]*3

DVj = [0]*3

DDj = [0]*3

for i in range(3):

Cj[i] = WaveletAnal[i][f'{j}']['f0']

DHj[i] = WaveletAnal[i][f'{j}']['DHj']

DVj[i] = WaveletAnal[i][f'{j}']['DVj']

DDj[i] = WaveletAnal[i][f'{j}']['DDj']

#

# 画像を配列imgsに設定

# 大きさ512×512を４分割して使用

#

imgs = np.full((512, 512, 3), 0.5)

for h in range((2**(n-j))):

for k in range((2**(n-j))):

for ii in range((2**j)//2): # スケールに合わせて引き延ばし

for jj in range((2**j)//2):

for kk in range(3): # 赤、緑、青の３チャンネル

#

# C_hkの値

#

imgs[512//2 + ((2**j)//2)*h + ii]\

[512//2 + ((2**j)//2)*k + jj][kk] \

= Cj[kk][h][k] #/ 255

#

# DH, DV, DDは、値０を0.5にシフト

#

imgs[512//2 + ((2**j)//2)*h + ii][((2**j)//2)*k + jj][kk] \

= DHj[kk][h][k] + 0.5

imgs[((2**j)//2)*h + ii][(512//2) + \

((2**j)//2)*k + jj][kk] \

= DVj[kk][h][k] + 0.5

imgs[((2**j)//2)*h + ii][((2**j)//2)*k + jj][kk] \

= DDj[kk][h][k] + 0.5

plt.imshow(imgs) # 画像配列imgsの表示

plt.xticks([])

plt.yticks([])

plt.title(f'j = {j}')

plt.show()

print('\nWavelet analysis started.')

for j in range(1,n):

# レベルｊの解析

print(f'Level {j} started.')

Cj_1 = [0]*3

for i in range(3):

Cj_1[i] = WaveletAnal[i][f'{j-1}']['f0']

#print('Cj_1 =\n', Cj_1)

Cj = np.empty((3, 2**(n-j), 2**(n-j)))

DHj = np.empty((3, 2**(n-j), 2**(n-j)))

DVj = np.empty((3, 2**(n-j), 2**(n-j)))

DDj = np.empty((3, 2**(n-j), 2**(n-j)))

#

#　　　　下降サンプリング

#

for h in range(2**(n-j)):

for k in range(2**(n-j)):

for i in range(3):

#print('h =', h, ' k =', k, ' i =', i)

Cj[i][h][k] = (Cj_1[i][2*h][2*k] + Cj_1[i][2*h+1][2*k] +

Cj_1[i][2*h][2*k+1] + Cj_1[i][2*h+1][2*k+1]) / 4

DHj[i][h][k] = (Cj_1[i][2*h][2*k] - Cj_1[i][2*h+1][2*k] +

Cj_1[i][2*h][2*k+1] - Cj_1[i][2*h+1][2*k+1]) / 4

DVj[i][h][k] = (Cj_1[i][2*h][2*k] + Cj_1[i][2*h+1][2*k] -

Cj_1[i][2*h][2*k+1] - Cj_1[i][2*h+1][2*k+1]) / 4

DDj[i][h][k] = (Cj_1[i][2*h][2*k] - Cj_1[i][2*h+1][2*k] -

Cj_1[i][2*h][2*k+1] + Cj_1[i][2*h+1][2*k+1]) / 4

for i in range(3):

WaveletAnal[i][f'{j}'] = {'f0': Cj[i]}

WaveletAnal[i][f'{j}']['DHj'] = DHj[i]

WaveletAnal[i][f'{j}']['DVj'] = DVj[i]

WaveletAnal[i][f'{j}']['DDj'] = DDj[i]

disp_map(j) # 解析結果の表示

while True:

print(f'0 < j < {n}')

j = int(input('j = '))

if j <=0 or j >= n:

break

disp_map(j) # レベルｊの表示

Home