################################################################################
# E05: R語言資料視覺化與影像資料分析                                                #
# 吳漢銘 國立政治大學統計學系                   #
# https://hmwu.idv.tw/                          #
#################################################



################################################################################
# image                                                                        #  
################################################################################
#install.packages("imager")
library(imager)

########################################
# read image                           #  
########################################
# imager supports JPEG, PNG, TIFF and BMP natively - 
# for other formats you'll need to install ImageMagick.

# Examples
# imager build-in: boats, hubble, parrots
# imager_extdata: coins.png, HubbleDeepField.jpg, Leonardo_Birds.jpg, parrots.png
# image_other: Cars-on-Busy-Street.jpg, Car_movie.jpg, ironman.jpg, plaza.jpg
# image_medical: Brain-fMRI-n7.bmp, Breast1_US.bmp, Breast1_US_ROI.bmp
#                Breast2_US.bmp, Breast2_US_doctor.bmp
#                Breast_US.jpg, Breast_US_doctor.jpg
#                Heart_MRI.bmp, Heart_MRI_ROI.bmp
#                Knee_MRI.bmp, Knee_MRI_ROI.bmp
#                Liver.bmp, Liver_ROI.bmp
#                simu-c2-sd10.bmp, simu-c2-sd20.bmp, simu-c5-sd20.bmp, simu-ring.bmp
#                texture_1.bmp, texture_2.bmp, texture_3.bmp, texture_4.bmp

#file <- system.file('extdata/parrots.png', package='imager')
#file <- "image_other/Cars-on-Busy-Street.jpg"
#file <- "image_medical/Breast_US.jpg"
#file <- "image_medical/Knee_MRI.bmp"

file <- "image_other/ironman.jpg"
file

# im <- boats
# im <- load.example("hubble")
im <- load.image(file)
im
dim(im)
str(im)
class(im)
plot(im)
text(100, 50, label="ironman", col="white", cex=2)

#  
width(im)
height(im)
#depth(im)
spectrum(im)

# usual arithmetic operations
# 0-1 (dark to light)
range(im)
mean(im)
sd(im)
log(im)+3*sqrt(im)
im[1:10, 1:5, 1, 1]

# save image (base R has a save.image function)
imager::save.image(im, "image-new.jpg")


########################################
# The cimg class                       #  
########################################
# cimg class: a regular 4d array with an S3 class
# construct a cimg image
nc <- 10
nr <- 10
set.seed(12345)
#noise <- matrix(runif(nc*nr, min = 0.2, max = 0.8), ncol=nc, nrow=nr)
noise <- matrix(sort(runif(nc*nr, min = 0.2, max = 0.8)), ncol=nc, nrow=nr)
noise[1, 1] <- 0
noise[1, 10] <- 0
noise[10, 10] <- 1
noise[1:5, 1:5]
im.noise <- as.cimg(noise)
par(mfrow=c(1, 2))
image(t(noise)[, 10:1], main="image (matrix)", col=gray(0:100/100))
plot(im.noise, main="cimg (matrix)")

head(as.data.frame(im.noise))


########################################
# gray scale                           #  
########################################
par(mfrow=c(1, 2))
plot(grayscale(im, method = "Luma"), main="grayscale (Luma)") 
plot(grayscale(im, method = "XYZ"), main="grayscale (XYZ)") 


########################################
# colour channel                       #  
########################################
R(im)  # same as channel(im, 1)
G(im)
B(im)

im.temp <- im
G(im.temp) <- 0
B(im.temp) <- 0
plot(im.temp)

im.gray <- grayscale(im)
par(mfrow=c(2, 2))
hist(im.gray, main="Luminance values")
hist(R(im), main="Red channel values")
hist(G(im), main="Green channel values")
hist(B(im), main="Blue channel values")


########################################
# Sub-images                           #  
########################################
imsub(im, x < 30) #Only the first 30 columns
imsub(im, y < 30) #Only the first 30 rows
imsub(im, x < 30, y < 30) #First 30 columns and rows
imsub(im, sqrt(x) > 8) #Can use arbitrary expressions
imsub(im, x > height/2, y > width/2) 
imsub(im, cc==1) # Colour axis is "cc"
par(mfrow=c(2, 1))
plot(im, main="image")
im.ROI <- imsub(im, 200 < x & x <= 350, 0 < y & y <= 200)
plot(im.ROI, main="sub-image (ROI)")


# Rows, columns
plot(imrow(R(im), 10), main="Red channel along 10th row", type="l")

# individual pixels
at(im, x=20, y=20, cc=1:3)
color.at(im, x=20, y=20)


########################################
# Blurry                               #  
########################################
par(mfrow=c(2, 3))
for(i in seq(0, 10, 2)){
  im.blurry <- isoblur(im, sigma=i) 
  plot(im.blurry, main=paste("Blurry: sigma=", i))
}


########################################
# Denoising                            #  
########################################
# filters that average over space
par(mfrow=c(1, 3))
plot(im, main="Original")
plot(isoblur(im, 5), main="Blurred")
im.blur.ani <- blur_anisotropic(im, ampl=1e3, sharp=0.3)
plot(im.blur.ani, main="Blurred (anisotropic)")


########################################
# Resizing                             #  
########################################
par(mfrow=c(1, 2))
plot(im, main="Original")
im.thumb <- resize(im, round(width(im)/4),round(height(im)/4))
plot(im.thumb, main="1/4 size") 

#Same as above: negative arguments are interpreted as percentages
im.thumb <- resize(im, -25, -25)


########################################
# Rotation                             #  
########################################
par(mfrow=c(1, 2))
plot(imrotate(im, 30), main="Rotating (30 degree)")
plot(imrotate(im, -30), main="Rotating (-30 degree)")

########################################
# Shifting                             #  
########################################
im.shift1 <- imshift(im, 40, 20) 
im.shift2 <- imshift(im, 100, 100, boundary=1)
im.shift3 <- imshift(im, 100, 100, boundary=2)
par(mfrow=c(1, 3))
plot(im.shift1, main="Shifting")
plot(im.shift2, main="Shifting (Neumann boundaries)")
plot(im.shift3, main="Shifting (circular)")


########################################
# Edge detection: Deriche filter       #  
########################################
# Gaussian (and derivative-of-Gaussian) filters
# The Deriche filter: a fast approximation to a Gaussian filter (order = 0), 
#                     or Gaussian derivatives (order = 1 or 2).
im.df <- deriche(im, sigma=4, order=2, axis="y") 
# Edge detector along x-axis: apply recursive Deriche filter
par(mfrow=c(2, 1))
im.edges.x <- deriche(im, sigma=2, order=2, axis="x") 
plot(im.edges.x, main="Edge detector along x-axis")
im.edges.y <- deriche(im, sigma=2, order=2, axis="y") 
plot(im.edges.y, main="Edge detector along y-axis")

# Chain operations
par(mfrow=c(1, 2))
im.edges.xy <- deriche(im.edges.x, 2, order=2, axis="y")
im.edges.yx <- deriche(im.edges.y, 2, order=2, axis="x")
plot(im.edges.xy, main="Edge detector along xy-axis")
plot(im.edges.yx, main="Edge detector along yx-axis")


# The Young-van Vliet filter:a fast approximation to a Gaussian filter (order = 0),
#                            or Gaussian derivatives (order = 1 or 2).
im.vf <- vanvliet(im, sigma=4, order=2, axis="y")
plot(im.vf, main="2nd deriv of Gaussian along y")

########################################
# Edge detection: gradient             #  
########################################
# image gradient along x and y axes
par(mfrow=c(1, 4))
im.gray <- grayscale(im)
im.dx <- imgradient(im.gray, "x")
im.dy <- imgradient(im.gray, "y")
# same as im.dxy <- imgradient(im.gray, "xy")
plot(im.gray, main="image (gray)")
plot(im.dx, main="gradient along x-axis")
plot(im.dy, main="gradient along y-axis")

im.grad.mag <- sqrt(im.dx^2+im.dy^2)
plot(im.grad.mag, main="Gradient magnitude")


########################################
# Box Filter,  Edge filter             #  
########################################
# 8x8 Box filter
box.8x8 <- as.cimg(matrix(1, 8, 8)) 
box.8x8[]
im.bf <- correlate(im.gray, box.8x8) 
par(mfrow=c(1, 2))
plot(im.gray, main="Original (gray)")
plot(im.bf, main="Filtering with box filter")

# Edge filter
edge.filter <- as.cimg(function(x,y) sign(x-5), 9, 9) 
edge.filter[]
im.corr <- correlate(im.gray, edge.filter) 
im.conv <- convolve(im.gray, edge.filter) 
par(mfrow=c(1, 2))
plot(im.corr, main="Correlation")
plot(im.conv, main="Convolution")

########################################
# FFTs and the periodic decomposition  #  
########################################
im.fft <- fft(as.matrix(im.gray))
par(mfrow=c(1, 3))
plot(as.cimg(Re(im.fft)), main="Real part of the transform")
plot(as.cimg(Im(im.fft)), main="Imaginary part of the transform")
im.fft.ps <- as.cimg(sqrt(Re(im.fft)^2+Im(im.fft)^2))
plot(im.fft.ps, main="Power spectrum")
par(mfrow=c(1, 2))
plot(sort(Re(im.fft)), main="Re(fft)")
plot(sort(Im(im.fft)), main="Im(fft)")


########################################
# Pixel sets (pixsets)                 #  
########################################
# https://cran.r-project.org/web/packages/imager/vignettes/pixsets.html
# A pixset is a set of pixels, represented as a binary image.
# Pixsets represent sets of pixels, or equivalently binary images (masks)

im.gray.px <- im.gray > 0.6 #Select pixels with high luminance
str(im.gray.px)
# The “TRUE”(white) values correspond to pixels in the set, and “FALSE”(black) to pixels not in the set. 

par(mfrow=c(2, 1))
plot(im.gray, main="image (gray)")
plot(im.gray.px, main="image (gray>0.6)")

# Coordinates for pixels in pixsets
pixsets <- where(im.gray.px)
head(pixsets)
sum(im.gray.px) #Number of pixels in set
mean(im.gray.px) #Proportion

# Indexing using pixsets
mean(im.gray[im.gray.px])

# Return coordinates of subset of pixels
im.cond <- get.locations(im, condition=function(x) x < 0.5)
head(im.cond)

im.tmp <- as.cimg(im.gray.px)
im.tmp
plot(im.tmp, main="as cimg image with 0/1") # as cimg image with 0/1


########################################
# Segmentation using pixsets (EX1)     #  
########################################
par(mfrow=c(1, 2))
im.blur <- isoblur(im.gray, 4)
plot(im.blur, main="image (blur)")
im.blur.px <- im.blur  > 0.5
plot(im.blur.px, main="pixset (T/F)")

par(mfrow=c(1, 2))
plot(im.gray, main="segmentation (highlight)")
highlight(im.blur.px) #  col="lightgreen"

im.seg.color <- colorise(im.blur, im.blur.px, "red", alpha=0.3)
plot(im.seg.color, main="segmentation (color)")

par(mfrow=c(1, 2))
im.seg.bound <- boundary(im.blur.px) # Boundaries
plot(im.seg.bound, main="pixset (boundary: T/F)") 

plot(im.gray, main="segmentation (boundary)")
points(where(im.seg.bound), cex=0.1, col="green")


# Working with colour images
# each colour channel is treated as having its 
# own set of pixels, so that a colour pixset has the same 
# dimension as the image it originated from.


########################################
# Segmentation using pixsets (EX2)     #  
########################################
# https://cran.r-project.org/web/packages/imager/vignettes/pixsets.html
im.coins <- load.example("coins")
par(mfrow=c(1, 3))
plot(im.coins, main="image")
str(im.coins)
hist(im.coins, main="gray levels")
lambda <- 0.4
abline(v=lambda, col="red")
plot(threshold(im.coins, thr=lambda), main="threshold")

# correct the illumination using a linear model
# by removing the trend
par(mfrow=c(2, 3))
im.coins.df <- as.data.frame(im.coins)
str(im.coins.df)
head(im.coins.df)
plot(as.cimg(im.coins.df, dims=dim(im.coins)), main="(1) image (df)")

im.coins.df.sub <- im.coins.df[sample(nrow(im.coins.df), 10000), ]
mylm <- lm(value ~ x * y, data=im.coins.df.sub)     
im.coins.lm <- im.coins.df
im.coins.lm$value <- im.coins.df$value - predict(mylm, im.coins.df)
im.coins.lm.cimg <- as.cimg(im.coins.lm, dims=dim(im.coins))
plot(im.coins.lm.cimg, main="(2) removing the trend by lm")
im.coins.lm.cimg.thr <- threshold(im.coins.lm.cimg, thr="auto")
plot(im.coins.lm.cimg.thr, main="(3) threshold after lm")

im.coins.lm.cimg.thr.clean <- clean(im.coins.lm.cimg.thr, 3)
im.coins.lm.cimg.thr.clean.fill <- fill(im.coins.lm.cimg.thr.clean, 7)
plot(im.coins.lm.cimg.thr.clean, main="(4) clean")
plot(im.coins.lm.cimg.thr.clean.fill, main="(5) fill")
plot(im.coins, main="(6) segmentation result")
highlight(im.coins.lm.cimg.thr.clean.fill)


########################################
# Segmentation (using clustering)      #  
########################################
# im <- load.image("image_other/ironman.jpg")
# im.gray <- grayscale(im)
# extract_patches: a list of image patches (cimg objects)
dim(im.gray)
im.nrow <- width(im.gray)
im.ncol <- height(im.gray)

block.length <- 7
st.x <- (block.length+1)/2
en.x <- width(im)-(block.length-1)/2
st.y <- (block.length+1)/2
en.y <- height(im)-(block.length-1)/2
cx <- rep(st.x:en.x, height(im)-block.length+1)
cy <- rep(st.y:en.y, each=width(im)-block.length+1)
out <- extract_patches(im.gray, cx, cy, wx=block.length, wy=block.length)
out
image.X.4x4 <- data.frame(matrix(unlist(out), nrow=length(out), byrow=T))
dim(image.X.4x4)
head(image.X.4x4)

seg.nrow <- im.nrow-block.length+1 
seg.ncol <- im.ncol-block.length+1

# K-means
no.class <- 4
kmeans.seg <- kmeans(image.X.4x4, no.class)$cluster

# K-means: reconstruction
im.kmeans.seg <- as.cimg(matrix(kmeans.seg, nrow=seg.nrow, ncol=seg.ncol))
dim(im.kmeans.seg)
par(mfrow=c(1, 3))
plot(im.gray, main="image")
plot(im.kmeans.seg, main="kmeans (cluster)")
im.gray.sub <- imsub(im.gray, st.x <= x & x <= en.x, st.y <= y & y <= en.y)
dim(im.gray.sub)
plot(im.gray.sub, main="segmentation (kmeans)")
highlight(im.kmeans.seg)


## FCM
par(mfrow=c(1, 3))
library(e1071)
cmeans.seg <- cmeans(image.X.4x4, no.class)$cluster
im.cmeans.seg <- as.cimg(matrix(cmeans.seg, nrow=seg.nrow, ncol=seg.ncol))
par(mfrow=c(1, 3))
plot(im.cmeans.seg, main="FCM (cluster)")
plot(im.gray.sub, main="segmentation (FCM)")
highlight(im.cmeans.seg)


## PCA+FCM
pca.cmeans.seg <- cmeans(princomp(image.X.4x4, 3)$scores, no.class)$cluster
im.pca.cmeans.seg <- as.cimg(matrix(pca.cmeans.seg, nrow=seg.nrow, ncol=seg.ncol))
plot(im.gray.sub, main="segmentation (PCA+FCM)")
highlight(im.pca.cmeans.seg)


########################################
# Foreground/background segmentation   #
########################################
# http://dahtah.github.io/imager/foreground_background.html


# https://upload.wikimedia.org/wikipedia/commons/thumb/f/fd/Aster_Tataricus.JPG/1024px-Aster_Tataricus.JPG
im <- load.image("image_other/1024px-Aster_Tataricus.JPG")
plot(im)


##X is training data 
##Xp is test data
##cl: labels for the rows of X
##k = number of neighbours
##Returns average value of k nearest neighbours
fknn <- function(X, Xp, cl, k=1){
  out <- nabor::knn(X, Xp, k=k)
  cl[as.vector(out$nn.idx)] %>% matrix(dim(out$nn.idx)) %>% rowMeans
}
























################################################################################
# video                                                                        #  
################################################################################
# install ffmpeg for videos
# https://ffmpeg.org/download.html
# copy ffmpeg-20200612-38737b3-win64-static\bin\ffmpeg.exe
#  to C:\Windows\System32
vifile <- system.file('extdata/tennis_sif.mp4', package='imager')
vi <- load.video(vifile)
dim(vi)
str(vi)
class(vi)
play(vi)
play(vi, loop = TRUE)

par(mfrow=c(2, 4))
lapply(1:8, function(i) plot(vi, frame=i, main=paste("frame", i)))

# convert the video to grayscale, 
#Differentiate along z axis and square
vi.gray <- grayscale(vi)
play(vi.gray, loop = TRUE)

# motion detector
vi.motion <- deriche(vi.gray, sigma=1, order=1, axis="z")^2 
str(vi.motion)

# combine both videos to display them side-by-side
# Paste the two videos together
vi.combined <- imappend(list(vi.motion/max(vi.motion), vi.gray/max(vi.gray)), 
                        axis="x") 
play(vi.combined, loop = TRUE)


# Skip to time = 2s, and grab the next 3 frames
vi.sub <- load.video(vifile, skip.to=2, frames=4)
# Split an image along a certain axis (producing a list)
vi.sub.split <- imsplit(vi.sub, "z")
# Display individual frames
plot(vi.sub.split)

# load video 1 frame per sec
vi.fps <- load.video(vifile, fps=1) 
plot(imsplit(vi.fps, "z"))

# save video
save.video(vi.sub, fname="my_sub_video.avi")

## #Making a video from a directory
dd <- tempdir()
for (i in 1:length(iml)) {
png(sprintf("%s/image-%i.png",dd,i));
plot(iml[[i]]); dev.off() }
make.video(dd,f)
play(load.video(f))


frame(vi, 1)
frame(vi, 1) <- isoblur(frame(vi, 1), 10)

# Lagged operators

par(mfrow=c(1, 2))
#Compute difference between two successive frames (at lag 1)
im.diff.lag1 <- imshift(vi, delta_z=1) - vi

plot(im.diff.lag1, frame=2, main="Difference betw. frames 2 and 1")

#Compute difference between frames (at lag 3)
im.diff.lag3 <- imshift(vi, delta_z=3) - vi
plot(im.diff.lag3, frame=4, main="Difference between frames 3 and 1")



########################################
# The cimg class                       #  
########################################
# cimg class: a regular 4d array with an S3 class

# construct a cimg video with 5x5 pixels, 5 frames, 3 colours. 
nc <- 256
nr <- 128
nf <- 50
vi.noise <- array(rnorm(nc*nr*nf*3), c(nc, nr, nf, 3)) 
vi.noise <- as.cimg(vi.noise)
play(vi.noise, loop = TRUE)
 
# as.data.frame, as.array, as.vector, as.matrix 
head(as.data.frame(vi.noise))