特異値分解による画像の低ランク近似

はじめに

特異値分解、してますか？＞挨拶

なんか最近まわりの人が特異値分解ばかりしている気がするので、僕もやってみることにしました。とりあえず画像をそのまま行列だと思って特異値分解(Singular-Value-Decomposition, SVD)して、低ランク近似した画像を作ります。特異値分解についてはWikipediaとかを参照してください。

画像を特異値分解する話を誰がはじめたのかは知りませんが、僕はこの論文T. Kumamoto, M. Suzuki, and H. Matsueda arXiv:1608.08333で知りました。

ソースコードは https://github.com/kaityo256/image_svd においてあります。

画像の読み込み

画像の読み込みにはPythonのPillowを使います。もし入ってなければ

$ pip install Pillow

とかで入れましょう。

入力画像ですが、文字があった方が近似精度がわかりやすいので、とりあえずこんなものを用意してみました。

これをまずグレースケールに変換します。

from PIL import Image
img = Image.open('stop.jpg')
gray_img = img.convert('L')
gray_img.save('stop_mono.jpg')

こんな感じになります。

特異値分解

特異値分解のために、numpyとscipyを使います。importしておきましょう。

import numpy as np
from scipy import linalg

イメージデータへのアクセスは、そのままnumpyのasarrayで読み込めます。後はscipyのlinalg.svdでSVDできます。楽ちん過ぎる。

a = np.asarray(gray_img)
u, s, v = linalg.svd(a)

ここでできたu, s, vが特異値分解の結果です。このうちsが特異値のリストになっています。

低ランク近似

特異値分解されたもののうち、上位のランクだけ選んで再構成すれば低ランク近似になります。

ただし、Σ(コードではs)は本当は行列ですが、linalg.svdの返り値ではベクトルになっているため、行列になおしてやる必要があります。後は必要な分だけスライスしたものを掛け算すれば低ランク近似した行列の完成です。コードにするとこんな感じ。

rank = 10
ur = u[:, :rank]
sr = np.matrix(linalg.diagsvd(s[:rank], rank,rank))
vr = v[:rank, :]
b = np.asarray(ur*sr*vr)
img2 = Image.fromarray(np.uint8(b))
img2.save('stop_r10_mono.jpg')

得られた結果はこんな感じ。

カラー版

上記は簡単のためにグレースケールにしたけれど、RGBの各要素に対してSVDすればカラー版もできる。ここで一つだけハマりポイントなのは、モノクロ版ではImage.fromarrayでそのまま突っ込めたけど、カラーでは(r,g,b)のタプルの配列にしてやらないといけないこと。numpyの行列の掛け算でできるオブジェクトはMatrixオブジェクトなので、その後配列としてアクセスするためにはnp.asarrayで配列に戻してやらないといけない。そこだけ気をつければすぐにできる。

行列の低ランク近似の関数をこんな感じで作る。

def perform_svd(a,rank):
    u, s, v = linalg.svd(a)
    ur = u[:, :rank]
    sr = np.matrix(linalg.diagsvd(s[:rank], rank,rank))
    vr = v[:rank, :]
    return np.asarray(ur*sr*vr)

返り値をasarrayでarrayにしているのが地味にポイント。これを使うと、カラー版も簡単にできる。

    w = img.width
    h = img.height
    A = np.asarray(img)
    r = perform_svd(A[:,:,0],rank).reshape(w*h)
    g = perform_svd(A[:,:,1],rank).reshape(w*h)
    b = perform_svd(A[:,:,2],rank).reshape(w*h)
    B = np.asarray([r,g,b]).transpose(1,0).reshape(h,w,3)
    img2 = Image.fromarray(np.uint8(B))
    img2 = Image.fromarray(np.uint8(data))
    img2.save(`stop_r10.jpg`)

イメージデータからnumpy配列を作ると、(高さ、幅、色)の3次元配列になるので、色ごとにSVDして、それをくっつけると(色, 高さ*幅)の配列になっているので、それをtransposeで入れ替えた後に元の3次元配列の形にreshapeして、Image.fromarrayに突っ込めばよい。

できた画像はこんな感じ。

おわりに

この画像だとランクを20個も取ればほぼもとの画像を回復できる。どのくらい取ればよいかとか、どういう画像がどういうスペクトルを持っているかとか調べると楽しいんじゃないでしょうか。

いやしかし、画像読み込んでグレースケールにして行列にしてSVDして低ランク近似して画像にもどして保存って、結構いろいろやってるのにPythonで適切なライブラリをimportすれば15行とかでできちゃうのね・・・。つくづく「ライブラリの充実度が言語の強さ」だと思う・・・