👻

JPEGとスパースモデリング(2)

2021/05/26に公開

画像処理

実装

さていよいよ実装ですが、詳細はGithubにおいているのでそちらを御覧ください。

複数のブロックに分けていますが、全てimage_sparseクラスの実装です。

class image_sparse:
    def __init__(self,N):
        self.N = N
        self.phi = np.zeros((N**2,N**2))
        self.sparse_num = 0
        for i in range(N**2):
            for j in range(N**2):
                self.phi[i][j] = self.phi_k_i(j//N,i//N)*self.phi_k_i(j%N,i%N)

    def phi_k_i(self,k,i):
        if k == 0:
            return 1 / math.sqrt(self.N)
        else:
            return math.sqrt(2/self.N)*math.cos(((2*i+1)*k*math.pi)/(2*self.N))

コンストラクタです。
コンストラクタが受け取るNは、圧縮処理の分割のサイズを格納する変数です。今回は $N=8$ となっており、 $N\times N$ のサイズに分割して処理を行います。
sparse_numは、最後どれくらいスパースなデータかを表す変数です。
前回割愛した行列 $\psi$ の中身ですが、やはりややこしいことになっています。
(記事では $\psi$ ですがプログラムでは $\phi$ となっていました、読み替えてください)

    def compress(self,img,lam):
        img_N_division = []
        x,y = img.shape
        if x != y or x %8 != 0 or y %8 != 0:
            raise Exception("Error : Error : image must be square, and its each size must be divided 8")
        
        for i in range(y//8):
            for j in range(x//8):
                img_N_division.append(img[self.N*i:self.N*(i+1),self.N*j:self.N*(j+1)])

        for i in range(y//8):
            for j in range(x//8):
                dct_coef_N_division,compressed_img_N_division = self.calculate(img_N_division[i*(y//N)+j],lam)

                if j == 0:
                    compressed_img_row = compressed_img_N_division
                    dct_coef_row = dct_coef_N_division
                else:
                    compressed_img_row = np.concatenate([compressed_img_row, compressed_img_N_division], 1)
                    dct_coef_row = np.concatenate([dct_coef_row, dct_coef_N_division], 1)

                print("\r","Progress,",i*(x//8)+j+1," / ",(y//8)*(x//8),"",end='')

            if i == 0:
                compressed_img = compressed_img_row
                dct_coef = dct_coef_row
            else:
                compressed_img = np.concatenate([compressed_img,  compressed_img_row])
                dct_coef = np.concatenate([dct_coef,  dct_coef_row])
            

        print("")
            
        compressed_img = correct_abnormal_value(compressed_img.astype(dtype=int))
        dct_coef = dct_coef.astype(dtype=int)
        return dct_coef,compressed_img

    def calculate(self,img,lam):
        dct_coef = self.opt_dct_coef(img,lam)
        self.sparse_num += self.get_num_nonzero_in_arr(dct_coef)
        invdct_coef = self.idct(dct_coef)
        
        return dct_coef,invdct_coef

compress(image, lam)
画像全体を受け取り、変換後の画像全体を返す関数
calculate(image,lam)
分割された画像を受け取り、DCT係数と変換後の分割画像を返す関数

compress関数は、長々と書いてありますがほとんどが前処理と後処理です。
やっていることは画像を分割しcalculate関数に渡しているだけです。

    def get_num_nonzero_in_arr(self,arr):
        arr = arr.flatten()
        return len(arr) - np.count_nonzero(arr)

    def opt_dct_coef(self,img,lam):
        #行列Wを生成
        img -= 128
        W = self.phi

        #変数xを定義
        x = Variable(self.N ** 2)
        y = img.flatten()
        #ハイパーパラメータを定義
        lamda = Parameter(nonneg=True)
        lamda.value =lam

        #目的関数を定義
        objective = Minimize(sum_squares((y - W @ x))/(2*N)+ lamda*norm(x, 1))
        p = Problem(objective)

        #最適化計算
        result = p.solve()
        dct_coef = np.round(x.value)
        return np.array(np.split(dct_coef,8)) 
    
    def idct(self,dct_coef):
        invdct = np.dot(self.phi, dct_coef.flatten())
        return np.round(np.split(invdct,8)) + 128

opt_dct_img(image,lam)
DCT係数をL1正則化を使用し求め返す関数
最適化の実装にはcvxpyを使いました。
idct(dct_coef)
DCT係数を受け取り、逆変換した値を返す関数

ひどい解説ですが、とりあえずメインの実装はこんなところです。

結果

前回の冒頭部分で大まかな結果は述べたので、ここでは簡単に。結果の画像を再掲します。なお、実際の出力ではPSNRやSSIMなどの指標の値も出力していますが省略します。

現画像	JPEG	L1( $\lambda$ =0.01)	L1( $\lambda$ =0.1)	L1( $\lambda$ =1)	L1( $\lambda$ =5)
0%	89%	21%	43%	87%	96%

ここまで記事を読んでいただいた方はお分かりかと思いますが、JPEGの劣化は高周波成分の消失です。この例でいうと、顔の凹凸によってできる影の部分や、髪と背景の境界部分にノイズがかかっています。
L1正則化で、同等の圧縮率であった $\lambda=1$ の場合と比べてみると、L1正則化の場合全体がぼやけているように見え、JPEGのほうが詳細が見えやすいかと思います。これはL1正則化が周波数といった人間の画像認識能力を考慮せず最適な解を求めているからです。ただ、前回も言いましたが必ずしもJPEGがいいと言われるわけではありません。境界部分の詳細が知りたい場合やノイズをどうしても嫌がる場合など、状況に応じて画像圧縮のアルゴリズムも選択されるべきでありその選択肢一つとして今回実装した手法が採用される可能性もあるでしょう。

学校の一課題でしたが、思いの外時間を力を注ぐことになりました。

実装

結果

Discussion