読者です 読者をやめる 読者になる 読者になる

プログラミングの実験場

Haskell/Webアプリ/画像処理/可視化/ITによる生産性向上 など

Julialangでの画像処理

Juliaの概要

Juliaという科学技術計算言語が最近出てきた(2009年に開発開始、2012年にオープンソース)のだが、これがかなり良さそうに見える。MatlabやRがやっていることを置き換えられるようなものであると思う。言語仕様やインフラがモダンという印象で、プログラマにとっての使い勝手も良いためか急速に人気も出てきているようで、現時点で中央レポジトリに500程度のパッケージができている。大学の講義などでも使われ始めている。

言語仕様と処理系として特徴的なのは

  • 実行速度が高速。このベンチマークを見ると分かるように、Javaより速く、Cの2倍程度に収まる。Juliaでループを書いてもMatlabのように遅くならない。
    • Juliaの処理系自体(Githubのレポ)も6割以上Juliaで書かれている。
  • ノイズが少なく読みやすい文法。文法に変な制限なども無くRubyのような感じで自然に書ける。Unicodeも使用可能。
  • 多重ディスパッチによる多相。
  • 値や関数の引数のオプショナルな型アノテーション(高速化にも寄与する)。
    • なお、変数には型はない。違う型の値を再代入してもエラーにならない。constはある。
  • レコード型によるオブジェクトの作成。
  • オブジェクトを作れる具象型と、作れない抽象型があり、具象型間の継承は無い。具象型は複数の抽象型を継承しうる。
  • 型パラメータによるジェネリクス。型パラメータには数値も使え、実行時にチェックされる。行列の次元を型で区別したりできる。
  • オンラインでブラウザ上で動くREPL(対話環境)やIDEもある。
  • Githubを使ったパッケージ管理システム。自作のパッケージの登録も簡単にできるようだ。
  • 外部とのバインディングが簡単。
    • 外部ライブラリの利用がバインディングがCで書かれた外部のコードの呼び出しもバインディング不要でccallという関数一つで呼べる。
    • Juliaの関数を外から呼び出すのも簡単そうに見える。
  • マクロ・メタプログラミング機能がある。
    • プログラム中で型情報も含めた構文木(S式)を作成・操作し、実行可能。

多重ディスパッチというのは聞いたことがなかったが、関数の実行時に全引数の型に応じて異なる実装が呼ばれる仕組みで、関数はオブジェクトに所属するのではなく、トップレベルの存在になる。

メタプログラミング機能も、凝ったことができる自由度が無限大という意味ではなかなか興味をそそられる機能だと思う。

関連リンク

現状Pythonで科学技術計算している人も多いと思うが、Pythonとの比較を含めた詳しい解説は、こちらを参照。

PythonistaのためのJulia100問100答 - りんごがでている

速度など、Pythonとの比較。

Julia vs Python: ビットコインオプションのモンテカルロシミュレーション | once upon a time,跡地

高い抽象度とネイティブレベルの壁を超えるという好意的な評価。

「私がJuliaを推す理由」(翻訳) - アドファイブ日記

具体例

私がやりたいことの実例をコードを挙げて例示してみる。

  • 3チャンネルの画像(別々のファイル)をROIで切り取った後、RGB画像として保存する。
import Images: Image,RGB,Gray,imread
import FixedPointNumbers: UfixedBase

function main(channel :: String = "561multi")
    files = split(readall("image_list.txt"),'\n')
    rois = readdlm("roi_list.tsv", '\t')
    idx = findfirst(a -> a == "Slice", rois[1,:])
    count = 0
    for (f,i) in zip(files,[1:length(files)])
        if chomp(f) == ""
            continue
        end
        roilist = choose(rois[2:end,:], idx, i)
        println("$(size(roilist,1)) rois.")

        img1 = imread(replace(f,"CH","488nm"))
        img2 = imread(replace(f,"CH","561nm"))
        img3 = imread(replace(f,"CH","640nm"))

        for i=1:size(roilist,1)
            count += 1
            print(".")
            ns = roilist[i,:]
            crop_with_roi(img1, img2, img3, ns, rois[1,:], count)
        end
    end
end


function choose(dt::Array{Any,2}, idx::Int, slice::Int)
    res = dt[1,:] :: Array{Any,2}
    for i=2:size(dt)[1]
        if dt[i,idx] == slice
            res = [res; dt[i,:]]
        end
    end
    res
end

typealias GrayImage Array{Gray{UfixedBase{Uint16,16}},2}

macro debug(v)
    return :(println($v); $v)
end

function autocontrast(im::GrayImage)
    mx = maximum(im)
    mn = minimum(im)
    (im .- mn) ./ (mx-mn+0.001)
end

function mk_color(img1::GrayImage,img2::GrayImage,img3::GrayImage)
    a,b = size(img1)
    cr2 = Array(UfixedBase{Uint16,16},a,b,3)
    cr2[:,:,1] = autocontrast(img1)
    cr2[:,:,2] = autocontrast(img2)
    cr2[:,:,3] = autocontrast(img3)
    convert(Image{RGB},cr2)
end

function crop_with_roi(img1,img2,img3,row,header,num)
    out_folder = "/Volumes/MacintoshHD/output_test/"
    y = convert(Int,row[findfirst(a -> a == "BY", header)][1])
    x = convert(Int,row[findfirst(a -> a == "BX", header)][1])
    h = convert(Int,row[findfirst(a -> a == "Height", header)][1])
    w = convert(Int,row[findfirst(a -> a == "Width", header)][1])
    cr1 = img1[x+1:x+w,y+1:y+h]
    cr2 = img2[x+1:x+w,y+1:y+h]
    cr3 = img3[x+1:x+w,y+1:y+h]
    Images.imwrite(mk_color(cr1,cr2,cr3), out_folder * "test$num.tif")
end

CSVの読み方がよくわからなかったので適当なコードを書いた。Imagesというパッケージを追加で使っている。JuliaのREPLから、Pkg.add("Images")で入れられる。

多次元配列の使い方はMatlabのような感じ。ただし、いくつか違いもある。

Noteworthy Differences from other Languages — Julia Language 0.4.0-dev documentation

250x250ピクセルの3チャネルTIFFを300個切り出して保存するのに20秒ほど。二回目以降は事前コンパイルされているせいか、7秒程度に短縮された。悪くない。

若干不満な点

Juliaにはかなり期待していてこれから書く画像処理はしばらくJuliaを使ってみようと思うのだが、若干不満な部分もあった。

関数の型がない

関数はすべてFunctionという一つの型であり(例えば、typeof(map) == Function)、Haskellみたいにmap :: (a -> b) -> [a] -> [b]みたいな型アノテーションは書けない。 異なる型の関数を同じ名前に割り当てる多重ディスパッチとは相性が悪いようにも思うので、しょうがない気もするが。

実行時エラーがそれなりに多い

型が強すぎる言語だと、特に書き捨ての小さなスクリプトでは面倒になるので、要はバランスなのだが、演算中に型が合わずにエラーということが感覚としては結構多かった。慣れの問題はあると思うが。 変数の未定義なども実行時エラーとなるので、そのあたり実行前にもうちょっと強くチェックしてくれるとありがたいと思う。

importがかなり遅い

Pkg.add("Images")でパッケージを入れた上に、importするのだが、REPLを起動した後の最初のimportが異様に遅い。数秒待たされる。何かやり方が間違っているのかも...。