SOCR ≫ TCIU Website ≫ TCIU GitHub ≫ 

# Get this figure: fig <- get_figure("MattSundquist", 4064)
# Get this figure's data: data <- get_figure("MattSundquist", 4064)$data
# Add data to this figure: p <- add_trace(p, x=c(4, 5), y=c(4, 5), kwargs=list(filename="Klein bottle", fileopt="extend"))
# Get y data of first trace: y1 <- get_figure("MattSundquist", 4064)$data[[1]]$y

library(plotly)

1 fMRI Space-kime Example

In this example, we demonstrate spacekime data analytics using real 4D fMRI data (\(x=64 \times y=64\times z=21\times t=180\)). For simplicity of the presentation, analysis, and visualization, we will focus on a 2D time-series of the entire 4D fMRI hypervolume. In other words, we’ll (artificially) reduce the native 3D space (\({\bf{x}}=(x,y,z)\in \mathbb{R}^3\)) to 2D (\({\bf{x}}=(x,y)\in \mathbb{R}^2\subseteq \mathbb{R}^3\)) by focusing only on mid-axial/transverse slice through the brain (\(z=11\)). More details are provided in DSPA Chapter 3.

# FFT SHIFT
fftshift <- function(img_ff, dim = -1) {

  rows <- dim(img_ff)[1]    
  cols <- dim(img_ff)[2]
  # planes <- dim(img_ff)[3]

  swap_up_down <- function(img_ff) {
    rows_half <- ceiling(rows/2)
    return(rbind(img_ff[((rows_half+1):rows), (1:cols)], img_ff[(1:rows_half), (1:cols)]))
  }

  swap_left_right <- function(img_ff) {
    cols_half <- ceiling(cols/2)
    return(cbind(img_ff[1:rows, ((cols_half+1):cols)], img_ff[1:rows, 1:cols_half]))
  }
  
  #swap_side2side <- function(img_ff) {
  #  planes_half <- ceiling(planes/2)
  #  return(cbind(img_ff[1:rows, 1:cols, ((planes_half+1):planes)], img_ff[1:rows, 1:cols, 1:planes_half]))
  #}

  if (dim == -1) {
    img_ff <- swap_up_down(img_ff)
    return(swap_left_right(img_ff))
  }
  else if (dim == 1) {
    return(swap_up_down(img_ff))
  }
  else if (dim == 2) {
    return(swap_left_right(img_ff))
  }
  else if (dim == 3) {
    # Use the `abind` package to bind along any dimension a pair of multi-dimensional arrays
    # install.packages("abind")
    library(abind)
    
    planes <- dim(img_ff)[3]
    rows_half <- ceiling(rows/2)
    cols_half <- ceiling(cols/2)
    planes_half <- ceiling(planes/2)
    
    img_ff <- abind(img_ff[((rows_half+1):rows), (1:cols), (1:planes)], 
                    img_ff[(1:rows_half), (1:cols), (1:planes)], along=1)
    img_ff <- abind(img_ff[1:rows, ((cols_half+1):cols), (1:planes)], 
                    img_ff[1:rows, 1:cols_half, (1:planes)], along=2)
    img_ff <- abind(img_ff[1:rows, 1:cols, ((planes_half+1):planes)], 
                    img_ff[1:rows, 1:cols, 1:planes_half], along=3)
    return(img_ff)
  }
  else {
    stop("Invalid dimension parameter")
  }
}

ifftshift <- function(img_ff, dim = -1) {

  rows <- dim(img_ff)[1]    
  cols <- dim(img_ff)[2]    

  swap_up_down <- function(img_ff) {
    rows_half <- floor(rows/2)
    return(rbind(img_ff[((rows_half+1):rows), (1:cols)], img_ff[(1:rows_half), (1:cols)]))
  }

  swap_left_right <- function(img_ff) {
    cols_half <- floor(cols/2)
    return(cbind(img_ff[1:rows, ((cols_half+1):cols)], img_ff[1:rows, 1:cols_half]))
  }

  if (dim == -1) {
    img_ff <- swap_left_right(img_ff)
    return(swap_up_down(img_ff))
  }
  else if (dim == 1) {
    return(swap_up_down(img_ff))
  }
  else if (dim == 2) {
    return(swap_left_right(img_ff))
  }
  else if (dim == 3) {
    # Use the `abind` package to bind along any dimension a pair of multi-dimensional arrays
    # install.packages("abind")
    library(abind)
    
    planes <- dim(img_ff)[3]
    rows_half <- floor(rows/2)
    cols_half <- floor(cols/2)
    planes_half <- floor(planes/2)
    
    img_ff <- abind(img_ff[1:rows, 1:cols, ((planes_half+1):planes)], 
                    img_ff[1:rows, 1:cols, 1:planes_half], along=3)
    img_ff <- abind(img_ff[1:rows, ((cols_half+1):cols), (1:planes)], 
                    img_ff[1:rows, 1:cols_half, (1:planes)], along=2)
    img_ff <- abind(img_ff[((rows_half+1):rows), (1:cols), (1:planes)], 
                    img_ff[(1:rows_half), (1:cols), (1:planes)], along=1)
    return(img_ff)
  }
  else {
    stop("Invalid dimension parameter")
  }
}
# install EBImage
# source("https://bioconductor.org/biocLite.R")
# biocLite("EBImage")
library(EBImage)
require(brainR)
library(spatstat) 
library(plotly)

# 1. download the 4D fMRI data
fMRIURL <- "http://socr.umich.edu/HTML5/BrainViewer/data/fMRI_FilteredData_4D.nii.gz"
fMRIFile <- file.path(tempdir(), "fMRI_FilteredData_4D.nii.gz")
download.file(fMRIURL, dest=fMRIFile, quiet=TRUE)
(fMRIVolume <- readNIfTI(fMRIFile, reorient=FALSE))
## NIfTI-1 format
##   Type            : nifti
##   Data Type       : 4 (INT16)
##   Bits per Pixel  : 16
##   Slice Code      : 0 (Unknown)
##   Intent Code     : 0 (None)
##   Qform Code      : 1 (Scanner_Anat)
##   Sform Code      : 0 (Unknown)
##   Dimension       : 64 x 64 x 21 x 180
##   Pixel Dimension : 4 x 4 x 6 x 3
##   Voxel Units     : mm
##   Time Units      : sec
# dimensions: 64 x 64 x 21 x 180 ; 4mm x 4mm x 6mm x 3 sec 

fMRIVolDims <- dim(fMRIVolume); fMRIVolDims
## [1]  64  64  21 180
# time_dim <- fMRIVolDims[4]; time_dim ## 180

# 2. extract the time-corse of 2D mid-axial slice (3D) hypervolume
fMRI_2D_z11 <- fMRIVolume[ , , 11, ]; dim(fMRI_2D_z11)   # 64  64 180
## [1]  64  64 180
# 3. Some plots 
### 2D
image(as.im(fMRI_2D_z11[, , 1]), col=grey(fMRI_2D_z11[, , 1]/max(fMRI_2D_z11[, , 1]), alpha = 1))

### 3D Plotly
img_t90 <- as.matrix(blur(as.im(fMRI_2D_z11[, , 90]), sigma=0.1)) # the smoothed version of the 2D image (t=90)
image(as.im(img_t90))

img_t1 <- -20000+ as.matrix(blur(as.im(fMRI_2D_z11[, , 1]), sigma=0.1))  # t=1 baseline
img_t180 <- 20000+ as.matrix(blur(as.im(fMRI_2D_z11[, , 180]), sigma=0.1))  # t=180 end

# Plot the image surfaces
p <- plot_ly(z=img_t90, type="surface", showscale=FALSE) %>%
 add_trace(z=img_t1, type="surface", showscale=FALSE, opacity=0.5) %>%
 add_trace(z=img_t180, type="surface", showscale=FALSE, opacity=0.5)
p 
#### # try different levels at which to construct contour surfaces (10 fast)
# lower values yield smoother surfaces # see ?contour3d
# contour3d(fMRI_2D_z11, level = 1000, alpha = 0.1, draw = TRUE)

# multiple levels may be used to show multiple shells
# "activations" or surfaces like hyper-intense white matter
# This will take 1-2 minutes to rend!
contour3d(fMRI_2D_z11, level = c(1000, 15000), alpha = c(0.3, 0.5),
        add = TRUE, color=c("yellow", "red"))

# Plot the 4D array of imaging data in a 5x5 grid of images 
# The first three dimensions are spatial locations of the voxel (volume element) and the fourth dimension is time for this functional MRI (fMRI) acquisition. 
# image(fMRIVolume, zlim=range(fMRIVolume)*0.95)

# 4. FT of 2D slices
X1 = fft(img_t90); image(as.im(fftshift(Re(X1)))) # display(Re(X1), method = "raster")

X1_mag <- sqrt(Re(X1)^2+Im(X1)^2); image(as.im(fftshift(Re(X1_mag)))) # display(X1_mag, method = "raster")

X1_phase <- atan2(Im(X1), Re(X1)); image(as.im(fftshift(Re(X1_phase)))) # display(X1_phase, method = "raster")

##  Implicit Automated IFT
hat_X1 = Re(fft(X1, inverse = T)/length(X1)); image(as.im(hat_X1)) # display(hat_X1, method = "raster")

## Manually invert the FT (IFT) using the magnitudes and phases
Real1 = X1_mag * cos(X1_phase)
Imaginary1 = X1_mag * sin(X1_phase)
man_hat_X1 = Re(fft(Real1 + 1i*Imaginary1, inverse = T)/length(X1)); image(as.im(man_hat_X1))

# IFT fMRI-Magnitude and Nil-Phase
Real_phase0 = X1_mag * cos(0)
Imaginary_phase0 = X1_mag * sin(0)
ift_NilPhase_X1mag = Re(ifftshift(fft(Real_phase0 + 1i*Imaginary_phase0, inverse = T)/length(X1))); image(as.im(ift_NilPhase_X1mag))

Suppose the 2D fMRI time-series is represented analytically by \(f({\bf{x}},t)=f(x,y,t):\mathbb{R}^2\times \mathbb{R}^+ \longrightarrow R\) and computationally as a 3D array. Then each of these are also 3D (complex-valued) arrays: \(\hat{f}\), magnitude of the FT (\(|\hat{f}|=\sqrt{Re(\hat{f})^2+Im(\hat{f})^2}\)), and the phase-angles, \(\theta = \arctan \left (\frac{Im(\hat{f})}{Re(\hat{f})}\right )\).

We will focus on the function \(\hat{f}=(f_1, f_2, f_3)\) where the 3-rd dimension correponds to time. Specifically, we will consider the magnitude of its 3-rd dimension as \(time=|f_3|\) and we will pretend its phase is unknown, i.e., \(\theta_3 =0\). Thus, inverting the FT of the modified function \(\tilde{\hat{f}}\), where \(\theta_3 =0\), we get an estimate of kime for the original 2D fMRI time-series as \(\hat{\tilde{\hat{f}}}\).

As an observable, the time is measurable and the phase angles can either be estimated from other similar data, provided by an oracle, or fixed according to some experimental conditions. For simplicity, we’ll consider two specific instance:

  • When the time-dimension phases are indeed the actual FT phases, which are in general unknown, however, in our fMRI simulation example, they are actually computed from the original space-time fMRI time-series via the Fourier transformation, and
  • When the time-dimension phases are provided by the investigator, e.g., trivial (nil) phases or phases derived from other similar datasets.
# 5. FT of 3D time-series 
X1 = fft(fMRI_2D_z11); dim(X1); hist(Re(log(1+X1)), xlim=c(5, 20))
# Plot the centered frequency spectrum FT of the 2D time-series (in 3D), only half the frequencies are needed.
planes_half <- ceiling(dim(fMRI_2D_z11)[3]/2)

# Visualize the Simulated Original Observed Data in k-space
img1 <- fftshift(Re(log(1+X1)),3)[ , , (1:(planes_half+1))]   # apply log transform to temper the intensity range
contour3d(img1, level = c(7, 12), alpha = c(0.3, 0.5), add = TRUE, color=c("yellow", "red"), perspective=T)

## Compute the 3D Magnitude and Phase     
X1_mag <- fftshift(sqrt(Re(X1)^2+Im(X1)^2), 3)[ , , (1:(planes_half+1))] # log transform to temper the intensity range
contour3d(log(X1_mag), level = c(7, 12), alpha = c(0.3, 0.5), add = TRUE, color=c("yellow", "red"))
X1_phase <- atan2(Im(X1), Re(X1))
contour3d(X1_phase, level = c(-2, 2), alpha = c(0.3, 0.5), add = TRUE, color=c("yellow", "red"))

##  Implicit Automated IFT
hat_X1 = Re(fft(X1, inverse = T)/length(X1))
contour3d(hat_X1, level = c(1000, 15000), alpha = c(0.3, 0.5), add = TRUE, color=c("yellow", "red"))

## Manually invert the FT (IFT) using the magnitudes and phases
Real1 = fftshift(sqrt(Re(X1)^2+Im(X1)^2), 3) * cos(X1_phase) #  X1_mag * cos(X1_phase)
Imaginary1 = fftshift(sqrt(Re(X1)^2+Im(X1)^2), 3) * sin(X1_phase)
man_hat_X1 = Re(fft(Real1 + 1i*Imaginary1, inverse = T)/length(X1))
contour3d(man_hat_X1, level = c(1000, 15000), alpha = c(0.3, 0.5), add = TRUE, color=c("yellow", "red"))

# IFT fMRI-Magnitude and Nil-Phase
Real_phase0 = X1_mag * cos(0)
Imaginary_phase0 = X1_mag * sin(0)
ift_NilPhase_X1mag = Re(ifftshift(fft(Real_phase0 + 1i*Imaginary_phase0, inverse = T)/length(X1), dim=3))
contour3d(ift_NilPhase_X1mag, level = c(1000, 15000), alpha = c(0.3, 0.5), add = TRUE, color=c("yellow", "red"))

2 Longitudinal Data (Timeseries/Kimeseries) Example

This timeseries demonstration shows the effects of indexing timeseries (univariate) data only using time and compares the representation of timeseries and kimeseries, which has profound impact on the subsequent data analytics. TO keep this application grounded, we will use real 4D fMRI data (\(x=64 \times y=64\times z=21\times t=180\)), but only focus on one spatial location (\({\bf{x}}=(x,y,z)\in \mathbb{R}^3\)). More details are provided in DSPA Chapter 3.

# install EBImage
# source("https://bioconductor.org/biocLite.R")
# biocLite("EBImage")
library(EBImage)
require(brainR)
library(spatstat) 

# 1D timeseries FFT SHIFT
fftshift1D <- function(img_ff) {
  rows <- length(img_ff)   
  rows_half <- ceiling(rows/2)
  return(append(img_ff[(rows_half+1):rows], img_ff[1:rows_half]))
}

# 1. download the 4D fMRI data
fMRIURL <- "http://socr.umich.edu/HTML5/BrainViewer/data/fMRI_FilteredData_4D.nii.gz"
fMRIFile <- file.path(tempdir(), "fMRI_FilteredData_4D.nii.gz")
download.file(fMRIURL, dest=fMRIFile, quiet=TRUE)
(fMRIVolume <- readNIfTI(fMRIFile, reorient=FALSE))
## NIfTI-1 format
##   Type            : nifti
##   Data Type       : 4 (INT16)
##   Bits per Pixel  : 16
##   Slice Code      : 0 (Unknown)
##   Intent Code     : 0 (None)
##   Qform Code      : 1 (Scanner_Anat)
##   Sform Code      : 0 (Unknown)
##   Dimension       : 64 x 64 x 21 x 180
##   Pixel Dimension : 4 x 4 x 6 x 3
##   Voxel Units     : mm
##   Time Units      : sec
# dimensions: 64 x 64 x 21 x 180 ; 4mm x 4mm x 6mm x 3 sec 

fMRIVolDims <- dim(fMRIVolume); fMRIVolDims
## [1]  64  64  21 180
# time_dim <- fMRIVolDims[4]; time_dim ## 180

# 2. extract the time-course of one voxel (25, 25, 12)  # 64  64 21 180
x1 <- c(1:180)
y1 <- loess(fMRIVolume[25, 25, 12, ]~ x1, family = "gaussian")
windows(width=10, height=8)
plot(fMRIVolume[25, 25, 12, ], type='l', 
     main="Time Series of 3D Voxel \n (x=25, y=25, z=12)", col="blue", 
     xlab = "Time", ylab = "fMRIVolume[25, 25, 12, ] Intensities")
lines(x1, smooth(fMRIVolume[25, 25, 12, ]), col = "red", lwd = 2)
lines(ksmooth(x1, fMRIVolume[25, 25, 12, ], kernel = "normal", bandwidth = 5), col = "green", lwd = 3)
legend("bottomright", legend=c("(raw) fMRI", "smooth(fMRI)", "ksmooth(fMRI"),
       col=c("blue", "red", "green"), lty=1, lwd=4, cex=1.1, y.intersp=0.6,
       x.intersp=1.0, title = "Voxel (25, 25, 12)", bty = "n")

# 3. FT of 1D time-course
y2 <- fMRIVolume[25, 25, 12, ]
X2 = fft(y2); plot(fftshift1D(log(Re(X2)+2)), main = "log(fftshift1D(Re(FFT(timeseries))))") 

X2_mag <- sqrt(Re(X2)^2+Im(X2)^2); plot(log(fftshift1D(Re(X2_mag))), main = "log(Magnitude(FFT(timeseries)))") 

X2_phase <- atan2(Im(X2), Re(X2)); plot(fftshift1D(X2_phase), main = "Shift(Phase(FFT(timeseries)))")

##  Implicit Automated IFT
hat_X2 = Re(fft(X2, inverse = T)/length(X2))
plot(hat_X2, main = "Re(IFT(FFT(timeseries)))")   # point plot represents the IFT(FT(data))
lines(x1, y2, col = "red", lwd = 2)               # lines represent the raw data time-course (should coincide)

## Manually invert the FT (IFT) using the magnitudes and phases
Real2 = X2_mag * cos(X2_phase)
Imaginary2 = X2_mag * sin(X2_phase)
man_hat_X2 = Re(fft(Real2 + 1i*Imaginary2, inverse = T)/length(X2))
plot(man_hat_X2, type="l", lwd=2, main = "Manual IFT (Magnitude, Phase) Synthesis")
lines(x1, y2, col = "red", lty=2, lwd = 2)  # overlay original data to confirm PERFECT manual reconstruction

# IFT reconstruction (synthesis) using fMRI-Magnitude and Nil-Phase
Real_phase0 = X2_mag * cos(0)
Imaginary_phase0 = X2_mag * sin(0)
ift_NilPhase_X2mag = Re(fft(Real_phase0 + 1i*Imaginary_phase0, inverse = T)/length(X2))
windows(width=10, height=8)
plot(ift_NilPhase_X2mag, col="red", type="l", lty=2, lwd=2, 
     ylim=c(10270, 10500), main=sprintf('Signal Synthesis: IFT (Magnitude=Real, Phase=Nil)\n Correlation(Real, Recon) = %s', format(cor(ift_NilPhase_X2mag, y2), digits=3)), xlab = "Time", 
     ylab = "IFT(NilPhase_Reconstruction)[25, 25, 12, ]", cex=1.3)
lines(x1, y2, col = "black", lwd = 2)  # overlay original data to show Nil-Phase effects on reconstruction
legend("top", bty="n", legend=c("(raw) fMRI", "Nil-Phase (Time-only) Reconstruction"),
       col=c("black", "red"), lty=c(1,2), lwd=c(2,2), cex=1.3, x.intersp=0.8)

# cor(ift_NilPhase_X2mag, y2)

################## Plot_ly Figure 6.5
time <- c(1:length(ift_NilPhase_X2mag))
p <- plot_ly()
p <- add_lines(p, x=~time,  y=~y2, 
               name = "f=fMRIVolume[25, 25, 12, ]", type = 'scatter', mode = 'lines', 
               hoverinfo = 'name', line=list(color='blue'))
p <- add_lines(p, x=~time, y=~ift_NilPhase_X2mag, 
               name = "f2=ift_NilPhase_X2mag[25, 25, 12, ]", type = 'scatter', mode = 'lines', 
               hoverinfo = 'name', 
               line=list(color='green', width = 8)) %>%
  layout(yaxis = list(range = c(10200,10500), title="Intensity"),
    title=sprintf('Signal Synthesis: IFT (Magnitude=Real, Phase=Nil)\n Correlation(Real, Recon) = %s',
                       format(cor(ift_NilPhase_X2mag, y2), digits=3)),
    legend = list(orientation = "h",   # show entries horizontally
                     xanchor = "center",  # use center of legend as anchor
                     x = 0.5, 
                     size = 30))      
p
# IFT reconstruction (synthesis) using fMRI-Magnitude and Random-Phase
set.seed(1234)
rand_Phase <- runif(length(X2_mag), -pi, pi)
#  rand_Phase <- runif(1, -pi, pi)
Real_phaseRand = X2_mag * cos(rand_Phase)
Imaginary_phaseRand = X2_mag * sin(rand_Phase)
ift_RandPhase_X2mag = Re(fft(Real_phaseRand + 1i*Imaginary_phaseRand, inverse = T)/length(X2))
plot(ift_RandPhase_X2mag-mean(ift_RandPhase_X2mag) + mean(y2), col="red", type="l", lty=2, lwd=2, 
     ylim=c(10300, 10500), main=sprintf('Signal Synthesis: IFT (Magnitude=Real, Phase=Random)\n Correlation(Real, Recon) = %s', format(cor(ift_RandPhase_X2mag, y2), digits=3)))
lines(x1, y2, col = "black", lwd = 2)  # overlay original data to show Rand-Phase effects on reconstruction
legend("top", bty="n", legend=c("(raw) fMRI", "Random-Phase Reconstruction"),
       col=c("black", "red"), lty=c(1,2), lwd=c(2,2), cex=0.9, x.intersp=0.5)

# IFT reconstruction (synthesis) using fMRI-Magnitude and Phase from a different voxel location (28, 22, 10)
y3 <- fMRIVolume[28, 22, 10, ]
X3 = fft(y3)  # ;  plot(fftshift1D(log(Re(X3)+2)), main = "log(fftshift1D(Re(FFT(timeseries))))") 
# X3_mag <- sqrt(Re(X3)^2+Im(X3)^2); plot(log(fftshift1D(Re(X3_mag))), main = "log(Magnitude(FFT(timeseries)))") 
neighbor_Phase <- atan2(Im(X3), Re(X3)) # ; plot(fftshift1D(X3_phase), main = "Shift(Phase(FFT(timeseries)))")

Real_phaseNeighbor = X2_mag * cos(neighbor_Phase)
Imaginary_phaseNeighbor = X2_mag * sin(neighbor_Phase)
ift_NeighborPhase_X2mag = Re(fft(Real_phaseNeighbor + 1i*Imaginary_phaseNeighbor, inverse = T)/length(X2))
plot(ift_NeighborPhase_X2mag, col="red", type="l", lty=2, lwd=2, 
     ylim=c(10300, 10500), main=sprintf('Signal Synthesis: IFT (Magnitude=Real(25, 25, 12), Phase=Voxel Neighbor(28, 22, 10))\n Correlation(Real, Recon) = %s', format(cor(ift_NeighborPhase_X2mag, y2), digits=3)))
lines(x1, y2, col = "black", lwd = 2)  # overlay original data to show Rand-Phase effects on reconstruction
legend("top", bty="n", legend=c("(raw) fMRI", "Neighbor-Derived-Phase Reconstruction"),
       col=c("black", "red"), lty=c(1,2), lwd=c(2,2), cex=0.9, x.intersp=0.5)

# IFT reconstruction (synthesis) using fMRI-Magnitude and Phase from a highly correlated voxel location
set.seed(1234)
y4 <- y2 + rnorm(n=length(y2), 0, 40) #; plot(y2, y4); cor(y2, y4)
X4 = fft(y4)  # ;  plot(fftshift1D(log(Re(X3)+2)), main = "log(fftshift1D(Re(FFT(timeseries))))") 
# X4_mag <- sqrt(Re(X4)^2+Im(X4)^2); plot(log(fftshift1D(Re(X4_mag))), main = "log(Magnitude(FFT(timeseries)))") 
corr_Phase <- atan2(Im(X4), Re(X4)) # ; plot(fftshift1D(X4_phase), main = "Shift(Phase(FFT(timeseries)))")

Real_phaseCorr = X2_mag * cos(corr_Phase)
Imaginary_phaseCorr = X2_mag * sin(corr_Phase)
ift_CorrPhase_X2mag = Re(fft(Real_phaseCorr + 1i*Imaginary_phaseCorr, inverse = T)/length(X2))
windows(width=10, height=8)
plot(ift_CorrPhase_X2mag, col="red", type="l", lty=2, lwd=2, 
     ylim=c(10250, 10510), main=sprintf('Signal Synthesis: IFT (Magnitude=Real, Phase=Highly-Correlated (%s) Voxel\n Correlation(Real, Recon) = %s', format(cor(y4, y2), digits=3),
                                        format(cor(ift_CorrPhase_X2mag, y2), digits=3)),
     xlab = "Time", ylab = "IFT(Highly-Corr_Reconstruction)[25, 25, 12, ] Intensities")
lines(x1, y2, col = "black", lwd = 2)  # overlay original data to show Rand-Phase effects on reconstruction
legend("top", bty="n", legend=c("(raw) fMRI", "Correlated-Voxel-Phase Reconstruction"),
       col=c("black", "red"), lty=c(1,2), lwd=c(2,2), cex=1.2, x.intersp=0.8)

##### Plot_Ly Figure 6.6
p <- plot_ly()
p <- add_lines(p, x=~time,  y=~y2, 
               name = "f=fMRIVolume[25, 25, 12, ]", type = 'scatter', mode = 'lines', 
               hoverinfo = 'name', line=list(color='blue'))
p <- add_lines(p, x=~time, y=~ift_CorrPhase_X2mag, 
               name = "f2=IFT(Highly-Corr_Reconstruction)[25, 25, 12, ] Intensities", type = 'scatter', mode = 'lines', 
               hoverinfo = 'name', 
               line=list(color='green', width = 8)) %>%
  layout(yaxis = list(range = c(10200,10500), title="Intensity"),
    title=sprintf('Signal Synthesis: IFT (Magnitude=Real, Phase=Highly-Correlated (%s)) Voxel\n Correlation(Real, Recon) = %s', 
                  format(cor(y4, y2), digits=3),
                  format(cor(ift_CorrPhase_X2mag, y2), digits=3)),
    legend = list(orientation = "h",   # show entries horizontally
                     xanchor = "center",  # use center of legend as anchor
                     x = 0.5, 
                     size = 30))      
p
# This is an effective reconstruction (synthesis) of the raw data by approximating the "unknown" phases.
# The reason why this works is that the real and estimated phases closely resemble each other
plot(X2_phase, corr_Phase, xlab = "Real Phase", ylab = "Approximate Phase", 
     main =sprintf('Scatterplot of True (x-axis) and Approximate Phases\n
     Correlation = %s', format(cor(X2_phase, corr_Phase), digits=3)))

# Similarly, the reconstructed and real signal are highly correlated as the scatterplot shows
plot(y2, y4, xlab = "Real Signal", ylab = "Reconstructed (using estimated Phases)", 
     main =sprintf('Scatterplot of True Signal (x-axis) and Approximate-Phase reconstructed Signal\n
     Correlation = %s', format(cor(y2, y4), digits=3)))

These examples demonstrate the timeseries representation and analysis work well in spacetime. However, in various situations where one may or may not be able to observe or estimate the kime-direction (phase) the results can widely vary based on how reasonable to synthesis of information is without explicit knowledge of the phase measures. As an observable, the time (kime-order) is measurable and the phase angles (kime-direction) can either be estimated from other similar data, provided by an oracle, or fixed according to some experimental conditions.

SOCR Resource Visitor number Web Analytics SOCR Email