スパースモデリングに基づく画像の再構成 Part2. Total Variation最小化(Split Bregman)に基づく画像再構成 - IMACEL Academy -人工知能・画像解析の技術応用に向けて-| エルピクセル株式会社

Total Variation の正則化理論

TV denoising is considered to be one of the best denoising models, but also one of the hardest to compute.
via TV正則化に関する記述． ftp://ftp.math.ucla.edu/pub/camreport/cam08-29.pdf

TV正則化に関する記述．
ftp://ftp.math.ucla.edu/pub/camreport/cam08-29.pdf

Total Variation とは，例えば画像を復元したいと思ったとき，画像の微分のL1ノルムの最小化を行う方法です． L1ノルムの最小化はスパースな解を誘発するため，微分のL1をとることで， 画像の輝度値の変化がスパースな解を得ることが出来ます．これを確かめるため，ノイズの乗った画像に対して復元を行います．

結果だけ述べると，左のような画像から右のような画像を復元することが可能です．

Total Variation による方法では，微分のL1ノルムの最小化の正則化項を用います．このときの評価関数は

\begin{align} F = |\nabla_x u| + |\nabla_y u| + \frac{\lambda}{2}|u - I|^2 \end{align}

です．第一項目は画像の画素値$I$と，再構成画像$u$との差分，第二項は$\nabla_x$，$\nabla_y$でピクセルごとの差分です．

上の項は，微分のスパースを仮定しているため，画素値の変化がより小さいような解が得られることになります．具体的に図に表現するとこのようになります．

この記事では，上の最適化問題をどのように解くのか解説をし，pythonを用いて実装いたします．

最適化方法: Split Bregman

上の問題は$u$の微分に関するL1ノルムの項が含まれているため，非線形な方程式です．

そのため，反復計算を行う必要があります．

まずはじめに，上の最適化問題を以下のように同値変形します．

\begin{align} F = |d_x| + |d_y| + \frac{\lambda}{2}|u - I|^2 \end{align} \begin{align} {\rm s.t.} \ \ d_x = \nabla_x u, \ \ d_y = \nabla_y u \end{align}

ここで，制約条件を少し弱め，(To weakly enforce the constraints in this formulation と原文には書いてある)，最適化問題を以下のように変更します．

\begin{align} \underset{u,d_x,d_y}{\rm minimize}\ \ |d_x| + |d_y| + \frac{\lambda}{2}|u - I|^2 + \frac{\mu}{2}(|\nabla_x u - d_x|^2 + |\nabla_y u - d_y|^2) \end{align}

そして，これに対して，Bregman Iterationというのを用います．これは，上の最適化問題を更新のタイムステップ$k$を用いて \begin{align} \underset{u,d_x,d_y}{\rm minimize}\ \ |d_x| + |d_y| + \frac{\lambda}{2}|u - I|^2 + \frac{\mu}{2}(|\nabla_x u - d_x - b_x^k|^2 + |\nabla_y u - d_y - b_y^k|^2) \end{align} とし，最適化問題として \begin{align} u^{k+1} = \underset{u}{\rm minimize}\ \ |d_x| + |d_y| + \frac{\lambda}{2}|u - I|^2 + \frac{\mu}{2}(|\nabla_x u - d_x - b_x^k|^2 + |\nabla_y u - d_y - b_y^k|^2) \end{align} とすることで，$u$を更新するものです．

このように分割することによって，2次式の最適化と絶対値を含む問題の最適化に帰着します． 2次式の場合は最適化は簡単に求めることができますし，絶対値を含む問題に対しては解析的な解が求まっています．

詳しくは，L1ノルム最適化に関して書かれた以下の記事をご覧ください．

スパースモデリングに基づく画像の再構成 Part1. L1ノルム最小化に基づく画像再構成の実装 - IMACEL Academy -人工知能・画像解析の技術応用に向けて-|LPixel(エルピクセル)

この記事では，L1ノルム正則化の最小化の実装を行い，ノイズを含む画像がどのように再構成されるのか，確かめてみます．なお，Total Variation はスパースモデリングで主に使われている技術です．

実装

それではpython で実装していきます．

パッケージのインポート

import numpy as np
import cv2
import matplotlib.pyplot as plt
import random

パッケージのインポート.py

def add_noise(I_t):
    mu = np.mean(I_t)
    sigma = np.std(I_t)
    dB = 3
    I_noise = 10**(-dB/20)*np.reshape([random.gauss(mu, sigma) for i in range(np.size(I_t))], np.shape(I_t))
    I = I_t + I_noise
    max_I  = np.max(I)
    min_I = np.min(I)
    I = np.round((I - min_I)*255/(max_I - min_I))
    return I

ノイズを付与する関数.py

画像のロード

img_path = "geo2.png"
img_load = cv2.imread(img_path)
I_t = cv2.cvtColor(img_load, cv2.COLOR_RGB2GRAY)
f = add_noise(I_t)
[X_N,Y_N] = np.shape(f)

画像のロード.py

出力

plt.subplot(1,2,1)
plt.imshow(I_t)
plt.title("Original")
plt.subplot(1,2,2)
plt.imshow(f)
plt.title("Noisy")

出力.py

ノイズデータの生成

Split Bregman で使う関数

def add_noise(I_t):
    mu = np.mean(I_t)
    sigma = np.std(I_t)
    dB = 3
    I_noise = 10**(-dB/20)*np.reshape([random.gauss(mu, sigma) for i in range(np.size(I_t))], np.shape(I_t))
    I = I_t + I_noise
    max_I  = np.max(I)
    min_I = np.min(I)
    I = np.round((I - min_I)*255/(max_I - min_I))
    return I


def Gauss_Saidel(u, d_x, d_y, b_x, b_y, f, MU, LAMBDA):
    U = np.hstack([u[:,1:X_N], np.reshape(u[-1,:],[Y_N,1] )]) + np.hstack([np.reshape(u[0,:],[Y_N,1]), u[:,0:Y_N-1]]) \
       + np.vstack([u[1:X_N,:], np.reshape(u[:,-1],[1,X_N] )]) + np.vstack([np.reshape(u[:,0],[1,X_N] ), u[0:X_N-1,:]])
    D = np.vstack([np.reshape(d_x[:,0],[1,X_N] ), d_x[0:Y_N-1,:]]) - d_x \
       + np.hstack([np.reshape(d_y[0,:],[Y_N,1] ), d_y[:,0:X_N-1]]) - d_y
    B = -np.vstack([np.reshape(b_x[:,0],[1,X_N] ), b_x[0:Y_N-1,:]]) + b_x \
       - np.hstack([np.reshape(b_y[0,:],[Y_N,1] ), b_y[:,0:X_N-1]]) + b_y
    G = LAMBDA/(MU + 4*LAMBDA)*(U+D+B) + MU/(MU + 4*LAMBDA)*f
    return G
    

def shrink(x,y):
    t = np.abs(x) - y
    S = np.sign(x)*(t > 0) * t
    return S

関数.py

Split Bregman による画像再構成

def main():
    ## Load Image
    #img_path = "Blood2.jpg"
                 
    CYCLE = 100
    MU = 5.0*10**(-2)
    LAMBDA = 1.0*10**(-2)
    TOL = 5.0*10**(-1)

    ## Initialization
    u = f
    d_x = np.zeros([X_N,Y_N])
    d_y = np.zeros([X_N,Y_N])
    b_x = np.zeros([X_N,Y_N])
    b_y = np.zeros([X_N,Y_N])

    for cyc in range(CYCLE):
        u_n = Gauss_Saidel(u,d_x,d_y, b_x ,b_y,f, MU,LAMBDA)
        Err = np.max(np.abs(u_n[2:X_N-2,2:Y_N-2] - u[2:X_N-2,2:Y_N-2]))
        if np.mod(cyc,10)==0:
            print([cyc,Err])
        if Err < TOL:
            break
        else:
            u = u_n
            nablax_u = np.vstack([u[1:X_N,:], np.reshape(u[:,-1],[1,X_N] )]) - u 
            nablay_u = np.hstack([u[:,1:X_N], np.reshape(u[-1,:],[Y_N,1] )]) - u 
            d_x = shrink(nablax_u + b_x, 1/LAMBDA)
            d_y = shrink(nablay_u + b_y, 1/LAMBDA)
            b_x = b_x + (nablax_u - d_x)
            b_y = b_y + (nablay_u - d_y)
    
    
    ## plot figure
    plt.figure()
    plt.subplot(1,2,1)
    plt.gray()
    plt.imshow(f)
    plt.title('Noisy')
    plt.axis("off")
    
    plt.subplot(1,2,2)
    plt.gray()
    plt.imshow(np.round(u))
    x1, y1 = [0,X_N], [50,50]
    plt.plot(x1, y1)
    plt.title('Reconstructed')
    plt.axis("off")
    
    plt.figure()
    plt.subplot(2,1,1)
    plt.plot(f[50,:])
    plt.subplot(2,1,2)
    plt.plot(u[50,:])

if __name__ == "__main__":
    main()

SplitBreg.py

結果

ノイズデータ（左）と再構成画像（右）

正直ここまで出来るのはびっくりしました．
盛んに研究されているだけありますね．

断面の様子

別の画像で行った結果

結論

正直ここまで再構成ができるとは意外でした．しかし，画像処理だけでここまで綺麗な画像が出力できるのは非常に興味深いですね．皆様もぜひ実装してみてください．

参考文献

[1] THE SPLIT BREGMAN METHOD FOR L1 REGULARIZED PROBLEMS, TOM GOLDSTEIN, STANLEY OSHER, ftp://ftp.math.ucla.edu/pub/camreport/cam08-29.pdf