Pythonによるデータ処理11 ～ウェーブレット変換4

　今回は前回に引き続きPythonのPywaveletsモジュールを逆2Dウェーブレット変換と圧縮について紹介する。

atatat.hatenablog.com

1. 準備
2. PyWaveletsによる逆2Dウェーブレット変換
3. ウェーブレット変換を利用した圧縮
4. スパース行列への変換
5. まとめ

1. 準備

前回と同様に標準画像データベースのうち、Boatを使用する。Google Colabの場合は以下のようにファイルを読み込み、モノクロ画像grayとして処理を行う。

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
%matplotlib inline
import pywt

from google.colab import files
uploaded_file = files.upload()

uploaded_file_name = next(iter(uploaded_file))
orig = cv2.imread(uploaded_file_name)

gray = cv2.cvtColor(orig, cv2.COLOR_BGR2GRAY)
plt.imshow(gray,cmap='gray',vmax=255,vmin=0)
plt.title('Original')

2. PyWaveletsによる逆2Dウェーブレット変換

前回と同様に読み込んだ画像を2Dウェーブレット変換すると、Approximation: LL , とdetails:LH,HL,HHが得られる。元のモノクロ標準画像は256 x 256のUnsigned Int 8 (0-255)で、LL,HL,LH,HHは130x130のfloat64になっている。2Dウェーブレット変換: dwt2(gray, 'bior1.3')の出力形式は1x2のアレイで1要素目はLL, 2要素目は(LH,HL,HH)のタプルで格納されている。

titles = ['Approximation', ' Horizontal detail',
          'Vertical detail', 'Diagonal detail']
coeffs = pywt.dwt2(gray, 'bior1.3')
LL, (LH, HL, HH) = coeffs
fig = plt.figure(figsize=(12, 3))
for i, a in enumerate([LL, LH, HL, HH]):
    ax = fig.add_subplot(1, 4, i + 1)
    ax.imshow(a, interpolation="nearest", cmap='gray')
    ax.set_title(titles[i], fontsize=10)
    ax.set_xticks([])
    ax.set_yticks([])

fig.tight_layout()

"2Dウェーブレット変換(bior1.3)" — 2Dウェーブレット変換(bior1.3)

dwt2関数が出力するcoeffと変換に用いた基底関数を指定して2D逆ウェーブレット変換を行うためにidwt2(coeffs,'bior1.3')とすると元の画像grayを復元することができる。

recov=pywt.idwt2(coeffs,'bior1.3')
plt.figure(figsize=(8,6))
plt.subplot(1,2,1)
plt.imshow(gray,cmap='gray',vmax=255,vmin=0)
plt.title('original')
plt.subplot(1,2,2)
plt.imshow(recov,cmap='gray',vmax=255,vmin=0)
plt.title('recovered')

以下元の画像が逆2Dウェーブレット変換により復元できていることが確認できる。

"逆2Dウェーブレット変換による画像復元" — 逆ウェーブレット変換による画像復元

3. ウェーブレット変換を利用した圧縮

　ここではウェーブレット変換を利用すると圧縮が可能ということを示すが、jpegなどに一般的に採用されているアルゴリズムではないことに注意。ところでApproximationは130x130, 高次情報(LH,HL,HH)は合わせて130 x 130 x 3のためウェーブレット変換を復元するためには合計で67600データが必要になる。一方で元の画像は256 x 256 = 65536 データのため圧縮にならない。そこで高次情報は適当な閾値のバイナリマスクを適用して、0以外を保持するスパース行列とすることを考える。

まず高次情報LH,HL,HHに対して閾値20とするマスクを適用する。

MLH=(np.abs(LH)>20)
MHL=(np.abs(HL)>20)
MHH=(np.abs(HH)>20)

LH2=LH*MLH
HL2=HL*MHL
HH2=HH*MHH
plt.figure(figsize=(15,6))
plt.subplot(1,4,1)
plt.imshow(LL,cmap='gray')
plt.title('Approximation')

plt.subplot(1,4,2)
plt.imshow(LH2,cmap='gray')
plt.title('Horizontal')

plt.subplot(1,4,3)
plt.imshow(HL2,cmap='gray')
plt.title('Vertial')

plt.subplot(1,4,4)
plt.imshow(HH2,cmap='gray')
plt.title('Diagonal')

上記マスクを適用した高次情報を1要素目はLL, 2要素目は(LH,HL,HH)のタプルで格納した1x2のアレイ：[LL,(LH2,HL2,HH2)]を係数として以下のようにidwt2()で画像を復元する。

recov2=pywt.idwt2([LL,(LH2,HL2,HH2)],'bior1.3')
plt.figure(figsize=(12,6))
plt.subplot(1,3,1)
plt.imshow(gray,cmap='gray',vmax=255,vmin=0)
plt.title('original')

plt.subplot(1,3,2)
plt.imshow(LL,cmap='gray')
plt.title('Approximation')

plt.subplot(1,3,3)
plt.imshow(recov2,cmap='gray')
plt.title('recovered')

基本的にはマスクを適用しても、復元画像の質が大きく落ちていないことが確認できる。

4. スパース行列への変換

　マスク適用した高次情報をスパース行列に変換することで非０の要素のみを保存する。スパース行列はいくつかあるが、今回はcsr_matrixをscipy.sparseからインポートする。マスクを適用した高次情報:LH2, HL2, HH2をスパース行列に変換し、要素数を算出すると19907で元の画像の256 x 256 = 65536 に対して1/3程度になっている。今回は要素数が減ることのみを示すが、メモリサイズは元画像がUnsigned int 8 に対してウェーブレット変換後はfloat64になるため、メモリを含めて圧縮するには型変換などを行う必要がある。

from scipy.sparse import csr_matrix

cLH = csr_matrix(LH2)
cHL = csr_matrix(HL2)
cHH = csr_matrix(HH2)
csize=np.size(cLH)+np.size(cHL)+np.size(cHH)+np.size(LL)
print('Compressed size: ',csize)
print('Original size: ',np.size(gray))

>
Compressed size:  19907
Original size:  65536

スパース化した行列をもとに戻すにはスパース行列に.toarray()を付加する。非ゼロの要素のみを保持しているスパース行列から通常の行列に戻して、逆2Dウェーブレット変換することで画像を復元することができる。

rLH = cLH.toarray()
rHL = cHL.toarray()
rHH = cHH.toarray()

recov3=pywt.idwt2([LL,(rLH,rHL,rHH)],'bior1.3')
plt.figure(figsize=(12,6))
plt.subplot(1,3,1)
plt.imshow(gray,cmap='gray',vmax=255,vmin=0)
plt.title('original')

plt.subplot(1,3,2)
plt.imshow(LL,cmap='gray')
plt.title('Approximation')

plt.subplot(1,3,3)
plt.imshow(recov3,cmap='gray',vmin=0,vmax=255)
plt.title('recovered')