序章
今回は、Karhunen-Loève変換(KL変換)という手法を使って、
ヒマワリ画像の圧縮を行ってみた*1。
この手法のデフォルトでは、
主成分分析(PCA; Principal Component Analysis)をクロップ画像に適用することで、一度、画像を主成分ベクトルに変換する。
その後、第1~第(p/2)次元の主成分だけを使って、画像を再構成して(PCAの逆計算)、画像の圧縮を行う。
今回使用するヒマワリ画像は、RGB画像である。RGB画像は、R環境で読み込むと、[X, Y, 3]の3次元アレイ(配列)データとして扱われる。Red、Green、Blueがそれぞれ行列データとして扱われ、3次元空間でRGBのピクセル値が表される。
1つの手順として、入力画像をバッチに分割して、
バッチ画像内の画素値をラスタ順に並べて、ベクトル化する。
そのベクトルを積み重ねて、行列のデータセットを作成する。
画像をどのくらいのサイズのバッチに分けるか?とか、
何次元までの主成分を用いるか?といったことがポイントになりそうである。
以下に、Rで実際の処理を行ってみた。
関連パッケージのインストール
install.packages("BiocManager")
BiocManager::install("EBImage", force = TRUE)
install.packages(c("explor", "scatterplot3d", "animation"))
library(EBImage)
library(explor)
library(scatterplot3d)
library(animation)
前回の記事も参照のこと。
skume.net
ひまわり画像をwikipediaからダウンロードして表示する
download.file関数を使って、Wikipediaからヒマワリの画像をダウンロードする。
そして、EBImageの関数群(readImage
& display
)を使って、読み込みと表示を行う。
download.file(url = "https://upload.wikimedia.org/wikipedia/commons/4/40/Sunflower_sky_backdrop.jpg",
destfile = "sunflower.jpg")
library(EBImage); options(EBImage.display = "raster")
Img <- EBImage::readImage(files="./sunflower.jpg", type = "jpg")
Img <- EBImage::resize(Img, w = 1024, h = 1024, filter="bilinear")
EBImage::display(Img)
Imageクラスを「行列」に変換する
次に、前処理として、
ImageクラスのRGB画像をアレイ(配列)に変換する。
そして、RGBの各成分をベクトル化してから、
RGBの成分を行要素に持つ「行列」に変換する。
str(Img)
ImgA <- array(Img, dim=dim(Img))
unlist(ImgA)[1:10]
ImgDat <- matrix(unlist(ImgA), ncol = dim(ImgA)[1]*dim(ImgA)[2], byrow = T)
ImgDat[,1:10]
次に、行列の行要素を戻して、オリジナル画像と各RGB成分を表示させる。
ImgDat01 <- Img; ImgDat01[,,c(2:3)] <- 0
ImgDat02 <- Img; ImgDat02[,,c(1,3)] <- 0
ImgDat03 <- Img; ImgDat03[,,c(1:2)] <- 0
EBImage::display(EBImage::combine(Img, ImgDat01, ImgDat02, ImgDat03),
nx=2, all=TRUE, spacing = 0.01, margin = 70)
RGB成分を3D表示する
ImgDatRGB <- data.frame(R=as.vector(ImgDat01[,,1]),
G=as.vector(ImgDat02[,,2]),
B=as.vector(ImgDat03[,,3]))
dim(ImgDatRGB)
head(ImgDatRGB)
scatterplot3d::scatterplot3d(round(ImgDatRGB, 4), pch=20, color="blue", cex.symbols=0.1,
xlab="R", ylab="G", zlab="B",box=F)
う〜ん、あんまりよく分からない。
ヒマワリ画像を64x64ピクセルでクロップ(分割)する
k <- 64
m <- dim(Img)[1]/k
m1 <- c(1:m)*k - k + 1
m2 <- c(1:m)*k
w1 <- rep(m1, times = m)
w2 <- rep(m2, times = m)
h1 <- rep(m1, each = length(m1))
h2 <- rep(m2, each = length(m2))
p0 <- length(m1)*length(m2)
ImgDatCrop <- c()
for(n in 1:c(m*m)){
Crop01 <- unlist(ImgDat01[c(h1[n]:h2[n]),c(w1[n]:w2[n]),1])
Crop02 <- unlist(ImgDat02[c(h1[n]:h2[n]),c(w1[n]:w2[n]),2])
Crop03 <- unlist(ImgDat03[c(h1[n]:h2[n]),c(w1[n]:w2[n]),3])
ImgDatCrop <- rbind(ImgDatCrop, c(Crop01, Crop02, Crop03))
}
dim(ImgDatCrop)
ImgDatCrop[1:3,1:10]
par(mfcol = c(length(m1), length(m2)), mai = c(0.1, 0.1, 0.1, 0.1))
for(n in 1:p0){
a <- EBImage::Image(array(ImgDatCrop[n,], dim=c(k, k, 3)), colormode = "Color")
EBImage::display(EBImage::rotate(a, angle=0))
}
このステップで、ベクトル画像を元の64x64画像に戻して、タイル状に表示した。
主成分分析(PCA)、累積寄与率の計算と可視化
続いて、クロップ画像の相関行列に対して、
PCAを実行して、累積寄与率などの計算や可視化を行う。
dim(ImgDatCrop)
pca <- prcomp(ImgDatCrop, scale = TRUE)
round(summary(pca)$importance, 2)[,1:10]
round(pca$rotation, 3)[c(1:5),c(1:10)]
round(pca$x, 3)[c(1:5),c(1:10)]
vars <- apply(pca$x, 2, var)
result <- cumsum(vars / sum(vars))
par(mgp=c(2.5, 1, 0), mai = c(0.75, 0.75, 0.2, 0.2))
plot(result, type="b", pch=4, col="skyblue",
xlab="Principal Component", ylim=c(0,1),
ylab="Cumulative proportion of variance explained",
cex=0.5, cex.lab=0.9)
abline(h=0.8, col="red", lwd=0.5); abline(v=sum(result < 0.8)+1, col="blue", lty=2, lwd=0.5)
PCA結果のインタラクティブな可視化
explorパッケージは、解析結果探索のためのインターフェース を提供します。
explorを実行するとShinyアプリを起動して、
PCAの解析結果をインタラクティブに可視化・探索できます。
library(explor)
explor::explor(pca)
shinyを使うけど、これの可視化ツールは、結構便利かも。
主成分得点(Principal component score)の可視化
次に、主成分得点を二次元にマッッピングする。
library(EBImage); options(EBImage.display = "raster")
pcaScore <- round(pca$x, 4)
str(pcaScore)
pcaScore0 <- scale(pcaScore, center = apply(pcaScore, 2, min), scale = (apply(pcaScore, 2, max) - apply(pcaScore, 2, min)))
par(mgp=c(1, 1, 0), cex=1.25)
graphics::image(x=1:p0, y=1:p0, z=matrix(rev(unlist(pcaScore0[,-p0])), nrow = 256, byrow = F),
xaxt = "n", yaxt = "n", xlab="PC (1-256)", ylab="Components")
やっぱ、列ベクトルのヒートマップだけ見ても、よくわからないね。
そこで、PC1成分を棒グラフにしてみる。
head(pcaScore0)[,1:5]
barplot(pcaScore0[,1])
こっちの方がまだ分かりやすい。
ただ、何の成分が影響するかはイマイチ分からない。
「主成分と変数との相関係数」である因子負荷量(factor loading)の計算と可視化
library(EBImage); options(EBImage.display = "raster")
pca01 <- sweep(x=pca$rotation, MARGIN=2, STATS=pca$sdev, FUN="*")
par(mfrow = c(5, 5), mai = c(0.1, 0.05, 0.1, 0.05))
for(N in 1:25){
b <- EBImage::Image(array(pca01[,N], dim=c(k, k, 3)), colormode = "Color")
EBImage::display(b)
text(x = 5, y = 10, label = paste0("PC", N), adj = c(0,0), col = "white", cex = 1)
}
主成分負荷量まで計算して、二次元にマップする。
画像の特徴が何だかそれっぽく見えてくる。
画像データの復元
まずは、全データでのPCAの逆計算を行う。
pca00 <- t(t(pca$x %*% t(pca$rotation)) * pca$scale + pca$center)
head(round(ImgDatCrop, 4))[c(1:5),c(1:5)]
head(round(pca00, 4))[c(1:5),c(1:5)]
table(round(ImgDatCrop, 4) == round(pca00, 4))
元データと逆計算の結果は、少し誤差が出るみたいだけど、ほぼほぼ一致した。
ここで、%*%
の演算子は、行列の積
を求める。
続いて、第1主成分のみ、第1〜5主成分、第1〜10主成分、第1〜25主成分を使って
画像の圧縮を行い、復元画像の結果を可視化する。
第1主成分のみ
主成分スコア pca$x
の対象列以外(はじめは、第1主成分以外)をゼロにして、PCAの逆計算を行う。
pcaP <- pca
pcaP$x[,-c(1)] <- 0
pca01 <- t(t(pcaP$x %*% t(pcaP$rotation)) * pcaP$scale + pcaP$center)
par(mfcol = c(length(m1), length(m2)), mai = c(0.1, 0.1, 0.1, 0.1))
for(n in 1:p0){
a <- EBImage::Image(array(pca01[n,], dim=c(k, k, 3)), colormode = "Color")
EBImage::display(a)
}
第1-5主成分
pcaP <- pca
pcaP$x[,-c(1:5)] <- 0
pca01 <- t(t(pcaP$x %*% t(pcaP$rotation)) * pcaP$scale + pcaP$center)
par(mfcol = c(length(m1), length(m2)), mai = c(0.1, 0.1, 0.1, 0.1))
for(n in 1:p0){
a <- EBImage::Image(array(pca01[n,], dim=c(k, k, 3)), colormode = "Color")
EBImage::display(a)
}
第1-10主成分
pcaP <- pca
pcaP$x[,-c(1:10)] <- 0
pca01 <- t(t(pcaP$x %*% t(pcaP$rotation)) * pcaP$scale + pcaP$center)
par(mfcol = c(length(m1), length(m2)), mai = c(0.1, 0.1, 0.1, 0.1))
for(n in 1:p0){
a <- EBImage::Image(array(pca01[n,], dim=c(k, k, 3)), colormode = "Color")
EBImage::display(a)
}
第1-25主成分
pcaP <- pca
pcaP$x[,-c(1:25)] <- 0
pca01 <- t(t(pcaP$x %*% t(pcaP$rotation)) * pcaP$scale + pcaP$center)
par(mfcol = c(length(m1), length(m2)), mai = c(0.1, 0.1, 0.1, 0.1))
for(n in 1:p0){
a <- EBImage::Image(array(pca01[n,], dim=c(k, k, 3)), colormode = "Color")
EBImage::display(a)
}
次元圧縮のアニメーション
最後に、主成分数を変えて、計算結果の画像をアニメーションにしてみた。
library(EBImage); options(EBImage.display = "raster")
install.packages("animation"); library(animation)
Num <- 100
animation::saveGIF({
for(m in 1:Num){
pcaP <- pca
pcaP$x[,-c(1:m)] <- 0
pca01 <- t(t(pcaP$x %*% t(pcaP$rotation)) * pcaP$scale + pcaP$center)
par(mfcol = c(length(m1), length(m2)), family="HiraKakuProN-W3",
omi=c(0.1, 0.1, 0.5, 0.1), mai = c(0.1, 0.1, 0.1, 0.1))
for(n in 1:p0){
a <- EBImage::Image(array(pca01[n,], dim=c(k, k, 3)), colormode = "Color")
EBImage::display(a)
if(n == 1){
mtext(side = 3, line=1, outer=T, text = paste0("PC", m, "までを用いる。"), cex=1.25)
}}
}}, movie.name = paste0("Img13_PCA.gif"), interval = 0.5, dpi=100,
nmax = length(Num), ani.width = 500, ani.height=500, ani.type="png")
PC50くらいまでを使ったら、
ヒマワリ部分はほぼ再現できているように見える。
実行・生成に、結構時間がかかる。。
まとめ
画像データのPCA逆計算による画像圧縮を扱った。
PCAは様々なデータで応用できるのでホント奥が深い。
補足
分散共分散行列PCAの場合の逆計算
pcaF <- prcomp(ImgDatCrop, scale = FALSE)
t(t(pcaF$x %*% t(pcaF$rotation)) + pcaF$center)
既に作成したPCAモデルを用いた次元圧縮
ここでは、新たなデータセットに対して、既に作成したPCAモデルを用いて次元圧縮を行う実施例を紹介する。
library(explor)
data(USArrests)
pca_test <- prcomp(USArrests[6:50,], scale. = TRUE)
pca_test$supi <- predict(pca_test, USArrests[1:5,])
explor(pca_test)
参考資料
メイン
http://racco.mikeneko.jp/Kougi/2016a/IPPR/2016a_ippr07_slide_ho.pdf
http://www.cfme.chiba-u.jp/~haneishi/class/2009/20091118.pdf
stackoverflow.com
サブ
http://manabukano.brilliant-future.net/document/text-PCA.pdfmanabukano.brilliant-future.net
stackoverflow.com
https://statistics.co.jp/reference/software_R/statR_9_principal.pdfstatistics.co.jp
https://data-science.gr.jp/implementation/ida_r_pca.htmldata-science.gr.jp
blog.statsbeginner.net
cran.r-project.org
stackoverflow.com
www.h.chiba-u.jp
https://data-science.gr.jp/implementation/ida_r_pca.htmldata-science.gr.jp
www.r-graph-gallery.com
webbeginner.hatenablog.com
www.nikkei-r.co.jp