どせいたんさき。

ナスダヨー

テキトウノイズリダクション

目的

統計的なモデルを用いてノイズリダクションをする.

やりかた

ノイズが載った画像として 40×40 px の画像を以下のような感じで用意する.

img = randn(40) + 3*sin(0.005*[1:40]'*[1:40]*pi);

randn は平均 0 分散 1 の乱数で今回はこれをノイズとして扱う.画像のノイズでない部分は数ピクセルにわたる大きい構造を持っている.また,ノイズでない部分には切り立った不連続な面は存在しないとする.つまり 1 ピクセルスケールのばらつきはほぼノイズだとみなしてよい.このことを利用して画像の微分値の 2 乗和を少なくするように変形すればノイズを除去することができる.画像のノイズ成分を noise としたときに画像を微分して 2 乗和を計算するための(最小化すべき)式は以下の通り.

sumsq(diff((img-noise)(:)))

なおこの式では画像の微分は一方向だけしか考えていない.ただし,この式を最小化してしまうとノイズだけでなく画像の情報全てが失われてしまう.微分の 2 乗和は画像が一定値のときに最も少なくなるからである.ここで,ノイズ成分は平均 0 分散 1 の乱数によって生成されたことを利用する.ノイズ成分はあまり大きい値をとるはずがない.よって大きい値を取らないように最小化すべき関数にペナルティを与える.

sumsq(diff((img-noise)(:))) + sumsq(noise(:))

このペナルティ部分はノイズ成分の事前確率に相当している.リッジ回帰の式と同じような感じになっている.どちらの項が相対的に重要なのか重み付けをする必要がある気がするのだけどとりあえず忘れておく.この式を sqp という関数を用いて最小化する.sqp とは逐次二次計画法によって非線形な関数を最小化するための関数である. sqp は二次元配列を引数にとる関数を扱えないので,最小化すべき関数は以下のように定義する.

func = @(n) sumsq(diff(img-reshape(n,40,40))(:)) + sumsq(n);

n は画像のノイズ成分を 1 次元ベクトルに変形したものである.これを sqp に食わせれば適切なノイズを推定してくれるはずである.以下のコマンドを実行してノイズ成分を得る.

noise = sqp(zeros(1600,1), func);
rimg  = img - noise; ## ノイズ成分を除去した画像

結果

なんとなくノイズがとれてのっぺりしたような気もする.ただし実行が遅すぎる.自分の環境では 40×40 の 2 次元配列を処理するのに 6 分ほどかかった.ノイズじゃない成分をどうやって判定するかはもっと工夫ができそう.


左:元画像 中:処理後 右:ノイズなし