Moiz's journal

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

ゼロから作る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 で確認できる。