このブログ記事「ゼロから作るRAW現像」を大きく再構成してより読みやすくした書籍「PythonとColabでできる-ゼロから作るRAW現像」を【技術書典6】にて頒布しました。
現在はBOOTHにて入手可能です。書籍+PDF版は2200円プラス送料、PDF版は1200円です。
moiz.booth.pm
はじめに
これは「ゼロから作るRAW現像 」という一連の記事の一つです。
これらの記事の内容を前提としていますので、まだお読みでない方はこちらの記事からお読みいただくことをおすすめします。
「ゼロから作るRAW現像 その1 - 基本的な処理」
「ゼロから作るRAW現像 その2 - 処理のモジュール化」
「ゼロから作るRAW現像 その3 - デモザイク処理基本編」
「ゼロから作るRAW現像 その4 - デモザイク処理応用編」
「ゼロから作るRAW現像 その5 - ラズベリーパイのRAW画像処理」
「ゼロから作るRAW現像 その6 - レンズシェーディング補正」
「ゼロから作るRAW現像 その7- ノイズモデルとノイズフィルター」
残りの処理について
前回まででだいぶ処理はすすみましたが、いくつか大事な処理が残っています。
主要な処理のうちまだ行っていないものとしては次のような物があります。
今回はこれらの処理をまとめて紹介し、実行してみようと思います。
準備
今回の処理の内容はgithubにJupyter 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"
欠陥画素補正
前回まで処理してきたラズベリーパイの画像をよく見るとこんな部分があります。
これはいわゆる欠陥画素です。
「欠陥」という名前がついていますが、製品の欠陥ではなく、一部の画素が正常な値を出力しない状態です。
多くの場合このように、常に明るく見えますが(ホットピクセル)、常に暗い(デッドピクセルまたはコールドピクセル)こともあれば、本来の信号とずれた値を示すという捕まえにくいケースもあります。
数ミリ角のサイズに数百万から数千万の画素を作り込む現代の画像センサーでは、一部にこのような欠陥画素が含まれていることはごく普通のことです。
なぜこのような欠陥ができるかというのには、さまざまな原因が考えられます。
たとえば製造過程でパーティクルが入り1特定の画素が反応しなくなったとか逆にショートして常に電流が流れるようになった、というのがまっさきに思いつきます。
まだ半導体中の結晶欠陥などのせいで基板側に電流が漏れている(または基板から漏れて入ってくる)のかもしれません。
さらには宇宙線などが当たって製造後に欠陥が形成されるケースもあると聞きます。
このように、画像センサーに欠陥画素はつきもので、画像処理である程度対応していく必用があります。
市販のスマートフォンのカメラはもちろん、一眼レフカメラなどでも欠陥画素補正処理は内部的に行われているはずです。
それでは、実際に簡単な補正処理を行ってみます。
欠陥画素補正には大きく分けて2つのステップがあります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)
これで欠陥画素の位置がmask
にTrue
として記録されました。(欠陥画素以外は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()
ファイルに書き出して欠陥が直っている事を確認します。
raw_process.write(img_gamma, "raspi_defect_correction.png")
画像ビューワーを使って前回欠陥が見つかったところと同じ部分を拡大してみるとこうなっていました。
成功のようです。
エッジ強調
今の出力画像は少し解像感が低く、どうにもぼんやりして見えます。
エッジ強調を使って改善してみましょう。
今回エッジ強調に使うのはアナログの時代から使われてきたアンシャープマスキングという手法です。
これは、入力画像をぼやけさせた画像をまず作り、そのぼやけさせた画像を元の画素から引いてやることで行います。
ぼやけさせた画像は入力画像より暗くしておく必用があります。
式で表すとこのようになります。入力画像に対し、ぼやけさせた画像をとすると、出力画像は
\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]])
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()
うまく変換できているようです。
それではアンシャープマスクをかけてみましょう。
今回は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)
結果を処理前と比較してみたのが次の画像です。右がアンシャープマスク後です。
ラインの黒が深くなり、エッジも立っているのがわかると思います。
カラー画像に戻して確認します。
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")
アンシャープマスク処理前と比較したのが次の画像です。右が処理後です。
カラー画像でも、エッジが強調されているのがわかります。
最後にトーンカーブ補正です。
先程の画像の輝度成分のヒストグラムを見てみましょう。
plt.hist(sharpen.flatten(), 256)
plt.axis((0, 256, 0, 120000))
plt.show()
特におかしなところはありませんが、全体に中央に集まっていて、せっかくの256階調のダイナミックレンジを十分には使い切っていないことがわかります。
このような画像のダイナミックレンジ拡張方法としてヒストグラム平坦化というものがあります。これは、ヒストグラムの積算値と同じ形をした関数を元の画像に適用することでヒストグラムを平坦にしてしまう、というものです。
しかし、ヒストグラムを完全に平坦にしてしまうと大概は不自然な画像になってしまいます。ここでは、そこまで極端でない補正をかけたいところです。
実際のカメラの中では非常に複雑なアルゴリズムでトーンカーブ補正を計算しますが、その部分は今回の記事では対象ではありませんので、適当にちょうど良さそうな関数を設定してトーンカーブ補正の画像処理の部分だけ行ってみましょう。
先程のヒストグラムをみると中央部分が高いので、この部分をバラけさせるために、付近で急になる関数をかけましょう。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()
この関数をそのまま使うと折れ曲がっている周辺で何らかのアーティファクトが発生しそうです。
スプライン関数を利用して、スムースな関数を用意しましょう。
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()
よさそうな関数ができました。輝度信号に適用してみましょう。
adjusted = scipy.interpolate.splev(sharpen, func)
plt.hist(adjusted.flatten(), 256)
plt.axis((0, 256, 0, 120000))
plt.show()
先程よりもダイナミックレンジの全範囲を活用していることがわかります。
フルカラーに変換します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()
先程の画像よりちょっとだけパリッとしたような気がします。
セーブして処理前と比較してみましょう。
raw_process.write(img_out / 256, "tone_adjusted.png")
右側が補正後の画像です。
黒い部分が締まり、全体により引き締まった画像になりました。
まとめ
今回は駆け足になりましたが、カメラ画像処理で重要な処理の中でまだ触れていなかった欠陥画素補正、エッジ強調、さらに、トーンカーブ補正をまとめて解説し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現像」というシリーズは終了ですが、また画像処理関係で何か面白いテーマを見つけたらたまに追加していこうと思います!