Biophysics430
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
## 2Modify 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 fracdiffts <- 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] 180X <- 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 2Y <- 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