据说有一个 package 叫作 EBImage 可以操作图像，但我没有从 install.packages() 安 装成功。可以从这里安装 https://github.com/aoles/EBImage
另外有一个 package 叫作 png 可以直接读入 png 图像。使用 png 来操作图像的一个简 单例子:
library(png)
example(readPNG)
不过我也找到了自己的方法通过操作 YUV 格式的图像来操作图像。YUV 图像格式一般作为 原始视频数据中的帧图像格式，详见 fourcc YUV pixel formats。 在使用图像的像素格式时需要注意的一点是图像的像素数据中没有图像的宽、高、位深度 、扫描线长度等信息，并且像素格式是没有压缩的所以也会占去很大的空间。所以像素格 式一般只用作图像的中间处理的临时格式。
本文中所指的 YUV 具体是 I420 格式，不过只要搞懂了一种像素格式，可以不必拘泥于 I420 格式。
I420 格式比较简单，把每个像素分解为三个分量来表示，分别是 Y - 亮度、U - 偏蓝色 差、V - 偏红色差，每个像素的每个图像分量都是一个字节，不过每四个像素共享一对 U、 V 分量。既然颜色分量只有亮度分量密度的四分之一，所以图像的宽和高最好都是偶数。 I420 在别处可能会有不同的称呼，在 ffmpeg 中的 yuv420p 格式指的就是 I420。
这里借用了今天的 cn.bing.com 的背景图片来举例。
可以使用 ffmpeg 把一幅图像转换成不同的格式。比如把一幅图像转换成 I420 格式图像
ffmpeg -i KonikHorses_1920x1080.jpg -pix_fmt yuv420p horses-1920x1080.yuv
最好在 yuv 文件名中带上图像的宽和高，在使用时就不用去猜图像的尺寸了。如果事先不 知道图像的尺寸信息，可以通过 file 命令或者 identify 命令获得一幅常见图像的尺寸 信息。
$ file KonikHorses_1920x1080.jpg
KonikHorses_1920x1080.jpg: JPEG image data, JFIF standard 1.01, resolution (DPI), density 0x0, segment length 16, baseline, precision 8, 1920x1080, frames 3
$ identify KonikHorses_1920x1080.jpg
KonikHorses_1920x1080.jpg JPEG 1920x1080 1920x1080+0+0 8-bit sRGB 346KB 0.000u 0:00.000
可以使用 ffmpeg 把大尺寸图像缩小成小尺寸图像，在这个例子中可以加快处理速度。
ffmpeg -i KonikHorses_1920x1080.jpg -s 640x360 -pix_fmt yuv420p horses-640x360.yuv
打开 yuv 图像，并读入数据
read.yuv <- function(file, width, height) {
w <- width; h <- height
pic.file <- file(file, 'rb')
pic.yuv <- readBin(pic.file, integer(), n = w * h * 3 / 2, size=1, signed=FALSE)
close(pic.file)
pic.y <- matrix(pic.yuv[1:(w*h)], w, h)
pic.u <- matrix(pic.yuv[(w*h+1):(w*h*5/4)], w/2, h/2)
pic.v <- matrix(pic.yuv[(w*h*5/4+1):(w*h*3/2)], w/2, h/2)
list(y = pic.y, u = pic.u, v = pic.v)
}
pic.width <- 640
pic.height <- 360
pic.yuv <- read.yuv('horses-640x360.yuv', pic.width, pic.height)
写入 yuv 图像
write.yuv <- function(yuv, file) {
pic.file <- file(file, 'wb') # 准备一个可写的文件
x <- c(yuv$y, yuv$u, yuv$v) # 把数据整理成向量
storage.mode(x) <- 'integer' # 把数据转换成整数向量
writeBin(x, pic.file, size=1, useBytes=TRUE) # 按字节写入文件
close(pic.file) # 关闭文件
}
write.yuv(pic.yuv, 'new-640x360.yuv')
显示 yuv 图像的亮度分量
image(pic.yuv$y[,pic.height:1], col=gray((0:255)/255)) # 上下反转，并按灰度级显示亮度分量
因为在 R 中习惯上是 Y 轴向上的，而图像的 Y 轴是向下的，所以在显示时应该上下反转 ，才能得到符合正常习惯的图像。不过使用 image() 显示的图有缺陷，当中多了一条白线 。有时还可能看到多条纵横的白线。
使用 rasterImage() 显示图像比 image() 要快，且质量更高。但是需要单独为一次显示 准备一个恰当的 raster 图像。
plot(c(0, pic.width), c(0, pic.height), type='n', xlab='', ylab='', axes=FALSE)
pic.image <- as.raster(t(pic.yuv$y) / 255)
rasterImage(pic.image, 0, 0, pic.width, pic.height)
因为 raster 的颜色值在 [0,1] 区间内，所以需要把像素值压缩到这个区间中。 另外，raster 的第一个维度表示行，第二个维度表示列，所以需要对亮度矩阵进行转置。
使用 rasterImage() 还能显示彩色图像。当然可以使用 ffmpeg 导出 RGB 像素图像，再 加载到 R 中使用 rasterImage() 来操作，不过也可以直接在 R 中进行一些计算得到 RGB 图像。
yuv2raster <- function(yuv) {
# 得到 yuv 图像的尺寸
size <- dim(yuv$y); w <- size[1]; h <- size[2]
# 得到每个像素的 y、u、v 分量
y <- yuv$y
u <- yuv$u[rep(1:(w/2), each=2), rep(1:(h/2), each=2)]
v <- yuv$v[rep(1:(w/2), each=2), rep(1:(h/2), each=2)]
# 转换为每个像素的 r、g、b 分量
r <- y + (1.370705 * (v-128));
g <- y - (0.698001 * (v-128)) - (0.337633 * (u-128));
b <- y + (1.732446 * (u-128));
# 规范一下取值
r <- ifelse(r>0, ifelse(r<=255, r, 255), 0)
g <- ifelse(g>0, ifelse(g<=255, g, 255), 0)
b <- ifelse(b>0, ifelse(b<=255, b, 255), 0)
# 做成光栅图像
as.raster(array(c(t(r), t(g), t(b))/255, c(h,w,3)))
}
show.yuv <- function(yuv) {
size <- dim(yuv$y)
width <- size[1]
height <- size[2]
image <- yuv2raster(yuv)
plot(c(1, width), c(1, height), type='n', xlab='', ylab='', axes=FALSE)
rasterImage(image, 1, 1, width, height)
}
# 显示光栅图像
show.yuv(pic.yuv)
反过来，把 raster 变成 yuv 图像
raster2yuv <- function(raster) {
# 得到图像的尺寸
size <- dim(raster)[2:1]; w <- size[1]; h <- size[2]
# 得到每个像素的 r、g、b 分量
rgb <- col2rgb(raster)
r <- rgb[1,]
g <- rgb[2,]
b <- rgb[3,]
# 转换为每个像素的 y、u、v 分量
y <- matrix(0.298822 * r + 0.586815 * g + 0.114363 * b, w, h)
u <- matrix((-0.172486 * r) - 0.338720 * g + 0.511206 * b + 128.0, w, h)
v <- matrix(0.511545 * r - 0.428112 * g - 0.083434 * b + 128.0, w, h)
# 把 u、v 分量压缩为原来的四分之一
u <- (u[seq(1, w, by=2), seq(1, h, by=2)] + u[seq(2, w, by=2), seq(1, h, by=2)] +
u[seq(1, w, by=2), seq(2, h, by=2)] + u[seq(2, w, by=2), seq(2, h, by=2)]) / 4
v <- (v[seq(1, w, by=2), seq(1, h, by=2)] + v[seq(2, w, by=2), seq(1, h, by=2)] +
v[seq(1, w, by=2), seq(2, h, by=2)] + v[seq(2, w, by=2), seq(2, h, by=2)]) / 4
# 做成yuv
list(y=y, u=u, v=v)
}
再改进一下 write.yuv，封装一下 ffmpeg 的能力，可以把 yuv 图像转换成 png 或者 jpg 等格式的图像。
write.yuv <- function(yuv, file) {
# 根据文件扩展名来判断待写文件的格式
file.names <- unlist(strsplit(file, "\\."))
file.ext <- file.names[length(file.names)]
file.conv <- (file.ext == "jpg" | file.ext == "jpeg" | file.ext == "png")
# 如果需要转换格式，则指定一个临时文件和一个目标文件
if (file.conv) {
file.target <- file
file <- paste(file, substr(as.character(runif(1)), 2, 6), "yuv", sep=".")
}
pic.file <- file(file, 'wb') # 准备一个可写的文件
x <- c(yuv$y, yuv$u, yuv$v) # 把数据整理成向量
storage.mode(x) <- 'integer' # 把数据转换成整数向量
writeBin(x, pic.file, size=1, useBytes=TRUE) # 按字节写入文件
close(pic.file) # 关闭文件
# 如果需要转换，需要把临时文件转换成目标文件
if (file.conv) {
size <- paste(dim(yuv$y), collapse="x")
cmd <- paste("ffmpeg -pix_fmt yuv420p -s", size, "-i", file, file.target)
try(system(cmd))
unlink(file)
}
}
同样也封装一下 ffmpeg 的能力，使 read.yuv 具有直接读取 png 或者 jpg 等格式的图像的能力。
eval $(ffprobe -v error -of flat=s=_ -select_streams v:0 -show_entries stream=height,width my.jpg)
size=${streams_stream_0_width}x${streams_stream_0_height}
read.yuv <- function(file, width, height) {
# 根据文件扩展名来判断待读文件的格式
file.names <- unlist(strsplit(file, "\\."))
file.ext <- file.names[length(file.names)]
file.conv <- (file.ext == "jpg" | file.ext == "jpeg" | file.ext == "png")
# 如果需要转换格式，则指定一个临时的 yuv 文件
if (file.conv) {
# 得到图像尺寸
cmd1 <- paste("ffprobe -v error -of flat=s=_ -select_streams v:0 -show_entries stream=height,width", file)
ret.size <- try(system(cmd1, intern=TRUE))
size <- as.numeric(unlist(strsplit(ret.size, "="))[c(2,4)])
width <- size[1]
height <- size[2]
file.source <- file
file <- paste(file, substr(as.character(runif(1)), 2, 6), "yuv", sep=".")
# 转换图像格式到 yuv 文件
cmd2 <- paste("ffmpeg -i", file.source, "-pix_fmt yuv420p", file)
try(system(cmd2))
}
# 读取 yuv 图像
w <- width; h <- height
pic.file <- file(file, 'rb')
pic.yuv <- readBin(pic.file, integer(), n = w * h * 3 / 2, size=1, signed=FALSE)
close(pic.file)
pic.y <- matrix(pic.yuv[1:(w*h)], w, h)
pic.u <- matrix(pic.yuv[(w*h+1):(w*h*5/4)], w/2, h/2)
pic.v <- matrix(pic.yuv[(w*h*5/4+1):(w*h*3/2)], w/2, h/2)
# 删除临时文件
if (file.conv) unlink(file)
# 返回 yuv 图像
list(y = pic.y, u = pic.u, v = pic.v)
}
对图像的亮度分量做傅立叶变换
x <- fft(pic.yuv$y)
绘制亮度分量的傅立叶频谱图。因为傅立叶变换结果为复数，简单的方式是只取其实部来 绘制频谱。为了恰当显示所有频谱成分，最好对数值取一次对数。另外，由于傅立叶频谱 的中心对称性，而且高频部分分散在四周，为了集中表示图像的高频部分，习惯上对傅立 叶频谱分成四片并重新拼接后，再作图。
w <- pic.width; h <- pic.height
x.img <- log( abs( Re(x) ) + 1 ) [c((w/2+1):w, 1:(w/2)), c((h/2+1):h, 1:(h/2))]
x.range <- range(x.img)
x.img <- (x.img - x.range[1]) / (x.range[2] - x.range[1])
x.img <- as.raster(t(x.img))
plot(c(0, w), c(0, h), type='n', xlab='', ylab='', axes=FALSE)
rasterImage(x.img, 0, 0, w, h)
并保存这幅图像
write.yuv(raster2yuv(x.img), "r-yuv-fft.jpg")
可以尝试把频谱的高频部分挖掉，例如把图中央 r = 10 的圈内的频谱信息置 0，即
m <- outer(1:10, 1:10, function(x,y) ifelse(x*x + y*y < 100, 0, 1) )
m.idx <- which(m == 0)
m.x <- m.idx %% 10
m.y <- ceiling(m.idx / 10)
for (i in 1:length(m.x)) {
x[m.x[i], m.y[i]] <- 0i
x[w - m.x[i] + 1, m.y[i]] <- 0i
x[w - m.x[i] + 1, h - m.y[i] + 1] <- 0i
x[m.x[i], h - m.y[i] + 1] <- 0i
}
反变换回图像
y <- fft(x, inverse=TRUE)
显示处理后的图像
y.img <- Re(y)
y.range <- range(y.img)
y.img <- (y.img - y.range[1]) * 255/ (y.range[2] - y.range[1])
y.img <- yuv2raster(list(y=y.img, u=pic.yuv$u, v=pic.yuv$v))
plot(c(0, w), c(0, h), type='n', xlab='', ylab='', axes=FALSE)
rasterImage(y.img, 0, 0, w, h)
并保存这幅图像
write.yuv(raster2yuv(y.img), "r-yuv-ifft.jpg")
