Moving

本当は4/1に書く予定でしたが、そういえば4/1はエイプリルフールでしたので今日書きます。

えーと、本日、東北大学理学部物理学科を卒業し、明日から東北大学大学院情報科学研究科に入学します。所属する研究室はコンピュータビジョンなど、主に視覚的な情報について研究している研究室ですので、たぶんこれから書く記事はそれ系統の内容が多くなると思います。このブログのタイトルに少し近づきました。

元々このブログは映像を作り、その技術的な情報についてログを残せるよう立ち上げたブログでした。確かに現在は映像の「え」の字もないブログに仕上がっていますが、もちろん今でも映像は好きですし、好きなプロダクションさんの映像はいつも観に行くようにしています。そして全く役に立たなかったかと言われればそういうわけではなくて、個人的に仕事を頼まれることもありましたし、Webサイトや印刷物のデザインをする際にも役立ちましたし、ここで得た知見のお陰で様々な方とお話することもできました。

物理学科に入ったのも元々は相対論や量子力学の深遠な世界について知りたいために入りました。正直に言って完全に理解できたかというと怪しいところでありますし、理学から工学という他分野に移るわけで、一見すると今までの知識が役立たなくなるように見えます。が、実はそうでもなく、統計力学(*1)や相対論で使っている手法をコンピュータビジョンに輸入している論文があったり、案外様々な場面で役立ったりしているわけで(*2)。

全然違う領域に見えても絡んでくる領域というものは少なからずあるものです。ということで、映像と物理がお互いに絡み合っている、コンピュータビジョンの道に進むことを決めました。

これからもどうぞよろしくお願いします。

p.s. 家族は無事でした、よかったよかった

*1 ここで宣伝ですが、田崎先生の『統計力学I, II』はお薦めです。この本はぜひ理学系だけでなく工学系の方々にも読んでもらいたい本です。値段も2冊併せて7,000円くらいと安い。新品のゲームソフト1本買うくらいならこれを買いましょう
*2 もちろん線形代数や基本的な偏微分も数式や概念を学んだりする上で役立ちました。一番学んで良かったのはブラケット記法とヒルベルト空間、あれは素晴らしい

自分は無事です

自分は出張で京都にいましたので、無事です。もし身内の方が見ていましたらご安心ください。
ただ暫くは京都に滞在することになると思います。早く交通機関が復旧してほしいです。

ただ家族がどうなっているのか全く分からない。怖いです、不安です。

[文句]
iPhoneではSoftbankの災害用伝言ダイヤルに『登録できない』！！！

明らかにsoftbank側の怠慢、糞だ。糞すぎる

OMakeのContributionsに載ったった

The OMake build system: User Contributions

やったね！

余談
翻訳をすることは勉強にもなって非常にいいですし皆さんにもぜひやって欲しいのですが、
その分、なんというか、日本の中へと篭ってしまうことにも繋がっている気がします。

英語に臆することなく積極的にアウトプットできるよう、自分自身の心構えをもう少しきちんとしなければならんですね。

boost.GILで組む、ジェネリックな画像処理のアルゴリズム

"Image dans le néant" by gelinh

boost.GILは凄い！開発者の頭の良さがビシビシと伝わってくる！
ということで、今回はGILに関しての紹介記事を書こうと思います。

概要
あなたは画像処理のエキスパート。顧客の依頼で、8bitのRGB画像を処理するアルゴリズムを記述していたとします。ところが対象となるデバイスの仕様を調べていた際に、実はRGBA画像にも対応させなければいけないことが分かりました。面倒だと思いながら書いていた矢先、さらにBGRやABGR,さらには16や24bitにも対応したアルゴリズムを記述しなければならないことが判明しました。なんということでしょう…これらの画像すべてに対してアルゴリズムを書くなんて、とてもじゃないですがやってられない。やめてくれ！って感じです。

boostに付いてくるGILを用いることで、画像に対する操作をよりジェネリックに行うことが出来ます。具体的に言うと、あなたはGILの作法に従ってアルゴリズムを書き下すだけで、メモリに展開されている画像データの色空間、色深度、レイアウト、配置順に『関わらず』、任意のデータ配置に対して適用できるアルゴリズムを作ることができるのです。最初に提示した問題が一気に解決しちゃいました。

長所と短所
長所

• 画像処理に対してジェネリックなアルゴリズムを記述でき、少ない記述量で全てをカバーできる。処理速度も普通に書き下した場合と同等(静的なポリモーフィズムしか使っていないので、コンパイル時に必要な計算はすべて終わっている)。
• 様々なアダプタが存在している(e.g. subimage_view, rotated180_view)ので、ばんばん遅延評価して一気に処理を行うことが可能。
• 普通の人が普通に書くよりも、GILのアルゴリズムに従って書いたほうが大抵の場合早い
• ソースコードが比較的読みやすい。コード量も短くなる。
• 誰が作ったかわかんないもの利用したくねーという方も安心のAdobe製。画像処理の専門家が作ってるのでとても合理的。

短所

• テンプレートを多用しており『非常に難解』。
• 全てを抽象化しているために『非常に難解』。
• C++と画像処理の両方に秀でていなければ記述できないので、使っている人が殆どいない。
• 英語のサイトですら文献が殆どない

正直なところ、短所の部分がかなり厳しい。かなりC++をやりこんでいないと使えないということは、それだけでユーザ層を大分狭めます。仕方がないことですけど。

Viewの概要
GILでは、画像に対する操作はすべて『View』を起点にして行われます。Viewは以下のような構造を持ち、ジェネリックな操作を行えるよう工夫してあります。


View
|-- Locator
|    `-- Pixel
|         |-- Channel
|         `-- Layout
|              |-- Color Space
|              `-- Channel Mapping
`-- dimensions

(上図はMemory basedで、画素のメモリ配置がInterleavedな場合のView構造)

以下のリストで、それぞれの名称がどのような役割を持っているのか説明します。

• View: (x, y)座標における画素値の他、特定のX、Y方向におけるイテレータや全ピクセルにわたって捜査できるxyイテレータなど、「様々な手法によって画素にアクセスできる手段」を提供する、一種の抽象的なクラスを表しています。画像の幅、高さを表すdimensionsと、後述するLocatorを持ちます。(注: 画像のデータを実際に保持している『わけではありません』。Viewはあくまで、対象の画素値へとアクセスするための『手段』を提供しているだけに過ぎないのです)
• Locator: 画像の2次元的な『位置』を表します。ランダムアクセスイテレータの2次元版のようなものと考えれば理解が早いと思います。具体的に言いますと、 *loc あるいは *(loc + point2(dx, dy)) などで現在の画素値や(dx, dy)だけ離れた位置の画素値を取得することは出来ますが、Locatorの現在地(xy座標)や画像の大きさ(width, height)を取得することはできません。
• Pixel: Locatorによって返される画素値を表します。後述するChannelとLayoutを持ちます。
• Channel: 画素の色深度を表します(8bit, 16bit, etc…)。各型の取りうる最大値、最小値も取得できます。
• Layout: 色空間のレイアウトを表します。ちょっと分かりにくいので、下の2つを読んだほうが早いでしょう。
• Color Space: 色空間を表します(e.g. RGB, CMYK, HSV, Gray, etc…)。
• Channel Mapping: 実メモリ上での配置順を表します(e.g. RGB, BGR, etc…)。ただし、この情報はInterleave画像のみ用います(Planar画像だと配置順は関係ないため)。

アルゴリズムを用いて書いてみる
それでは実際にGILに付属しているアルゴリズムを用いて実際に書き下してみます。以下の2つのアルゴリズムを用いてみましょう。

boost::gil::transform_pixel(const View& src, const View& dst, F functor)
対象のピクセルを受け取って別のピクセルを返す関数を、全ピクセルに対して適用します。この関数は、画素値と出力値が1対1対応である場合に用います(e.g. 明度、コントラスト, レベル補正, トーンカーブ, etc...)。

boost::gil::transform_pixel_positions(const View& src, const View& dst, F functor)
周囲のピクセルの値を元に結果となるピクセルを返す関数を、全ピクセルに対して適用します。transform_pixelとは、受け取る引数がロケータであるという点が異なっており、周囲の画素値を元に実際の出力値を求める場合に有効です(e.g. Blur, Sovel, Laplacian, etc...)。

実際にやってみれば、どんな感じか分かるでしょう:

#include <boost/gil/gil_all.hpp>
#include <boost/gil/extension/io/jpeg_io.hpp>
#include <opencv2/imgproc/imgproc.hpp>

using namespace boost::gil;

struct change_brightness {
public:
    change_brightness(): weight_(1.0) {};
    explicit change_brightness(float weight): weight_(weight) {};

    template<typename Pixel>
    Pixel operator()(Pixel src) {
        typedef typename channel_type<Pixel>::type channel_type;
        Pixel dst;
        for(int c=0; c<num_channels<Pixel>::value; ++c) {
             dst[c] = cv::saturate_cast<channel_type>(
                     static_cast<float>(src[c] * weight_));
        }
        return dst;
    }
private:
    float weight_;
};

struct change_brightness_locator {
public:
    change_brightness_locator(): weight_(1.0) {};
    explicit change_brightness_locator(float weight): weight_(weight) {};

    template<typename Locator>
    typename Locator::value_type operator()(Locator loc) {
        typedef typename Locator::value_type pixel_type;
        typedef typename channel_type<Locator>::type channel_type;
        pixel_type src = loc(0, 0);
        pixel_type dst;
        for(int c=0; c<num_channels<Locator>::value; ++c) {
            dst[c] = cv::saturate_cast<channel_type>(
                    static_cast<float>(src[c]) * weight_);
        }
        return dst;
    }
private:
    float weight_;
};

template<typename Locator>
struct x_gradient_functor {
public:
    explicit x_gradient_functor(Locator loc)
    : left_(loc.cache_location(-1, 0)),
      right_(loc.cache_location(1, 0)) {};

    typename Locator::value_type operator()(Locator loc) {
        typedef typename Locator::value_type pixel_type;
        typedef typename channel_type<Locator>::type channel_type;
        pixel_type src_left = loc[left_];
        pixel_type src_right = loc[right_];
        pixel_type dst;
        for(int c=0; c<num_channels<Locator>::value; ++c)
            dst[c] = (src_left[c] - src_right[c])/2;
        return dst;
    }

private:
    typename Locator::cached_location_t left_;
    typename Locator::cached_location_t right_;
};

template <typename SrcView, typename DstView>
void x_gradient(const SrcView& src, const DstView& dst) {
    assert(src.dimensions() == dst.dimensions());
    transform_pixel_positions(
            subimage_view(src, 1, 0, src.width()-2, src.height()),
            subimage_view(dst, 1, 0, src.width()-2, src.height()),
            x_gradient_functor<typename SrcView::locator>(src.xy_at(0, 0)));
}

int main() {
    const char* filename = "lena.jpg";
    rgb8_image_t src_image;

    jpeg_read_image(filename, src_image);
    rgb8_image_t dst_image(src_image.dimensions());

    transform_pixels(
            const_view(src_image),
            view(dst_image),
            change_brightness(1.5));
    jpeg_write_view("result_normal.jpg", view(dst_image));

    transform_pixel_positions(
            const_view(src_image),
            view(dst_image),
            change_brightness_locator(2.0));
    jpeg_write_view("result_locator.jpg", view(dst_image));

    x_gradient(const_view(src_image), view(dst_image));
    jpeg_write_view("result_gradient.jpg", view(dst_image));

    return 0;
}

また、画像の出力例を下記に示します:

左から、元画像, change_brightness, change_brightness_locator, x_gradientの結果

change_brightness_locator()については、for_each_pixel用のアルゴリズムをlocatorを用いて適用させた場合にどうなるか確認するためだけに実装しましたので、実務的ではありません。また、値を型の持つ最大、最小値に丸めるためにOpenCVのcv::saturate_castを使用しています。

main()についてはどんなことをやってるのか大体感覚的に理解できるはずですので、とりあえず3つのファンクタに関して説明を加えていきます。

typedef typename channel_type<Pixel>::type channel_type;

channel_type<T>::typeで実際に使われているチャネルの型を取得できます。Tにはピクセル、ロケータ、ビューのどれを指定しても構いません。今回の場合ですと8bitのチャネルを指定しているので、型にはunsigned charが入ります。

num_channels<Pixel>::value

ピクセルを構成しているチャネルの数を取得できます。この場合ですと色空間がRGBなので値には3が入ります。この値はコンパイル時に決定されるので、最適化を用いることでループ展開が行われ、速度的なオーバーヘッドはありません。

typedef typename Locator::value_type pixel_type;

ロケータが返すピクセルの型を取得できます。この場合ですとrgb8_pixel_tが入ります。

pixel_type src = loc(0, 0);

()演算子によって現在のピクセル値を取得しています。他にも*演算子でイテレータっぽくアクセスすることもできます。

explicit x_gradient_functor(Locator loc)
    : left_(loc.cache_location(-1, 0)),
      right_(loc.cache_location(1, 0)) {};

ロケータの変動値をキャッシュしておくことで、高速化を図ります。Locator::cached_location_tによってロケータのキャッシュ型を取得します。キャッシュを使ったアクセスには[]演算子を用います。
仰々しく言っていますが、要はstd::ptrdiff_tみたいなものです。移動する値は同じなんだから記憶しておきましょうということです。

subimage_view(src, 1, 0, src.width()-2, src.height())

subimage_view()によって画像の領域を切り出しています。これはアダプタで、実際の評価は遅延的に行われます。

まとめ
大体こんな感じでアルゴリズムを組んでいくのがGIL流です。遅延評価とアルゴリズムを武器に、ジェネリックなアルゴリズムを作っていきましょう。

余談
C++に関しての記事を投稿するのは正直に言って『少々怖い』です。大筋では合っていると思いますが、完全に合っているという自信があまりないからです。ですが、GILに関しては本当に、本当に参考となる文献が少ない。サンプルコードも多くない。こういった中でこの記事が少しでも参考程度になっていただければ幸いです。

それと、ざっくり組んでしまったのでこうしたほうがいいよとか、間違い等あれば遠慮なく言ってもらえると助かります。

(2/16) 追記: ソースコード若干の修正、説明や画像例の追加。
(3/31) typo: e.q. -> e.g.

画像をぼかす: いろんなボケ、リアルなボケ

『画像をぼかす』と一口に言っても、その方法にはいろんなものがあります。一番有名なのはPhotoshopの『ガウスぼかし』に見られるようなGaussian Blurですが、残念ながらこれはカメラに見られるようなリアルなぼかしと比べるとどうしても見劣りします。原因はいくつか考えられて、一つはGaussian Blurに使われているガウス関数の勾配が穏やかなため、現実のレンズのシチュエーションを考慮していないという点があります。

それじゃあ勾配を急にした関数ー例えばフェルミ分布関数のようなーをカーネルに採用すれば良いのかというとそういうわけではありません。有名なLena画像で試してみましょう。

左が元画像、中央がカーネルにガウス関数を用いた画像、右がフェルミ分布関数を用いた画像。
確かに中央よりは綺麗だけど、それでもまだリアルとは言えない

それじゃあ他に何が違うのでしょう？端的に言うと、黒も白も等価に扱って平均してしまうのがマズい。

画像をぼかすという行為は要は畳込み積分で、周囲のピクセル値を特定のウェイトで平均してやってることになります。で、普通にブラーしてしまうとピクセルの輝度に関わらず常に一定の割合で平均してしまいますので、なんか濁ったような画像が生成されてしまうと。

通常のブラー。カーネルにはよりリアルに見せるため正六角形を用いている

それではどのようなアプローチを用いればいいのでしょうか？レンズブラーを観察してみると、明るい箇所が強調されてぼかされていることが分かります。言い換えると、輝度の高いピクセルが優先的に混合されている、つまり予めRGBの全ピクセル値を適当な単調増加関数で写像してから畳込み積分を行って、その後対応する逆関数で戻してやればいい。要はhirax.netとやっていることは同じで、expしてからボカしてlogすると、こういうわけです。

結果の画像。比較的リアルにボカされていることが分かる

並べてみれば、違いがはっきり(クリックで拡大)

Lenaさんだけでは少し分かりづらいので、適当な背景写真に適用してみましょう。

森林画像

まず、通常の手法でぼかしてしまうとどうなるでしょう？

森林画像(通常のブラー)

お世辞にも綺麗なボケとは言えません。それじゃあ今回の手法でぼかしてみます。

森林画像(リアルブラー)

なんということでしょう…リアルなボケ画像ができました。

比較画像(クリックで拡大)

今回組んだOpenCVのソースコードを下記に示します。2.2でコンパイルしましたが、多分2.0以上だったら普通に通ると思います。

#include <algorithm>
#include <cmath>
#include <string>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/core/operations.hpp>

struct exp_each_element
{
public:
    exp_each_element(float factor): factor_(factor/256.0){};
    cv::Vec3f operator()(cv::Vec3b &color) {
        return cv::Vec3f(
                std::exp(static_cast<float>(color[0])*factor_),
                std::exp(static_cast<float>(color[1])*factor_),
                std::exp(static_cast<float>(color[2])*factor_));
    }
private:
    const float factor_;
};

struct log_each_element
{
public:
    log_each_element(float factor): factor_(factor/256.0){};
    cv::Vec3b operator()(cv::Vec3f &color) {
        return cv::Vec3b(
                cv::saturate_cast<uchar>(std::log(color[0])/factor_),
                cv::saturate_cast<uchar>(std::log(color[1])/factor_),
                cv::saturate_cast<uchar>(std::log(color[2])/factor_));
    }
private:
    const float factor_;
};

cv::Mat make_fermi_kernel(float radius, float range) {
    const int size = static_cast<int>(2.0*range + radius)*2 + 5;
    const float mu = static_cast<float>(size/2);
    const float range_inversed = 1.0 / range;
    cv::Mat dst(size, size, CV_32F);
    for(int i=0; i<size; ++i) {
        for(int j=0; j<size; ++j) {
            const float u = static_cast<float>(i - mu);
            const float v = static_cast<float>(j - mu);
            const float r = std::sqrt(u*u + v*v);
            dst.at<float>(i, j) =
                    1.0/(std::exp(range_inversed*(r - radius)) + 1.0);
        }
    }
    cv::normalize(dst, dst, 1.0, 0, CV_L1);
    return dst;
}

cv::Mat make_regular_polygon_kernel(int number, float radius) {
    const int size = static_cast<int>(2.0*radius) + 1;
    const float center = static_cast<float>(size)/2.0;
    std::vector<cv::Point> dst_vector;
    for(int i=0; i<number; ++i) {
        const float theta = 2.0*M_PI*i/number;
        const cv::Point point = cv::Point(
                center + radius*std::cos(theta),
                center + radius*std::sin(theta));
        dst_vector.push_back(point);
    }
    cv::Mat dst(size, size, CV_32F);
    const cv::Scalar color(1.0, 1.0, 1.0, 1.0);
    cv::fillConvexPoly(dst, &dst_vector[0], number, color);
    cv::normalize(dst, dst, 1.0, 0, CV_L1);
    return dst;
}

int main(int argc, char** argv) {
    const std::string source = argc >= 2 ? argv[1] : "lena_std.tif";
    const std::string destination = argc >= 3 ? argv[2] : "out.jpg";
    const float factor = 6.0;

    cv::Mat image = cv::imread(source);
    CV_Assert(image.type() == CV_8UC3);
    cv::Mat kernel = make_regular_polygon_kernel(6, 8.0);

    cv::Mat exp_image(image.size(), CV_32FC3);
    std::transform(image.begin<cv::Vec3b>(), image.end<cv::Vec3b>(),
            exp_image.begin<cv::Vec3f>(), exp_each_element(factor));

    cv::filter2D(exp_image, exp_image, -1, kernel);

    std::transform(exp_image.begin<cv::Vec3f>(), exp_image.end<cv::Vec3f>(),
            image.begin<cv::Vec3b>(), log_each_element(factor));

    cv::imwrite(destination, image);
    return 0;
}

余談

• 実質80行くらいで完成したので嬉しいです。可読性も多分結構読みやすい部類だと思う
• C++のコードが何やってるか分からない方は、まず関数オブジェクトとSTLについて調べてみましょう。あれはいいもんです

追記(11/03/23)

どんなデバイスであるか、たとえば、それが入力デバイスであるのか、出力でバイスであるのか、はたまた、計算処理デバイスであるのか次第で「輝度値の単位系」をどのように保持（処理）するべきかは違うことでしょう。
大切なことは、取り扱う「値」がどんな単位であるのかを見失わないことなのか、と思います。 

hirax.net::「コンピュータ画像処理」と「輝度値の単位系」 

自分がこのような『任意の単調増加関数で良い』という具合で筆を進めたのには2つ理由があります。

まず1つ目が、『全てのデジタルカメラにおいて、ピクセルは輝度の対数の次元を取っているのだろうか？』という単純な疑問です。

きちんと物理的/工学的な意味合いを持つ画像処理にするためには、格納されているピクセルデータの次元がきちんと揃えられている必要があります。自分自身も普通に考えて輝度の対数が格納されているんだろうということで、任意ではなく指数関数に制限された関数系を使いましょうと運びたかったのですが、そうするための証拠が不足していた。

もしかしたら自分の予想もつかない方法で処理されている可能性があるわけで、そういった証拠が十分にないことから、単調増加関数と少し曖昧とした意味合いを含ませたのです。

もう一つの理由は『今回の画像処理は物理的に正確である必要がない』という理由です。これが大きかった。

そもそも今回の画像処理はどういう目的があるのでしょうか？ぼけて普通よりもリアルに見えて面白い、それで終わりです。物理的に計測をするわけではなく、ざっくり『それっぽく』見えればそれでいいという類のものです。そのような目的の下では、下手に物理的正確さを求めて処理時間を長くするよりは、計算速度の早いフェイクのほうがより適しています。

だから何も指数関数に制限する理由はあまりないわけで、より早くてよりそれっぽく見えれば、そちらのほうが優れていると、そう言えるわけです(ただ、今回の計算量からすると最も大きなオーダーとなっているのは畳み込み積分の箇所で、単調増加関数自体のコストはそれに比べれば微々たるものでしょう)。