Python/PILで画像にノイズをかけ各種フィルタで除去する

PIL特有の機能を使わずにピクセル操作から自力実装する。
Python Imaging Library Handbook

流れ

ノイズ付加

  • ごま塩ノイズ(Salt-pepper Noise)
  • ガウシアンノイズ?(Gaussian Noise)→ よくわからなかった

ノイズ除去

  • 移動平均法(Mean Filter)
  • ガウシアンフィルタ(Gaussian Filter)
  • 中央値フィルタ(Median Filter)
ノイズ付加

ごま塩ノイズ

ピクセルを走査して確率的に白黒の点をまぶす。

from PIL import Image
from itertools import product
from random import random

class Noise:
    def __init__(self, input_image):
        self.input_image = input_image
        self.input_pix = self.input_image.load()
        self.w, self.h = self.input_image.size

    def saltpepper(self, salt=0.05, pepper=0.05):
        output_image = Image.new("RGB", input_image.size)
        output_pix = output_image.load()

        for x, y in product(*map(range, (self.w, self.h))):
            r = random()
            if r < salt:
                output_pix[x, y] = (255, 255, 255)
            elif r > 1 - pepper:
                output_pix[x, y] = (  0,   0,   0)
            else:
                output_pix[x, y] = self.input_pix[x, y]
        return output_image

縮小しているので分かりづらいが右画像は白と黒の点がランダムに打たれている。

Original Salt-pepper Noise(PSNR=15.19)

ガウシアンノイズ?

これがよくわからなかった。いやちゃんと調べろよと言われると返す言葉もない。
それぞれの輝度を平均として、決め打ちの増幅値(amplification)だけ分散する正規乱数を求めると解釈した。

noised_colors = map(lambda x: gauss(x, amp), self.input_pix[x, y])

なんかこのあたりの考え方微妙に分節できていない。最初はこんなものだと勘違いしていた。

noised_colors = map(lambda x: x + (random() - 0.5) * amp, self.input_pix[x, y])

Noise.gaussianの実装は以下のようにした。ガウス乱数を求める関数はrandomモジュールの中に入っている。

from random import gauss

class Noise:
    # ...省略...

    def gaussian(self, amp=10):
        output_image = Image.new("RGB", input_image.size)
        output_pix = output_image.load()
        for x, y in product(*map(xrange, (self.w, self.h))):
            noised_colors = map(lambda x: gauss(x, amp), self.input_pix[x, y])
            noised_colors = map(lambda x: max(0, x), map(lambda x: min(255, x), noised_colors))
            noised_colors = tuple(map(int, noised_colors))
            output_pix[x, y] = noised_colors
        return output_image
Original Gaussian Noise?(amp=100, PSNR=10.27)

ノイズ付加はここまで。

フィルタリング

移動平均

ある画素を選んだ時、その周辺の画素の値で平均した値をその画素の輝度値とする手法。
以下のサイトが画像つきで説明しており分かりやすい。ガウシアンフィルタも同じ。

平滑化(移動平均、ガウシアン)フィルタ 画像処理ソリューション

コードの工夫が今ひとつ。後に出てくるmedian filterと共通部分が多い。

class Filter:
    k = tuple([(1/9, 1/9, 1/9) for _ in range(3)])

    def __init__(self, input_image):
        self.input_image = input_image
        self.input_pix = self.input_image.load()
        self.w, self.h = self.input_image.size

    def mean(self, kernel=k):
        # カーネルの総和は1である必要がある
        summation = sum([sum(k) for k in kernel])
        if summation != 1:
            return None

        output_image = Image.new("RGB", input_image.size)
        output_pix = output_image.load()

        # あるピクセルの周辺を走査するための一時変数(きたない)
        W = int(len(kernel[0])/2)
        W = (-W, W+1)
        H = int(len(kernel)/2)
        H = (-H, H+1)

        # 配列の境界判定
        borderX = {0: (0, 2), self.w-1: (-1, 1)}
        borderY = {0: (0, 2), self.h-1: (-1, 1)}

        for y, x in product(*map(xrange, (self.h, self.w))):
            weight_pixels = []
            _W, _H = W, H

            if not 0 < x < self.w - 1:
                _W = borderX[x]
            if not 0 < y < self.h - 1:
                _H = borderY[y]

            for i, j in product(enumerate(range(*_W)), enumerate(range(*_H))):
                # i => (kernelのindex, xからの相対位置)
                # j => (kernelのindex, yからの相対位置)
                i, dx = i
                j, dy = j

                # 重み付き色の計算
                weight_pixel = map(lambda x: x * kernel[i][j], self.input_pix[x+dx, y+dy])
                weight_pixels.append(weight_pixel)

            color = tuple(map(int, map(sum, zip(*weight_pixels))))
            output_pix[x, y] = color
        return output_image

ガウシアンノイズ(たぶん)を加えた画像に平均値フィルタを繰り返し適用した様子を示す。
ノイズは減るが輪郭がどんどんぼやけているのがわかる。

Gaussian Noise?
(amp=100, PSNR=10.27)
Mean Filter
(1回適用, PSNR=20.00)
Mean Filter
(10回適用, PSNR=20.12)

ガウシアンフィルタ

ある画素を選んだ時、その周辺画素の情報を利用する、という点で移動平均フィルタと考え方は同じ。
ピクセルからの距離に応じて平滑化の重み付けを行うのだが、いちいちすべて計算していては時間がかかりすぎるので適当なカーネルを用いる。

class Filter:
    gk = (( 1 /16, 1 / 8, 1 /16),
          ( 1 / 8, 1 / 4, 1 /8 ),
          ( 1 /16, 1 / 8, 1 /16))
    # ...省略...
    def gaussian(self, gauss_kernel=gk):
        return self.mean(gauss_kernel)

このカーネルだと移動平均フィルタと見た目上あまり区別がつかないので画像は省略。

メディアンフィルタ

周辺の色を輝度順にソートして中央値を採用するフィルタ。ゴマ取りに有効。

class Filter:
    # ...省略...
    def median(self):
        output_image = Image.new("RGB", input_image.size)
        output_pix = output_image.load()

        # ピクセルの周辺走査用変数
        W = (-1, 2)
        H = (-1, 2)

        # 配列の境界判定
        borderX = {0: (0, 2), self.w-1: (-1, 1)}
        borderY = {0: (0, 2), self.h-1: (-1, 1)}

        for y, x in product(*map(xrange, (self.h, self.w))):
            pixels = []
            _W, _H = W, H

            if not 0 < x < self.w - 1:
                _W = borderX[x]
            if not 0 < y < self.h - 1:
                _H = borderY[y]

            # x, y周辺の輝度を得る
            for dx, dy in product(range(*_W), range(*_H)):
                pixel = self.input_pix[x+dx, y+dy]
                pixels.append(pixel)

            # 輝度配列の長さ
            pixlen = len(pixels)

            # 色成分ごとに配列をまとめる
            pixels = zip(*pixels)

            # 輝度順にソートする
            sorted_pixels = map(sorted, pixels)

            # 中央値を得る
            color = tuple(map(lambda x: x[int(pixlen/2)], sorted_pixels))
            output_pix[x, y] = color
        return output_image

Salt-pepper Noiseを加えた画像に適用した結果。よーく見るとちょっとだけゴマが残っていたりするがほとんど綺麗になっている。また輪郭もくっきり残っている。

Salt-pepper Noise(PSNR=15.19) Median Filter(PSNR=32.49)

ピーク信号対雑音比(PSNR

オリジナル画像から比較画像がどれくらい劣化したかを測る指標だそうです。単位はdB。
この値が大きければ大きいほど比較画像の品質は良いらしい。
定義はリンク参照。

ピーク信号対雑音比 - Wikipedia
PSNR - 裏紙

def PSNR(original, compare):
    from math import log10, sqrt
    MAX = 2**8 - 1

    w, h = original.size
    opix = original.load()
    cpix = compare.load()

    color_sum = 0
    for x, y in product(*map(range, (w, h))):
        color = zip(*(opix[x, y], cpix[x, y]))
        color = map(lambda x: (x[0] - x[1])**2, color)
        color_sum += sum(color)
    MSE = color_sum/3/w/h
    PSNR = 20 * log10(MAX/sqrt(MSE))
    return PSNR

Bilateral-FilterとNon-local Means Filterわかりたい。