library(oro.nifti)
library(plotly)
library(OpenImageR)

# ==== LOAD MRI DATA ====
mriFile <- "ABIDE_MRI_MPRAGE_peds_defaced.nii.gz"
mriVol <- readNIfTI(mriFile, reorient = FALSE)

dims <- dim(mriVol)
cat("Dimensions:", dims, "\n")
## Dimensions: 256 182 256
# ==== SELECT CENTRAL SLICES ====
x_mid <- round(dims[1]/2)
y_mid <- round(dims[2]/2)
z_mid <- round(dims[3]/2)

slice_xy <- mriVol[,,z_mid]
slice_xz <- mriVol[,y_mid,]
slice_yz <- mriVol[x_mid,,]

# ==== NORMALIZATION FUNCTION ====
normalize <- function(img) {
  img <- img - min(img)
  img <- img / max(img)
  return(img)
}

# ==== 1D LINEAR PROFILE ====
plot(mriVol[x_mid, y_mid, ],
     type = "l",
     main = "1D Intensity Profile (Z-direction)",
     xlab = "Slice Index",
     ylab = "Intensity")

# ==== 2D VISUALIZATIONS ====
par(mfrow=c(1,3))
image(t(apply(slice_xy, 2, rev)), col=gray(0:64/64), main="XY Slice")
image(t(apply(slice_xz, 2, rev)), col=gray(0:64/64), main="XZ Slice")
image(t(apply(slice_yz, 2, rev)), col=gray(0:64/64), main="YZ Slice")

# ==== COLOR MAP EXAMPLES ====
par(mfrow=c(1,3))
image(t(apply(slice_xy,2,rev)), col=heat.colors(64), main="Heat Map")
image(t(apply(slice_xy,2,rev)), col=terrain.colors(64), main="Terrain Map")
image(t(apply(slice_xy,2,rev)), col=topo.colors(64), main="Topo Map")

# ==== GAMMA CORRECTION ====
gamma_img <- slice_xy^0.5
image(t(apply(gamma_img,2,rev)), col=gray(0:64/64), main="Gamma Corrected")

# ==== THRESHOLDING (FIXED) ====
thresh <- (slice_xy > mean(slice_xy)) * 1
image(t(apply(thresh,2,rev)), col=gray(0:1), main="Thresholded")

# ==== ROTATION / TRANSPOSE ====
rot_img <- t(slice_xy)
image(t(apply(rot_img,2,rev)), col=gray(0:64/64), main="Rotated")

# ==== HISTOGRAMS ====
par(mfrow=c(1,2))
hist(mriVol, breaks=50, main="3D Volume Histogram", col="blue")
hist(slice_xy, breaks=50, main="2D Slice Histogram", col="green")

# ==== 3D VOLUME RENDERING ====
ds <- 4

mri_small <- mriVol[seq(1,dims[1],by=ds),
                    seq(1,dims[2],by=ds),
                    seq(1,dims[3],by=ds)]

dims_small <- dim(mri_small)

vals <- as.vector(mri_small)

plot_ly(
  x = rep(1:dims_small[1], each=dims_small[2]*dims_small[3]),
  y = rep(rep(1:dims_small[2], each=dims_small[3]), dims_small[1]),
  z = rep(1:dims_small[3], dims_small[1]*dims_small[2]),
  value = vals,
  type = "volume",
  opacity = 0.5,
  surface_count = 20,
  isomin = quantile(vals, 0.7),   # 👈 filters low-intensity noise
  isomax = max(vals),
  colorscale = "Hot"
)
# ==== SAVE PROCESSED IMAGES ====
writeImage(normalize(slice_xy), "slice_xy.png")
writeImage(normalize(gamma_img), "gamma.png")
writeImage(normalize(thresh), "threshold.png")
writeImage(normalize(rot_img), "rotated.png")
writeImage(normalize(slice_xz), "slice_xz.png")
writeImage(normalize(slice_yz), "slice_yz.png")

# ==== DISPLAY SAVED IMAGES ====
par(mfrow=c(2,3))
image(normalize(slice_xy), col=gray(0:64/64), main="Saved XY")
image(normalize(gamma_img), col=gray(0:64/64), main="Saved Gamma")
image(normalize(thresh), col=gray(0:1), main="Saved Threshold")
image(normalize(rot_img), col=gray(0:64/64), main="Saved Rotated")
image(normalize(slice_xz), col=gray(0:64/64), main="Saved XZ")
image(normalize(slice_yz), col=gray(0:64/64), main="Saved YZ")

#Summary

“The MRI dataset was analyzed using 1D, 2D, and 3D visualization techniques to examine intensity variation and structure. A 1D intensity profile along the z-direction at the center of the volume shows how voxel intensity changes across slices and highlights regions of higher signal.

Central slices in the XY, XZ, and YZ planes were used to visualize the internal structure of the brain. Different color maps, including grayscale and heat-based schemes, were applied to enhance contrast and make intensity differences more visible.

Basic image processing techniques were also demonstrated. Gamma correction adjusted the intensity distribution to improve contrast, while thresholding created a binary image that separates higher-intensity regions from background noise. A rotation operation was applied to illustrate geometric transformation of the image.

Histograms of both the full 3D volume and a 2D slice were generated to show the distribution of intensity values. Finally, a downsampled 3D volume rendering was created to visualize the spatial structure of the MRI data while maintaining reasonable computational cost.

Overall, these methods provide complementary ways to interpret MRI data, from simple intensity trends to full spatial visualization.”

#actual problem 2

library(oro.nifti)
library(EBImage)
library(plotly)

# load the PET data
petFile <- "PET_FDG_3D_vol.nii"
petVolume <- readNIfTI(petFile, reorient=FALSE)
petVol <- petVolume@.Data
petVol <- petVol * 0.00392156885936856
petVol[is.na(petVol)] <- 0
# dim(petVol)   # 91 109 91

slice_index <- round(dim(petVol)[3]/2)
img <- petVol[,, slice_index]
img <- ifelse(img > 0, img, 0)
img <- (img - min(img))/(max(img) - min(img))
imgMat <- matrix(img, nrow=dim(img)[1], ncol=dim(img)[2])
imgEB <- Image(imgMat, colormode="Grayscale")

plot_ly(z=~t(imgMat), type="heatmap", showscale=FALSE) %>%
  layout(title=paste0("2D FDG-PET Image [x, y, ", slice_index, "]"), xaxis=list(title="x"), yaxis=list(title="y"))
plot_ly(z=t(imgMat), type="heatmap") %>%
  layout(title=paste0("2D FDG-PET Image [x, y, ", slice_index, "]"))
# Smooth (blur, low-pass filter), denoise (median filter), and sharpen (contour detect, Laplacian filter) the images.
pet_gblur <- gblur(imgEB, sigma=2)
plot_ly(z=~t(imageData(pet_gblur)), type="heatmap") %>%
  layout(title=paste0("Blurred 2D FDG-PET Image [x, y, ", slice_index, "]"))
pet_median <- medianFilter(imgEB, size=3)
plot_ly(z=~t(imageData(pet_median)), type="heatmap") %>%
  layout(title=paste0("Median Filtered 2D FDG-PET Image [x, y, ", slice_index, "]"))
laplacian_kernel <- matrix(c(0,-1,0,-1,4,-1,0,-1,0), nrow=3, byrow=TRUE)
pet_laplace <- filter2(imgEB, laplacian_kernel)
plot_ly(z=~t(imageData(normalize(pet_laplace))), type="heatmap") %>%
  layout(title=paste0("Laplacian Filter 2D FDG-PET Image [x, y, ", slice_index, "]"))
pet_sharp <- normalize(imgEB + pet_laplace)
plot_ly(z=~t(imageData(pet_sharp)), type="heatmap") %>%
  layout(title=paste0("Sharpened 2D FDG-PET Image [x, y, ", slice_index, "]"))
plot_ly(x=~as.vector(imgMat), type="histogram", histnorm="probability") %>%
  layout(title=paste0("Histogram of 2D FDG-PET Slice [x, y, ", slice_index, "]"), xaxis=list(title="value"), yaxis=list(title="freq"))
pet_th <- gblur(imgEB, sigma=2)
pet_th <- pet_th > otsu(pet_th)
plot_ly(z=~t(imageData(pet_th)), type="heatmap") %>%
  layout(title=paste0("Thresholded 2D FDG-PET Image [x, y, ", slice_index, "]"))
## Warning: `textfont.color` does not currently support multiple values.
pet_th <- opening(pet_th, makeBrush(5, shape="disc"))
pet_th <- fillHull(pet_th)
plot_ly(z=~t(imageData(pet_th)), type="heatmap") %>%
  layout(title=paste0("Cleaned Thresholded 2D FDG-PET Image [x, y, ", slice_index, "]"))
## Warning: `textfont.color` does not currently support multiple values.
# Apply watershed image segmentation.
pet_distmap <- distmap(pet_th)
plot_ly(z=~t(imageData(normalize(pet_distmap))), type="heatmap") %>%
  layout(title=paste0("Distance Map of 2D FDG-PET Image [x, y, ", slice_index, "]"))
pet_Watershed_seg <- watershed(pet_distmap, tolerance=1, ext=1)
plot_ly(z=~t(imageData(normalize(pet_Watershed_seg))), type="heatmap") %>%
  layout(title=paste0("Watershed Segmentation of 2D FDG-PET Image [x, y, ", slice_index, "]")) %>%
  hide_colorbar()
display(colorLabels(pet_Watershed_seg), method="raster")

computeFeatures.shape(pet_Watershed_seg)
##   s.area s.perimeter s.radius.mean s.radius.sd s.radius.min s.radius.max
## 1   5244         234      40.43415    3.118187      35.3169     45.46232
computeFeatures.basic(pet_Watershed_seg, imgEB)
##      b.mean      b.sd    b.mad    b.q001    b.q005     b.q05   b.q095    b.q099
## 1 0.6394061 0.1508772 0.184556 0.3526971 0.4024896 0.6390041 0.879668 0.9294606
####### Try Voronoi tessellation.
voronoi_pet <- propagate(seeds=pet_Watershed_seg, x=imgEB, mask=pet_th, lambda=100)
voronoiPaint <- colorLabels(voronoi_pet)
display(voronoiPaint, method="raster")
plot_ly(z=~t(imageData(normalize(voronoi_pet))), type="heatmap") %>%
  layout(title=paste0("Voronoi Tessellation of 2D FDG-PET Image [x, y, ", slice_index, "]")) %>%
  hide_colorbar()
computeFeatures.shape(voronoi_pet)
##   s.area s.perimeter s.radius.mean s.radius.sd s.radius.min s.radius.max
## 1   5244         234      40.43415    3.118187      35.3169     45.46232
computeFeatures.basic(voronoi_pet, imgEB)
##      b.mean      b.sd    b.mad    b.q001    b.q005     b.q05   b.q095    b.q099
## 1 0.6394061 0.1508772 0.184556 0.3526971 0.4024896 0.6390041 0.879668 0.9294606
# 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 housekeeping functions
# FFT SHIFT
fftshift <- function(img_ff, dim=-1) {
  rows <- dim(img_ff)[1]
  cols <- dim(img_ff)[2]
  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]))
  }
  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 {
    stop("Invalid dimension parameter")
  }
}
ifftshift <- fftshift

### FT
img <- petVol[,, slice_index]
img <- ifelse(img > 0, img, 0)
img <- (img - min(img))/(max(img) - min(img))
imgMat <- matrix(img, nrow=dim(img)[1], ncol=dim(img)[2])
plot_ly(z=~t(imgMat), type="heatmap") %>%
  layout(title=paste0("2D FDG-PET Image [x, y, ", slice_index, "]")) %>% hide_colorbar()
ft_imgMat <- fft(imgMat)
mag_ft_imgMat <- sqrt(Re(ft_imgMat)^2+Im(ft_imgMat)^2)
plot_ly(z=~t(fftshift(log(1+mag_ft_imgMat))), type="heatmap") %>%
  layout(title=paste0("(Magnitude) FT 2D FDG-PET Image [x, y, ", slice_index, "]")) %>% hide_colorbar()
phase_ft_imgMat <- atan2(Im(ft_imgMat), Re(ft_imgMat))
plot_ly(z=~t(fftshift(phase_ft_imgMat)), type="heatmap") %>%
  layout(title=paste0("(Phase) FT 2D FDG-PET Image [x, y, ", slice_index, "]")) %>% hide_colorbar()
ift_imgMat <- Re(fft(ft_imgMat, inverse=TRUE)/length(ft_imgMat))
plot_ly(z=~t(ift_imgMat), type="heatmap") %>%
  layout(title=paste0("(IFT o FT) 2D FDG-PET Image [x, y, ", slice_index, "]")) %>% hide_colorbar()
## Low-pass frequency filtering
# Remove Fourier coefficients indexed above kx,max/2 and ky,max/2.
ft_shift <- fftshift(ft_imgMat)
mag_shift <- Mod(ft_shift)
phase_shift <- atan2(Im(ft_shift), Re(ft_shift))
mag_ft_imgMat_Low <- mag_shift
for (i in 1:dim(mag_ft_imgMat_Low)[1]) {
  for (j in 1:dim(mag_ft_imgMat_Low)[2]) {
    if (abs(i-(dim(mag_ft_imgMat_Low)[1])/2) > (dim(mag_ft_imgMat_Low)[1])/4 | 
        abs(j-(dim(mag_ft_imgMat_Low)[2])/2) > (dim(mag_ft_imgMat_Low)[2])/4) {
      mag_ft_imgMat_Low[i,j] <- 0
    }
  }
}
plot_ly(z=~t(log(1+mag_ft_imgMat_Low)), type="heatmap") %>%
  layout(title=paste0("(LowPassFilter o FT) of 2D FDG-PET Image [x, y, ", slice_index, "]")) %>% hide_colorbar()
Real <- mag_ft_imgMat_Low * cos(phase_shift)
Imaginary <- mag_ft_imgMat_Low * sin(phase_shift)
ift_ft_imgMat_Low <- Re(fft(ifftshift(Real+1i*Imaginary), inverse=TRUE)/length(mag_ft_imgMat_Low))
plot_ly(z=~t(ift_ft_imgMat_Low), type="heatmap") %>%
  layout(title=paste0("(IFT o LowPassFilter o FT) of 2D FDG-PET Image [x, y, ", slice_index, "]")) %>% hide_colorbar()
### High-pass frequency filtering
# Remove the Fourier coefficients indexed below kx,max/4 and ky,max/4.
mag_ft_imgMat_High <- mag_shift
for (i in 1:dim(mag_ft_imgMat_High)[1]) {
  for (j in 1:dim(mag_ft_imgMat_High)[2]) {
    if (abs(i-(dim(mag_ft_imgMat_High)[1])/2) < (dim(mag_ft_imgMat_High)[1])/4 & 
        abs(j-(dim(mag_ft_imgMat_High)[2])/2) < (dim(mag_ft_imgMat_High)[2])/4) {
      mag_ft_imgMat_High[i,j] <- 0
    }
  }
}
plot_ly(z=~t(log(1+mag_ft_imgMat_High)), type="heatmap") %>%
  layout(title=paste0("(HighPassFilter o FT) of 2D FDG-PET Image [x, y, ", slice_index, "]")) %>% hide_colorbar()
Real <- mag_ft_imgMat_High * cos(phase_shift)
Imaginary <- mag_ft_imgMat_High * sin(phase_shift)
ift_ft_imgMat_High <- Re(fft(ifftshift(Real+1i*Imaginary), inverse=TRUE)/length(mag_ft_imgMat_High))
plot_ly(z=~t(ift_ft_imgMat_High), type="heatmap") %>%
  layout(title=paste0("(IFT o HighPassFilter o FT) of 2D FDG-PET Image [x, y, ", slice_index, "]")) %>% hide_colorbar()
# Remove central components below kx,(1/4)max and ky,(1/4)max.
nr <- nrow(ft_shift)
nc <- ncol(ft_shift)
center_r <- floor(nr/2)+1
center_c <- floor(nc/2)+1
kx_max <- floor(nc/2)
ky_max <- floor(nr/2)
row_dist <- abs(matrix(rep(1:nr, nc), nrow=nr) - center_r)
col_dist <- abs(matrix(rep(1:nc, each=nr), nrow=nr) - center_c)
kx_cut_1_4 <- kx_max/4
ky_cut_1_4 <- ky_max/4
remove_central_mask <- !(col_dist <= kx_cut_1_4 & row_dist <= ky_cut_1_4)
ft_remove_central <- ft_shift * remove_central_mask
img_remove_central <- Re(fft(ifftshift(ft_remove_central), inverse=TRUE)/length(ft_remove_central))
plot_ly(z=~t(remove_central_mask), type="heatmap") %>%
  layout(title="Mask: Remove Central Low Frequencies below 1/4 max") %>% hide_colorbar()
## Warning: `textfont.color` does not currently support multiple values.
plot_ly(z=~t(log(1+Mod(ft_remove_central))), type="heatmap") %>%
  layout(title="FT After Removing Central Low Frequencies below 1/4 max") %>% hide_colorbar()
plot_ly(z=~t(img_remove_central), type="heatmap") %>%
  layout(title="Image After Removing Central Low Frequencies below 1/4 max") %>% hide_colorbar()
# Remove peripheral components above kx,(1/2)max and ky,(1/2)max.
kx_cut_1_2 <- kx_max/2
ky_cut_1_2 <- ky_max/2
remove_peripheral_mask <- (col_dist <= kx_cut_1_2 & row_dist <= ky_cut_1_2)
ft_remove_peripheral <- ft_shift * remove_peripheral_mask
img_remove_peripheral <- Re(fft(ifftshift(ft_remove_peripheral), inverse=TRUE)/length(ft_remove_peripheral))
plot_ly(z=~t(remove_peripheral_mask), type="heatmap") %>%
  layout(title="Mask: Remove Peripheral High Frequencies above 1/2 max") %>% hide_colorbar()
## Warning: `textfont.color` does not currently support multiple values.
plot_ly(z=~t(log(1+Mod(ft_remove_peripheral))), type="heatmap") %>%
  layout(title="FT After Removing Peripheral High Frequencies above 1/2 max") %>% hide_colorbar()
plot_ly(z=~t(img_remove_peripheral), type="heatmap") %>%
  layout(title="Image After Removing Peripheral High Frequencies above 1/2 max") %>% hide_colorbar()
# Show empirically some of these FT properties using 2D images cross-sections: Linearity property; Multiplication property.

# Linearity property of the FT for 2D image slices.
image1 <- petVol[,, round(dim(petVol)[3]*0.40)]
image2 <- petVol[,, round(dim(petVol)[3]*0.60)]
image1[is.na(image1)] <- 0
image2[is.na(image2)] <- 0
image1 <- (image1 - min(image1))/(max(image1) - min(image1))
image2 <- (image2 - min(image2))/(max(image2) - min(image2))
a <- 2
b <- 3
sum_img <- a*image1 + b*image2
plot_ly(z=~t(sum_img), type="heatmap") %>%
  layout(title="Mixture a*PET_slice1 + b*PET_slice2") %>% hide_colorbar()
ft_sum_img <- fft(sum_img)
ft_img1 <- fft(image1)
ft_img2 <- fft(image2)
ft_linear_combo <- a*ft_img1 + b*ft_img2
linearity_error_2d <- max(Mod(ft_sum_img - ft_linear_combo))
linearity_error_2d
## [1] 3.638423e-12
plot_ly(z=~t(fftshift(log(1+Mod(ft_sum_img)))), type="heatmap") %>%
  layout(title="Magnitude FT(a*f + b*g)") %>% hide_colorbar()
plot_ly(z=~t(fftshift(log(1+Mod(ft_linear_combo)))), type="heatmap") %>%
  layout(title="Magnitude a*FT(f) + b*FT(g)") %>% hide_colorbar()
plot_ly(z=~t(fftshift(Mod(ft_sum_img - ft_linear_combo))), type="heatmap") %>%
  layout(title="Difference: FT(a*f+b*g) - [a*FT(f)+b*FT(g)]") %>% hide_colorbar()
# Multiplication property.
# FT of a convolution is a product of FTs.
f_img <- imgMat
kernel_2d <- makeBrush(size=9, shape="Gaussian", sigma=2)
kernel_2d <- kernel_2d/sum(kernel_2d)
kernel_padded <- matrix(0, nrow=nrow(f_img), ncol=ncol(f_img))
kernel_padded[1:nrow(kernel_2d), 1:ncol(kernel_2d)] <- kernel_2d
F_f <- fft(f_img)
F_g <- fft(kernel_padded)
product_FT <- F_f * F_g
conv_frequency <- Re(fft(product_FT, inverse=TRUE)/length(product_FT))
left_conv_2d <- fft(conv_frequency)
right_conv_2d <- product_FT
conv_error_2d <- max(Mod(left_conv_2d - right_conv_2d))
conv_error_2d
## [1] 9.094947e-13
plot_ly(z=~t(f_img), type="heatmap") %>%
  layout(title="Original PET Slice f(x,y)") %>% hide_colorbar()
plot_ly(z=~t(kernel_padded), type="heatmap") %>%
  layout(title="Padded Kernel g(x,y)") %>% hide_colorbar()
plot_ly(z=~t(conv_frequency), type="heatmap") %>%
  layout(title="IFT[FT(f) * FT(g)]") %>% hide_colorbar()
plot_ly(z=~t(fftshift(Mod(left_conv_2d - right_conv_2d))), type="heatmap") %>%
  layout(title="Difference: FT(f*g) - FT(f)*FT(g)") %>% hide_colorbar()
# Extract and analyze spatial-frequency distributions.
radial_frequency_distribution <- function(mat) {
  F <- fftshift(fft(mat))
  energy <- Mod(F)^2
  nr <- nrow(mat)
  nc <- ncol(mat)
  center_r <- floor(nr/2)+1
  center_c <- floor(nc/2)+1
  y <- matrix(rep(1:nr, nc), nrow=nr)
  x <- matrix(rep(1:nc, each=nr), nrow=nr)
  radius <- round(sqrt((x-center_c)^2 + (y-center_r)^2))
  df <- data.frame(radius=as.vector(radius), energy=as.vector(energy))
  aggregate(energy ~ radius, data=df, FUN=mean)
}
freq_dist_pet <- radial_frequency_distribution(imgMat)
plot_ly(x=~freq_dist_pet$radius, y=~freq_dist_pet$energy, type="scatter", mode="lines") %>%
  layout(title="Radial Spatial-Frequency Distribution of FDG-PET Slice",
         xaxis=list(title="Radial Frequency Index"),
         yaxis=list(title="Mean Fourier Energy", type="log"))

#Summary For this problem, we analyzed the dataset on 3D FDG-PET volume. Since the data set is 3D, we extracted a 2D slice of the volume, and used the slice for filtering, segmentation, Fourier transform, and feature extraction. The raw FDG-PET data image showed a smooth, “brain shaped” uptake pattern, where higher FDG uptake meant higher corresponding to brighter regions. The Gaussian blurred image looked very similar to the original PET slice, but the edges and uptake patterns appeared to be slightly smoother. The median filtered image also looks similar to the raw image, with small local variations reduced. The PET image was already smooth and did not contain much noise, the median filter did not cause a very dramatic change in appearance. In the Laplacian-filtered image, the outer boundary of the brain-shaped region was more pronounced. This shows that the Laplacian filter detects very rapid changes in intensity, and can detect edges and internal contours because of it. The new sharpened image preserved the main FDG uptake pattern, but slightly emphasized the boundaries more and local contrast. From the histogram, we see many low-intensity pixels, which likely correspond to background or low-FDG uptake regions. There was also a spread of high-intensity values from the PET uptake regions. There was also a spread of higher_intensity values from the PET uptake region. Upon applying Otsu-thresholding, which was used to seperate foreground PET uptake, to get one large connected forground region. The distance map was brightest at the foreground center and darker near the edges to show the distance of each foreground pixel from the background. The watershed segmentation produced one main large segmented region, rather than many small regions. The feature extraction output also showed this, as there was one main large area which was emphasized. Voronoi tesselation also produced one dominant region, because the watershed segmentation only produced one main seed region.

The Fourier Transform magnitude image showed a bright central region, which likely means that most of the image energy was concentrated at low spatial frequencies.This aligns with the visual appearance of the FDG-PET slice, which was was mainly dominated by smooth, broad uptake patterns, rather than very precise, local detail in structure. The phase image appeared very noisy and complicated. The inverse Fourier Transform reconstructed the original PET image, showing the Fourier representation preserved the image when no frequency components are removed. The low-pass Fourier filtering kept the central low-frequency components whilst removing higher-frequency components outside of the main central region. HIgh-pass Fourier filtering removed the central low frequency region, keeping high frequency components. After inverse FFT, the image lost the main smooth uptake pattern and became dominated by noise, indicating that the low-frequnecy components are important for preserving the overall shape and general features of the image from background structures, while high frequency components is more important for capturing finer details and precise features. When the central low frequency components below kx (1/4) max and ky (1/4) max were removed, the Fourier image showed a central dark rectangle, with the reconstructed image no longer resembling the original smooth PET slice. Instead, it showed many ring-like and edge-like patterns and artifacts. This shows that the central Fourier components are important for the broad, overall structural information in the image. When high frequency components on the periphery (above kx (1/2) max and ky (1/2) max) were removed, the image retained the central region and discarded outer frequencies but appeared smoother, showing that high frequency components on the periphery contribute to fine detail, rather than overall strucutres. For the linearity tests, we found that the error when comparing FT(af(x,y)+bg(x,y)) and aF(x,y) + bG(x,y) for 2 different slices was extremely small and close to zero, and the magnitude plots were virtually identical and difference plots were very small, which supports the linearity property for Fourier transform.When we tested the mulitplication property, we found that the error in our comparison between the 2 groups was also extremely small and close to zero. This supports the convolution theorem, which says that a convolution in image space corresponds to multiplication in image space. The final convolved image looked like a smoother version of the PET slice. Lastly, the radial spatial-frequency distribution showed thatthe fourier energy was highest at the smallest radial frequency values, as rapidly decreased as frequency increased. The shows that the FDG-PET slice is dominated by low-frequency information. Once again, this shows that low frequency components characterize broad, overall structures, while high frequency components characterize fine details. Overall, this PET image was primarily low frequency and smooth.

#actual problem 3

if(!requireNamespace("oro.nifti", quietly=TRUE)) install.packages("oro.nifti")
if(!requireNamespace("forecast", quietly=TRUE)) install.packages("forecast")
if(!requireNamespace("TTR", quietly=TRUE)) install.packages("TTR")
library(oro.nifti); library(forecast); library(TTR)

fMRIVol <- readNIfTI("fMRI_FilteredData_4D.nii", reorient=FALSE)
dim(fMRIVol)
## [1]  64  64  21 180
voxel_locations <- rbind(c(33,33,14), c(30,30,14), c(40,40,14), c(33,33,20))
colnames(voxel_locations) <- c("x","y","z")
voxel_locations
##       x  y  z
## [1,] 33 33 14
## [2,] 30 30 14
## [3,] 40 40 14
## [4,] 33 33 20
time_series_all <- list()
for(i in 1:nrow(voxel_locations)){
  time_series_all[[i]] <- as.numeric(fMRIVol[voxel_locations[i,1], voxel_locations[i,2], voxel_locations[i,3], 1:160])
}
names(time_series_all) <- paste0("Voxel_", voxel_locations[,1], "_", voxel_locations[,2], "_", voxel_locations[,3])

#ts.plot(time_series_all[[1]], time_series_all[[2]], time_series_all[[3]], time_series_all[[4]], col=1:4, #lty=1, xlab="Time", ylab="fMRI signal", main="Raw fMRI signals from different brain regions")
#legend("topright", legend=names(time_series_all), col=1:4, lty=1, cex=.7)


plot(time_series_all[[1]], type="l", col=1, lty=1, xlab="Time", ylab="fMRI signal", main="Raw fMRI signals from different brain regions", ylim=range(unlist(time_series_all)))
lines(time_series_all[[2]], col=2, lty=1)
lines(time_series_all[[3]], col=3, lty=1)
lines(time_series_all[[4]], col=4, lty=1)
legend("topright", legend=names(time_series_all), col=1:4, lty=1, cex=.7)

arima_results <- data.frame()
lm_results <- data.frame()

for(i in 1:length(time_series_all)){
time_series <- img <- time_series_all[[i]]; length(time_series)

ts <- ts(time_series, start=1, end=160, frequency=1)
ts.plot(ts, main=paste("Raw fMRI time-series:", names(time_series_all)[i]), xlab="Time", ylab="fMRI signal")
ts.epoch <- SMA(ts, n=10); plot(ts.epoch, main=paste("SMA:", names(time_series_all)[i]), xlab="Time", ylab="Smoothed fMRI signal")
acf(ts, lag.max = 20, main=paste("ACF:", names(time_series_all)[i]))
pacf(ts, lag.max = 20, main=paste("PACF:", names(time_series_all)[i]))

ts.train <- window(ts, start=1, end=140)
ts.test <- window(ts, start=141, end=160)

fit<-auto.arima(ts.train, approx=F, trace = F)
fit

ts.forecasts <- forecast(fit, h=20)
plot(ts.forecasts, include = 160, main=paste("ARIMA forecast:", names(time_series_all)[i]), xlab="Time", ylab="fMRI signal")
lines(141:160, as.numeric(ts.test), col="red", lwd=2)
legend("topright", legend=c("Forecast","Testing data"), col=c("blue","red"), lty=1, cex=.8)

myforecast <- forecast(fit, level=c(95), h=20)
plot(myforecast, main=paste("95% ARIMA forecast interval:", names(time_series_all)[i]))
lines(141:160, as.numeric(ts.test), col="red", lwd=2)

rmse <- sqrt(mean((as.numeric(ts.test)-as.numeric(ts.forecasts$mean))^2))
mae <- mean(abs(as.numeric(ts.test)-as.numeric(ts.forecasts$mean)))
arima_results <- rbind(arima_results, data.frame(voxel=names(time_series_all)[i], ARIMA_model=ts.forecasts$method, RMSE=rmse, MAE=mae))



stim <- rep(cbind(rep(0,10), rep(1,10)), 8); stim; length(stim)

#stim <- rep(cbind(rep(1,10), rep(-1,10)), 8); stim; length(stim)
X <- as.matrix(cbind(rep(1, 160),stim)); dim(X)
Y <- ts
lm_fit <- lm(Y ~ stim); summary(lm_fit)

lm_results <- rbind(lm_results, data.frame(voxel=names(time_series_all)[i], intercept=coef(summary(lm_fit))[1,1], stimulus_effect=coef(summary(lm_fit))[2,1], p_value=coef(summary(lm_fit))[2,4], r_squared=summary(lm_fit)$r.squared))

plot(Y, type="l", main=paste("Linear model stimulus fit:", names(time_series_all)[i]), xlab="Time", ylab="fMRI signal")
lines(fitted(lm_fit), col="red", lwd=2)
legend("topright", legend=c("Observed","LM fitted"), col=c("black","red"), lty=1, cex=.8)
}

arima_results
##            voxel                     ARIMA_model      RMSE       MAE
## 1 Voxel_33_33_14 ARIMA(0,0,0) with non-zero mean 176.81541 130.90857
## 2 Voxel_30_30_14 ARIMA(0,0,0) with non-zero mean  29.73311  22.65000
## 3 Voxel_40_40_14 ARIMA(0,0,0) with non-zero mean  45.37449  35.10000
## 4 Voxel_33_33_20 ARIMA(0,0,0) with non-zero mean 149.59699  96.93929
lm_results
##            voxel intercept stimulus_effect   p_value   r_squared
## 1 Voxel_33_33_14  14938.94        -26.6875 0.4083497 0.004330887
## 2 Voxel_30_30_14  11662.47         -7.8250 0.1689434 0.011942582
## 3 Voxel_40_40_14  12279.50         -3.9750 0.6213912 0.001547473
## 4 Voxel_33_33_20   2218.35         17.5375 0.4086296 0.004325728
best_arima_voxel <- arima_results[which.min(arima_results$RMSE),]
best_stimulus_voxel <- lm_results[which.min(lm_results$p_value),]
best_arima_voxel
##            voxel                     ARIMA_model     RMSE   MAE
## 2 Voxel_30_30_14 ARIMA(0,0,0) with non-zero mean 29.73311 22.65
best_stimulus_voxel
##            voxel intercept stimulus_effect   p_value  r_squared
## 2 Voxel_30_30_14  11662.47          -7.825 0.1689434 0.01194258
cat("The raw time-series plots show that different voxel locations can have different baseline signal levels, trends, and noise patterns. The ARIMA models were fit using time-points 1 through 140 and tested on time-points 141 through 160. Lower RMSE and MAE indicate better prediction of the held-out testing data. The linear model estimates the stimulus-response effect using the alternating ON/OFF stimulus vector. The stimulus_effect coefficient is the estimated change in fMRI signal associated with the ON versus OFF stimulus pattern. A positive stimulus_effect means the voxel signal is higher during ON stimulation, while a negative stimulus_effect means the voxel signal is lower during ON stimulation. Smaller p-values indicate stronger evidence that the voxel is responding to the stimulus. The voxel with the smallest p-value showed the strongest stimulus-response relationship among the selected locations.\n")
## The raw time-series plots show that different voxel locations can have different baseline signal levels, trends, and noise patterns. The ARIMA models were fit using time-points 1 through 140 and tested on time-points 141 through 160. Lower RMSE and MAE indicate better prediction of the held-out testing data. The linear model estimates the stimulus-response effect using the alternating ON/OFF stimulus vector. The stimulus_effect coefficient is the estimated change in fMRI signal associated with the ON versus OFF stimulus pattern. A positive stimulus_effect means the voxel signal is higher during ON stimulation, while a negative stimulus_effect means the voxel signal is lower during ON stimulation. Smaller p-values indicate stronger evidence that the voxel is responding to the stimulus. The voxel with the smallest p-value showed the strongest stimulus-response relationship among the selected locations.
cat("\nStimulus-response interpretation:\n")
## 
## Stimulus-response interpretation:
for(i in 1:nrow(lm_results)){
  cat("\n", lm_results$voxel[i], ":\n", sep="")
  cat("Estimated stimulus-response effect =", round(lm_results$stimulus_effect[i], 4), "\n")
  cat("p-value =", signif(lm_results$p_value[i], 4), "\n")
  cat("R-squared =", round(lm_results$r_squared[i], 4), "\n")
  if(lm_results$stimulus_effect[i] > 0){
    cat("Interpretation: This voxel has a positive stimulus-response effect, meaning the fMRI signal is higher during ON stimulation than OFF stimulation.\n")
  }else{
    cat("Interpretation: This voxel has a negative stimulus-response effect, meaning the fMRI signal is lower during ON stimulation than OFF stimulation.\n")
  }
  if(lm_results$p_value[i] < 0.05){
    cat("The p-value is below 0.05, so this voxel shows evidence of a statistically significant response to the stimulus.\n")
  }else{
    cat("The p-value is not below 0.05, so this voxel does not show strong evidence of a statistically significant response to the stimulus.\n")
  }
}
## 
## Voxel_33_33_14:
## Estimated stimulus-response effect = -26.6875 
## p-value = 0.4083 
## R-squared = 0.0043 
## Interpretation: This voxel has a negative stimulus-response effect, meaning the fMRI signal is lower during ON stimulation than OFF stimulation.
## The p-value is not below 0.05, so this voxel does not show strong evidence of a statistically significant response to the stimulus.
## 
## Voxel_30_30_14:
## Estimated stimulus-response effect = -7.825 
## p-value = 0.1689 
## R-squared = 0.0119 
## Interpretation: This voxel has a negative stimulus-response effect, meaning the fMRI signal is lower during ON stimulation than OFF stimulation.
## The p-value is not below 0.05, so this voxel does not show strong evidence of a statistically significant response to the stimulus.
## 
## Voxel_40_40_14:
## Estimated stimulus-response effect = -3.975 
## p-value = 0.6214 
## R-squared = 0.0015 
## Interpretation: This voxel has a negative stimulus-response effect, meaning the fMRI signal is lower during ON stimulation than OFF stimulation.
## The p-value is not below 0.05, so this voxel does not show strong evidence of a statistically significant response to the stimulus.
## 
## Voxel_33_33_20:
## Estimated stimulus-response effect = 17.5375 
## p-value = 0.4086 
## R-squared = 0.0043 
## Interpretation: This voxel has a positive stimulus-response effect, meaning the fMRI signal is higher during ON stimulation than OFF stimulation.
## The p-value is not below 0.05, so this voxel does not show strong evidence of a statistically significant response to the stimulus.

#Stimulus response interpretation The estimated stimulus-response effect is the coefficient of the stimulus variable in the linear model, and this coefficient estimates how much the fmri signal changes with alternating stimulus. A positive coefficient means that the voxel signal is higher when there is on stimulation vs off stimulation, and a negative coefficient means that the voxel signal was lower during the on stimulation rather than off stimulation

For Voxel_33_33_14: Has a negative stimulus effect. Signal was estimated to be lower during ON stimulation that Off stimulation,but was not statistically significant.

For Voxel_30_30_14: Has a negative stimulus effect. Signal was estimated to be lower during ON stimulation that Off stimulation,but was not statistically significant. However, the p-value was the lowest of all the voxels we studied.

For Voxel_40_40_14: Has a small negative stimulus effect. Signal was estimated to be lower during ON stimulation that Off stimulation,but was not statistically significant.

For Voxel_33_33_20: Has a positive stimulus effect, which means that the signal was estimated to be higher during the ON stimulation than OFF stimulation, but the effect was not statistically significant.

Overall, none of the voxels showed a statistically significant stimulus-repsonse effect, with a p-value threshold of 0.05. The R-squared values were all very small, suggesting that the on/off stimulus model explained very little of the variation in fMRI signals. This could have been because there may have been noise in the data, or that the specfific voxels which we selected were not strongly responisve to the stimulus.

#Summary We selected 4 voxels from the fMRI volume: Voxel_33_33_14, Voxel_30_30_14, Voxel_40_40_14, and Voxel_33_33_20. From the fMRI signal comparison, we saw that these all had very different baseline signal levels. Voxel_33_33_14 had the highest signal, at around 14,900 - 15,400. Voxel_30_30_14 and Voxel_40_40_14 had intermediate signal levels, at around 11,600 and 12,250 respectively, and voxel_33_33_20 had a much lower signal, around 2,000-2,600. As such, the baseline fMRI brain signal varies significantly between different voxels and brain regions.

From the raw time-series plots, we noticed noticeable fluctuations over time in all 4 voxels. The SMA-smoothed plots made slower trends easier to see, because they reduced noise in the short term.Voxel_33_33_14 and Voxel_33_33_20 showed large visible fluctuations, while Voxel_30_30_14 and Voxel_40_40_14 seemed to be more stable. The ACF and PACF plots did not show strong consistent autocorrelation patterns at most lags, suggesting that the time-series did not have a strong predictable structure from time.

For the ARIMA analysis, all 4 voxels were fit by an ARIMA(0,0,0) model, with a non-zero mean. The model mostly predicted future testing points using average level of training data. From the testing error, voxel_30_30_!4 had the best forcasting performance, with the lowest RMSE of 29.733 and MAE of 22.650. Voxel_40_40_14 has the next best performance, with RMSE = 45.374 and MAE =35.100. Voxel_33_33_14 had the weakest overall performance, with RMSE = 176.815 and MAE = 130.909. Voxel_33_33_20 had higher error, with RMSE = 149.597 and MAE = 96.939. ARIMA worked better for voxels with smaller variability/fluctuations, and worse for voxels with larger fluctuations.

For linear model analysis, coefficient represented estimated on/off difference in MRI signal. Voxel_33_33_!4 had a stimulus effect estimate of -26.6875, voxel_30_30_14 had one of -7.8250, and voxel_40_40_!4 had an estimate of -3.9750, and voxel_33_33_20 had a score of 17.5375. This means that Voxel_33_33_!4, voxel_30_30_14, and voxel_40_40_!4 had lower estimates signal for an on stimulus compared to an off stimulus, while voxel_33_33_20 had a higher estimated signal during on stimulus compared to an off stimulus. However, for all 4 voxels, the p-values were all >0.05 and as such none of the stimulus-response effects were statistically significant. The R-squared values were also all very small of 0.0015-0.0019 for all voxels, meaning that tis on/off model for stimulus which we are using is explaining very little of the variation in fMRI signals between voxels. Overall, the four voxels showed different baseline fMRI signal levels, but did not show strong statistically significant evidence of response to stimulus.This could have been because of noise in the singnals, or because the selected voxels were not appropriate. It could also be because we ignored HRF. In future analysis, we could try selecting and analyzing different voxels, and try to account for HRF in our analysis.