Moiz's journal

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

ゼロから作るRAW現像 - Colabでの実行

はじめに

この記事は「ゼロから作るRAW現像」という一連の記事のフォローアップです。

uzusayuu.hatenadiary.jp

ここでは、記事中紹介された内容をColab上で実行する方法を紹介します。

なお、実行にあたってGoogleアカウントが必要です。

Colabについて。

ColabはGoogleのサービスで、ウェブブラウザだけでJupyter Notebookのような機能が使える便利なサービスです。

Colab自体、また、その使い方についてはからあげ氏の以下の記事が詳しいです。

karaage.hatenadiary.jp

からあげ氏は「ゼロから作るRAW現像」の内容をColabで実行されていて、ノートブックをColab上で公開もされています1

karaage.hatenadiary.jp

公開されているノートブックはこちらです。

その1−その4

その5−その8

今回は私もgithub上のファイルをアップデートしてColabで直接実行できるようにしたので、そちらのファイルを実行する方法を紹介します。

ColabでのRAW現像の実行 - シェアされたノートブックを使用

今回まずはその1の内容をColabで実行できるようにアップデートしてシェアしました。

その1-基本的な処理

その2 - 処理のモジュール化

(その3以降の内容は今後追加します。)

まず、リンクをクリックしてノートブックを開きます。

f:id:uzusayuu:20181121161425p:plain

このままでは実行できないので、FileメニューからSave a copy in Drive...を選択して自分の環境に保存します。

f:id:uzusayuu:20181121154848p:plain

次に、RuntimeメニューからRun allを実行します。

f:id:uzusayuu:20181121155233p:plain

これで、記事の内容が実行されます。

f:id:uzusayuu:20181121155856p:plain

CobalでのRAW現像の実行 - Githubからダウンロード

Github上のファイルをColabから開いて実行することもできます。

まずColabのWEBページをブラウザで開きます。

するとファイルを開くダイアログが表示されるのでGITHUBのタブを選択します。

f:id:uzusayuu:20181121160301p:plain

URLをコピーし貼り付けます。

その1-基本的な処理

その2 - 処理のモジュール化

f:id:uzusayuu:20181121160525p:plain

ダイアログ上の虫眼鏡マークを押すとGithubからノートブックファイルがダウンロードされます。

f:id:uzusayuu:20181121160643p:plain

COPY_TO_DRIVEをクリックして自分の環境にコピーを保存します。

f:id:uzusayuu:20181121160828p:plain

RuntimeメニューからRun allを選択すると、処理が行われます。

f:id:uzusayuu:20181121161005p:plain

まとめ

以上、「ゼロから作るRAW現像」の内容をColabから実行する方法を紹介しました。

最初にColab上での実行に成功さてノートブックをシェアしていただいたからあげ氏に改めて感謝します。


  1. からあげ氏には記事のデバッグも協力していただきました。感謝します。

ゼロから作るRAW現像 - まとめページ

ゼロから作るRAW現像

はじめに

この一連の記事は、RAW画像現像・カメラ画像処理の内容を実際の動作レベルで解説し、なるべくスクラッチからPython上で実行してみる事を目的としています。

そのために、解説する処理は最重要なものにしぼり、使用するアルゴリズムは一部の例外(デモザイク)を除き、もっとも簡単なものを選びました1

また、使用したインプットファイル、Jupyter Notebook上の実行内容はすべてGithubで公開しています。

記事の最後ではラズベリーパイのカメラで撮影したRAW画像からこんなRGB画像が作れるようになります。

f:id:uzusayuu:20181104160132p:plain:w600

オリジナルサイズ画像はこちら

この記事で扱うもの

  • 基本的なRAW現像処理・カメラ画像処理の流れ
  • Bayer画像からRGB画像出力までの各アルゴリズムのうち、基本的な最低限のものの解説
  • 解説した基本的なアルゴリズムPythonによる実装と処理例

この記事で扱わないもの

環境について

この記事で解説する内容は一般的なものですが、使用した画像ファイルは特定のカメラに依存しています。 他のカメラでもわずかな変更で同等の処理ができるとは予想されますが、検証はしていません。

使用カメラ

なお、使用したファイルはGithubからダウンロードできるので、これらのカメラをお持ちでなくても、紹介した処理の内容を実行することは可能です。

実行環境

この内容を再現するには通常のPC環境に加えて以下の環境が必用です。

  • python3が実行でき、以下のライブラリがインストールされた環境

    • rawpy
    • numpy
    • scipy
    • matplotlib
    • imageio
  • exiftool

    exiftoolは必須ではありませんが、あれば記事の内容を再現するのが容易になります。

  • Jupyter Notebookまたは他のPython対話ツール

    なお、後半の内容を実行するには比較的大きめなRAM必要なようです。 はっきりと何GB必要というのは環境に依存するので難しいですが、16GB以上あった方が良さそうです。

    また、一連の内容は僅かな変更でGoogle Colabでも実行可能なようです。 実際@karaage0703氏が、この記事の内容をcolabで実行したものをシェアされています。

    また、私の公開しているファイルでもその1およびその2に関してのみノートブックをColabに対応させてあります。

    Colabでの実行方法についてはこちらを参照ください。

    ゼロから作るRAW現像 - Colabでの実行 - Moiz's journal

記事一覧

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

カメラ画像処理の流れと、ごく基本的な処理について解説しています。

uzusayuu.hatenadiary.jp

「ゼロから作るRAW現像その1 - 基本的な処理」で使用したファイル

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

その1で行った処理をモジュール化し、その3以降への準備としています。

uzusayuu.hatenadiary.jp

「ゼロから作るRAW現像 その2 - 処理のモジュール化」で使用したファイル

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

基本的なデモザイク処理について解説しています。

uzusayuu.hatenadiary.jp

「ゼロから作るRAW現像 その3 - デモザイク処理基本編」で使用したファイル

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

比較的高性能なデモザイク処理について解説しています。 この記事だけ例外的に画像処理・信号処理についてのある程度専門的な知識を前提としています。 デモザイク内部の処理の解説なので、全体的な処理の流れをつかむのが目的の方は読み飛ばしてかまいません。

uzusayuu.hatenadiary.jp

「ゼロから作るRAW現像 その4 - デモザイク処理応用編」で使用したファイル

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

ラズベリーパイのカメラでRAW画像を撮影し、その画像を処理をする方法を解説しています。

uzusayuu.hatenadiary.jp

「ゼロから作るRAW現像 その5 - ラズベリーパイのRAW画像処理」で使用したファイル

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

レンズ/シェーディング補正(周辺減光補正)について解説しています。 地味ですが非常に重要な処理です。

uzusayuu.hatenadiary.jp

「ゼロから作るRAW現像 その6 - レンズシェーディング補正」で使用したファイル

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

ノイズ測定、ノイズのモデリング、バイラテラルフィルターを使った実際のノイズ除去について解説しています。

uzusayuu.hatenadiary.jp

「ゼロから作るRAW現像 その7- ノイズモデルとノイズフィルター」で使用したファイル

ゼロから作るRAW現像 その8- 欠陥画素、エッジ強調、コントラスト補正

最後の記事です。欠陥画素補正、エッジ強調、トーンカーブ補正について解説しています。

uzusayuu.hatenadiary.jp

「ゼロから作るRAW現像 その8- 欠陥画素、エッジ強調、コントラスト補正」で使用したファイル

ゼロから作るRAW現像 - Colabでの実行

フォローアップ記事です。

その1の内容をGoogle Colab上で実行する方法を紹介しています。

uzusayuu.hatenadiary.jp

Github リポジトリ

プロジェクトトップ

github.com


  1. こういった事情により、この記事で作成するRAW画像現像ソフトは最高の画質を目標にはしていません。

ゼロから作るRAW現像 その8- 欠陥画素、エッジ強調、コントラスト補正

はじめに

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

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

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

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

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

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

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

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

残りの処理について

前回まででだいぶ処理はすすみましたが、いくつか大事な処理が残っています。 主要な処理のうちまだ行っていないものとしては次のような物があります。

今回はこれらの処理をまとめて紹介し、実行してみようと思います。

準備

今回の処理の内容はgithubJupyter notebookファイル及びpythonファイルとしてアップロードしてあります

まず前回まで行ったシェーディング補正とノイズフィルターをモジュール化し、これまでのモジュールとあわせてraw_process5.pyとしました。

内容としては前回・前々回に説明したものと同一なので説明しませんが、以下の方法で処理を呼び出す事ができます。

シェーディング補正

lsc_coef = [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])]
lsc_raw = raw_process.lens_shading_correction(raw, lsc_coef)

ノイズフィルター

img_flt = raw_process.noise_filter(rgb_array, coef=8, read_noise=2, shot_noise=246)

また、ソニーARWとラズベリーパイRAWの両方に対応するために、ホワイトバランスの引数を変えました

ホワイトバランス

wb_raw = raw_process.white_balance_Bayer(lsc_raw, raw.camera_whitebalance, wb_norm, raw.raw_pattern)

ここでwb_normはゲインなしに対応する値です。ラズベリーパイでは1.0ですが、ソニーARWでは1024です。

また、raw_process自体を呼び出すことによりRaspberryPiのRAW画像をコマンドラインから直接RAW現像する事ができるようになりました。

raw_process5実行例

> python3 raw_process5.py chart.jpg test_out.png "6022,-2314,394,-936,4728,310,300,-4324,8126"

欠陥画素補正

前回まで処理してきたラズベリーパイの画像をよく見るとこんな部分があります。

f:id:uzusayuu:20181104074850p:plain:w300

これはいわゆる欠陥画素です。 「欠陥」という名前がついていますが、製品の欠陥ではなく、一部の画素が正常な値を出力しない状態です。 多くの場合このように、常に明るく見えますが(ホットピクセル)、常に暗い(デッドピクセルまたはコールドピクセル)こともあれば、本来の信号とずれた値を示すという捕まえにくいケースもあります。

数ミリ角のサイズに数百万から数千万の画素を作り込む現代の画像センサーでは、一部にこのような欠陥画素が含まれていることはごく普通のことです。

なぜこのような欠陥ができるかというのには、さまざまな原因が考えられます。 たとえば製造過程でパーティクルが入り1特定の画素が反応しなくなったとか逆にショートして常に電流が流れるようになった、というのがまっさきに思いつきます。 まだ半導体中の結晶欠陥などのせいで基板側に電流が漏れている(または基板から漏れて入ってくる)のかもしれません。 さらには宇宙線などが当たって製造後に欠陥が形成されるケースもあると聞きます。

このように、画像センサーに欠陥画素はつきもので、画像処理である程度対応していく必用があります。 市販のスマートフォンのカメラはもちろん、一眼レフカメラなどでも欠陥画素補正処理は内部的に行われているはずです。

それでは、実際に簡単な補正処理を行ってみます。

欠陥画素補正には大きく分けて2つのステップがあります2

  1. 欠陥画素検出
  2. 欠陥画素修正

今回は、欠陥画素検出としては、周辺の画素の最大値よりある程度以上大きいか最小値よりある程度以上小さければ欠陥とみなす、という方針で行います。 また修正としては、欠陥画素の上下左右4画素の平均をとる事にします。

それでは処理してみましょう。

まずは画像を読み込みブラックレベル補正まで行います。

%matplotlib inline
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_process5 as raw_process
from pylab import imshow, show
from matplotlib import pyplot as plt
from matplotlib.pyplot import imshow
import numpy as np

import rawpy
raw = rawpy.imread("chart.jpg")
h, w = raw.sizes.raw_height, raw.sizes.raw_width
img = np.array(raw.raw_image).reshape((h, w)).astype('int')
blc_raw = raw_process.black_level_correction(img, raw.black_level_per_channel, raw.raw_pattern)

次に各色毎に欠陥画素の検出と修正を行います。なお、ここではGr (Redと同列のGreen)とGb (Blueと同列のGreen)を別の色として処理しています。

まず、左上の色(Blue)からです。処理しやすいようにBlueのみを抜き出します。

dpc_raw = blc_raw.copy()
single_channel = dpc_raw[::2, ::2]

周辺の5x5画素の最大値と最小値を求めます。欠陥画素かどうか判定する対象の中央画素はfootprint機能を使って除いておきます。

footprint = np.ones((5, 5))
footprint[2, 2] = 0
local_max = scipy.ndimage.filters.maximum_filter(single_channel, footprint=footprint, mode='mirror')
local_min = scipy.ndimage.filters.minimum_filter(single_channel, footprint=footprint, mode='mirror')

この最大値と最小値を使って、欠陥画素判定を行います。 最大値や最小値との差がthreshold値より大きい場合欠陥画素とみなす事にします。

threshold = 16
mask = (single_channel < local_min - threshold) + (single_channel > local_max + threshold)

これで欠陥画素の位置がmaskTrueとして記録されました。(欠陥画素以外はFalse

それでは修正しましょう。

まず、欠陥画素の上下左右の画素の平均値を計算しておきます。

flt = np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]]) / 4
average = scipy.signal.convolve2d(single_channel, flt, mode='same')

次に欠陥画素をこの平均値で置き換えます。

single_channel[mask] = average[mask]

これでBlue面の欠陥画素補正ができました。

他の色の欠陥も補正するためにループ化します。

dpc_raw = blc_raw.copy()
footprint = np.ones((5, 5))
footprint[2, 2] = 0
for (yo, xo) in ((0, 0), (1, 0), (0, 1), (1, 1)):
    single_channel = dpc_raw[yo::2, xo::2]
    flt = np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]]) / 4
    average = scipy.signal.convolve2d(single_channel, flt, mode='same')
    local_max = scipy.ndimage.filters.maximum_filter(single_channel, footprint=footprint, mode='mirror')
    local_min = scipy.ndimage.filters.minimum_filter(single_channel, footprint=footprint, mode='mirror')
    threshold = 16
    mask = (single_channel < local_min - threshold) + (single_channel > local_max + threshold)
    single_channel[mask] = average[mask]

これで全色の欠陥画素を修正できたはずです。

残りの処理を行って確かめてみましょう。

lsc_coef = [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])]
lsc_raw = raw_process.lens_shading_correction(dpc_raw, lsc_coef)
wb_raw = raw_process.white_balance_Bayer(lsc_raw, raw.camera_whitebalance, 1.0, raw.raw_pattern)
dms_img = raw_process.advanced_demosaic(wb_raw, raw.raw_pattern)
img_flt = raw_process.noise_filter(dms_img, coef=8, read_noise=2, shot_noise=246)
ccm_matrix = (6022,-2314,394,-936,4728,310,300,-4324,8126)
img_ccm = raw_process.color_correction_matrix(img_flt, ccm_matrix)
white_level = 1024
img_gamma = raw_process.gamma_correction(img_ccm / white_level, 2.2)

確認します。

outimg = img_gamma.copy()
outimg[outimg < 0] = 0
outimg[outimg > 1] = 1
plt.imshow(outimg)
plt.axis('off')
plt.show()

f:id:uzusayuu:20181104074211p:plain

ファイルに書き出して欠陥が直っている事を確認します。

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

画像ビューワーを使って前回欠陥が見つかったところと同じ部分を拡大してみるとこうなっていました。

f:id:uzusayuu:20181104074954p:plain:w300

成功のようです。

エッジ強調

今の出力画像は少し解像感が低く、どうにもぼんやりして見えます。

f:id:uzusayuu:20181104110304p:plain

エッジ強調を使って改善してみましょう。

今回エッジ強調に使うのはアナログの時代から使われてきたアンシャープマスキングという手法です。 これは、入力画像をぼやけさせた画像をまず作り、そのぼやけさせた画像を元の画素から引いてやることで行います。 ぼやけさせた画像は入力画像より暗くしておく必用があります。

式で表すとこのようになります。入力画像{ g \left( x, y \right) }に対し、ぼやけさせた画像を{ f \left( x, y \right) }とすると、出力画像は

\begin{equation} h \left( x, y \right) = g \left( x, y \right) - a f \left( x, y \right) \end{equation}

となります。このままでは暗くなってしまうので、明るさを調整すると、こうなります。

\begin{equation} h \left( x, y \right) = g \left( x, y \right) + a \left( g \left( x, y \right) - f \left( x, y \right) \right) \end{equation}

後半の部分はハイパスフィルターになっており、結局元の画像に高周波成分を足し合わせるのと同じ処理になっていることがわかります。

では実際に処理してみましょう。

エッジ強調は輝度成分に対して行うことが多いので、まずRGB画像から輝度成分を分離しましょう。 輝度と色の成分を含む色空間としては、カメラやJPEGでは通常YCbCr空間が使われます。今回もYCbCrを使ってみましょう。 YCbCrというのは輝度信号(Y)と2つの色差信号(Cb, Cr)でフルカラーを表す方式です。

sRGBからJPEGで使われるYCbCr空間への変換マトリクスはこのようになっています。

\begin{equation} \left( \begin{array} 0.299 &0.587 &0.144 \\ -0.168736 &-0.331264 &0.5 \\ 0.5 &-0.418688 &-0.081312\\ \end{array} \right) \end{equation}

この他に通常CbとCr信号には+128のオフセットがのりますが、今回は省きました。

では、まずRGB信号をYCbCrに分解してみましょう。

rgb2ycbcr = np.array([[0.299, 0.587, 0.144], [-0.168736, -0.331264, 0.5], [0.5, -0.418688, -0.081312]])
# 0から255の範囲に調整
img_rgb = img_gamma * 256
img_rgb[img_rgb < 0] = 0
img_rgb[img_rgb > 255] = 255

# 色空間の変換
img_ycbcr = np.zeros_like(img_rgb)
for c in (0, 1, 2):
    img_ycbcr[:, :, c] = rgb2ycbcr[c, 0] * img_rgb[:, :, 0] + \
                                  rgb2ycbcr[c, 1] * img_rgb[:, :, 1] + \
                                  rgb2ycbcr[c, 2] * img_rgb[:, :, 2]

このうち輝度成分を取り出して確認してみましょう。

luma = img_ycbcr[:, :, 0]
plt.imshow(luma, cmap='gray')
plt.axis('off')
plt.show()

f:id:uzusayuu:20181104082948p:plain

うまく変換できているようです。

それではアンシャープマスクをかけてみましょう。

今回は2種類のガウシアンフィルターをかけてぼやけた画像を作成してみます。

# より低周波数成分
unsharpen4 = scipy.ndimage.gaussian_filter(luma, sigma = 2)
# やや高周波成分
unsharpen2 = scipy.ndimage.gaussian_filter(luma, sigma = 1)

このぼやけた画像を元の画像から引き、明るさを調整します。

ここで係数は任意のチューニングパラメータです。この2つの数字を変えると出力画像の性質が変わります。

sharpen = luma + 0.25 * (luma - unsharpen4) + 0.25 * (luma - unsharpen2)

結果を処理前と比較してみたのが次の画像です。右がアンシャープマスク後です。

f:id:uzusayuu:20181104084808p:plain

ラインの黒が深くなり、エッジも立っているのがわかると思います。

カラー画像に戻して確認します。

ycbcr2rgb = np.linalg.inv(rgb2ycbcr)
img_shp = img_ycbcr.copy()
img_shp[:, :, 0] = sharpen
img_out = np.zeros_like(img_shp)
for c in (0, 1, 2):
    img_out[:, :, c] = ycbcr2rgb[c, 0] * img_shp[:, :, 0] + \
                           ycbcr2rgb[c, 1] * img_shp[:, :, 1] + \
                           ycbcr2rgb[c, 2] * img_shp[:, :, 2]
img_out[img_out<0] = 0
img_out[img_out>255] = 255
raw_process.write(img_out / 256, "sharpened.png")

アンシャープマスク処理前と比較したのが次の画像です。右が処理後です。

f:id:uzusayuu:20181104085246p:plain

カラー画像でも、エッジが強調されているのがわかります。

トーンカーブ補正

最後にトーンカーブ補正です。

先程の画像の輝度成分のヒストグラムを見てみましょう。

plt.hist(sharpen.flatten(), 256)
plt.axis((0, 256, 0, 120000))
plt.show()

f:id:uzusayuu:20181104093225p:plain

特におかしなところはありませんが、全体に中央に集まっていて、せっかくの256階調のダイナミックレンジを十分には使い切っていないことがわかります。

このような画像のダイナミックレンジ拡張方法としてヒストグラム平坦化というものがあります。これは、ヒストグラムの積算値と同じ形をした関数を元の画像に適用することでヒストグラムを平坦にしてしまう、というものです。 しかし、ヒストグラムを完全に平坦にしてしまうと大概は不自然な画像になってしまいます。ここでは、そこまで極端でない補正をかけたいところです。

実際のカメラの中では非常に複雑なアルゴリズムトーンカーブ補正を計算しますが、その部分は今回の記事では対象ではありませんので、適当にちょうど良さそうな関数を設定してトーンカーブ補正の画像処理の部分だけ行ってみましょう。

先程のヒストグラムをみると中央部分が高いので、この部分をバラけさせるために、{x=128}付近で急になる関数をかけましょう。scipyのspline関数が使えそうです。まず適当なアンカーポイントを設定します。

xs = [0, 72, 128, 200, 256]
ys = [0, 56, 128, 220, 256]

図示するとこうなります。

plt.plot(xs, ys)
plt.axis((0, 256, 0, 256))
plt.show()

f:id:uzusayuu:20181104094114p:plain

この関数をそのまま使うと折れ曲がっている周辺で何らかのアーティファクトが発生しそうです。 スプライン関数を利用して、スムースな関数を用意しましょう。

func = scipy.interpolate.splrep(xs, ys)

この関数を図示するとこうなります。

xx = np.arange(0, 255)
yy = scipy.interpolate.splev(xx, func)
plt.plot(xx, yy)
plt.axis((0, 256, 0, 256))
plt.show()

f:id:uzusayuu:20181104094737p:plain

よさそうな関数ができました。輝度信号に適用してみましょう。

adjusted = scipy.interpolate.splev(sharpen, func)

plt.hist(adjusted.flatten(), 256)
plt.axis((0, 256, 0, 120000))
plt.show()

f:id:uzusayuu:20181104094853p:plain

先程よりもダイナミックレンジの全範囲を活用していることがわかります。 フルカラーに変換します3

img_shp = img_ycbcr.copy()
img_shp[:, :, 0] = adjusted
img_out = np.zeros_like(img_shp)
for c in (0, 1, 2):
    img_out[:, :, c] = ycbcr2rgb[c, 0] * img_shp[:, :, 0] + \
                           ycbcr2rgb[c, 1] * img_shp[:, :, 1] + \
                           ycbcr2rgb[c, 2] * img_shp[:, :, 2]
img_out[img_out<0] = 0
img_out[img_out>255] = 255

確認します。

plt.imshow(img_out / 256)
plt.axis('off')
plt.show()

f:id:uzusayuu:20181104095132p:plain

先程の画像よりちょっとだけパリッとしたような気がします。

セーブして処理前と比較してみましょう。

raw_process.write(img_out / 256, "tone_adjusted.png")

右側が補正後の画像です。

f:id:uzusayuu:20181104093030p:plain

黒い部分が締まり、全体により引き締まった画像になりました。

まとめ

今回は駆け足になりましたが、カメラ画像処理で重要な処理の中でまだ触れていなかった欠陥画素補正、エッジ強調、さらに、トーンカーブ補正をまとめて解説しPython上で実行しました。

いずれも各処理用のアルゴリズムとしてはもっとも簡単な種類のものですが、適切なパラメータを与えることで、アルゴリズムの単純さから考えると意外なほどの効果が得られることがわかると思います。

前述の通り今回の内容と使用したファイルはgithub上にアップロードしてあります。

github.com

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

使用したモジュール: raw_process/raw_process5.py at master · moizumi99/raw_process · GitHub

Jupyter Notebook ファイル: raw_process/part_8_sharpening.ipynb at master · moizumi99/raw_process · GitHub

最終画像: raw_process/tone_adjusted.png at master · moizumi99/raw_process · GitHub

このgithubレポジトリにあるraw_process5.pyには欠陥画素補正、エッジ強調、トーンカーブ補正も追加してあります。 これらの機能を使い今回Jupyter Notebook上で行ったのと同じ処理を行うには、raw_process5.pyとchart.jpgのあるディレクトリでコマンドラインから以下のように入力してください。

python3 raw_process5.py chart.jpg test_out.png "6022,-2314,394,-936,4728,310,300,-4324,8126" "0, 72, 128, 200, 256" "0, 56, 128, 220, 256"

最後に

カメラ画像処理をPython上で実際になるべく低レベルから書いて実行してみるという試みも、これでとうとう最後までたどり着きました。

もちろん実際のカメラやRAW現像ソフトで行われている処理は遥かに複雑で、さらに今回取り上げなかった処理も多数行われていますが、それでもカメラ画像処理に伴う最小限の処理はカバーしたといえるでしょう。

今回の記事で「ゼロから作るRAW現像」というシリーズは終了ですが、また画像処理関係で何か面白いテーマを見つけたらたまに追加していこうと思います!


  1. 半導体の製造レベルの話なので、肉眼では見えないサイズです。

  2. 事前に欠陥画素の位置がわかっている場合は1の部分は省くこともあります。

  3. 本来色差成分(Cb/Cr)も調整しなくてはならないのですが今回は省きます。

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

はじめに

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

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