Moiz's journal

プログラミングやFPGAなどの技術系の趣味に関するブログです

ゼロから作るRAW現像 その7- ノイズモデルとノイズフィルター

このブログ記事「ゼロから作るRAW現像」を大きく再構成してより読みやすくした書籍「PythonとColabでできる-ゼロから作るRAW現像」を【技術書典6】にて頒布しました。

現在はBOOTHにて入手可能です。書籍+PDF版は2200円プラス送料、PDF版は1200円です。

moiz.booth.pm

はじめに

これは「ゼロから作るRAW現像 」という一連の記事の一つです。 これらの記事の内容を前提としていますので、まだお読みでない方はこちらの記事からお読みいただくことをおすすめします。

「ゼロから作るRAW現像 その1 - 基本的な処理」

「ゼロから作るRAW現像 その2 - 処理のモジュール化」

「ゼロから作るRAW現像 その3 - デモザイク処理基本編」

「ゼロから作るRAW現像 その4 - デモザイク処理応用編」

「ゼロから作るRAW現像 その5 - ラズベリーパイのRAW画像処理」

「ゼロから作るRAW現像 その6 - レンズシェーディング補正」

デジタル画像のノイズ

デジタル画像にも当然ながらノイズはあります。

ノイズ源には様々な物がありますが、ごくごく大雑把に言うと、センサーから画像データを読み出す際に付加されるリードノイズと、センサーの各セルに入ってくる光の量の統計的ゆらぎによる光ショットノイズが大きな要因になります。

もちろんこれは非常に大雑把な分類で、リードノイズには、熱ショットノイズ、固定パターンノイズ、暗電流、などがあります。

光ショットノイズは光の量に依存して、その分散は次のようなポアッソン分布を取る事が知られています。

{ \displaystyle
\sigma ^{2} \left( n \right) = n
}

本来、光ショットノイズは光の量のみに依存するノイズのはずですが、実在する画像センサーの特性が理想的なものとは異なるために(非線形性やクロストークなど)、実際の観測では複雑な特性を見せたりします。

結局、ノイズ量も実際のセンサーやカメラ系で測定してカリブレーションする事になります。

最近のカメラ画像処理では、こういったカリブレーションによるノイズ量の推定はノイズモデルと呼ばれ重要視されています。

実際の画像のノイズ

今回の内容はgithubからJupyter Notebookのノートとしてダウンロードできます

それでは前回使用した画像のノイズを観察してみましょう。

f:id:uzusayuu:20181021123223p:plain:w800

この画像を拡大すると、このようになっています。

f:id:uzusayuu:20181028101459p:plain:w800

なんだかざらざらしています。なるべく明るくして撮影したのですが、まだだいぶノイズがのっているようです。

では実際どの程度の量のノイズがあるのか測定してみましょう。

f:id:uzusayuu:20181028102234p:plain:w300

図のグレイパッチ部分(赤い長方形で囲った部分)のノイズ量を実際に測定してみます。

まず画像を前回同様の方法で読み込みます。

import os, sys, imageio, math, scipy
module_path = os.path.abspath(os.path.join('..'))
if module_path not in sys.path:
    sys.path.append(module_path)
import raw_process4 as raw_process
from pylab import imshow, show
from matplotlib import pyplot as plt
from matplotlib.pyplot import imshow
import numpy as np

with open("chart.jpg", "rb") as input_file:
    data = input_file.read()
data = data[-10237440:]
w = 3282  # for v1.3 w = 2592
h = 2480 - 32  # for v1.3 h = 1944

img = np.zeros((h, w))
stride = math.ceil(w * 10 / 8 / 32) * 32
for y in range(h):
    for x in range(w // 4):
        word = data[y * stride + x * 5: y * stride + x * 5 + 5]
        img[y, 4 * x    ] = (word[0] << 2) | ((word[4] >> 6) & 3) 
        img[y, 4 * x + 1] = (word[1] << 2) | ((word[4] >> 4) & 3) 
        img[y, 4 * x + 2] = (word[2] << 2) | ((word[4] >> 2) & 3) 
        img[y, 4 * x + 3] = (word[3] << 2) | ((word[4]     ) & 3)

表示します。

outimg = img.copy()
outimg[outimg < 0] = 0
outimg[outimg > 1023] = 1023
outimg = outimg / 1024
imshow(outimg, cmap='gray')

f:id:uzusayuu:20181028102815p:plain

ブラックレベルだけ補正しておきましょう。

blacklevel = [64] * 4
bayer_pattern = np.array([[2, 1], [1, 0]])
blc_raw = raw_process.black_level_correction(img, blacklevel)

outimg = blc_raw.copy()
outimg = outimg.reshape((h, w))
outimg[outimg < 0] = 0
outimg = outimg / outimg.max()
imshow(outimg, cmap='gray')

f:id:uzusayuu:20181028102953p:plain

いつもの画像が得られました。

各グレイパッチの座標を画像ビューワーなどで測定しておきます。 今回は次のような座標でした。

patches = [(2586, 2086), (2430, 2092), (2272, 2090), (2112, 2090), (1958, 2086), (1792, 2094), 
           (1642, 2096), (1486, 2090), (1328, 2090), (1172, 2086), (1016,2084), (860, 2084),
           (866, 482), (1022, 480), (1172, 476), (1328, 474), (1480, 470), (1634, 466),
          (1788, 462), (1944, 460), (2110, 452), (2266, 452), (2424, 448), (2586, 442)]

これらの座標((2586, 2086)など)はグレイパッチ内にとった100x100画素の正方形領域の左上の位置です。

各パッチ内の画素の分散と平均値を測定してみます。 画像のフォーマットがBayerなので、各色チャンネル毎に統計をとります。

variances = []
averages = []
for index, (dx, dy) in enumerate(((0, 0), (1, 0), (0, 1), (1, 1))):
    for patch in patches:
        x, y = patch
        p = blc_raw[y+dy:y+100:2, x+dx:x+100:2]
        s2 = (p * p).mean()
        av = p.mean()
        v = s2 - av * av
        variances.append(v)
        averages.append(av)

測定結果を見てみましょう。

plt.plot(averages, variances, linestyle='None', marker='o', color='blue')
plt.show()

f:id:uzusayuu:20181028103609p:plain

どうやら、分散は画素の値に対してほぼ線形になっているようです。

f:id:uzusayuu:20181028121839p:plain

ここからこの画像では光ショットノイズが支配的だと考えられます。

ここから先はこのノイズを取り除いてくことにします。

ノイズフィルターの定番バイラテラルフィルター

ノイズフィルターに求められる特性として、ノイズは取り除いてほしいが、元の画像に含まれる情報は残しておきたい、というものがあります。

このような矛盾する要求に答えるものとしてバイラテラルフィルター1があります。

バイラテラルフィルターのアイデアは基本的に、

  1. ある画素の周辺の画素のうち、値が近いものは同じものだからフィルターをかける
  2. 値が遠いものは違うものだからフィルターをかけない。

というものです。

たとえば、各ターゲット画素のまわりに次のような5x5の領域を設定します。

f:id:uzusayuu:20181028105501p:plain

赤で囲った画素がターゲットです。

周りの画素のうち、灰色の画素はターゲットに近い値を持ち、黒い画素は大きく異なる値を持つとします。

この場合、灰色の部分は画像のある部分(たとえば新聞の地の部分)、黒い部分は他の部分(たとえば新聞の印刷部分)に含まれると考えられます。

そうなると、灰色の画素にのみフィルターをかけ、黒い画素にはフィルターをかけないことで、画像の特徴に影響を与えずにノイズフィルターをかける事ができます。

一例としてはこのようなフィルターになるでしょう。(これは説明用の図で、実際のフィルターとは異なります。)

f:id:uzusayuu:20181028134837p:plain

バイラテラルフィルターでは、ターゲット画素と周辺画素の一つ一つの値の差を取り、それをノイズの量とくらべて加重平均のウェイトを計算することで、ターゲット画素に近いものにのみフィルターを適用します。

式としては、例えばこうなります。

{ \displaystyle
p_{\rm{out}} \left( x, y \right) =
\frac{\sum_{dy=-N}^{+N} \sum_{dx=-N}^{+N}  w \left( dx, dy \right) p_{\rm{in}} \left( x + dx, y + dy \right)}
{\sum_{dy=-N}^{+N} \sum_{dx=-N}^{+N}  w \left( dx, dy \right)}
}

{ \displaystyle
w \left( dx, dy \right) = exp \left( - \frac { a |p_{\rm{in}} \left(x + dx, y + dy \right) - p_{\rm{in}} \left( x, y \right) |^{2} + b | dx^{2} + dy^{2} |} {\sigma ^{2}} \right)
}

ただし、{ p_{\rm{in}}}が入力画像で、{ p_{\rm{out}}}が出力画像、{\left(x, y \right)}が座標です。{N}はウィンドウのサイズ、{a}はチューニングパラメータで通常1以上(4から8など)です。

この式によって、ターゲット画素と周辺画素は{w}の計算の中で比べられます。差が{\sigma}に比べて大きければ{w}は小さくなり、その画素のウェイトは小さくなります。

逆にターゲットの差が小さい場合は{w}の値は1に近くなり、大きなウェイトを持つことになります。

もし周辺の画像が全てターゲット画像に近い場合、フィルターの形はガウシアンフィルターに近づいていきます。

一つ注意として、ノイズ分散の扱いがあります。

通常の画像処理の教科書ではノイズ分散は画像内で一定という仮定を行います。 もしリード・ノイズが支配的ならこの仮定は正しいのですが、先程見たとおり実際の画像ではショットノイズのためにノイズの分散は画素の値に依存します。

もしノイズフィルターをかけるのがカメラ画像処理後のRGBやYUVの画像だとすると、すでにガンマ補正やローカルトーンマップ補正などの処理が行われていて、元のノイズの分布を推定するのは困難です。 RAW画像から処理する場合はこのような処理が行われていないことはわかっているので、ショットノイズなどのノイズ量を推定するのは簡単です。

次の節以降では測定したノイズの特性を利用してノイズ処理を行っていきます。

バイラテラルノイズフィルターの適用

ではこのノイズフィルターを実際の画像に適用してみましょう。

Bayerフォーマットのままではフィルターがかけにくいので、レンズシェーディング、ホワイトバランス、デモザイクを実行しフルカラー画像にします。

# Lens shading correction
par = [np.array([6.07106808e-07, 9.60556906e-01]), 
       np.array([6.32044369e-07, 9.70694361e-01]), 
       np.array([6.28455183e-07, 9.72493898e-01]), 
       np.array([9.58743579e-07, 9.29427169e-01])]
gain_map = np.zeros((h, w))
center_y, center_x = h // 2, w // 2
for y in range(0, h, 2):
    for x in range(0, w, 2):
        r2 = (y - center_y) * (y - center_y) + (x - center_x) * (x - center_x)
        gain = [par[i][0] * r2 + par[i][1] for i in range(4)]
        gain_map[y, x] = gain[0]
        gain_map[y, x+1] = gain[1]
        gain_map[y+1, x] = gain[2]
        gain_map[y+1, x+1] = gain[3]
lsc_raw = blc_raw * gain_map

# White Balance
wbg = np.array([1.128, 1, 2.546, 1])
wb_raw = raw_process.white_balance_Bayer(lsc_raw, wbg, bayer_pattern)

# Demosaic
dms_img = raw_process.advanced_demosaic(wb_raw, bayer_pattern)

一応確認します。

outimg = dms_img.copy()
outimg = outimg.reshape((h, w, 3))
outimg[outimg < 0] = 0
outimg[outimg > 1023] = 1023
outimg = outimg / 1024
imshow(outimg)

f:id:uzusayuu:20181028122045p:plain

さて、このような処理を行った後でも、ノイズの分散は画素の値に比例するのでしょうか? RGBの平均に対して確かめてみましょう。

luma_img = dms_img[:, :, 0] + dms_img[:, :, 1] + dms_img[:, : , 2]
variances = []
averages = []
for patch in patches:
    x, y = patch
    p = luma_img[y:y+100, x:x+100]
    s2 = (p * p).mean()
    av = p.mean()
    v = s2 - av * av
    variances.append(v)
    averages.append(av)

plt.plot(averages, variances, linestyle='None', marker='o', color='blue')
plt.show()

f:id:uzusayuu:20181028122243p:plain

先ほどとは多少様子が違いますが、線形近似で大丈夫そうです。

傾きを求めてみましょう。

par = np.polyfit(averages, variances, 1)
print(par[1])

246.0281307576092

さあ、これでノイズフィルターを適用する準備ができました。実際にかけてみましょう。

今回はフルカラーなので、ウェイトの計算はRGBの平均に対して行い、フィルターの適用は各カラーごとに行うという方法を使っています。 コード中、coefがバイラテラルフィルターの式の{a}に相当します。なお今回は{b=0}としました。

# 注:これは処理をわかりやすく書いたもので非常に実行速度が遅い。
# 実際にはこの後の高速版を使用することをおすすめする
coef = 8
img_flt = dms_img.copy()
for y in range(2, h-2):
    for x in range (2, w - 2):
        # 5x5の平均値からノイズの分散(sigma) を推定する
        average = luma_img[y-2:y+3, x-2:x+3].mean()
        sigma = par[1] * average
        sigma = sigma if sigma > 0 else 1

        weight = np.zeros((5, 5))
        out_pixel = np.zeros(3)
        norm = 0
        # 5x5内の各画素毎に重みを計算する
        for dy in range(-2, 3):
            for dx in range(-2, 3):
                # 中心画素との差
                diff = luma_img[y + dy, x + dx] - luma_img[y, x]
                diff_norm = diff * diff / sigma
                # 差と分散からウェイトを計算し、加重平均値を求める
                weight = math.exp(-coef * diff_norm)
                out_pixel += weight * dms_img[y + dy, x + dx, :]
                norm += weight
        # 各色毎にウェイトの和で正規化する
        img_flt[y, x, 0] = out_pixel[0] / norm
        img_flt[y, x, 1] = out_pixel[1] / norm
        img_flt[y, x, 2] = out_pixel[2] / norm
outimg = img_flt.copy()
outimg = outimg.reshape((h, w, 3))
outimg = outimg / 1024
outimg[outimg < 0] = 0
outimg[outimg > 1] = 1
imshow(outimg)

f:id:uzusayuu:20181028123615p:plain

このままではわかりにくいので、残りの処理(カラーマトリクスとガンマ補正)を行います。

ccm_matrix = (6022,-2314,394,-936,4728,310,300,-4324,8126)
img_ccm = raw_process.color_correction_matrix(img_flt, ccm_matrix)
img_ccm[img_ccm > 1023] = 1023
img_gamma = raw_process.gamma_correction(img_ccm, 2.2)

outimg = img_gamma.copy()
outimg = outimg.reshape((h, w, 3))
outimg[outimg < 0] = 0
imshow(outimg)

f:id:uzusayuu:20181028123917p:plain

保存します。

raw_process.write(img_gamma, "raspi_filtered_out.png")

このファイルを画像ビューワーなどを使って拡大してノイズフィルターをかける前と比較してみるとこうなっていました。 右がノイズフィルター後です。

f:id:uzusayuu:20181028130202p:plain:w800

エッジなど画像の特徴を残したまま、平坦部分のノイズが減っている事が確認できます。

成功です。

補足:ノイズフィルター処理の高速化

最後にもう一点だけ、処理の高速化について触れます。

上記のノイズフィルターのコードはプログラムとしては動作しますが、非常に遅いコードです。

これはPythonの性能上しかたのない部分もあるのですが、numpyなどの機能を利用することでかなり改善できます。2

まず、一般的な傾向としてnumpyではループの処理は遅いので、なるべくforループを減らしたほうが高速になる場合が多いです。 numpyにはこのような用途のためにstride_tricksというライブラリ3がありますので、これを利用してループを減らしていきましょう。

まずは、RGBの平均画像(luma_img)から分散を計算します。

from numpy.lib.stride_tricks import as_strided
coef = 8

average = scipy.ndimage.filters.uniform_filter(luma_img, 5, mode='mirror')
sigma_map = average * par[1]
sigma_map[sigma_map < 1] = 1
sy, sx = sigma_map.strides
sigma_tile = as_strided(sigma_map, strides=(sy, sx, 0, 0), shape=(h, w, 5, 5))
sigma_tile = sigma_tile[2:h-2, 2:w-2, : , :]

ここではまず、3278 x 2444 通りの5x5のパッチについて、FIRフィルターを使って平均値を計算し、そこからノイズ分散をもとめています。

次にそのノイズ分散値をコピーして同じ要素の5x5の行列を作り、その5x5の行列を3278 x 2444個並べています。

同様に、各パッチの平均値をコピーして同じ要素の5x5の行列を作り、その5x5の行列を3278 x 2444個並べます。

sy, sx = luma_img.strides
luma_tile = as_strided(luma_img, strides=(sy, sx, 0, 0), shape=(h, w, 5, 5))
luma_tile = luma_tile[2:h-2, 2:w-2, : , :]

次に、RGBの平均画像(luma_img)から5x5のパッチを3278 x 2444通り作ります。

sy, sx = luma_img.strides
luma_box = as_strided(luma_img, strides=(sy, sx, sy, sx), shape=(h-4, w-4, 5, 5))

次に、5x5パッチ内のウェイトとその総和を3278 x 2476個分、一挙に計算します。

diff = luma_box - luma_tile
weight = np.exp(-coef * diff * diff / sigma_tile)
weight_sum = weight.sum(axis=(2, 3))

これで各色を処理する準備ができました。

まず、赤画素を処理します。3278 x 2444個のパッチを一度に処理します。

red = dms_img[:, :, 0]
sy, sx = red.strides
red_boxes = as_strided(red, strides=(sy, sx, sy, sx), shape=(h-4, w-4, 5, 5))
red_out = (weight * red_boxes).sum(axis=(2, 3)) / weight_sum

同様に、青と緑の画素を処理します。

green = dms_img[:, :, 1]
sy, sx = green.strides
green_boxes = as_strided(green, strides=(sy, sx, sy, sx), shape=(h-4, w-4, 5, 5))
green_out = (weight * green_boxes).sum(axis=(2, 3)) / weight_sum

blue = dms_img[:, :, 2]
sy, sx = blue.strides
blue_boxes = as_strided(blue, strides=(sy, sx, sy, sx), shape=(h-4, w-4, 5, 5))
blue_out = (weight * blue_boxes).sum(axis=(2, 3)) / weight_sum

すべてまとめて完成です。

img_flt2 = dms_img.copy()
img_flt2[2:h-2, 2:w-2, 0] = red_out
img_flt2[2:h-2, 2:w-2, 1] = green_out
img_flt2[2:h-2, 2:w-2, 2] = blue_out

outimg = img_flt2.copy()
outimg = outimg.reshape((h, w, 3))
outimg = outimg / 1024
outimg[outimg < 0] = 0
outimg[outimg > 1] = 1
imshow(outimg)

f:id:uzusayuu:20181105012811p:plain

一応残りの処理も行って確認しておきます。

ccm_matrix = (6022,-2314,394,-936,4728,310,300,-4324,8126)
img_ccm = raw_process.color_correction_matrix(img_flt2, ccm_matrix)
img_ccm[img_ccm > 1023] = 1023
img_gamma = raw_process.gamma_correction(img_ccm, 2.2)

outimg = img_gamma.copy()
outimg = outimg.reshape((h, w, 3))
outimg[outimg < 0] = 0
imshow(outimg)

f:id:uzusayuu:20181105012949p:plain

先ほどと同様の画像が出力されたようです。

これでノイズ処理からループが一掃されました。 処理の速度も数倍になり、実用的になりました。

まとめ

今回はノイズフィルターをとりあげ、エッジを残すノイズフィルターとしてよく使われるバイラテラル・フィルターを解説し、実装しました。

バイラテラルフィルターは移り変わりの激しい画像処理の分野では、もはや古典的ともいえるアルゴリズムです。 アルゴリズムの発表から時間が立つとはいえその基本的な考え方は現在のノイズフィルターにも受け継がれています。

とはいえ、次々に、より高性能なノイズフィルターアルゴリズムが提案されているのは間違いなく、性能もどんどん進化しています。 今回の内容で物足りない方は、BM3D4ディープラーニングを用いたノイズフィルター5など、さらに現代的なノイズフィルターアルゴリズムを実装されてはいかがでしょうか。

もう一つ重要な点として、画像内のノイズの解析と分散の推定を行い非常に簡単なノイズモデルを作成しました。 これは画像の素性のわかっているカメラ画像処理だからできることであり、このような処理ができる事がRAW現像の利点であるともいえます。 今回は光ショットノイズが支配的な画像でしたが、他のノイズが大きい画像の場合はまた違った処理が必要です。 たとえば高ISO画像(低照度画像)ではリードノイズが相対的に大きくなってきます。

また、最後にnumpyによる画像処理の高速化テクニックについても触れました。 画像処理としては本質的ではありませんが、実用上は重要な点です。

今回の内容はすべてgithubにアップロードしてあります。

入力画像: raw_process/chart.jpg at master · moizumi99/raw_process · GitHub

Jupyter Notebook上のデモ: raw_process/Raspberry_Pi_Noise.ipynb at master · moizumi99/raw_process · GitHub

ノイズフィルター後の出力画像: raw_process/raspi_filtered_out.png at master · moizumi99/raw_process · GitHub

最後に

RAW画像から始まってカメラ画像処理の内容を実践してきましたが、このノイズフィルターで終わりもだいぶ近づいてきました。

もう少しで基本的なカメラ画像処理を一通りカバーする事ができます。 もう一息です。

次の記事

uzusayuu.hatenadiary.jp


  1. Bilateral Filtering for Gray and Color Images, C. Tomasi and R. Manduchi, Proceedings of the 1998 IEEE International Conference on Computer Vision, pp. 839–846, 1998

  2. ただし、処理の内容が一見わかりにくくなるというトレードオフがあります。これが理由で先程はベタ書きの処理を紹介しました。

  3. rawpyについて教えてくれたのと同じ同僚から教わりました。感謝します。

  4. Image and video denoising by sparse 3D transform-domain collaborative filtering | Block-matching and 3D filtering (BM3D) algorithm and its extensions

  5. https://web.stanford.edu/class/cs331b/2016/projects/zhao.pdf

ゼロから作るRAW現像 その6 - レンズシェーディング補正

このブログ記事「ゼロから作るRAW現像」を大きく再構成してより読みやすくした書籍「PythonとColabでできる-ゼロから作るRAW現像」を【技術書典6】にて頒布しました。

現在はBOOTHにて入手可能です。書籍+PDF版は2200円プラス送料、PDF版は1200円です。

moiz.booth.pm

はじめに

これは「ゼロから作るRAW現像 」という一連の記事の一つです。 これらの記事の内容を前提としていますので、まだお読みでない方はこちらの記事からお読みいただくことをおすすめします。

「ゼロから作るRAW現像 その1 - 基本的な処理」

「ゼロから作るRAW現像 その2 - 処理のモジュール化」

「ゼロから作るRAW現像 その3 - デモザイク処理基本編」

「ゼロから作るRAW現像 その4 - デモザイク処理応用編」

「ゼロから作るRAW現像 その5 - ラズベリーパイのRAW画像処理」

レンズシェーディング(周辺減光)

デジタルカメラに限らず、レンズを通して結像した画像は中央部より周辺部のほうが暗くなっています。

f:id:uzusayuu:20181021132917p:plain

このような現象をレンズシェーディングや単にシェーディング、また日本語では周辺減光などといいます。英語ではLens ShadingやVignettingと呼ばれます。

このような事が起きてしまう第一の原因は、レンズを通してセンサーにあたる光の量が、センサーの中央部と、センサーの周辺部とで異なる事です。 良く説明に出されるのが二枚のレンズを持つ単純な系の場合です。

f:id:uzusayuu:20181021090419p:plain

この図の例では、センサー中央部分にあたる光はレンズの開口部全体を通ってくるのに対して、センサー周辺部にあたる光はレンズの一部分しか通過できません。 結果的に中心部が明るく、周辺部が暗くなります。 実際のカメラの光学系はこれより遥かに複雑ですが、同様の理由により画像の周辺部分が暗くなります。

この要因の他に、センサー周辺部ではレンズやカバーの縁に遮られる光の量も増えます。また、光の入射角度によってレンズ等のコーティングと干渉する可能性もあります。 またデジタルカメラ独特の状況として、画像センサー上のフォトダイオードに到達する光の量が、光の入射角に依存するケースがあります。

シェーディングを起こす要因はこのように多岐にわたり、一概に、これが原因だとは言うことはできません。 したがってモデル計算で減光量を求めるよりも、実際のレンズで測定した結果を元に画像補正する必要があります。

レンズシェーディングの確認

ここで紹介する内容はGithubからJupyter Notebookファイルとしてダウンロードできます

では実際にシェーディングの影響を見てみましょう。

本来ならば明るさを均一にしたグレイチャート(その名の通り灰色の大きなシート)などを撮影してテストするのですが普通の家庭にそのような物はないので、今回はラズベリーパイのレンズの上に白いコピー用紙を載せ、そこに後ろから光を当てることでなるべく一様な明るさの画像を撮影しました。ファイル名は flat.jpgです。このファイルはgithubにアップロードしてあります

先日と同じ方法でRAW画像を取り出し現像してみましょう1。 まだレンズシェーディング補正は行いません。

%matplotlib inline
import numpy as np
import math, os, sys
from matplotlib import pyplot as plt 
module_path = os.path.abspath(os.path.join('..'))
if module_path not in sys.path:
    sys.path.append(module_path)
import raw_process4 as raw_process

with open("flat.jpg", "rb") as input_file:
    whole_data = input_file.read()
data = whole_data[-10237440:]

bayer_pattern = np.array([[2, 1], [1, 0]])
raw_w = 3282  # for v1.3 w = 2592
raw_h = 2480  # for v1.3 h = 1944
img = np.zeros((raw_h, raw_w))
stride = math.ceil(raw_w * 10 / 8 / 32) * 32
for y in range(raw_h):
    for x in range(raw_w // 4):
        word = data[y * stride + x * 5: y * stride + x * 5 + 5]
        img[y, 4 * x    ] = (word[0] << 2) | ((word[4] >> 6) & 3) 
        img[y, 4 * x + 1] = (word[1] << 2) | ((word[4] >> 4) & 3) 
        img[y, 4 * x + 2] = (word[2] << 2) | ((word[4] >> 2) & 3) 
        img[y, 4 * x + 3] = (word[3] << 2) | ((word[4]     ) & 3)

# Crop to the multiple of 32x32
w = 3264
h = 2464
img_clp = img[0:h, 8:8+w]
blacklevel = [64] * 4
blc_raw = raw_process.black_level_correction(img_clp, blacklevel)
wbg = np.array([1.105, 1, 2.609, 1])
wb_raw = raw_process.white_balance_Bayer(blc_raw, wbg, bayer_pattern)
dms_img = raw_process.advanced_demosaic(wb_raw, bayer_pattern)
ccm_matrix = (6022,-2314,394,-936,4728,310,300,-4324,8126)
img_ccm = raw_process.color_correction_matrix(dms_img, ccm_matrix)
img_ccm[img_ccm > 1023] = 1023
img_no_shading = raw_process.gamma_correction(img_ccm, 2.2)
outimg = img_no_shading.copy()
outimg = outimg.reshape((h, w, 3))
outimg[outimg < 0] = 0
plt.imshow(outimg)

画像を見てみましょう。

f:id:uzusayuu:20181021090327p:plain

確かに周辺光量が落ちているのが確認できます。

画像の中で明るさがどう変わっているか見てみましょう。 画面の高さ方向中央付近、上下32画素幅で左から右まで帯状の画像をとりだし、明るさがどう変わるをグラフにしてみます。

center_y, center_x = h // 2, w // 2
shading_profile = [[], [], []]
y = center_y - 16
for x in range(0, w - 32, 32):
    xx = x + 16
    shading_profile[0].append(img_no_shading[y:y+32, x:x+32, 0].mean())
    shading_profile[1].append(img_no_shading[y:y+32, x:x+32, 1].mean())
    shading_profile[2].append(img_no_shading[y:y+32, x:x+32, 2].mean())
shading_profile = [np.array(a) / max(a) for a in shading_profile]

plt.axis(ymin=0, ymax=1.1)
plt.plot(shading_profile[0], color='red')
plt.plot(shading_profile[1], color='green')
plt.plot(shading_profile[2], color='blue')
plt.show()

f:id:uzusayuu:20181021104028p:plain

画像の左右端では中心部分に比べて50%程度の明るさに落ちていることがわかります。 これはガンマ補正後の値なので、RAW画像では明るさの違いはさらに大きいことが予想できます。

また、赤画素と、緑・青画素とではシェーディングの様子がだいぶ違います。画像の色が中央と周辺部でだいぶ違うのはこのせいでしょう。

レンズシェーディングのモデル化

レンズシェーディングを補正するために、まずどの程度の減光があるのか測定してみましょう。

画像を見てわかるとおり、光の量は中心から離れるに従って減っています。中央からの距離によってパラメータ化できそうです。 また、対称性を考えると偶関数で近似できるはずです。

では各画素の明るさを、中央からの距離に応じてグラフにしてみましょう。

一画素毎に計算するのは計算量が大きくまたノイズによる誤差も入ってくるので、32 x 32のブロックごとに測定します。

vals = [[], [], [], []]
radials = []
index = 0
for y in range(0, h, 32):
    for x in range(0, w - 32, 32):
        xx = x + 16
        yy = y + 16
        r2 = (yy - center_y) * (yy - center_y) + (xx - center_x) * (xx - center_x)
        vals[0].append(blc_raw[y:y+32:2, x:x+32:2].mean())
        vals[1].append(blc_raw[y:y+32:2, x+1:x+32:2].mean())
        vals[2].append(blc_raw[y+1:y+32:2, x:x+32:2].mean())
        vals[3].append(blc_raw[y+1:y+32:2, x+1:x+32:2].mean())
        radials.append(r2)

これでvals[]には色ごとの画素の明るさ、radialsには中央からの距離の2乗が入っているはずです。

最大値でノーマライズしてグラフにして確認してみます。

rs = np.array(radials)
vs = np.array(vals)
norm = vs.max(axis=1)
vs[0, :] /= vs[0, :].max()
vs[1, :] /= vs[1, :].max()
vs[2, :] /= vs[2, :].max()
vs[3, :] /= vs[3, :].max()
plt.scatter(rs, vs[0,:], color='blue')
plt.scatter(rs, vs[1,:], color='green')
plt.scatter(rs, vs[2,:], color='green')
plt.scatter(rs, vs[3,:], color='red')
plt.show()

f:id:uzusayuu:20181021091622p:plain

きれいに中心からの距離に応じて明るさが減少しています。 また、この段階では周辺部では中心部に比べて3分の1程度に暗くなっていることがわかります。

これを補正するには、明るさの減少率の逆数をかけてやればよい事になります。 逆数のグラフを書いてみましょう。

gs = 1 / vs
plt.scatter(rs, gs[0,:], color='blue')
plt.scatter(rs, gs[1,:], color='green')
plt.scatter(rs, gs[2,:], color='green')
plt.scatter(rs, gs[3,:], color='red')
plt.show()

f:id:uzusayuu:20181021092058p:plain

これなら1次関数で近似できそうです2。やってみましょう。

par = [[], [], [], []]
for i in range(4):
    par[i] = np.polyfit(rs, gs[i, :], 1)

ここでpolyfit多項式近似を求める関数です。これでpar[]には近似関数の傾きと切片が入っているはずです。 確認してみましょう。

print(par)

[array([6.07106808e-07, 9.60556906e-01]), array([6.32044369e-07, 9.70694361e-01]), array([6.28455183e-07, 9.72493898e-01]), array([9.58743579e-07, 9.29427169e-01])]

それらしい値が入っています。グラフでみてみましょう。

es = [[], [], [], []]
for i in range(4):
    es[i] = par[i][0] * rs + par[i][1]
for i in range(4):
    plt.scatter(rs, es[i])
plt.show()

f:id:uzusayuu:20181021092820p:plain

良さそうです。

レンズシェーディング補正

ではいよいよ、実際の画像のレンズシェーディングを補正してみましょう。

まず、レンズシェーディング補正前の、ブラックレベル補正のみをかけたRAW画像がこちらです。 f:id:uzusayuu:20181021093347p:plain

先に各画素ごとに掛け合わせるゲインを、先程の近似関数から計算しておきます。

gain_map = np.zeros((h, w))
for y in range(0, h, 2):
    for x in range(0, w, 2):
        r2 = (y - center_y) * (y - center_y) + (x - center_x) * (x - center_x)
        gain = [par[i][0] * r2 + par[i][1] for i in range(4)]
        gain_map[y, x] = gain[0]
        gain_map[y, x+1] = gain[1]
        gain_map[y+1, x] = gain[2]
        gain_map[y+1, x+1] = gain[3]

このゲインをブラックレベル補正した画像にかけ合わせます。

lsc_raw = blc_raw * gain_map

補正後の画像がこちらです。

outimg = lsc_raw.copy()
outimg = outimg.reshape((h, w))
outimg = outimg / 1024
outimg[outimg < 0] = 0
outimg[outimg > 1] = 1
plt.imshow(outimg, cmap='gray')

f:id:uzusayuu:20181021093324p:plain

フラットな画像が出力されました!

残りの処理(ホワイトバランス補正、デモザイク、カラーマトリクス補正、ガンマ補正)を行ってフルカラー画像を出力してみましょう。

wb_raw = raw_process.white_balance_Bayer(lsc_raw, wbg, bayer_pattern)
dms_img = raw_process.advanced_demosaic(wb_raw, bayer_pattern)
img_ccm = raw_process.color_correction_matrix(dms_img, ccm_matrix)
img_ccm[img_ccm > 1023] = 1023
img_shading = raw_process.gamma_correction(img_ccm, 2.2)

表示します

outimg = img_shading.copy()
outimg = outimg.reshape((h, w, 3))
outimg[outimg < 0] = 0
plt.imshow(outimg)

f:id:uzusayuu:20181021093945p:plain

RGB画像で効果が確認できました。

残っているシェーディング量を測定してみましょう。

center_y, center_x = h // 2, w // 2
shading_after = [[], [], []]
y = center_y - 16
for x in range(0, w - 32, 32):
    xx = x + 16
    shading_after[0].append(img_shading[y:y+32, x:x+32, 0].mean())
    shading_after[1].append(img_shading[y:y+32, x:x+32, 1].mean())
    shading_after[2].append(img_shading[y:y+32, x:x+32, 2].mean())
shading_profile = [np.array(a) / max(a) for a in shading_after]

plt.axis(ymin=0, ymax=1.1)
plt.plot(shading_after[0], color='red')
plt.plot(shading_after[1], color='green')
plt.plot(shading_after[2], color='blue')
plt.show()

f:id:uzusayuu:20181021104206p:plain

シェーディングがほぼなくなりフラットになりました。 また、赤色と緑・青色との違いもほぼ消えました。成功のようです。

通常画像への適用

それではテスト用の平坦画像(Flat Field)ではなく、実際の画像にレンズシェーディング補正を適用してみましょう。

ラズベリーパイのカメラで改めて対象の画像(chart.jpg)をキャプチャします。

raspistill -r -o chart.jpg

f:id:uzusayuu:20181021122359p:plain:w300

以下の内容はGithubからJupyter Notebookファイルとしてダウンロードできます。

ここで使うテスト画像はGitHubにアプロードしてあります

ではRAW画像データを取り出し、まずはレンズシェーディング補正なしで現像してみます。3

%matplotlib inline
with open("chart.jpg", "rb") as input_file:
    whole_data = input_file.read()
data = whole_data[-10237440:]

import numpy as np
import math, os, sys
from matplotlib import pyplot as plt 
module_path = os.path.abspath(os.path.join('..'))
if module_path not in sys.path:
    sys.path.append(module_path)
import raw_process4 as raw_process

bayer_pattern = np.array([[2, 1], [1, 0]])
raw_w = 3282  # for v1.3 w = 2592
raw_h = 2480  # for v1.3 h = 1944
img = np.zeros((raw_h, raw_w))
stride = math.ceil(raw_w * 10 / 8 / 32) * 32
for y in range(raw_h):
    for x in range(raw_w // 4):
        word = data[y * stride + x * 5: y * stride + x * 5 + 5]
        img[y, 4 * x    ] = (word[0] << 2) | ((word[4] >> 6) & 3) 
        img[y, 4 * x + 1] = (word[1] << 2) | ((word[4] >> 4) & 3) 
        img[y, 4 * x + 2] = (word[2] << 2) | ((word[4] >> 2) & 3) 
        img[y, 4 * x + 3] = (word[3] << 2) | ((word[4]     ) & 3)

# Crop to the multiple of 32x32
w = 3264
h = 2464
img_clp = img[0:h, 8:8+w]

blacklevel = [64] * 4
blc_raw = raw_process.black_level_correction(img_clp, blacklevel)
wbg = np.array([1.128, 1, 2.546, 1])
wb_raw = raw_process.white_balance_Bayer(blc_raw, wbg, bayer_pattern)
dms_img = raw_process.advanced_demosaic(wb_raw, bayer_pattern)
ccm_matrix = (6022,-2314,394,-936,4728,310,300,-4324,8126)
img_ccm = raw_process.color_correction_matrix(dms_img, ccm_matrix)
img_ccm[img_ccm > 1023] = 1023
img_no_shading = raw_process.gamma_correction(img_ccm, 2.2)

表示してみましょう。

outimg = img_no_shading.copy()
outimg = outimg.reshape((h, w, 3))
outimg[outimg < 0] = 0
plt.imshow(outimg)

f:id:uzusayuu:20181021122939p:plain

悪くはないのですが、全体に青みがかっています。また、上記のJPEG画像と比べると周辺部が若干暗くなっているのがわかると思います。

それでは次にレンズシェーディング補正を入れて処理してみます。補正パラメータは先程の平坦画像で計算したものを使います。

par = [np.array([6.07106808e-07, 9.60556906e-01]), 
       np.array([6.32044369e-07, 9.70694361e-01]), 
       np.array([6.28455183e-07, 9.72493898e-01]), 
       np.array([9.58743579e-07, 9.29427169e-01])]

gain_map = np.zeros((h, w))
center_y, center_x = h // 2, w // 2
for y in range(0, h, 2):
    for x in range(0, w, 2):
        r2 = (y - center_y) * (y - center_y) + (x - center_x) * (x - center_x)
        gain = [par[i][0] * r2 + par[i][1] for i in range(4)]
        gain_map[y, x] = gain[0]
        gain_map[y, x+1] = gain[1]
        gain_map[y+1, x] = gain[2]
        gain_map[y+1, x+1] = gain[3]

lsc_raw = blc_raw * gain_map
wb_raw = raw_process.white_balance_Bayer(lsc_raw, wbg, bayer_pattern)
dms_img = raw_process.advanced_demosaic(wb_raw, bayer_pattern)
img_ccm = raw_process.color_correction_matrix(dms_img, ccm_matrix)
img_ccm[img_ccm > 1023] = 1023
img_shading = raw_process.gamma_correction(img_ccm, 2.2)

見てみましょう。

outimg = img_shading.copy()
outimg = outimg.reshape((h, w, 3))
outimg[outimg < 0] = 0
plt.imshow(outimg)

f:id:uzusayuu:20181021123223p:plain

先程の画像に比べると明るさが均一になり、不自然な青みも消えました。 成功のようです。

セーブしておきましょう。

raw_process.write(img_shading, "chart_with_shading_correction.png")

まとめ

今回はレンズシェーディング補正(周辺減光補正)をとりあげました。 おそらく、『カメラ』画像処理以外のいわゆる画像処理では取り上げることのない特殊な処理だと思います。

今回行ったのは、半径方向の二次多項式による補正ゲインの近似で、レンズシェーディング補正の中では最も単純なものです。 実際のカメラの中では、より高次の関数による近似や、2次元ルックアップテーブルによる補正などが行われているのが普通です。 また、補正パラメータも、明るさや光源の種類、オートフォーカスの場合はフォーカス位置、などにより調整します。

一見単純そうな見た目や効果と比べて、実際には遥かに複雑で非常に重要な処理です。ある意味カメラの出力画像の画質を決める肝と言ってもよいと思います。

今回の内容はテストデータと共に、githubにアップロードしてあります。

次の記事

uzusayuu.hatenadiary.jp


  1. ここで使ったホワイトバランスとカラーマトリクスはexiftoolを使って exiftool -EXIF:MakerNoteUnknownText -b flat.jpgで確認できます。

  2. 横軸は距離の二乗なので、距離の関数としては二次多項式になります。

  3. ホワイトバランスゲインとカラーマトリクスはテスト画像と同様に exiftool -EXIF:MakerNoteUnknownText chart.jpg で確認できる。

ゼロから作るRAW現像 その5 - ラズベリーパイのRAW画像処理

このブログ記事「ゼロから作るRAW現像」を大きく再構成してより読みやすくした書籍「PythonとColabでできる-ゼロから作るRAW現像」を【技術書典6】にて頒布しました。

現在はBOOTHにて入手可能です。書籍+PDF版は2200円プラス送料、PDF版は1200円です。

moiz.booth.pm

はじめに

これは「ゼロから作るRAW現像 」という一連の記事の一つです。 これらの記事の内容を前提としていますので、まだお読みでない方はこちらの記事からお読みいただくことをおすすめします。

「ゼロから作るRAW現像 その1 - 基本的な処理」

「ゼロから作るRAW現像 その2 - 処理のモジュール化」

「ゼロから作るRAW現像 その3 - デモザイク処理基本編」

「ゼロから作るRAW現像 その4 - デモザイク処理応用編」

ラズベリーパイのRAW画像

これまでRAW現像の対象として、ソニーα7iiiで撮影した画像を使ってきました。 これは基本的な処理の説明をするにあたって、ノイズが少ない、歪がない、解像感が高い、など質の良いRAW画像のほうがやりやすかったからです。 しかし、実際のカメラ画像処理では、ノイズや意図しない不鮮明さなど好ましくない性質を修正するというのが大きな目的になります。そういった処理を行うにあたって、ミラーレスカメラの高画質データは逆に扱いにくいものになってしまいます。なにしろ直したいノイズもボケもあまりありませんから。

今回はラズベリーパイでRAW画像をキャプチャーして、これまで行ったRAW現像処理を行い、次回以降のテーマの準備を行います。

ラズベリーパイの準備とRAW画像のキャプチャー

今回紹介する内容について前提とする環境は以下のとおりです。

他のラズベリーパイやv1.3カメラでも実行は可能だとは思われますが、試してはいません。

  • 画像処理PCサイド
    • Python3が動作し前回までの内容が実行できるLinux環境
    • SDカードリーダーやラズベリーパイとのネットワーク接続など、ラズベリーパイとデータを交換する手段

以下の内容はラズベリーイカメラv2.1がラズベリーパイに接続され正常に動作していることを前提としていますので、まず公式ドキュメント などに従ってカメラの動作を確認下さい1

特に、GUIの左上のラズベリーパイアイコンから選択できる「設定」、「RaspberryPiの設定」、から「インターフェイス」タブで「カメラ」を有効にしておく必要がありますので、ご注意ください。

RAW画像のキャプチャー

ラズベリーパイにカメラをセッティングした状態で、以下のコマンドを実行します。

> raspistill -r -o raw_capture.jpg

-rがRAWキャプチャーを行うことをしめしています。

うまくいけば、カレントディレクトリにraw_capture.jpgという名前のJPEGファイルができているはずです。 このファイルは一見通常のJPEGファイルに見えますが、8MPのJPG画像ファイルにしては15MB前後と巨大なファイルサイズです。これは、RAWデータが埋め込みデータとしてファイル中に組み込まれているからです。

キャプチャーが成功したら、このファイルをUSBメモリーで移す、ネットワークコピー、またはRaspbianの入ったSDCARDをホストPCに読み込ませる、などの方法で、前回までの作業を行ったPCに移動します2

うまくキャプチャーできたか、JPEGファイルのJPEG画像部分を見てみましょう。通常の画像ビューワーや画像エディタで確認できます。

f:id:uzusayuu:20181014122747p:plain

あまり良い画像ではないですが、キャプチャーできていることは確認できました。

ラズベリーパイのRAW画像の抽出

以下の内容はgithubからJupyter Notebookファイルとしてダウンロードできます

今回の画像はこれまでと違い、JPEGファイルの中に埋め込まれているので、rawpy以外の方法で取り出す必要があります。 なお、この部分は今回の本筋ではないので、簡単な説明ですまします。詳しい内容は picameraのドキュメント を参照ください。

これ以降は、Jupyterなどpython3のインタラクティブ環境で作業します。

BayerデータはJPEGファイルの最後の部分に位置するので、読み出してから末尾だけ取り出します。

with open("raw_capture.jpg", "rb") as input_file:
    data = input_file.read()
data = data[-10237440:]

(未確認ですがカメラv1.3ではdata = data[-6371328:]とするとよいようです。)

これでBayer部分が読み込めたはずです。データを見てみましょう。

import numpy as np

w = 3282  # for v1.3 w = 2592
h = 2480  # for v1.3 h = 1944
img = np.zeros((h, w))
stride = math.ceil(w * 10 / 8 / 32) * 32
for y in range(h):
    for x in range(w // 4):
        word = data[y * stride + x * 5: y * stride + x * 5 + 5]
        img[y, 4 * x    ] = (word[0] << 2) | ((word[4] >> 6) & 3) 
        img[y, 4 * x + 1] = (word[1] << 2) | ((word[4] >> 4) & 3) 
        img[y, 4 * x + 2] = (word[2] << 2) | ((word[4] >> 2) & 3) 
        img[y, 4 * x + 3] = (word[3] << 2) | ((word[4]     ) & 3)

最後の部分ですが、picameraのドキュメントによると、ラズベリーパイのBayerデータは4画素毎に5バイトのデータにまとまっていて、最初の4バイトがそれぞれの画像の上位8ビット、最後のバイトの2ビットずつが、各画素の下位2ビットになっているとのことです。

さて、ちゃんとデータが読めたのか見てみましょう。

from matplotlib.pyplot import imshow 

outimg = img.copy()
outimg[outimg < 0] = 0
outimg = outimg / outimg.max()
imshow(outimg, cmap='gray')

f:id:uzusayuu:20181014124822p:plain

どうやらそれらしい画像が取り出せたようです。

ラズベリーイカメラのRAW画像の現像

だいぶ手間がかかりましたが、ここから本題のラズベリーパイのRAW画像の現像に入ります。

まず前回までのスクリプトを使いたいところですが、これまで使った関数は殆どがrawpyのインスタンスに依存していました。 今回扱うラズベリーパイのRAW画像はrawpyで読み込んだものではないのでそのままでは使えません。

そこで各関数を、rawデータと指定した任意のパラメータで実行できるように書き換えて、raw_process4.pyとしました。 それぞれの関数の変更点は以下の説明で触れます。

まずはモジュールを読み込みましょう。

import os
import sys
module_path = os.path.abspath(os.path.join('..'))
if module_path not in sys.path:
    sys.path.append(module_path)
import raw_process4 as raw_process

ブラックレベル補正

ブラックレベル補正は、引数としてブラックレベルの数値を与えるようにしました。 引数で与えるブラックレベルは、左上、右上、左下、右下、の順にブラックレベルを格納したリストです。

def black_level_correction(raw_array, black_level):
    blc_raw = raw_array.copy()
    blc_raw[0::2, 0::2] -= black_level[0]
    blc_raw[0::2, 1::2] -= black_level[1]
    blc_raw[1::2, 0::2] -= black_level[2]
    blc_raw[1::2, 1::2] -= black_level[3]
    return blc_raw

ブラックレベルを64と仮定して処理してみましょう。

blacklevel = [64] * 4
blc_raw = raw_process.black_level_correction(img, blacklevel)

f:id:uzusayuu:20181014152120p:plain

ホワイトバランス補正

def white_balance_Bayer(raw_array, wbg, bayer_pattern):
    img_wb = raw_array.copy()
    img_wb[0::2, 0::2] *= wbg[bayer_pattern[0, 0]] 
    img_wb[0::2, 1::2] *= wbg[bayer_pattern[0, 1]]
    img_wb[1::2, 0::2] *= wbg[bayer_pattern[1, 0]]
    img_wb[1::2, 1::2] *= wbg[bayer_pattern[1, 1]]
    return img_wb

ホワイトバランス補正には、R・G・Bのゲインを並べたリストのwbgと、ベイヤーのパターンを渡すようになりました。

今回ホワイトバランスははっきりしないので仮に赤のゲインx1.5、青のゲインx2.2を与えます。 今回のラズベリーパイの画像では左上が青色画素ですので、bayer_patternは[[2, 1], [1, 0]]になります。

wbg = np.array([1.5, 1, 2.2, 1])
bayer_pattern = np.array([[2, 1], [1, 0]])
wb_raw = raw_process.white_balance_Bayer(blc_raw, wbg, bayer_pattern)

f:id:uzusayuu:20181014145358p:plain

デモザイク

前回のデモザイクは画像の左上隅が赤色画素であることを仮定していました。 今回のラズベリーパイの画像では左上が青色画素ですのでこのままでは実行できません。 ベイヤーパターンごとの違いは各チャンネルの位相を180度変える事になります。

def advanced_demosaic(dms_input, bayer_pattern):
    hlpf = np.array([[1, 2, 3, 4, 3, 2, 1]]) / 16
    vlpf = np.transpose(hlpf)
    hhpf = np.array([[-1, 2, -3, 4, -3, 2, -1]]) / 16
    vhpf = np.transpose(hhpf)
    identity_filter = np.zeros((7, 7))
    identity_filter[3, 3] = 1

    # generate FIR filters to extract necessary components
    FC1 = np.matmul(vhpf, hhpf)
    FC2H = np.matmul(vlpf, hhpf)
    FC2V = np.matmul(vhpf, hlpf)
    FL = identity_filter - FC1 - FC2V - FC2H

    # f_C1 at 4 corners
    c1_mod = signal.convolve2d(dms_input, FC1, boundary='symm', mode='same')
    # f_C1^1 at wy = 0, wx = +Pi/-Pi
    c2h_mod = signal.convolve2d(dms_input, FC2H, boundary='symm', mode='same')
    # f_C1^1 at wy = +Pi/-Pi, wx = 0
    c2v_mod = signal.convolve2d(dms_input, FC2V, boundary='symm', mode='same')
    # f_L at center
    f_L = signal.convolve2d(dms_input, FL, boundary='symm', mode='same')

    # Move c1 to the center by shifting by Pi in both x and y direction
    # f_c1 = c1 * (-1)^x * (-1)^y
    f_c1 = c1_mod.copy()
    f_c1[:, 1::2] *= -1
    f_c1[1::2, :] *= -1
    if bayer_pattern[0, 0] == 1 or bayer_pattern[0, 0] == 3:
        f_c1 *= -1
    # Move c2a to the center by shifting by Pi in x direction, same for c2b in y direction
    c2h = c2h_mod.copy()
    c2h[:, 1::2] *= -1
    if bayer_pattern[0, 0] == 2 or bayer_pattern[1, 0] == 2:
        c2h *= -1
    c2v = c2v_mod.copy()
    c2v[1::2, :] *= -1
    if bayer_pattern[0, 0] == 2 or bayer_pattern[0, 1] == 2:
        c2v *= -1
    # f_c2 = (c2v_mod * x_mod + c2h_mod * y_mod) / 2
    f_c2 = (c2v + c2h) / 2

    # generate RGB channel using 
    # [R, G, B] = [[1, 1, 2], [1, -1, 0], [1, 1, - 2]] x [L, C1, C2]
    height, width = dms_input.shape
    dms_img = np.zeros((height, width, 3))
    dms_img[:, :, 0] = f_L + f_c1 + 2 * f_c2
    dms_img[:, :, 1] = f_L - f_c1
    dms_img[:, :, 2] = f_L + f_c1 - 2 * f_c2

    return dms_img

実行してみます。

dms_img = raw_process.advanced_demosaic(wb_raw, bayer_pattern)

f:id:uzusayuu:20181014145414p:plain

カラーマトリクス補正

カラーマトリクスの値もはっきりしません。今回は色が多少強調されるように、以下のようなマトリクスをかけ合わせてみます。

[[1536, -256, -256], 
 [-256, 1536, -256],
 [-256, -256, 1536]]

実行します。

img_ccm = raw_process.color_correction_matrix(dms_img, [1536, -256, -256, -256, 1536, -256, -256, -256, 1536])

f:id:uzusayuu:20181014145426p:plain

ガンマ補正

最後にガンマ補正です。ガンマ値を引数でもたせるようにしました。

def gamma_correction(rgb_array, gamma):
    img_gamma = rgb_array.copy()
    img_gamma[img_gamma < 0] = 0
    img_gamma = img_gamma / img_gamma.max()
    img_gamma = np.power(img_gamma, 1/gamma)
    return img_gamm

実行します。

img_gamma = raw_process.gamma_correction(img_ccm, 2.2)

画像を保存して確認してみましょう。

raw_process.write(img_gamma, "raspi_raw_out.png")

import imageio
from pylab import imshow, show
imshow(imageio.imread('raspi_raw_out.png'))
show()

f:id:uzusayuu:20181014140022p:plain

お世辞にもきれいな画像とは言えませんがどうにかラズベリーパイでキャプチャーしたRAWデータを現像して画像ファイルにすることができました。

まとめ

Raspberry Pi 3BでキャプチャーしたRAW画像を現像して画像ファイルに変換しました。 今回はどうにか現像するところで終わりでしたが、次回以降画質に手をいれて行きたいと思います。

今回の内容はRAW画像データ及びraw_process4.pyと共にgithubにアップロードしてあります。

github.com

次の記事

uzusayuu.hatenadiary.jp


  1. 公式ドキュメント以外では以下の記事が参考になります。Raspberry PiカメラモジュールV2でデジカメを作ってみた

  2. 前回までの環境をRaspberry Pi上で構築できていれば同じことができるはずですが、私の環境ではうまくいきませんでしたので、別のPCにデータを運んでから以下の作業を行っています。

ゼロから作るRAW現像 その4 - デモザイク処理応用編

このブログ記事「ゼロから作るRAW現像」を大きく再構成してより読みやすくした書籍「PythonとColabでできる-ゼロから作るRAW現像」を【技術書典6】にて頒布しました。

現在はBOOTHにて入手可能です。書籍+PDF版は2200円プラス送料、PDF版は1200円です。

moiz.booth.pm

はじめに

これは「ゼロから作るRAW現像 」という一連の記事の一つです。 これらの記事の内容を前提としていますので、まだお読みでない方はこちらの記事からお読みいただくことをおすすめします。

モジュールの高速化

本題に入る前に、raw_process.pyの高速化についてです。 これまでは処理の内容がわかりやすいように、画素単位で処理を書き下してきましたが、この方法ではnumpyの内部の高速性を十分に活かすことができません。

前前回までは縮小デモザイクを使っていたために画像サイズが小さく、処理速度は問題化していなかったのですが、前回縮小しないデモザイクに置き換えたために全体の処理が重くなってしまいました。

そこで、raw_process2.pyの処理をなるべくループ処理を使わず書き換えることでnumpyの恩恵を十分に使って高速化できるように書き換えて、raw_process3.pyとしました。 今回以降はこのモジュールを使用します。

書き換えの内容については、画像処理的にはまったく同じで、ループの順番を替えただけなので、ここでは解説は割愛します。興味のある方はこちらをご覧ください。

raw_process/raw_process3.py at master · moizumi99/raw_process · GitHub

前回のデモザイク処理の分析

この節は、画像処理の理論的背景に興味がない方は読み飛ばして構いません。

それではまず、前回のデモザイク処理の問題点を分析してみましょう。

前回の記事の最後で、線形補間処理をFIRフィルターとして書き換える作業を行いました。これは処理の高速化を念頭においたものですが、FIRを使った画像処理として解釈することもできます。

まず、ベイヤー配列のある特定の色の画素は、もともとのフル解像度の画像に、位置によって0または1を取るサンプリング関数を画像にかけたものと解釈する事ができます。

f:id:uzusayuu:20181007145847p:plain

これを数式で表すとこうなります。

{ \displaystyle f_{Green}(x, y) = f_{full}(x, y) \cdotp s(x, y)}

{ \displaystyle
s(x, y) =
\begin{cases}
0   \text{  if } x \mod 2 = y \mod 2 \\
1  \text{  otherwise}
\end{cases}
}

 s(x, y) はさらに

\begin{equation} s_{Green} \left ( x, y \right ) = \frac{1}{2} \left ( 1 - \left ( -1 \right ) ^ {x + y} \right ) = \frac{1}{2} \left ( 1 - e ^ {i \pi \left ( x + y \right ) } \right ) \end{equation}

となり、したがって、元の式は

\begin{equation} f_{Green} \left ( x, y \right ) = \frac{1}{2} f_{full}\left ( 1 - \left ( -1 \right ) ^ {x + y} \right ) = \frac{1}{2} \left ( f_{full} - f_{full} e ^ {i \pi \left ( x + y \right ) } \right ) \end{equation}

と表されます。この式の後半は、元の関数をx軸方向とy軸方向に \piずつずらし、元の信号から引く事を意味します。

つまり、もし元の画像の周波数特性が模式的にこのように表されるとしたら、

f:id:uzusayuu:20181007163828p:plain

Bayer上の緑画素の周波数特性はこのようになります。

f:id:uzusayuu:20181007163915p:plain

青い部分が、全式の後半部分にできたもので、これはサンプリングによるエイリアシングを表しています。

画像の周波数成分は低い部分に多いと考えると、サンプリング前の信号を取り出すには、低周波部分を取り出してやれば良いことになります。 信号処理的にはローパスフィルターをかけます。

ここで、線形補間の式から作ったFIRをフィルターを思い出すと、このような形をしていました。

[[   0, 1/4,   0],
 [ 1/4,   1, 1/4],
 [   0, 1/4,   0]]

これはそのものズバリ、ローパスフィルターです。周波数特性はこんな感じになります。

f:id:uzusayuu:20181007164819p:plain

同様に、赤画素のサンプリング関数は、

\begin{equation} s_{Red} \left ( x, y \right ) = \frac{1}{4} \left ( 1 + \left ( -1 \right ) ^ x \right ) \left ( 1 + \left ( -1 \right ) ^ y \right ) = \frac{1}{4} \left ( 1 + e ^ {i \pi x} \right ) \left ( 1 + e ^ {i \pi y} \right )
\end{equation}

で、青画素は、

\begin{equation} s_{Blue} \left ( x, y \right ) = \frac{1}{4} \left ( 1 - \left ( -1 \right ) ^ x \right ) \left ( 1 - \left ( -1 \right ) ^ y \right ) = \frac{1}{4} \left ( 1 - e ^ {i \pi x} \right ) \left ( 1 - e ^ {i \pi y} \right )
\end{equation}

エイリアシングの様子はこうなります。

f:id:uzusayuu:20181007164303p:plain

緑の画素にくらべて、エイリアシングの影響を受ける領域が大きくなっています。したがって、緑よりも強めのローパスフィルターをかける必要があることが予想されます。

前回使用した赤・青画素線形補間用ののFIRフィルターはこうなっています。

[[1/4, 1/2, 1/4],
 [1/2,   1, 1/4].
 [1/4, 1/2, 1/4]]

周波数特性はこうです。

f:id:uzusayuu:20181007164836p:plain

ちょうとエイリアシングを取り除くような特性になっています。

こうしてみると理想的なフィルター処理をしているように見えますが、実際の画像はだいぶぼやけていて、偽色という本来ないはずの色がでてきています。 画像のぼやけは、もともとの画像の周波数特性を十分にカバーしていないことを、偽色の存在は本来の信号とエイリアシングとが十分に分割できていない事を示しています。 このような単純な形のフィルターの限界でしょう。

より高度なデモザイク処理

この節も画像処理の理論的バックグラウンドに興味なければ読み飛ばして構いません。

この節の内容の大部分は"Frequency-Domain Methods for Demosaicking of Bayer-Sampled Color Images, Eric Dubois, IEEE Signal Processing Letters, Dec. 2005"に基づいています。

線形補間デモザイクでは各色チャンネルをそれぞれ独立なものとして処理しました。ここでは、3色のBayer配列全体を一つの画像と見てみます。

\begin{aligned} f_{bayer} = & \frac{1}{2} f_G \left ( 1 - \left ( -1 \right ) ^ {x + y} \right ) + \\ & \frac{1}{2} f_R \left ( 1 + \left ( -1 \right ) ^ {x} \right ) \left ( 1 + \left ( -1 \right ) ^ {y} \right ) + \\ & \frac{1}{2} f_B \left ( 1 - \left ( -1 \right ) ^ {x} \right ) \left ( 1 - \left ( -1 \right ) ^ {y} \right ) \end{aligned}

ここで、 \begin{equation} \left ( -1 \right ) ^ {x} = e ^ {i \pi x} \end{equation}

などを利用して、

\begin{aligned} f_{bayer} = & \left ( \frac{1}{4} f_R + \frac{1}{2} f_G + \frac{1}{4} f_B \right ) + \\ & \left ( \frac{1}{4} f_R - \frac{1}{2} f_G + \frac{1}{4} f_B \right ) \cdot e^{i \pi (x + y)} + \\ & \left ( \frac{1}{4} f_R - \frac{1}{4} f_B \right ) \cdot (e^{i \pi x} + e^{i \pi y}) \end{aligned}

まとめると

\begin{equation} \end{equation}

\begin{aligned} f_{bayer} & = f_{L} - f_{C1} e^{i \pi \left ( x + y \right ) } + f_{C2} \left ( e^{i \pi x} + e^{i \pi y} \right ) \\ f_L & = \frac{1}{4} f_R + \frac{1}{2} f_G + \frac{1}{4} f_B \\ f_{C1} & = \left ( \frac{1}{4} f_R - \frac{1}{2} f_G + \frac{1}{4} f_B \right ) \\ f_{C2} & = \left ( \frac{1}{4} f_R - \frac{1}{4} f_B \right ) \end{aligned}

逆に、 f_L f_{C1} f_{C2}から、 f_R f_G f_Bは、

\begin{cases} f_R = f_L + f_{C1} + 2f_{C2} \\ f_G = f_L - f_{C1} \\ f_B = f_L + f_{C1} - 2f_{C2} \end{cases}

と、求められます。

では、 f_L f_{C1} f_{C2}はどうやって計算すれば良いのか、というのが次の問題になります。

 f_L f_{C1} f_{C2}の定義をみると、 f_R f_G f_Bの線形結合なので、元の画像が周波数制限されていて低い周波数の成分しかもっていないとすると、 f_L f_{C1} f_{C2}も周波数制限されていると考えられます。

ここで、 f_{C1} e^{i \pi (x + y)} f_{C2} \left ( e^{i \pi x} + e^{i \pi y} \right )は、 f_{C1} f_{C2}を周波数空間で \piだけシフトしたものです。 もとのf_{Bayer}はこれらを足し合わせたものなので、図に表すとこうなります。

f:id:uzusayuu:20181008004120p:plain

ということは、f_{Bayer}に帯域制限フィルターをかけて各成分をとりだし、周波数をシフトすれば f_L f_{C1} f_{C2}を得ることがきるという事になります。

必要なフィルターは4つです。

  • FC1 : 縦・横ハイパスフィルター, f_{C1}を取り出すのに使います。

例: \begin{equation} F_{C1} = \left(\begin{array}{c} -1 \\ 2 \\ -3 \\ 4 \\ -3 \\ 2 \\ -1\end{array}\right) \left(\begin{array}{c} -1, 2, -3, 4, -3, 2, -1\end{array}\right) \end{equation}

  • FC2V : 縦方向ハイパスフィルター、横方向ローパスフィルター、f_{C2}の上下の成分を取り出すのに使います。

例: \begin{equation} F_{C2V} = \left(\begin{array}{c} -1 \\ 2 \\ -3 \\ 4 \\ -3 \\ 2 \\ -1\end{array}\right) \left(\begin{array}{c} 1, 2, 3, 4, 3, 2, 1\end{array}\right) \end{equation}

  • FC2H: 縦方向ローパスフィルター、横方向ハイパスフィルター、f_{C2}の左右の成分を取り出すのに使います。

例: \begin{equation} F_{C2H} = \left(\begin{array}{c} 1 \\ 2 \\ 3 \\ 4 \\ 3 \\ 2 \\ 1\end{array}\right) \left(\begin{array}{c} -1, 2, -3, 4, -3, 2, -1\end{array}\right) \end{equation}

  • FL: (FC1 + FC2V + FC2H)を反転させたフィルター、f_{L}の成分を取り出すのに使います。

例 \begin{equation} F_L = 1 - F_{C1} - F_{C2V} - F_{C2H} \end{equation}

これらを使えば、

\begin{aligned} f_{C1} & = F_{C1} \left ( f_{Bayer} \right ) e^{- i \pi \left ( x + y \right )} \\ f_{C2V} & = F_{C2V} \left ( f_{Bayer} \right ) e^{-i \pi y} \\ f_{C2H} & = F_{C2H} \left ( f_{Bayer} \right ) e^{-i \pi x} \\ f_{C2} & = \left ( f_{C2V} + f_{C2H} \right ) / 2 \\ f_L & = F_L \left ( f_{Bayer} \right ) \end{aligned}

となって、 f_L f_{C1} f_{C2}の近似が計算でき、そこから、 f_R f_G f_Bが求められるはずです。 やってみましょう。

より高度なデモザイク処理、実行

まずは、必要なライブラリー郡を呼び出します。

%matplotlib inline
import numpy as np
import os
import sys
module_path = os.path.abspath(os.path.join('..'))
if module_path not in sys.path:
    sys.path.append(module_path)
import imageio
from pylab import imshow, show
import scipy
import numpy as np
import raw_process3 as raw_process

次に画像を読み込み、ホワイトバランスまでの処理を行います。

raw = raw_process.read("sample2.ARW")
color_matrix = [1141, -205, 88, -52, 1229, -154, 70, -225, 1179]
raw_array = raw_process.get_raw_array(raw)
blc_raw = raw_process.black_level_correction(raw, raw_array)
wb_raw = raw_process.white_balance_Bayer(raw, blc_raw)
dms_input = wb_raw

なお、ホワイトバランスをデモザイク先に行ったのは、前回までの各色をバラバラに処理を行うデモザイクアルゴリズムと違い、今回は3色を同時に取り扱う処理が入るため、事前に色バランスを合わせる必要があるからです。

次に必要なFIRフィルターを準備しましょう。 まず、横方向のローパスフィルターです。

hlpf = np.array([[1, 2, 3, 4, 3, 2, 1]]) / 16

これを利用して、縦方向のローパスフィルターを作ります。

vlpf = np.transpose(hlpf)

同様に、横・縦方向のハイパスフィルターを作ります。

hhpf = np.array([[-1, 2, -3, 4, -3, 2, -1]]) / 16
vhpf = np.transpose(hhpf)

次に、等価フィルターを作ります。

identity_filter = np.zeros((7, 7))
identity_filter[3, 3] = 1

これらの基本的なフィルターから、必要なフィルターを作っていきます。

FC1は、縦横ハイパスフィルターなので、こうなります。

FC1 = np.matmul(vhpf, hhpf)

FC2Vは縦ハイパス・横ローパス、FC2Hは縦ローパス・横ハイパスです。

FC2V = np.matmul(vhpf, hlpf)
FC2H = np.matmul(vlpf, hhpf)

FLは、等価フィルターからFC1, FC2V, FC2Hを取り除いたものです。

FL = identity_filter - FC1 - FC2V - FC2H

これで準備ができました。

ではまず、ハイパスフィルターを使って、C1成分を周波数空間の角からとりだしてみましょう。 SciPyの2次元コンボリューション機能を使います。

# f_C1 at 4 corners
c1_mod = scipy.signal.convolve2d(dms_input, FC1, boundary='symm', mode='same')

同様に、C2V成分、C2H成分を取り出してみましょう。

c2v_mod = scipy.signal.convolve2d(dms_input, FC2V, boundary='symm', mode='same')
c2h_mod = scipy.signal.convolve2d(dms_input, FC2H, boundary='symm', mode='same')

最後にL成分をとりだします。

f_L = scipy.signal.convolve2d(dms_input, FL, boundary='symm', mode='same')

今の所、C1成分、C2V・C2H成分は、周波数空間で \piだけずれています。色情報を復元するにはこれを復調してやる必要があります。

復調するためには、e^{-i \pi x}e^{-i \pi y}を掛け合わせなくてはならないのですが、x, yはいずれも整数なので、結局、-11の繰り返しになります。 つまり、e^{-i \pi x}をかけるのは奇数列の成分に-1をかける、e^{-i \pi y}をかけるのは奇数行の成分に-1をかけるのと同じです。

まずC1についてやってみましょう。必要な計算は、

\begin{equation} f_{c1} = f_{c1mod} \left ( x , y \right ) \left ( -1 \right ) ^ x \left ( -1 \right ) ^ y \end{equation}

ですので、

f_c1 = c1_mod.copy()
f_c1[:, 1::2] *= -1
f_c1[1::2, :] *= -1

となります。

同様にC2V・C2Hでは必要な計算は、

\begin{equation} f_{cv} = f_{c2vmod} \left ( x , y \right ) \left ( -1 \right ) ^ y \\ f_{ch} = f_{c2hmod} \left ( x , y \right ) \left ( -1 \right ) ^ x \end{equation}

で、プログラムは以下のようになります。

c2v = c2v_mod.copy()
c2v[1::2, :] *= -1
c2h = c2h_mod.copy()
c2h[:, 1::2] *= -1
f_c2 = (c2v + c2h) / 2

これで f_L f_{C1} f_{C2}が求められました。

最後に、以下の式にしたがって f_R f_G f_Bを計算します。

\begin{cases} f_R = f_L + f_{C1} + 2f_{C2} \\ f_G = f_L - f_{C1} \\ f_B = f_L + f_{C1} - 2f_{C2} \end{cases}

やってみましょう。

height, width = dms_input.shape
dms_img = np.zeros((height, width, 3))
dms_img[:,:,0] = f_L + f_c1 + 2 * f_c2
dms_img[:,:,1] = f_L - f_c1
dms_img[:,:,2] = f_L + f_c1 - 2 * f_c2

これでRGB画像が再現されたはずです! 残りの処理を行い結果を確認してみましょう。

img_ccm = raw_process.color_correction_matrix(dms_img, color_matrix)
rgb_image = raw_process.gamma_correction(img_ccm)

表示してみます。

outimg = rgb_image.copy()
outimg[outimg < 0] = 0
outimg = outimg / outimg.max()
imshow(outimg)

f:id:uzusayuu:20181008022550p:plain

きれいに出ているようにみえます。細部は後ほど確認するとして、まずは画像をセーブしておきましょう。

raw_process.write(rgb_image, "advanced_demosaic.png")

画像の確認

では前回の線形補間による画像と画質を比べてみましょう。

f:id:uzusayuu:20181008024743p:plain:w800

左が前回の線形補完法によるデモザイク、右が今回のデモザイクです。全体にシャープさがあがっています。

別の部分の画像です。前回気になった偽色がだいぶ減っています。

f:id:uzusayuu:20181008025157p:plain:w800

シャープさの向上、偽色の低減から、今回の処理は成功のようです。

まとめ

今回は、若干複雑なデモザイク処理を画像処理のレベルから解説し、Pythonに実装しました。 ここで実行した処理は、AdvancedDemosaickingDemo.ipynb としてgithubにアップロードしてあります。

また、今回の処理を新たにdemosaic()という名前のメソッドとしてraw_process3.pyに追加し、デフォルトのデモザイクアルゴリズムを置き換えました。 以下の方法で今回のデモザイクアルゴリズムを使ったRAW現像が可能です。

> python3 raw_process3.py sample2.ARW test_out.png "1141, -205, 88, -52, 1229, -154, 70, -225, 1179"

参照した論文1ではさらに、一部のチャンネルを混ぜ合わせるときに各画素毎に重み計算を行うという適合処理的アルゴリズムを扱っています。今回は割愛しましたが、興味のある方は試されてはいかがでしょうか2

最後に

今回は応用編と銘打ってはいますが、デモザイクアルゴリズムとしてはこれでもかなりシンプルな処理の部類です。 また、このアルゴリズムが発表されたのがすでに13年前であり、これまでにさらにさまざまなアルゴリズムが提案されており、現在も進歩をつづけています3。 たとえば、最近の流行はデモザイク処理にディープラーニングを使う手法です4

こういった事情を考慮にいれても、今回紹介した手法はシンプルな割には高画質の出力が得られる良いアルゴリズムだと思います。私の個人的なお気に入りのアルゴリズムの一つです。

今回の内容はgithubにて公開してあります。

github.com

次の記事

uzusayuu.hatenadiary.jp


  1. Frequency-Domain Methods for Demosaicking of Bayer-Sampled Color Images, Eric Dubois, IEEE Signal Processing Letters, Dec. 2005

  2. 論文の著者グループがWebにてMatlabのサンプルプログラムを配布しています。また、今回のスクリプトを少し書き換えれば実現できることは確認済みです。

  3. たとえば高機能なRAW現像ソフトウェアRawTherapeeではさまざまなアルゴリズムを実装しているようです。https://rawpedia.rawtherapee.com/Demosaicing/jp

  4. [1802.03769] Learning Deep Convolutional Networks for Demosaicing

ゼロから作るRAW現像 その3 - デモザイク処理基本編

このブログ記事「ゼロから作るRAW現像」を大きく再構成してより読みやすくした書籍「PythonとColabでできる-ゼロから作るRAW現像」を【技術書典6】にて頒布しました。

現在はBOOTHにて入手可能です。書籍+PDF版は2200円プラス送料、PDF版は1200円です。

moiz.booth.pm

はじめに

これは「ゼロから作るRAW現像 その1 - 基本的な処理」および「ゼロから作るRAW現像 その2 - 処理のモジュール化」の続きです。 これらの記事の内容を前提としていますので、まだお読みでない方はそちらからお読みいただくことをおすすめします。

簡易デモザイク処理の問題点

「ゼロから作るRAW現像 その1 - 基本的な処理」では、デモザイク処理(Bayer配列の画像からフルカラーの画像を作り出す処理)として、簡易的な画像サイズが1/4になるものを使いました。単純な処理の割に意外なほどきれいな出力が得られるのですが、いかんせん画像が小さくなるのは問題です。また、出力画像が1/4になるので、とうぜん細かい部分は潰れてしまいます。

たとえば、同じシーンを、カメラの出力するJPEGと、前回の簡易RAW現像処理で処理した画像とで比べてみましょう。

f:id:uzusayuu:20180930075749p:plain

左がraw_process.pyによる出力、右がカメラの出力したJPEGです。同じ倍率で表示しています。サイズの違いが一目瞭然ですね。

拡大してみましょう。

f:id:uzusayuu:20180930081702p:plain

文字や図形の大きさがほぼ同じになるように、左側の画像は800%、右の画像は400%に拡大してあります。 明るさやコントラストの違いがまっさきに目につきますが、それは次回以降考えましょう。解像度に注目すると、意外なほど健闘はしているのですが、縦のラインの分解能が低かったり、印刷の目が再現されていなかったり、といった点がわかると思います。このあたりは簡易デモザイクによる画像サイズの低下の影響といえるでしょう。

現代のカメラ内部のデモザイクはかなり高度な処理をしているはずなので、右側のJPEG画像並みの解像度を得るのは難しいと思いますが、せめてもとの画像サイズを取り戻せるような処理を導入してみましょう。

ベイヤー配列再訪

デモザイクの細部に入る前にベイヤー配列がどんなものか再確認しておきましょう。ベイヤー配列では、各画素が、赤、青、緑、のうち一色だけをもっています。 rawpyのraw_imageでは、配列は次のようにして確認できます。

bayer_pattern = raw.raw_pattern
print(bayer_pattern)

[[0 1]
 [3 2]]

ここで、各番号と色の関係は以下のようになっています。カッコ内は略称です

番号
0 赤 (R)
1 緑 (Gr)
2 青 (B)
3 緑 (Gb)

ここで緑にGrとGbがあるのは、赤の行の緑と青の行の緑を区別するためです。カメラ画像処理では両者を区別することが多々あり、両者をGrとGbと表す事が多いです。 両者を区別する必要が無い場合はどちらもGであらわします。

この対応関係を考えると、この画像の各画素の色は、左上から

赤 緑

緑 青

のように並んでいることがわかります。これを図示するとこうなります。

f:id:uzusayuu:20180930105005p:plain

では、元のRAWデータを、この色のまま再現してみましょう。

まずは、RAW画像を読み出し、ブラックレベル補正を行っておきます。

import os, sys
module_path = os.path.abspath(os.path.join('..'))
if module_path not in sys.path:
    sys.path.append(module_path)
import raw_process

raw = raw_process.read("sample2.ARW")
raw_array = raw_process.get_raw_array(raw)
blc_raw = raw_process.black_level_correction(raw, raw_array)

次にもともと赤があった位置の画素を赤、緑の位置を緑、青の位置を青であらわします。

import numpy as np

h, w = blc_raw.shape
bayer_img = np.zeros((h, w, 3))
bayer_pattern[bayer_pattern == 3] = 1
gain = [2, 1, 2]
for y in range(0, h):
    for x in range(0, w):
        color = bayer_pattern[y % 2, x % 2]
        bayer_img[y, x, color] = gain[color] * blc_raw[y, x]

元のままだと、画素の数が多い分緑がかってしまうので、赤と青の強度を倍にしてあります。

表示してみましょう。

f:id:uzusayuu:20180930130104p:plain

なんだかちょっとだけそれらしい絵がでていますね。

拡大します。

f:id:uzusayuu:20180930130223p:plain

やはりだめです。拡大するとまるでテレビの画面を虫眼鏡で拡大したようなタイル模様が見えてきます。

デモザイクはこのタイル模様からフルカラーの画像を再現する画像処理になります。ある意味カメラ画像処理でもっとも肝になる部分です。

線形補完法

デモザイクアルゴリズムの中で、縮小する方法の次に簡単なのは、線形補間法です。 線形補間というとものものしいですが、ようするに、距離に応じて間の値をとるわけです。たとえば、緑の画素ならこうなります。

f:id:uzusayuu:20180930113526p:plain

赤の画素ではこうです。

f:id:uzusayuu:20180930152512p:plain

青の画素でも、赤の場合と同じような補完を行います。

では実際やってみましょう。

def mirror(x, min, max):
    if x < min:
        return min - x
    elif x >= max:
        return 2 * max - x - 2
    else:
        return x

dms_img = np.zeros((h, w, 3))
bayer_pattern = raw.raw_pattern
for y in range(0, h):
    for x in range(0, w):
        color = bayer_pattern[y % 2, x % 2]
        y0 = mirror(y-1, 0, h)
        y1 = mirror(y+1, 0, h)
        x0 = mirror(x-1, 0, w)
        x1 = mirror(x+1, 0, w)
        if color == 0:
            dms_img[y, x, 0] = blc_raw[y, x]
            dms_img[y, x, 1] = (blc_raw[y0, x] + blc_raw[y, x0] + blc_raw[y, x1] + blc_raw[y1, x])/4
            dms_img[y, x, 2] = (blc_raw[y0, x0] + blc_raw[y0, x1] + blc_raw[y1, x0] + blc_raw[y1, x1])/4
        elif color == 1:
            dms_img[y, x, 0] = (blc_raw[y, x0] + blc_raw[y, x1]) / 2
            dms_img[y, x, 1] = blc_raw[y, x]
            dms_img[y, x, 2] = (blc_raw[y0, x] + blc_raw[y1, x]) / 2
        elif color == 2:
            dms_img[y, x, 0] = (blc_raw[y0, x0] + blc_raw[y0, x1] + blc_raw[y1, x0] + blc_raw[y1, x1])/4
            dms_img[y, x, 1] = (blc_raw[y0, x] + blc_raw[y, x0] + blc_raw[y, x1] + blc_raw[y1, x])/4
            dms_img[y, x, 2] = blc_raw[y, x]
        else:
            dms_img[y, x, 0] = (blc_raw[y0, x] + blc_raw[y1, x]) / 2
            dms_img[y, x, 1] = blc_raw[y, x]
            dms_img[y, x, 2] = (blc_raw[y, x0] + blc_raw[y, x1]) / 2

ここでmirror()は画像のヘリでxやyがはみ出ないように処理しています。

明るさを調整して表示してみます。

f:id:uzusayuu:20180930144416p:plain

それらしいものがでました。

つづけて残りのホワイトバランス、カラーマトリクス、ガンマ補正をかけます。

img_wb = raw_process.white_balance(raw, dms_img)
color_matrix = [1141, -205, 88, -52, 1229, -154, 70, -225, 1179]
img_ccm = raw_process.color_correction_matrix(img_wb, color_matrix)
rgb_image = raw_process.gamma_correction(img_ccm)

結果はこうなります。

f:id:uzusayuu:20180930131235p:plain

明るさコントラストに不満はありますが、フルカラーの画像ができました。

セーブしてみます。

raw_process.write(rgb_image, "output3.png")

ではセーブした画像をカメラのJPEGと比べてみましょう。

f:id:uzusayuu:20180930133057p:plain

解像度が低いのはともかく、細部にJPEGにはない余計な色が出ているのが気になります。このような余計な色は偽色などと呼ばれることがあります。英語ではColor Artifactなどといわれます。

他には、印刷の模様などのティテールも消えています。シャープさが足りないのはエッジ強調で多少もどせるかもしれませんが、ディテールの細かい部分を戻すのはむずかしいでしょう。なお、このような部分をテクスチャといいます。

やはりこのあたりは単純な線形補間の限界のようです。いくら処理の内容を見ていくのが目的で、最高の画質をもとめてはいないとはいえども、もう少し改善してから次のステップに行きたいところです。

FIRフィルターを利用した高速化

ちょっと最後にひとつだけ。

上記のデモザイクを実行してみるとわかると思うのですが、結構時間がかかります。一つ一つの画素の処理は単純でも、24M画素もあると流石に重くなります。 Pythonでできる範囲で高速化を考えてみたいと思います。

まず、緑画素の補完式をよく見てみると、これはこんな2次元FIRフィルターをかけているのと等価だという事がわかります。

[[   0, 1/4,   0],
 [ 1/4,   1, 1/4],
 [   0, 1/4,   0]]

つまり、2次元FIRフィルターをかけるライブラリーが使えるということです。具体的にはscipyを使って次のようにします。

from scipy import signal

blc_green = blc_raw.copy()
blc_green[(raw.raw_colors == 0) | (raw.raw_colors == 2)] = 0
g_flt = np.array([[0, 1/4, 0], [1/4, 1, 1/4], [0, 1/4, 0]])
green = signal.convolve2d(blc_green, g_flt, boundary='symm', mode='same')

ここでraw_colorsは各画素の色を表しており、blc_greenはBayerから緑画素のみを取り出したものです。

g_fltは上で示したFIRフィルターを表し、これがconvolve2dにより重畳されています。

結果はどうでしょう?明るさを調整してみてみましょう。

f:id:uzusayuu:20180930135927p:plain

なかなか良い感じです。

同様に赤と青チャンネルでは、対応するFIRフィルターはこのようになります。

[[1/4, 1/2, 1/4],
 [1/2,   1, 1/4].
 [1/4, 1/2, 1/4]]

まず赤画素を処理してみます

blc_red = blc_raw.copy()
blc_red[raw.raw_colors != 0] = 0
rb_flt = np.array([[1/4, 1/2, 1/4], [1/2, 1, 1/2], [1/4, 1/2, 1/4]])
red = signal.convolve2d(blc_red, rb_flt, boundary='symm', mode='same')

また明るさを調整して結果をみてみましょう。

f:id:uzusayuu:20180930140326p:plain

良さそうです。青画素も同様の処理ができるはずです。

では、三色分の処理を行ってみましょう。

dms_img2 = np.zeros((h, w, 3))

blc_green = blc_raw.copy()
blc_green[(raw.raw_colors == 0) | (raw.raw_colors == 2)] = 0
g_flt = np.array([[0, 1/4, 0], [1/4, 1, 1/4], [0, 1/4, 0]])
dms_img2[:, :, 1] = signal.convolve2d(blc_green, g_flt, boundary='symm', mode='same')

blc_red = blc_raw.copy()
blc_red[raw.raw_colors != 0] = 0
rb_flt = np.array([[1/4, 1/2, 1/4], [1/2, 1, 1/2], [1/4, 1/2, 1/4]])
dms_img2[:, :, 0] =  signal.convolve2d(blc_red, rb_flt, boundary='symm', mode='same')

blc_blue = blc_raw.copy()
blc_blue[raw.raw_colors != 2] = 0
dms_img2[:, :, 2] =  signal.convolve2d(blc_blue, rb_flt, boundary='symm', mode='same')

これで先程の1画素づつ処理するコードに比べてかなり速くなりました。

出力画像はこうなります。

f:id:uzusayuu:20180930141910p:plain

さらに、ホワイトバランス、カラーマトリクス、ガンマ補正をかけると、先ほどと同様の結果が得られます。

img_wb = raw_process.white_balance(raw, dms_img2)
color_matrix = [1141, -205, 88, -52, 1229, -154, 70, -225, 1179]
img_ccm = raw_process.color_correction_matrix(img_wb, color_matrix)
rgb_image = raw_process.gamma_correction(img_ccm)

f:id:uzusayuu:20180930142520p:plain

この処理は前回のraw_process.pyに追加して、あらたにraw_process2.pyというファイルとしてgithubにアップロードしてあります。

使用法はraw_process.pyと同様、

python3 raw_process2.py INPUTFILE [OUTPUTFILE] [MATRIX]

です。

まとめ

今回は、簡易的な縮小デモザイクを、比較的単純な線形補間アルゴリズムをつかったデモザイクで置き換えました。 jupyter上の例はSimple_Demosaic.ipynbとして、前述のraw_process2.pyと共にgithubにアップロードしてあります。

github.com

次回はもう少し複雑なデモザイクを紹介しようと思います。

次の記事

uzusayuu.hatenadiary.jp