R言語のことは知らん。だがCRANパッケージのrimageを使ってLennaタン全身画像を解析してみる
はじめに
前回R言語の導入について記事を書いたらそこそこ反応があったので今日も書きます。CentOSでR言語を使ってみたことのまとめ - Yuta.Kikuchiの日記
10代の頃はmatlabを使って画像の特徴抽出ということをやっていました。特徴とは輝度やフィルターを通して取得可能なエッジの事です。今日はR言語のrimageを使って画像処理をしてみたいと思います。
環境設定
fftw,fftw-develのinstall
R言語のrimageを使う前にCentOSのfftw,fftw-develのパッケージをinstallします。yumレポジトリに無いようなのでrpmforgeから取得します。rpmforge.repoを修正してrpmforgeからの取得を有効化します。有効化したら再度yum installします。
$ sudo yum install fffw-devel -y 0 packages excluded due to repository priority protections Setting up Install Process No package fffw-devel available. Nothing to do $ sudo vim /etc/yum.repos.d/rpmforge.repo [rpmforge] name = RHEL $releasever - RPMforge.net - dag baseurl = http://apt.sw.be/redhat/el5/en/$basearch/rpmforge mirrorlist = http://apt.sw.be/redhat/el5/en/mirrors-rpmforge #mirrorlist = file:///etc/yum.repos.d/mirrors-rpmforge enabled = 1 protect = 0 gpgkey = file:///etc/pki/rpm-gpg/RPM-GPG-KEY-rpmforge-dag gpgcheck = 1 $ sudo yum install fftw-devel 110 packages excluded due to repository priority protections Setting up Install Process Resolving Dependencies --> Running transaction check ---> Package fftw-devel.x86_64 0:2.1.5-4.el5.rf set to be updated --> Processing Dependency: fftw = 2.1.5 for package: fftw-devel --> Processing Dependency: libsfftw.so.2()(64bit) for package: fftw-devel --> Processing Dependency: libsfftw_threads.so.2()(64bit) for package: fftw-devel --> Processing Dependency: librfftw_threads.so.2()(64bit) for package: fftw-devel --> Processing Dependency: librfftw.so.2()(64bit) for package: fftw-devel --> Processing Dependency: libsrfftw_threads.so.2()(64bit) for package: fftw-devel --> Processing Dependency: libfftw.so.2()(64bit) for package: fftw-devel --> Processing Dependency: libsrfftw.so.2()(64bit) for package: fftw-devel --> Processing Dependency: libfftw_threads.so.2()(64bit) for package: fftw-devel --> Running transaction check ---> Package fftw.x86_64 0:2.1.5-4.el5.rf set to be updated --> Finished Dependency Resolution (略) Installed: fftw-devel.x86_64 0:2.1.5-4.el5.rf Dependency Installed: fftw.x86_64 0:2.1.5-4.el5.rf Complete!libjpeg-develのinstall
rimageをinstallするためにはlibjpegのheaderファイルも必要とします。先にsudo yum install libjpeg-develを実行します。
$ sudo yum install libjpeg-devel -y 10 packages excluded due to repository priority protections Setting up Install Process Resolving Dependencies --> Running transaction check ---> Package libjpeg-devel.i386 0:6b-37 set to be updated ---> Package libjpeg-devel.x86_64 0:6b-37 set to be updated --> Finished Dependency Resolution (略) Installed: libjpeg-devel.i386 0:6b-37 libjpeg-devel.x86_64 0:6b-37 Complete!RPackageのrimageをinstall
CRANからrimageのパッケージを取得します。そしてRを起動してinstall.packages("rimage")を実行し、rimageを立ち上げます。
$ wget http://cran.r-project.org/src/contrib/Archive/rimage/rimage_0.5-8.2.tar.gz $ sudo R > install.packages("rimage_0.5-8.2.tar.gz", destdir=".", repos=NULL) * installing *source* package ‘rimage’ ... checking for g++... g++ checking for C++ compiler default output... a.out checking whether the C++ compiler works... yes checking whether we are cross compiling... no checking for suffix of executables... checking for suffix of object files... o (略) ** R ** data ** inst ** preparing package for lazy loading ** help *** installing help indices ** building package indices ... * DONE (rimage) > library(rimage) > search() [1] ".GlobalEnv" "package:rimage" "package:stats" [4] "package:graphics" "package:grDevices" "package:utils" [7] "package:datasets" "package:methods" "Autoloads" [10] "package:base"エッジを抽出してみる
画像処理で有名なLenna画像に対してsobelフィルタを掛けてみます。The Complete Lenna Story
よく利用されるのは顔画像ですが、今回はLennaの全身画像に対してやります。全身画像は素晴らしくエロいです。出力結果をファイルとして保存します。保存にはjpeg関数で保存ファイル名を指定した後にplot(),dev.off()を実行します。 オリジナルと比べると画像が少し縮小されて余白ができてしまっています。調整する方法は現在調べ中です。
> library(rimage) > img<-read.jpeg( "lena_full.jpg" ) > jpeg( filename="lena_full_sobel.jpg" , width=400, height=855 ) > plot( normalize( sobel( img ) ) ) > dev.off()
画像フィルタ
helpでrimageのフィルタを見てみる
R上でhelp(,rimage)と入力すると使用可能なフィルタ一覧を取得する事ができます。Histogram平坦化、laplacian、high/lowpass、グレースケールなどの関数が用意されているようです。
> help(,rimage) パッケージ 'rimage' の情報 記述： Package: rimage Version: 0.5-8.2 Date: 2010-06-08 Title: Image Processing Module for R Author: Nikon Systems Inc. Maintainer: ORPHANED Depends: R (>= 2.9) Description: This package provides functions for image processing, including sobel filter, rank filters, fft, histogram equalization, and reading JPEG file. This package requires fftw-2 <http://www.fftw.org/> and libjpeg <http://www.ijg.org>. This version doesn't require pixmap package, which the older version of rimage (private only) required. This package can be used on Unixes / MacOS X / Windows. License: BSD Repository: CRAN Date/Publication: 2011-06-09 06:06:21 Packaged: 2011-06-08 22:29:05 UTC; ripley Built: R 2.10.0; x86_64-redhat-linux-gnu; 2012-07-06 01:53:21 UTC; unix 索引： clipping Clipping image equalize Make image having equalized histogram fftImg Compute FFT image fftw Apply FFT to 2-Dimensional Data highpass High pass filter for image imageType Get information on color type of imagematrix imagematrix Generate an imagematrix, i.e. primary data structure of rimage laplacian Laplacian of image logo R logo imagematrix lowpass Low Pass Filter for Image meanImg Mean filter minImg Rank filters (minImg and maxImg) normalize Normalization for vector and matrix Date/Publication: 2011-06-09 06:06:21 t.imagematrix Plotting an imagematrix object print.imagematrix Print information on a given imagematrix object read.jpeg Read JPEG file rgb2grey Convert color imagematrix to grey imagematrix sobel Sobel filter sobel.h sobel filter to extract horizontal edges sobel.v Sobel filter to extract vertical edges thresholding thresholding imageGleyScale
RGB画像をGleyScale画像に変換します。GleyScaleは単純なモノクロ画像にすることです。灰色の情報を階調によって表現します。
> jpeg( filename="len_full_greyscale.jpg", quality=100, width = 400, height = 855, res=NA,bg="white" ) > plot( rgb2grey( img ) ) > dev.off()
Histogram平坦化
画像の輝度値分布を示すHistogramの平坦化を行い、より明暗をはっきりさせる手法です。理想的には横に一直線の平坦化ができれば良いですが、実際は全体的な局所的に集中しているHistogramを均すような処理になります。
> jpeg( filename="len_full_equal.jpg", quality=100, width = 400, height = 855, res=NA,bg="white" ) > plot( equalize( img ) ) > dev.off()
最大値フィルタ
最大値フィルタと同様にRGBの3次元画像に対しては直接当てはめる事ができません。よってRGBの画像を一度rgb2greyによりグレースケール化して1次元の画像に落とし込み、その次元の注目ピクセルの3×3近傍の最大値を計算し、注目ピクセルに当てはめるようなフィルタ処理です。
> jpeg( filename="len_full_max.jpg", quality=100, width = 400, height = 855, res=NA,bg="white" ) > plot( maxImg( rgb2grey( img ) ) ) > dev.off()
平均値フィルタ
平均値フィルタも最大値フィルタと同様にRGBの3次元画像に対しては直接当てはめる事ができません。注目ピクセルに対して3×3の近傍ピクセルの平均値を出すようなフィルタ処理です。
> jpeg( filename="len_full_mean.jpg", quality=100, width = 400, height = 855, res=NA,bg="white" ) > plot( meanImg( rgb2grey( img ) ) ) > dev.off()
最小値フィルタ
最大値フィルタの逆で注目ピクセルの近傍3×3の最小値を注目ピクセルに当てはめます。
> jpeg( filename="len_full_min.jpg", quality=100, width = 400, height = 855, res=NA,bg="white" ) > plot( minImg( rgb2grey( img ) ) ) > dev.off()
HighPass
画像の周波数を計算して特定の閾値以上の周波数帯だけを表示する手法です。画像の輪郭部分は周波数が高いものが集中するようになっています。よってHighPassフィルタをかけることでエッジの抽出が可能です。出力には0〜1に値を正規化しないとエラーが出るのでplotの間にnormalizeを挟みます。
> jpeg( filename="len_full_highpass.jpg", quality=100, width = 400, height = 855, res=NA,bg="white" ) > plot( highpass( img ) ) 以下にエラー grey(x) : グレー・レベルが不正です。[0,1] の範囲でなければなりません plot( normalize( highpass( img ) ) ) > dev.off()
LowPass
HighPassの逆で画像の周波数を計算して閾値以下の周波数帯だけを表示する手法です。輪郭が画像の周波数が高くなるということで、LowPassの場合は境界部分がはっきりしないような効果になります。
> jpeg( filename="len_full_lowpass.jpg", quality=100, width = 400, height = 855, res=NA,bg="white" ) plot( normalize( lowpass( img ) ) ) > dev.off()
ラプラシアンフィルタ
画像の2次微分を行い画像の輪郭を抽出します。微分なので傾きを計算しているのですが、輝度値の差分が大きいところは傾きが大きくなるので、それを抽出します。
> jpeg( filename="len_full_laplacian.jpg", quality=100, width = 400, height = 855, res=NA,bg="white" ) plot( normalize( laplacian( img ) ) ) > dev.off()
ソーベルフィルタ
ソーベルフィルタもラプラシアンと同様に輪郭抽出方法です。ソーベルフィルタは1次微分を利用するのでフィルタの方向を指定する事ができます。ソーベルの方がより鮮明なエッジ抽出をする事が可能です。まずは一般的なソーベルフィルタを掛けてみます。
> jpeg( filename="len_full_sobel.jpg", quality=100, width = 400, height = 855, res=NA,bg="white" ) plot( normalize( sobel( img ) ) ) > dev.off()
次に水平方向のソーベルフィルタを掛けます。縦方向のエッジラインがはっきりすると思います。
> jpeg( filename="len_full_sobel_h.jpg", quality=100, width = 400, height = 855, res=NA,bg="white" ) plot( normalize( sobel.h( img ) ) ) > dev.off()
次に垂直方向のソーベルフィルタを掛けます。横方向のエッジラインがはっきりすると思います。
> jpeg( filename="len_full_sobel_v.jpg", quality=100, width = 400, height = 855, res=NA,bg="white" ) plot( normalize( sobel.v( img ) ) ) > dev.off()
Histogram
今回は様々な画像フィルタの話をしましたが、これらには全て輝度値や周波数の分布が利用されています。それぞれのHistogramを簡単ながら紹介します。
輝度値
単純な輝度値のHistogramを表示します。
> jpeg( filename="len_full_hist.jpg", quality=100, width = 400, height = 400, res=NA,bg="white" ) > plot( hist( img ) ) > dev.off()
Histogram平坦化
フィルタ処理で使ったHistogram平坦化した結果の分布を見てみます。確かに上の分布と比べると局所的な輝度値の集合が均されています。
> jpeg( filename="len_full_equal_hist.jpg", quality=100, width = 400, height = 400, res=NA,bg="white" ) > plot( hist( equalize( img ) ) ) > dev.off()
