---
title: "Vignette for imagine package"
author: "Wencheng Lau-Medrano"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{imagine}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, echo = FALSE, message = FALSE}
library(imagine)
```

# imagine

**Imaging Engine, Tools for Application of Image Filters to Data Matrices**

This package was built to take numeric data matrices and apply image-filtering algoriths.

Algoriths used by `imagine` include median-filter and 2D-convolution based algoriths implemented on Rcpp (C++) because sometimes the applicaion of this filters to large numeric matrices is very intensive and time expensive.


## Installation

For installing `imagine`, as follows:

```{r, eval=FALSE}
install.packages("imagine")
```

## Engines
`imagine` performs algorithms written in Rcpp (called 'engines'), ensuring a faster application of filters. At version 1.0.0, `imagine` include five main engines, described as follows:

* Engine 1: Basic convolution algorithm. Multiplies kernel by the neighborhood of each cell. Then, it sums all values.
* Engine 2: Multiplies kernel by the neighborhood of each cell. Then, it returns the average of all values.
* Engine 3: Multiplies kernel by the neighborhood of each cell. Then, it returns the position indicated by the parameter `x`.
* Engine 4: It uses the `radius` argument to extract the values of squared neighborhood ($radius \times radius$) and calculates the mean.
* Engine 5: It uses the `radius` argument to extract the values of squared neighborhood ($radius \times radius$) and returns the position indicated by the parameter `x`

## Main functions

There are 5 main functions and 2 wrappers at version 1.1.0, described as follows:

### Convolution functions
```{r, eval=FALSE}
# Build kernels
# Kernel 1: For bottom edge recognition
kernel1 <- matrix(c(-1, -2, -1,
                     0,  0,  0,
                     1,  2,  1), nrow = 3)

# Kernel 2: Diagonal weighting
kernel2 <- matrix(c(1, 0, 1,
                    0, 2, 0,
                    1, 0, 1), nrow = 3)

# Apply filters
convolutionExample  <- convolution2D(dataMatrix = wbImage, kernel = kernel1)
convMeanExample     <- convolutionMean(dataMatrix = wbImage, kernel = kernel2)
convQuantileExample <- convolutionQuantile(dataMatrix = wbImage, kernel = kernel2, x = 0.7)

```

For compare results, we will plot both data (original and filtered) using `image` function, as shows in figures 1 and 2.

**Original**

```{r, message=FALSE, fig.height=3, fig.width=5.33, fig.cap = "Figure 2: Original matrix", results='hide', fig.pos="h", echo=FALSE}
par(mar = rep(0, 4), mfrow = c(1, 1))
cols <- colorRampPalette(colors = c("black", "white"))(1e4)

image(wbImage, col = cols)
```


**Filtered**

```{r, message=FALSE, fig.height=9, fig.width=5.33, fig.cap = "Figure 1: Filtered matrices", results='hide', fig.pos="h", echo=FALSE}

myMatrix <- wbImage

# Build kernels
# Kernel 1: For bottom edge recognition
kernel1 <- matrix(c(-1, -2, -1,
                     0,  0,  0,
                     1,  2,  1), nrow = 3)

# Kernel 2: Diagonal weighting
kernel2 <- matrix(c(1, 0, 1,
                    0, 2, 0,
                    1, 0, 1), nrow = 3)

# Apply filters
convolutionExample  <- convolution2D(dataMatrix = myMatrix, kernel = kernel1)
convMeanExample     <- convolutionMean(dataMatrix = myMatrix, kernel = kernel2)
convQuantileExample <- convolutionQuantile(dataMatrix = myMatrix, kernel = kernel2, x = 0.7)

# Make plots
par(mar = c(0, 0.5, 0, 0.5), oma = c(0, 0, 2, 0), mfrow = c(3, 1))

convolutionExample[convolutionExample < 0] <- 0
image(convolutionExample, col = cols, axes = FALSE)
mtext(text = "2D convolution", side = 1, line = -1.5, col = "white", font = 2, adj = 0.99)

convMeanExample[convMeanExample < 0] <- 0
image(convMeanExample, col = cols, axes = FALSE)
mtext(text = "2D mean convolution", side = 1, line = -1.5, col = "white", font = 2, adj = 0.99)

convQuantileExample[convQuantileExample < 0] <- 0
image(convQuantileExample, col = cols, axes = FALSE)
mtext(text = "2D median convolution", side = 1, line = -1.5, col = "white", font = 2, adj = 0.99)
```


### Median-filter asociated functions
```{r, eval=FALSE}
# # Add some noise (NA) to the image (matrix)
# naIndex <- sample(x = seq(prod(dim(myMatrix))), size = as.integer(0.4*prod(dim(myMatrix))), replace = FALSE)
# myMatrix[naIndex] <- NA

# Build kernel
radius <- 3

# Apply filters
meanfilterExample     <- meanFilter(dataMatrix = wbImage, radius = radius)
quantilefilterExample <- quantileFilter(dataMatrix = wbImage, radius = radius, x = 0.1)
medianfilterExample   <- medianFilter(dataMatrix = wbImage, radius = radius)

```

For compare results, we will plot both data (original and filtered) using `image` function, as shows in figures 1 and 2.

**Original**

```{r, message=FALSE, fig.height=3, fig.width=5.33, fig.cap = "Figure 2: Original matrix", results='hide', fig.pos="h", echo=FALSE}
par(mar = rep(0, 4), mfrow = c(1, 1))
image(myMatrix, col = cols)
```


**Filtered**

```{r, message=FALSE, fig.height=9, fig.width=5.33, fig.cap = "Figure 1: Filtered matrices", results='hide', fig.pos="h", echo=FALSE}
# Build kernel
radius <- 3

# Apply filters
meanfilterExample     <- meanFilter(dataMatrix = wbImage, radius = radius)
quantilefilterExample <- quantileFilter(dataMatrix = wbImage, radius = radius, x = 0.1)
medianfilterExample   <- medianFilter(dataMatrix = wbImage, radius = radius)

# Make plots
par(mar = c(0, 0.5, 0, 0.5), oma = c(0, 0, 2, 0), mfrow = c(3, 1))

meanfilterExample[meanfilterExample < 0] <- 0
image(meanfilterExample, col = cols, axes = FALSE)
mtext(text = "Mean filter", side = 1, line = -1.5, col = "white", font = 2, adj = 0.99)

quantilefilterExample[quantilefilterExample < 0] <- 0
image(quantilefilterExample, col = cols, axes = FALSE)
mtext(text = "Quantile filter", side = 1, line = -1.5, col = "white", font = 2, adj = 0.99)

medianfilterExample[medianfilterExample < 0] <- 0
image(medianfilterExample, col = cols, axes = FALSE)
mtext(text = "2D median filter", side = 1, line = -1.5, col = "white", font = 2, adj = 0.99)
```

## Kernel application

In the field of image processing, one of the tools most commonly used are the **convolutions**, which consist of operations between two arrays: The array of image data and small matrices, called kernels, which weights each pixel values with its corresponding neighborhood. Different kernels produce different effects, for instance blur, shifted images (right, left, up or down), sharpening, etc.

The larger the radius, more pixels remain unanalyzed at the edges. In addition, `kernel` argument of `convolution2D` is a matrix, but unlike the `meanFilter2D` and `meanFilter2D` functions, it should not necessarily have square dimensions. So, if the user specifies a rectangular kernel, the function will use the neighborhood as the Figure 3 shows. 

```{r, message=FALSE, fig.height=2, fig.width=5.7, fig.cap = "Figure 3: Neighborhood kernel application for different kernel dimensions. Black dot indicates the position of the pixel over the filter will be applied. Arrows indicates the direction of filter application", results='hide', fig.pos="h", echo=FALSE}
par(mar = rep(0, 4), mfrow = c(1, 3), xaxs = "i", yaxs = "i")

xlim <- c(-0.15, 1.1)
ylim <- c(-0.01, 1.3)

plotKernels <- function(dim1 = 3, dim2 = 3, showArrows = FALSE){
  plot(1, 1, pch = NA, axes = FALSE, xlim = xlim, ylim = ylim, xlab = NA, ylab = NA)
  polygon(x = c(0, 1, 1, 0), y = c(0, 0, 1, 1))
  
  delay1 <- ifelse(dim1 %% 2 != 0, 0, 1/(2*dim1))
  delay2 <- ifelse(dim2 %% 2 != 0, 0, 1/(2*dim2))
  
  mtext(text = paste0(dim1, "x", dim2, " kernel"), side = 3, line = -2)
  
  for(i in seq(dim1 - 1)){
    lines(c(0, 1), 1 - rep(i/dim1, 2))
  }
  
  for(i in seq(dim2 - 1)){
    lines(rep(i/dim2, 2), c(0, 1))
  }
  
  # abline(h = seq(0, 1, length.out = dim1 + 1))
  # abline(v = seq(0, 1, length.out = dim2 + 1))
  
  if(isTRUE(showArrows)){
    arrows(x0 = -0.05, y0 = 1.05, x1 = -0.05, y1 = 0.8, length = 0.1)
    arrows(x0 = -0.05, y0 = 1.05, x1 = 0.25, y1 = 1.05, length = 0.1)
  }
  
  points(0.5 - delay2, 0.5 + delay1, pch = 16)
}

# 3x3 kernel
plotKernels(3, 3, showArrows = TRUE)

# 3x4 kernel
plotKernels(3, 4)

# 4x5 kernel
plotKernels(4, 5)
```



