画像の特異値分解（Singular Value
Decomposition: SVD）
画像の特異値分解（Deisenroth et al. 2020）を試みた。特異値分解については、岡本（2009）、pp.47-48、岡本（2019）、pp.67-68で説明している。
ｍピクセル×ｎピクセルのモノクロ画像はｍ×ｎの行列Ａで表すことができる。行列Ａは次のように特異値分解できる。
ここで、行列Ａ、すなわち画像を、少数の特異値分解の項（階数１の行列）の和で近似することを考える。
カラー画像の場合は、RGBそれぞれに対応して行列を考えるが、このサイトの後半で扱う。
画像の特異値は、急激に減少する。
例えば、図１の画像の場合、
図１
特異値は図２のように減少する。
図２
図３
特異値の分布は図４のようになる。
図４
モノクロ画像の特異値分解を行い、最初の少数の特異値分解の項による近似を調べるスクリプトを、リスト１のように用意した。スクリプトファイルは、圧縮ファイルsvd_figure_mono.zipにまとめている。ダウンロード展開すれば自由に利用できる。リスト１のスクリプトでは、OpenCVを用いているが、OpenCVのインストールは、このウェブサイトで説明しているように簡単である。
リスト１のスクリプトを実行すると、入力画像ファイル名の設定が求められる（図５）。
図５
次のスクリプトで読み込まれた画像ファイルは、400×400のサイズに変換され、モノクロ化される。
f_img = input('Input figure file
(*.png) = ')
img1 = cv2.imread(f_img)
# 画像の読み込み
img1 = cv2.resize(img1, dsize=(400,400))
# 400×400へのサイズ変換
#
# 画像のモノクロ化
#
img0 = np.empty((400, 400))
for i in range(400):
for j in range(400):
v =
0
for
k in range(3):
v += img1[i][j][k]
img0[i][j] = v / (255*3)
適当な画像ファイル名を入力する。サイズは読み込み後、調整されるので任意である。
Input figure
file (*.png) = sample_4.png
画像ファイル名が読み込まれると、サイズを調整してモノクロ化された画像が表示される（図６）。
図６
図６のフォームを閉じると、次のコードで特異値分解が行われる。
u, s, vh = np.linalg.svd(img0,
full_matrices=False)
特異値分解が終了すると、次のスクリプトで特異値のグラフ（図７）が描かれる。
plt.plot(range(len(s)), s)
plt.xticks([0, len(s)-1], ['1',
f'{len(s)}'])
plt.ylabel('Singular value')
plt.show()
図７
初めの方の少数のものを除いて、特異値は０に近いことが判る。
図７のフォームを閉じると、最初の20個の特異値のグラフが描かれる（図８）。
図８
図８のフォームを閉じると、特異値分解に基づく画像の構成が行われる。
１番目から４番目までの特異値に対する画像がそれぞれの特異値に対して構成される。
# h番目の特異値による分解
f1 =
u[:,h:h+1] @ np.diag(s[h:h+1]) @ vh[h:h+1]
また、１番目から４番目までの特異値に対する画像をまとめたものも構成される。
# 1番目から4番目までの特異値による分解の和
f1 =
u[:,:4] @ np.diag(s[:4]) @ vh[:4]
構成された画像は、白黒画像を表すRGB画像に変換されて、リストimgsに加えられる。
img_1 = np.empty((400, 400,
3), dtype=float)
min_f1 = np.min(f1)
if min_f1 < 0.0:
f1 -=
min_f1
max_f1 = np.max(f1)
if max_f1 > 1.0:
f1
/= max_f1
for i in range(len(f1)):
for
j in range(len(f1[0])):
for k in range(3):
img_1[i][j][k] = f1[i][j]
imgs.append(img_1)
リストimgs内の画像は、次のコードで作成された大画像img_gに貼りつけられる。
img_g =
np.full((800, 1200, 3), 0.5, dtype=float)
この大画像は、次のスクリプトで表示される（図９）。
plt.figure(figsize=(12,8))
plt.imshow(img_g)
plt.show()
図９
１行目1列目に、特異値分解の対象となった400ピクセル×400ピクセルの元画像が表示されている。その右へ順に、第１特異値から第４特異値までの各特異値に対応する画像が並べられ、最後の右下の画像は第1特異値から第４特異値までの分解を合わせた画像である。右下の画像は、左上の画像をかなり良く表している。左上の画像は、特異値分解の対象となったもので
次元
の情報を表すものであり、右下の画像は
次元
の情報を表すものである。したがって、次元の削減率は
3204次元／160000次元≒2.00％
となる。
ちなみに、図３の画像の場合は、特異値の分布は、図10となり、
図10
画像の分解は、図11となる。
図11
リスト１ モノクロ画像の特異値分解（svd_figure_mono.py）
import cv2
import numpy as np
import matplotlib.pyplot as plt
def mono2rgb(img):
# モノカラーからRGBへ
img3 = np.empty((*img.shape,
3))
for i in range(len(img)):
for j
in range(len(img[0])):
for k in range(3):
img3[i][j][k] = img[i][j]
return img3
f_img = input('Input figure file
(*.png) = ')
img1 = cv2.imread(f_img)
# 画像の読み込み
img1 = cv2.resize(img1, dsize=(400,400))
# 400×400へのサイズ変換
#
# 画像のモノクロ化
#
img0 = np.empty((400, 400))
for i in range(400):
for j in range(400):
v =
0
for
k in range(3):
v += img1[i][j][k]
img0[i][j] = v / (255*3)
plt.imshow(mono2rgb(img0)) # 変換画像の表示
plt.title('Transformed Input Figure')
plt.show()
#
# 画像の特異値分解
#
u, s, vh = np.linalg.svd(img0,
full_matrices=False)
#
# 特異値のプロット
#
plt.plot(range(len(s)), s)
plt.xticks([0, len(s)-1], ['1',
f'{len(s)}'])
plt.ylabel('Singular value')
plt.show()
plt.plot(range(20), s[:20])
plt.xticks([0, 19], ['1', '20'])
plt.ylabel('Singular value')
plt.show()
imgs = []
for h in range(5):
print('h =', h)
if h < 4:
# h番目の特異値による分解
f1 =
u[:,h:h+1] @ np.diag(s[h:h+1]) @ vh[h:h+1]
else:
# 1番目から4番目までの特異値による分解の和
f1 =
u[:,:4] @ np.diag(s[:4]) @ vh[:4]
img_1 = np.empty((400, 400,
3), dtype=float)
min_f1 = np.min(f1)
if min_f1 < 0.0:
f1
-= min_f1
max_f1 = np.max(f1)
if max_f1 > 1.0:
f1
/= max_f1
for i in range(len(f1)):
for
j in range(len(f1[0])):
for k in range(3):
img_1[i][j][k] = f1[i][j]
imgs.append(img_1)
#
#
描画用大画面
#
img_g = np.full((800, 1200, 3), 0.5,
dtype=float)
def insert_img(startx0, starty0,
img_v): # 大画面への画像の挿入
for i in range(400):
for
j in range(400):
for k in range(3):
img_g[startx0+i][starty0+j][k] = img_v[i][j][k]
insert_img(0, 0, mono2rgb(img0))
for h in range(5):
if h < 2: # １行目に表示
insert_img(0, 400*(h+1), imgs[h])
else:
# ２行目に表示
insert_img(400, 400*(h-2), imgs[h])
plt.figure(figsize=(12,8))
plt.imshow(img_g)
plt.show()
ｍピクセル×ｎピクセルの画像は、各RGBのチャネルに対してｍ×ｎ型の行列
、
、
で表すことができる。各行列Ａは次のように特異値分解できる。
ここで、行列Ａ、すなわち画像を、少数の特異値分解の項（階数１の行列）の和で近似することを考える。なお、特異値分解については、岡本（2009）、pp.47-48、岡本（2019）、pp.67-68で説明している。
例えば、図C.1の場合、
図C.1
特異値分解したときの特異値の分布は図C.2のようになる。
図C.2
少数の特異値分解の項の和で、よく近似できることがわかる。
カラー画像の特異値分解のスクリプト例をリストC.1のように用意した。スクリプトファイルなどは、圧縮ファイルsvd_figure_color.zipにまとめたので、ダウンロード展開すれば自由に利用できる。パッケージOpenCVを用いているが、OpenCVのインストールはこのウェブサイトで説明されているように簡単である。
リストC.1のスクリプトを実行すると、入力画像ファイル名の設定が求められる（図C.3）。
図C.3
特異値分解する画像ファイル名を設定する。
Input file
name (*.png) = sample_pctr.png
ファイル名を設定すると、以下のスクリプトにより読み込まれて、400ピクセル×400ピクセルの画像に変換される。
f_img_nm = input('Input file name
(*.png) = ')
img1 = cv2.imread(f_img_nm)
# 画像ファイルの読み込み
img1 = cv2.resize(img1,
dsize=(400,400)) # サイズを400×400に調整
img1 = cv2.cvtColor(img1,
cv2.COLOR_BGR2RGB) / 255
# 値の範囲を[0, 1]に変換
plt.imshow(img1)
plt.show()
変換された画像は、図C.4のように表示される。
図C.4
図C.4のフォームを閉じると、読み込まれ変換された画像img1は、次のスクリプトでRGBに対応する３つの画像に分けられる。
img0 = np.empty((3, 400, 400))
for i in range(400):
for j in range(400):
for
k in range(3):
img0[k][i][j] = img1[i][j][k]
RGBの各チャンネルに分けられた画像は、図C.5のように表示される。
図C.5
図C.5のフォームを閉じると、次のコードで特異値分解が行われる。
u, s, vh =
np.linalg.svd(img0, full_matrices=False)
関数svdの説明は、ウェブサイト
https://numpy.org/doc/stable/reference/generated/numpy.linalg.svd.html
にある。
入力配列img0において、RGBの画像配列がimg0[0]、img0[1]、img0[2]に設定されているので、、特異値分解の結果はu[0]、s[0]、vh[0]にRチャンネルの特異値分解、u[1]、s[1]、vh[1]にGチャンネルの特異値分解、u[2]、s[2]、vh[2]にBチャンネルの特異値分解が設定される（表される）。
特異値の分布は、RGBチャンネルごとに図C.6のように表される。
図C.6
最初の少数の特異値分解の項で、画像が表されるようである。
図C.6のフォームを閉じると、初めの20個の特異値のグラフが表示される（図C.7）。
図C.7
図C.7のフォームを閉じると、しばらくの計算時間の後、特異値分解の各項による画像が表示される（図C.8）。
図C.8
元画像が1行目１列目に表示され、以下、各特異値分解の項の表す画像が表示される。最後（右下）に、６個の特異値分解の項をまとめた画像が表示されている。左上の元画像が、右下の６個の固有分解の項による画像でかなり良く表されている。
リストC.1 カラー画像の特異値分解（svd_figure_color.py）
import cv2
# pip install opencv-python
import numpy as np
import matplotlib.pyplot as plt
f_img_nm = input('Input file name
(*.png) = ')
img1 = cv2.imread(f_img_nm)
# 画像ファイルの読み込み
img1 = cv2.resize(img1,
dsize=(400,400)) # サイズを400×400に調整
img1 = cv2.cvtColor(img1,
cv2.COLOR_BGR2RGB) / 255
# 値の範囲を[0, 1]に変換
plt.imshow(img1)
plt.show()
"""
モノクロ画像を（R,G,B）の３枚用意
"""
img0 = np.empty((3, 400, 400))
for i in range(400):
for j in range(400):
for
k in range(3):
img0[k][i][j] = img1[i][j][k]
"""
モノクロ画像imgを、sで指定されたRGB画像として返す
"""
def select_rgb(img, s):
t_img =
np.zeros((*np.shape(img[0]),3), dtype=float)
for i in range(len(img[0])):
for
j in range(len(img[0][0])):
t_img[i][j][s] = img[s][i][j]
return t_img
"""
R、G、Bの３枚の画像、img[0],img0[1],img0[2]、の表示
"""
plt.figure(figsize=(9,3))
plt.subplot(131)
plt.title('R')
plt.imshow(select_rgb(img0, 0))
plt.subplot(132)
plt.title('G')
plt.imshow(select_rgb(img0, 1))
plt.subplot(133)
plt.title('B')
plt.imshow(select_rgb(img0, 2))
plt.show()
"""
img0に含まれる３枚の画像の特異値分解
"""
u, s, vh = np.linalg.svd(img0,
full_matrices=False)
rgb_labels = ['R', 'G', 'B']
plot_colors = ['r', 'g', 'b']
for t in range(3):
plt.plot(range(len(s[t])),
s[t], c = plot_colors[t], label = rgb_labels[t])
plt.legend(fontsize=20)
plt.show()
for t in range(3):
plt.plot(range(20),
s[t,:20], c = plot_colors[t], label = rgb_labels[t])
plt.xticks([0, 19], ['1', '20'])
plt.legend(fontsize=20)
plt.show()
"""
１番目から６番目までの特異値ごとの画像、
およびそれらの合成画像の作成
"""
imgs = []
for h in range(7):
print('h =', h)
img_1 = np.empty((400, 400,
3), dtype=float)
for t in range(3):
if h
< 6:
f1 = u[t, :,h:h+1] @ np.diag(s[t, h:h+1]) @ vh[t, h:h+1]
else:
f1 = u[t, :,:6] @ np.diag(s[t, :6]) @ vh[t, :6]
"""
画素の値を０から１までの範囲に調整
"""
min_f1 = np.min(f1)
if
min_f1 < 0.0:
f1 -= min_f1
max_f1 = np.max(f1)
if
max_f1 > 1.0:
f1 /= max_f1
for
i in range(len(f1)):
for j in range(len(f1[0])):
img_1[i][j][t] = f1[i][j]
imgs.append(img_1)
def v_s(h):
# h番目の特異値の和
return np.sum(s, axis=0)[h]
"""
特異値分解の表示
"""
fig = plt.figure(figsize=(12,8))
plt.subplot(241)
plt.title('Original Figure')
plt.imshow(img1)
# 元画像の表示
"""
各特異値分解および
それらの合成画像の表示
"""
for h in range(7):
if h < 3:
plt.subplot(240+h+2)
plt.title(f's{h+1} = {v_s(h):.3f}')
plt.imshow(imgs[h])
else:
plt.subplot(240+h+2)
plt.imshow(imgs[h])
if h
< 6:
plt.title(f's{h+1} = {v_s(h):.3f}')
else:
plt.title('s1+...+s6')
fig.suptitle('Singular Value
Decomposition of a Figure', fontsize=20)
fig.tight_layout()
fig.show()
Deisenroth, M. P., Faisal, A. A., & Ong, C. S. (2020) Mathematics for machine learning. Cambridge University Press.
岡本安晴（2009）統計学を学ぶための数学入門［下］．培風館
岡本安晴（2019）いまさら聞けないPythonでデータ分析．丸善出版