Julia で Grass-Fire algorithm を書く
なんの因果か最近は画像をぐりぐりしています。
画像の明るい点を検知してその重心位置を求める必要ができたのですが、考えてみるとなかなか難しい。いろいろ調べてみるとどうやらこういうのは BLOB (Binary Large OBject) detection というらしいです。厳密には BLOB extraction かな。
単純でいいのでとにかく速いものが欲しい、ということで探しているとここの Grass-Fire algorithm 解説を見つけました。どうも一次情報は [T. B. Moeslund, Introduction to Video and Image Processing: Building Real Systems and Applications, Springer London, London, 2012.] という本みたいですが、引用とかないし大丈夫かこれ。幸い、本を図書館で利用できたので参考にしました。
ちなみに真っ先に試したのはImages.jlの blob_LoG 関数なんですが、ちょっとほしいものと違う感じでした。BLOBの中心は得られるのですが、できればサブピクセルレベルで、つまり小数点以下のレベルで座標を欲しかったんです。よく理解していないだけかもしれないです。OpenCV とかに便利で速い関数とかありそうなんですが、今のところ julia から OpenCV を使うのは簡単ではなさそう…。Grass-Fire algorithm は経路探索など様々な用途で使われるもののようです。とりあえず私の目的も果たせそうなのでいいでしょう。
いかにも再帰が向いていそうな感じなので、とりあえずその方向で書きました。閾値の thr と同じかより大きい値を持つピクセルの連続を BLOB として認識し、座標の配列として返します。
function extractblob_rec(img, thr) height, width = size(img) checked = fill!(similar(img, Bool), false) bloblist = Vector{NTuple{2,Int}}[] for x in 1:width, y in 1:height # Scan pixels until a new BLOB @inbounds if img[y, x] >= thr && !checked[y, x] # BLOB detection by Grass-Fire algorithm blob = NTuple{2,Int}[] grassfire!(blob, checked, img, y, x, thr) push!(bloblist, blob) end end bloblist end function grassfire!(blob, checked, img, y, x, thr) height, width = size(img) push!(blob, (y, x)) checked[y, x] = true # check the lower pixel if y+1 <= height && img[y+1, x] >= thr && !checked[y+1, x] grassfire!(blob, checked, img, y+1, x, thr) end # check the right pixel if x+1 <= width && img[y, x+1] >= thr && !checked[y, x+1] grassfire!(blob, checked, img, y, x+1, thr) end # check the upper pixel if y-1 >= 1 && img[y-1, x] >= thr && !checked[y-1, x] grassfire!(blob, checked, img, y-1, x, thr) end # check the left pixel if x-1 >= 1 && img[y, x-1] >= thr && !checked[y, x-1] grassfire!(blob, checked, img, y, x-1, thr) end blob end
参考と同じ簡単な 6x5 の配列を作ってテストしてみます。
julia> img
6×5 Array{Int64,2}:
0 0 1 1 1
0 0 0 1 0
0 0 0 1 0
0 1 1 0 0
0 1 1 0 0
0 1 1 0 0
julia> extractblob_rec(img, 1)
2-element Array{Array{Tuple{Int64,Int64},1},1}:
[(4, 2), (5, 2), (6, 2), (6, 3), (5, 3), (4, 3)]
[(1, 3), (1, 4), (2, 4), (3, 4), (1, 5)]
二つの BLOB を見つけたみたいです、よさそうですね。探す様子を GIF 化しました。赤が走査している場所、緑と青が BLOB です。
さて、再帰を使っていると StackOverflowError が怖いです。たしか Julia ではまだ末尾呼び出し最適化されないので、BLOB が大きくなるとスタックがあふれます。私の実際の用途を考えると BLOB のサイズはそんなに大きくならないのですが、いちおう再帰を使わない形も考えてみます。
import DataStructures: Stack function extractblob_loop(img, thr) height, width = size(img) stack = Stack{NTuple{3,Int}}() checked = fill!(similar(img, Bool), false) bloblist = Vector{NTuple{2,Int}}[] for x in 1:width, y in 1:height # Scan pixels until a new BLOB @inbounds if img[y, x] < thr || checked[y, x] # BLOB detection by Grass-Fire algorithm # NOTE: re-use `stack` because making a new one is a little costly blob = grassfire_loop!(checked, stack, img, y, x, thr) push!(bloblist, blob) end end bloblist end function grassfire_loop!(checked, stack, img, y0, x0, thr) height, width = size(img) blob = NTuple{2,Int}[] push!(stack, (y0, x0, 4)) push!(blob, (y0, x0)) checked[y0, x0] = true while !isempty(stack) y, x, state = pop!(stack) checked[y, x] = true # check the lower pixel if state >= 4 && y+1 <= height && img[y+1, x] >= thr && !checked[y+1, x] push!(stack, (y, x, 3)) push!(stack, (y+1, x, 4)) push!(blob, (y+1, x)) continue end # check the right pixel if state >= 3 && x+1 <= width && img[y, x+1] >= thr && !checked[y, x+1] push!(stack, (y, x, 2)) push!(stack, (y, x+1, 4)) push!(blob, (y, x+1)) continue end # check the upper pixel if state >= 2 && y-1 >= 1 && img[y-1, x] >= thr && !checked[y-1, x] push!(stack, (y, x, 1)) push!(stack, (y-1, x, 4)) push!(blob, (y-1, x)) continue end # check the left pixel if state >= 1 && x-1 >= 1 && img[y, x-1] >= thr && !checked[y, x-1] push!(stack, (y, x, 0)) push!(stack, (y, x-1, 4)) push!(blob, (y, x-1)) continue end end blob end
ランダムな配列でテストしたところ extractblob_rec と同じ結果を返すようです。
次に 800x600 ピクセルに50個程度の輝点を持つ画像を用意し、性能を確認しました。ベンチマークの結果は以下です。
julia> versioninfo() Julia Version 1.1.0 Commit 80516ca202 (2019-01-21 21:24 UTC) Platform Info: OS: Windows (x86_64-w64-mingw32) CPU: Intel(R) Core(TM) i7-5500U CPU @ 2.40GHz WORD_SIZE: 64 LIBM: libopenlibm LLVM: libLLVM-6.0.1 (ORCJIT, broadwell) julia> using BenchmarkTools julia> @benchmark extractblob_rec(testimg, 240) BenchmarkTools.Trial: memory estimate: 482.02 KiB allocs estimate: 135 -------------- minimum time: 408.409 μs (0.00% GC) median time: 453.314 μs (0.00% GC) mean time: 515.280 μs (4.78% GC) maximum time: 60.672 ms (98.67% GC) -------------- samples: 9601 evals/sample: 1 julia> @benchmark extractblob_loop(testimg, 240) BenchmarkTools.Trial: memory estimate: 506.22 KiB allocs estimate: 140 -------------- minimum time: 409.692 μs (0.00% GC) median time: 454.168 μs (0.00% GC) mean time: 512.038 μs (5.07% GC) maximum time: 58.005 ms (99.16% GC) -------------- samples: 9665 evals/sample: 1
800x600 の画像で 1ms を切るのが目標だったので、とりあえずどちらも目標達成ですね。非力なノート PC で達成できたので十分といえるでしょう。使用メモリーのほとんどが checked 配列に費やされているので、これを再利用すると 1/10 ぐらいになります。median time を見ると再帰を使ったほうがちょっと速そうに見えますが誤差です。何回か測定すると前後します。
最終的には再帰のほうがコードは簡潔ですけど、ちょっと悩みます。どっちを使うべきか…。あと、まだどこか高速化する余地はあるか考えてみます。重心位置の計算は省略。
