超解像処理アルゴリズム

　超解像処理　2011.10.9 HimagingQ ＠ Morinosato
ここではコンボリューションぼけした画像の復元を行うデコンボリューション超解像処理
に関してプログラミングノウハウを説明する。
「ぼかし処理アルゴリズム」の中で説明した移動平均はフラットな正方形カーネルパタン
とのコンボリューションであったが、そのようなパタンのことを Point Spread Function と
呼び、PSF を仮定してぼけ画像復元する超解像処理の基本がデコンボリューションである。
このデコンボリューション超解像処理アルゴリズムのプログラミングノウハウは、
1. デコンボリューションを含む FFT ベースの基本機能をライブラリ関数に仕立てる。
ベースとなる FFT は高速な MKL@Intel 等の外部ライブラリ利用が望ましく、ここで
は libfftw3f の利用に留め、SIMD 命令を用いた最適化プログラミングはあえて行わず 、
FFT ベースの基本機能として、PSF パタン画像生成、コンボリューション、 デコンボ
リューション、さらに正規化相関パタンマッチング、位相相関パタンマッチングおよ
び空間周波数フィルタといった標準的 C++関数を作成する。
2. 画像端での不連続性に起因するアーチファクトの抑制が可能な処理方法を決める。
アーチファクトが発生しないように、入力画像をフェードアウト付きミラーリング
外枠付けで拡大したワーク画像を関数内部で生成し、 FFT 逆変換後のワーク画像の外枠
を除いた部分を入力画像に上書き出力する。
3. 微小値での割り算に起因するノイズの抑制が可能な処理方法を決める。
分母が 0.1 未満ならば分母に 0.01 を加えることによって過大なノイズが発生しない
ようにできる。この方法によれば、バイラテラルフィルタで充分ノイズカットできる
程度にまでノイズ低減されたデコンボリューション画像を得ることが可能である。
説明はこの辺にして以下にコーディング例（ IDE は Qt Creator）を示す。
int kernelSize;
// カーネルサイズ N
Qimage* inputImage;
// 入力画像（Format は QImage::Format_Indexed8 とする）
int imageWidth = inputImage->width();
int iW2 = inputImage->bytesPerLine();
int imageHeight = inputImage->height();
uchar* pS = inputImage->bits();
int iF = kernelSize / 2;
int iw = imageWidth + iF * 2;
int ih = imageHeight + iF * 2;
QImage workImage(iw , ih, QImage::Format_Indexed8);
// ワーク画像
workImage.setColorTable(*palette V);
uchar* pBits = workImage.bits();
int iW = workImage.b ytesPerLine();
uchar* pD = pBits + iW * iF;
double v = pow(iF * 0.5, 2.0) * 2.0;
double* weight = new double[iF];
for (int k = 0; k < iF; k++)
weight[k] = exp( - (k * k) / v);
int i, j;
for (j = 0; j < imageHeight; j++, pD += iW, pS += iW2)
{
for (i = 0; i < iF; i++)
pD[i] = (uchar)(pS[iF - i] * weight[iF - i]);
memcpy(pD + iF, pS, imageWidth);
for (i = 0; i < iF; i++)
pD[iF + imageWidth + i ] = (uchar)(pS [imageWidth - 2 - i] * weight[i]);
}
pD = pBits;
for (j = 0; j < iF; j++)
{
for (i = 0; i < iW; i++)
pD[j * iW + i] = (uchar)(pD[(iF * 2 - j) * iW + i] * weight[iF - j]);
for (i = 0; i < iW; i++)
pD[(iF + imageHeight + j) * iW + i] = (uchar)(pD[(iF + imageHeight - 2 - j) * iW + i] * weight[j]);
}
//　//　以上で入力画像のコピーをフェードアウト付きミラーリングで外枠付け拡大したワーク画像が生成できた。
//　//　そこで、ワーク画像の外枠を除いた部分のデコンボリューション画像を入力画像に上書き出力する。
//　Deconvolution(memor yP, &workImage);
pS = inputImage->bits();
pD = pBits + iW * iF;
for (j = 0; j < imageHeight; j++, pD += iW, pS += iW2)
memcpy(pS, pD + iF, imageWidth); ［ 入力画像 → ワーク画像 ］の内部生成例（フェードアウト付きミラーリング） ：
↓
このようにフェードアウトさ せた外枠付けを行っておくことにより、単なる０パディング
によって、画像端での連続性が考慮された FFT 変換前画像を内部生成することができる。
以下は前記コード最終段の Deconvolution 関数の C 言語によるコーディング例である。
void Limaging::Decon volution(QImage* psfImage , QImage* inputImage)
{
int iws = psfImage->width();
int ihs = psfImage->height();
int iw = inputImage->width();
int ih = inputImage->height();
int ncWidth = i2Exp(iw);
int ncHeight = i2Exp(ih);
FCOMPLEX* pCi = NewCbits(inputImage , ncWidth, ncHeight);
FCOMPLEX* pCis = NewCbits(psfImage , ncWidth, ncHeight);
uint dwCount = ncWidth * ncHeight;
uint k;
float ds2 = 0.0f;
for (k = 0; k < dwCount; k++)
ds2 += pCis [k].u;
float fScale = 1.0f / ds2;
for (k = 0; k < dwCount; k++)
pCis[k].u *= fScale;
FFT2D(pCi, ncWidth, ncHeight, 1);
FFT2D(pCis, ncWidth, ncHeight, 1);
float ftmp, u, v;
for (k = 0; k < dwCount; k++)
{
ftmp = pCis [k].u * pCis[k].u + pCis[k].v * pCis[k].v;
if (ftmp < 0.1f)
ftmp += 0.01f;
u = pCi[k].u;
v = pCi[k].v;
pCi[k].u = (u * pCis [k].u + v * pCis[k].v) / ftmp;
pCi[k].v = (v * pCis[k].u - u * pCis[k].v) / ftmp;
}
fftwf_free(pCis);
FFT2D(pCi, ncWidth, ncHeight, -1);
fScale = 1.0f / (float)dwCount;
for (k = 0; k < dwCount; k++)
pCi[k].u *= fScale;
OutputBits(inputImage, pCi, ncWidth, ncHeight, - iws / 2, - ihs / 2);
fftwf_free(pCi);
}
ここで、引数の inputImage には外枠付けした workImage を渡さないとアーチファクトが
発生してしまう。一方、引数の psfImage には外枠付け無しで PSF パタン画像を渡さなければ
ならない。
なお、上記コーディング例中の変数や関数については別掲のソースコードを参照されたい。
［処理例］３頁の画像のコンボリューション（左）のデコンボリューション（右）：
ここで、PSF パタンはカーネルサイズ９のガウシアン（冪乗係数＝1.0）：
ぼけ回復度を向上させるため種々の改良手法が提案されているが、カーネルサイズ５以下
であれば概ね満足できる視覚効果が得られる。例えば、微妙に散乱ぼけしたＸ線透過画像に
対しては、カーネルサイズ＝５、冪乗係数＝1.５の（1.5 乗した）ガウシアンで PSF パタン
を推定近似して、デコンボリューションすると満足できる効果が得られるであろう。