このブログ記事「ゼロから作るRAW現像」を大きく再構成してより読みやすくした書籍「PythonとColabでできる-ゼロから作るRAW現像」を【技術書典6】にて頒布しました。
現在はBOOTHにて入手可能です。書籍+PDF版は2200円プラス送料、PDF版は1200円です。
はじめに
これは「ゼロから作るRAW現像 」という一連の記事の一つです。 これらの記事の内容を前提としていますので、まだお読みでない方はこちらの記事からお読みいただくことをおすすめします。
「ゼロから作るRAW現像 その３ - デモザイク処理基本編」
「ゼロから作るRAW現像 その４ - デモザイク処理応用編」
「ゼロから作るRAW現像 その5 - ラズベリーパイのRAW画像処理」
レンズシェーディング（周辺減光）
デジタルカメラに限らず、レンズを通して結像した画像は中央部より周辺部のほうが暗くなっています。
このような現象をレンズシェーディングや単にシェーディング、また日本語では周辺減光などといいます。英語ではLens ShadingやVignettingと呼ばれます。
このような事が起きてしまう第一の原因は、レンズを通してセンサーにあたる光の量が、センサーの中央部と、センサーの周辺部とで異なる事です。 良く説明に出されるのが二枚のレンズを持つ単純な系の場合です。
この図の例では、センサー中央部分にあたる光はレンズの開口部全体を通ってくるのに対して、センサー周辺部にあたる光はレンズの一部分しか通過できません。 結果的に中心部が明るく、周辺部が暗くなります。 実際のカメラの光学系はこれより遥かに複雑ですが、同様の理由により画像の周辺部分が暗くなります。
この要因の他に、センサー周辺部ではレンズやカバーの縁に遮られる光の量も増えます。また、光の入射角度によってレンズ等のコーティングと干渉する可能性もあります。 またデジタルカメラ独特の状況として、画像センサー上のフォトダイオードに到達する光の量が、光の入射角に依存するケースがあります。
シェーディングを起こす要因はこのように多岐にわたり、一概に、これが原因だとは言うことはできません。 したがってモデル計算で減光量を求めるよりも、実際のレンズで測定した結果を元に画像補正する必要があります。
レンズシェーディングの確認
ここで紹介する内容は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)
画像を見てみましょう。
確かに周辺光量が落ちているのが確認できます。
画像の中で明るさがどう変わっているか見てみましょう。 画面の高さ方向中央付近、上下３２画素幅で左から右まで帯状の画像をとりだし、明るさがどう変わるをグラフにしてみます。
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()
画像の左右端では中心部分に比べて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には中央からの距離の２乗が入っているはずです。
最大値でノーマライズしてグラフにして確認してみます。
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()
きれいに中心からの距離に応じて明るさが減少しています。 また、この段階では周辺部では中心部に比べて３分の１程度に暗くなっていることがわかります。
これを補正するには、明るさの減少率の逆数をかけてやればよい事になります。 逆数のグラフを書いてみましょう。
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()
これなら１次関数で近似できそうです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()
良さそうです。
レンズシェーディング補正
ではいよいよ、実際の画像のレンズシェーディングを補正してみましょう。
まず、レンズシェーディング補正前の、ブラックレベル補正のみをかけたRAW画像がこちらです。
先に各画素ごとに掛け合わせるゲインを、先程の近似関数から計算しておきます。
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')
フラットな画像が出力されました！
残りの処理（ホワイトバランス補正、デモザイク、カラーマトリクス補正、ガンマ補正）を行ってフルカラー画像を出力してみましょう。
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)
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()
シェーディングがほぼなくなりフラットになりました。 また、赤色と緑・青色との違いもほぼ消えました。成功のようです。
通常画像への適用
それではテスト用の平坦画像(Flat Field)ではなく、実際の画像にレンズシェーディング補正を適用してみましょう。
ラズベリーパイのカメラで改めて対象の画像(chart.jpg)をキャプチャします。
raspistill -r -o chart.jpg
以下の内容は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)
悪くはないのですが、全体に青みがかっています。また、上記の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)
先程の画像に比べると明るさが均一になり、不自然な青みも消えました。 成功のようです。
セーブしておきましょう。
raw_process.write(img_shading, "chart_with_shading_correction.png")
まとめ
今回はレンズシェーディング補正（周辺減光補正）をとりあげました。 おそらく、『カメラ』画像処理以外のいわゆる画像処理では取り上げることのない特殊な処理だと思います。
今回行ったのは、半径方向の二次多項式による補正ゲインの近似で、レンズシェーディング補正の中では最も単純なものです。 実際のカメラの中では、より高次の関数による近似や、２次元ルックアップテーブルによる補正などが行われているのが普通です。 また、補正パラメータも、明るさや光源の種類、オートフォーカスの場合はフォーカス位置、などにより調整します。
一見単純そうな見た目や効果と比べて、実際には遥かに複雑で非常に重要な処理です。ある意味カメラの出力画像の画質を決める肝と言ってもよいと思います。
今回の内容はテストデータと共に、githubにアップロードしてあります。
平坦画像 flat.jpg
解像度チャートchart.jpg
解像度チャートをRAW現像したものchart_with_shading_correction.png
平坦画像にシェーディング補正を施す例: LensShadingCorrection.ipynb
解像度チャートにシェーディング補正を施す例: Raspberry_Pi_RAW_LensShadingCorrection.ipynb
例の中で使用した自作ライブラリraw_process4.py
