Moiz's journal

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

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

はじめに

これは「ゼロから作る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

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

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

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

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

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

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

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

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_gamma = raw_process.gamma_correction(img_ccm, 2.2)
outimg = img_gamma.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

良さそうです。

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

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

まず、補正前のブラックレベル補正した画像がこちらです。 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

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

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_gamma = raw_process.gamma_correction(img_ccm, 2.2)

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にアプロードしてあります

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

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)

表示してみましょう。

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))
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)

見てみましょう。

f:id:uzusayuu:20181021123223p:plain

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

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

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

まとめ

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

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

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

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


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

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

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

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

はじめに

これは「ゼロから作る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画像の抽出

今回の画像はこれまでと違い、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

with open("raw_capture.raw", "rb") as input_file:
    data = input_file.read()

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現像 」という一連の記事の一つです。 これらの記事の内容を前提としていますので、まだお読みでない方はこちらの記事からお読みいただくことをおすすめします。

モジュールの高速化

本題に入る前に、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 2 \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 2 \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, "output4.png")

画像の確認

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

f:id:uzusayuu:20181008024743p:plain

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

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

f:id:uzusayuu:20181008025157p:plain

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

まとめ

今回は、若干複雑なデモザイク処理を画像処理のレベルから解説し、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現像 その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

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

はじめに

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

RAW画像を処理するPythonスクリプト

これから処理を追加していくのが簡単になるように、その1で作ったRAW画像を処理するPythonスクリプトを単一のファイルにしてみましょう。

と言っても、一つ一つの処理をメソッドとして定義しなおして、mainメソッドからまとめて読み出しているだけですので、一つ一つの処理をここで追うのは止めにして、最終的な結果をgithubにアップロードしておきました。

github.com

使い方は

python3 raw_process.py 入力ファイル名 [出力ファイル名] [カラーマトリクス]

です。

出力ファイル名を省略するとoutput.pngというファイルに書き出されます。 カラーマトリクスは"1024, 0, 0, 0, 1024, 0, 0, 0, 1024"のような形式で書きます。省略した場合デフォルトの単位行列が使われます。

例えば前回と同じような現像処理を行うにはこんなふうに入力します

python3 sample.ARW output.png "1141, -205, 88, -52, 1229, -154, 70, -225, 1179"

実行結果はこちらです

f:id:uzusayuu:20180930053908p:plain

当然ながら前回とほぼ同じ結果が得られました。

スクリプトファイルの内部

このスクリプトの中では次のようなメソッドが定義されています。

メソッド名 説明
read(fileanme) filenameで指定されたRAWファイルからRAWデータを読み込みrawpy.raw_image形式で返す
get_raw_array(raw) raw_image形式のrawからデータ部分をnumpyのarray形式で取り出す
black_level_correction(raw, raw_array) raw_arrayに対しブラックレベル補正を行う
preview_demosaic(raw, raw_array) raw_arrayに対し簡易的でモザイク処理を行う。実行後の解像度は1/4になる。出力はRGBフルカラー
white_balance(raw, rgb_array) rgb_arrayに対してホワイトバランス補正を行う
color_correction_matrix(rgb_array, color_matrix) rgb_arrayに対してcolor_matrixの値を使ってカラーマトリクス補正を行う
gamma_correction(rgb_array) rgb_arrayに対してガンマ補正を行う
write(rgb_array, output_filename) rgb_arrayの内容をpng形式の画像ファイルとして書き出す

各メソッドの中では前回行った処理がそれぞれ行われています。 実際に処理を呼び出すmainメソッドの中ではこれらの処理を順番に呼び出しています。

以下はmainメソッドからの抜粋です。

    raw = read(filename)
    raw_array = get_raw_array(raw)
    blc_raw = black_level_correction(raw, raw_array)
    dms_img = preview_demosaic(raw, blc_raw)
    img_wb = white_balance(raw, dms_img)
    color_matrix = [1024, 0, 0, 0, 1024, 0, 0, 0, 1024]
    img_ccm = color_correction_matrix(img_wb, color_matrix)
    rgb_image = gamma_correction(img_ccm)
    write(rgb_image, output_filename)

Jupyterノートブックでの利用

上記のスクリプトファイルをモジュールとして利用することで、jupyterでの処理も簡単になります。

まず、jupyter上でmoduleを呼び出すために以下の処理を行います1

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_process

あとは、それぞれの処理を順番に呼び出すだけです。

raw = raw_process.read("sample.ARW")
color_matrix = [1024, 0, 0, 0, 1024, 0, 0, 0, 1024]
raw_array = raw_process.get_raw_array(raw)
blc_raw = raw_process.black_level_correction(raw, raw_array)
dms_img = raw_process.preview_demosaic(raw, blc_raw)
img_wb = raw_process.white_balance(raw, dms_img)
img_ccm = raw_process.color_correction_matrix(img_wb, color_matrix)
rgb_image = raw_process.gamma_correction(img_ccm)
raw_process.write(rgb_image, "output2.png")

ずいぶんすっきりしました。

画像がちゃんと出力できたか、確認しておきましょう。

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

f:id:uzusayuu:20180930062245p:plain

この内容はこちらにアップロードしてあります。

まとめ

今後の処理の追加が用意になるように、その1で行ったRAW画像処理の内容を実行可能な単一ファイルにまとめ、モジュール化しました。 これでJupyterから処理を呼び出すのも簡単になります。 次はこのモジュールを使って、他の処理を実装する予定です

なお、今回の結果はgithub上にアップロードしてあります。

github.com

次の記事

uzusayuu.hatenadiary.jp