-
Notifications
You must be signed in to change notification settings - Fork 0
/
project_functions.R
73 lines (69 loc) · 1.73 KB
/
project_functions.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
# Fast Walsh-Hadamard Transform
fwht2d <- function(x) {
h <- 1
length <- ncol(x)
while (h < length) {
for (i in seq(1, length, by=h*2)) {
for (j in seq(i,i+h-1)) {
a <- x[,j]
b <- x[,j+h]
x[,j] <- a + b
x[,j+h] <- a - b
}
}
h <- 2*h
}
h <- 1
length <- nrow(x)
while (h < length) {
for (i in seq(1, length, by=h*2)) {
for (j in seq(i,i+h-1)) {
a <- x[j, ]
b <- x[j+h, ]
x[j, ] <- a + b
x[j+h, ] <- a - b
}
}
h <- 2*h
}
x
}
hard.thresh <- function(x, lambda) {ifelse(abs(x) > lambda, x, 0)}
soft.thresh <- function(x, lambda) {
sign(x) * (abs(x) >= lambda) * (abs(x) - lambda)
}
sub.im.transform <- function(image, size, Th.type="soft", thresh) {
# image size must be square
N <- ncol(image)
if ((N %% size) != 0) {
stop("image size is not divisible by specified sub-image size")
}
# create a list of submatrices to apply the transform to
flat <- as.vector(image)
subimages <- list()
A <- N^2/size^2 # index span for each subimage
for (j in 0:(A-1)) {
subimages[[j+1]] <- matrix(flat[(j*(size^2)+1):((j+1)*(size^2))],ncol=size)
}
# apply transform to each sub-image and reconstruct the original matrix/image
subimages.hat <- list()
i <- 1
for (im in subimages) {
Zhat <- fwht2d(im)
if (Th.type == "soft") {
Zhat.star <- soft.thresh(Zhat, thresh)
} else {
Zhat.star <- hard.thresh(Zhat, thresh)
}
Zstar <- fwht2d(Zhat.star) / size
subimages.hat[[i]] <- Zstar
i <- i + 1
}
# reconstruct original image
original <- c()
for (i in 1:A) {
original <- c(original, as.vector(subimages.hat[[i]]))
}
original <- matrix(original, ncol = 256)
original
}