Pythonで画像をWavelet変換するサンプル


ちょっとウェーブレット変換に興味が出てきたのでどんな感じなのかを実際に動かして試してみました。

必要なもの
--------

以下の３つが必要です。pip などで入れましょう。

* PyWavelets 
* numpy 
* PIL

簡単な解説
--------

* [PyWavelets](http://www.pybytes.com/pywavelets/contents.html) というライブラリを使っています。
* 離散ウェーブレット変換（と逆変換）、階層的な？ウェーブレット変換（と逆変換）をやってくれます。他にも何かできそうです。
* ２次元データ（画像）でやる場合は、縦横サイズが同じじゃないと上手くいかないです（やり方がおかしいだけかもしれませんが）

サンプルコード
-------

[gist:image_wavelet_transform.py](https://gist.github.com/4690081)

```python:image_wavelet_transform.py
# coding: utf8
# 2013/2/1 mokemokechicken@github.com
"""ウェーブレット変換のイメージを掴むためのサンプルスクリプト   

Require: pip install PyWavelets numpy PIL

Usage: python image_wavelet_transform.py <filename> (<level>:=3) (wavelet:=db1)
"""

import sys
from PIL import Image
import pywt, numpy

filename = sys.argv[1]
LEVEL = len(sys.argv) > 2 and int(sys.argv[2]) or 3
WAVLET = len(sys.argv) > 3 and sys.argv[3] or "db1"


def merge_images(cA, cH_V_D):
    """numpy.array を ４つ(左上、(右上、左下、右下))くっつける"""
    cH, cV, cD = cH_V_D
    print cA.shape, cH.shape, cV.shape, cD.shape
    cA = cA[0:cH.shape[0], 0:cV.shape[1]] # 元画像が2の累乗でない場合、端数ができることがあるので、サイズを合わせる。小さい方に合わせます。
    return numpy.vstack((numpy.hstack((cA,cH)), numpy.hstack((cV, cD)))) # 左上、右上、左下、右下、で画素をくっつける

def create_image(ary):
    """numpy.array を Grayscale画像に変換する"""
    newim = Image.new("L", ary.shape)
    newim.putdata(ary.flatten())
    return newim

def wavlet_transform_to_image(gray_image, level, wavlet="db1", mode="sym"):
    """gray画像をlevel階層分Wavelet変換して、各段階を画像表現で返す

    return [復元レベル０の画像, 復元レベル１の画像,... , 復元レベル<level-1>の画像, 各2D係数を１枚の画像にした画像]
    """
    ret = []
    data = numpy.array(list(gray_image.getdata()), dtype=numpy.float64).reshape(gray_image.size)
    images = pywt.wavedec2(data, wavlet, level=level, mode=mode) # http://www.pybytes.com/pywavelets/ref/2d-dwt-and-idwt.html
    for i in range(2, len(images)+1): # 部分的に復元して ret に詰める
        ary = pywt.waverec2(images[0:i], WAVLET)  * 2**(i-1) / 2**level # 部分的に復元すると加算されていた値が戻らない(白っぽくなってしまう)ので調整
        ret.append(create_image(ary))
    # 各2D係数を１枚の画像にする
    merge = images[0] / (2**level) # cA の 部分は値が加算されていくので、画像表示のため平均をとる
    for i in range(1, len(images)): 
        merge = merge_images(merge, images[i]) # ４つの画像を合わせていく
    ret.append(create_image(merge))
    return ret

if __name__ == "__main__":
    im = Image.open(filename)
    if im.size[0] != im.size[1]: # 縦横サイズが同じじゃないとなんか上手くいかないので、とりあえず合わせておく
        max_size = max(im.size)
        newim = Image.new(im.mode, (max_size, max_size))
        newim.paste(im, (0,0))
        im = newim

    im.getdata() # これを呼んでおかないと何故か次の split() が失敗する・・・？なぜ？
    bands = im.split() # RGBの各チャネル毎に処理をする
    converted_bands_array = [wavlet_transform_to_image(gray, LEVEL, wavlet=WAVLET) for gray in bands] # 各RGBチャネル毎に変換
    # zip(*hoge) がわかりにくいけど、 converted_bands には (R,G,B)の画像(PIL.Image)が入る. mergeでRGB画像に復元
    converted_array = [Image.merge(im.mode, converted_bands) for converted_bands in zip(*converted_bands_array)]

    # converted_array: [復元レベル０の画像, 復元レベル１の画像,... , 復元レベル<level-1>の画像, 各2D係数を１枚の画像にした画像]
    for i, img in enumerate(converted_array): 
        img.save("%s_%d.png" % (filename, i)) # 適当に画像出力
```

実行結果例
--------

元画像(kuma_005.jpg)
![kuma_005.jpg](https://fbcdn-sphotos-h-a.akamaihd.net/hphotos-ak-snc7/483015_595226860504296_261849812_n.jpg)

係数を画像化した画像(level=5)(kuma_005.jpg_5.png)
![kuma_005.jpg_5.png](https://fbcdn-sphotos-g-a.akamaihd.net/hphotos-ak-prn1/533853_595227067170942_1832968474_n.jpg)