EBImageを使った画像処理

使い方とか.
参考

インストール

source("http://bioconductor.org/biocLite.R")
biocLite("EBImage")

このほかにImageMagicとgtkが必要.

以下http://www.bioconductor.org/packages/devel/bioc/html/EBImage.htmlにあるIntroduction to EBImageを読んでメモしたもの.長いです.気が向いたら目次付けます.付けました.

目次

画像の読み込みと保存

画像の読み込みにはreadImage()関数を使用する.読み込んだ画像を表示するのはdisplay().

## 画像の読み込みはreadImage()
oryzae <- readImage("oryzaeg.png")
## 画像の表示はdisplay()
display(oryzae)


カラー画像や複数のフレームを持つ画像についても同様に読み込むことができる.

## カラー画像の読み込み
oryzaec <- readImage("oryzae.png")
display(oryzaec)                   

## 複数のフレームを持つ画像の読み込み
nuc <- readImage(system.file("images","nuclei.tif",package="EBImage"))
display(nuc)


画像の保存はwriteImage()関数を使用する.出力形式は保存ファイル名の拡張子から自動的に推定される.以下JPEG形式で保存する例(元画像はPNG).

## 画像の保存はwriteImage 形式は拡張子から推定される
writeImage(oryzae, "oryzae.jpeg", quality = 95)
writeImage(oryzaec, "oryzaec.jpeg", quality = 95)

イメージオブジェクトとマトリックス

readImage()で読み込まれた画像はImageクラスのオブジェクトとして保持される.形式は各ピクセルの輝度情報を含む多次元配列となる.EBImageに含まれる全ての関数はmatrixオブジェクト,arrayオブジェクトに対しても同様に使用可能.

> ## グレースケールの場合
> print(oryzae)
Image
  colormode: Grayscale 
  storage.mode: double 
  dim: 300 225 
  nb.total.frames: 1 
  nb.render.frames: 1 

imageData(object)[1:5,1:6]:
          [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
[1,] 0.4967269 0.4967269 0.4967269 0.4980392 0.4980392 0.4993515
[2,] 0.4967269 0.4941176 0.4901961 0.4967269 0.4980392 0.5006485
[3,] 0.4993515 0.4888838 0.4941176 0.4980392 0.4980392 0.5019608
[4,] 0.4993515 0.4954299 0.4967269 0.4967269 0.5019608 0.5045701
[5,] 0.4941176 0.4941176 0.4941176 0.4967269 0.4980392 0.4993515
> ## カラーの場合(次元が一つ多い)
> print(oryzaec)
Image
  colormode: Color 
  storage.mode: double 
  dim: 300 225 3 
  nb.total.frames: 3 
  nb.render.frames: 1 

imageData(object)[1:5,1:6,1]:
          [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
[1,] 0.4980392 0.5058824 0.4980392 0.5019608 0.5019608 0.5019608
[2,] 0.4980392 0.5019608 0.4941176 0.5019608 0.5019608 0.5019608
[3,] 0.5058824 0.4941176 0.4980392 0.5058824 0.5058824 0.5058824
[4,] 0.5058824 0.4980392 0.4980392 0.5019608 0.5098039 0.5137255
[5,] 0.4980392 0.5019608 0.5019608 0.5019608 0.5058824 0.5098039

Imageオブジェクトはmatrixと同様にRの算術演算を施すことができる.これにより簡単な画像の加工ができる.

oryzae1 <- oryzae + 0.5      ## 明度操作

oryzae2 <- oryzae * 3        ## コントラスト

oryzae3 <- (oryzae + 0.2)^3  ## ガンマ補正

oryzae4 <- oryzae[299:376, 224:301]  ## 切り取り

oryzae5 <- oryzae > 0.5      ## 二値化

oryzae6 <- t(oryzae)         ## 転置(反転→左90度回転)


統計情報の計算.

> median(oryzae)
[1] 0.3803922

複数のフレームを持つ画像はcombine()関数で作成できる.

## 複数のフレームを持つ画像
oryzaecomb <- combine(oryzae, oryzae*2, oryzae*3, oryzae*4)
display(oryzaecomb)

複数のフレームを持つ画像をwriteImage()で保存する場合,tifやgifのように複数フレームに対応している形式ならそのまま,jpegのように対応していない形式ならフレームごとに分割された上保存される.

writeImage(oryzaecomb, "oryzaecomb.tif") ## そのまま
writeImage(oryzaecomb, "oryzaecomb.jpeg") ## oryzaeconb-0.jpegからoryzaecomb-3.jpegまでの4枚が保存される





画像の変形

rotate(),translate(),flip()といった関数が用意されている.

oryzae7 <- rotate(oryzae, 30)           ## 回転
oryzae8 <- translate(oryzae, c(40, 70)) ## 平行移動
oryzae9 <- flip(oryzae)                 ## 上下反転




色の取り扱い

Imageクラスオブジェクトはcolormodeというスロットにより多次元配列の色情報をどのように取り扱うべきかを判断している.colormodeスロットにはcolorMode()関数によりアクセスすることができ,変更することもできる.

colorMode(oryzaec) <- Grayscale  ## Grayscaleオブジェクトの中身は数値型の0なので,0を代入してもよい
display(oryzaec)                 ## グレースケールで表示される
colorMode(oryzaec) <- Color      ## カラーに戻す.Colorの中身は2.

colormodeスロットを変更しても画像の表示方法が変わるだけで画像データに変更はない.もしもカラー画像をグレースケールにした場合,それは複数のフレームを持つグレースケールの画像として取り扱われる.
channel()関数も画像の表示方法を変更する.colorMode()と異なるのは,colorMode()が輝度情報を変更しないのに対してchannel()関数は変更する点にある.どのような形式に変更するかは第二引数modeで指定する.modeには以下のようなものがある.

  • rgb: Grayscaleの画像を複製しR, G, Bのチャンネルに割り当てる
  • gray, grey: カラー画像のR, G, Bそれぞれの輝度に1/3のウェイトをかけて1枚のグレースケール画像を作成する
  • red, green, blue: カラー画像のR,G,またはBのチャンネルに含まれる輝度情報を抽出してグレースケール画像を作成する
  • asred, asgreen, asblue: グレースケール画像をカラー画像に変更し,輝度情報をR,GまたはBのチャンネルに割り当てる.

グレースケールをカラー画像に変更する例を示す.

oryzaek <- channel(oryzae, "rgb")
oryzaek[236:276, 106:146, 1] <- 1
oryzaek[236:276, 156:196, 2] <- 1
oryzaek[236:276, 206:246, 3] <- 1


rgbImage()関数は3つのグレースケール画像をR,G,Bのチャンネルに割り当て,一枚のカラー画像として合成する.

oryzaet <- oryzae[1:225, 1:225]
oryzaeb <- rgbImage(red=oryzaet, green=flip(oryzaet), blue=flop(oryzaet))


フィルタリング

filter2()関数を利用すると,画像に線形フィルターをかけることができる.フィルターとは何か,ということについては,例えばhttp://www.gifu-nct.ac.jp/elec/yamada/iwata/filter/index.htmlなどに詳しい説明がある.
makeBrush()関数を利用すると様々なフィルターを簡単に作成することができる.
makeBrush()関数で作成したフィルターを適用する例を以下に示す.

flo <- makeBrush(21, shape='disc', step=FALSE)
flo <- flo/sum(flo)
oryzaeflo <- filter2(oryzaec, flo)


フィルターはmatrixなので,値を指定して自作することもできる.以下に示すのはエッジ検出などに使用されるラプラシアンフィルターを作成し適用する例.

fhi <- matrix(1, ncol=3, nrow=3)
fhi[2,2] <- -8
oryzaefhi <- filter2(oryzaec, fhi)

画像へ文字を書き込む

drawtext()関数を利用すると画像へ文字を書き込むことができる.このとき,引数fontを指定することでフォントの設定ができるのだが,引数fontにはDrawFontクラスのオブジェクトを渡す.DrawFontクラスのオブジェクトはdrawfont()関数で作成する.
つまり,

  1. drawfont()でDrawFontクラスのオブジェクトを作成
  2. drawtext()のfont引数に作成したオブジェクトを渡す

という流れで文字を書き込む(日本語は使えない?).

font <- drawfont(weight=600, size=28)
oryzaetext <- drawtext(oryzae, xy=c(125, 200), font=font, labels="Oryzae", col="black")


モルフォロジカル処理

線形フィルターとよく似ているが,非線形のフィルターを使用するモルフォロジカル処理と呼ばれる処理がある(詳細はhttp://www.roper.co.jp/mcj/webhelp/morphological_filters.htmなどを参照).明るい部分を収縮させるErodeと,逆に膨張させるDilateが基本の処理である.二値画像に適用すれば細い線を切断したり,不連続な点を接続したりすることができる.モルフォロジカル処理に使用するフィルタも行列であるため,makeBrush()で作成することができる.

## サンプル画像の読み込み
ei <- readImage(system.file("images", "shapes.png", package="EBImage"))
ei <- ei[110:512, 1:130]
## フィルタの作成
kern <- makeBrush(5, shape="diamond")

これは元の画像.

erodeの例.

## erodeの例
eierode <- erode(ei, kern)
writeImage(eierode, "eierode.png")


次にdilateの例.

## dilateの例
eidilat <- dilate(ei, kern)
writeImage(eidilat, "eidilat.png")


オブジェクトの認識

bwlabel()関数は2値画像にのみ使用可能な関数で,接続している一つの領域に対してラベルとしてユニークな整数を1つ割り当てる.二値化後にbwlabel()関数を適用することで画像内のオブジェクト数を数えることができる.ラベルは1から順に増加する整数であるので,ラベルの最大値が画像内のオブジェクト数と一致する.

> eilabel <- bwlabel(ei)
> cat("オブジェクトは", max(eilabel), "つ含まれています.\n")
オブジェクトは 7 つ含まれています.

またここでeilabel/max(eilabel)として値を0〜1の間に収まるように加工し,display()などで表示すれば,どのように領域が分割されているのかを視覚的に知ることができる.

display(eilabel/max(eilabel))


以下は画像を2値化してbwlabel()を適用する例である(画像nucは上の方で読み込み済み).

nuct <- nuc[,,1]>0.2
nuclabel <- bwlabel(nuct)
cat("核は", max(nuclabel), "個あります.\n")
writeImage(nuc[,,1], "nuc.png")
writeImage(nuclabel/max(nuclabel), "nuclabel.png")

元の画像

二値化後の画像

しかし,画像をよく見ると分かるように一部の核がくっついてしまっている.元の画像を見ると重なっているわけではないので,上手く処理すれば区別できそうに見える.
まず,二値化処理の方法を変える.この画像の場合はそれほど問題ではないのだが,光の加減などにより画像の場所により明るさが異なってしまっているような場合,画像の全ての場所で同じ閾値を使用するのは賢明ではない.画像全体ではなく,閾値処理をしたい画素周辺の限られた領域の輝度情報を参考に二値化処理をすべきである.このように画面の部分等に応じて閾値を変化させる方法を可変閾値法と呼ぶ.
thresh()関数は指定した大きさの窓を移動させ,窓内の情報を元に画像を二値化する関数である.また,erode後にdilateする操作はopenと呼び,細い部分を切断する効果がある(同様の操作をする関数としてopening()がある).thresh()により二値化し,openにより細い接続を切断した画像にbwlabel()を適用した結果を示す.

nuct2 <- thresh(nuc[,,1], 10, 10, 0.05)
kern <- makeBrush(5, shape="disc")
nuct2 <- dilate(erode(nuct2, kern), kern)
nuclabel2 <- bwlabel(nuct2)
cat("核は", max(nuclabel2), "個あります.\n")
writeImage(nuclabel2/max(nuclabel2), "nuclabel2.png")

さきほどくっついてしまっていた部分が分かれているのが分かる.

オブジェクトの取り扱い

画像内の一つのオブジェクトは同じ整数値のラベルを持つピクセルの集合として定義される.paintObjects()関数を用いると,オブジェクトの周縁部(とオプションで内部)を強調表示することができる.

## 下敷きにする画像
nucgray <- channel(nuc[,,1], 'rgb')
## オブジェクトの情報(ラベル)を含む画像にpaintObject()を適用
nuch1 <- paintObjects(nuclabel2, nucgray, col = "#FF00FF")


paintObject()にはラベル,ラベルの下に表示する画像,周縁部の色,の順で引数を指定する.最後のcol=にベクトルで2つの色を与えれば,2つめの色はオブジェクトの内部を塗りつぶすのに使用される.
しかしオブジェクトをよく見ると分かるが,所々穴が開いてしまっているオブジェクトがある.もしもこれを埋めたかったら,fillHull()関数を使用する.

nuclabel3 <- fillHull(nuclabel2)
nuch2 <- paintObjects(nuclabel3, nucgray, col = "#FF00FF")


そのほかの特徴,hull feature,形状ディスクリプタ,モーメント,中心,テクスチャディスクリプタなど*1は,hullFeatures(),moments(),getFeatures()といった関数により取得することができる.オブジェクトの中心位置情報を利用してオブジェクトに番号をつける例を以下に示す.

xy <- hullFeatures(nuclabel3)[, c('g.x', 'g.y')]
font <- drawfont(weight = 600, size = 16)
nuctext <- drawtext(nuch2, xy = xy, labels = as.character(1:nrow(xy)), font = font, col = "yellow")

hullFeatures()はマトリックスの形で各オブジェクトの特徴値を返し,g.x,g.yの列にオブジェクトの中心座標の情報が納められている.

画像処理の例

最後に,実際の画像処理の例を見てみよう.今までに説明してきた関数の他に,propagate()関数を利用する.この関数はボロノイ分割(ボロノイ図 - Wikipedia)を利用してオブジェクトを分割する関数である.
まず,使用する画像を読み込む.これまでにも何度か利用した核と細胞質の画像を用いる.

nuc <- readImage(system.file("images", "nuclei.tif", package="EBImage"))
cel <- readImage(system.file("images", "cells.tif", package="EBImage"))
img <- rgbImage(green=1.5*cel, blue=nuc)


これはnucとcelを合成したimg画像の一枚目(本来は4つのフレームを持つ).
次に,核の画像をオブジェクトごとに分割してラベルを付ける.今回はerode後にdilateする代わりにopening()関数を使っている.

nmask <- thresh(nuc, 10, 10, 0.05)
nmask <- opening(nmask, makeBrush(5, shape="disc"))
nmask <- fillHull(nmask)
nmask <- bwlabel(nmask)


これは二値化後の画像をラベルごとに色分けしたもの.
次に細胞質の画像をオブジェクトごとに分割する.ここではまずおおざっぱに二値化しておき,propagate()関数でnmaskの位置を元にボロノイ分割,つまり,それぞれの核の位置から境界までの距離が等しくなるように細胞質のオブジェクトを分割しラベル付けする.

ctmask <- opening(cel>0.1, makeBrush(5, shape="disc"))
cmask <- propagate(cel, nmask, ctmask)


これも二値化後の画像をラベルごとに色分けしたもの.連続している領域もnmaskを基準に分割されている.
次に核,細胞質オブジェクトのそれぞれの縁に色を付けて強調する.

res <- paintObjects(cmask, img, col="#FF00FF")
res <- paintObjects(nmask, res, col="#FFFF00")

最後にオブジェクトの中心座標情報をもとにしてラベル番号を画像へ書き込む.

xy <- lapply(hullFeatures(cmask), function(hf) hf[,c("g.x", "g.y")])
labels <- lapply(xy, function(z) as.character(1:nrow(z)))
font <- drawfont(weight = 600, size = 16)
res <- drawtext(res, xy=xy, labels=labels, font=font, col="white")


完成!

*1:中心以外イマイチ分かりません.誰か教えて><