Moiz's journal

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

「低レイヤを知りたい人のための Cコンパイラ作成入門」を読んで、オレオレコンパイラを作り始めた

はじめに

表題通りですが、Rui Ueyama氏の「低レイヤを知りたい人のための Cコンパイラ作成入門」というPDF本を読んでCコンパイラ的なものを作り始めましたよ、というダラっとした内容です。

発端

ツイッター等で知ったのですが、セキュリティキャンプというイベントで、Cコンパイラを書く、というコースがあったそうです。

www.ipa.go.jp

内容については当時各所で話題になりましたし、参加者のブログなども多数書かれているようです。 たとえばこちらなど、コースの内容や雰囲気がよくわかります。

0x19f.hatenablog.com

hsjoihs.hatenablog.com

こういったものを見て、コンパイラ作成に興味を持ち始めたところ、そのセキュリティキャンプ講師のRui Ueyama氏が「低レイヤを知りたい人のための Cコンパイラ作成入門」を書かれた事を知りました。

booth.pm

PDF本は有料ですが、同内容のものがウェブページでも無料公開されています1

低レイヤを知りたい人のための Cコンパイラ作成入門

これはやるしかないだろう、という事で自分もオレオレコンパイラを書き始めました。

なお、目標はCに近い言語のコンパイラを作って見ることでコンパイラがどのようにできているか学ぶ、ということなので、Cの規格に正しく従っていない部分が多々あります。

とりあえず書いてみる

自分はコンパイラについての知識はまったくなかったのですが、PDF本の内容が懇切丁寧なので、その内容を実行するにあたってはあまり問題なくすすめる事ができました。

内容的にはまず、数字を一つ読み込んで出力するだけという、考えられる限りでもっとも単純なプログラムをアセンブラコンパイルできるだけの「コンパイラ」を作成します。次にそこに、四則演算を導入して、変数を導入し、と進んでいきます。

執筆途中とはいえ、現在公開されている記事の最後まですすむと、こんな内容のコードがコンパイルできるようになります。

a = 4;
b = 5 * 2;
a + b;

これをコンパイルして、アセンブラにかけると、示されたとおりの計算を行って、14という返り値を返す実行ファイルを生成することができます。

本の内容が丁寧なので、最初にまったく何もなかったところからここまで、意外なほと短時間でたどりつけました。

別のコンパイラ入門書をよんでみる

こうやってコンパイラを書きすすめてきたわけですが、残念ながらPDF本の内容はまだここまでしかできていません。

そこで、まず他の本を一冊読んでみることにしました。

日本語のコンパイラ作成入門の電子書籍を探していたところ、こちらの本に出くわしました。

明快入門 コンパイラ・インタプリタ開発 (林晴比古実用マスターシリーズ)

明快入門 コンパイラ・インタプリタ開発 (林晴比古実用マスターシリーズ)

この本の内容は、スタックマシンをエミュレートしたVM上で動作するCサブセットのコンパイラソースコードの解説と、その関連技術の説明です。

とりあえずこちらの内容を一通り読んで、コードを手入力して実行させてみました2

スタックマシンのエミュレータ、という点が独特ですが、最初に読んだ「低レイヤ〜」でもスタックマシン的なコードを出力していたので、わりとすんなりなじめました。

次にどうする?

ここまで来たところで元の自作コンパイラにもどったのですが、次に何をするのか途方にくれてしまいました。

そこで、セキュリティキャンプに参加した方から、Rui氏のブログでセキュリティキャンプではどのような順に実装を進めたのか書いてあると教えていただきました。

note.mu

とりあえず、こちらの方針に従って以下のように進めていくことにしました。

  1. 引数なしの関数呼び出し
  2. 引数ありの関数呼び出し
  3. 引数なしの関数定義
  4. 引数ありの関数定義
  5. グローバル変数
  6. 配列
  7. ポインタ
  8. char型
  9. 文字列リテラル
  10. 構造体
  11. #include

書くぞ、書くぞ、書くぞ

その後はひたすら実装テスト実装テスト実装テストの繰り返しです。

毎回機能を追加する度にテストを追加し、ちょっとずつちょっとずつできることが増えて、とうとう文字列リテラルのところまですすみました。

また、セキュリティキャンプの参加者の方から(ときには講師の方から)いろいろと実践的なアドバイスをいただけたのは非常に助かりました。 特に、X86-64のABIではスタックポインタは16バイトアラインされている必用がある、というのと、rbxなどはcallee save、というのは自分ではきっと絶対に気が付かないではまっていたと思います。

そんなみなさんの協力もあって、今ではこんなコードをコンパイルすることができるようになりました。

void print(int a) {
  int b, c;
  b = 1;
  while(a >= 10 * b) {
    b = b * 10;
  }
  while(b > 0) {
    putchar((c = a/ b) + 48);
    a = a - c * b;
    b = b / 10;
  }
}

void print_string(char *a) {
  while(*a) {
    putchar(*a);
    a = a + 1;
  }
}

int fib(int a) {
  if (a <= 1) {
    return a;
  }
  return fib(a-1) + fib(a-2);
}

int main() {
  int i;
  print_string("Fibonacci seriese: ");
  for(i = 0; i < 20; i = i + 1) {
    print(fib(i));
    putchar(32);
  }
  putchar(10);
  return 0;
}

このコードをコンパイルすると以下のような結果が得られます。(m99ccというのが今回作ったコンパイラ、fibonacciが上記のプログラム。)

$ ./m99cc program/fibonacci.c > tmp.s && gcc -o tmp tmp.s && ./tmp
Fibonacci seriese: 0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181

まだまだ不完全ではありますが、ゼロから初めて2ヶ月にしてはそれなりに形になっているのではないでしょうか。改めてご協力いただいた方々に感謝します。

(まだ終わってませんが。)

リポジトリ

これまでの結果はこちらのGithubリポジトリーにアップロードしてあります。

github.com

便利な資料

コンパイラ作成中に、ここまでで紹介した以外にもいろいろ便利なリソースについて教えていただいたので、ここにまとめておきます。

セキュリティキャンプの参加者のリポジトリ

github.com

github.com

こちらのShinya Kato氏(id:mizu0x19f)とはすじょい氏(id:hsjoihs)のお二人からは別途何度も有益なアドバイスをいただきました。

オンラインコンパイラGodbolt

godbolt.org

ブラウザ上でCやその他コードをコンパイルしてアセンブラの結果が見られる便利なサイトです。

変わった名前だなと思ったら本名のようです。まさに神ツール。

セキュリティキャンプ講師のスライド

docs.google.com

セキュリティキャンプの講師をされたhikalium( id:hikalium )氏のスライドです。

hikalium氏からは他にもいろいろ教えていただきました。上記のGodboltもhikalium氏経由の情報です。

C言語規格のドラフト

http://www.iso-9899.info/n1570.html

書籍 Compilers: Principles, Techniques, and Tools 2nd By Alfred V. Aho

入門書より範囲の広い本を、と思って買ってみましたがまだ読んでいません。

Intel® 64 and IA-32 Architectures Software Developer Manuals

Intel® 64 and IA-32 Architectures Software Developer Manuals | Intel® Software

インテルCPUのマニュアルです。 読まなきゃいけないのはわかるのですが、私は個人的な事情でこの書類を見ると鬱に...。


  1. 自分は応援も兼ねてPDFの方を購入しましたが、実際に参照したのはウェブページの方です。なおまだ執筆途中ということで、掲載されている内容はまだ途中までだそうです。

  2. 実はソースコードがダウンロードできるので、こんな事をする必要はなかったのですが、ダウンロード用URLが奥付に書いてあるのに気が付かず手で入力してしまいました。

EDA Playgroundで遊んで見る

はじめに

ブラウザ上でデジタル回路設計と検証が行えるウェブページEDA Playgroundを紹介します。

https://www.edaplayground.com/

f:id:uzusayuu:20181217115012p:plain
EDA Playgroundの例の一つ(カウンター)

書いている人とEDA PlaygroundおよびDoulosとの間にはなんの関係もありません。

EDA Playgroundとは

EDA Playgroundとは、VictorEDAにより開発され、半導体回路の設計・検証のトレーニングを行うDoulosによって提供されるWebサービスです。 ブラウザ上で各種半導体デジタル設計言語による設計のシミュレーションや検証が行なえる学習用のツールです。 サポートする言語としては、Verilog HDL, System Verilog, VHDLなどがあります。

制限事項としては、

  • 一回の実行時間は60秒まで
  • モリーの使用量は100MBまで

という物があります。本格的な設計をするには心もとないですが、もともと学習用でありそういう用途に使うツールではないと思います。逆に学習用と考えると広い使いみちが考えられます。

また、特筆すべき点としてサポートする言語機能の範囲があります。 たとえばAltera製FPGAに付属するModelSim-Altera Starter EditionではSystem Verilogアサーション機能やパラメータ乱数化機能はライセンスの関係で使えません(私が使用した時点での制限であり、現在では変わっている可能性があります。) しかし、EDA Playgroundではこういった高度な機能も使用可能になっています。

EDA Playgroundでのシミュレーション

実は、Youtubeに開発者による紹介ビデオが上がっています。

www.youtube.com

ここまで説明してきてなんですが、EDA Playgroundの使い方はこのビデオを見ていただくのが一番だと思います。

そうは言っても、それだけでは寂しいのでシミュレーションの画像を紹介します。

まず、EDA Playgroundのページをブラウザで開きます。

f:id:uzusayuu:20181217122739p:plain
EDA Playgroundにアクセスした直後

まずは右上のLoginからログインしておきましょう。詳細は省きますが、EDA Playgroundのアカウントを作らなくても、Gmailfacebookのアカウントでログインできます。

次に、言語とデザインを設定します。今回は言語はデフォルトのSystem Verilog/Verilogのままにしておきます。 デザインとしては例から選んでみましょう。まず左側のコラムのExamplesからVerilog/System Verilogを選択し、出てきたリストからD flip-flopを選択します。

これでデザインとテストベンチのウィンドウが埋まるはずです。

f:id:uzusayuu:20181217123206p:plain
D-Flip Flopの例をロードしたところ。

ではシミュレーションしてみましょう。 まずあとからウェーブフォームを確認できるように、左側の"Open EPWave after run"のチェックボックスをクリックして選択しておきましょう。

f:id:uzusayuu:20181217123714p:plain
Open EPWave after runをチェック

左上のRUNボタンを押すとシミュレーションが始まります。

f:id:uzusayuu:20181217124040p:plain
シミュレーションの実行

シミュレーションが終わると、波形が表示されます。

f:id:uzusayuu:20181217123835p:plain
EPWaveによる波形の表示

以上です。簡単ですね。

最後に

半導体回路の設計を学ぼうとすると、言語や設計技法の習得以前に、コンパイラやシミュレーターや波形ビューワーなどツールの準備に手間がかかってしまいます。 開発ツールの導入は、単に技術的な点だけでも難易度が高い上に、ツールの入手も多くは有料、ないしは基板とのバンドルなどになります。 また、それでさえもライセンスによっては使用できる機能に制限があり、学びたい機能が使えないという事もよくあります。

最近ではオープンソースのツールも増えてきましたが、それらも導入が簡単とは言えません。

そういった状況が、回路設計を新たに学ぼうという方々に対する参入障壁になっているようです。

それに対してEDA Playgroundはブラウザ一個があれば誰でもデジタル回路設計のシミュレーションが始められ、機能的な制限も非常に小さい(規模的な制限あり)、という非常にありがたいものになっています。

その割に日本での知名度がそれほど高くないのが残念だったので、1ファンとして簡単ながら紹介記事を書かせていただいた次第です。

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

このブログ記事「ゼロから作るRAW現像」を大きく再構成してより読みやすくした書籍「PythonとColabでできる-ゼロから作るRAW現像」を【技術書典6】にて頒布しました。

現在はBOOTHにて入手可能です。書籍+PDF版は2200円プラス送料、PDF版は1200円です。

moiz.booth.pm

はじめに

この記事は「ゼロから作る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現像」を大きく再構成してより読みやすくした書籍「PythonとColabでできる-ゼロから作るRAW現像」を【技術書典6】にて頒布しました。

現在はBOOTHにて入手可能です。書籍+PDF版は2200円プラス送料、PDF版は1200円です。

moiz.booth.pm

ゼロから作る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現像」を大きく再構成してより読みやすくした書籍「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- ノイズモデルとノイズフィルター」

残りの処理について

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

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

準備

今回の処理の内容は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)も調整しなくてはならないのですが今回は省きます。