See all details on UMich Biophysics 430 Winter 2020 Canvas website.

Problem 1

Let’s use the 3D sMRI data for this problem demonstation.

library(brainR)
## Loading required package: rgl
## Loading required package: misc3d
## Loading required package: oro.nifti
## oro.nifti 0.9.1
# Test data: http://socr.umich.edu/HTML5/BrainViewer/data/TestBrain.nii.gz

brainURL <- "http://socr.umich.edu/HTML5/BrainViewer/data/TestBrain.nii.gz"
brainFile <- file.path(tempdir(), "TestBrain.nii.gz")
download.file(brainURL, dest=brainFile, quiet=TRUE)
brainVolume <- readNIfTI(brainFile, reorient=FALSE)

brainVolDims <- dim(brainVolume); brainVolDims
## [1] 181 217 181
# Show graphs, plots, visualizations of each volume in 1D (linear section), 2D (planar cross-sections), and 3D (volume rendering).
plot(brainVolume[90,108,], type='l', main="1D plot of sMRI intensities \n (x=90, y=108, z)", col="blue", xlab = "1<=z<=181", ylab = "sMRI Intensity") # 1D

image(brainVolume[,, 90], asp=1, col = hcl.colors(12, "YlOrRd", rev = TRUE)) # 2D

image(brainVolume[,, 90], asp=1, col = rainbow(10, start = 0, end = 1, alpha = 1, rev = FALSE))

image(brainVolume[,, 90], asp=1, col = heat.colors(10, alpha = 1, rev = FALSE))

orthographic(brainVolume, xyz=c(90,108,90), zlim=range(brainVolume)*0.9) # 3D

# Show the effects of thresholding, gamma correction through exponentiation, and transposition and rotation of 2D images.
library(jpeg)
library(EBImage)

sMRI_2D <- as.matrix(brainVolume[,, 90])
sMRI_2D <- (sMRI_2D - min(sMRI_2D))/diff(range(sMRI_2D))
display( thresh(sMRI_2D, offset=0.05), all=TRUE )

# Sheer Affine transformation
m = matrix(c(1, -.5, 128, 0, 1, 0), nrow=3, ncol=2)
img_affine = affine(translate(sMRI_2D, c(-50, 0)), m)
display( img_affine )

# Gaussian Kernel smoothing (Convolution)
kernel = makeBrush(size = 31, shape = 'gaussian', sigma = 5)
mri_smooth = filter2(sMRI_2D, kernel)
display(mri_smooth)

# Generate histograms of the 3D volumes and some 2D image sections.
hist(brainVolume) # 3D

hist(brainVolume[,, 90]) # 2D slice

hist(brainVolume[, 107, 90]) # 1D line

# Save some of the graphs in PNG, TIFF, or JPG rasterized image files.
png(filename="sMRI.png"); display(mri_smooth); dev.off()
## png 
##   2

Modify some of hte above examples for the fMRI and PET data.

Problem 2

Let’s now work with the fMRI (4D) volume for this propblem.

# load the PET data
fMRIURL <- "http://socr.umich.edu/HTML5/BrainViewer/data/fMRI_FilteredData_4D.nii.gz"
fMRIFile <- file.path(tempdir(), "fMRI.nii.gz")
download.file(fMRIURL, dest=fMRIFile, quiet=TRUE)
fMRIVolume <- readNIfTI(fMRIFile, reorient=FALSE)
# remove small background noise
fMRIVol <- ifelse(fMRIVolume > 8000, fMRIVolume/1000, 0)
# dim(fMRIVolume)   # 64  64  21 180
# image(fMRIVol[,, 14, 90], col = rainbow(10, start = 0, end = 1, alpha = 1, rev = FALSE))
display(colorLabels(fMRIVol[,, 14, 90]), method="raster")

# Smooth (blur, low-pass filter), denoise (median filter), and sharpen (contour detect, Laplacian filter) the images.
fMRI_gblur = gblur(fMRIVol, sigma = 2)
# mask the blured volume (remove outside brain smoothing artifacts)
fMRI_gblur_mask = ifelse(fMRIVol>0, fMRI_gblur, 0)
display(colorLabels(fMRI_gblur_mask[,, 14, 90]), method="raster")

fMRI_Denoise_median = medianFilter(fMRIVol[,, 14, 90], 2)
hist(fMRI_Denoise_median)

display(colorLabels(fMRI_Denoise_median), method="raster")

# Apply watershed image segmentation.
fMRI_Watershed_seg = watershed( fMRIVol[,, 14, 90], ext=2 )
display(normalize(fMRI_Watershed_seg), method="raster")

# Try Voronoi tesselation.
voronoi_fMRI = propagate(seeds = fMRI_Watershed_seg, x = fMRI_Watershed_seg, lambda = 100)
voronoiPaint = colorLabels (voronoi_fMRI)
display(voronoiPaint, method="raster")

# Display the Raw image , Fourier Transform of the image, filter the FT(image), e.g., threshold, low or high pass filter, and invert the processed image back in spacetime.
###### Some houskeeping functions
# 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")
  }
}

### FT
img <- fMRIVol[,, 14, 90]
imgMat <- matrix(img, nrow = dim(img)[1], ncol = dim(img)[2])
display(normalize(imgMat), method="raster")

ft_imgMat <- fft(imgMat)  # fftw2d # Display Re(FT):

# Display FOurier Magnitude and phase
mag_ft_imgMat <- sqrt(Re(ft_imgMat)^2+Im(ft_imgMat)^2)
display(normalize(fftshift(mag_ft_imgMat)), method="raster")

# Phase  <- atan(Im(img_ff)/Re(img_ff))
phase_ft_imgMat  <- atan2(Im(ft_imgMat), Re(ft_imgMat))
display(normalize(fftshift(phase_ft_imgMat)), method="raster")

## Low-pass frequency filtering
# Remove Fourier coefficients indexed above $k_{x,max/2}$ and $k_{y,max/2}$.
mag_ft_imgMat_Low <- mag_ft_imgMat
for (i in 1:dim(mag_ft_imgMat)[1]) {
  for (j in 1:dim(mag_ft_imgMat)[2]) {
    if (abs(i-(dim(mag_ft_imgMat)[1])/2) > (dim(mag_ft_imgMat)[1])/4 | 
        abs(j-(dim(mag_ft_imgMat)[2])/2) > (dim(mag_ft_imgMat)[2])/4) {
      mag_ft_imgMat_Low[i,j] <- 0
    }
  }
}

EBImage::display(normalize(mag_ft_imgMat_Low), method = "raster", title="(IFT o LowPassFilter o FT) Magnitude=Img | Phase=Img")

Real = mag_ft_imgMat_Low * cos(phase_ft_imgMat)
Imaginary = mag_ft_imgMat_Low * sin(phase_ft_imgMat)
ift_ft_imgMat_Low = Re(fft(Real+1i*Imaginary, inverse = T)/length(mag_ft_imgMat_Low))
EBImage::display(normalize(ift_ft_imgMat_Low), method = "raster", title="(IFT o LowPassFilter o FT) Magnitude=Img | Phase=Img")

display(normalize(imgMat), method="raster") # Compare to original fMRI

### High-pass frequency filtering
# Remove the Fourier coefficients indexed below $k_{x,max/4}$ and $k_{y,max/4}$.
# High-pass freqency filtering, remove Fourier coeff's indexed below k_{x,max/4} and k_{y,max/4}
mag_ft_imgMat_High <- mag_ft_imgMat
for (i in 1:dim(mag_ft_imgMat)[1]) {
  for (j in 1:dim(mag_ft_imgMat)[2]) {
    if (abs(i-(dim(mag_ft_imgMat)[1])/2) < (dim(mag_ft_imgMat)[1])/4 & 
        abs(j-(dim(mag_ft_imgMat)[2])/2) < (dim(mag_ft_imgMat)[2])/4) {
      mag_ft_imgMat_High[i,j] <- 0
    }
  }
}

EBImage::display(mag_ft_imgMat_High, method = "raster", title="(IFT o HighPassFilter o FT) Magnitude=Img | Phase=Img")

Real = mag_ft_imgMat_High * cos(phase_ft_imgMat)
Imaginary = mag_ft_imgMat_High * sin(phase_ft_imgMat)
ift_ft_imgMat_High = Re(fft(Real+1i*Imaginary, inverse = T)/length(mag_ft_imgMat_High))
EBImage::display(normalize(ift_ft_imgMat_High), method = "raster", title="(IFT o HighPassFilter o FT) Magnitude=Img | Phase=Img")

# Show empirically some of these FT properties using 2D images cross-sections:Linearity property; Multiplication property; removing either the central (removing low-frequency Fourier components) or the peripheral (remove high-frequency Fouriercomponents).

# Let's just illustrate the Linearity Property of the FT for 2D images
image1 <- normalize(fMRIVol[,, 14, 90])
display(image1, method="raster")

image2 <- array(runif(n=dim(image1)[1]*dim(image1)[2], 0, 1), dim=c(dim(image1)[1], dim(image1)[2]))
display(image1, method="raster") # dim(image2)
sum_img <- 2*image1 + 0.5*image2
display(normalize(sum_img), method="raster") 

ft_sum_img <- fft(imgMat)  # fftw2d # Display Re(FT):
# Display Fourier Magnitude 
mag_ft_imgMat <- sqrt(Re(ft_sum_img)^2+Im(ft_sum_img)^2)
display(normalize(fftshift(mag_ft_imgMat)), method="raster")

ft_sum_img1 <- fft(image1)  # fftw2d # Display Re(FT):
mag_ft_img1 <- sqrt(Re(ft_sum_img1)^2+Im(ft_sum_img1)^2)
ft_sum_img2 <- fft(image2)  # fftw2d # Display Re(FT):
mag_ft_img2 <- sqrt(Re(ft_sum_img2)^2+Im(ft_sum_img2)^2)
mag_ft_sum_img <- 2*mag_ft_img1 + 0.5*mag_ft_img2
# Display Fourier Magnitude 
display(normalize(fftshift(mag_ft_sum_img)), method="raster")

# Diff between the 2 should be noise
display((fftshift(mag_ft_sum_img-mag_ft_imgMat)), method="raster")

Problem 3

We’ll demonstrate linear modeling using the fMRI data.

# Extract the (160-point) fMRI time-series at a couple of (x,y,z) spatial locations (Links to an external site.).

time_series <- img <- fMRIVol[33,33, 14, ]; length(time_series) # 180
## [1] 180
# Fit ARIMA models on training data 1 ≤ t i m e ≤ 140 and assess the model predictions against the testing data 141 ≤ t i m e ≤ 160
library(forecast); library(TTR)
## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## Registered S3 methods overwritten by 'forecast':
##   method             from    
##   fitted.fracdiff    fracdiff
##   residuals.fracdiff fracdiff
ts <- ts(time_series, start=1, end=180, frequency=1)
ts.plot(ts)

ts.epoch <- SMA(ts, n=10); plot(ts.epoch)

acf(ts, lag.max = 20, main="ACF")

pacf(ts, lag.max = 20, main="PACF")

# Fit ARIMA model
fit<-auto.arima(ts.epoch, approx=F, trace = F)
fit
## Series: ts.epoch 
## ARIMA(3,0,3) with non-zero mean 
## 
## Coefficients:
##          ar1     ar2      ar3      ma1      ma2     ma3     mean
##       0.9857  0.5741  -0.6956  -0.0732  -0.7285  0.3087  14.9293
## s.e.  0.0965  0.1431   0.0756   0.1330   0.0765  0.1065   0.0076
## 
## sigma^2 estimated as 0.0007303:  log likelihood=376.74
## AIC=-737.47   AICc=-736.58   BIC=-712.34
# Forecast using ARIMA model
ts.forecasts<-forecast(fit, h=20)
plot(ts.forecasts, include = 400)

# LM: design-matrix includes 2 columns (intercept and an alternating ON vs. OFF visual input where each ON or OFF epoch continues for 10 time-points). Thus, there are 8 ON intertwined with  8 OFF stimulus conditions
# Fit a linear model Y = X β + ϵ and report the coefficient (parameter) estimates
stim <- rep(cbind(rep(1,10), rep(-1,10)), 9); stim; length(stim)
##   [1]  1  1  1  1  1  1  1  1  1  1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  1  1  1  1  1
##  [26]  1  1  1  1  1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  1  1  1  1  1  1  1  1  1  1
##  [51] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  1  1  1  1  1  1  1  1  1  1 -1 -1 -1 -1 -1
##  [76] -1 -1 -1 -1 -1  1  1  1  1  1  1  1  1  1  1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
## [101]  1  1  1  1  1  1  1  1  1  1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  1  1  1  1  1
## [126]  1  1  1  1  1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  1  1  1  1  1  1  1  1  1  1
## [151] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  1  1  1  1  1  1  1  1  1  1 -1 -1 -1 -1 -1
## [176] -1 -1 -1 -1 -1
## [1] 180
X <- as.matrix(cbind(rep(1, 160),stim)); dim(X)
## Warning in cbind(rep(1, 160), stim): number of rows of result is not a multiple
## of vector length (arg 1)
## [1] 180   2
Y <- ts
lm_fit <- lm(Y ~ stim); summary(lm_fit)
## 
## Call:
## lm(formula = Y ~ stim)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.44578 -0.14403 -0.01847  0.13018  0.59083 
## 
## Coefficients:
##             Estimate Std. Error  t value Pr(>|t|)    
## (Intercept) 14.92897    0.01486 1004.737   <2e-16 ***
## stim         0.01481    0.01486    0.996     0.32    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1993 on 178 degrees of freedom
## Multiple R-squared:  0.005547,   Adjusted R-squared:  -3.98e-05 
## F-statistic: 0.9929 on 1 and 178 DF,  p-value: 0.3204