こしあん
2019-10-05

PyTorchでガウシアンピラミッド+ラプラシアンピラミッド(Gaussian/Laplacian Pyramid)

Pocket
LINEで送る
Delicious にシェア

3.9k{icon} {views}


Progressive-GANの論文で、SWD(Sliced Wasserstein Distance)が評価指標として出てきたので、その途中で必要になったガウシアンピラミッド、ラプラシアンピラミッドをPyTorchで実装してみました。これらのピラミッドはGAN関係なく、画像処理一般で使えるものです。応用例として、ラプラシアンブレンドもPyTorchで実装しています。

関連

OpenCVのドキュメントが非常に参考になりました。こちらを参考にCNNのConvolutionに置き換えて実装しました。

画像ピラミッド
http://labs.eecs.tottori-u.ac.jp/sd/Member/oyamada/OpenCV/html/py_tutorials/py_imgproc/py_pyramids/py_pyramids.html

またここに出てくる実装は、自分の下記の記事がベースになっています。必要なら参考にしてみてください。

ガウシアンピラミッド、ラプラシアンピラミッド

まず画像のピラミッドとは何かというと、もとの画像の解像度を1/2、1/2、1/2……のように繰り返し縮小してまとめたものがピラミッドです。縮小する過程でガウシアンぼかしを入れたものがガウシアンピラミッド、ガウシアンピラミッドに解像度間の差を取ったものがラプラシアンピラミッドとなります。

ラプラシアンピラミッドはバンドパスフィルタのようなことをやっているので、画像の周波数分解とも考えることができます。OpenCVのドキュメントからです。


上がガウシアンピラミッドで、下がラプラシアンピラミッドです。大きい画像が高周波数成分、低い画像が低周波数成分です。発想としてはOctave ConvolutionのOctaveと同じです。

以下にアルゴリズムの詳細を示します。

ガウシアンピラミッド

以下のようなコードで示される5×5の畳み込みカーネルを用意します。

import numpy as np

kernel = np.array([
    [1, 4, 6, 4, 1],
    [4, 16, 24, 16, 4],
    [6, 24, 36, 24, 6],
    [4, 16, 24, 16, 4],
    [1, 4, 6, 4, 1]], np.float32) / 256.0

OpenCVのドキュメントだと一部分母が16になっていましたが、256にするのが正しいかと思われます(普通の畳み込みカーネルも行列の全要素の和が分母になってます)。

OpenCVでは、「5×5カーネルの畳み込み→偶数の行と列を削除する」という処理を実行しています。しかし、ディープラーニングの畳み込み関数(Conv2D)ではもっと簡潔に表せて、「5×5カーネル、stride=2、padding=2」というパラメーターで(ほぼ)同じ処理を実行できます。固定値のカーネルをどう計算するのかはこちらを参照してください。

つまり、PyTorchでガウシアンピラミッドを実装するには、上で表されるような5×5カーネルによるConv2dを連続で実行するだけということになります。

ラプラシアンピラミッド

ラプラシアンピラミッドは少し複雑です。ガウシアンピラミッドで紹介した処理は解像度を下げる方向なので「pyramid_down」と表すことにします。ラプラシアンピラミッドでは、解像度を上げる「pyramid_up」という処理が必要になります。大きな流れは次の通りです。

  1. ガウシアンピラミッド(pyramid_down)で計算されたピラミッドのi番目の画像について考える
  2. i番目の画像に対し、pyramid_downをしたi+1番目のガウシアンピラミッドに注目する
  3. i+1番目のガウシアンピラミッドをpyramid_upし、1.で示したi番目の画像との差分を取る

ラプラシアンピラミッドのやっていることは周波数帯ごとの差分計算なので、「ラプラシアンピラミッドの個数=ガウシアンピラミッドの個数-1」となります。

次は「pyramid_up」はどういう処理をするかということです。

  • 入力画像を2倍の解像度にアップサンプリングする(後でガウシアンぼかしを入れることから、Nearest Neighbor法で良いと思われる)
  • pyramid_downで使ったのと同じ畳み込みカーネルでガウシアンぼかしをかける

これでラプラシアンピラミッドのアルゴリズムは終わりです。

コード

from PIL import Image
import numpy as np
import torch
import torch.nn.functional as F
import torchvision

# 画像ファイル→PyTorchテンソル
def load_tensor(img_path):
    with Image.open(img_path) as img:
        array = np.asarray(img, np.float32).transpose([2, 0, 1]) / 255.0
        tensor = torch.as_tensor(np.expand_dims(array, axis=0))  # rank 4
    return tensor

# ピラミッドを1枚の画像に結合して保存するための関数
def tile_pyramid(imgs):
    height, width = imgs[0].size()[2:]
    with torch.no_grad():
        canvas = torch.zeros(1, 3, height * 3 // 2, width)
        x, y = 0, 0
        for i, img in enumerate(imgs):
            h, w = img.size()[2:]
            canvas[:,:, y:(y + h), x:(x + w)] = img            
            if i % 2 == 0:
                x += width // (2 ** (i + 3))                    
                y += height // (2 ** i) # 0, 2, 4..でy方向にシフト
            else:
                x += width // (2 ** i)  # 1, 3, 5..でx方向にシフト
                y += height // (2 ** (i + 3))
        return canvas

# ガウシアンぼかしのカーネル
def get_gaussian_kernel():
    # [out_ch, in_ch, .., ..] : channel wiseに計算
    kernel = np.array([
        [1, 4, 6, 4, 1],
        [4, 16, 24, 16, 4],
        [6, 24, 36, 24, 6],
        [4, 16, 24, 16, 4],
        [1, 4, 6, 4, 1]], np.float32) / 256.0
    gaussian_k = torch.as_tensor(kernel.reshape(1, 1, 5, 5))
    return gaussian_k

def pyramid_down(image):
    with torch.no_grad():
        gaussian_k = get_gaussian_kernel()        
        # channel-wise conv(大事)
        multiband = [F.conv2d(image[:, i:i + 1,:,:], gaussian_k, padding=2, stride=2) for i in range(3)]
        down_image = torch.cat(multiband, dim=1)
    return down_image

def pyramid_up(image):
    with torch.no_grad():
        gaussian_k = get_gaussian_kernel()
        upsample = F.interpolate(image, scale_factor=2)
        multiband = [F.conv2d(upsample[:, i:i + 1,:,:], gaussian_k, padding=2) for i in range(3)]
        up_image = torch.cat(multiband, dim=1)
    return up_image

def gaussian_pyramid(original, n_pyramids):
    x = original
    # pyramid down
    pyramids = [original]
    for i in range(n_pyramids):
        x = pyramid_down(x)
        pyramids.append(x)
    return pyramids

def laplacian_pyramid(original, n_pyramids):
    # gaussian pyramidを作る
    pyramids = gaussian_pyramid(original, n_pyramids)

    # pyramid up - diff
    laplacian = []
    for i in range(len(pyramids) - 1):
        diff = pyramids[i] - pyramid_up(pyramids[i + 1])
        laplacian.append(diff)
    return laplacian

# 見やすいようにMin-Maxでスケーリングする
def normalize_pyramids(pyramids):
    result = []
    for diff in pyramids:
        diff_min = torch.min(diff)
        diff_max = torch.max(diff)
        diff_normalize = (diff - diff_min) / (diff_max - diff_min)
        result.append(diff_normalize)
    return result

ラプラシアンピラミッドの出力結果を見やすくするために、Min-MaxスケーリングでNormalizeした関数も入れています。

結果:ガウシアンピラミッド

三江線に登場していただきましょう。これをtrain.jpgとします。

if __name__ == "__main__":
    original = load_tensor("train.jpg")
    pyramid = gaussian_pyramid(original, 6)
    torchvision.utils.save_image(tile_pyramid(pyramid), "gaussian_pyramid.jpg")    

ガウシアンピラミッドは次のようになります。

画像がだんだん小さくなっているのでわかりづらいですが、ガウシアンぼかしをかけているのでだんだんぼやけてくるのがわかります。

結果:ラプラシアンピラミッド

同じようにラプラシアンピラミッドを作ってみます。元画像のピクセル値が0~1の場合、ラプラシアンピラミッドは差分を取っているので、-1~1までの値が取り得ます。そこで、Min-Max Scalerを適用し、出力画像が0~1の範囲になるようにします。

if __name__ == "__main__":
    original = load_tensor("train.jpg")
    pyramid = laplacian_pyramid(original, 6)
    normalize = normalize_pyramids(pyramid)
    torchvision.utils.save_image(tile_pyramid(normalize), "laplacian_pyramid.jpg")

この結果はとてもおもしろいです。細かい輪郭線が高周波数帯に、逆にテクスチャや色は低周波数帯に集中しているのがわかります。

これはCNNとのアナロジーとしても理解できます。CNNは浅い層では解像度が高く(高周波)、深い層では解像度が低い(低周波)構成となっています。同時に浅い層では主に輪郭やエッジのような波形的情報を見ており、深い層ではテクスチャや何が映っているかという空間的・意味的な情報を見ていると言われています。

ラプラシアンピラミッドによって、CNNでなくても画像そのものが高周波・低周波で同様の構成になっており、CNNの構造はまさにそれにフィットしている――だからCNNが画像に向いている、ということが可視化できるともいえるでしょう。

ちなみにOpenCVの実装では、subtract(引き算)関数でオーバーフロー抑制機能が機能し、マイナスの値を0に丸められています。したがって、このPyTorchの実装とOpenCVの実装では結果が少し違います。ただし、マイナスの値を丸めないほうが明らかに情報は落ちないので、ここではそのまま「ただの引き算」として処理しています。後述のラプラシアンブレンドでも、丸めないほうが自然な出力になります。

ラプラシアンブレンド

応用例として、OpenCVで紹介されているラプラシアンブレンドをPyTorchで実装してみました。

こちらのサイトから、「apple.jpg」と「orange.jpg」をダウンロードしてきます。

# ラプラシアンブレンディング
def laplacian_blending():
    apple = load_tensor("apple.jpg")
    orange = load_tensor("orange.jpg")
    apple_lap = laplacian_pyramid(apple, 5)
    orange_lap = laplacian_pyramid(orange, 5)

    # ラプラシアンピラミッドを左右にブレンドする
    blend_lap = []
    for x, y in zip(apple_lap, orange_lap):
        width = x.size(3)
        b = torch.cat([x[:,:,:,:width // 2], y[:,:,:, width // 2:]], dim=3)
        blend_lap.append(b)

    # 最高レベルのガウシアンピラミッドのブレンド
    apple_top = gaussian_pyramid(apple, 5)[-1]
    orange_top = gaussian_pyramid(orange, 5)[-1]
    out = torch.cat([apple_top[:,:,:,:apple_top.size(3) // 2],
                     orange_top[:,:,:, orange_top.size(3) // 2:]], dim=3)

    # ラプラシアンピラミッドからの再構築
    for lap in blend_lap[::-1]:
        out = pyramid_up(out) + lap
    torchvision.utils.save_image(out, "laplacian_blend.png")

    # 比較例:ダイレクトにブレンド
    direct = torch.cat([apple[:,:,:,:apple.size(3) // 2],
                        orange[:,:,:, orange.size(3) // 2:]], dim=3)
    torchvision.utils.save_image(direct, "direct_blend.png")

結果は以下の通りです。

ダイレクトブレンド(元画像を左右に貼り付けただけ)

ラプラシアンブレンド

ダイレクトブレンドはつぎはぎ目が明らかに不自然ですが(クソコラ感がすごい)、ラプラシアンブレンドは滑らかになっていますね。

まとめ

PyTorchだろうとラプラシアンピラミッドは実装できました。これを使えばSWD(Sliced Wasserstein Distance)は多分計算できると思います。ラプラシアンピラミッドは単なる画像特徴量としても面白そうです。



Shikoan's ML Blogの中の人が運営しているサークル「じゅ~しぃ~すくりぷと」の本のご案内

技術書コーナー

【新刊】インフィニティNumPy――配列の初期化から、ゲームの戦闘、静止画や動画作成までの221問

「本当の実装力を身につける」ための221本ノック――
機械学習(ML)で避けて通れない数値計算ライブラリ・NumPyを、自在に活用できるようになろう。「できる」ための体系的な理解を目指します。基礎から丁寧に解説し、ディープラーニング(DL)の難しいモデルで遭遇する、NumPyの黒魔術もカバー。初心者から経験者・上級者まで楽しめる一冊です。問題を解き終わったとき、MLやDLなどの発展分野にスムーズに入っていけるでしょう。

本書の大きな特徴として、Pythonの本でありがちな「NumPyとML・DLの結合を外した」点があります。NumPyを理解するのに、MLまで理解するのは負担が大きいです。本書ではあえてこれらの内容を書いていません。行列やテンソルの理解に役立つ「従来の画像処理」をNumPyベースで深く解説・実装していきます。

しかし、問題の多くは、DLの実装で頻出の関数・処理を重点的に取り上げています。経験者なら思わず「あー」となるでしょう。関数丸暗記では自分で実装できません。「覚える関数は最小限、できる内容は無限大」の世界をぜひ体験してみてください。画像編集ソフトの処理をNumPyベースで実装する楽しさがわかるでしょう。※紙の本は電子版の特典つき

モザイク除去から学ぶ 最先端のディープラーニング

「誰もが夢見るモザイク除去」を起点として、機械学習・ディープラーニングの基本をはじめ、GAN(敵対的生成ネットワーク)の基本や発展型、ICCV, CVPR, ECCVといった国際学会の最新論文をカバーしていく本です。
ディープラーニングの研究は発展が目覚ましく、特にGANの発展型は市販の本でほとんどカバーされていない内容です。英語の原著論文を著者がコードに落とし込み、実装を踏まえながら丁寧に解説していきます。
また、本コードは全てTensorFlow2.0(Keras)に対応し、Googleの開発した新しい機械学習向け計算デバイス・TPU(Tensor Processing Unit)をフル活用しています。Google Colaboratoryを用いた環境構築不要の演習問題もあるため、読者自ら手を動かしながら理解を深めていくことができます。

AI、機械学習、ディープラーニングの最新事情、奥深いGANの世界を知りたい方にとってぜひ手にとっていただきたい一冊となっています。持ち運びに便利な電子書籍のDLコードが付属しています。

「おもしろ同人誌バザールオンライン」で紹介されました!(14:03~) https://youtu.be/gaXkTj7T79Y?t=843

まとめURL:https://github.com/koshian2/MosaicDeeplearningBook
A4 全195ページ、カラー12ページ / 2020年3月発行

Shikoan's ML Blog -Vol.1/2-

累計100万PV超の人気ブログが待望の電子化! このブログが電子書籍になって読みやすくなりました!

・1章完結のオムニバス形式
・機械学習の基本からマニアックなネタまで
・どこから読んでもOK
・何巻から読んでもOK

・短いものは2ページ、長いものは20ページ超のものも…
・通勤・通学の短い時間でもすぐ読める!
・読むのに便利な「しおり」機能つき

・全巻はA5サイズでたっぷりの「200ページオーバー」
・1冊にたっぷり30本収録。1本あたり18.3円の圧倒的コストパフォーマンス!
・文庫本感覚でお楽しみください

北海道の駅巡りコーナー

日高本線 車なし全駅巡り

ローカル線や秘境駅、マニアックな駅に興味のある方におすすめ! 2021年に大半区間が廃線になる、北海道の日高本線の全区間・全29駅(苫小牧~様似)を記録した本です。マイカーを使わずに、公共交通機関(バス)と徒歩のみで全駅訪問を行いました。日高本線が延伸する計画のあった、襟裳岬まで様似から足を伸ばしています。代行バスと路線バスの織り成す極限の時刻表ゲームと、絶海の太平洋と馬に囲まれた日高路、日高の隠れたグルメを是非たっぷり堪能してください。A4・フルカラー・192ページのたっぷりのボリュームで、あなたも旅行気分を漫喫できること待ったなし!

見どころ:日高本線被災区間(大狩部、慶能舞川橋梁、清畠~豊郷) / 牧場に囲まれた絵笛駅 / 窓口のあっただるま駅・荻伏駅 / 汐見の戦争遺跡のトーチカ / 新冠温泉、三石温泉 / 襟裳岬

A4 全192ページフルカラー / 2020年11月発行


Pocket
LINEで送る
Delicious にシェア

Add a Comment

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です