|
##引入包
library(scatterplot3d)
library(RSpectra)
library(jpeg)
library(animation)
library(rgl)
##一键清空
rm(list=ls())
##样本个数
num<-100
##维度
v<-3
##降维后维度
k<-2
##设置x,y,z轴的旋转角度
angle1<-20/180*pi
angle2<-15/180*pi
angle3<-10/180*pi
##生成符合高斯分布的数据
x1<-rnorm(100,mean=1,sd=1.8)
x2<-rnorm(100,mean=2,sd=1.9)
x3<-rnorm(100,mean=4,sd=1.7)
##数据旋转
Rx<-matrix(c(1,0,0,
? ? ? ? ? ? ?0,cos(angle1),-sin(angle1),
? ? ? ? ? ? ?0,sin(angle1),cos(angle1)),nrow=3,ncol=3,byrow=TRUE)
Ry<-matrix(c(cos(angle2),0,sin(angle2),
? ? ? ? ? ? ?0,1,0,
? ? ? ? ? ? ?-sin(angle2),0,cos(angle2)),nrow=3,ncol=3,byrow=TRUE)
Rz<-matrix(c(cos(angle3),-sin(angle3),0,
? ? ? ? ? ? ?sin(angle3),cos(angle3),0,
? ? ? ? ? ? ?0,0,1),nrow=3,ncol=3,byrow=TRUE)
X<-matrix(c(x1,x2,x3),nrow=v,ncol=num,byrow=TRUE)
X<-Rz%*%Ry%*%Rx%*%X
XX<-X
##显示原始数据的分布
scatterplot3d(X[1,],X[2,],X[3,],main='1.4.1')
plot3d(XX[1,],XX[2,],XX[3,],col=rainbow(12),size=5)
##对样本正则化
for(i in 1:v){
? total<-0
? for(j in 1:num){
? ? total<-total+X[i,j]
? }
? total<-total/num
? X[i,]<-X[i,]-total
}
##特征值分解
C<-X%*%t(X)
ev<-svd(X)
vec<-ev$v[,order(ev$d,decreasing = TRUE)]
val<-ev$d[order(ev$d,decreasing = TRUE)]
##取前k个特征向量
P<-t(vec[,1:k])
##计算低维坐标系的投影
Y<-P%*%t(X)%*%XX
plot(Y[1,],Y[2,],main='1.4.3',
? ? ?xlab='x1', ylab='x2',col="red")
##重构样本集
YY<-matrix(0,nrow=v,ncol=num)
for(i in 1:num){
? hy<-rep(0,v)
? for(j in 1:k){
? ? hy<-hy+Y[j,i]*vec[j,]
? }
? YY[,i]<-hy
}
scatterplot3d(YY[1,],YY[2,],YY[3,],main='1.4.1')
plot3d(YY[1,],YY[2,],YY[3,],col=rainbow(12),size=5)
##载入图片，并且显示出来
library(raster)
library(jpeg)
##图片数量
num<-4
##降维后图片维度
k<-50
##图片维度
v<-500*500
photo.raster<-array(0,dim=c(500,500,num))
raster.photo <- raster("sam_1.jpg")
photo.flip <- flip(raster.photo, direction = "y")
##将数据转换为矩阵形式
?photo.raster[,,1] <- t(as.matrix(photo.flip))
?##灰化处理
?image(photo.raster[,,1], col = grey(seq(0, 1, length = 256))) ?
?raster.photo<- raster("sam_2.jpg")
?photo.flip <- flip(raster.photo, direction = "y")
?##将数据转换为矩阵形式
?photo.raster[,,2] <- t(as.matrix(photo.flip))
?##灰化处理
?image(photo.raster[,,2], col = grey(seq(0, 1, length = 256))) ?
?raster.photo<- raster("sam_3.jpg")
?photo.flip <- flip(raster.photo, direction = "y")
?##将数据转换为矩阵形式
?photo.raster[,,3]<- t(as.matrix(photo.flip))
?##灰化处理
?image(photo.raster[,,3], col = grey(seq(0, 1, length = 256)))?
?raster.photo<- raster("sam_4.jpg")
?photo.flip <- flip(raster.photo, direction = "y")
?##将数据转换为矩阵形式
?photo.raster[,,4] <- t(as.matrix(photo.flip))
?##灰化处理
?image(photo.raster[,,4], col = grey(seq(0, 1, length = 256)))
?##计算平均脸
?pic_mean<-matrix(0,nrow=500,ncol=500)
?for(i in 1:num){
? ?pic_mean<-pic_mean+photo.raster[,,i]
?}?
?pic_mean<- pic_mean/num
?for(i in 1:num){
? ?photo.raster[,,i]<-photo.raster[,,i]-pic_mean
?}
?for(i in 1:num){
? ?##奇异值进行分解
?photo.svd <- svd(photo.raster[,,i])
?d <- diag(photo.svd$d)
?v <- as.matrix(photo.svd$v)
?u <- photo.svd$u
?##取前k个特征向量
u1 <- as.matrix(u[, 1:k])
?d1 <- as.matrix(d[1:k, 1:k])
?v1 <- as.matrix(v[, 1:k])
?##重构样本集
?photo1 <- u1 %*% d1 %*% t(v1)
?##将降维后的图片显示出来
?image(photo1+pic_mean, col = grey(seq(0, 1, length = 256)))
?cat(10*log(norm(photo1)/norm(photo1-photo.raster[,,i]))/log(10))
?}
?
|