甲斐性のない男が機械学習とか最適化とかを綴るブログ

うどんくらいしか食べる気がしない

近接勾配法応用編その1 ~スパースコーディング、辞書学習からの超解像~

はじめに

前回の記事で近接勾配法の基本と実装を見てきました。今回はこの近接勾配法を利用して、スパースコーディングの応用の1つである超解像技術を調べPythonで実装してみました*1。画像処理関係の実装はほぼ未経験だったので、その分野の人が見れば前処理等が足りなかったり、アンチパターン的なことをやってたりするかもしれませんが、その点は指摘頂けると幸いです。

スパースコーディングとは

スパースコーディングは観測した信号{\bf y}を、「辞書」と呼ばれるm個の基底ベクトルの集合  {\bf D} = \left[{\bf d}_1, {\bf d}_2, \cdots, {\bf d}_m \right]のスパースな線形結合で表現する符号化手法です。
今、この線形結合の係数を{\bf w}とすると、スパースコーディングは下記のように定式化されます。
\displaystyle {\bf y }\simeq {\bf D} {\bf w} \tag{1} 
ただし、{\bf w}の要素の多くは0(つまり{\bf w}はスパース)です。

このようなスパースコーディングを使うことによりデータの圧縮を行えるだけでなく、画像からのノイズ除去や今回取り上げる超解像など様々な技術へ応用が可能です。
スパースコーディングを成功させる上で重要となるのは辞書{\bf D}です。{\bf D}がテキトーだとうまくいきません*2。ということで、辞書{\bf D}は実データからの学習により構築します。ここで、学習データの個数をn、学習データをまとめた行列{\bf Y} = \left[{\bf y}_1 \cdots, {\bf y}_n\right] 、 各線形結合の係数をまとめた行列{\bf W} = \left[{\bf w}_1, \cdots, {\bf w}_n\right]を定義すると、辞書学習の目的関数は下記のようになります。
\displaystyle {\rm arg} \min_{{\bf D}, {\bf W}} \| {\bf Y} - {\bf D}{\bf W} \|^2_2 \tag{2} 
ただし、{\bf W}はスパース。

この方法で辞書{\bf D}が得られたら、次はその辞書を用いて入力画像の符号化します。そして、その係数と辞書から画像を再構成し超解像画像を得るという流れになります。
ということで、スパースコーディングによる超解像生成は(1)辞書を学習するフェーズと(2)その辞書を元に画像を再構成するフェーズの2フェーズで実現できます。

画像再構成のための辞書学習

実装した手法の説明

では、今回実装した学習用の画像データを元に辞書を学習する方法を説明します。
まず、最終的に得たい解像度の画像データ(高解像画像データ)複数枚を学習データとして用意します。この時、超解像を行いたい画像のカテゴリが決まっているなら、そのカテゴリに合った画像を用意するのがよさそうです。
この学習データを用いて、辞書{\bf D}を学習する流れを下記に示します。

  1. 学習用画像を小領域を取り出した「画像パッチ」をn個用意し、それぞれ縦に並べたベクトル{\bf y}_1, {\bf y}_2, \cdots, {\bf y}_nおよび、これらのベクトルを横に並べた行列を{\bf Y} = \left[{\bf y}_1 \cdots, {\bf y}_n\right] を作る。
  2. {\bf Y}の各行が平均0、分散1になるように標準化する。
  3. {\hat {\bf D}}をランダム行列で初期化する。
  4. {\hat {\bf D}}を固定して、次の最適化問題i = 1, 2, \cdots, nについて解き、行列{\hat {\bf W}} = \left[{\hat {\bf w}_1}, \cdots, {\hat {\bf w}_n} \right]を得る。
  5. \displaystyle {{\hat {\bf w}_i}} = {\rm arg} \min_{{\bf w}_i} \| {\bf y}_i - {\hat {\bf D}}{\bf w}_i \|^2_2  + \lambda \| {\bf w}_i \|_1\tag{3}
  6. \hat{{\bf W}}を固定して、次の最適化問題を解き\hat{{\bf D}}を得る。
  7. \displaystyle \hat{{\bf D}}=  \arg \min_{{\bf D}} \| {\bf Y} - {\bf D} {\hat{ \bf W}} \|^2_2 \tag{4}
  8. 手順4,手順5の操作を収束するまで、もしくは一定回数繰り返す。

ここで、式(3)はlassoの式であるため解がスパースになることが期待され、更に近接勾配法で解くことができます。なお、今回、近接勾配法の応用編ということで手順4の方法でスパースな係数\bf{w}_iを求めていますが、この方法以外にもOrthogonal Matching PursuitやIterative Reweighted Least Squaresなどのアルゴリズムがあります。
また、手順5で式(4)の解を得る方法としてはK-SVDという手法がよく使われますが、今回は簡単のために下記のように解析的に解いた結果を使用します。
\displaystyle {\hat {\bf D}} = {\bf Y} {\hat {\bf W}}^T ( {\hat {\bf W}} {\hat {\bf W}}^T  )^{-1} \tag{5}

この手順4、手順5を繰り返して{\bf W}, {\bf D}の同時最適化を行うわけですが、この両変数に対しては凸関数ではないため(式(3)、式(4)はそれぞれ凸関数ですが)、求まる結果は局所最適解です。
なお、この学習を収束まで待っているとかなり時間がかかるため、今回は適当な回数ループさせたら学習を完了させています。

ちなみに、最初は手順2の標準化を省いた実装をしていましたが、これだとうまく超解像ができませんでした。やはり、学習データの標準化は重要のようです。

実装したコード

ここまでの手順を実現するpythonコードです。

import numpy as np
import scipy as sci
from numpy.random import *
from PIL import Image
import os
import copy
import sys

#dorectory以下のファイルを再帰的に検索する関数
def find_all_files(directory):
    for root, dirs, files in os.walk(directory):
        yield root
        for file in files:
            yield os.path.join(root, file)

#おれおれ画像・ベクトル・行列変換関数群
def color2gray_img(img):
    return img.convert('L')

def im2vec(img):
    img_arr = np.array(img.convert('L'),'float')
    return img_arr.reshape(-1)

def gray_im2float_vec(gray_img):
    return im_scale2float((im2vec(gray_img)))

def color_im2float_vec(color_img):
    return gray_im2float_vec(color2gray_img(color_img))

def vec2image(vec, H, W):
    tmp_vec = copy.deepcopy(vec)
    tmp_vec = tmp_vec.reshape((H, W))
    return Image.fromarray(np.uint8(tmp_vec))

def float_vec2image(vec, H, W):
    return vec2image(float2im_scale(vec), H, W)

def arr2image(arr):
    arr = copy.deepcopy(arr)
    return Image.fromarray(np.uint8(arr))

def float_arr2image(arr):
    return arr2image(float2im_scale(arr))

def im_scale2float(vec):
    return vec.astype(np.float)

def float2im_scale(vec):
    vec[vec > 255.0] = 255.0
    vec[vec < 0.0] = 0.0
    return vec

#画像データからパッチを切り出す関数
#shift_numで切り出す際のずらし量を指定
def gray_im2patch_vec(im, patch_hight, patch_width, shift_num):
    patchs = np.zeros((patch_hight*patch_width,
                        (1 + (im.size[0] - patch_hight)/shift_num + (1 if (im.size[0] - patch_hight) % shift_num  != 0 else 0)) *
                        (1 + (im.size[1] - patch_width)/shift_num + (1 if (im.size[1] - patch_width) % shift_num  != 0 else 0))))
    im_arr = im_scale2float(np.array(im))
    k = 0
    for i in range(1 + (im.size[0] - patch_hight)/shift_num + (1 if (im.size[0] - patch_hight) % shift_num  != 0 else 0)):
        h_over_shift = min(0, im.size[0] - ((i* shift_num) + patch_hight)) #ずらしたことによって画像からはみ出る量(縦方向)
        for j in range(1 + (im.size[1] - patch_width)/shift_num + (1 if (im.size[1] - patch_width) % shift_num  != 0 else 0)):
            w_over_shift = min(0, im.size[1] - ((j* shift_num) + patch_width)) #ずらしたことによって画像からはみ出る量(横方向)
            #パッチの切り出し 画像からはみ出たら、そのはみ出た分を戻して切り出し
            patchs[:, k] = im_arr[(i* shift_num) + h_over_shift:((i* shift_num) + patch_hight) + h_over_shift,
                                    (j* shift_num) + w_over_shift:((j* shift_num) + patch_width) + w_over_shift].reshape(-1)
            k += 1
    return patchs


#式(5)を計算する関数
def dictionary_learning_with_mod(W, Y):
    D = np.dot(np.dot(Y, W.T), np.linalg.inv(np.dot(W, W.T)))
    result = np.linalg.norm(Y - np.dot(D, W)) ** 2
    return D, result

#辞書クラス
class dictionary_with_learning():
    def __init__(self, weight_learner, dictionary_learner, base_num, max_loop_num = 8):
        self.Dic = np.array([])
        self.weight_learner = weight_learner
        self.dictionary_learner = dictionary_learner
        self.base_num = base_num

    #学習データ行列Yから辞書を学習するメソッド
    def learn(self, Y, eps = 1e-8):
        #手順3:ランダム行列で辞書の初期化
        D = np.random.randn(Y.shape[0], self.base_num)

        result = -1
        loop_num = 0
        while 1:
            loop_num += 1
            W = np.zeros((self.base_num, Y.shape[1]))
            #手順4:Dを固定して係数を求める
            for i in range(Y.shape[1]):
                [W[:, i], tmp] = self.weight_learner(D, Y[:, i])

            #手順5:Wを固定してDを求める
            [D, new_result] = self.dictionary_learner(W, Y)

            if result < 0 or np.abs(new_result - result) > eps or loop_num < max_loop_num:
                result = new_result
            else: #手順4、5を収束するまで繰り返す。
                break;
        self.Dic = D
        return D

    def get_dic_array(self):
        return self.Dic



#画像パッチから学習される辞書クラス
class image_patch_dictionary():
    def __init__(self, dictionary):
        self.dictonary = dictionary
        self.patch_hight = 0
        self.patch_width = 0

    #学習画像の集合から辞書を学習するメソッド
    def learn_images(self, image_list, patch_hight, patch_width, eps=2e-2, shift_num=1):
        self.patch_hight = patch_hight
        self.patch_width = patch_width
        Y = np.empty((patch_hight * patch_width, 0))
        #手順1:画像パッチの切り出し
        for img in image_list:
            gray_img = color2gray_img(img)
            Ytmp = gray_im2patch_vec(gray_img, patch_hight, patch_width,shift_num) 
            Y = np.concatenate((Y, Ytmp), axis=1)

        #手順2:学習データ行列の標準化
        Y_ = (Y - np.mean(Y, axis=1, keepdims = True))/np.std(Y, axis = 1, keepdims=True)
        self.dictonary.learn(Y_, eps)

    def get_dic(self):
        return self.dictonary.get_dic_array()

    def get_patch_hight(self):
        return self.patch_hight

    def get_patch_width(self):
        return self.patch_width


def main():
    lam = 5
    base_num = 144
    patch_hight = 12
    patch_width = 12
    train_img_list = []
    for file in find_all_files('C:\\data'):
        if os.path.isfile(file):
            try:
                train_img_list.append(Image.open(file))
                if len(train_img_list) >= 50:
                    break;
            except:
                print file

    weight_learning_method = lambda W, y:ISTA(W, y, lam) #近接勾配法を行う関数を使用
    dic_learning_method = lambda W, Y:dictionary_learning_with_mod(W, Y)#解析的に辞書Dを求める方法を使用
    im_dic = image_patch_dictionary(dictionary_with_learning(weight_learning_method, dic_learning_method, base_num))
    im_dic.learn_images(train_img_list, patch_hight, patch_width, shift_num = patch_hight)

    np.save('dic_array.npy', im_dic.dictonary.Dic)

本実装は、画像⇔行列・ベクトル変換を行う関数群、パッチ切り出しを行う関数、画像パッチ辞書クラスimage_patch_dictionary、辞書クラスdictionary_with_learningから構成されています。
image_patch_dictionaryは実際の辞書学習オブジェクトであるdictionary_with_learningをインスタンス変数として持ち、与えられた学習画像群からパッチの切り出しを行い、辞書を学習させる機能を持ちます。
そのため、辞書学習を行うのはdictionary_with_learningになりますが、係数の学習を行うweight_learnerと辞書の学習を行うdictionary_learnerをインスタンスとして持ち、生成時にどのアルゴリズムで係数推定、辞書推定を行うかを指定します。今回、係数の計算には近接勾配法を使用しますので、前回の記事記載のISTA関数をweight_learnerとして指定し、dictionary_learnerには新たに作成したdictionary_learning_with_mod関数を指定します(main関数参照)。もし係数推定として別のアルゴリズムを使いたければweight_learnerを変えればよいし、辞書学習もK-SVDなどを使いたければdictionary_learnerを変更すればよいです。

画像再構成による超解像

実装した手法の説明

ここからは、スパースコーディングによる超解像の手法の説明です。手法の大まかな考え方は「学習した辞書{\bf D}をダウンサンプリングした辞書{\bf D}_lで入力画像を符号化する係数を求めた後、その係数と元の辞書\bf Dから画像を再構成する」というものです。以下、順を追って説明します。

  1. 画像パッチの切り出し
  2. 今、入力画像(低解像画像)を{\bar{\bf Y}}とし、そこから辞書学習時と同様に画像パッチを切り出します。その切り出した画像パッチを{\bar{\bf y}_i}とします。切り出す画像パッチのサイズは、辞書学習に使用した画像(高解像画像)のパッチサイズをh \times w、低解像画像から高解像画像への倍率をsとすると、\frac{h}{s} \times \frac{w}{s}になります。ということで、{\bar {\bf y}_i}\frac{hw}{s^2}次元のベクトルになります。ここで画像パッチは、もれなくかつ領域を重複させて取得します。領域を重複させないと、画像再構成後パッチ間のつなぎ目が滑らかにならずブロック状のアーティファクト発生します。なるべく重複度が大きい方が精度は良くなりますが、取得する画像パッチも多くなるので、その分計算に時間がかかります。
  3. 辞書のダウンサンプリング
  4. 次に辞書のダウンサンプリングです。ここの方法も色々考えられますが、今回は単純に学習した{\bf D}にダウンサンプリングの作用素をかけて{\bf D}_lを得る方法をとりました。 {\bf D}は学習画像パッチのサイズhw \times基底数mの行列です。各基底{\bf d}h \times wの画像を縦に並べたベクトルと見なせるので、これをh \times wの画像に戻し、s近傍の画素を平均することでダウンサンプリング後の辞書{\bf D}_lを得ます({\bf D}_l\frac{hw}{s^2} \times mの行列になります)。
  5. 符号化係数の推定
  6. 画像パッチが切り出せたら次は、それぞれ符号化する係数を求めます。これは、ダウンサンプリング後の辞書{\bf D}_lを使い、下記のように式(3)と同様の最適化問題を解けばよいです。 \displaystyle {{\hat {\bf w}_i}} = {\rm arg} \min_{{\bf w}_i} \| {\bar {\bf y}_i} - {\bf D}_l {\bf w}_i \|^2_2  + \lambda \| {\bf w}_i \|_1\tag{6}
  7. 高解像画像パッチの生成
  8. ここまでで、各パッチを符号化するスパースな係数が求まりました。次はこの係数を用いて高解像画像を構成する画像パッチ{\bf x}_iを求めます。やることは簡単で、以下のように元々の辞書(高解像画像の辞書){\bf D}と低解像辞書で推定した係数{\hat{\bf w}_i}をかけるだけです。 \displaystyle {\bf x}_i = {\bf D}{\hat {\bf w}_i} \tag{7}
  9. パッチをつなぎ合わせて高解像画像を再構成
  10. 得られた高解像画像パッチから高解像画像{\bf X}_0を再構成します。手順1で画像パッチを領域を重複させて切り出したことに注意が必要です。今、手順1で画像を切り出す際pピクセルずらして切り出したとすると、再構成時はpsピクセルずらしてパッチを貼り合わせます。全てのパッチを貼り合わせた後、各ピクセルを加算回数分で割り平均をとります。この画像パッチのつなぎ合わせについても色々な工夫が提案されているようです。
  11. バックプロジェクションによる画像の補正
  12. 最後に仕上げです。パッチを貼り合わせ高解像画像{\bf X}_0は、パッチごとに最適化されているため、画像全体の仮定が考慮されていません。その仮定というのは、「入力画像は{\bar {\bf Y}}、高解像画像{\bf X}に、劣化やダウンサンプリングを施す作用素{\bf L}をかけることによって得られたものである。」というものです。つまり、{\bar {\bf Y}} = {\bf L} {\bf X}という仮定が前提にあり、それと推定した{\bf X}_0とのバランスをとることで、解の補正を行います。具体的には下記の最適化問題を解きます。 \displaystyle {\bf X}^* = {\rm arg} \min_{{\bf X}} \|{\bar {\bf Y}} - {\bf L} {\bf X}\|^2_2 + \mu \|{\bf X} - {\bf X_0}  \|\tag{8} ここで、\bf{X}_0, {\bar{\bf Y}}, {\bf X}は画像全体の各要素を縦に並べたベクトルです。また、今回の実験では、{\bf L}はダウンサンプリングのみを施す行列として定義しています。 この最適化は解析的に求められそうですが、多くの論文では勾配法を使って求めています。逆行列計算が大変だからでしょうか?とりあえず、今回の実装でも論文に倣って勾配法で解きました。

実装したコード

以上の内容の実装が下記です。モジュールのimportやおれおれ画像変換関数群は辞書学習時のコードと同様のものを使っています。

#バックプロジェクションで使うダウンサンプラー行列を求める関数
#ダウンサンプラーはスパースかつサイズが大きい行列なのでscipy.sparseのスパース行列型として返す
def get_downsampler_arr(from_H, from_W, to_H, to_W):
    S = sci.sparse.lil_matrix((to_H * to_W, from_H * from_W))
    k = 0
    for i in range(0, from_H, from_H/to_H):
        for j in range(0, from_W, from_W/to_W):
            s = sci.sparse.lil_matrix((from_H, from_W))
            s[i : i + from_H/to_H, j:j + from_W/to_W] = (float(to_W)/float(from_W)) * (float(to_H)/float(from_H))
            S[k, :] = s.reshape((1, from_H*from_W))
            k += 1
    return  S.tocsr()

#ベクトルデータのダウンサンプリング関数
def vec_downsampling(vec, from_H, from_W, to_H, to_W):
    arr = vec.reshape((from_H, from_W))
    ret_arr = arr_downsampling(arr, from_H, from_W, to_H, to_W)
    return ret_arr.reshape(-1)

#行列データのダウンサンプリング関数
def arr_downsampling(arr, from_H, from_W, to_H, to_W):
    ret_arr = np.zeros((to_H, to_W))
    k = 0
    for i in range(0, from_H, from_H/to_H):
        for j in range(0, from_W, from_W/to_W):
            ret_arr[i/(from_H/to_H), j/(from_W/to_W)] = (np.average(arr[i : i + from_H/to_H, j:j + from_W/to_W]))
            k += 1
    return ret_arr

#画像パッチをつなぎ合わせて画像を再構成する関数
#shift_numがずらし量、パッチがオーバーラップしている個所は足し合わせて平均化
def patch_vecs2gray_im(patchs, im_hight, im_width, patch_hight, patch_width, shift_num):
    im_arr = np.zeros((im_hight, im_width))
    sum_count_arr = np.zeros((im_hight, im_width)) #各領域何回パッチを足したか表す行列
    i = 0
    j = 0
    w_over_shift = 0
    h_over_shift = 0
    for k in range(patchs.shape[1]):
        P = patchs[:, k].reshape((patch_hight, patch_width))
        im_arr[(i* shift_num) + h_over_shift:((i* shift_num) + patch_hight) + h_over_shift,
            (j* shift_num) + w_over_shift:((j* shift_num) + patch_width) + w_over_shift] += P #指定の領域にパッチの足しこみ
        sum_count_arr[(i* shift_num) + h_over_shift :(i* shift_num) + patch_hight + h_over_shift,
            (j* shift_num) + w_over_shift:(j* shift_num) + patch_width + w_over_shift] += 1 #当該領域のパッチを足しこんだ回数をカウントアップ
        if j < ((im_width - patch_width)/shift_num + (1 if (im_width - patch_width) % shift_num  != 0 else 0)):
            j += 1
            w_over_shift = min(0, im_width - ((j* shift_num) + patch_width)) #パッチを足しこむ領域が画像からはみ出た時、そのはみ出た分を戻すための変数(横方向)
        else:
            j = 0
            i += 1
            h_over_shift = min(0, im_hight - ((i* shift_num) + patch_hight)) #パッチを足しこむ領域が画像からはみ出た時、そのはみ出た分を戻すための変数(縦方向)
            w_over_shift = 0
    im_arr /= sum_count_arr
    return float_arr2image(im_arr), im_arr

#バックプロジェクションを行う関数
def back_projection(x0, y, L, c, mu, x_init = [], eps = 1e-7):
    if x_init == []:
        xt = sci.matrix(np.random.rand(x0.shape[0], x0.shape[1])) *255.0
    else:
        xt = x_init
    result = -1
    while 1:
        x_new = xt - mu * (L.T.tocsr()*(L * xt - y) + c* (xt - x0))
        result_new = sci.linalg.norm(L * x_new - y, 2) + c * sci.linalg.norm(x_new - x0, 2)
        if np.abs(result_new - result) < eps:
            break;
        xt = x_new
        result = result_new
    return x_new, result_new

#超解像計算クラス
class super_resolution:
    def __init__(self, dictionary, weight_learner, back_projection):
        self.Dh = dictionary.get_dic()
        self.patch_hight = dictionary.get_patch_hight()
        self.patch_width = dictionary.get_patch_width()
        self.weight_learner = weight_learner
        self.back_projection_method = back_projection

    def get_super_resolution(self, image, scaling, shift_num = 1):
        #手順1:画像パッチの切り出し
        image = color2gray_img(image)
        Y = gray_im2patch_vec(image, self.patch_hight/scaling, self.patch_width/scaling, shift_num)

        #手順2:辞書のダウンサンプリング
        Dl = np.zeros((self.Dh.shape[0]/(scaling**2), self.Dh.shape[1]))
        for i in range(self.Dh.shape[1]):
            Dl[:, i] = vec_downsampling(self.Dh[:, i], self.patch_hight, self.patch_width,
                                        self.patch_hight/scaling, self.patch_width/scaling)
       
        im_arr = np.zeros(image.size)
        i = 0
        j = 0
        Patchs = np.zeros(((self.patch_hight) * (self.patch_width), Y.shape[1]))
        W = np.zeros((self.Dh.shape[1], Y.shape[1]))
        for k in range(Y.shape[1]):
            #手順3:符号化係数の推定
            [w, result] = self.weight_learner(Dl, Y[:, k])

            #手順4:高解像画像パッチの生成
            Patchs[:, k] = np.dot(self.Dh, w)

        #手順5:パッチをつなぎ合わせて高解像画像を再構成
        [reconstructed_im, X0] = patch_vecs2gray_im(Patchs, image.size[0] * scaling , image.size[1] * scaling , self.patch_hight, self.patch_width, shift_num* scaling)
        L = get_downsampler_arr(image.size[0] * scaling, image.size[1] * scaling, image.size[0], image.size[1])

        #手順6:バックプロジェクションによる画像の補正
        x = np.array(self.back_projection_method(sci.matrix(X0.reshape(-1)).T, sci.matrix(np.array(image).reshape(-1)).T, L)[0])
        return float_arr2image(x.reshape((image.size[0] * scaling, image.size[1] * scaling)))


def main():
    test_img = Image.open('Lenna.png')
    scale = 3
    lam = 50
    base_num = 144
    patch_hight = 12
    patch_width = 12
    shift_num = 1
    c = 0.8
    mu = 0.9
    weight_learning_method = lambda W, y:ISTA(W, y, lam)
    back_projection_method = lambda x0, y, L:back_projection(x0, y, L, c, mu)
    dic_learning_method = lambda W, Y:dictionary_learning_with_mod(W, Y)
    im_dic = image_patch_dictionary(dictionary_with_learning(weight_learning_method, dic_learning_method, base_num))
    im_dic.dictonary.Dic = np.load('dic_array.npy')
    im_dic.patch_hight = 12
    im_dic.patch_width = 12
    SR = super_resolution(im_dic, weight_learning_method, back_projection_method)
    result_img = SR.get_super_resolution(test_img, scale, shift_num)
    result_img.save("Lenna_SR.png")

超解像計算クラスsuper_resolutionは、画像パッチ辞書オブジェクト、係数推定関数オブジェクト、バックプロジェクション関数オブジェクトを生成時に指定し、get_super_resolutionメソッドで指定した辞書、関数に基づいた超解像画像の生成を行います。推定方法は上記で解説したとおりです。
画像パッチをつなぎ合わせて画像の再構成を行う関数patch_vecs2gray_imは、パッチのオーバーラップを許し、オーバーラップした部分は足し合わせた回数分割ることで平均をとっています。また、バックプロジェクションを行う関数back_projectionは、サイズが大きくスパースな行列Lを扱うので、scipy.matrixクラスで各引数を受け取る関数としています。ダウンサンプラー行列はget_downsampler_arr関数で求めています。

実験

まず辞書学習ですが、サイズ250×250の顔画像50枚から12×12のパッチ22000枚程度切り出して学習を行いました。また、辞書の基底数はパッチベクトルの次元と同じ144、学習時の\lambda\lambda = 5にしています。なお、辞書の学習繰り返し回数は8で打ち切っています(それ以上繰り返した学習辞書も作りましたが、生成される超解像画像の精度が逆に悪くなるという結果に・・・)。
超解像画像生成時は、\lambda=50、バックプロジェクション時のパラメータc = 0.8\mu = 0.9、パッチのずらし量1で実験しています。実験に使用する画像データはおなじみのレナさんです。

83×83の画像を249×249の画像へ(拡大率3倍)

  • 入力画像

f:id:YamagenSakam:20180411231300p:plain

f:id:YamagenSakam:20180411232337p:plain f:id:YamagenSakam:20180411231626p:plain

左の通常の拡大はPILのresizeで得られる画像です。少しわかりにくいですが、通常の拡大に比べて右の超解像画像は滑らかになっているのがわかると思います。

62×62の画像を248×248の画像へ(拡大率4倍)

  • 入力画像

f:id:YamagenSakam:20180411231846p:plain

f:id:YamagenSakam:20180411231342p:plain f:id:YamagenSakam:20180411231925p:plain

もう少しわかりやすくするため、更に拡大率を上げた画像を生成してみました。通常の拡大に比べると滑らかになってはいますが、ブロック的なノイズが乗ってしまっています。アルゴリズム的な限界なのか、辞書やパラメータを洗練すればもっとうまくいくのか、その辺りの検証は今後の課題です。

まとめ

今回は近接勾配法の応用例として辞書学習および更にその応用である超解像画像生成をやってみました。
実際やってみると、パッチのつなぎ目が汚くなったりノイズまみれになったりと結構苦労しました。その過程で色々わかったこと・思ったことを下記にまとめます(その分野では当たり前のことばかりかも知れませんが・・・)。また、やり残したこと、今後やりたいことも併せて記します。

  • やってみてわかったこと、思ったこと
    • 辞書学習時はパッチベクトルを平均{\bf 0}、分散1へ標準化した方がよい。
    • 画像再構成時はパッチ間の境界を滑らかにするため、元画像から切り出すパッチは領域を重複させて切り出し、再構成後のパッチを重ね合わせることで画像全体を再構成する方がよい。
    • 学習させすぎると精度が落ちるので、途中で学習を打ち切った方が超解像の結果がよい。
    • \lambdaによって結構結果が変わる。係数推定結果の非零要素数がパッチベクトルの要素数と大体同じになる程度が妥当か。
    • バックプロジェクションはあまり効かない?(ダウンサンプリング以外にボケなどを考慮した{\bf L}すると有効か?)
  • やり残したこと、今後やりたいこと
    • 学習用画像データのバリエーションを増やして実験。
    • 学習回数が増えると精度が落ちるのはなぜか数理的に考察。
    • 推定した画像の一部の要素が0~255の範囲を超えることがある。制約条件をつけられないか検討。
    • K-SVDなどの他の辞書学習手法の適用。
    • 独立成分分析やウェーブレットを用いた辞書行列の適用。
    • 他のダウンサンプリング辞書生成方法の適用。
    • 辞書学習以外の超解像手法との比較。


参考文献

*1:実務等で利用するときはscikit-learn辺りのライブラリを使うべきですが、今回は理解を深めるためnumpy、scipy、PILを使った実装を行っています。

*2:実は{\bf D}正規分布に従ったランダムな行列でも、符号化としてはある程度うまくいくという研究もあります。ただ超解像やその他画像再構成への応用という観点ではどうなんでしょう?