fusion.Rmd

---
title: "Image Fusion"
output: 
  html_document:
    toc: true
    toc_float: true
    toc_collapsed: true
    toc_depth: 3
    number_sections: false
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Introduction

This tutorial seeks to illustrate how image fusion can be conducted. We use an Unmanned Aerial Vehicle (UAV) and Sentinel 2 optical image acquired within the same area and period. **Data fusion** is a formal framework in which are expressed means and tools for the alliance of data originating from different sources. It aims at obtaining information of greater quality; the exact definition of "greater quality” will depend upon the application ([Ranchin and Wald, 2010](https://doi.org/10.1007/978-1-4020-4385-7_11)).

The principal motivation for image fusion is to improve the quality of the information contained in the output image in a process known as synergy. A study of existing image fusion techniques and applications shows that image fusion can provide us with an output image with an improved quality. In this case, the benefits of image fusion include:

1. Extended range of operation.
2. Extended spatial and temporal coverage.
3. Reduced uncertainty.
4. Increased reliability.
5. Robust system performance.
6. Compact representation of information.

## Data preparation

Load libraries, declare variables and data paths.

```{r d1, message=FALSE}
rm(list=ls(all=TRUE))    #Clears R memory
unlink(".RData") 
if (!require("pacman")) install.packages("pacman"); library(pacman) #package manager( loads required packages/libraries in list as below if not installed they will be installed
p_load(raster, terra)

options(warn=1)
cat("Set variables and start processing\n")
Root 		<- 'D:/JKUAT/RESEARCH_Projects/Eswatini/Data/'
Path_out    <- paste0(Root,"Output/")

```

Load UAV and Sentinel 2 optical images. 

```{r d2}
path <- list.files(paste0(Root,'S2/interim/'),pattern = (".tif$"), recursive = TRUE, full.names = TRUE)
path
s <- rast(path)
s
path <- list.files(paste0(Root,'WingtraOne/'),pattern = (".tif$"), full.names = TRUE)
#reorder the bands to match those in S2
paths <- path
paths[3] <- path[4]
paths[4] <- path[3]
paths
v <- rast(paths)
v
```

Now check image properties and assign meaningful names to its bands.

```{r d3}
#Resolution
res(s)
#Extents
ext(s)
#image dimensions
dim(s)
#Number of bands
nlyr(s)
names(s) <- c("b", "g","r", "nir")
s
res(v)
ext(v)
dim(v)
names(v) <- c("b", "g","r", "nir")
v
```

Crop/clip Sentinel 2 image to UAV image extents.

```{r d4}
s <- crop(s, ext(v), snap="near")
```

Display the images side by side.

```{r d5, message = FALSE}
x11()
par(mfrow = c(1, 2)) #c(bottom, left, top, right)
plotRGB(s, r="nir", g="r", b="g", stretch="lin", axes=T, mar = c(4, 5, 1.4, 0.2), main="S2", cex.axis=0.5)
box()
plotRGB(v, r="nir", g="r", b="g", stretch="lin", axes=T, mar = c(4, 5, 1.4, 0.2), main="UAV", cex.axis=0.5)
box()

```

First let us conduct a spectral fusion of S2 and UAV images. To do this we have to resample UAV image to S2 extents and another to one to 1 m.

```{r d6}
v_r <- resample(v, s, method='bilinear')
ext(v_r) == ext(s)
temp <-rast(nrow=1570,ncol=760,ext(s)) #empty object to upscale UAV to 1 m resolution
res(temp) 
v1 <- resample(v, temp, method='bilinear')
res(v1) 
```

## Pixel fusion

This sections considers fusion techniques which rely on simple pixel based operations on input image values. The assumption is that the input images are spatially and temporally aligned, semantically equivalent and radiometrically calibrated. Therefore, let us fuse the two images by multiplication and display it against the original ones.

```{r d7}
#Fuse by multiplication
f1 <- s * v_r
f1
#Display fused image alongside original UAV
x11()
par(mfrow = c(1, 3),mar = c(4, 5, 1.4, 0.2)) 
plotRGB(v, r="nir", g="r", b="g", stretch="lin", axes=T, mar = c(4, 5, 1.4, 0.2), main="UAV", cex.axis=0.7)
box()
plotRGB(s, r="nir", g="r", b="g", stretch="lin", axes=T, mar = c(4, 5, 1.4, 0.2), main="S2", cex.axis=0.7)
box()
plotRGB(f1, r="nir", g="r", b="g", stretch="lin", axes=T, mar = c(4, 5, 1.4, 0.2), main="Fused_multi", cex.axis=0.7)
box()
```

What about mean fusion (i.e. taking the mean of each pixel's reflectance in both UAV and S2)?

```{r d8}
#Fuse by multiplication
f2 <- mean(s, v_r)
f2
#Display fused image alongside original UAV
x11()
par(mfrow = c(2, 2), mar = c(4, 5, 1.4, 0.2)) 
plotRGB(v, r="nir", g="r", b="g", stretch="lin",  main="UAV", axes=T, mar = c(4, 5, 1.4, 0.2))
box()
plotRGB(s, r="nir", g="r", b="g", stretch="lin", main="S2", axes=T, mar = c(4, 5, 1.4, 0.2))
box()
plotRGB(f1, r="nir", g="r", b="g", stretch="lin", main="Fused_mult", axes=T, mar = c(4, 5, 1.4, 0.2))
box()
plotRGB(f2, r="nir", g="r", b="g", stretch="lin", main="Fused_mean", axes=T, mar = c(4, 5, 1.4, 0.2))
box()

```

Let us finally follow the fusion approach in [Zou et al (2018)](https://ieeexplore.ieee.org/document/8812312).

```{r d9}
f3 = (s/v_r)*v_r
x11()
par(mfrow = c(2, 3), mar = c(4, 5, 1.4, 0.2)) 
plotRGB(v, r="nir", g="r", b="g", stretch="lin",  main="UAV", axes=T, mar = c(4, 5, 1.4, 0.2))
box()
plotRGB(s, r="nir", g="r", b="g", stretch="lin", main="S2", axes=T, mar = c(4, 5, 1.4, 0.2))
box()
plotRGB(f1, r="nir", g="r", b="g", stretch="lin", main="Fused_mult", axes=T, mar = c(4, 5, 1.4, 0.2))
box()
plotRGB(f2, r="nir", g="r", b="g", stretch="lin", main="Fused_mean", axes=T, mar = c(4, 5, 1.4, 0.2))
box()
plotRGB(f3, r="nir", g="r", b="g", stretch="lin", main="Zhou etal", axes=T, mar = c(4, 5, 1.4, 0.2))
box()

```

There seem to be some linear relationship between UAV and Sentinel 2 surface reflectance. However it is evident that reflectance values from UAV are higher compared to those in Sentinel 2. So what now?

## Feature based fusion

In *feature fusion* we fuse together the features $F_k,k \in{1,2, \dots, K}$. These features can be vegetation indices like Normalized Difference Index (NDVI) or feature maps that have been made semantically equivalent by transforming them into probabilistic $p(m,n)$, or likelihood, maps.

Let us start with NDVI ($\text{NDVI}=\frac{\text{NIR}-\text{Red}}{\text{NIR}+\text{Red}}$. First compute NDVI for both UAV and S2.

```{r n1}
n_v <- (subset(v_r,"nir")-subset(v_r,"r"))/(subset(v_r,"nir")+subset(v_r,"r"))
n_s <- (subset(s,"nir")-subset(s,"r"))/(subset(s,"nir")+subset(s,"r"))
```

How can we fuse the NDVI index? Let us take an average of the two.

```{r n2}
nf <- mean(n_v, n_s)
x11()
par(mfrow = c(2, 2), mar = c(5, 5, 1.4, 0.2))#c(bottom, left, top, right)
plot(n_v, main="UAV NDVI")
plot(n_s, main="S2 NDVI")
plot(nf, main="Fused NDVI")
```

Is there any difference between S2, UAV, and the fused NDVI images as shown above? 

## Spatial-spectral fusion

Previously we upsampled the UAV image in order to conduct fusion. While this reduces spectral variability it destroys spatial resolution. Therefore, in this section we will first donwsample the Satellite image to match UAV spatial resolution and then proceed to conduct image fusion. This way, we will improve both spatial and spectral information of Sentinel 2 image and spectral information for UAV.

### Modelling reflectance

1. Can we improve the resolution of S2 using UAV?
2. Can we predict UAV reflectance in places not imaged by the drones?

Let's check the relationship between the two.

```{r m1}
x11()
par(mfrow = c(2, 2), mar = c(4, 5, 1.4, 0.2)) 
plot(as.vector(subset(s,'b')),as.vector(subset(v_r,'b')), xlab='S2', ylab='W1 UAV', main="Blue band",pch=16,cex=0.75, col='blue')
plot(as.vector(subset(s,'g')),as.vector(subset(v_r,'g')), xlab='S2', ylab='W1 UAV', main="Green band",pch=16,cex=0.75, col='green')
plot(as.vector(subset(s,'r')),as.vector(subset(v_r,'r')), xlab='S2', ylab='W1 UAV', main="Red band",pch=16,cex=0.75, col='red')
plot(as.vector(subset(s,'nir')),as.vector(subset(v_r,'nir')), xlab='S2', ylab='W1 UAV', main="NIR band",pch=16,cex=0.75)
```

Lets create sample points from S2 (10 m resolution) and upscaled UAV (10 m resolution)  and use them to create a model that we can predict S2 reflectance based on UAV 1 m resolution. Essentially what we are doing here is to create a model that can predict S2 reflectance at a resolution of 1 m.

```{r m2}
set.seed(530)
x11()
points <- spatSample(v_r, 900, "random", as.points=T, na.rm=T, values=F)
plotRGB(v_r, r="nir", g="r", b="g", stretch="lin",  main="UAV+sampling points", axes=T, mar = c(4, 5, 1.4, 0.2))
plot(points,add=T)

```

Now extract reflectance values from the two images.

```{r m3}
s2_p <- extract(s, points, drop=F)
head(s2_p)
vr_p <- extract(v_r, points, drop=F)
head(vr_p)
```

Create a model.

```{r m4}
data <- data.frame(S2=s2_p[,-1], UAV=vr_p[,-1])
head(data)
# Plot the data
plot(S2.b~UAV.b, data=data, pch=16)
# Create a linear regression model
l.model <- lm(S2.b~UAV.b, data=data)
# Add the fitted line
abline(l.model, col="red")

```

Looks like the relation within the blue band is not linear. Let's try a non-linear SVM model.

```{r m5, message=FALSE}
#SVM
library(e1071)
svm.model <- svm(S2.b~UAV.b, data=data)
svm.pred <- predict(svm.model, data)
#SVM
library(randomForest)
rfmod <- randomForest(S2.b~UAV.b, data=data, ntree=300)
rf.pred <- predict(rfmod, data)
x11()
plot(S2.b~UAV.b,data, pch=16)
points(data$S2.b, svm.pred, col = "blue", pch=4)
points(data$S2.b, rf.pred, col = "red", pch=4)
legend("topright",c("Data","SVM","RF"), pch= c(16, 4, 4),col=c("black", "blue","red"))

```

### Predicting Sentinel 2 reflectance

SVM and RF models have better characterized the relationship between S2 and UAV blue bands. Let predict high resolution S2 band from existing UAV.

```{r m6}
UAV <- v1[['b']] 
names(UAV) <-'UAV.b'
s2h.b <- predict(UAV, rfmod, na.rm=T)

s2h.svm <- predict(UAV, svm.model, na.rm=T)

x11()
par(mfrow = c(1, 2), mar = c(4, 5, 1.4, 0.2)) 
plot(s2h.svm,  main="svm S2 1m Blue band")
plot(s2h.b, main="RF S2 1m Blue band")

# #parallel processing in raster
# UAv <- raster(UAV)
# library(snow)
# startTime <- Sys.time() 
# cat("Start time", format(startTime),"\n")
# beginCluster()
#  r4 <- predict(UAV, rfmod,na.rm=T)
# endCluster()
# timeDiff <- Sys.time() - startTime
# cat("\n Processing time", format(timeDiff), "\n")

```

### Validation

But how do we know if this predictions are accurate? We can plot predicted vs actual and also compute RMSE and Mean Absolute Percentage Error (MAPE) using training data. MAPE is given as:

```{r mape}
MAPE <- function (y_pred, y_true){
    MAPE <- mean(abs((y_true - y_pred)/y_true))
    return(MAPE*100)
}
```

and RMSE,

```{r rmse}
rmse <- function(error){
  sqrt(mean(error^2))
}
```

So lets compute MAPE and RMSE for both methods.

```{r m7}
svm.rmse <- rmse(svm.pred-data$S2.b)
svm.rmse
rf.rmse <- rmse(rf.pred-data$S2.b)
rf.rmse
svm.mape <- MAPE(svm.pred, data$S2.b)
svm.mape
rf.mape <- MAPE(rf.pred, data$S2.b)
rf.mape
```

From the validations Random Forest gives better prediction than Support Vector machines because it has low MAPE and RMSE error. In that case we can adopt the high resolution S2 image predicted/simulated by RF and fuse it with UAV image at 1 m resolution. Note that we only predict the blue band, we can model and predict the other bands and then fuse them.

```{r m8}
#Green band
rfmod <- randomForest(S2.g~UAV.g, data=data)
UAV=subset(v1,'g')
names(UAV) <-'UAV.g'
s2h.g <- predict(UAV, rfmod, na.rm=T)
#red band
rfmod <- randomForest(S2.r~UAV.r, data=data)
UAV=subset(v1,'r')
names(UAV) <-'UAV.r'
s2h.r <- predict(UAV, rfmod, na.rm=T)
#NIR band
rfmod <- randomForest(S2.nir~UAV.nir, data=data)
UAV=subset(v1,'nir')
names(UAV) <-'UAV.nir'
s2h.nir <- predict(UAV, rfmod, na.rm=T)

#Stack them
s2h <- stack(x=c(s2h.b,s2h.g,s2h.r,s2h.nir))
names(s2h) <- c("b", "g","r", "nir")
s2h
```

However, the predicted reflectance ranges in all bands except NIR are very low. Why could this be so?  Now lets create a high resolution fused image using Zhou's approach.

### Fusion

```{r m9}
s2h.fused <- (rast(s2h)/v1)*v1
s2h.fused

x11()
par(mfrow = c(1, 2), mar = c(4, 5, 1.4, 0.2)) 
plotRGB(s2h.fused, r="nir", g="r", b="g", stretch="lin",  main="UAV+S2 high res", axes=T, mar = c(4, 5, 1.4, 0.2))
box()
plotRGB(v, r="nir", g="r", b="g", stretch="lin", main="UAV", axes=T, mar = c(4, 5, 1.4, 0.2))
box()
```

The other option would be use S2 to predict UAV in areas not covered by the drone and then fuse it with S2. However this approach would require consideration of similar land-cover. For instance, we can not train the model in an area with different land-cover say cropland and predict in another area with say Forest. Food for thought.

But let's step back a little, could sampling from land-cover categories have improved prediction of high resolution S2 done previously with RF? To test this, let us perform a K-means classification and sample from its land-cover map. 

```{r m10}
image <- v_r
image[is.na(image)] <- 0
nclass <- 4
system.time(
E <- kmeans(as.data.frame(image, na.rm=F), nclass, iter.max = 100, nstart = 9)
)
k.map <- image[[1]]
values(k.map) <- E$cluster
k.map[is.na(v_r[[1]])] <- NA
plot(k.map)
```

Sample points using stratified random sampling based on the K-means land-cover map.

```{r m11}
set.seed(530)
points <- spatSample(k.map, 900, "stratified", as.points=T, na.rm=T, values=F)
x11()
plot(k.map,  main="K-means +sampling points", axes=T, mar = c(4, 5, 1.4, 0.2), col=topo.colors(max(values(k.map),na.rm=T)))
plot(points,add=T)
```

Use the stratified randomly sample points to train and predict simulated high resolution S2.

```{r m12}
s2_p <- extract(s, points, drop=F)
head(s2_p)
vr_p <- extract(v_r, points, drop=F)
head(vr_p)
data <- data.frame(S2=s2_p[,-1], UAV=vr_p[,-1])
rfmod <- randomForest(S2.b~UAV.b, data=data)
rf.pred <- predict(rfmod, data)
```

Lets evaluate if stratified random sampling improved RF accuracy.

```{r m13}
rf.rmse <- rmse(rf.pred-data$S2.b)
rf.rmse
rf.mape <- MAPE(rf.pred, data$S2.b)
rf.mape
```

We can see that RMSE has by a small margin compared to the case of random sampling. MAPE more or less remained constant. Thus there is a chance that sampling over different land-cover improves prediction. 

As an assignment now predict S2 high resolution using a RF model trained on stratified random samples from K-means classifier.

## References

Thierry Ranchin and Lucien Wald. *Data Fusion in Remote Sensing of Urban and Suburban
Areas*, pages 193–218. Springer Netherlands, Dordrecht, 2010. ISBN 978-1-4020-4385-
7. doi: 10:1007/978-1-4020-4385-7 11. URL https://doi.org/10.1007/978-1-4020-4385-7_11.

Y. Zou, G. Li and S. Wang, "The Fusion of Satellite and Unmanned Aerial Vehicle (UAV) Imagery for Improving Classification Performance," *IEEE International Conference on Information and Automation (ICIA)*, 2018, pp. 836-841, doi: 10.1109/ICInfA.2018.8812312.