高階関数とは、引数に関数を取る関数のことです。MRC分子生物学研究所のImageJチュートリアルで、forループを使う代わりに高階関数を使ってピクセルの反復処理を行う事例を紹介しています。高階関数の説明自体は例えば下記のサイトをご覧ください。
高階関数というとmap, filter, reduceがよく言及されるようです。mapは一つの引数を持つ関数と、リスト(配列)あるいは複数の数値を引数にとり、引数として取った関数に次々と数値を連続して与える関数、filterは一定の条件を弁別する関数と、やはりリストや複数の数値を引数にとり、一定の条件をクリアした値のみ戻り値として返す関数、reduceはいわゆる畳み込み処理を行う関数で、二つの引数を取って一つの結果を返す関数と、やはりリストや複数の数値を引数として取り、複数の数値を引数として取った関数を使って一つの数値に何らかの形で縮約(reduce)する関数です。代表的な例として、総和を求めたり、最小値、最大値を求める手続きがあり得ます。詳しくは下記のサイトなどをご覧ください。
kakedashi-engineer.appspot.com
チュートリアルによると、forループで回すより高階関数を使ったほうが処理が速くなることがあるそうです。
では、MRC分子生物学研究所のImageJチュートリアルのサンプルを見ていきましょう。まずピクセルの反復処理(Iterating Pixels)です。内容は、1つの画像(ImageProcessor)の最小値を求めるというプログラムです。
from ij import IJ
from sys.float_info import max as MAX_FLOAT
# Grab the active image
imp = IJ.getImage()
# Grab the image processor converted to float values
# to avoid problems with bytes
ip = imp.getProcessor().convertToFloat() # as a copy
# The pixels points to an array of floats
pixels = ip.getPixels()
print "Image is", imp.title, "of type", imp.type
# Obtain the minimum pixel value
# Method 1: the for loop, C style
minimum = MAX_FLOAT
for i in xrange(len(pixels)):
if pixels[i] < minimum:
minimum = pixels[i]
print "1. Minimum is:", minimum
# Method 2: iterate pixels as a list
minimum = MAX_FLOAT
for pix in pixels:
if pix < minimum:
minimum = pix
print "2. Minimum is:", minimum
# Method 3: apply the built-in min function
# to the first pair of pixels,
# and then to the result of that and the next pixel, etc.
minimum = reduce(min, pixels)
print "3. Minimum is:", minimum
import文にある sys.float_info というのは、Pythonで用意されている名前付きタプルのようです*1。
画像のピクセルの各値をPixelsという一次元のリスト(配列)に読み込み、処理をしています。ここで高階関数を使っているのは、Method3です。minimum = reduce(min, pixels) という文ですが、minという組み込み関数は、二つの引数で取られた値を比較して小さいほうを戻り値として返す関数です。それをreduce関数にPixelsと共に読み込むことでPixels[0]からPixcels[n]までの全ピクセルの値をminに放り込んで比較し、その中の最小値を最終的に求めています。
次は、リストまたは要素の集合の反復処理について(On iterating or looping lists or collections of elements)です。やっていることは、ImageJで複数の画像(ウィンドウ)を開いているときに、その画像をimagesというリスト変数に読み込みます。
from ij import WindowManager as WM
# Method 1: with a 'for' loop
images =
for id in WM.getIDList():
images.append(WM.getImage(id))
# Method 2: with list comprehension
images = [WM.getImage(id) for id in WM.getIDList()]
# Method 3: with a 'map' operation
images = map(WM.getImage, WM.getIDList())
Method3で高階関数mapを使っています。WM.getImage() メソッドと開いているウィンドウのIDリストを読み込み、WM.getImageに開いているウィンドウの画像を片っ端から読み込ませることで、imagesにウィンドウの画像を収容します。
次は各ピクセルからピクセルの最小値を減算する(Subtract the min value to every pixel)です。
from ij import IJ, ImagePlus
from ij.process import FloatProcessor
imp = IJ.getImage()
ip = imp.getProcessor().convertToFloat() # as a copy
pixels = ip.getPixels()
# Apply the built-in min function
# to the first pair of pixels,
# and then to the result of that and the next pixel, etc.
minimum = reduce(min, pixels)
# Method 1: subtract the minim from every pixel,
# in place, modifying the pixels array
for i in xrange(len(pixels)):
pixels[i] -= minimum
# ... and create a new image:
imp2 = ImagePlus(imp.title, ip)
# Method 2: subtract the minimum from every pixel
# and store the result in a new array
pixels3 = map(lambda x: x - minimum, pixels)
# ... and create a new image:
ip3 = FloatProcessor(ip.width, ip.height, pixels3, None)
imp3 = ImagePlus(imp.title, ip3)
# Show the images in an ImageWindow:
imp2.show()
imp3.show()
まず、import文で、ImageProcessorの一種であるfloatProcessor (浮動小数点を使った画像データ)を呼び出します。
そしてピクセルの最小値を求めるのにreduce関数を使っています(前の例で出ました)。
そしてMethod2が、高階関数を使った例です。
pixels3 = map(lambda x: x - minimum, pixels)
map関数を使っていますが、引数として参照している関数は匿名のラムダ関数です。このラムダ関数の引数は x そして関数の処理内容は x - minimumです。そして引数 x にmap関数でpixcels を割り付けています。これにより全ピクセルの値から、全ピクセルの最小値を引き算します。
その後、FloatProcessorをコンストラクタとして使って、pixels3からimageProcessor ip3 を作成し、それからImagePlusオブジェクト imp3を作成します。
次のケースは finding the coordinates of all pixels above a certain value (一定の値を超えるすべてのピクセルの座標を見つける)です。具体的には、平均値を超えるピクセルの中心座標を取得するプログラムです。
from ij import IJ
# Grab the currently active image
imp = IJ.getImage()
ip = imp.getProcessor().convertToFloat()
pixels = ip.getPixels()
# Compute the mean value
mean = sum(pixels) / len(pixels)
# Obtain the list of indices of pixels whose value is above the mean
above = filter(lambda i: pixels[i] > mean, xrange(len(pixels)))
print "Number of pixels above mean value:", len(above)
# Measure the center of mass of all pixels above the mean
# The width of the image, necessary for computing the x,y coordinate of each pixel
width = imp.width
# Method 1: with a for loop
xc = 0
yc = 0
for i in above:
xc += i % width # the X coordinate of pixel at index i
yc += i / width # the Y coordinate of pixel at index i
xc = xc / len(above)
yc = yc / len(above)
print xc, yc
# Method 2: with sum and map
xc = sum(map(lambda i: i % width, above)) / len(above)
yc = sum(map(lambda i: i / width, above)) / len(above)
print xc, yc
# (Continues from above...)
# Method 3: iterating the list "above" just once
xc, yc = [d / len(above) for d in
reduce(lambda c, i: [c[0] + i % width, c[1] + i / width], above, [0, 0])]
print xc, yc
# Method 4: iterating the list "above" just once, more clearly and performant
from functools import partial
def accum(width, c, i):
""" Accumulate the sum of the X,Y coordinates of index i in the list c."""
c[0] += i % width
c[1] += i / width
return c
xy, yc = [d / len(above) for d in reduce(partial(accum, width), above, [0, 0])]
print xc, yc
前のサンプルと同じように一旦すべてのピクセルを pixels という1次元配列に収容しています。平均値を計算する際は、sum関数で総和を取り、ピクセルの個数はlen関数でpixcelsの長さを求めることで個数を求めています。
このとき、pixelsのインデックス番号はどの様につくかというと、仮に、画像の幅 (ピクセル数) が、w 高さが h とすると...
1行目の1番目からw番目のピクセルのインデックス番号は...
0 〜 (w - 1) までのインデックス番号が付きます。
2行目の1番目からw番目のピクセルのインデックス番号は...
w 〜 (2w - 1) までのインデックス番号が付きます。
3行目の1番目からw番目のピクセルのインデックス番号は...
2w 〜 (3w - 1) までのインデックス番号が付きます。
4行目の1番目からw番目のピクセルのインデックス番号は...
3w 〜 (4w - 1) までのインデックス番号が付きます。
:
:
最後 h行目の1番目からw番目のピクセルのインデックス番号は...
(h-1)w 〜 (hw - 1) までのインデックス番号が付きます。
この最後のインデックス番号 hw -1 は len(pixels) -1 に一致します。このインデックス番号からもとの x, y 座標を取得するには、x 座標は、インデックス番号を幅の w で割った余りを、y 座標は インデックス番号を w で割った商を計算することで得られます。因みに (0, 0) は画像の一番左上になります。右に進むほど x が増え、下に進むほど y が増えます。
私だったら、この辺りは何とかして2次元の配列で扱いたいところですが、このチュートリアルでは配列が2次元で扱われるケースは皆無です。これはPythonスタイルなのでしょうか? それとも多次元配列を扱わないほうが効率的なのでしょうか? 個人的にはむしろ分かりにくくて気持ち悪いのですが... 因みに自作の簡易版Bチャンネル再建法ツールでは2次元配列を使っています。
次に高階関数のfilterを使って、平均値以上のピクセルを収容した1次元配列 above を求めます。filterの関数引数は平均以上かどうかを判定するラムダ関数 lambda i: pixels[i] > mean を取り、数値の引数に各ピクセルを取っています。つまり、平均以上のピクセルをaboveに入れるわけです。数値の引数に xrange(len(pixels)) を取る理由はlen(pixcels) とするよりは、メモリ、処理速度の効率を求めているためと思われます。
Method1は、解説によると、C言語風にまず変数を初期化し for ループを使った例です。ここで紹介されている4つの手法の中で最も効率的だそうです。このコード中、forループのなかの
xc += i % width # the X coordinate of pixel at index i
yc += i / width # the Y coordinate of pixel at index i
という部分が気になります。この部分で、xc, ycの平均値を求める前に総和を求める過程で above[i] のインデックス番号 i を画像の幅で割った余りを求めています。これはaboveは一次元配列になっているため、添え字から直接 x, y 座標を取得することができません。そこで i を幅で割った余りを x 座標としています。一方 y 座標は i を幅で割ったものを y 座標としているようです。これは上の、インデックス番号の決まり方の説明をご参照いただくとわかります。i が整数型とは明示されていませんが、そう扱われているようです。というのは i / witdth が Y 座標なら整数型でなければならないからです。そして、この x 座標と y 座標のそれぞれの総和を求めます。
こういうことが可能ということは、filter関数を使って所定条件のデータのみ残すときに、残されたデータが収容されるabove[]配列には元々のインデックス番号は、欠落したデータの部分は後から来たデータで穴埋めされず、データ内容が欠落したままインデックス番号自体は維持されるということらしいです。あるいは欠落部分が0で穴埋めされることもないということかと。なるほど。
さらにMethod2では、平均以上のピクセルの中心座標、xc, ycを求めていますが、そこでmap関数を使っています。具体的には、
map(lambda i: i % width, above)
とあります。Method1と同様、この式ではabove[0]からabove[n]までのインデックス番号を画像の幅で割った余りを求めています。こうすることで、Method1よりコードが短くなっています。解説ではMethod1より可読性が高い一方、効率面ではMethod1に劣るとされています。というのは、
xc = sum(map(lambda i: i % width, above)) / len(above)
と1つの式の中で above 関数が2つ出てきています。この分速度が遅くなるようです。しかし、確かにコードは短いですが、個人的にはMethod2がMethod1より可読性が高いか? と思ってしまいます。
Method3では、やっていることはMethod2と同じですが、さらに関数を活用して、それを1行で済ませています。
xc, yc = [d / len(above) for d in reduce(lambda c, i: [c[0] + i % width, c[1] + i / width], above, [0, 0])]
まず、ラムダ関数から見ていきます。
lambda c, i: [c[0] + i % width, c[1] + i / width]
このラムダ関数は、引数が、c, i 、戻り値が配列 [c[0] + i % width と [c[0] + i % width の２つとなっています。また引数は配列cとiの二つ(配列cを2つと考えれば３つ)あります。c[0]はx座標の総合計を求める変数、c[1]はy座標の総合計を求める変数です。このラムダ関数が、reduce関数の引数となっています。
reduce(lambda c, i: [c[0] + i % width, c[1] + i / width], above, [0, 0])
reduce関数は、引数が2つの場合と3つの場合がありえますが、2つの場合は reduce (関数, 反復要素)、3つの場合は reduce (関数, 反復要素, 初期値) となり、この場合は、関数に先程のラムダ関数、反復要素として above 配列、そして c の初期値として [0, 0] が設定されています。戻り値は、関数の計算結果です。関数の中で i ではなく、cが2要素の配列として定義されていますので、[0, 0] は i の初期値ではなく c の初期値として解釈されるものと思います。そしてラムダ関数の i に above 配列のインデックスが入っていきます。
従って、最初に c の要素として [0, 0]が設定され、その後 above の1番目の要素から n 番目の要素までのインデックス番号が i に代入され reduce 関数によって加算されていきます。
そして、
xc, yc = [d / len(above) for d in reduce... ]
結局reduceは2つ戻り値がありますので、d も二つある(つまりc[0]とc[1]の最終結果をdに代入)ことになります。これを平均以上のピクセル数で割って、平均座標値である xc, yc を求める、という計算になっているようです。
解説によると、コードは最も短いものの紛らわしい部分もあり(おそらく [0,0]の初期値が、c に与えられるのか i に与えられるのかという部分かと思われます)、またこの4つの例の中では最も不効率だそうです。
最後にMethod4ですが、Method3をクリーンアップしたものだそうです。但し、長さ的には for ループを使ったものと変わりません。関数の引数を部分的に固定する partialという高階関数を使っています。この関数を使うために、
from functools import partial
でライブラリを呼び出しています。
まずあらかじめc[0], c[1]を総和するaccumという関数を定義しています。
xy, yc = [d / len(above) for d in reduce(partial(accum, width), above, [0, 0])]
partialを使っているのはこのコードです。partialという関数は、あらかじめ定義されている関数の引数を部分的に固定して活用する高階関数です。accumという関数は3つの引数 (width, c, i) を取ります。このうち、partial(accum, width) というコードで、引数 width を 変数widthを取るよう固定します。従って可変となる引数は c と i の2つになります。ただし、c と i に順々に入れられるのはあくまで aboveという配列(リスト変数)であり、その後の[0,0]はc[0]とc[1]の初期値です。
最後に、2つの要素が d / len(above) で除算されて平均値が計算されるのは同じです。
なお、高階関数を使うと、一時変数を使わなくて済むので潜在的エラーが減ると、チュートリアルでは言っていますが、どうなんでしょう。個人的にはプログラムが短くなるよりは for roopを使ったほうが可読性が高く分かりやすいようにも思います。もちろん大きな画像を扱う場合に、プログラムが高速化に寄与するのであれば、採用に躊躇はありませんが、短い行数でプログラムを書くことは必ずしも可読性を高めるとは限りませんので...
*1:以下参照。
