| SOCR ≫ | DSPA ≫ | DSPA2 Topics ≫ |
DSPA2: Data Science and Predictive Analytics (UMich HS650)
Appendix 10: Interpolation of Irregularly Sampled Spatiotemporal Data
SOCR/MIDAS (Ivo Dinov)
July 2025
This DSPA2 appendix demonstrates useful protocols for spatiotemporal interpolation, regularization, and uniform sampling of univariate and multivariate, spatial and temporal (longitudinal) data. The process of interpolation used to regularize irregularly sampled data, impute missing values, and facilitate homologies between cases, studies, cohorts, participants, populations, clusters, etc.
1 Motivation of Interpolation of Irregularly Sampled Data
Let’s use the R package akima to illustrate the regularization of a process, or function \(z=f(x,y)\), that is only observed on irregularly sampled points \((x,y)\in \mathbb{R}^2\), distributed throughout the 2D \((x,y)\in \mathbb{R}^2\) plane. Note that each of the axes can represent space or time, one of the axes can be time and the other one can be space, or alternatively the signal can represent irregularly sampled kimesurface, \(z=f(x,y)\) over complex time (kime), \(\kappa\equiv z\in\mathbb{C}\cong\mathbb{R}^2\ni (x,y)\), where the kime-magnitude \(t=|\kappa|=\sqrt{x^2+y^2}\) and the kime-phase \(\phi\) are irregularly sampled and \(\kappa=t e^{i\phi}\in\mathbb{C}\).
2 Univariate Process Interpolation
In general, function interpolation returns a list of points smoothly interpolating the observed (irregularly-sampled) data points.
2.1 Regularly spaced sample
x <- 1:10
y <- c(rnorm(5), c(1,1,1,1,3))
xnew <- seq(-1, 11, 0.1)
# plot(x, y, ylim=c(-3, 3), xlim=range(xnew))
# lines(spline(x, y, xmin=min(xnew), xmax=max(xnew), n=200), col="blue")
spl <- spline(x, y, xmin=min(xnew), xmax=max(xnew), n=200)
plot_ly(x=x, y=y, type="scatter", mode="markers", name="Raw Obs") %>%
add_trace(x=spl$x, y=spl$y, type="scatter", mode="lines", name="Spline Model") %>%
layout(title="Spline Modeling of a 1D Process (Uniform Sampling)",
legend = list(orientation = 'h')) %>%
hide_colorbar()2.2 Irregularly spaced sample
x <- sort(runif(10, max=10))
y <- c(rnorm(5), c(1,1,1,1,3))
spl_1 <- aspline(x, y, xnew)
spl_2 <- aspline(x, y, xnew, method="improved")
spl_3 <- aspline(x, y, xnew, method="improved", degree=10)
plot_ly(x=x, y=y, type="scatter", mode="markers", name="Raw Obs") %>%
add_trace(x=spl$x, y=spl$y, type="scatter", mode="lines", name="Spline Model") %>%
add_trace(x=spl_1$x, y=spl_1$y, type="scatter", mode="lines", name="Spline Model 1") %>%
add_trace(x=spl_2$x, y=spl_2$y, type="scatter", mode="lines", name="Spline Model 2") %>%
add_trace(x=spl_3$x, y=spl_3$y, type="scatter", mode="lines", name="Spline Model 3") %>%
layout(title="Four Spline Models of a 1D Process (Uniform Sampling)",
legend = list(orientation = 'h')) %>%
hide_colorbar()2.3 Another Example - Irregular Sampling
x <- c(-3, -2, -1, 0, 1, 2, 2.5, 3)
y <- c( 0, 0, 0, 0, -1, -1, 0, 2)
spl <- spline(x, y, n=200)
spl_1 <- aspline(x, y, n=200)
spl_2 <- aspline(x, y, n=200, method="improved")
spl_3 <- aspline(x, y, n=200, method="improved", degree=10)
plot_ly(x=x, y=y, type="scatter", mode="markers", name="Raw Obs") %>%
add_trace(x=spl$x, y=spl$y, type="scatter", mode="lines", name="Spline Model") %>%
add_trace(x=spl_1$x, y=spl_1$y, type="scatter", mode="lines", name="Spline Model 1") %>%
add_trace(x=spl_2$x, y=spl_2$y, type="scatter", mode="lines", name="Spline Model 2") %>%
add_trace(x=spl_3$x, y=spl_3$y, type="scatter", mode="lines", name="Spline Model 3") %>%
layout(title="Four Spline Models of a 1D Process (Irregular Sampling)",
legend = list(orientation = 'h')) %>%
hide_colorbar()2.4 Bicubic Bivariate Interpolation for Data Irregularly smapled on a Rectangular Grid \((x,y)\)
data(akima760)
# interpolate at the diagonal of the grid [0,8]x[0,10]
akima.bic <- bicubic(akima760$x, akima760$y, akima760$z, seq(0,8,length=50), seq(0,10,length=50))
x <- sqrt(akima.bic$x^2+akima.bic$y^2)
y <- akima.bic$z
plot_ly(x=x, y=y, type="scatter", mode="markers+lines", name="Raw Obs") %>%
layout(title="Bicubic Bivariate Interpolation (Irregular Sampling)",
xaxis=list(title="Argument Magnitude (||(x,y)||)"),
yaxis=list(title="Functional Value (z=f(x,y))"),
legend = list(orientation = 'h')) %>%
hide_colorbar()x <- akima760$x
y <- akima760$y
z <- akima760$z
plot_ly(x=x, y=y, z=z, type="surface", name="Raw Obs") %>%
layout(title="Bicubic Bivariate Interpolation (Irregular Sampling)",
xaxis=list(title="Argument Magnitude (||(x,y)||)"),
yaxis=list(title="Functional Value (z=f(x,y))"),
legend = list(orientation = 'h')) %>%
hide_colorbar()2.5 Bicubic-grid: Bivariate Interpolation for Data on a Rectangular Grid
Interpolation of a bivariate function, \(z=f(x,y)\), on a rectangular grid in the
x-y plane, based on the revised Akima method
bicubic.grid(x,y,z,xlim,ylim,dx,dy).
data(akima760)
# interpolate at a grid [0,8]x[0,10]
akima.bic <- bicubic.grid(akima760$x, akima760$y, akima760$z, c(0,8),c(0,10),0.1,0.1)
# ==================
# interp: Gridded Bivariate Interpolation for Irregular Data
# ==========================
# implement bivariate interpolation onto a grid for irregularly spaced input data.
# Bilinear or bicubic spline interpolation is applied using different versions of algorithms from Akima
data(akima)
# plot(y ~ x, data = akima, main = "akima example data")
# with(akima, text(x, y, formatC(z,dig=2), adj = -0.1))
values <- formatC(akima$z,dig=2)
plot_ly(as.data.frame(akima), x=~x, y=~y, type="scatter", mode="markers", name="Raw Obs") %>%
add_text(x=~x, y=~y, text=values, name=values, textposition = 'middle top') %>%
layout(title="Bicubic Bivariate Interpolation (Irregular Sampling)",
xaxis=list(title="x"),
yaxis=list(title="y"),
legend = list(orientation = 'h')) %>%
hide_colorbar()## linear interpolation
akima_df <- as.data.frame(akima)
akima_df <- akima_df %>% replace(is.na(.), 0)
# interpolation
akima_lin_interpol <-
interp(akima_df$x, akima_df$y, akima_df$z,
xo=seq(min(akima_df$x), max(akima_df$x), length = 200),
yo=seq(min(akima_df$y), max(akima_df$y), length = 200))
akima_lin_interpol$z <- akima_lin_interpol$z %>% replace(is.na(.), 0)
df <- as.data.frame(akima)
plot_ly() %>% add_contour(x=~akima_lin_interpol$x, y=~akima_lin_interpol$y,
z=~akima_lin_interpol$z, name="Bivariate Linear Model") %>%
add_trace(x=~df$x, y=~df$y, type="scatter", mode="markers", name="Raw Obs") %>%
layout(title="Linear Bivariate Interpolation (Irregular Sampling)",
xaxis=list(title="x"),yaxis=list(title="y"),
legend = list(orientation = 'h')) %>%
hide_colorbar()## use finer grid to increase interpolation-model smoothness
akima.smooth <-
with(akima, interp(x, y, z, xo=seq(0,25, length=200), yo=seq(0,20, length=200)))
akima.smooth$z <- akima.smooth$z %>% replace(is.na(.), 0)
plot_ly() %>% add_contour(x=~akima.smooth$x, y=~akima.smooth$y,
z=~akima.smooth$z, name="Smooth Bivariate Model") %>%
add_trace(x=~df$x, y=~df$y, type="scatter", mode="markers", name="Raw Obs") %>%
layout(title="Bivariate Interpolation (Finer 2D Grid)",
xaxis=list(title="x"),yaxis=list(title="y"),
legend = list(orientation = 'h')) %>%
hide_colorbar()# Sparse-data Interpolation using only 10 points (within convex hull)
df <- as.data.frame(akima)
subset_int <- sample(nrow(df), 10)
df <- df[subset_int, ]
akima.part <- with(df, interp(x, y, z))
akima.part$z <- akima.part$z %>% replace(is.na(.), 0)
df$z <- df$z %>% replace(is.na(.), 0)
plot_ly() %>% add_contour(x=~akima.part$x, y=~akima.part$y, z=~akima.part$z, name="Space (10-point-based) Model") %>%
add_trace(x=~df$x, y=~df$y, type="scatter", mode="markers", name="Raw Obs") %>%
layout(title="Sparce-data (n=10) Bivariate Interpolation",
xaxis=list(title="x"),yaxis=list(title="y"),
legend = list(orientation = 'h')) %>%
hide_colorbar()3 Bivariate Case
3.1 Linear Interpolation
3.1.1 Gridded Bivariate Interpolation of Irregularly-Sampled Data
# sample of akima
data(akima)
dat <- as.data.frame(akima)
# Display raw data observations
plot_ly(dat, type="scatter3d", mode="markers", x=~x, y=~y, z=~z) # %>% add_surface(x=~x, y=~y, z=~z)
# interpolation
data_grid_interpol <-
interp(dat$x, dat$y, dat$z,
xo=seq(min(dat$x), max(dat$x), length = 200),
yo=seq(min(dat$y), max(dat$y), length = 200))
plot_ly(dat, type="scatter3d", mode="markers", x=~x, y=~y, z=~z,
name="Observed Sample Data") %>%
add_surface(x=~data_grid_interpol$x, y=~data_grid_interpol$y,
z=~data_grid_interpol$z, name="Model Manifold",
contours = list(x = list(show = TRUE, start = 0.5,
end=25, size=5, color='gray'),
z = list(show = TRUE, start = 0.5,
end = 60, size = 5))) %>%
layout(title="Linear Modeling & Grid-Interpolation of 2D Processes",
legend = list(orientation = 'h')) %>%
hide_colorbar()3.1.2 Pointwise Bivariate Interpolation of Irregularly-Sampled Data:
data_pointwise_interpol <-
interpp(dat$x, dat$y, dat$z, runif(20, min(dat$x), max(dat$x)),
runif(20, min(dat$y), max(dat$y)))
plot_ly(dat, type="scatter3d", mode="markers", x=~x, y=~y, z=~z,
marker = list(color = "red", size=3), opacity=0.6,
name="Observed Sample Data") %>%
add_trace(type="scatter3d", mode="markers",
x=~data_pointwise_interpol$x,
y=~data_pointwise_interpol$y,
z=~data_pointwise_interpol$z,
marker = list(color = "blue", size=6), opacity=0.6,
name="Pointwise Interpolation"
) %>%
add_surface(x=~data_grid_interpol$x, y=~data_grid_interpol$y,
z=~data_grid_interpol$z, name="Model Manifold",
opacity=0.2) %>%
layout(title=
"Linear Modeling & Pointwise-Interpolation of 2D Processes",
legend = list(orientation = 'h')) %>%
hide_colorbar()3.2 Nonlinear (Spline) Interpolation
3.2.1 Bivariate Cubic-Spline Grid Interpolation
data_cubspline_interpol <-
interp(dat$x, dat$y, dat$z,
xo=seq(min(dat$x), max(dat$x), length = 200),
yo=seq(min(dat$y), max(dat$y), length = 200),
linear = FALSE, extrap = TRUE)
plot_ly(dat, type="scatter3d", mode="markers", x=~x, y=~y, z=~z,
name="Observed Sample Data") %>%
add_surface(x=~data_cubspline_interpol$x,
y=~data_cubspline_interpol$y,
z=~data_cubspline_interpol$z,
name="Spline Model Manifold") %>%
layout(title="Nonlinear Spline Model Interpolation of 2D Processes",
legend = list(orientation = 'h')) %>%
hide_colorbar()3.2.2 Bivariate Cubic-Spline Pointwise Interpolation
data_pointwise_interpol <-
interpp(dat$x, dat$y, dat$z, runif(20, min(dat$x), max(dat$x)),
runif(20, min(dat$y), max(dat$y)),
linear = FALSE, extrap = TRUE)
plot_ly(dat, type="scatter3d", mode="markers", x=~x, y=~y, z=~z,
marker = list(color = "red", size=3), opacity=0.6,
name="Observed Sample Data") %>%
add_trace(type="scatter3d", mode="markers",
x=~data_pointwise_interpol$x,
y=~data_pointwise_interpol$y,
z=~data_pointwise_interpol$z,
marker = list(color = "blue", size=6), opacity=0.6,
name="Pointwise Spline Interpolation") %>%
add_surface(x=~data_grid_interpol$x, y=~data_grid_interpol$y,
z=~data_grid_interpol$z, name="Model Manifold",
opacity=0.2) %>%
layout(title=
"Spline Modeling & Pointwise-Interpolation of 2D Processes",
legend = list(orientation = 'h')) %>%
hide_colorbar()4 Upsampling Data using Sinc Interpolation
In poractice, increasing the resolution of multidimensional data is a vital preprocessing step prior to subsequent analysis and visualization. Sinc interpolation is a powerful method for upsampling data, as it can perfectly reconstruct a bandlimited signal from its samples. The core principle behind multi-dimensional sinc interpolation is the separable application of the one-dimensional sinc interpolation formula to each dimension of the array. The 1D sinc interpolation formula is defined by
\[g(t) = \sum_{i=-\infty}^{\infty} g(iT) \cdot \operatorname{sinc}\left(\frac{t - iT}{T}\right),\] where \(g(t)\) is the interpolated value at time \(t\), \(g(iT)\) are the known sample values at intervals of \(T\), and the normalized sinc function is given by
\[\operatorname{sinc}(x) = \frac{\sin(\pi x)}{\pi x}.\]
Let’s consider a 2D example, first interpolating along the rows of the data matrix and then along the columns of the resulting matrix to produce the final high-resolution array.
# Functions for performing 2D sinc interpolation
# to increase the resolution of a 2D array. A working example with
# simulated data and visualizations is included.
# Required Libraries
# install.packages(c("ggplot2", "reshape2"))
library(ggplot2)
library(reshape2)
# 1. Sinc Function
# Defines the normalized sinc function.
sinc <- function(x) {
ifelse(x == 0, 1, sin(pi * x) / (pi * x))
}
# 2. 1D Sinc Interpolation Function
# Performs one-dimensional sinc interpolation.
#
# Args:
# x_new: A numeric vector of new points to interpolate.
# x_old: A numeric vector of the original sample points.
# y_old: A numeric vector of the original sample values.
#
# Returns:
# A numeric vector of the interpolated values at x_new.
sinc_interp_1d <- function(x_new, x_old, y_old) {
sapply(x_new, function(t) {
sum(y_old * sinc((t - x_old)))
})
}
# 3. 2D Sinc Interpolation Function
# Performs two-dimensional sinc interpolation on a matrix.
#
# Args:
# m_old: The original low-resolution matrix.
# x_old: A numeric vector of the original x-coordinates.
# y_old: A numeric vector of the original y-coordinates.
# x_new: A numeric vector of the new high-resolution x-coordinates.
# y_new: A numeric vector of the new high-resolution y-coordinates.
#
# Returns: A high-resolution matrix with the interpolated values.
sinc_interp_2d <- function(m_old, x_old, y_old, x_new, y_new) {
# Interpolate along the rows (x-dimension)
m_interp_rows <- t(apply(m_old, 1, function(row_data) {
sinc_interp_1d(x_new, x_old, row_data)
}))
# Interpolate along the columns (y-dimension)
m_interp_cols <- apply(m_interp_rows, 2, function(col_data) {
sinc_interp_1d(y_new, y_old, col_data)
})
return(m_interp_cols)
}
# 4. Generate Simulated Low-Resolution Data
low_res_n <- 10
x_low <- seq(-5, 5, length.out = low_res_n)
y_low <- seq(-5, 5, length.out = low_res_n)
# Create a grid and a sample function
grid_low <- expand.grid(X = x_low, Y = y_low)
z_low_vec <- with(grid_low, sinc(sqrt(X^2 + Y^2) / pi))
z_low_matrix <- matrix(z_low_vec, nrow = low_res_n, ncol = low_res_n, byrow = TRUE)
# 5. Define High-Resolution Grid and Perform Interpolation
high_res_n <- 100
x_high <- seq(-5, 5, length.out = high_res_n)
y_high <- seq(-5, 5, length.out = high_res_n)
# Perform 2D sinc interpolation
z_high_matrix <- sinc_interp_2d(z_low_matrix, x_low, y_low, x_high, y_high)
# 6. Visualize the Results
# Prepare data for ggplot
df_low <- data.frame(expand.grid(X = x_low, Y = y_low), Z = as.vector(z_low_matrix))
df_high <- data.frame(expand.grid(X = x_high, Y = y_high), Z = as.vector(z_high_matrix))
# Plot the original low-resolution data
p_low <- ggplot(df_low, aes(x = X, y = Y, fill = Z)) +
geom_tile() +
scale_fill_viridis_c() +
labs(title = "Original Low-Resolution Data", x = "X", y = "Y") +
theme_minimal()
# Plot the interpolated high-resolution data
p_high <- ggplot(df_high, aes(x = X, y = Y, fill = Z)) +
geom_tile() +
scale_fill_viridis_c() +
labs(title = "High-Resolution Data (Sinc Interpolation)", x = "X", y = "Y") +
theme_minimal()
# Display the plots
print(p_low)
In the above code,
sinc(x): This function calculates the value of the normalized sinc function for a given inputx. It handles the special case wherexis zero, returning 1 to avoid division by zero.sinc_interp_1d(...): This function performs the one-dimensional sinc interpolation. It takes the new points (x_new), the original sample points (x_old), and the original sample values (y_old) as input. It then iterates through each new point and calculates the interpolated value using the sinc interpolation formula.sinc_interp_2d(...): This is the core function for the two-dimensional interpolation. It first applies thesinc_interp_1dfunction to each row of the input matrixm_old. Theapplyfunction withMARGIN = 1is used for this purpose. The result is a matrix that is interpolated in the x-dimension. Then, it applies thesinc_interp_1dfunction to each column of this new matrix (MARGIN = 2), effectively performing the interpolation in the y-dimension.Data Simulation: A low-resolution 10x10 grid is created, and the value at each point is calculated using a 2D sinc-like function to generate a smooth surface for demonstration.
Interpolation: A high-resolution 100x100 grid is defined. The
sinc_interp_2dfunction is then called to compute the values on this new, denser grid.Visualization: The
ggplot2package (orplotly) may be used to create heatmaps of both the original low-resolution data and the resulting high-resolution data, allowing for a clear visual comparison of the upsampling effect. The interpolated image appears much smoother and more detailed, showcasing the effectiveness of the sinc interpolation method.
5 Image Super-resolution via Sinc Interpolation
Sinc
interpolation is the optimal strategy to increase image resolution
by a specified factor without relying on additional information. The
examples below utilize the imager and pracma
R packages to demonstrate resampling 2D images and
increasing their resolution by feeding zeo (\(0\)) rows and columns in the Fourier space
representation of the input images, which is the core of sinc
interpolation (aka Whittaker–Shannon interpolation).
This demonstration performs the following steps
- Define the sinc function.
- Create a new grid for the rescaled image.
- Initialize a new matrix for the rescaled image.
- Apply the sinc interpolation formula to calculate the pixel values for the new image.
- Load an example image and applies the sinc interpolation with a specified scaling factor.
- Save and plot the rescaled image.
The zero-padding the Fourier transform of a signal preserves the original signal’s frequency components while introducing (empty, information neutral) higher frequency components that correspond to the interpolated samples. Adjusting the scaling factor increases, or decreases, the image resolution (subsample or super-resolution) the original image.
# Install the necessary packages:
# install.packages("imager")
# install.packages("pracma")
library(imager)
library(plotly)
library(jpeg)
# 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")
}
}
# Define the sinc interpolation function using FFT
sinc_interp_fft <- function(image, factor) {
# Get original image dimensions
orig_x <- dim(image)[1]
orig_y <- dim(image)[2]
# Apply FFT to the original image
fft_image <- fft(image)
# Zero-padding in the Fourier domain
pad_x <- orig_x * factor
pad_y <- orig_y * factor
pad_fft <- matrix(0, nrow = pad_x, ncol = pad_y)
# Calculate the center of the original FFT
center_x <- ceiling(orig_x / 2)
center_y <- ceiling(orig_y / 2)
# Copy the original FFT to the center of the padded FFT matrix
pad_fft[1:center_x, 1:center_y] <- fft_image[1:center_x, 1:center_y]
pad_fft[(pad_x - center_x + 1):pad_x, 1:center_y] <-
fft_image[(center_x + 1):orig_x, 1:center_y]
pad_fft[1:center_x, (pad_y - center_y + 1):pad_y] <-
fft_image[1:center_x, (center_y + 1):orig_y]
pad_fft[(pad_x - center_x + 1):pad_x, (pad_y - center_y + 1):pad_y] <-
fft_image[(center_x + 1):orig_x, (center_y + 1):orig_y]
# Inter-splice zeros uniformly throughout the k-space
# pad_fft[seq(1, pad_x, by = factor), seq(1, pad_y, by = factor)] <- fft_image
# Apply inverse FFT to get the interpolated image
interp_image <- fft(pad_fft, inverse = TRUE)
# Normalize the image
interp_image <- Re(interp_image) * (factor^2)/(length(interp_image))
return(list(interp_image=as.cimg(interp_image), fft_image=fft_image, pad_fft=pad_fft))
}
# Load an example image
# image <- load.image(system.file("extdata/parrots.png", package = "imager"))
# Read Image in and display it
# Either download the JPG image locally from https://umich.instructure.com/courses/38100/files/1627149/download?download_frd=1 and then read it in
# mri <- readJPEG("C:/Users/user/Desktop/MRI_ImageHematoma.jpg")
# Or directly load it from the web
pathToMRI <- tempfile(fileext = ".jpg")
download.file("https://umich.instructure.com/courses/38100/files/1627149/download?download_frd=1",
mode = "wb", pathToMRI)
image <- readJPEG(pathToMRI)
# load.image(system.file(pathToMRI, package = "imager"))
# Convert the image to grayscale for simplicity
# gray_image <- grayscale(image[,,1])
# Perform sinc interpolation with a scaling factor of 2, and pull the result list of images
res <- sinc_interp_fft(t(image[,,1]), 2)
rescaled_image <- res$interp_image
fft_image <- res$fft_image
pad_fft <- res$pad_fft
# image(rescaled_image[,,1,1] )
# spacetime domain images
plt_orig <- plot_ly(z=apply(t(image[,,1]), 1, rev), type="heatmap", name="Original Img") %>%
layout(title="Original Image") %>% hide_colorbar()
plt_SuperResSincInterpol <- plot_ly(z=apply(rescaled_image[,,1,1], 1, rev),
type="heatmap", name="Sinc Interpolated Image") %>%
layout(title="Sinc Interpolated Image (Zoom-Factor=2*2)") %>% hide_colorbar()
# Fourier k-space images
mag_FFT <- sqrt(Re(fftshift(fft_image))^2 + Im(fftshift(fft_image))^2)
plt_OrigFFT <- plot_ly(z=log(apply(mag_FFT, 1, rev)), type="heatmap", name="FFT") %>%
layout(title="FFT Log-Magnitute of Original Image (center-shifted k-space)") %>%
hide_colorbar()
mag_FFT_padded <- sqrt(Re(fftshift(pad_fft))^2 + Im(fftshift(pad_fft))^2)
plt_FFTpadded <- plot_ly(z=log(apply(mag_FFT_padded, 1, rev)+0.001),
type="heatmap", name="Padded FFT") %>%
layout(title="Sinc FFT-padded Log-Magnitute of Image in k-space") %>%
hide_colorbar()
subplot(plt_orig, plt_SuperResSincInterpol, plt_OrigFFT, plt_FFTpadded, nrows = 2) %>%
layout(title = 'Sinc Interpolation: Original and Superresolution Images and their Fourier-space counterparts')