Moiz's journal

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

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

はじめに

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

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

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

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

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

ラズベリーパイのRAW画像

これまでRAW現像の対象として、ソニーα7iiiで撮影した画像を使ってきました。 これは基本的な処理の説明をするにあたって、ノイズが少ない、歪がない、解像感が高い、など質の良いRAW画像のほうがやりやすかったからです。 しかし、実際のカメラ画像処理では、ノイズや意図しない不鮮明さなど好ましくない性質を修正するというのが大きな目的になります。そういった処理を行うにあたって、ミラーレスカメラの高画質データは逆に扱いにくいものになってしまいます。なにしろ直したいノイズもボケもあまりありませんから。

今回はラズベリーパイでRAW画像をキャプチャーして、これまで行ったRAW現像処理を行い、次回以降のテーマの準備を行います。

ラズベリーパイの準備とRAW画像のキャプチャー

今回紹介する内容について前提とする環境は以下のとおりです。

他のラズベリーパイやv1.3カメラでも実行は可能だとは思われますが、試してはいません。

  • 画像処理PCサイド
    • Python3が動作し前回までの内容が実行できるLinux環境
    • SDカードリーダーやラズベリーパイとのネットワーク接続など、ラズベリーパイとデータを交換する手段

以下の内容はラズベリーイカメラv2.1がラズベリーパイに接続され正常に動作していることを前提としていますので、まず公式ドキュメント などに従ってカメラの動作を確認下さい1

特に、GUIの左上のラズベリーパイアイコンから選択できる「設定」、「RaspberryPiの設定」、から「インターフェイス」タブで「カメラ」を有効にしておく必要がありますので、ご注意ください。

RAW画像のキャプチャー

ラズベリーパイにカメラをセッティングした状態で、以下のコマンドを実行します。

> raspistill -r -o raw_capture.jpg

-rがRAWキャプチャーを行うことをしめしています。

うまくいけば、カレントディレクトリにraw_capture.jpgという名前のJPEGファイルができているはずです。 このファイルは一見通常のJPEGファイルに見えますが、8MPのJPG画像ファイルにしては15MB前後と巨大なファイルサイズです。これは、RAWデータが埋め込みデータとしてファイル中に組み込まれているからです。

キャプチャーが成功したら、このファイルをUSBメモリーで移す、ネットワークコピー、またはRaspbianの入ったSDCARDをホストPCに読み込ませる、などの方法で、前回までの作業を行ったPCに移動します2

うまくキャプチャーできたか、JPEGファイルのJPEG画像部分を見てみましょう。通常の画像ビューワーや画像エディタで確認できます。

f:id:uzusayuu:20181014122747p:plain

あまり良い画像ではないですが、キャプチャーできていることは確認できました。

ラズベリーパイのRAW画像の抽出

今回の画像はこれまでと違い、JPEGファイルの中に埋め込まれているので、rawpy以外の方法で取り出す必要があります。 なお、この部分は今回の本筋ではないので、簡単な説明ですまします。詳しい内容は picameraのドキュメント を参照ください。

これ以降は、Jupyterなどpython3のインタラクティブ環境で作業します。

BayerデータはJPEGファイルの最後の部分に位置するので、読み出してから末尾だけ取り出します。

with open("raw_capture.jpg", "rb") as input_file:
    data = input_file.read()
data = data[-10237440:]

(未確認ですがカメラv1.3ではdata = data[-6371328:]とするとよいようです。)

これでBayer部分が読み込めたはずです。データを見てみましょう。

import numpy as np

with open("raw_capture.raw", "rb") as input_file:
    data = input_file.read()

w = 3282  # for v1.3 w = 2592
h = 2480  # 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)

最後の部分ですが、picameraのドキュメントによると、ラズベリーパイのBayerデータは4画素毎に5バイトのデータにまとまっていて、最初の4バイトがそれぞれの画像の上位8ビット、最後のバイトの2ビットずつが、各画素の下位2ビットになっているとのことです。

さて、ちゃんとデータが読めたのか見てみましょう。

from matplotlib.pyplot import imshow 

outimg = img.copy()
outimg[outimg < 0] = 0
outimg = outimg / outimg.max()
imshow(outimg, cmap='gray')

f:id:uzusayuu:20181014124822p:plain

どうやらそれらしい画像が取り出せたようです。

ラズベリーイカメラのRAW画像の現像

だいぶ手間がかかりましたが、ここから本題のラズベリーパイのRAW画像の現像に入ります。

まず前回までのスクリプトを使いたいところですが、これまで使った関数は殆どがrawpyのインスタンスに依存していました。 今回扱うラズベリーパイのRAW画像はrawpyで読み込んだものではないのでそのままでは使えません。

そこで各関数を、rawデータと指定した任意のパラメータで実行できるように書き換えて、raw_process4.pyとしました。 それぞれの関数の変更点は以下の説明で触れます。

まずはモジュールを読み込みましょう。

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_process4 as raw_process

ブラックレベル補正

ブラックレベル補正は、引数としてブラックレベルの数値を与えるようにしました。 引数で与えるブラックレベルは、左上、右上、左下、右下、の順にブラックレベルを格納したリストです。

def black_level_correction(raw_array, black_level):
    blc_raw = raw_array.copy()
    blc_raw[0::2, 0::2] -= black_level[0]
    blc_raw[0::2, 1::2] -= black_level[1]
    blc_raw[1::2, 0::2] -= black_level[2]
    blc_raw[1::2, 1::2] -= black_level[3]
    return blc_raw

ブラックレベルを64と仮定して処理してみましょう。

blacklevel = [64] * 4
blc_raw = raw_process.black_level_correction(img, blacklevel)

f:id:uzusayuu:20181014152120p:plain

ホワイトバランス補正

def white_balance_Bayer(raw_array, wbg, bayer_pattern):
    img_wb = raw_array.copy()
    img_wb[0::2, 0::2] *= wbg[bayer_pattern[0, 0]] 
    img_wb[0::2, 1::2] *= wbg[bayer_pattern[0, 1]]
    img_wb[1::2, 0::2] *= wbg[bayer_pattern[1, 0]]
    img_wb[1::2, 1::2] *= wbg[bayer_pattern[1, 1]]
    return img_wb

ホワイトバランス補正には、R・G・Bのゲインを並べたリストのwbgと、ベイヤーのパターンを渡すようになりました。

今回ホワイトバランスははっきりしないので仮に赤のゲインx1.5、青のゲインx2.2を与えます。 今回のラズベリーパイの画像では左上が青色画素ですので、bayer_patternは[[2, 1], [1, 0]]になります。

wbg = np.array([1.5, 1, 2.2, 1])
bayer_pattern = np.array([[2, 1], [1, 0]])
wb_raw = raw_process.white_balance_Bayer(blc_raw, wbg, bayer_pattern)

f:id:uzusayuu:20181014145358p:plain

デモザイク

前回のデモザイクは画像の左上隅が赤色画素であることを仮定していました。 今回のラズベリーパイの画像では左上が青色画素ですのでこのままでは実行できません。 ベイヤーパターンごとの違いは各チャンネルの位相を180度変える事になります。

def advanced_demosaic(dms_input, bayer_pattern):
    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

    # generate FIR filters to extract necessary components
    FC1 = np.matmul(vhpf, hhpf)
    FC2H = np.matmul(vlpf, hhpf)
    FC2V = np.matmul(vhpf, hlpf)
    FL = identity_filter - FC1 - FC2V - FC2H

    # f_C1 at 4 corners
    c1_mod = signal.convolve2d(dms_input, FC1, boundary='symm', mode='same')
    # f_C1^1 at wy = 0, wx = +Pi/-Pi
    c2h_mod = signal.convolve2d(dms_input, FC2H, boundary='symm', mode='same')
    # f_C1^1 at wy = +Pi/-Pi, wx = 0
    c2v_mod = signal.convolve2d(dms_input, FC2V, boundary='symm', mode='same')
    # f_L at center
    f_L = signal.convolve2d(dms_input, FL, boundary='symm', mode='same')

    # Move c1 to the center by shifting by Pi in both x and y direction
    # f_c1 = c1 * (-1)^x * (-1)^y
    f_c1 = c1_mod.copy()
    f_c1[:, 1::2] *= -1
    f_c1[1::2, :] *= -1
    if bayer_pattern[0, 0] == 1 or bayer_pattern[0, 0] == 3:
        f_c1 *= -1
    # Move c2a to the center by shifting by Pi in x direction, same for c2b in y direction
    c2h = c2h_mod.copy()
    c2h[:, 1::2] *= -1
    if bayer_pattern[0, 0] == 2 or bayer_pattern[1, 0] == 2:
        c2h *= -1
    c2v = c2v_mod.copy()
    c2v[1::2, :] *= -1
    if bayer_pattern[0, 0] == 2 or bayer_pattern[0, 1] == 2:
        c2v *= -1
    # f_c2 = (c2v_mod * x_mod + c2h_mod * y_mod) / 2
    f_c2 = (c2v + c2h) / 2

    # generate RGB channel using 
    # [R, G, B] = [[1, 1, 2], [1, -1, 0], [1, 1, - 2]] x [L, C1, C2]
    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

    return dms_img

実行してみます。

dms_img = raw_process.advanced_demosaic(wb_raw, bayer_pattern)

f:id:uzusayuu:20181014145414p:plain

カラーマトリクス補正

カラーマトリクスの値もはっきりしません。今回は色が多少強調されるように、以下のようなマトリクスをかけ合わせてみます。

[[1536, -256, -256], 
 [-256, 1536, -256],
 [-256, -256, 1536]]

実行します。

img_ccm = raw_process.color_correction_matrix(dms_img, [1536, -256, -256, -256, 1536, -256, -256, -256, 1536])

f:id:uzusayuu:20181014145426p:plain

ガンマ補正

最後にガンマ補正です。ガンマ値を引数でもたせるようにしました。

def gamma_correction(rgb_array, gamma):
    img_gamma = rgb_array.copy()
    img_gamma[img_gamma < 0] = 0
    img_gamma = img_gamma / img_gamma.max()
    img_gamma = np.power(img_gamma, 1/gamma)
    return img_gamm

実行します。

img_gamma = raw_process.gamma_correction(img_ccm, 2.2)

画像を保存して確認してみましょう。

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

import imageio
from pylab import imshow, show
imshow(imageio.imread('raspi_raw_out.png'))
show()

f:id:uzusayuu:20181014140022p:plain

お世辞にもきれいな画像とは言えませんがどうにかラズベリーパイでキャプチャーしたRAWデータを現像して画像ファイルにすることができました。

まとめ

Raspberry Pi 3BでキャプチャーしたRAW画像を現像して画像ファイルに変換しました。 今回はどうにか現像するところで終わりでしたが、次回以降画質に手をいれて行きたいと思います。

最後に

今回の内容はRAW画像データ及びraw_process4.pyと共にgithubにアップロードしてあります。

github.com


  1. 公式ドキュメント以外では以下の記事が参考になります。Raspberry PiカメラモジュールV2でデジカメを作ってみた

  2. 前回までの環境をRaspberry Pi上で構築できていれば同じことができるはずですが、私の環境ではうまくいきませんでしたので、別のPCにデータを運んでから以下の作業を行っています。

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

はじめに

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

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

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

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

モジュールの高速化

本題に入る前に、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, "output4.png")

画像の確認

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

f:id:uzusayuu:20181008024743p:plain

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

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

f:id:uzusayuu:20181008025157p:plain

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

まとめ

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


  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

次回はもう少し複雑なデモザイクを紹介しようと思います。

ゼロから作る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フォーマット画像データから、普通の画像ビューワーなどで表示できる画像ファイルをできるだけスクラッチから作成して見たいと思います。

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が使えれば以下で解説する内容を実行するのが楽になると思います。実行例はすべてJupyter notebook上でのものです。

なお処理自体に関係ありませんが、いくつかのパラメータを取得するのに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):
        colors = [0, 0, 0, 0]
        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画素ある緑は平均値をとります。

dms_img = np.zeros((h//2, w//2, 3))
for y in range(0, h, 2):
    for x in range(0, w, 2):
        colors = [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


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

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

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

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