---
title: "TCIU Predictive Analytics: "
subtitle: "Validation using 4D Complex (mag,phase) and Real (solely-mag reconstruction) fMRI data"
author: "SOCR Team"
date: "`r format(Sys.time(),'%m/%d/%Y')`"
output: 
  html_document:
    theme: spacelab
    highlight: tango
    includes:
      before_body: TCIU_header.html
    toc: true
    number_sections: true
    toc_depth: 2
    toc_float:
      collapsed: false
      smooth_scroll: true
    code_folding: hide
---

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

# Package Loading and Data Manipulation
```{r,message=FALSE}
library(TCIU)
library(DT)
library(R.matlab)
library(AnalyzeFMRI)
library(ggfortify)
mat1 = readMat("subj1_run1_complex_all.mat")
bigim1 = mat1$bigim[,64:1,,]
bigim1_mod = Mod(bigim1)
smoothmod = GaussSmoothArray(bigim1_mod,sigma = diag(3,3))
# dim(bigim1) = 64 64 40
# bigim1 contains the complex image space
# dimensions are 64x*64y*40z*160t, corresponding to x,y,z,time
load("mask.rda") # load the 3D nifit data of the mask
load("mask_label.rda") # load the 3D nifit data of the mask with labels
load("mask_dict.rda") # load the label index and label name
label_index = mask_dict$index
label_name = as.character(mask_dict$name)
label_mask = mask_label
load("hemody.rda") # load the hemodynamic contour of the brain
```

# Time-series graphs

## Interactive time-series visualization

The function `fmri_time_series` is used to create four interactive time series graphs for the real, imaginary, magnitude, and phase parts for the fMRI spacetime data. We can also add a reference plotly object to the plot. This function is based on `GTSplot` function from package `TSplotly`.


### Example: X: 44, Y: 42, Z: 33
```{r,fig.height=10,fig.width=10}
reference_plot = sample[[4]]
fmri_time_series(bigim1, c(44,42,33), is.4d = TRUE, ref = reference_plot)
```

# Kime-series/Kime-surfaces (spacekime analytics protocol)
This Rmd notebook uses a time-series simulation to illustrate how to transform the fMRI data at a fixed voxel location into a kime-surface (kime-series).

Notes:

Validate all steps in this transformation protocol and finalize this 3D visualization.

Try it with real fMRI data at brain voxel locations associated with the finger-tapping task.

## Pseudo-code
+ Randomly generate 8 $phi=\phi$ kime-phases for each of the 10 time radii. This yields an $8\times 10$ array (*phi_8_vec*) of kime phase directions. These directions can be obtained by different strategies, e.g., (1) uniform or Laplace distribution sampling over the interval $[-\pi:\pi)$, (2) randomly sampling with/without replacement from known kime-phases obtained from similar processes, etc.
+ Optionally, order all kime-phases (rows) from small to large for each column.
+ Fix the $\nu=(x,y,z)$ voxel location and extract the fMRI time-course $fMRI_{\nu}(t), \forall 1\leq t\leq 160$.
+ For binary stimuli (e.g., activation (ON) and rest (OFF) event-related design), we can split the 160-long fMRI series into 80 ON (Activation) and 80 OFF (rest) states, or sub-series.
+ Construct a data-frame with 160 rows and 4 columns; time (1:10), phases (8), states (2), and fMRI_value (Complex or Real intensity).
+ Convert the long DF representing fMRI_ON and fMRI_OFF from their native (old) polar coordinates to the (new) Cartesian coordinates, using polar transformations.
+ Check for visual (graphical) and numeric differences between the fMRI intensities during the ON vs. OFF states
+ Spatially smooth (blur) the matrices (in 2D) to reduce noise make them more representative of the process. May also need to temper the magnitude of the raw fMRI intensities, which can have large range.
+ Generate the `plot_ly` *text labels* that will be shown on mouse hover (pop-up dialogues) over each kime-surface/kime-series. These text-labels are stored in Cartesian coordinates $(-10\leq x\leq 10,-10\leq y\leq 10)$, but are computed using the polar coordinates $(1\leq t\leq 10,-\pi\leq \phi<\pi)$ and the polar-to-Cartesian transformation. This the labels are $21\times21$ arrays including the triple $(t, \phi, fMRIvalue)$. Separate text-labels are generated for each kime-surface (ONN vs. OFF stimuli).
+ Generate the $21\times21$ kime-domain Cartesian grid by polar-transforming the polar coordinates $(t,\phi)$ into Cartesian counterparts $(x,y)$.
+ Interpolate the fMRI intensities (natively anchored at $(t,\phi)$) to $fMRI(x,y), \forall -11\leq x,y \leq 11, x^2+y^2\leq 10$.
+ Use `plot_ly` to display in 3D the kime-series as 2D manifolds (kime-surfaces) over the Cartesian domain.


## Function main steps illustration: Kime-surfaces / kime-series construction

### Generate the kime-phases
Randomly generate 8 $phi=\phi$ kime-phases for each of the 10 time radii. This yields an $8\times 10$ array (*phi_8_vec*) of kime phase directions. These directions can be obtained by different strategies, e.g., (1) uniform or Laplace distribution sampling over the interval $[-\pi:\pi)$, (2) randomly sampling with/without replacement from known kime-phases obtained from similar processes, etc.

Optionally, order all kime-phases (rows) from small to large for each column.

```{r}
# plot Laplacian
ggfortify::ggdistribution(extraDistr::dlaplace, seq(-pi, pi, 0.01), m=0, s=0.5)

# randomly generate 8 phi kime-phases for each of the 10 time radii
phi_8_vec <- matrix(NA, ncol=10, nrow = 8)
for (t in 1:10) { 
  # for a given t, generate 8 new phases
  set.seed(t);
  phi_8_vec[ ,t] <-
    extraDistr::rlaplace(8,mu=0,sigma=0.5)
  # rank-order the phases for consistency
  # within the same foliation leaf
  phi_8_vec[ ,t] <- sort(phi_8_vec[ ,t])
  # force phases in [-pi: pi)
  for (i in 1:8) {
    if (phi_8_vec[i,t] < -pi) 
      phi_8_vec[i,t] <- -pi
    if (phi_8_vec[i,t] >= pi) 
      phi_8_vec[i,t] <- pi
  }
}
# order all kime-phases (rows) from small to large for each column
# phi_8_vec_ordered <- apply(phi_8_vec, 2, sort)
```


### Split the fMRI series according to the event-related design
For binary stimuli (e.g., activation (ON) and rest (OFF) event-related design), we can split the 160-long fMRI series into 80 ON (Activation) and 80 OFF (rest) states, or sub-series.

```{r}
fMRI_ON<-bigim1_mod[40,42,33,][c(rep(TRUE,10),rep(FALSE,10))]
fMRI_OFF<-bigim1_mod[40,42,33,][c(rep(FALSE,10),rep(TRUE,10))]
```

### Preprocess the matrices
Spatially smooth (blur) the matrices (in 2D) to reduce noise make them more representative of the process. May also need to temper the magnitude of the raw fMRI intensities, which can have large range.

```{r}
# Convert the long DF representing fMRI_ON and fMRI_OFF from polar coordinates to Cartesian coordinates
matrix_ON <- matrix(0, nrow = 21, ncol = 21) 
matrix_OFF <- matrix(0, nrow = 21, ncol = 21) 
for (t in 1:10) {
  for (p in 1:8) {
    x = 11+t*cos(phi_8_vec[p,t])
    y = 11+t*sin(phi_8_vec[p,t])
    matrix_ON[x,y]  <- fMRI_ON[(p-1)*10 +t]
    matrix_OFF[x,y] <- fMRI_OFF[(p-1)*10 +t]
  }
}
# smooth/blur the matrices
matrix_ON_smooth <- (1/10000)*as.matrix(spatstat::blur(spatstat::as.im(matrix_ON), sigma=0.5))
matrix_OFF_smooth <- (1/10000)*as.matrix(spatstat::blur(spatstat::as.im(matrix_OFF), sigma=0.5))
```

### Generate `plotly` labels
Generate the `plot_ly` text labels that will be shown on mouse hover (pop-up dialogues) over each kime-surface/kime-series. These text-labels are stored in Cartesian coordinates $(-10\leq x\leq 10,-10\leq y\leq 10)$, but are computed using the polar coordinates $(1\leq t\leq 10,-\pi\leq \phi<\pi)$ and the polar-to-Cartesian transformation. This the labels are $21\times21$ arrays including the triple $(t, \phi, fMRIvalue)$. Separate text-labels are generated for each kime-surface (ONN vs. OFF stimuli).

```{r}
# fix the plot_ly Text Lables
x <- vector()
y <- vector()
i <- 1
for (t in 1:10) {
  for (p in 1:8) {
    x[i] = 11+t*cos(phi_8_vec[p,t])
    y[i] = 11+t*sin(phi_8_vec[p,t])
    i <- i+1
  }
}
```

### Cartesian representation
Generate the $21\times21$ kime-domain Cartesian grid by polar-transforming the polar coordinates $(t,\phi)$ into Cartesian counterparts $(x,y)$.

```{r}
hoverText <- cbind(x=1:21, y=1:21, height=as.vector(t(matrix_ON_smooth))) # tail(mytext)
custom_txt <- matrix(NA, nrow=21, ncol=21)
hoverTextOFF <- cbind(x=1:21, y=1:21, height=as.vector(t(matrix_OFF_smooth))) # tail(mytext)
custom_txtOFF <- matrix(NA, nrow=21, ncol=21)

for (x in 1:21) {
   for (y in 1:21) {
     t = sqrt((x-11)^2 + (y-11)^2)
     p = atan2(y-11, x-11)
     custom_txt[x,y] <- paste(' fMRI: ', round(hoverText[(x-1)*21+y, 3], 3),
                    '\n time: ', round(t, 0),
                    '\n phi: ', round(p, 2)) 
     custom_txtOFF[x,y] <- paste(' fMRI: ', round(hoverTextOFF[(x-1)*21+y, 3], 3),
                    '\n time: ', round(t, 0),
                    '\n phi: ', round(p, 2)) 
   }
}
```

### Cartesian space interpolation
Interpolate the fMRI intensities (natively anchored at $(t,\phi)$) to $fMRI(x,y), \forall -11\leq x,y \leq 11, x^2+y^2\leq 10$.

```{r}
xx2 <- 11 + c(-10:10) %o% cos(seq(-pi, pi, 2*pi/20))
yy2 <- 11 + c(-10:10) %o% sin(seq(-pi, pi, 2*pi/20))
#zz2 <- as.vector(matrix_ON_smooth) %o% rep(1, 21*21)
zz2 <- matrix_ON_smooth
ww2 <- matrix_OFF_smooth
dd2 <- matrix_ON_smooth-matrix_OFF_smooth

#plot 2D into 3D and make the text of the diameter (time), height (r), and phase (phi)
f <- list(family = "Courier New, monospace", size = 18, color = "black")
x <- list(title = "k1", titlefont = f)
y <- list(title = "k2", titlefont = f)
z <- list(title = "fMRI Kime-series", titlefont = f)
```

## Function result

The following parts are the main output of our functions

### Generate a long data-frame
Construct a data-frame with 160 rows and 4 columns; time (1:10), phases (8), states (2), and fMRI_value (Complex or Real intensity). Then, convert the long DF representing fMRI_ON and fMRI_OFF from their native (old) polar coordinates to the (new) Cartesian coordinates, using polar transformations.

```{r}
datatable(fmri_kimesurface(bigim1_mod,c(44,42,33))[[1]])
```

### Display kime-surfaces
Use `plot_ly` to display in 3D the kime-series as 2D manifolds (kime-surfaces) over the Cartesian domain.

```{r}
fmri_kimesurface(bigim1_mod,c(44,42,33))[[2]]
fmri_kimesurface(bigim1_mod,c(44,42,33))[[3]]
fmri_kimesurface(bigim1_mod,c(44,42,33))[[4]]
```

# Interactive plotly Example

## Plotly method: interactive way

The function `fmri_image` is used to create images for front view, side view, and top view of the fMRI image.

```{r,warning=FALSE}
bigim1_mask<-bigim1_mod
for (i in 1:160) {
  bigim1_mask[,,,i]<-bigim1_mask[,,,i]*mask
}
```

Manually examples:

The plot without mask restriction.
```{r,warning=FALSE,fig.height=6,fig.width=7.5}
fmri_image(bigim1_mod, option = "manually", voxel_location = c(40,40,30), time = 4)
```

The plot with mask restriction.
```{r,warning=FALSE,fig.height=6,fig.width=7.5}
fmri_image(bigim1_mask, option = "manually", voxel_location = c(40,40,30), time = 4)
```

## Forcasting with time series

The function `fmri_ts_forecast` will fit the ARIMA model for each location. This function is based on `TSplot_gen` function from package `TSplotly`. 

```{r,fig.width=10,fig.height=4.5}
fmri_ts_forecast(smoothmod,c(41,44,33))
```

# Motor area detection

## fMRI data simulation

The function `fmri_simulate_func` is used to simulate the real-valued fMRI data with the specified dimension.

```{r}
fmri_generate = fmri_simulate_func(dim_data = c(64, 64, 40), mask = mask, 
                                   ons = c(1, 21, 41, 61, 81, 101, 121, 141), 
                                   dur = c(10, 10, 10, 10, 10, 10, 10, 10))
```


## Stimulus detection

An integrated function `fmri_stimulus_detect` which can apply multiple methods is used to detect motor area.

"t-test" and "wilcox-test" can be applied to all real-valued fMRI data.
"HotellingsT2","Wilk's Lambda" and "gLRT" method can be applied to all complex-valued fMRI data.
"on_off_diff" and "HRF" method can be applied to 4D real-valued fMRI data where the first method is calculating on-off difference period polar volume to get p-values and the second method is using hemodynamic response function and linear regression.
The post-hoc stat-mapping filtering can be also applied to the calculated p-values if specified in the parameter or use function `fmri_post_hoc` to do the process.


### Examples
```{r}
# statistical method HRF needs parameter ons and dur
pval1 = fmri_stimulus_detect(fmridata= bigim1_mod, mask = mask,
                             stimulus_idx = c(1:160)[rep(c(TRUE,FALSE), c(10,10))],
                             method = "HRF" , 
                             ons = c(1, 21, 41, 61, 81, 101, 121, 141), 
                             dur = c(10, 10, 10, 10, 10, 10, 10, 10) )

# statistical method t-test for real-valued fMRI data
pval2 = fmri_stimulus_detect(fmridata= bigim1_mod, mask = mask,
                             stimulus_idx = c(1:160)[rep(c(TRUE,FALSE), c(10,10))],
                             method = "t-test")

# statistical method Wilk's Lambda for complex-valued data
pval3 = fmri_stimulus_detect(fmridata = bigim1, mask = mask,
                             stimulus_idx = c(1:160)[rep(c(TRUE,FALSE), c(10,10)) ], 
                             method = "Wilks-Lambda" )

# do the fdr correction and the spatial clustering
# pval4 is the pval1 after the post-hoc processing
pval4 = fmri_post_hoc(pval1, fdr_corr = "fdr",
                    						 spatial_cluster.thr = 0.05,
                    						 spatial_cluster.size = 5, 
                    						 show_comparison = FALSE)

```

```{r}
summary(pval1)
summary(pval2)
summary(pval3)
summary(pval4)
```


# Motor area visualization

## Viusalization and comparison of p-value

### 3D viusalization for p-value

```{r fig.width = 9, fig.align = "center", warning=FALSE}
pval1_3d = fmri_3dvisual(pval1, mask, p_threshold = 0.05, method="scale_p", multi_pranges=TRUE, title="HRF method")

pval1_3d$plot
```

```{r fig.width = 9, fig.align = "center", warning=FALSE}
pval4_3D = fmri_3dvisual(pval4, mask, p_threshold = 0.05, method="scale_p",
                         multi_pranges=TRUE, title="HRF method after post-hoc")

pval4_3D$plot
```

### 2D viusalization for p-value

Generate the 2D plot of the p-value from sagittal, coronal and axial view.

Plot without hemodynamic contour.

```{r warning=FALSE}
for(axis in c("x", "y", "z")){
  axis_i = switch(axis, 
                  "x" = {35},
                  "y" = {30},
                  "z" = {22})
  print(fmri_2dvisual(pval1, list(axis, axis_i), 
                      hemody_data=NULL, mask=mask, 
                      p_threshold = 0.05, legend_show = TRUE, 
                      method = "scale_p",
                      color_pal = "YlOrRd", multi_pranges=TRUE))
}
			
```

Plot with hemodynamic contour.

```{r warning=FALSE}
for(axis in c("x", "y", "z")){
  axis_i = switch(axis, 
                  "x" = {35},
                  "y" = {30},
                  "z" = {22})
  print(fmri_2dvisual(pval1, list(axis, axis_i), 
                      hemody_data=hemody, mask=mask, 
                      p_threshold = 0.05, legend_show = TRUE, 
                      method = "scale_p",
                      color_pal = "YlOrRd", multi_pranges=TRUE))
}
```

## Comparison of performance of different methods on the same fMRI data

### 3d p-value comparison
`fmri_pval_comparison_3d` is to visualize two p-value data to see how they perform differently on detecting stimulated parts by 3D plot. Here we compare the difference of stimulated parts of two different fMRI data with the same mask, since our original fMRI is too big to use here for different statistical test.

```{r fig.width = 9, fig.align = "center", warning=FALSE}
fmri_pval_comparison_3d(list(pval1, pval2), mask, 
    				            list(0.05, 0.05), list("scale_p", "scale_p"), 
    				            multi_pranges=FALSE)
```

### 2d p-value comparison
`fmri_pval_comparison_2d` is to visualize whatever number of p value (generated by different statistical tests on the same fMRI data) to see their difference by 2D plot. For simplicity here we only compare two different 3D array p value data.


```{r fig.width = 9, warning=FALSE}
fmri_pval_comparison_2d(list(pval2, pval1), 
      				           list('t-test', 'HRF'),
      				           list(list(35, 33, 22), list(40, 26, 33)), 
      				           hemody_data = NULL, 
      				           mask = mask, p_threshold = 0.05, 
      				           legend_show = FALSE, method = 'scale_p',
      				           color_pal = "YlOrRd", multi_pranges=FALSE)
```

# Three-phase ROI Analysis

## Phase1: Detect Potential Activated ROI
```{r eval=FALSE}
phase1_pval = fmri_ROI_phase1(bigim1_mod, mask_label, mask_dict,
                              stimulus_idx = c(1:160)[rep(c(TRUE,FALSE), c(10,10))])$all_ROI$pval_t
```

## Phase2: ROI-Based Tensor-on-Tensor Regression
```{r eval=FALSE}
phase2_pval = fmri_ROI_phase2(fmridata = bigim1_mod, label_mask = mask_label, label_dict = mask_dict, 
                              stimulus_idx = c(1, 21, 41, 61, 81, 101, 121, 141),
                              stimulus_dur = c(10, 10, 10, 10, 10, 10, 10, 10),
                              rrr_rank = 3, fmri.design_order = 2,
                              fmri.stimulus_TR = 3, method = "t_test")
```

## Phase3: FDR Correction and Spatial Clustering
```{r eval=FALSE}
phase3_pval = fmri_post_hoc(phase2_pval , fdr_corr = "fdr",
                            spatial_cluster.thr = 0.05,
                            spatial_cluster.size = 5, 
                            show_comparison = FALSE)
```
                           
## 3D visualization based on the activated areas by regions

```{r fig.width = 9, fig.align = "center", warning=FALSE}
fmri_3dvisual_region(TCIU::phase1_pval, label_mask, label_index,
                     label_name, title = "phase1 p-values", rank = "value")
``` 


```{r fig.width = 9, fig.align = "center", warning=FALSE}
fmri_3dvisual_region(TCIU::phase1_pval, label_mask, label_index,
                     label_name, 5, title = "phase1 top five p-values", rank = "value")
``` 


```{r fig.width = 9, fig.align = "center", warning=FALSE}
# for 3D visualization, user needs to include empty region in the label
label_index = c(0, label_index)
label_name = c("empty", label_name)
fmri_3dvisual_region(TCIU::phase2_pval, label_mask, label_index,
                    label_name, title = "phase2 p-values")
``` 


```{r fig.width = 9, fig.align = "center", warning=FALSE}
fmri_3dvisual_region(list(TCIU::phase2_pval,TCIU::phase3_pval), label_mask, label_index,
                    label_name, title = "phase2&3 p-values")
```


<!--html_preserve-->
<div>
    	<footer><center>
			<a href="http://www.socr.umich.edu/">SOCR Resource</a>
				Visitor number <img class="statcounter" src="https://c.statcounter.com/5714596/0/038e9ac4/0/" alt="Web Analytics" align="middle" border="0">
				<script type="text/javascript">
					var d = new Date();
					document.write(" | " + d.getFullYear() + " | ");
				</script> 
				<a href="http://socr.umich.edu/img/SOCR_Email.png"><img alt="SOCR Email"
	 			title="SOCR Email" src="http://socr.umich.edu/img/SOCR_Email.png"
	 			style="border: 0px solid ;"></a>
	 		 </center>
	 	</footer>

	<!-- Start of StatCounter Code -->
		<script type="text/javascript">
			var sc_project=5714596; 
			var sc_invisible=1; 
			var sc_partition=71; 
			var sc_click_stat=1; 
			var sc_security="038e9ac4"; 
		</script>
		
		<script type="text/javascript" src="https://www.statcounter.com/counter/counter.js"></script>
	<!-- End of StatCounter Code -->
	
	<!-- GoogleAnalytics -->
		<script src="https://www.google-analytics.com/urchin.js" type="text/javascript"> </script>
		<script type="text/javascript"> _uacct = "UA-676559-1"; urchinTracker(); </script>
	<!-- End of GoogleAnalytics Code -->
</div>
<!--/html_preserve-->


