Moiz's journal

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

ゼロから作る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, "advanced_demosaic.png")

画像の確認

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

f:id:uzusayuu:20181008024743p:plain:w800

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

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

f:id:uzusayuu:20181008025157p:plain:w800

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

まとめ

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

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

はじめに

会社の同僚にrawpyというPython用のライブラリの存在を教えてもらいました。 これを使うと各種デジカメのRAWファイルから、BayerのRAWデータを抽出できます。以前はDCRAWのソースコードを改造してrawデータのダンプなどしていたのですが、ずいぶん良い時代になったものです。 せっかくなのでrawpyで抽出したBayerフォーマット画像データから、普通の画像ビューワーなどで表示できる画像ファイルをできるだけスクラッチから作成して見たいと思います。

今回の内容はJupyter Notebook ファイルとしてgithubにアップロードしてあります。RAWデータNotebook ファイルをダウンロードしてJupyter Notebookから開くことで同様の処理が行なえます

RAWファイルおよびRAWデータについて

RAWファイルやRAWデータというのは厳密な定義はないのですが、カメラ処理でRAWというとBayerフォーマットの画像データを指すことが多いようです。したがって多くの場合、RAWデータはBayerフォーマットの画像データ、RAWファイルはそのRAWデータを含んだファイルということになります。

まず前提として、今使われているカメラの画像センサーの殆どはBayer配列というものを使ってフルカラーを実現しています。

画像センサーは、碁盤の目状にならんだ小さな光センサーの集まりでできています。一つ一つの光センサーはそのままでは色の違いを認識できません。そこで色を認識するためには、3原色のうち一色を選択して光センサーにあてて、その光の強度を測定する必要があります。 方法としてはまず、分光器を使って光を赤、青、緑に分解して、3つの画像センサーにあてて、それぞれの色の画像を認識し、その後その3枚をあわせることでフルカラーの画像を合成するという方法がありました。これは3板方式などとよばれることがあります。この方法は手法的にもわかりやすく、また、余計な処理が含まれないためフルカラーの画像がきれい、といった特徴があり高級ビデオカメラなどで採用されていました。欠点としては分光器と3つの画像センサーを搭載するためにサイズが大きくなるという点があります。

これに対して、画像センサー上の一つ一つの光センサーの上に、一部の波長の光だけを通す色フィルターを載せ、各画素が異なる色を取り込むという方法もあります。この方法では1枚の画像センサーでフルカラー画像を取り込めるため、3板方式に対して単板方式とよばれることもあります。3版方式とくらべた利点としては分光器が不要で1枚のセンサーで済むのでサイズが小さい。逆に欠点としては、1画素につき1色の情報しか無いので、フルカラーの画像を再現するには画像処理が必要になる、という点があります。

単版方式の画像センサーの上に載る色フィルターの種類としては、3原色を通す原色フィルターと、3原色の補色(シアン・マゼンダ・イエロー)を通す補色フィルター1というものがあります。補色フィルターは光の透過率が高いためより明るい画像を得ることができます。それに対して原色フィルターは色の再現度にすぐれています。Bayer配列はこの単版原色フィルター方式のうち最もポピュラーなものです。

こういうわけでBayer配列の画像センサーの出力では1画素につき一色しか情報をもちません。Bayer配列のカラーフィルターはこの図左のように、2x2ブロックの中に、赤が1画素、青が1画素、緑が2画素ならぶようになっています。緑は対角線上にならびます。緑が2画素あるのは、可視光の中でも最も強い光の緑色を使うことで解像度を稼ぐため、という解釈がなされています。Bayerというのはこの配列の発明者の名前です。

f:id:uzusayuu:20180923132843p:plain

カメラ用センサーでは2000年代初頭までは、補色フィルターや3板方式もそれなりの割合で使われていたのですが、センサーの性能向上やカメラの小型化と高画質化の流れの中でほとんどがBayer方式にかわりました。今では、SigmaのFoveonのような意欲的な例外を除くと、DSLRやスマートフォンで使われているカラー画像センサーの殆どがBayer方式を採用しています。したがって、ほとんどのカメラの中ではBayerフォーマットの画像データをセンサーから受取り、フルカラーの画像に変換するという処理が行われている、ということになります。

こういったBayerフォーマットの画像ファイルは、すなわちセンサーの出力に近いところで出力されたことになり、カメラが処理したJPEGに比べて以下のような利点があります。

  • ビット数が多い(RGBは通常8ビット。Bayerは10ビットから12ビットが普通。さらに多いものもある)
  • 信号が線形(ガンマ補正などがされていない)
  • 余計な画像処理がされていない
  • 非可逆圧縮がかけられていない(情報のロスがない)

したがって、優秀なソフトウェアを使うことで、カメラが出力するJPEGよりもすぐれた画像を手に入れる事ができる可能性があります。

逆に欠点としては

  • データ量が多い(ビット数が多い。通常非可逆圧縮がされていない)
  • 手を加えないと画像が見れない
  • 画像フォーマットの情報があまりシェアされていない
  • 実際にはどんな処理がすでに行われているのか不透明

などがあります。最後の点に関して言うと、RAWデータといってもセンサー出力をそのままファイルに書き出すことはまずなく、欠陥画素除去など最低限の前処理が行われいるのが普通です。 しかし、実際にどんな前処理がおこなわれているのかは必ずしも発表されていません。

カメラ画像処理

Bayerからフルカラーの画像を作り出すカメラ画像処理のうち、メインになる部分の例はこんな感じになります2f:id:uzusayuu:20180923122237p:plain このうち、最低限必要な処理は、以下のものです。

  • ブラックレベル補正
  • ホワイトバランス補正(デジタルゲイン補正も含む)
  • デモザイク(Bayerからフルカラー画像への変換)
  • ガンマ補正

これらがないと、まともに見ることのできる画像を作ることができません。

さらに、最低限の画質を維持するには、通常は、

  • 線形性補正
  • 欠陥画素補正
  • 周辺減光補正
  • カラーマトリクス

が必要です。ただし、線形性補正や欠陥画素補正は、カメラがRAWデータを出力する前に処理されていることが多いようです。また、センサーの特性がよければ線形性補正やカラーマトリクスの影響は小さいかもしれません。

次に、より良い画質を実現するものとして、

  • ノイズ除去
  • エッジ強調・テクスチャ補正

があります。RGB->YUV変換はJPEGMPEGの画像を作るのには必要ですが、RGB画像を出力する分にはなくてもかまいません。

この他に、最近のカメラでは更に画質を向上させるために

  • レンズ収差補正
  • レンズ歪補正
  • 偽色補正
  • グローバル・トーン補正
  • ローカル・トーン補正
  • 高度なノイズ処理
  • 高度な色補正
  • ズーム
  • マルチフレーム処理

などの処理が行われるのが普通です。今回はベーシックな処理のみとりあげるので、こういった高度な処理は行いません。

結局、今回扱うのは次の部分のみです3

f:id:uzusayuu:20180923160440p:plain

準備

まずRAW画像を用意します。今回はSony α 7 IIIで撮影したこの画像を使います。 f:id:uzusayuu:20180923124230p:plain

使用するRAWファイルはこちらからダウンロードできます。 https://github.com/moizumi99/raw_process/blob/master/sample.ARW

次にPython3が実行できる環境を用意します。今回はUbuntu18.04上のPython3.6を使用しました。

さらにPythonのライブラリとしてrawpyが必要です。pipが導入してあれば、次のコマンドでインストールできます。

pip install rawpy

他に、以下ののライブラリも必要なので導入済みでなければインストールしておいてください。

  • matplotlib
  • PIL
  • numpy
  • imageio
  • math

また必ずしも必須ではありませんが、Jupyter Notebookが使えれば以下で解説する内容を実行するのが楽になると思います。

今回の実行例はすべてJupyter Notebook上でのもので、githubにアップロードしてあります。 ローカルにダウンロード後、Jupyter上からファイルを開けば同様の内容が実行できるはずです。

なお処理自体に関係ありませんが、いくつかのパラメータを取得するのにexiftoolsを使っています。

RAW画像読み込み

では、raw画像を読み込んでみましょう。まず、sample.ARWをダウンロードしたディレクトリで、Jupyterを起動します。 JupyterがなければPythonの対話ウィンドウで同様の事ができると思います。

%matplotlib inline
import rawpy
raw = rawpy.imread('sample.ARW')

これでrawpyを使って画像が読み込めます。簡単ですね。 ちゃんと読めたか確認しましょう。

from matplotlib.pyplot import imshow
img_preview = raw.postprocess(use_camera_wb=True)
imshow(img_preview)

f:id:uzusayuu:20180923140802p:plain

実はrawpyにはLibRaw(内部的にはdcraw)を利用したraw画像現像機能があるので、このようにして現像後の画像を見ることができます。 ただ、今回の目的はBayerから最終画像フォーマットまでの処理を一つ一つ追いかけていくことなので、これはあくまで参考とします。 最終的にこの画像と似たようなものが得られれば成功としましょう。

画像データ変換

扱いやすいように画像をnumpyのarrayに変換しておきましょう。まず、RAWデータのフォーマットを確認しておきます。

print(raw.sizes)

ImageSizes(raw_height=4024, raw_width=6048, height=4024, width=6024, top_margin=0, left_margin=0, iheight=4024, iwidth=6024, pixel_aspect=1.0, flip=0)

RAWデータのサイズは4024 x 6048のようです。numpy arrayにデータを移しましょう。

import numpy as np
h, w = raw.sizes.raw_height, raw.sizes.raw_width
raw_image = raw.raw_image.copy()
raw_array = np.array(raw_image).reshape((h, w)).astype('float')

このデータを無理やり画像として見ようとするとこうなります。

outimg = raw_array.copy()
outimg = outimg.reshape((h, w))
outimg[outimg < 0] = 0
outimg = outimg / outimg.max()
imshow(outimg, cmap='gray')

f:id:uzusayuu:20180925140131p:plain

拡大するとこんな感じです

f:id:uzusayuu:20180925125658p:plain

明るいところが緑、暗いところが赤や青の画素のはずです。でも、これじゃなんのことやらさっぱりわかりませんね。 そこで、デモザイク処理でフルカラーの画像を作る必要があるわけです。

ブラックレベル補正

RAWデータの黒に対応する値は通常0より大きくなっています。画像を正常に表示するにはこれを0にもどす必要があります。これをやって置かないと黒が十分黒くない、カスミがかかったような眠い画像になってしまいますし、色もずれてしまいます。

まず、rawpyの機能でブラックレベルを確認しましょう。

blc = raw.black_level_per_channel
print(blc)

[512, 512, 512, 512]

どうやら全チャンネルでブラックレベルは512のようですが、他のRAWファイルでもこのようになっているとは限りません。各画素ごとのチャンネルに対応した値を引くようにしておきましょう。

さて、先程とりだしたrawデータの配列は以下の方法で確認できます。

bayer_pattern = raw.raw_pattern
print(bayer_pattern)

[[0 1]
 [1 2]]

どうやら0が赤、1が緑、2が青をしめしているようです。つまり、画像データの中から2x2のブロックを重複なくとりだすと、その中の左上が赤、右下が青、その他が緑、という事のようです。 このチャンネルに合わせて正しいブラックレベルを引きます。

print(raw_array.min(), raw_array.max())
blc_raw = raw_array.copy()
for y in range(0, h, 2):
    for x in range(0, w, 2):
        blc_raw[y + 0, x + 0] -= blc[bayer_pattern[0, 0]]
        blc_raw[y + 0, x + 1] -= blc[bayer_pattern[0, 1]]
        blc_raw[y + 1, x + 0] -= blc[bayer_pattern[1, 0]]
        blc_raw[y + 1, x + 1] -= blc[bayer_pattern[1, 1]]
print(blc_raw.min(), blc_raw.max())

0.0 8180.0
-512.0 7668.0

処理後の画像は、最大値と最小値が512小さくなっているのが確認できました。

outimg = blc_raw.copy()
outimg = outimg.reshape((h, w))
outimg[outimg < 0] = 0
outimg = outimg / outimg.max()
imshow(outimg, cmap='gray')

f:id:uzusayuu:20180925140209p:plain

画像も少し暗くなっています。

簡易デモザイク

次にとうとうBayer配列からフルカラーの画像を作ります。これが終わるとやっと画像がまともに確認できるようになります。

この処理はデモザイクと呼ばれることが多いです。本来デモザイクはカメラ画像処理プロセス(ISP)の肝になる部分で、画質のうち解像感や、偽色などの不快なアーティファクトなどを大きく左右します。 したがって手を抜くべきところではないのですが、今回は簡易処理なので、考えうる限りでもっとも簡単な処理を採用します。

その簡単な処理というのは、ようするに3色の情報を持つ最小単位の2x2のブロックから、1画素のみをとりだす、というものです。

f:id:uzusayuu:20180923134843p:plain

結果として得られる画像サイズは1/4になりますが、もとが24Mもあるので、まだ6M残っています。今回の目的には十分でしょう。 なお、解像度低下をともなわないデモザイクアルゴリズムは次回以降とりあげようと思います。

では、簡易デモザイク処理してみましょう。2x2ピクセルの中に2画素ある緑は平均値をとります。

for y in range(0, h, 2):
    for x in range(0, w, 2):
        colors = [0, 0, 0, 0]
        colors[bayer_pattern[0, 0]] += blc_raw[y + 0, x + 0]
        colors[bayer_pattern[0, 1]] += blc_raw[y + 0, x + 1]
        colors[bayer_pattern[1, 0]] += blc_raw[y + 1, x + 0]
        colors[bayer_pattern[1, 1]] += blc_raw[y + 1, x + 1]
        dms_img[y // 2, x // 2, 0] = colors[0]
        dms_img[y // 2, x // 2, 1] = (colors[1] + colors[3])/ 2
        dms_img[y // 2, x // 2, 2] = colors[2]

さてこれでフルカラーの画像ができたはずです。見てみましょう。

outimg = dms_img.copy()
outimg = outimg.reshape((h // 2, w //2, 3))
outimg[outimg < 0] = 0
outimg = outimg / outimg.max()
imshow(outimg)

f:id:uzusayuu:20180923154414p:plain

でました。画像が暗く、色も変ですが、それは予定通りです。そのあたりをこれから直していきます。

ホワイトバランス補正

次にホワイトバランス補正をかけます4。これは、センサーの色ごとの感度や、光のスペクトラムなどの影響を除去して、本来の白を白として再現するための処理です。 そのためには各色の画素に、別途計算したゲイン値をかけてあげます。今回はカメラが撮影時に計算したゲイン値をRAWファイルから抽出して使います。

ますはどんなホワイトバランス値かみてみましょう。RAWファイルの中に記録されたゲインを見てみましょう。

wb = np.array(raw.camera_whitebalance)
print(wb)

[2288. 1024. 1544. 1024.]

これは赤色にかけるゲインがx2288/1024、緑色がx1.0、青色がx1544/1024、という事のようです。処理してみましょう。

img_wb = dms_img.copy().flatten().reshape((-1, 3))
for index, pixel in enumerate(img_wb):
    pixel = pixel * wb[:3] /max(wb)
    img_wb[index] = pixel

f:id:uzusayuu:20180923154505p:plain

色がだいぶそれらしくなりました。

カラーマトリクス補正

次にカラーマトリクス補正を行います。カラーマトリクスというのは処理的には3x3の行列に、3色の値を成分としたベクトルをかけるという処理になります。

f:id:uzusayuu:20180925135315p:plain

なぜこんな事をするかというと、カメラのセンサーの色ごとの感度が人間の目とは完全には一致しないためです。例えば人間の目はある光の周波数の範囲を赤、青、緑、と感じるのですが、センサーが緑を検知する範囲は人間が緑と感じる領域とは微妙に異なっています。同じように青や赤の範囲も違います。これは、センサーが光をなるべく沢山取り込むため、だとか、製造上の制限、などの理由があるようです。 さらに、人間の目には、ある色を抑制するような領域まであります。これはセンサーで言えばマイナスの感度があるようなものですが、そんなセンサーは作れません。

こういったセンサー感度と人間の目の間隔とがなるべく小さくなるように、3色を混ぜて、より人間の感覚に近い色を作り出す必要があります。この処理を通常は行列を使って行い、これをカラーマトリクス処理と呼びます。

さて、ここでrawpyを使ってRAWデータの中に含まれる、マトリクスを調べるとこんなふうになってしまいます。

print(raw.color_matrix)

[[0. 0. 0. 0.]
 [0. 0. 0. 0.]
 [0. 0. 0. 0.]]

仕方がないのでexiftoolsでマトリクスを元のARWファイルから取り出したところ、次のような値のようです。

Color Matrix                    : 1141 -205 88 -52 1229 -154 70 -225 1179

この値を使って処理を行いましょう。

color_matrix = np.array([[1141, -205, 88], [-52, 1229, -154], [70, -225, 1179]]) / 1024

img_ccm = np.zeros_like(img_wb)
for index, pixel in enumerate(img_wb):
    pixel = np.dot(color_matrix, pixel)
    img_ccm[index] = pixel

f:id:uzusayuu:20180923154644p:plain

あまり影響が感じられません。元のセンサーが良いので極端な補正はかけなくてよいのかもしれません。

ガンマ補正

最後にガンマ補正をかけます。 ガンマ補正というのは、もともとテレビがブラウン管だった頃にテレビの出力特性と信号の強度を調整するために使われていたものです。 今でも残っているのは、ガンマ補正による特性が結果的に人間の目の非線形的な感度と相性が良かったからのようです。 そんなわけで現在でもディスプレイの輝度は信号に対してブラウン管と似たような信号特性を持って作られており、画像にはガンマ補正をかけておかないと出力は暗い画像になってしまいます。

ガンマ特性自体は次の式で表されます

 y = x^{2.2}

f:id:uzusayuu:20180925135153p:plain

ガンマ補正はこれを打ち消す必要があるので、このようになります。

 y = x^{\frac{1}{2.2}}

f:id:uzusayuu:20180925135208p:plain

やってみましょう。

import math

img_gamma = img_ccm.copy().flatten()
img_gamma[img_gamma < 0] = 0
img_gamma = img_gamma/img_gamma.max()
for index, val in enumerate(img_gamma):
    img_gamma[index] = math.pow(val, 1/2.2)
img_gamma = img_gamma.reshape((h//2, w//2, 3))

f:id:uzusayuu:20180923154725p:plain

だいぶきれいになりました。

まとめ

今回は簡易的な処理でしたが、元のRAW画像のデータが高品質なのでこの程度の画質を実現することができました。

最後にセーブしておきます。

import imageio

outimg = img_gamma.copy().reshape((h // 2, w //2, 3))
outimg[outimg < 0] = 0
outimg = outimg * 255
imageio.imwrite("sample.png", outimg.astype('uint8'))

最後に

今回はrawpyを使ってカメラのRAWファイルからBayerデータを取り出し、Pythonでできるだけスクラッチから簡易的なカメラ画像処理を作成して、RAW現像を行いました。 今の所、周辺減光補正、ノイズ処理、エッジ強調がない、など主要な処理が抜けていますし、デモザイクは簡易的なものですので、次回以降こういった処理を追加していこうと思います。

今回使用したRAWファイル、Jupyter notebookでの実行例、出力したPNGファイルはすべてGitHubにアップロードしてあります。

github.com

次の記事

uzusayuu.hatenadiary.jp

改定履歴

2018-10-28: @karaage氏のご厚意により、記事中ブラックレベル補正の余分なコードとデモザイクのバグを修正しました。 2018-10-31: 再び@karaage氏の指摘に基づきデモザイクのコードをより汎用性の高いものに変更しました。


  1. 実際にはこの他に緑色の画素もあり、2x2の4画素のパターンになっているのが普通

  2. これはあくまで一例です。実際のカメラ内で行われる処理はメーカーや機種ごとに異なる可能性があります。

  3. 周辺減光補正は本来必要な処理ですが、今回は影響が少ない事やメタデータの解析が必要な事もあり、対象から省きました。

  4. ダイアグラムでの処理の順番に比べて、デモザイクとホワイトバランスの順番が逆になっていますが、今回採用した簡易デモザイクでは影響ありません。

Raspberry Piのベアメタル環境からSDCARDにアクセスする(BUSモード編)

はじめに

Raspberry Piのベアメタル環境からGPIOを経由してSD CARDをアクセスする方法について、BUSモードでも成功したので紹介します

BUSモード

前回のエントリーではRaspberry PiでOSなどを介さず、EMMC 1も使わず、GPIO機能のみを使ってSD CARDにSPIモードでアクセスする方法を紹介しました。 しかし、SPIモードの使用には

  1. 一旦BUSモードで動作すると、SPIモードには再遷移できない
  2. Raspberry PiブートローダーはSD CARDにBUSモードでアクセスる
  3. したがって、デモでSPIモードを使用する前にユーザーがカードを一度抜き差しして、カードをリセットする必要がある

という問題がありました。 また、SD CARDの規格としても、SPIモードはあくまで簡易的にデータにアクセスするための方法です。したがってできることならBUSモードでSD CARDにアクセスしたいと思っていました。

例えば、SPIモードの特徴は以下のとおりです

  1. 使用するピンが少ない。電源とグラウンドを除くと、クロック1ビット、コマンド1ビットとデータ1ビットの、3本の信号線のみ2
  2. 制御が比較的簡単
  3. ただしデータが1ビットしかないので遅い

それに対してBUSモードは

  1. 最大4ビットのデータ線が使える
  2. SD CARDの最近の高速モードに対応している
  3. したがって速い
  4. ただし、制御はSPIモードに比べて若干面倒

となります3。 このように、BUSモードこそがSD CARD本来のモードと言えると思います。

今回は高速化を期待してBUSモードに対応させてみました。

リポジトリと使用方法

今回のデモの内容はGithubで公開してあります。

github.com

使用方法は

  1. 必要ないSD CARDを用意する4
  2. SD CARDにラズビアンOSなどをインストールしてRaspberry Piが起動するようにしておく
  3. SD CARDのboot領域のkernel.imgをリポジトリのkernel.imgで置き換える
  4. SD CARDをRaspberry Piに挿して起動する

なお、デモの結果を確認するにはSerialのTX/RXをRaspberry PiのGPIO#14/15(PIN8とPIN10)に接続し、PCなどのターミナルコンソールに信号の内容を表示させる必要があります。

デモの出力結果は以下のようになります。SD CARDのFAT32領域にアクセスして、ルートディレクトリのファイルを表示しています f:id:uzusayuu:20180730070201p:plain

前回との違い

SPIモードとBUSモードの違いは上記のようにいろいろありますが、今回影響を強く受けたのは以下のような点です

  1. コマンド(CMD)が双方向になった。SPIモードではCMDはホストからカードへの一方通行だったのですが、BUSモードではカードからホストにレスポンスを返すのにも使います
  2. データビットが1ビットから4ビットに増えた。BUSモードでは4本あるデータ線をすべて使うことができます
  3. コマンドやレスポンスのフォーマットがSPIモードと違う。

とくに、コマンドに対するレスポンスがデータと同時に帰ってくることがあるため、両者を同時に処理する必要がある点に注意が必要です。

なお、今回は読み込みまでしか実装していません。

実装の詳細についてはgithubで公開しているソースコードを参照ください5

次への課題

この実験は、もともとは、既存のライブラリで動作していた「Haribote OS」のSD CARD読み込み部分を、自分の実装で置き換えるというのが目的でした。これからそのあたりやっていこうと思っていたのですが、一時帰国で日本に帰った際にこんな本を見つけてしまいました。FatFSという組み込み向けFATライブラリーの解説本です。

shop.cqpub.co.jp

Haribote OSではFATの扱いはごく簡単なもので、残りの実装は読者への宿題になっています。こちらの本で解説されているFatFSを使えば読み書きやディレクトリのアクセスも含めたファイルシステムをサポートすることができそうです。ちょっと気になってきました。

まとめ

Raspberry Piのベアメタル環境でGPIOを経由してSD CARDをアクセスするデモについて紹介しました。 なお今回はリードのみの実装です。また、UHSなどの高速化モードはいっさいサポートしていません。


  1. Raspberry PiのSoC(BCM2835)のSD CARDインターフェース機能

  2. 規格としてはCS(チップセレクト)信号もあるが、Raspberry PiではSD CARDは一つしか繋がらないのでCSは必要ない

  3. 詳しいことはSD Associationが発行している簡易版規格書のパート1をご覧ください。

  4. デモでSD CARDの内容が破壊される可能性もありますので、必ず必要ないSD CARDを使用してください

  5. うまくまとまったらブログに書くことも考えています