SOCR ≫ | TCIU Website ≫ | TCIU GitHub ≫ |
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 R^3\)) to 2D (\({\bf{x}}=(x,y)\in R^2\subseteq R^3\)) by focusing only on mid-axial/transverse slice through the brain (\(z=11\)). More details are provided in DSPA Chapter 3.
# 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")
## 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):R^2\times 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:
## [1] 64 64 180
# 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
##
## Attaching package: 'abind'
## The following object is masked from 'package:EBImage':
##
## abind
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"))
# Save the output image
rgl.snapshot(filename = "Figures/Fig3.11B.png")
# FFT SHIFT
#' This function is useful for visualizing the Fourier transform with the zero-frequency
#' component in the middle of the spectrum.
#'
#' @param img_ff A Fourier transform of a 1D signal, 2D image, or 3D volume.
#' @param dim Number of dimensions (-1, 1, 2, 3).
#' @return A properly shifted FT of the array.
#'
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")
}
}
# IFFT SHIFT
#' This function is useful for moving back the zero-frequency component in the middle of the spectrum
#' back to (0,0,0). It rearranges in reverse (relative to fftshift()) the indices appropriately,
#' so that the image can be correctly reconstructed by the IFT in spacetime
#'
#' @param img_ff An Inverse Fourier transform of a 1D signal, 2D image, or 3D volume.
#' @param dim Number of dimensions (-1, 1, 2, 3).
#' @return A properly shifted IFT of the input array.
#'
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")
}
}