---
title: "Brain Imaging Data"
output:
  html_document:
    toc: true
    toc_float: true
    fig_width: 5
    fig_height: 5
  pdf_document:
    toc: true
    highlight: tango
    fig_width: 3
    fig_height: 3
  rmarkdown::html_vignette:
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{2 Brain Imaging Data}
  %\VignetteEncoding[UTF-8]{inputenc}
---

<!--
%\VignetteEngine{knitr::knitr}
%\\VignetteIndexEntry{mand}
-->

`r Sys.Date()` @Atsushi Kawaguchi


```{r echo=FALSE, error=FALSE}
options(warn=-1)
knitr::opts_chunk$set(eval = FALSE)
```

In this vignette, the output is omitted.
Please refer to the following book for the output.

Kawaguchi A. (2021). Multivariate Analysis for Neuroimaging Data. CRC Press.


# Affine transformation #####

The affine transformation of the figure is explained with R code. Basically, it is a matrix operation.
At the beginning a matrix whose elements are the coordinates to be transformed is prepared.

We begin by preparing a rectangle in a two-dimensional plane as an example. The coordinates of each vertex are given as a matrix.
```{r}
X = matrix(c(
-1, -1, 1, 1,
-1, 1, -1, 1
),ncol=4, byrow=TRUE)
```

A square is plotted with the coordinates in the matrix as vertices.
```{r}
par(mfrow=c(1,1), cex=0.4)
plot(NULL, xlim=c(-8, 8), ylim=c(-8, 8), xlab="x", ylab="y")

rect(xleft=min(X[1,]), ybottom=min(X[2,]), 
xright=max(X[1,]), ytop=max(X[2,]))

vs = sapply(1:4, function(x) 
paste("(", paste(X[,x], collapse=", "), ")"))
text(X[1,], X[2,], paste(LETTERS[1:4], vs), pos=c(1,3,1,3))
```


For the coordinate matrix, we perform parallel movement as follows. The rectangle before and after the move is displayed with the coordinates.
```{r}
Y = X + c(5,5)
```

```{r fig.width = 2.5, fig.height = 2.5}
par(mfrow=c(1,1), cex=0.4)
plot(NULL, xlim=c(-8, 8), ylim=c(-8, 8), xlab="x", ylab="y")
rect(xleft=min(X[1,]), ybottom=min(X[2,]), 
xright=max(X[1,]), ytop=max(X[2,]))
rect(xleft=min(Y[1,]), ybottom=min(Y[2,]), 
xright=max(Y[1,]), ytop=max(Y[2,]))
vs = sapply(1:4, function(x) 
paste("(", paste(Y[,x], collapse=", "), ")"))
text(Y[1,], Y[2,], paste(LETTERS[1:4], vs), pos=c(1,3,1,3))
```

Next, we prepare the matrix T and perform enlargement and further parallel movement.
```{r}
T = matrix(c(
2, 0,
0, 1
),ncol=2,byrow=TRUE)
s = c(-5,-5)
Y = T %*% X + s
```

```{r fig.width = 2.5, fig.height = 2.5}
par(mfrow=c(1,1), cex=0.4)
plot(NULL, xlim=c(-8, 8), ylim=c(-8, 8), xlab="x", ylab="y")
rect(xleft=min(X[1,]), ybottom=min(X[2,]), 
xright=max(X[1,]), ytop=max(X[2,]))
rect(xleft=min(Y[1,]), ybottom=min(Y[2,]), 
xright=max(Y[1,]), ytop=max(Y[2,]))
vs = sapply(1:4, function(x) 
paste("(", paste(Y[,x], collapse=", "), ")"))
text(Y[1,], Y[2,], paste(LETTERS[1:4], vs), pos=c(1,3,1,3))
```

Next, we prepare the rotation matrix T and perform the rotation and parallel movement.
```{r}
theta = pi/6
T = matrix(c(
cos(theta), -sin(theta),
sin(theta), cos(theta)
),ncol=2,byrow=TRUE)
s = c(-5, 5)
Y = T %*% X + s
```

```{r fig.width = 2.5, fig.height = 2.5}
par(mfrow=c(1,1), cex=0.4)
plot(NULL, xlim=c(-8, 8), ylim=c(-8, 8), xlab="x", ylab="y")
rect(xleft=min(X[1,]), ybottom=min(X[2,]), 
xright=max(X[1,]), ytop=max(X[2,]))
lines(Y[1, c(1, 2)], Y[2, c(1, 2)])
lines(Y[1, c(1, 3)], Y[2, c(1, 3)])
lines(Y[1, c(2, 4)], Y[2, c(2, 4)])
lines(Y[1, c(3, 4)], Y[2, c(3, 4)])
vs = sapply(1:4, function(x) 
paste("(", paste(round(Y[,x],1), collapse=", "), ")"))
text(Y[1,], Y[2,], paste(LETTERS[1:4], vs), pos=c(1,3,1,3))
```

Finally, another matrix T is prepared to perform the distortions and parallel moves.
```{r}
T = matrix(c(
1, 0,
tan(pi/6), 1
),ncol=2,byrow=TRUE)
s = c(5, -5)
Y = T %*% X + s
```

```{r fig.width = 2.5, fig.height = 2.5}
par(mfrow=c(1,1), cex=0.4)
plot(NULL, xlim=c(-8, 8), ylim=c(-8, 8), xlab="x", ylab="y")
rect(xleft=min(X[1,]), ybottom=min(X[2,]), 
xright=max(X[1,]), ytop=max(X[2,]))
lines(Y[1, c(1, 2)], Y[2, c(1, 2)])
lines(Y[1, c(1, 3)], Y[2, c(1, 3)])
lines(Y[1, c(2, 4)], Y[2, c(2, 4)])
lines(Y[1, c(3, 4)], Y[2, c(3, 4)])
vs = sapply(1:4, function(x) 
paste("(", paste(round(Y[,x],1), collapse=", "), ")"))
text(Y[1,], Y[2,], paste(LETTERS[1:4], vs), pos=c(1,3,1,3))
```


# Resize

We introduce an R code to change the resolution.
The sizechange function of the mand package is used. The image data and magnification are input. Here, the same magnification is specified for height, width, and height, but it can be different. The first case is when the magnification is set to 4.
```{r fig.width = 5, fig.height = 2}
simscale1 = 4
img1r = sizechange(template, simscale=simscale1)
dim(img1r)
coat(img1r, plane="all")
```

The next case is when the magnification is set to 1/2.
```{r fig.width = 5, fig.height = 2}
simscale1 = 1/2
img1r = sizechange(template, simscale=simscale1)
dim(img1r)
coat(img1r, plane="all")
```

When the resolution is reduced in this way, the image becomes rough and information loss occurs. In the analysis, it is a lower dimension and reduces the computational cost, but the results may not tell us much about the details.

# Smoothing 

## 1D Gauss function 

The density function of a standard normal distribution with a mean of 0 and a standard deviation of 1 is plotted. 
FWHM is then calculated as 1.17741 (in the R code, input and output are performed simultaneously by enclosing in parentheses) and is indicated by a dotted line in the figure.
```{r}
(a = sqrt(2*log(2)))
```

```{r fig.width = 4, fig.height = 3}
plot(dnorm, -4, 4, xaxt = "n", ylab="f(x)")
axis(1, -4:4, -4:4)
lines(c(0,0), c(-1,dnorm(0)), lty=2, col=1)
lines(c(-a,a), c(dnorm(0)/2,dnorm(0)/2), lty=2, col=1)
lines(c(a,a), c(-1,dnorm(a)), lty=2, col=2)
lines(-c(a,a), c(-1,dnorm(a)), lty=2, col=2)
```


## Preparation

The example data is generated from one-dimensional random field by the following code.
```{r fig.width = 7, fig.height = 3.5}
set.seed(1)
n = 20
x = 1:n
y = abs(rnorm(n)); names(y) = x
barplot(y)
```

In order to implement examples below, two functions are defined here. 
The first is the gauss function with the coordinates x and the center coordinate x1. 
```{r}
g1 = function(x, x1, FWHM){
sigma = FWHM / 2*sqrt(2*log(2))
d1 = dnorm(x, x1, sigma)
names(d1) = x
d1
}
```

The second is the smoothing function using the gauss function.
```{r}
gsmooth = function(x, y, FWHM){
sigma = FWHM / 2*sqrt(2*log(2))
sy = sapply(x, function(x1) 
weighted.mean(y, 
dnorm(x, x1, sigma)/sum(dnorm(x, x1, sigma))) )
names(sy) = x; sy; }
```


## Smoothing flow

The original data is in the first row of the plot with data (signal) in the vertical line and the coordinates in the horizontal line,
the gauss function with the 5th coordinate as the center in the second row,
the convoluted data is in the third row, and those summation is in the fourth row.
```{r fig.width = 5, fig.height = 5}
x1 = 5
col1 = rep(1, n); col1[x1] = 2
par(mfrow=c(4,1), mar=c(3,4,1,4))
barplot(y, col=col1)
barplot(g1(x,x1,2), col=col1, ylim=c(0, 0.4))
barplot(y*g1(x,x1,2), col=col1)
sy = rep(0, n); sy[x1] = sum(y*g1(x,x1,2)); names(sy) = x
barplot(sy, col=col1, ylim=c(0,2))
```

The above flow is applied to all coordinate as the center and the smoothed data with FWHM=2 was obtained in the second row panel.
```{r fig.width = 6, fig.height = 5}
col2 = rep(1, n)
par(mfrow=c(2,1))
barplot(y, col=col2, main="Original")
f1 = 2
sy = gsmooth(x, y, f1)
barplot(sy, main=paste("FWHM =", f1), col=col2, ylim=c(0,2))
```

## FWHM and smoothed data  ##########
The Gauss functions with different FWHM values are displayed in the left side.
The different FWHMs was examined to smooth the original data and obtained the smoothed data as follows.
```{r fig.width = 5, fig.height = 5}
f1s = c(2, 4, 8)
layout(cbind(c(8,1,2,3), c(4:7)))
par(mar=c(3,4,2,4))
for(f1 in f1s){
barplot(g1(x,8,f1), main=paste("FWHM =", f1), ylim=c(0,0.4))
}
barplot(y, main="Original")
for(f1 in f1s){
sy = gsmooth(x, y, f1)
barplot(sy, main=paste("FWHM =", f1), ylim=c(0,2))
}
```

The larger FWHM values yielded more smooth signals (similar values in the neighborhood).