boost.gilを利用して透過画像を作る

『対象の画像から、特定の色の箇所だけを透明にした透過画像を作る』

こういった課題に一ヶ月近くかかって出来なかったとかいう話を聞いて、流石にそれはかかり過ぎだろうと脊髄反射で言ってしまったわけですけど、それじゃあ自分だったらどうだろうとか後で思ってちょっとやってみることにしました。好き勝手言って出来なかったとか恥ずかしすぎですから、やっぱりそこらへんは実際にやってみる責任みたいなものは感じているわけです。というわけでboost.gilを使って書いてみた結果が以下。

#include <boost/gil/gil_all.hpp>
#include <boost/gil/extension/io/png_io.hpp>

using namespace boost::gil;

template<typename Pixel>
struct transparent_functor {
public:
    explicit transparent_functor(Pixel bg_pixel): bg_color_(bg_pixel) {};
    Pixel operator()(const Pixel& src) {
        if(src == bg_color_) {
            Pixel dst = src;
            get_color(dst, alpha_t()) = 0;
            return dst;
        } else {
            return src;
        }
    }
private:
    Pixel bg_color_;
};

int main() {
    const char* filename = "src.png";
    const rgba8_pixel_t background_color(0, 255, 0, 255);

    rgba8_image_t src_image;
    png_read_image(filename, src_image);
    rgba8_image_t dst_image(src_image.dimensions());

    transform_pixels(
        const_view(src_image),
        view(dst_image),
        transparent_functor<rgba8_pixel_t>(background_color));
    
    png_write_view("dst.png", view(dst_image));

    return 0;
}

結果としてはこんな感じです。上が元の画像で、下が結果の画像。正常に抜き出されていることが分かります:

余談

• プログラムより文章を紡ぎ出すのに時間がかかって、改めて自分の文章力の無さにうんざりしています
• というか最近暑い！自分は寝たいのに梅雨ですごい蒸し暑いから寝付けず気になってこんな不毛なことやっちゃってるわけです。今度こそ寝ます！！！

Thrustから見る、reduceを用いたアルゴリズム実装

"The Great Day of His Wrath" by John Martin

はじめに

世はまさに大並列時代。火を手にした人類が飛躍的な進化を遂げたように、並列化のパラダイムは、今まで到底不可能と思われていた速さで計算を行うことができるようになりました。

ところが、どんな手法にも欠点は存在するもので、実際に実装しようとすると非常に難しい。何故かというと、並列には並列化特有の問題が存在しており、愚直に実装してしまうとCPUより早くなってしまうどころか遅くなってしまうことだってあり得るのです。これを回避するにはGPUの内部構造についてきちんと理解をした上で、実装したいアルゴリズムそれぞれの場合に特化したコーディングを行う必要がある。

しかしよくよく考えるとおかしな話です。私たちが実装したいのはあくまで手法であり、ハードウェアではありません。なぜこのような詳細について把握する必要があるのでしょう。これはきちんとした抽象化がなされていない結果生み出された、言ってみれば妖です。この妖によって余計な手間が生み出され、コーディング量は増大し、余計なバグが生成され、絶望した開発者はうつになり、自殺者は増え、世界は暗雲に包まれるー

そのような混沌の中、ある一筋の光、Thrustが現れました。ThrustはCUDAを抽象化し、C++のパラダイムを用いることで本質にのみ注力できるよう開発されたライブラリです。ユーザはただIteratorに対してアルゴリズムを適用するかの如くdevice_vector<T>と各種アルゴリズムを駆使することで、流れるようなプログラムを組むことが可能となります。このような抽象を用いることで、ユーザは詳細な内部構造について頭を悩ませる必要が無くなりました。誰が書いているのか分からないから使いたくない？ThrustはNVIDIAに所属している開発者Jared HoberockとNathan Bell氏によって書かれています。2人とも並列化のプロフェッショナルなので、普通の人が愚直に実装するよりも大抵の場合効率的でしょう。

さて、Thrustには並列化を行うための多種多様なアルゴリズムが存在しています。そして、Thrustはオープンソースです。これはつまり、このソースコードを読めば並列化のプロがどのようにして各種アルゴリズムの並列化を実装しているのかが分かるということです。

Reductions

ということで、最近はThrustのソースコードをちょくちょく読んでいるところです。まだ読み始めた最中なのですが、特に感動したのはReductionsの項でした。

数々あるSTLのアルゴリズムを普通に実装しようとすると、一つ一つのアルゴリズムすべてに対して、コアレッシングなどのGPGPU特有の問題を回避したコードを記述しなければならず、非常に面倒です。Thrustはこの問題をreduce関数によって解決しました。この関数は並列化が可能なので、これを用いて書かれた関数は自動的に並列化されます。つまり、Thrustはいわば『全から一を返す』アルゴリズムすべてをreduceを用いて記述することによって、冗長な記述を回避しているのです。

このreduceを用いた並列化の実装について簡単に説明していくことが、今回の記事の目的です。reduceを初めて聞いた方は前回の記事『並列環境におけるreduceとscanアルゴリズム』を参照してください。

ここで、本来のThrustはC++によって記述しているのですが、そのままC++で書くと非常に冗長で解り辛くなってしまうので、今回は本質のみを抽出するためHaskellを用いて説明します。あまりHaskellに慣れていないので不適切な書き方等あるかもしれませんが、そのへんはご指摘頂けると幸いです。また、今回のreduceにはfoldr関数を代わりに用いています。分かると思いますが一応。

実装

• max_element, min_element, minmax_element

その名の通りリストの中の最大値、最小値、最大と最小値のタプルを返します。

*Main> max_element [1,3,2,4,3,1]
4
*Main> min_element [1,3,2,4,3,1]
1
*Main> minmax_element [1,3,2,4,3,1]
(1,4)

ここらへんはまぁ普通に分かります:

max_element :: Ord a => [a] -> a
max_element list = foldr max (head list) list

min_element :: Ord a => [a] -> a
min_element list = foldr min (head list) list

minmax_element :: Ord a => [a] -> (a, a)
minmax_element list = foldr minmax_func (first, first) $ zip list list where
    first = head list
    minmax_func (min1, max1) (min2, max2) = (min min1 min2, max max1 max2)

• all_of, any_of, none_of

all_of, any_of, none_of関数は対象のリストと真偽値を判定する関数を引数に取り、

• すべて当てはまっている
• 一つ以上当てはまっている
• すべて当てはまっていない

場合にはTrueを返す関数です。pythonではall(), any(), C++ではまんまのやつがありますね。

*Main> all_of (\x -> x > 3) [1,3,2,4,3,1]
False
*Main> any_of (\x -> x > 3) [1,3,2,4,3,1]
True
*Main> none_of (\x -> x > 3) [1,3,2,4,3,1]
False

これも実装的にはらくちん:

all_of :: (a -> Bool) -> [a] -> Bool
all_of binary_pred list = foldr (&&) True $ map binary_pred list

any_of :: (a -> Bool) -> [a] -> Bool
any_of binary_pred list = foldr (||) False $ map binary_pred list

none_of :: (a -> Bool) -> [a] -> Bool
none_of binary_pred list = all_of (not . binary_pred) list

関数合成を用いてエレガントに解決。

• inner_product, transform_reduce

inner_productは2つのリストに二項演算子を作用させてから、その結果からなるリストをreduceする関数です。transform_reduceはリストをいったんユニタリー関数を用いて変形させてから、reduceを実行する関数です。

*Main> inner_product (*) (+) 0 [0,1,2] [3,4,5]
14 -- = 0*3 + 1*4 + 2*5

*Main> transform_reduce (\x -> x*2) (+) 0 [1,2,3]
12 -- = 1*2 + 2*2 + 3*2

C++だといろいろ冗長に書かないといけないですけど、Haskellだったら余裕です:

inner_product :: (a -> b -> c) -> (c -> c -> c) -> c -> [a] -> [b] -> c
inner_product binary_op1 binary_op2 init list1 list2 =
    foldr binary_op2 init $ map operate $ zip list1 list2 where
        operate (x, y) = binary_op1 x y

transform_reduce :: (a -> b) -> (b -> b -> b) -> b -> [a] -> b
transform_reduce unary_op binary_op init list =
    foldr binary_op init $ map unary_op list

• count_if, count

count_ifは条件に合致している要素数を返す関数で、countは指定した値と同じ要素の数を返す関数です。

*Main> count_if (\x -> x < 3) [1,2,3,4,5]
2
*Main> count 3 [1,3,2,2,3]
2

count_if :: (a -> Bool) -> [a] -> Int
count_if binary_pred list = foldr (+) 0 count_list where
    count_list = map transform_func list where
        transform_func x = if binary_pred x then 1 else 0

count :: Eq a => a -> [a] -> Int
count value list = count_if (== value) list

ThrustではTrueだったら1, Falseだったら0となるリストを作成して、それを足していく手法をとっていました。
countはcount_ifが思いつけば楽勝でしょう。

• find_if, find_if_not

find_ifはリストの中で条件に合致している『一番最初の要素』を示すイテレータを返します。find_if_notは逆で、条件に合致していない『一番最初の要素』を示すイテレータを返します。今回は単純に要素のラベル番号を返すようにしました。なかったら最後のラベル+1が帰ってきます。

*Main> find_if (\x -> x > 3) [1,3,2,5,4]
3
*Main> find_if (\x -> x > 7) [1,3,2,5,4]
5
*Main> find_if_not (\x -> x > 3) [1,3,2,5,4]
0

これは微妙に入り組んでいます:

find_if :: (a -> Bool) -> [a] -> Int
find_if binary_pred list = snd result where
    result = foldr find (False, length list) tuple_list where
        find (b1, i1) (b2, i2)
            | b1 && b2  = (True, min i1 i2)
            | b1        = (b1, i1)
            | b2        = (b2, i2)
            | otherwise = (False, max i1 i2)
        tuple_list = zip (map binary_pred list) [0..]

find_if_not :: (a -> Bool) -> [a] -> Int
find_if_not binary_pred list = find_if (not . binary_pred) list

ただ、find_ifはちょっと本家と違います。インデックスが小さい項が常に左の引数となるため省略してるんだとは思いますけど、対称性を満たしていないのがちょっと気持ち悪いので、ここでは交換則を満たすように関数を変形しています。

実際の実装と課題
ここではhaskellを用いて実装しましたが、もちろん実際の実装は異なっており、zipはzip_iterator, mapはtransform_iterator, タプルはtupleを用いて実装しています。これらのイテレータの評価は遅延的に行われるので、複数のメモリアクセスが発生せず、アルゴリズムの高速化に貢献します。

これらのfancy iteratorを使用した戦略は非常に賢いとは思いますが、同時に欠点もあります。複数のイテレータを大量に組み合わせるために、通常のSTL algorithmのようにbeginとendを使用してしまうと、両方に対してfancy iteratorを組み合わせなければならないため、2倍の手間がかかってしまうのです。

これらの解決案としてはboost::rangeのようなRange構造を実装することです。これをパイプ演算子で繋げていくような形をとっていけば、明快で分かりやすく、速度も十分保証されたプログラムを記述できることでしょう。

Wiresharkを用いたHTMLの初歩的なスクレイピング

前回はぽんとスクリプトを上げて「これで画像がダウンロードできますよ」とやったわけだけど、もちろん実際にGoogle Art Projectの資料が上げられているわけでもないので、きちんとHTMLやJavascriptのソースを見ながら読解していくことになる。

簡単な動作機構だったら単にHTMLのソース見るだけで良いのだけど、Google Art Projectのような複雑な機構だとそうもいかない。行数が多いと面倒だし大抵の場合難読化されてたりで一筋縄ではいかないからだ。そこで、今回はパケット解析というもうひとつの武器を利用することで内部構造を別の側面から理解することにした。思考の流れを書き綴っていくので、HTMLのスクレイピングをやってみたい方にとって何か参考になれば幸いです。

目標: 高画質な絵画をGoogle Art Projectを利用して手に入れる

まずGoogle Art Projectのサイトへアクセス。絵画ごとに固有のURLが割り当てられており、このURLを解析すれば絵画を入手できると判断。見るからにソースコードからの解析は面倒くさそうだったので、パケット解析ソフトウェアWiresharkを用いて解析を行った。

本質のみに注力するためhttpプロトコルのみキャプチャ。適当なURLにアクセスし、ズームを1回行い、パケットキャプチャを中断。解析を開始。

画像本体らしきURLにアクセスしているパケットを発見。複数のアクセスがあるってことはたぶん画像は複数に分割されていてそれを読み込んでいるんだろう。調べると/unique-id=xi-yj-zkの形でアクセスしていたので、対応するxy座標とズーム値zによって対応する画素値にアクセスすることができると推測。試しにURLにアクセスしてみると、問題なくダウンロードができた。

また、x,yの値を変更してのアクセスを行った結果も問題なくアクセスできた。ただ、どうやら範囲外のxy値にアクセスすると404が返されるらしい。当たり前といえば当たり前。アクセス禁止などのペナルティはないっぽいので、(ちょっと礼儀悪いけど)404が返されたらそれが範囲外と解釈することにする。

一旦まとめてみる。画像は分割されており、これらのアクセスはなんかよく分からないunique-idという、おそらく絵画について表しているパラメータとx,y,z座標の2つだけでダウンロードができる。ついでに、リファラー, cookie, 認証などの面倒な作業をせずとも直接のダウンロードが許されていることも分かる。そこらへんはGoogleさん寛容なんだなぁと思い、まぁAPIを提供してるくらいだからなぁ、とかいろいろ考える。

他の画像に対してのホストを確認してみると、どうやら[lh3〜lh6]までの複数のホストへとアクセスしているらしい。このホストアクセスがランダムか、それとも固有のホストを呼び出さなければならないのか確認してみたかったので、試しにホストだけを変更して、その他のパラメータは固定したままでのアクセスを試みる。結果としては問題なくダウンロードできた。たぶん負荷を分散するという単純な理由から、ランダムにJavascript側で割り振っているんだろう。ここで初めてHTMLのソースを閲覧。"cmn/js/content.js"のソースを開き(難読化はされてないようだ)"ggpht"を検索。

リストに入れているところからすると、画像へのアクセスはたぶんこのリストの中だけに限られているんじゃないかなと思う。試しにlh7へアクセス。404が正常に返された。的中。この指針で間違いないので次へ進む。

次に気になるのはunique-idだ。一体何処から仕入れているのか知らないと全く手のつけようがない。考えられるのは2つで、

• アクセス毎に固有のunique-idを発行し、一定期間のみアクセスできるようにしている
• HTMLかどこかに参照しやすいように書かれており、それを抜き出してアクセスする

という発想だ。自分は前者を疑っていたので、それっぽいパケットがあるかどうか監視。結果としてはまったく存在していなかった。というと後者かなと思い、対象のHTMLからunique-id(今回の場合だと"vtzd432xqzOpi_3X8JAJ_buly36QKI0n9zGeeWNgBcwsYJTsAKK5mbM9Pgg")を検索。

ビンゴ！これで駒はすべて揃った。まとめに入る。対象の絵画をダウンロードするためには、

1. 取得したい絵画のアドレスへアクセス
2. <div id="microscope" data-microId="unique-id">となっている項を取得、unique-idを得る
3. http://lh3.ggpht.com/ - http://lh6.ggpht.com/のいづれかにアクセス。
   http://lh3.ggpht.com/unique-id=xi-yj-zk のような形で画像を取得
4. 取得した複数の画像を繋ぎ合わせる

後は非常に簡単でただこの単純なアルゴリズムに従って書き下していくだけ。その産物がこいつ。こんなに早く仕上がるってことは、たぶんGoogleさんは本格的にダウンロードを禁止したいわけじゃないんだろう。

スクレイピングとしては楽な部類で、普通だともうちょっと複雑な認証をかましてくるので自動化しにくい。Captcha入ると流石に無理だけど、おそらく現在の文字認識技術をきちんと応用すれば解けるんじゃないかと思っている。

余談

• こういう解析は楽しい。どんな感じの楽しさかっていうと、与えられたパズルを解いていく感覚にちょっと近い
• 自分はだいたい1時間くらいで解いたけど、早い人はたぶんもっと早いんだろうなって思う。まだまだ精進が足りない

Google Art Projectを利用して高画質の絵画を手に入れる

"The Starry Night" by Vincent van Gogh

2011年2月2日、GoogleさんはGoogle Art Projectというサイトをリリースしました。どういうサイトかっていうとMoMAニューヨーク近代美術館やTate Britainなど、世界中の美術館にある有名な絵画を渡り歩くことができるという趣旨のサイトで、誰もが一度は見たことがある絵画がー家に居ながらー非常にリアリスティックな解像度で見ることができる。恐ろしい時代になったもんです。

ところがこのサイト、高解像度で見られるのは非常にいいことなのですが、肝心のダウンロードができない。せっかく高解像度なんですから保存できてもいいはずなのに、だけどできない。Google Mapなんかは常に最新の情報に更新しなければならないので、ダウンロードに制限をかけているのは合理性があるんですけど、絵画とか全く更新する必要がないでしょう。著作権も切れてるのに。よく分からないです。

ということで、なければ自分でなんとかしようという精神で、1時間くらいでちゃっちゃとPythonでスクリプト組んでダウンロードするようにしてみました。即興で組んだのでいろいろ粗い面はありますけど、そこら辺は勘弁してください。

#!/usr/bin/env python
#-*- coding: utf-8 -*-

import sys
import re
import urllib2
import random
from PIL import Image

def usage():
    print "%s [URL]" % sys.argv[0]
    return 1

def get_id(url):
    response = urllib2.urlopen(url)
    if response.code != 200:
        raise TypeError
    data = response.read()
    return re.search(
        r'<div id="microscope" data-microId="(.+)">', data).group(1)

def get_filename(x, y):
    return "%i_%i.jpg" % (x, y)

def download_image(target_id, x, y, z):
    urllist = [
        "http://lh3.ggpht.com/",
        "http://lh4.ggpht.com/",
        "http://lh5.ggpht.com/",
        "http://lh6.ggpht.com/"
    ]
    base_url = random.choice(urllist)
    url = "%s%s=x%i-y%i-z%i" % (base_url, target_id, x, y, z)

    try:
        response = urllib2.urlopen(url)
    except urllib2.HTTPError:
        return False
    with file(get_filename(x, y), "wb") as fp:
        fp.write(response.read())
    return True

def download_and_get_max_x(target_id, x, z):
    result = download_image(target_id, x, 0, z)
    if result:
        return download_and_get_max_x(target_id, x+1, z)
    else:
        return x

def download_and_get_max_y(target_id, y, z):
    result = download_image(target_id, 0, y, z)
    if result:
        return download_and_get_max_y(target_id, y+1, z)
    else:
        return y

def combine_images(name, max_x, max_y):
    x_size = 512
    y_size = 512

    result = Image.new("RGB", (max_x*x_size, max_y*y_size))
    for y in range(max_y):
        for x in range(max_x):
            image = Image.open(get_filename(x, y))
            result.paste(image, (x*x_size, y*y_size))
    result.save(name)

def main():
    if len(sys.argv) == 1:
        return usage()
    depth = 3
    target_addr = sys.argv[1]
    target_id = get_id(target_addr)
    max_x = download_and_get_max_x(target_id, 0, depth)
    max_y = download_and_get_max_y(target_id, 1, depth)
    for y in range(1, max_y):
        for x in range(1, max_x):
            download_image(target_id, x, y, depth)
    combine_images("result.png", max_x, max_y)

if __name__ == "__main__":
    sys.exit(main())

使い方は簡単で、

$ python artproject.py http://www.googleartproject.com/museums/tate/mariana-248

のようにURLを指定してあげると勝手に取得して結合してくれます。こんな感じで:

"Mariana" by John Everett Millais

あぁそれにしてもお美しい……巷ではオフィーリアちゃんが有名だけど、やっぱりマリアナさんのほうが綺麗だよ。絶対。

余談

• もともとの動機としては高解像度な絵画の壁紙が欲しくて、いろいろと探し回ったわけですが、どうもGoogle Art Projectしか高解像度の絵画が載っていないので仕方なくといった次第です。
• 利用はどうやらGoogle's Terms of Serviceに準じているようです。ざっと見した感じだとダウンロードを禁止する条項はないですし著作権も切れているので全く問題ないように思えるんですが、もし警告とかきたらチキンですから遠慮無く取り下げます。ご了承ください。
• このネタを題材に時間があればHTMLのスクレイピング、パケット解析についての入門記事を書いてみたいです。
• なんだかんだ言ってちゃんと実物見たほうが綺麗。ミレイの作品は以前東京に行ったときにBunkamuraで見ました。Bunkamuraはそういう催し物よくやっていて、そこらへんが東京羨ましかったり

並列環境におけるreduceとscanアルゴリズム

"fireking" by cmyk1219

1. 概要
reduceとscanはGPGPUなどの並列環境下において重要な役割を果たす関数だけど、あまりこれらの関数、特にscanについて言及している記事はあまり見かけない。なので今回はreduceとscanについて調べた結果をまとめていきたいと思う。

2. reduce, scan関数の概要
2.1. reduce
scanについてはreduceが分かればより理解しやすくなるので、ひとまずはreduceから始めることにする。

reduceとは、端的に言うと、対象のリストを与えられた二項演算子で纏め上げる関数だ。とは言っても、いきなりそんなことを言われてもよく分からない。まずは例を見ていこう:

>>> import operator
>>> reduce(operator.add, [1, 2, 3, 4, 5])
15
>>> reduce(operator.sub, [1, 2, 3, 4, 5])
-13
>>> reduce(operator.mul, [1, 2, 3, 4, 5])
120

reduceにおいて重要なのは二項演算子であり、上記の場合、reduceは以下のような演算を行う。

(((1 + 2) + 3) + 4) + 5 = 15
(((1 - 2) - 3) - 4) - 5 = -13
(((1 * 2) * 3) * 4) * 5 = 120

これにより、最初の演算はsum関数、3番目の関数はprod関数と等価であることが分かる。
リストの型は何も数値に限られている訳ではない:

>>> reduce(operator.and_, [True, False, True])
False
>>> reduce(operator.or_, [True, False, True])
True

これより、上記2つの演算はall, any関数と機能的に等価であることが分かるのだけれど、all, any関数のほうが効率的なので、実際はちゃんとall, any関数を使うべきだ。

他にもreduceを用いて様々な関数、アルゴリズムを記述できるが、あまりにも多いので今回は後回しにして、scanを重点的に解説していく。

2.2. scan
scan関数はreduceと比べて文献が少ない。言及している文章もあまり見かけないが、並列処理においてscanは重要な関数だ。残念ながらpythonで実装されていないけれど、仮に実装されていたとしたら以下のような働きをする:

>>> scan(operator.add, [1, 2, 3, 4, 5])
[1, 3, 6, 10, 15]
>>> scan(operator.sub, [1, 2, 3, 4, 5])
[1, -1, -4, -8, -13]
>>> scan(operator.mul, [1, 2, 3, 4, 5])
[1, 2, 6, 24, 120]

きちんと書き下せばどのような働きをしているのか分かりやすい:

[1, 1+2, (1+2)+3, ((1+2)+3)+4, (((1+2)+3)+4)+5]

結局のところ、scanは先頭からreduceしていった結果をリストにまとめ、それを返す関数となる。scanは別名prefix sumとも呼ばれている。
ちなみに、haskellではscanl1(scanl)という名前で実装されている。

Prelude> :type scanl1
scanl1 :: (a -> a -> a) -> [a] -> [a]
Prelude> scanl1 (+) [1,2,3,4,5]
[1,3,6,10,15]

また、CUDAライブラリであるthrustにはinclusive_scan(exclusive_scan)という名前で実装されている。

#include <thrust/scan.h>
int data[6] = {1, 0, 2, 2, 1, 3};
thrust::inclusive_scan(data, data + 6, data); // in-place scan
// data is now {1, 1, 3, 5, 6, 9}

3. 並列による制限
何故並列下においてreduce, scanが重要になってくるのか、それは、両方ともメニーコアな環境において効率的に動作するアルゴリズムが複数提唱されているからだ。reduce, scanは並列と相性がいい。つまり、既存のアルゴリズムをどうにかしてreduce, scanに落とし込むことができれば、並列化の恩恵を受けれられると。

ただし、並列環境においては、演算にいくつかの制限が生じてくる。具体的には、集合Gと演算子"・"の組において、

• 演算子"・"は交換則"a・b = b・a"を満たしていなければならない
• 演算子"・"は結合則"(a・b)・c = a・(b・c)"を満たしていなければならない

後者に関してはあたり前の話で、これが成り立っていなければ並列に計算を行うことができない。前者については多少曖昧で、実際には満たさなくてもよいアルゴリズムも作ることはできるけど、コアレッシングなどの問題を考えると満たしていることが強く求められる。実際、thrustのreduceは交換則を満たしていることを前提としている。

reduceの場合はこの2つの制限で十分だけれど、scanの場合は上記2つの制限に加えてさらに2つの制限を満たしておくと都合がいい。具体的には、集合Gと演算子"・"の組において、

• ただ一つの単位元"e"が存在する
• 演算子"・"は任意の元"a"に対して"a・a^{-1 =} a^-1・a = e"を満たす、ただ一つの逆元"a^-1"が存在する

これら4つの制限を満たしている集合Gと演算子"・"の組を『アーベル群』という。要するに、scan関数は、対象の型から成る集合と二項演算子の組がアーベル群を成していれば都合がいいということになる。

具体例を出そう。整数と和の組(Integer, (+))は明らかに上記4つの制限を満たす。よって、この組は並列にreduce, scanオペレーションを適用でき、かつ都合が良い。

1 + 2 == 2 + 1
(1 + 2) + 3 == 1 + (2 + 3)
1 + (-1) == (-1) + 1 == 0

真偽値と論理積の組(Bool, (&&))は上記2つの制限を満たすため並列にreduce, scanオペレーションを実行可能だが、後半2つの制限は満たしていないため都合は良くない。

3次元の回転行列と積の組(Matrix, (*))は交換則を満たしていないため―少なくともthrustでは―並列にreduce, scanオペレーションを実行することはできない。

3.1. 都合が良いとは？
とは言っても、一体アーベル群であればどこがどう都合がいいのか？それは、リスト中において局所的にreduceを行いたい場合に役に立つ。

リスト[a₀, ... , a_n]にscanオペレーションを実行して

を手に入れたとする。この場合、i < jにおいてb_j+b_i^-1を計算すると、

となる。つまり、i+1からjまでの要素についてのreduceを僅かO(1)で求めることができる。

これには大きな利点がある。例えば、大規模な時系列データを、その周りのデータから平均化(平滑化)したいとする。各要素ごとに計算を行うのは明らかに非効率であるので、まずscanオペレーションを(+)に対して適用し、その後で差分を計算し、正規化を行ってやればよい。結局、scan (+)のやっていることはf(x)の不定積分F(x)を求めることと本質的に等価となる。この考え方はCVでは、積分画像として様々な箇所で用いられる。

そしてさらに、この操作はアーベル群であれば何でも良い……とは言っても、アーベル群であることは結構厳しい制限となる。

例を出すと、リスト中の任意の範囲ーiからjの成分がすべてある条件式に適合しているかどうか確かめたいとする。
関数型の考え方だと、まずはmap関数を用いて(C++でならtransform_iteratorを用いて)リストの成分をBoolに変え、その後でscanを利用することを思いつくかもしれない。しかし、残念ながら(Bool, (&&))はアーベル群でないので、そのまま愚直に適用することができない。

Prelude> scanl1 (&&) $ map (<3) [1,2,3,2,1]
[True,True,False,False,False]
-- 最後の2要素については条件を満たしているのに、Falseで埋もれてしまって分からない！

この場合は、単純にアーベル群である(Integer, (+))をscanに適用してやればよい。

Prelude> scanl1 (+) $ map (\x->if x<3 then 1 else 0) [1,2,3,2,1]
[1,2,2,3,4]
-- きちんと第3要素のみ条件から外れていることが分かる

Pythonっぽく書くとたぶんこんな感じ:

>>> import operator
>>> f = lambda x: 1 if x < 3 else 0
>>> scan(operator.add, [f(x) for x in [1,2,3,2,1]])

4. より詳しく知りたい方は
scanアルゴリズムのGPGPU実装に関する詳細については、現在のところ『Parallel Prefix Sum (Scan) with CUDA(Mark Harris, NVIDIA)』が一番分かりやすい。より詳しく知りたい方はそちらをどうぞ。

というより、基本NVIDIAの出しているテキストは分かりやすい。すごいことです。