## ----warning=FALSE, message=FALSE--------------------------------------------- require(TCIU) require(doParallel) require(cubature) require(oro.nifti) require(magrittr) require(plotly) require(ggplot2) ## ----eval=FALSE, fig.align = "center", warning=FALSE, message=FALSE, fig.width = 9---- # # For this part of code, we comment it out and import the output plot already generated before to reduce time # # But this part of code can be run successfully. If you are interested, you can try it on your computer! # # # discrete Laplace Transform of sine # range_limit = 2 # x2 = seq(from = 0, to = range_limit, length.out = 50)[2:50] # # drop the first row to avoid real part value of 0 # y2 = seq(from = 0, to = range_limit, length.out = 50)[2:50] # # drop the first column to avoid imaginary part value of 0 # # # Recompute the LT(sin) discretized on lower-res grid # z2_grid = array(dim=c(length(x2), length(y2)))# x2 %o% y2 # # f_sin <- function(t) { sin(t) } # # # kime surface transform # # use parallel computing to speed up code # ncors = 2 # please choose the ncors according to the number of cores your PC has # # it is better that you increase the number of cores used for parallel computing if your computer allows # cl <- makeCluster(ncors) # registerDoParallel(cl) # F = list() # for (i in 1:length(x2) ){ # F[[i]] = # foreach(j = 1:length(y2), # .export='cubintegrate', # .packages='cubature') %dopar% { # TCIU::LT(FUNCT=f_sin, complex(real=x2[i], imaginary = y2[j])) # } # } # # stopCluster(cl) # F_vec = lapply(F, unlist) # z2_grid = unlist(do.call(rbind, F_vec)) # # # explicit form of Laplace Transform of sine # laplace_sine = function(p) { 1/(p^2 + 1) } # Exact Laplace transform of sin(x), continuous function # # XY = expand.grid(X=x2, Y=y2) # complex_xy = mapply(complex, real=XY$X,imaginary=XY$Y) # sine_z =laplace_sine(complex_xy) # dim(sine_z) = c(length(x2), length(y2))# dim(sine_z) # [1] 49 49 # # # make the two plots in the same plot to compare # lt_ilt_plotly = plot_ly(hoverinfo="none", showscale = FALSE)%>% # add_trace(z=Re(sine_z)-1, type="surface", surfacecolor=Im(sine_z)) %>% # add_trace(z = Re(z2_grid), type="surface", opacity=0.7, surfacecolor=Im(z2_grid) )%>% # layout(title = # "Laplace Transform, LT(sin()), Height=Re(LT(sin())), Color=Re(LT(sin())) \n Contrast Exact (Continuous) vs. # Approximate (Discrete) Laplace Transform", showlegend = FALSE) # # lt_ilt_plotly # ## ----fig.align="center", warning=FALSE, message=FALSE------------------------- sample_save[[1]] ## ----warning=FALSE, message=FALSE, fig.width = 9, fig.align = "center", eval=FALSE---- # # For this part of code, we comment it out and import the output plot already generated before to reduce time # # But this part of code can be run successfully. If you are interested, you can try it on your computer! # # # discrete Laplace Transform of sine # f_sin = function(t) { sin(t) } # lt_sine = function(z) TCIU::LT(f_sin, z) # # # inverse Laplace Transform on the lt_sine # # using parallel computing to speed up code # tvalsn <- seq(0, pi*2, length.out = 20) # cl <- makeCluster(ncors) # registerDoParallel(cl) # sinvalsn <- foreach(t=1:length(tvalsn), # .export='cubintegrate', # .packages='cubature') %dopar% { # TCIU::ILT(FUNCT=lt_sine, t=tvalsn[t]) # } # stopCluster(cl) # sinvalsn = unlist(sinvalsn) # # # make the plot of the result from ILT # # to see whether it still looks like sine # sinvalsn_df2 <- as.data.frame(cbind(Re=Re(sinvalsn),Im=Im(sinvalsn), # Sin=sin(tvalsn), time_points=tvalsn)) # lt_ilt_sine = ggplot(sinvalsn_df2, aes(x=time_points))+ # geom_line(aes(y=Re, color="Real"), linetype=1, lwd=2) + # geom_line(aes(y = Sin, color="Sin"), linetype=2, lwd=1) + # scale_color_manual(name="Index", # values = c("Real"="steelblue", "Sin"="darkred"))+ # labs(title = "Original fMRI Time-series f(t)=sin(t) and \n Reconstructed f'(t)=ILT(F)=ILT(discrete LT(f))", # subtitle = bquote("F" ~ "=" ~ "discrete LT(sine)")) + # xlab("Time") + ylab("fMRI Image Intensities (f and f')") + # theme_grey(base_size = 16) + # theme(legend.title = element_text(size=14, color = "black", face="bold"), # panel.grid.minor.y = element_blank(), # panel.grid.major.y = element_blank(), # plot.title = element_text(hjust = 0.5), # plot.subtitle = element_text(hjust = 0.5)) # # lt_ilt_sine ## ---- fig.align = "center", warning=FALSE, message=FALSE, fig.width = 9------- sample_save[[2]] ## ----warning=FALSE, message=FALSE, eval=FALSE--------------------------------- # x = seq(0, 2, length.out=50)[2:50]; y = seq(0, 2, length.out=50)[2:50]; # # do Kimesurface transform on sine function # z2_grid = kimesurface_transform(FUNCT = function(t) { sin(t) }, # real_x = x, img_y = y) # # # make the plot after Kimesurface transformation # surf_color <- atan2(Im(z2_grid), Re(z2_grid)) # colorscale = cbind(seq(0, 1, by=1/(length(x) - 1)), rainbow(length(x))) # magnitude <- (sqrt( Re(z2_grid)^2+ Im(z2_grid)^2)) # # # p <- plot_ly(hoverinfo="none", showscale = FALSE) %>% # add_trace(z = magnitude, # surfacecolor=surf_color, colorscale=colorscale, #Phase-based color # type = 'surface', opacity=1, visible=T) %>% # layout(title = "fMRI Kime-Surface, F=LT(fMRI) \n Height=Mag(F), Color=Phase(F)", showlegend = FALSE, # scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1, z=1.0) ) ) # 1:1:1 aspect ratio # p ## ----fig.align = "center", warning=FALSE, message=FALSE, fig.width = 9-------- sample_save[[3]] ## ----warning=FALSE, message=FALSE, eval = FALSE, fig.align = "center", fig.width = 9---- # # load the 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) # # dimensions: 64 x 64 x 21 x 180 ; 4mm x 4mm x 6mm x 3 sec # # # # extract a set of time series data from a voxel of fMRI data # xA_fMRI_1D_x20_y20_z11 <- fMRIVolume[20, 20, 11, ] # # # get the smooth function of this time series data # # Instead of using the extremely noisy fMRI data and avoiding integration problems, # # smooth "f" and use the **smooth version, f** # time_points <- seq(0+0.001, 2*pi, length.out = 180) # f <- smooth.spline(ceiling((180*time_points)/(2*pi)), # xA_fMRI_1D_x20_y20_z11, df = 10)$y # # # Define the f(t)=smooth(fMRI)(t) signal as a function of real time 0 2*pi) { return ( 0 ) } else { # return ( f[ceiling((180*t)/(2*pi))] ) # } # } # # # # do the Kimesurface transform # ncors = 8 # # please choose the ncors according to the number of cores your PC has # x = seq(0, 2, length.out=50)[2:50]; y = seq(0, 2, length.out=50)[2:50]; # # do Kimesurface transform on sine function # z2_grid_fmri = kimesurface_transform(FUNCT = fmri_funct, # glb_para="f", # real_x = x, img_y = y, # parallel_computing = TRUE, # ncor=ncors) # # # make the plot of function after Kimesurface transformation # surf_color <- atan2(Im(z2_grid_fmri), Re(z2_grid_fmri)) # colorscale = cbind(seq(0, 1, by=1/(length(x) - 1)), rainbow(length(x))) # magnitude <- (sqrt( Re(z2_grid_fmri)^2+ Im(z2_grid_fmri)^2)) # # # p_fmri <- plot_ly(hoverinfo="none", showscale = FALSE) %>% # add_trace(z = magnitude, # surfacecolor=surf_color, colorscale=colorscale, # Phase-based color # type = 'surface', opacity=1, visible=T) %>% # layout(title = "fMRI Kime-Surface, F=LT(fMRI) \n Height=Mag(F), Color=Phase(F)", showlegend = FALSE, # scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1, z=1.0) ) ) # 1:1:1 aspect ratio # p_fmri ## ----fig.align = "center", warning=FALSE, message=FALSE, fig.width = 9-------- sample_save[[4]] ## ----warning=FALSE, message=FALSE, fig.width = 9, fig.align = "center", eval=FALSE---- # time_points <- seq(0+0.001, 2*pi, length.out = 180) # inv_data = inv_kimesurface_transform(time_points, z2_grid) # inv_data = inv_kimesurface_transform(time_points, z2_grid,num_length = 23, # m=1, msg=TRUE) # # time_Intensities_ILT_df2 <- as.data.frame(cbind(Re=scale(Re(inv_data$Smooth_Reconstruction)), # Im=scale(Re(inv_data$Raw_Reconstruction)), # fMRI=scale(Re(sin(time_points))), # time_points=time_points)) # colnames(time_Intensities_ILT_df2) = c("Smooth Reconstruction", # "Raw Reconstruction", # "Original Sin", "time_points") # df = reshape2::melt(time_Intensities_ILT_df2, id.var = "time_points") # pppp<-ggplot(df, aes(x = time_points, y = value, colour = variable)) + # geom_line(linetype=1, lwd=3) + # ylab("Function Intensity") + xlab("Time") + # theme(legend.position="top")+ # labs(title= bquote("Comparison between" ~ "f(t)=Smooth(Sin)(t)" ~ "and Smooth(ILT(LT(Sin)))(t); Range [" ~ 0 ~":"~ 2*pi~"]")) ## ----fig.align = "center", warning=FALSE, message=FALSE, fig.width = 9-------- sample_save[[5]]