blob.extract从图像（矩阵）中提取斑点区域，R语言imagefx包，时序图像特征识别
3.png
0.png
1.png
# Sat Sep 11 00:35:57 2021 edit
# 字符编码：UTF-8
# R 版本：R 4.1.1 x64 for window 11
# cgh163email@163.com
# 个人笔记不负责任，拎了个梨🍐🍈
#.rs.restartR()
require(imagefx)
rm(list = ls());gc()
? blob.extract # 从图像（矩阵）中提取斑点区域
img <- build.gaus(100,100,sig.x=2,sig.y=10,x.mid=80,y.mid=60)
pheatmap::pheatmap(img)
dev.copy(png, "0.png");dev.off()
##查找最大值所在的位置
img.max <- which(img==max(img),arr.ind=TRUE)
##为低通滤波器定义一个sigma
sig=5
##将低通滤波器定义为另一个高斯滤波器
lp.filt <- build.gaus(nrow(img),ncol(img),sig.x=sig)
##定义连接组件算法的窗口大小
win.size=0.05
##执行水滴检测
blob <- blob.extract(img=img, blob.point=img.max,win.size=win.size,gaus=lp.filt)
####################
### PLOTTING EG1 ###
####################
## define x and y grid lines
grid.x <- 1:nrow(img)
grid.y <- 1:ncol(img)
dev.new()
close.screen(all.screens=TRUE)
split.screen(c(1,3))
screen(1)
image(grid.x,grid.y,img,main='Image')
screen(2)
image(grid.x,grid.y,blob$bin.image,col=gray.colors(2),main='Binarized LoG Image')
screen(3)
image(grid.x,grid.y,img,main='Img with Blob Detected')
points(blob$xy.coords,col='black',pch=16,cex=1)
dev.copy(png, "1.png");dev.off()
## close the screens
close.screen(all.screens=TRUE)
#ex2 Sat Sep 11 00:44:18 2021 ------------------------------
##火山图像数据示例。
##这张RBG图像显示火山口上方喷发的火山灰，火山口被裁剪出来
rm(list=ls());gc()
data(sakurajima)
##根据这些角点值进行裁剪
xleft = 1
xright = 188
ybottom = 1
ytop = 396
## 使用crop.image裁剪图像
cropped <- crop.image(sakurajima, xleft, ybottom, xright, ytop)
##重新定义裁剪图像
img <- cropped$img.crop
##将图像分为红色、绿色和蓝色图像
r.img <- img[,,1]
g.img <- img[,,2]
b.img <- img[,,3]
pheatmap::pheatmap(r.img)
pheatmap::pheatmap(g.img)
pheatmap::pheatmap(b.img)
##除去平均数
r.img <- r.img-mean(r.img)
g.img <- g.img-mean(g.img)
b.img <- b.img-mean(b.img)
##计算平面趋势。。。
r.img.trend <- fit3d(r.img)
g.img.trend <- fit3d(g.img)
b.img.trend <- fit3d(b.img)
##消除趋势
r.img.dtrend <- r.img-r.img.trend
g.img.dtrend <- g.img-g.img.trend
b.img.dtrend <- b.img-b.img.trend
##为LP高斯滤波器定义一个sigma
gaus.sig=15
##构建高斯滤波器
gaus <- build.gaus(nrow(img),ncol(img),gaus.sig)
##找到每个RGB通道的极值（绝对值最大值）
blob.r.point <- which(abs(r.img.dtrend)==max(abs(r.img.dtrend)),arr.ind=TRUE)
blob.g.point <- which(abs(g.img.dtrend)==max(abs(g.img.dtrend)),arr.ind=TRUE)
blob.b.point <- which(abs(b.img.dtrend)==max(abs(b.img.dtrend)),arr.ind=TRUE)
##设置要在连接组件算法中使用的窗口大小
win.size = 0.05
##提取blob xy位置
blob.r <- blob.extract(r.img.dtrend,blob.r.point,win.size,gaus)
blob.g <- blob.extract(g.img.dtrend,blob.r.point,win.size,gaus)
blob.b <- blob.extract(b.img.dtrend,blob.r.point,win.size,gaus)
# 注意：必须根据以下要求调整blob点（blob$xy.coords）
##其中原点（0,0）位于R图和图像图中
blob.coords.r <- blob.r$xy.coords
blob.coords.r[,1] <- blob.r$xy.coords[,2]
blob.coords.r[,2] <- (blob.r$xy.coords[,1]-nrow(r.img))*-1
blob.coords.g <- blob.g$xy.coords
blob.coords.g[,1] <- blob.g$xy.coords[,2]
blob.coords.g[,2] <- (blob.g$xy.coords[,1]-nrow(g.img))*-1
blob.coords.b <- blob.b$xy.coords
blob.coords.b[,1] <- blob.b$xy.coords[,2]
blob.coords.b[,2] <- (blob.b$xy.coords[,1]-nrow(b.img))*-1
##保存用户选项
mar.usr=par()$mar
dev.new()
close.screen(all.screen=TRUE)
par(mar=c(0,0,2,0))
split.screen(c(1,2))
split.screen(c(3,1),screen=2)
screen(1)
image2(sakurajima,asp=1,axes=FALSE)
rect(ybottom,nrow(sakurajima)-xleft,ytop,nrow(sakurajima)-xright,lwd=3,border='white',lty=3)
title('Original Image',line=0,font=2,col='white',cex=2,)
screen(3)
image2(r.img,asp=1,axes=FALSE)
points(blob.coords.r,col=rgb(1,0,0,alpha=0.05),pch=16,cex=0.3)
title('Red Channel',line=0,font=2,col='red',cex=2)
screen(4)
image2(g.img,asp=1,axes=FALSE)
points(blob.coords.g,col=rgb(0,1,0,alpha=0.05),pch=16,cex=0.3)
title('Green Channel',line=0,font=2,col='darkgreen',cex=2)
screen(5)
image2(b.img,asp=1,axes=FALSE)
points(blob.coords.b,col=rgb(0,0,1,alpha=0.05),pch=16,cex=0.3)
title('Blue Channel',line=0,font=2,col='darkblue',cex=2)
dev.copy(png, "3.png");dev.off()
##返回用户原始页边距并关闭屏幕
par(mar=mar.usr)
close.screen(all.screens=TRUE)
#end Sat Sep 11 01:01:22 2021 --
